Abstract

LEAK2D (L2D) is a new, two-dimensional, dual-permeability model for the simulation of preferential water flow in the vadose zone, allowing for the continuous exchange of water between the matrix and the fracture domain. It is based on the two-dimensional Richards equation for the simulation of flow in the matrix domain and on the kinematic wave equation for the simulation of flow in the fracture domain. The Richards equation is solved by a combination of the Alternating Direction Implicit method and the Douglas-Jones predictor-corrector method. This combination leads to a very efficient, stable, and time-consuming method. A variable time step is used by which any instability of the numerical solution is avoided. The water transfer from the fracture to the matrix domain is estimated as a first-order approximation of the water diffusion equation. The model was used to satisfactorily simulate preferential flow under an extreme rainfall/irrigation event. The exchange of water between the two domains depends on parameters which have physical meaning; however, their exact values are difficult to be determined or measured. Based on the most common values of these parameters found in the literature, a sensitivity analysis was performed to define their effect on the output of the model.

1 Introduction

1.1 Background

The generic term ‘preferential flow’ was defined by Guo and Lin (2018) and Hendrickx and Flury (2001) as all phenomena where water and solutes move along certain pathways, while bypassing a fraction of the porous matrix. Preferential flow results in irregular wetting patterns in the soil profile as a direct consequence of water moving faster in certain parts of the soil profile than in others (Šimůnek et al. 2003; Nimmo 2021). The result of these irregular water patterns is the irregular concentration of the solutes in the soil. Preferential flow influences the leaching of contaminants to groundwater because the biologically and chemically-reactive surface soil layers can be quickly bypassed (Larsbo et al. 2005).

The preferential flow phenomenon is observed in aggregated fields, rocky, and homogeneous soils (Gerke 2006; Guo and Lin 2018; Pluer et al. 2020). Preferential flow is increased with irrigation intensity and is stronger in surface irrigation than in sprinkling.

The need for a preferential flow simulation led to the development of models in the last three decades, which are based on the theory of dual/multiple porosity (van Genuchten and Wierenga 1976) and on the theory of dual/multiple permeability (Gerke and van Genuchten 1993a). Dual-porosity and dual-permeability models both assume that the porous medium consists of two interacting regions, one associated with the inter-aggregate macropore, or fracture system, and one comprising micropores (or intra-aggregate pores) inside soil aggregates, or inside the rock matrix. In dual porosity models, it is assumed that the water flow is restricted only in the fracture domain, while the water in the matrix domain remains immobile (Šimůnek et al. 2003).

The dual permeability approach is the most appropriate to describe preferential flow phenomena avoiding the complexity of the multiple permeability models. The dual-permeability models differ mainly in the methodology for the determination of water movement in the fracture domain, and also in the description of the interaction process between the two domains (Šimůnek et al. 2003).

Another point that deserves further investigation in dual-permeability models is the interaction between the two domains as well as the equations that describe it. The equations that have been presented in the literature are empirical and in most cases contain parameters that have no physical meaning (Šimůnek et al. 2003). The role of these parameters in the description of the water dynamics in the soil profile is significant and needs to be clarified through their sensitivity analysis.

1.2 Purpose and scope

The primary objective of the present work is the development of a new, two-dimensional dual-permeability model for the preferential flow of the water in the vadose zone. The model adopts the advantages of the models of the above-mentioned two categories, while avoiding their drawbacks. It combines the two-dimensional Richards equation only in the matrix domain (where this equation is always valid) and the kinematic wave equation in the fracture domain (where the validity of Richards equation under certain circumstances can be questionable), allowing the water transfer between the two domains. The Richards equation is solved by a combination of Alternating Direction Implicit (A.D.I.) (An et al. 2011) method and the Douglas and Jones predictor-corrector method (Douglas and Jones 1963; Babajimopoulos 2000). This new method seems to be accurate and fast, without stability problems. Secondarily, a sensitivity analysis of the main input parameters of the fracture domain equations is accomplished to define their effect on the output of the model (water distribution in the soil slab and deep drainage).

2 Materials and methods

2.1 Model conceptualization

Leak2D (L2D) is a novel, two-dimensional, dual-permeability model to describe the preferential flow of water in the vadose zone. As in Gerke and van Genuchten (1993a) it is assumed that the soil can be separated into two pore domains. Each one of these domains has different hydraulic properties and different degrees of saturation. Also, water flow in these two domains is described by different equations. These two pore domains are: a) the soil matrix (subscript m), and b) the fracture (subscript f). During the execution of the model, the equations in each domain are solved sequentially in every time step. The interaction of the flow between the two domains takes place in the case where some appropriate conditions are met. This approach is an alteration of the Gerke and van Genuchten (1993a) concept for water flow in structured soils.

As in Gerke and van Genuchten (1993a) the relative volumetric proportion of the fracture domain is defined by the variable w (0<w<1) which is defined as:

(1)

Where:   

Vt,f	=	volume of fracture domain of the soil [L3], and
Vt	=	total volume of the soil [L3].

According to Equation (1) there is a different water content expression for each domain of the soil. A local one (θ) referring to the volumetric water content relative to the volume of each domain and an overall one (Θ) referring to the volumetric water content relative to the total volume of the soil. According to these concepts, the following relations apply:

(2)

(3)

(4)

(5)

(6)

(7)

(8)

Where: 

θf	=	local volumetric water content of the fracture domain [L3 L-3],
Θf	=	total water content of the fracture domain [L3 L-3],
θm	=	local water content of the matrix domain [L3 L-3],
Θm	=	total water content of the matrix domain [L3 L-3],
Θ	=	total water content of the soil [L3 L-3],
Θs,f	=	total saturated soil water content of the fracture domain [L3L-3],
Vw,f and Vw,m	=	volumes of water in the fracture and in the matrix domain respectively [L3],
Vt,m	=	volume of the matrix domain of the soil [L3], and
Vw	=	volume of water in both domains [L3] (matrix and fracture).

The degree of saturation of each domain is defined as:

(9)

(10)

(11)

Where:

Sem and Sef	=	degrees of saturation of the matrix and of the fracture domain, respectively,
Se	=	degree of saturation of the whole soil,
θs,m	=	local saturated soil water content of the matrix domain [L3 L-3],
θr,m	=	local residual soil water content of matrix the domain [L3 L-3],
Θs,m	=	total saturated soil water content of the matrix domain [L3 L-3],
Θr,m	=	total residual soil water content of the matrix domain [L3 L-3],
Θs,f	=	total saturated soil water content of the fracture domain [L3L-3],
θs,f	=	local saturated soil water content of the fracture domain [L3 L-3],
Θr	=	total residual soil water content of the soil [L3 L-3], and
Θs	=	total saturated soil water content of the soil [L3 L-3].

2.2 Model conceptualization

The soil water flow in the soil matrix domain is computed using the two-dimensional h-based Richards equation:

(12)

 Where:

z	=	vertical coordinate (positive downward) [L],
x	=	horizontal coordinate [L],
hm	=	pressure head at the matrix domain [L],
Km	=	hydraulic conductivity of the matrix domain [L/T],
Cm	=	specific water capacity (∂θm/∂hm) [L-1], and
Γw	=	water exchange term between the two domains [Τ-1].

According to Larsbo and Jarvis (2003), the water flow in the fracture domain is dominated by gravity and the capillarity is assumed to be negligible (∂h/∂z=0). Therefore, the governing water flow equation in the fracture domain is:

(13)

Where:

qf

=

flow in the fracture domain [L/T] which is described by the kinematic wave approach (Germann and Beven 1985) as:

(14)

Where:

Ks,f	=	saturated hydraulic conductivity in the fracture domain [L/T],
α (=~2)	=	a kinematic exponent (reflecting the fracture size distribution and tortuosity) (Larsbo and Jarvis 2003), and
Kf	=	unsaturated hydraulic conductivity of the fracture domain [L/T].

By replacing q_f in Equation (13) by (14), the following equation is obtained for the water flow in the fracture domain.

(15)

2.3 Hydraulic properties

In dual permeability models, the flow separation between the two domains (soil matrix and fracture) is achieved using a critical value of pressure head h_cr [L] (threshold). The flow in the fracture domain is activated when the pressure head in the soil matrix overcomes this critical value. Moreover, the local θ_cr water content, the total Θ_cr water content, and the hydraulic conductivity K_cr corresponding to h_cr are also defined.

Preferential flow, which is a dual domain problem, cannot be described using the simple form of the van Genuchten (1980) equation. A modification of the equations proposed by Vogel et al. (2000) is used in this work:

 

Where:

h	=	water pressure head [L],
αvg [L-1]	=	empirical fitting parameters (αvg >0),
nvg	=	empirical fitting parameters (nvg >1),
mvg	=	empirical fitting parameters (0<mvg<1), and
mvg	=	1-1/nvg.

The values of these fitting parameters are obtained by the water retention data for h≤h_cr. Therefore, the values of Θ_m for h>h_cr are obtained by an extrapolation of the retention curve for h≤h_cr (Figure 1).

The value of θ_s,m (and also of Θ_s,m) is a user-defined value representing the saturated conditions of the matrix domain. This approach is quite different from the approach of Larsbo and Jarvis (2003) in a MACRO model. In a MACRO model, the θ_cr and the hydraulic conductivity K_cr refer to saturated conditions, even though the degree of saturation of the matrix domain is not equal to one at this point. In Larsbo and Jarvis (2003), when the pressure head of the soil matrix exceeds the h_cr value, then an over-saturation condition takes place in the soil matrix. For this reason, a user-defined fictitious saturated soil moisture value is set (θ_s,m>θ_cr) as it is .

Figure 1 Modified van Genuchten soil water retention curve.

In Figure 1, the soil water retention curve is divided into two parts. The threshold of the two parts is the h_cr value. For pressure head values smaller than h_cr, water flow is taking place only in the matrix domain. For values of pressure head greater than h_cr, water flow is taking place in both domains (soil matrix and fracture). Water is transferred from the soil matrix to the fracture domain when h>h_cr and θ_m>θ_cr. This procedure is continued until the fracture domain is saturated (θ_f=θ_s,f). At this point, the water content in the matrix domain is allowed to exceed the value of θ_cr up to θ_s,m and the pressure head is allowed to exceed the value of h_cr up to the zero value, respectively.

The specific water capacity value, C_m is obtained by a differentiation of Equation 16.

(17)

The unsaturated hydraulic conductivity function of the soil matrix domain (K_m) is computed as in Luckner et al. (1989):

(18)

Where:

Kcr	=	hydraulic conductivity at near saturation conditions [L/T]
ℓ	=	the tortuosity factor of the soil matrix domain, and
Sem,θ cr	=	the value of degree of saturation of the soil matrix for pressure head equal to hcr (Equation 16).

According to Luckner et al. (1989) the use of K_cr (hydraulic conductivity at near saturation) instead of K_s,m (hydraulic conductivity at saturation) in Equation 18 is a preferable assumption. The near saturation value of hydraulic conductivity is easier to be determined in experimental conditions, instead of the problematic determination of the saturated hydraulic conductivity.

The hydraulic conductivity in the fracture domain (K_f) is computed as a power law expression of the degree of saturation (Se_f) (Larsbo et al. 2005):

(19)

2.4 Model conceptualization

The water transfer rate from the fracture to the matrix domain Γ_w [T^-1) is estimated as a first-order approximation of the water diffusion equation, with an assumption of rectangular-slab geometry for the aggregates (Šimůnek et al. 2003; Larsbo et al. 2005):

(20)

Where: 

d	=	an effective diffusion path-length related to soil’s aggregate size (the value of d is equal to the half of the aggregate width or half of the fracture spacing) [L],
Dw	=	water’s effective diffusivity [L2/T],
Gf	=	an aggregate’s geometry factor (for rectangular aggregate’s geometry is set equal to 3 (Gerke and van Genuchten 1996), and
γw	=	a correction factor in order to match the exact results with the approximated solutions of the diffusion problem (Gerke and van Genuchten 1993a).
Θcr	=	total soil water content at the critical value of pressure head hcr [L3L-3]

The value of Y_w depends on the initial soil moisture and the hydraulic properties of the soil. This dependence is not strong, so for simplicity, the value of Y_w is set equal to 0.4 (Gerke and van Genuchten 1993b).

Equation (20) describes the water flow from the fracture to the soil matrix domain, based on the soil water deficit of the soil matrix domain. The flow from the soil matrix to the fracture takes place in every time step of the numerical solution when the water content of the soil matrix exceeds the value of θ_cr (also when h>h_cr) following the soil physics principal that governs the water filling of pores when the water entry pressure value is exceeded (Larsbo and Jarvis 2003).

2.5 Numerical solution for the water flow equations in the matrix and the fracture domain

Richard’s equation (12) is solved in a cell-centered grid (Figure 2). The method of solution is a combination of the Alternating Direction Implicit (A.D.I.) method and the Douglas and Jones predictor-corrector method (An et al. 2011; Babajimopoulos 2000; Douglas and Rachford 1956).

Figure 2 Two-dimensional grid of Richards equation computational scheme.

The spatial subscript in the horizontal direction (x) is the letter j, and in the vertical direction (z) is the letter i. The maximum values of these two subscripts are the j_max and i_max, respectively. Finally, each time step Δt is symbolized with the superscript n (Δt=t_n+1-t_n).

The method of solution is described in the following procedure:

1^st step A.D.I. Predictor

 

1^st step A.D.I. Corrector

 

2^nd step A.D.I. Predictor

 

2^nd step A.D.I. Corrector

 

Where:

n	=	number of solving methos's time step, and
δx and δz	=	central difference operators in the x- and z- direction, respectively.

Equations 21a-d lead to tridiagonal systems of equations which are solved using the Thomas algorithm. The internodal hydraulic properties at i±1/2,j and I,j±1/2 can be computed by arithmetic, geometric, or harmonic mean.

Equation 15, which describes the one-dimensional water flow in the fracture domain, is solved at each column of the grid of Figure 2 explicitly:

 

The internodal values of saturated hydraulic conductivity K_s,f and saturated water content, θ_s,f can be computed by an arithmetic, geometric, or harmonic mean. The internodal values of soil water values θ_f (i±1/2,j) are computed using the arithmetic mean.

2.6 Boundary conditions

In single domain models, the atmospheric boundary conditions can be described by a simple equation as:

(23)

Where: 

ASF	=	actual surface flux [L/T],
r	=	applied water in soil surface (irrigation or rainfall) [L/T], and
E	=	evaporation rate from the soil surface [L/T].

In dual domain (permeability) models, the water inflow must be separated proportionally in each domain (soil matrix and fracture). In the Leak2D model, the total actual surface flux is separated according to Dusek, Gerke, and Vogel (2008) by the following equation:

(24)

Where:

ASFf  and ASFm

=

the actual surface fluxes in the fracture and in the soil matrix domain, respectively.

When r >E, the water first enters the soil matrix domain until a critical value of infiltration i_m,cr [L/T] is reached. After that, the water starts entering the fracture domain.

The critical value of infiltration i_m,cr [L/T] is defined as in Feddes et al. (1978):

(25)

Where:

Ksurface	=	hydraulic conductivity of the matrix domain at the soil surface [L/T] (=K1/2,j),
h1/2,j	=	pressure head at the surface nodal points [L], and
h1,j	=	pressure head at the first nodal points of the grid [L].

The value of h_1/2,j is set equal to h_cr when r >E and equal to h_lim when r<E. h_lim is the minimum pressure head allowed under air-dry conditions (Feddes et al. 1978):

(26)

Where:

R	=	universal gas constant [J/mole/K],
T	=	absolute temperature (K),
Μ	=	molecular weight of water [kg/mole],
g	=	gravity acceleration [m/sec2], and
f	=	relative air humidity.

In a similar way, a critical value of infiltration rate i_f,cr [L/T] is defined for the flow in the fracture domain.

(27)

Where:

Ks,f 1,j	=	saturated hydraulic conductivity [L/T],
θs,f 1,j	=	saturated water content,
θf 1,j	=	water content, at the first nodal points of the grid (i=1), and
α (=~2)	=	a kinematic exponent (reflecting the fracture size distribution and tortuosity).

According to Equation 27, the inflow in the fracture domain is defined by the following relation:

(28)

Using the critical values of infiltration defined in the equations above, inflow and outflow from the soil profile is described by the following equations:

When ASF>0, r >E, then:

(29a)

(29b)

When ASF<=0, r <=E, then:

(29c)

(29d)

Equations (29b) and (29c) describe the process, and the quantification of water entrance into the two domains (soil matrix and fracture). If water exceeds the limits of the fracture (as described in Equation 29b), then water continues to enter in the soil matrix domain and it reaches saturation at the soil surface (i.e., θ_1,j = θ_s,m). In this case, the value of h_1/2,j is changed, and is equal to 0 in the calculation of Equation (25).

Equations 29d and 29e describe the outflow of water from the soil. It is assumed that there is no evaporation from the fracture domain. The fracture domain is normally inactive, except during intense rainfall and irrigation events. Consequently, the fracture’s moisture content is at zero most of the time, and therefore cannot contribute to the evaporation process. But to avoid possible water balance errors, evaporation is calculated to the volume of the whole soil profile. This assumption is also used in dual-permeability flow models (Dusek et al. 2008).

Runoff occurs when the total moisture content at the surface nodal points reaches saturation (Θ_1,j=Θ_s). The excess water beyond that point is considered to stagnate on the soil surface and is transferred computationally to the next time step to infiltrate or evaporate after the rain–irrigation event.

2.7 Water balance computations

The results of the L2D model are evaluated through a number of criteria as proposed by Celia et al. (1990). The method proposed by these authors provides excellent results in terms of minimizing the mass balance error (Šimůnek et al. 2012).

The global mass balance error (e_w,tot) is computed by the following criterion. This criterion is important, but sometimes is not sufficient enough to guarantee good simulation results (Celia et al. 1990). The form of this criterion in a time step n+1 is:

 

Where:

CumQtop and CumQbot	=	total inflow and outflow values from the flow region [L3/T], and
CumdV	=	total water storage change both matrix and in the fracture domains [L3/T].

Another criterion for the evaluation of the simulation results is the relative mass balance error (%). This criterion is computed in every time step of the simulation and is given by:

 

Where:

Qtop and Qbot	=	inflow and outflow values from the flow region in every time step [L3/T], and
dV	=	total water storage change in the flow region in every time step [L3].

The total absolute mass balance error e_w,a (L³) in a time step n+1 is given by Šimůnek et al. (2012):

(32)

Where:

Vw n+1	=	volume of the water stored in the flow region in the time step n+1 [L3], and
Vw 0	=	volume of the water stored in the flow region in the start of the simulation procedure [L3].

2.8 Software development and interface of graphical visualization of the results of the model

The LEAK2D model was developed in a Fortran^© 2003/2008 environment. The graphical visualization of the results was developed in a Matlab^© 2021a environment.

3 Results

3.1 Flow domain and hydraulic parameters

The behaviour of the model was examined in a four-layer soil under extreme inflow conditions. Α rectangular slab, 1.0 x 0.875 m, which was discretized into 1,400 cells of 0.025 m by 0.025 m was chosen for the application. The hydraulic properties of the soil are shown in Table 1. These properties correspond to loamy soils (Schaap et al. 2001).

Table 1 Hydraulic properties of each layer.

Layer (cm)	αvg (1/m)	nvg	θs,m (cm3/cm3)	θr,m (cm3/cm3)	θs,f (cm3/cm3)	Ks,f (m/day)	W	hcr (m)	Kcr (m/day)	l
0-25	6.5	1.865	0.5	0.01	0.41	4.9	0.05	-0.012	0.15	0.5
25-50	5.5	1.853	0.51	0.01	0.42	4.91	0.06	-0.012	0.16	0.5
50-75	5.1	1.851	0.52	0.01	0.43	4.92	0.07	-0.012	0.17	0.5
75-100	5.3	1.85	0.53	0.01	0.44	4.93	0.08	-0.012	0.18	0.5

It was also considered that the value of the geometry factor G_f was 3, the effective diffusion path-length d was 45 mm, and the kinematic exponent α was equal to 2.2. An initial time step of 0.25 min was also chosen, and the duration of the simulations was 24 hours. The duration of 24 hours was chosen because at the end of this period, the flow in the fracture domain was practically zero.

3.2 Initial and boundary conditions of simulation

Slightly dry conditions (volumetric moisture content 24-26%) in most of the matrix domain were chosen for the beginning of the computations (Figure 3). In one spot in the soil slab, moisture content was set equal to 20-22%, and in three other spots, moisture content was a little higher, up to 36%. The initial moisture content of the fracture domain was considered to be zero.

A strong rainfall/irrigation event of 110 mm was considered for 0.5 hours, which was initiated 2 hours after the beginning of the simulation. This extreme value of precipitation was chosen to achieve ponding of water at the surface of the soil after the end of the precipitation. A constant evaporation of 3 mm/day was considered in all the simulations. A free drainage condition was applied at the bottom of the soil profile.

Figure 3 Initial condition of soil moisture in matrix domain.

3.3 Results of initial execution of the model

Figures 4a-d show the results of the basic execution of the model at four characteristic moments of the simulation. The first water profile (Figure 4a) is 5 minutes after the initiation of irrigation (2 hours after the beginning of the simulation). At this moment, the water starts entering the two domains of the soil. The second water profile (Figure 4b) is shown 3 hours and 45 minutes from the beginning of the simulation. At this moment, one half of the infiltration process has been completed. The third water profile (Figure 4c) is shown 5 hours and 25 minutes after the beginning of the simulation. At this moment, ponding of water at the soil surface has completed and evaporation from the soil surface is initiated. At this moment, water starts transferring from the fracture domain to the soil matrix. The last water profile shown (Figure 4d), and occurs 12 hours and 45 minutes after the beginning of the simulation. At this moment, the water has reached the bottom of the fracture domain and drainage out of the soil slab starts.

In all the figures, the local soil moisture is presented in both domains, and the total soil moisture content in the whole soil profile is presented as well.

Figure 4a Water profile at the matrix and fracture domain, and total moisture of the whole soil slab 2 hrs and 5 minutes after the start of simulation.

Figure 4b Water profile at the matrix and fracture domain, and total moisture of the whole soil slab in the middle of the irrigation event.

Figure 4c Water profile at the matrix and fracture domain, and total moisture of the whole soil slab 5 hrs and 25 minutes after the start of simulation (end of ponding).

Figure 4d Water profile at the matrix and fracture domain, and total moisture of the whole soil slab 12 hr and 45 min after the start of simulation.

Figures 4a-d show that the model responds very satisfactorily to an extreme irrigation/precipitation event. The role of the fracture domain is obvious in these figures. It is also shown that the initial distribution of the water in the soil slab plays an important role in the water profiles in all the domains in the subsequent moments.

It is also pointed out that the total relative water mass balance error (%) of this simulation was 3.55%, and the incremental water mass balance error was -0.0022178 (m³) which are considered very satisfactory in relation to the extreme irrigation/precipitation event which was simulated.

3.4 Sensitivity analysis of parameters affecting flow in the fracture domain

Modeling non-equilibrium water flow processes in soils is an inherently difficult task due to the unknown nature of the flow pathways, their large spatial and temporal variability, and our partial lack of knowledge of the transport processes (Jarvis and Larsbo 2012).

Flow in the fracture domain depends on the value of the kinematic exponent, α (Equations 14 and 15) which reflects the fracture size distribution and tortuosity (Larsbo and Jarvis 2003). The interaction of flow between the fracture and the matrix domain depends on the effective diffusion path length, d (Equation 20) which is related to the aggregate size and the influence of coatings on macropore and aggregate surfaces (Larsbo and Jarvis 2003), and also on the geometry factor G_f. According to Larsbo et al. (2005) the value of α varies between 2 and 8, the value of G_f. for a rectangular slab geometry is equal to 3, and 15 for spheres (Ray et al. 1997). The value of d varies between 1 and 50 mm (Jarvis and Larsbo 2012).

Although the above three parameters have a physical meaning, their exact values are difficult to be determined. Also, not all parameters are measurable, they vary spatially and also temporally (which technically means that they are not true parameters), and it is impossible to fully characterize a field site by measurements alone, especially in the subsoil (Jarvis and Larsbo 2012). Jarvis and Larsbo (2012) also recommend that these parameters be included in the calibration process. Therefore, a sensitivity analysis of these parameters is very important to define their effect on the output of the model. The results of this analysis are presented in the next paragraphs.

The initial values of these parameters found in the literature and their variation is shown in Table 2. Each parameter varied independently, keeping the other two parameters constant on their table values.

Table 2 Model’s input factors for sensitivity analysis.

Input factor	Initial value	Tested values
α	2.2	±10%, ±25%, ±50%
Gf	3.0	±10%, ±25%, ±50%
d (mm)	25	±10%, ±25%, ±50%

An extreme rainfall event of 110 mm for half an hour was chosen in all the simulations. Simulation results are presented in three instances: a) 3 hr and 45 min after the start of simulation, b) at the end of ponding (5 hr and 25 min after the start of simulation), and c) 7 hr and 5 min after the start of simulation, when the differences among the various simulations are more distinguishable.

In all the simulations, only the water profile in the fracture domain is presented since the major changes in the water content take place in this domain.

Sensitivity analysis of the kinematic exponent, α

Figures 5a–5c show the effect of the variation of α according to the values shown in Table 2. Water profiles of the fracture domain are shown at the middle of the irrigation event (3 hr and 45 min from the beginning of the simulation), at the end of the ponding (5 hr and 25 min), and 7 hr and 5 min from the beginning of the simulation.

Figure 5a Sensitivity analysis of kinematic exponent α. Water profiles of the fracture domain in the middle of the irrigation/rainfall event.

Figure 5b Sensitivity analysis of kinematic exponent α. Water profiles of the fracture domain at the end of ponding.

Figure 5c Sensitivity analysis of kinematic exponent α. Water profiles of the fracture domain 7 hr and 5 min after the start of simulation.

It is concluded that the kinematic exponent, α plays an important role in the water distribution in the fracture domain. The smaller the value of α, the more water enters the fracture domain. This is more obvious when comparing the water profiles at the end of ponding (Figure 5b), and at 7 hr and 5 min after the simulation start (Figure 5c). The moisture content of the lower nodal points of the fracture domain is especially higher when α takes its smallest value. However, higher water content leads to higher deep percolation, and higher water exchange with the matrix domain. Therefore, after all the rain/irrigation water has entered the soil profile, the profiles with the smaller values of α, and therefore with the higher water content, start losing more water.

Figure 6 shows the volume of water which is stored in the fracture domain up to the end of the simulation (24 hrs). It is shown that when α has its smallest value, the water content in the fracture domain is higher up to sometime after ponding. After that time, water content starts reducing at a higher rate in accordance with the smaller value of α. Table 3 shows the water that is stored in the fracture domain, and the water content at the end of the 24-hr simulation. It is shown that the smaller the value of α, the higher the water content of the domain. For example, when the value of α is reduced by -50%, the maximum value of water in the fracture domain is 1.26E-02 m³, which is 29.1% higher than the water stored in this domain when α has its literature value. When the value of α is at its largest value, then the maximum value of water in the domain is only 7.91E-03 m³, which is 19.1% less than the water stored when α has its literature value. On the contrary, the water remaining in the profile after the end of the 24-hr simulation is 0 when α is at its smallest value, and 1.72E-03 m³ (about 157% more than the water remaining in the domain with the literature value of α) when α is at its largest value.

Figure 6 Water volume in fracture domain (sensitivity analysis of kinematic exponent).

Table 3 Maximum volume of water stored in the fracture domain, and volume of water stored in the fracture domain, at the end of the 24-hr simulation.

α	-50%	-25%	-10%	0%	+10%	+25%	+50%
Maximum volume of water stored in the fracture domain (m3)	1.26E- 02	1.02E- 02	1.00E- 02	9.76E- 03	9.52E- 03	8.91E- 03	7.91E- 03
Difference (%)	29.1%	4.5%	2.5%	-	-2.5%	-8.7%	-19.1%
Volume of water stored in the fracture domain at the end of the 24-hr simulation (m3)	0.00E+ 00	4.93E- 05	2.87E- 04	6.68E- 04	1.19E- 03	1.70E- 03	1.72E- 03

It was shown that the larger the value of α, the more water enters the fracture domain, leading to a higher contribution of the fracture domain to deep percolation of water. This is shown in Figure 7 where the total cumulative volume of water leaving the soil profile from the bottom nodal points of the matrix, and the fracture domain is presented. Since this volume is leaving from the system, it is presented with negative values.

It is very clear that the larger the value of α, the more water is percolating out of the soil slab.

Figure 7 Cumulative volume of water (m³) percolating out of the soil slab in relation to the value of the kinematic exponent, α.

Sensitivity analysis of the geometry factor G_f

The effect of the variation of the geometry factor G_f to the water profiles of the fracture domain are shown in Figures 8a-c. The value of G_f is varied according to the values shown in Table 4. As it happened with the kinematic exponent α, water profiles of the fracture domain are shown at the middle of the irrigation event (3 hr and 45 min from the beginning of the simulation), at the end of the ponding (5 hr and 25 min), and 7 hr and 5 min from the beginning of the simulation.

Figure 8a Sensitivity analysis of geometry factor G_f. Water profiles of the fracture domain in the middle of the irrigation/rainfall event.

Figure 8b Sensitivity analysis of geometry factor G_f. Water profiles of the fracture domain at the end of ponding.

Figure 8c Sensitivity analysis of geometry factor G_f. Water profiles of the fracture domain 7 hrs and 5 mins after the start of simulation.

The effect of G_f on the water profiles in the fracture domain is like that of the kinematic exponent α, but to a lesser extent. As with α, the smaller the value of G_f, the more water is stored in the fracture domain. This is more distinguishable by observing the profiles at the end of ponding (Figures 8b) and at 7 hr and 5 min (Figure 8c) after the start of the simulation.

In a similar way to that of α, Figure 9 shows the volume of water which is stored in the fracture domain up to the end of the 24-hr simulation. The value of G_f affects the volume of water that is stored in the fracture domain, but to a lesser degree than that of α. Table 4 shows the values of the maximum volume of water stored in the fracture domain in relation to the values of G_f. When G_f is reduced by -50%, the maximum volume of water stored in the fracture domain is 1.04E-02 m³, whereby it is 9.08E-03 m³ when G_f is increased by +50% in relation to the literature value. These values are 6.6% more and 7.0% less, respectively, than the maximum volume of water stored in the fracture domain when G_f takes its literature value. These differences are smaller than those observed with the kinematic exponent α. Regardless of the value of G_f, the water remaining in the fracture domain is almost the same, to about 7E-04 m³.

Figure 9 Water volume in fracture domain (sensitivity analysis of geometry factor).

Table 4 Maximum volume of water stored in the fracture domain, and volume of water stored in the fracture domain, at the end of the 24-hr simulation.

Gf	-50%	-25%	-10%	0%	+10%	+25%	+50%
Maximum volume of water stored in the fracture domain (m3)	1.04E- 02	9.94E- 03	9.81E- 03	9.76E- 03	9.68E- 03	9.45E- 03	9.08E- 03
Difference (%)	6.6%	1.8%	0.5%	-	-0.8%	-3.2%	-7.0%
Volume of water stored in the fracture domain at the end of the 24-hr simulation (m3)	7.31E- 04	7.22E- 04	6.91E- 04	6.68E- 04	6.71E- 04	7.08E- 04	7.00E- 04

As observed with the kinematic exponent α, an increased volume of water in the fracture domain means higher contribution of this domain to deep percolation. Figure 10 shows the cumulative volume of water (m³) which is percolating out of the soil slab from each of the bottom nodal points, from both the domains (matrix and fracture). As it was shown with the kinematic exponent α, the smaller the value of G_f, the more water is removed from the soil slab.

Figure 10 Cumulative volume of water (m³) percolating out of the soil slab in relation to the value of the geometry factor, G_f.

Effective diffusion path-length d

The sensitivity analysis of the effective diffusion path-length d is presented in the following figures (Figures 11a-c):

Figure 11a Sensitivity analysis of effective diffusion path-length d. Fracture’s domain water profiles in the middle of irrigation event.

Figure 11b Sensitivity analysis of effective diffusion path-length d. Fracture’s domain water profiles at the end of the ponding.

Figure 11c Sensitivity analysis of effective diffusion path-length d. Fracture’s domain water profiles at 7 hr and 5 min after the simulation start.

The effect of the effective diffusion path-length, d to the water profiles in the fracture domain is opposite to what happens with α and G_f. The smaller the value of d, the slower the movement of water in the fracture domain, as can be seen in Figures 11a-c. The higher value of d implies faster movement of the water in the domain, and greater percolation out of the soil slab.

Figure 12 shows the volume of water stored in the fracture domain up to the end of the 24-hr simulation. It shows that the effective diffusion path length d has an opposite effect to those of α and G_f. The larger the value of d, the more water is stored in the fracture domain. Table 5 shows the values of the maximum volume of water stored in the fracture domain in relation to the values of d. When d is increased by 50%, the maximum volume of water in the domain is 1.06E-02 m³, which is 8.6% more than the volume stored with the literature value of d. The corresponding value of the volume of water when d is decreased by -50% is 7.96E-03 m³, which is 18.4% less than the value obtained with the literature value of d. As happened with G_f, in all the simulations, the volume of water at the end of the 24-hr period is almost 7.E-04 m³.

Figure 12 Water volume in fracture domain (sensitivity analysis of d).

Table 5 Maximum volume of water stored in the fracture domain and volume of water stored in the fracture domain at the end of the 24-hr simulation.

d	-50%	-25%	-10%	0%	+10%	+25%	+50%
Maximum volume of water stored in the fracture domain (m3)	7.96E- 03	8.81E- 03	9.47E- 03	9.76E- 03	9.86E- 03	1.01E- 02	1.06E- 02
Difference (%)	-18.4%	-9.7%	-3.0%	-	1.0%	3.5%	8.6%
Volume of water stored in the fracture domain at the end of the 24-hr simulation (m3)	6.84E- 04	6.62E- 04	7.04E- 04	6.68E- 04	7.07E- 04	7.49E- 04	7.43E- 04

In accordance with the results relating to d, Figure 13 shows the cumulative volume of water leaving from both the domains (matrix and fracture) from the bottom nodal points. As was expected, this figure differs from the corresponding Figures 6 and 9 relating to α and G_f.

The smaller the value of d, the lower the amount of water leaving the soil slab to deep percolation.

Figure 13 Cumulative volume of water (m³) percolating out of the soil slab in relation to the value of the effective diffusion path length, d.

4 Conclusions and discussion

The simulation of preferential flow in soils is a complex and complicated process, mainly for three reasons.

1. The velocity by which the water moves in the fracture domain is much faster than the velocity in the soil matrix domain. The validation of the Darcy equation in the fracture domain is questionable, and therefore special effort is required to simulate the flow in this domain;
2. The difficulty simulating the exchange of water between the two domains. This process depends on several parameters which, although they have physical meaning, are difficult to measure; and
3. The difficulty measuring and defining the hydraulic properties in the fracture domain.

The dual-permeability methodology, which allows the movement of water in both domains—matrix and fracture—has often been used in the literature and has proven very promising in describing the preferential flow of water in the vadose zone.

In this paper, a new model is presented to simulate preferential flow in the soil. The two-dimensional nonlinear Richards equation is used to describe the flow in the soil matrix domain. In the fracture domain, it is assumed that capillary forces are absent. Therefore, flow is described here by the kinematic wave equation, which does not require the definition of water retention parameters.

An analytical solution of the Richards equation is excluded because this equation is nonlinear.

The Richards equation is solved numerically, usually following the finite difference or the finite element methodology. Following the first methodology, the most common method of solution is the implicit method, which requires the solution of very big systems of equations, or several iterations in each time-step up to convergence. Both procedures are time consuming, and because of the huge number of computations, they cannot avoid the existence of many round-off errors. The method proposed in this paper is implicit but requires only two iterations in each line and column of the grid. It is based on the alternating direction implicit (A.D.I.) method combined with the Douglas-Jones predictor-corrector method leading to tri-diagonal systems of equations. These systems are solved by the Thomas algorithm, which is very efficient, and therefore avoids many round-off errors. The kinematic wave equation, which describes flow in the fracture domain, is one-dimensional and is solved very easily by an explicit method.

The model is developed in Fortran^© 2003/2008. Also, the graphical interface of Matlab^© for post processing of the data is used as a powerful tool for precise and dynamic presentation of the results of the model.

A sensitivity analysis of the parameters involved in the water exchange between the fracture and the matrix domain was performed. The kinematic exponent α, the geometry factor G_f, and the diffusion path length d, were varied by ±10%, ±25%, and ±50% of the most common values provided in the literature.

It was found that the kinematic exponent α, has the strongest influence on the exchange process. When the value of α was reduced, then more water entered the fracture domain. Opposite results, but to a lesser extent, were obtained when α values were larger than the literature value. For example, when the value of α was reduced to 50% less than the literature value, the maximum volume of water entering the fracture domain during the whole simulation period was 29.1% more than the maximum volume of water entering the fracture domain when α was equivalent to the literature value (standard volume) (Table 3). When α was increased by 50%, the maximum volume of water entering the fracture domain was reduced by 19.1% relative to the standard volume. In agreement with these results, more water is leaving the soil slab to deep percolation when α has a smaller value, and less water is percolating when α has a larger value (Figure 7).

The geometry factor G_f, has a similar effect on the exchange of water between the fracture and the matrix domain as α, but to a smaller degree, as it can be seen from Figures 8a-c and 5a-c. When G_f is reduced by 50% of its literature value, then the maximum volume of water entering the fracture domain is decreased by 6.6% of the standard volume (Table 4). In contrast, when G_f is increased by 50%, the maximum volume of water entering the fracture domain is decreased by 7%. Similarly, the water leaving the soil slab to deep percolation is greater when G_f is less than the literature values, and less when G_f has greater values (Figure 10).

In contrast to α and G_f, the diffusion path length d has an opposite effect on the exchange of water process between the fracture and the matrix domain. When the value of d is greater than the literature values, this results in more water entering the fracture domain, and smaller values result in less water entering the fracture domain. When d is reduced by 50%, then the maximum volume of water entering the fracture domain is 18.4% less than the standard value (Table 5). When the value of d is increased by 50%, then the maximum volume of water entering the fracture domain is increased by 8.6% compared to the standard value. In accordance with these results, the water leaving the soil slab to deep percolation is greater when d is higher than its literature value, and it is lower when d is less than the literature value (Figure 13). The unique approach of the L2D model to modeling preferential flow in both the matrix and fracture domain makes it a valuable tool for managing stormwater and groundwater in civil engineering applications. In urban areas, effective stormwater management is crucial to preventing flooding and environmental damage. The L2D model can assist in designing stormwater infrastructure by providing detailed information on water flow patterns and hydraulic properties, enabling the optimization of stormwater management systems for more efficient and effective water management.

Furthermore, in civil engineering applications such as groundwater modeling and contaminated site remediation, the L2D model's ability to accurately model preferential flow can aid in understanding and predicting water movement through complex subsurface environments. This information can be used to design more effective remediation strategies and optimize water resource utilization.

References

1. An, H., Y. Ichikawa, Y. Tachikawa, and M. Shiiba. 2011. “A New Iterative Alternating Direction Implicit (IADI) Algorithm for Multi-Dimensional Saturated–Unsaturated Flow.” Journal of Hydrology 408 (1–2): 127–39.
2. Babajimopoulos, C. 2000. “Revisiting the Douglas–Jones Method for Modelling Unsaturated Flow in a Cultivated Soil.” Environmental Modelling & Software 15 (3): 303–12.
3. Celia, M., E. Bouloutas, and R. Zarba. 1990. “A General Mass-Conservative Numerical Solution for the Unsaturated Flow Equation.” Water Resources Research 26 (7): 1483–96.
4. Douglas, J., and B. Jones. 1963. “One Predictor-Corrector Method for Non-Linear Parabolic Differential Equations.” Journal Society for Industrial and Applied Mathematics 11: 195–204.
5. Douglas, J., and H. Rachford. 1956. “On the Numerical Solution of Heat Conduction Problems in Two and Three Space Variables.” Transactions of the American Mathematical Society 82: 421. https://www.ams.org/journals/tran/1956-082-02/S0002-9947-1956-0084194-4/
6. Dusek, J., H. Gerke, and T. Vogel. 2008. “Surface Boundary Conditions in Two-Dimensional Dual-Permeability Modeling of Tile Drain Bromide Leaching.” Vadose Zone Journal 7 (4): 1287–1301.
7. Feddes, R., P. Kowalik, and H. Zaradny. 1978. Simulation of Field Water Use and Crop Yield. Simulation Monographs, PUDOC, Wageningen, Netherlands. Wageningen: Centre for Agricultural Publishing and Documentation, Simulation Monographs.
8. Gerke, H. 2006. “Preferential Flow Descriptions for Structured Soils.” Journal of Plant Nutrition and Soil Science 169 (3): 382–400.
9. Gerke, H., and M. van Genuchten. 1993a. “A Dual-Porosity Model for Simulating the Preferential Movement of Water and Solutes in Structured Porous Media.” Water Resources Research 29 (2): 305–19.
10. Gerke, H., and M.Th. van Genuchten. 1993b. “Evaluation of a First-Order Water Transfer Term for Variably Saturated Dual-Porosity Flow Models.” Water Resources Research 29 (4): 1225–38.
11. Gerke, H., and M.Th. van Genuchten. 1996. “Macroscopic Representation of Structural Geometry for Simulating Water and Solute Movement in Dual-Porosity Media.” Advances in Water Resources 19 (6): 343–57.
12. Germann, P., and K. Beven. 1985. “Kinematic Wave Approximation to Infiltration into Soils with Sorbing Macropores.” Water Resources Research 21 (7): 990–96.
13. Guo, L., and H. Lin. 2018. Addressing Two Bottlenecks to Advance the Understanding of Preferential Flow in Soils. Advances in Agronomy. 1st ed. Vol. 147. Elsevier Inc. https://doi.org/10.1016/bs.agron.2017.10.002
14. Hendrickx, J., and M. Flury. 2001. “Uniform and Preferential Flow Mechanisms in the Vadose Zone.” In Conceptual Models of Flow and Transport in the Fractured Vadose Zone, edited by Paul Hsieh, 149–87. Washington, D.C.: National Academy Press.
15. Jarvis, N., and M. Larsbo. 2012. “MACRO (v5.2): Model Use, Calibration, and Validation.” Transactions of ASABE 55 (4): 1413–23.
16. Larsbo, M., and N. Jarvis. 2003. MACRO 5.0. A Model of Water Flow and Solute Transport in Macroporous Soil. Technical Description. Department of Soil Sciences, Swedish University of Agricultural Sciences. Uppsala, Sweden.
17. Larsbo, M., S. Roulier, F. Stenemo, R. Kasteel, and N. Jarvis. 2005. “An Improved Dual-Permeability Model of Water Flow and Solute Transport in the Vadose Zone.” Vadose Zone Journal 4 (2): 398–406.
18. Luckner, L., M.Th. Van Genuchten, and D. Nielsen. 1989. “A Consistent Set of Parametric Models for the Two-Phase Flow of Immiscible Fluids in the Subsurface.” Water Resources Research 25 (10): 2187–93.
19. Nimmo, J.R. 2021. “The Processes of Preferential Flow in the Unsaturated Zone.” Soil Science Society of America Journal 85 (1): 1–27. https://doi.org/10.1002/saj2.20143
20. Pluer, W.T., M. Macrae, A. Buckley, and K. Reid. 2020. “Contribution of Preferential Flow to Tile Drainage Varies Spatially and Temporally.” Vadose Zone Journal 19 (1): 1–15. https://doi.org/10.1002/vzj2.20043
21. Ray, C., T. Ellsworth, A. Valocchi, and C. Boast. 1997. “An Improved Dual Porosity Model for Chemical Transport in Macroporous Soils.” Journal of Hydrology 193 (June): 270–92. https://doi.org/10.1016/S0022-1694(96)03141-1
22. Schaap, M., F. Leij, and M.Th. van Genuchten. 2001. “ROSETTA: A Computer Program for Estimating Soil Hydraulic Parameters with Hierarchical Pedotransfer Functions.” Journal of Hydrology 251 (October): 163–76.
23. Šimůnek, J., M.Th. van Genuchten, and M. Šejna. 2012. The HYDRUS Software Package for Simulating Two- and Three-Dimensional Movement of Water, Heat, and Multiple Solutes in Variably-Saturated Media, Technical Manual, Version 2.0. PC Progress. Prague, Czech Republic.
24. Šimůnek, J., N. Jarvis, M.Th. Van Genuchten, and A. Gärdenäs. 2003. “Review and Comparison of Models for Describing Non-Equilibrium and Preferential Flow and Transport in the Vadose Zone.” Journal of Hydrology 272 (1–4): 14–35.
25. van Genuchten, M.Th. 1980. “A Closed-Form Equation for Predicting the Hydraulic Conductivity of Unsaturated Soils.” Soil Science Society of America Journal 44 (5): 892–98. https://doi.org/10.2136/sssaj1980.03615995004400050002x
26. van Genuchten, M.Th., and P. Wierenga. 1976. “Mass Transfer Studies in Sorbing Porous Media I. Analytical Solutions.” Soil Science Society of America Journal 40 (4): 473.
27. Vogel, T., M.Th. van Genuchten, and M. Cislerova. 2000. “Effect of the Shape of the Soil Hydraulic Functions near Saturation on Variably-Saturated Flow Predictions.” Advances in Water Resources 24 (2): 133–44.