Simulating and Comparison of Non-Equilibrium Water and Potassium Transport in Sandy Loam Soil Using HYDRUS-1D Software Package

INTRODUCTION

The problem of non-equilibrium and/or preferential flow and transport has received much attention in the soil and agricultural sciences because of its implications in accelerating the movement of agricultural contaminants (fertilizers, pesticides,pathogenic microorganisms and toxic trace elements) through the unsaturated zone to underlying groundwater.

Non-equilibrium flow and transport are probably the most frustrating processes in terms of restrain accurate predictions of contaminant transport in soils and fractured rocks. Non-equilibrium or/and preferential flow, as opposed to uniform flow, results in irregular wetting of the soil profile as a direct consequence of water moving faster in certain parts of the soil profile than in others.

Uniform flow in granular soils and preferential flow in structured media (both macroporous soils andfractured rocks) can be described using a variety of single-porosity, dual-porosity, dual-permeability,multi-porosity, and multi-permeability models (Richards, 1931; Pruess and Wang, 1987; Gerke and vanGenuchten, 1993a; Gwoet al., 1995; Jarvis, 1998; Šimůneket al., 2003,2005, 2008). While single-porositymodels assume that a single pore system exists that is fully accessible to both water and solute, dualporosityand dual-permeability models both assume that the porous medium consists of two interactingpore regions, one associated with the inter-aggregate, macropore, or fracture system, and one comprisingthe micropores (or intra-aggregate pores) inside soil aggregates or the rock matrix.Whereas dual-porositymodels assume that water in the matrix is stagnant, dual-permeability models allow also for water flowwithin the soil or rock matrix (Jarvis, 1998; Šimůneket al., 2003, 2005 and 2008).

The dual-porosity modelassumes the presence of two pore regions, with water in one region being immobileand in the other region mobile. This model allows the exchange of both water and solute between thetwo regions (Šimůneket al., 2003 and 2008). Conceptually, this formulation views the soil as consisting of a soilmatrix containing grains/aggregates with a certain internal microporosity (intra-aggregate porosity) anda macropore or fracture domain containing the larger pores (inter-aggregate porosity). Dual- permeability models are those in which water can move in both the inter-and intra-aggregate pore regions (and matrix and fracture domains). These models in various forms are now also becoming increasingly popular (Pruess and Wang, 1987; Gerke and van Genuchten, 1993a; Jarvis, 1994; Pruess, 2004).

The HYDRUS software packages (Šimůneket al., 1998, 2005, 2007, and 2008) are among the most widely used models simulating water flow and solute transport in soils.Traditionally, the HYDRUS software, as well as many other models simulating variably saturated water flow and solute transport, either did not consider non-equilibrium flow and transport processes at all or only considered them separately. For example, previous versions of HYDRUS codes considered physical and chemical non-equilibrium separately. Physical non-equilibrium solute transport was accounted for by assuming a two-region, dual-porosity type formulation that partitions the liquid phase into mobile and immobile regions (van Genuchten andWierenga 1976). Chemical non-equilibrium solute transport was accounted for by assuming a two-site sorption model, which assumes that sorption sites can be divided into two fractions with sorption to different fractions of sorption sites being either instantaneous or kinetic (van Genuchten andWagenet, 1989).

There is increasing evidence that flow and transport processes in soils often cannot be described using classical models that assume uniform flow and transport (Nkedi- Kizzaet al., 1984; Hendrickx and Flury, 2001; Pot et al., 2005; Kohneet al., 2006). Many laboratory and field experiments have demonstrated the presence of non-equilibrium flow and transport conditions in soils. Non-equilibrium water flow and solute transport in the unsaturated zone can be simulated at present by means of a large number of models of various degrees of complexity and dimensionality. Modeling approaches range from relatively simple analytical solutions for solute transport (van Genuchten,1981; Torideet al., 1993) to complex numerical codes (Šimůneket al., 2005 and 2008; Jacques and Šimůnek, 2005).

The objectives of the present work were to evaluate and comparison the non-equilibrium processes importance in potassium transport under the present experimental conditions.

Theoretical basis

Equilibrium and Non-equilibrium flow and transport

Various Equilibrium and Non-equilibrium Solute Transport Models can be used. These models include (Šimůneket al., 2008):

Equilibrium model (EQ)

Numerical models for water flow in soils are usually based on the following equation:

or its extensions (e.g., for two- and three-dimensional systems). In Eq. (15), often referred to as the Richards equation, z is the vertical coordinate positive upward [L], tis time [T], h is the pressure head [L], θ is the water content [L³ L⁻³], S is a sink term representing root water uptake or some other source or sink [T⁻¹], and K(h) is the unsaturated hydraulic conductivity function, often given as the product of the relative hydraulic conductivity, K_r (dimensionless), and the saturated hydraulic conductivity, K_s [L T⁻¹].

Dual-Porosity Model (DPM)

Dual-porosity models assume that water flow is restricted to the macropores (or interaggregate pores and fractures), and that water in the matrix (intraaggregate pores or the rock matrix) does not move at all. This conceptualization leads to two-region type flow and transport models (van Genuchten and Wierenga, 1976) that partition the liquid phase into mobile (flowing, interaggregate), θ_mo, and immobile (stagnant, intraaggregate), θ_im, regions [L³ L⁻³]:

The dual-porosity formulation for water flow can be based on a mixed formulation of the Richards Eq. (15) to describe water flow in the macropores (the preferential flow pathways) and a mass balance equation to describe moisture dynamics in the matrix as follows (Šimůneket al., 2003):

Where:  S_mo and S_im are sink terms for the mobile and immobile regions, respectively [T⁻¹], and Γ_w is the transfer rate for water between the inter- and intra-aggregate pore domains [T⁻¹]. Šimůneket al. (2003) and Kohneet al. (2004) discussed different formulations that can be used to evaluate the mass transfer rate Γ_w.

Dual-Permeability Model (DP)

Different dual-permeability approaches may be used to describe flow and transport in structured media. While some models invoke similar equations for flow in the fracture and matrix regions, others use different formulations for the two regions. A typical example of the first approach, implemented in HYDRUS-1D, is the work of Gerke and van Genuchten (1993a,b and 1996), who applied the Richards equation to each of the two pore regions. The flow equations for the macropore (fracture), subscript f and matrix (subscript m) pore systems in their approach are given by:

Where: w is the ratio of the volumes of the macropore or fracture domain and the total soil system (dimensionless). This approach is relatively complicated in that the model requires characterization of water retention and hydraulic conductivity functions (potentially of different form) for both pore regions, as well as a hydraulic conductivity function of the fracture–matrix interface. Note that the water contents, θ_f and θ_m in Eq. (19), have different meanings than in Eq. (16), where they represented water contents of the total pore space (i.e., θ = θ_mo + θ_im), while here they refer to water contents of the two separate (fracture or matrix) pore domains such that:

Hence, lowercase subscripts in the dual-permeability model refer to the local (pore-region) scale, while uppercase subscripts refer to the global (total soil medium) scale.

Mobile-Immobile Water and Dual-Porosity Models (physical non-equilibrium), MIMDP

The concept of two-region, dual-porosity type solute transport was implemented already in earlier versions (1.0 and 2.0) of HYDRUS-1D to permit consideration of physical non-equilibrium transport. While the physical non-equilibrium transport model in the earlier versions was combined only with uniform water flow Eq. (12), Version 4.17 of HYDRUS-1D (Šimůneket al., 2013) was expanded to also consider the dual-porosity water flow model Eq. (20) with transient immobile water content. In both implementations, the governing solute transport equations are as follows:

in which solute exchange between the two liquid regions is modeled as the sum of an apparent first-order diffusion process and advective transport (where applicable). In Eq. (20), c_mo and c_im are concentrations of the mobile and immobile regions [M L⁻³], respectively; s_mo and s_im are sorbed concentrations of the mobile and immobile regions [M M⁻¹], respectively; D_mo is the dispersion coefficient in the mobile region [L² T⁻¹], q_mo is the volumetric fluid flux density in the mobile region [L T⁻¹], φ_mo and φ_im are sink–source terms that account for various zero- andfirst-order or other reactions in both regions [M L⁻³ T⁻¹]; f_mo is the fraction of sorption sites in contact with the mobile water content (dimensionless), ω_mim is the mass transfer coefficient [T⁻¹], and Γ_s is the mass transfer term for solutes between the mobile and immobile regions [M L⁻³ T⁻¹]. Equation (20a) describes solute transport in the mobile (macropore) zone, Eq. (20b) is a mass balance for the immobile (micropore) domain, while Eq. (20c) (Γs) describes the rate of mass transfer between the mobile and immobile domains. The second (advective) term of Γ_s in Eq. (20) is equal to zero for the Mobile–Immobile Model since that model does not consider water flow between the two regions. In the Dual-Porosity Model, c* is equal to c_mo for Γ_w> 0 and c_im for Γ_w< 0.

Dual-Permeability Model with Immobile Water (physical and chemical non-equilibrium), DPIMM

The Dual-Permeability Model with Immobile Water assumes that the liquid phase of the matrix can be further partitioned into mobile (flowing), θ_m,m[L³ L⁻³], and immobile (stagnant), θ_im,m[L³ L⁻³], regions as follows:

Where: θ_m is the volumetric water content of the matrix pore system [L³ L⁻³]. The governing advection–dispersion equation for transport in the matrix region (Eq. (22b)) is then replaced with the modified equations (Eq. (22) (e.g., Pot et al., 2005) to yield

where c_im,mand c_m,mare solute concentrations in the immobile and mobile zones of the matrix region [M L⁻³], respectively; φ_m,mand φ_im,mrepresent various reactions in the mobile and immobile parts of the matrix [M L⁻³ T⁻¹], respectively; f_m is again the fraction of sorption sites in contact with the mobile region of the matrix (dimensionless), ω_dpm is the mass transfer coefficient between mobile and immobile zones of the matrix region [T⁻¹], and Γ_s* is the mass transfer term for solutes between the mobile and immobile regions of the matrix domain [M L⁻³ T⁻¹]. Equation (22a) now describes solute transport in the fracture domain, Eq. [14b] transport in the mobile zone of the matrix domain, Eq. (22c) is a mass balance for the immobile zone of the matrix domain, Eq. (22d) describes mass transfer between the fracture and matrix domains, while Eq. (22e) describes mass transfer between the mobile and immobile zones within the matrix domain.

One Site Sorption Model (Chemical non-equilibrium), OSM

When sorption in the Uniform Transport Model is considered to be kinetic, Eq. (23) needs to be supplemented with an equation describing the kinetics of the sorption process. This is usually done by assuming a first-order process as follows:

where s_e^k is the sorbed concentration that would be reached at equilibrium with the liquid-phase concentration [M M⁻¹], s^k is the sorbed concentration of the kinetic sorption sites [M M⁻¹],α_k is a first-order rate constant describing the kinetics of the sorption process [T⁻¹], and φ_krepresents a sink–source term that accounts for various zero- and first-order or other reactions at the kinetic sorption sites [M L⁻³ T⁻¹].

Two-Site sorption Model (chemical non-equilibrium), TSM

Similarly to the mobile–immobile water concept, the concept of two-site sorption (Selimet al., 1976; van Genuchten and Wagenet, 1989) was implemented already in Versions 4.17 of HYDRUS-1D (Šimůneket al., 2013)to permit consideration of non-equilibrium adsorption–desorption reactions. The two-site sorption concept assumes that the sorption sites can be divided into two fractions:

Sorption s^e [M M⁻¹], on one fraction of the sites (Type 1 sites) is assumed to be instantaneous, while sorption s^k [M M⁻¹], on the remaining (Type 2) sites is considered to be a first-order kinetic rate process. The system of equations describing the Two-Site Model is given by:

Where:  f_e is the fraction of exchange sites assumed to be in equilibrium with the liquid phase (dimensionless), and α_k is a first-order rate constant [T⁻¹]. Equation (25a) describes solute transport in the total system, Eq. (25b) equilibrium sorption onto the instantaneous sorption sites, Eq. (25c) is a mass balance of the kinetic sorption sites (van Genuchten and Wagenet, 1989), while Eq. (25d) represents the sorbed concentration of the kinetic sites when equilibrium would be reached with the liquid-phase concentration.

MATERIALS AND METHODS

Soil

The soil used in the present experiment was collected from surface layer (0-30 cm depth) of El-Hammam region, Matrouh Governorate. The soil was air dried and passed through 2.0 mm sieve. Some physical and chemical properties of the soilsample are reported in Table (1). The soil properties were performed according to the procedures outlined in Carter and Gregorich (2008).

Table (1). Some physical and chemical characteristics of the experimental soil

Soil hydraulic parameters

Soil water retention, q(h) and soil hydraulic conductivity, K(h) functions are specified. The soil water retention curve was determined using the pressure plate extractor according to the method described in Klute (1986). The water retention data was described using the van Genuchten model (van Genuchten, 1980):

Where:

S_e    is effective saturation(-);

q(h) is the soil water content at matric potential, h (L³L^-3);

h      is the soil matric potential (L);

q_r    is the residual volumetric water content (L³L^-3);

q_s    is the saturated volumetric water content (L³L^-3);

n and m(m=1-1/n) are empirical coefficients; and

ais a fitting parameter (L^-1)

The hydraulic conductivity function was described using the capillary model (Mualem, 1976) as applied to the van Genuchten model (van Genuchten, 1980):

Where:

K(h)  is the unsaturated hydraulic conductivity(LT^-1);

K_s     is the saturated hydraulic conductivity (LT^-1); and

i                is the pore connectivity coefficient (assumed as 0.5, Mualem, 1976)

The saturated hydraulic conductivity (K_s) was determined using constant head permeameter method as described in Klute (1986).

The soil hydraulic parameters were fitted using the RETC model (van Genuchten et al., 1991).

Durner(1994) divided the porous medium into two (or more) overlapping regions and suggested to use for each of these regions a van Genuchten-Mualem type function (van Genuchten, 1980) of the soil hydraulic properties. Linear superposition of the functions for each particular region gives then the functions for the composite multimodal pore system (Durneret al., 1999):

Combining this retention model with Mualem [1976] pore-size distribution model leads now to:

Where: w_iare the weighting factors for the two overlapping regions, and α_i, n_i, m_i (=1-1/n_i), and l are empirical parameters of the separate hydraulic functions (i=1, 2), (Šimůneket al., 2008).

Potassium Sorption Isotherm

Triplicate 5 g soil sample were equilibrated in a 100 ml centrifuge tube with 50 ml of 0.01 M CaCl₂ solution containing 0, 5, 10, 20, 30, 50, 100, 150, 200, 300 and 400 mg Kl^-1 as KCl. The soil samples were then agitated on a rotary shaker to achieve equilibrium for 30 min. and K ion concentration was measured by Flame photometer and expressed as mgl^-1.

            Potassium sorption is assumed to be given by linear isotherm equation as follows (Matott, 2004; Matott and Rabideau, 2008):

Where: S is the sorbed K (mgkg^-1), C_e is the equilibrium concentration (mg/l) and K_d is the partition (sorption) coefficient (Lkg^-1) as:

The sorbed potassium was calculated according to the following formula:

Where: C₀ is the initial concentration (mgl^-1)

V is the volume of solution (ml)

W is the weight of soil (g)

Column experiment

Polyvinyl chloride columns (PVC), 65 cm long and 19 cm diameter (283.5 cm² surface area) with closed bottoms, were used in the study (Fig. 1). The base of the columns was tightly sealed with silicon adhesive. The bottom in each column was covered with soft tissue to facilitate the effluent.  A plastic tube of 0.5 cm diameter fitted into the drainage layer at 60 cm soil depth to collect the drainage water. The columns were hand-packed with air dried sandy loam soil for length of 60 cm at constant bulk density by gently tapping (1.42Mg m^-3). A soft tissue was placed on the top of the soil to decrease channeling. A Mariotte bottle is a device that allows constant flow of a fluid from a container, even when the fluid level in the container changes (Holden, 2005; Moore, 2004). The solute drained from the bottom of column was collected in glass bottle and then volume and K concentration were measured(Table 2).

Table (2). Soil column experimental data (sandy loam soil)

Fig. (1 ). Soil columns experiment layout with Marriott's device for

constant water flux

Water Flow

            The one dimensional water flow can be described by the Richards equation (Richards,1931):

Where:q is the volumetric water content (L³ L^-3),  h is the metric head (L), K(h) is the unsaturated hydraulic conductivity (LT^-1), t is the time (T) and Z is the vertical coordination (L) taken positively upward.

            The water retention characteristics q (h) and the unsaturated hydraulic conductivity function, K (h) are given by the Mualem-van Genuchten model (van Genuchten, 1980).

            The initial and boundary conditions of water flow are given as:

            The lower boundary conditions are:

Solute Transport

            The partial differential equation governing one-dimensional advective-dispersive solute transport equation (ADE) under transient water flow conditions in partially saturated porous medium is taken as:

Where: C is the total solute concentration in solution (ML^-3), S is the sorbed solute concentration (MM^-1), r is the soil bulk density (ML^-3), D is the effective dispersion coefficient (L² T^-1), q_wis the volumetric water flux  (LT^-1). The second term on the left side, Eq. (3-24) is equal to zero for non-reactive solute.

            The volumetric flux q_wis calculated with Darcy’s Law:

Where: K is the hydraulic conductivity (LT^-1) and h is the metric head (L) and z is the spatial distance (L).

            The dispersion coefficient (D) is calculated according to Bear (1972):

Where:  l_L is the longitudinal dispersivity (L),D_wis the aqueous ionic or molecular diffusion coefficient of solute in water (L²T^-1) q is the volumetric water content and τ is the tortuosity factor given by (Millington and Quirk, 1959):

Where: q_s is the saturated water content.

RESULTS AND DISCUSSION

Water retention function

Water retention function of sandy loam soil used in the present study was expressed using two hydraulic models; van Genuchten-Mualem type function model (van Genuchten, 1980) and Durner multimodal pore system (Durneret al., 1999). The results are reported in Tables (3 and 4) and Figures (2 and 3).

Table (3).van Genuchten soil parameters of  sandy loam soil

Table (4). Soil parameters of Durner multimodal pore system for sandy  loamsoils

Fig. (2). Soil water retention curve of sandy loam soil (van Genuchten model)

Fig.(3). Soil water retention curve of sandy loam soil (Durner multimodal pore system)

Potassium sorption

Potassium sorption isotherm was done according the method as described in (Gazoly, 2002) using linear sorption isotherm (Šimůneket al., 2013) and two site sorption model (Šimůneket al., 2008).  The results are reported in Table (5) and Figures (4 and 5).

Table (5). Potassium sorption parameters according linear and two-site sorption models

Fig. (4). Potassium sorption isotherm of sandy loam soil (linear model)

Fig. (5). Potassium sorption isotherm of sandy loam soil (two-site sorption model)

Potassium transport with equilibrium and non-equilibrium model

Potassium transport through the soil column was done and applying some equilibrium and non-equilibrium transport models.

Figure (6) shows the distribution of potassium concentration through the soil column. The results indicate that potassium was uniformly distributed through the soil profile with low input concentration (12.5 mg/l). While, with medium and high input concentration showed an increase in K concentration in the upper layer (down to 15 cm depth), then the concentration showed the same concentration with high magnitude of high K concentration.

Fig.(6). Potassium distribution with depth for different K input concentration (Experimental data)

The K distribution through the soil column was illustrated in Figures ( 7to 9) for different equilibrium and non-equilibrium transport models, i.e. equilibrium model (EQ), one-site sorption model (OSM), two-site sorption model (TSM), dual permeability model (DPM) and dual permeability with two-site sorption model (DPTS).

Fig. (7). Potassium distribution with depth for low K input concentration, 12.5 mg/l (modeling data) according to different transport models

Fig. (8). Potassium distribution with depth for medium K input concentration, 100 mg/l (modeling data) according to different transport models

Fig. (9). Potassium distribution with depth for high K input concentration, 200 mg/l (modeling data) according to different transport models

The predicted data according to the equilibrium and non-equilibrium transport model showed a symmetrically distribution for medium and high K input concentration. For low K input concentration, the results showed a symmetrically distribution for OSM (one-site sorption model) and DPTS (dual permeability with two-site sorption model). Also, a symmetrically distribution was found with DPM (dual permeability model) and TSM (two –site sorption model). The K distribution of EQ (equilibrium model) has a high values than other four models.

Distribution of sorbed K on soil column according to the equilibrium and non-equilibrium transport models is illustrated in Figures (10 to 12). The modeled data showed that the same trend was found in case of equilibrium model (EQ) at all K input concentration. Also, DPM and TSM were having the same trend and OSM has the higher value of sorbed K in all K input concentrations. The DPTS behave the mid trend of sorbed K. The differences between modelsmay be due to the fraction of exchange sites assumed to be in equilibrium with the liquid phase.

Fig. (10). Distribution of sorbed K on soil matrix for low input concentration

 (12.5 mg/l) according to different transport models

Fig. (11). Distribution of sorbed K on soil matrix for medium input concentration (100 mg/l) according to different transport models

Fig. (12). Distribution of sorbed K on soil matrix for medium input concentration (200 mg/l) according to different transport models

Figure (13) shows the distribution of sorbed potassium through the soil column according to the experimental data. The results indicate that sorbed potassium was uniformly distributed through the soil profile with medium and high input concentration (100 and 200 mg/l). While, with low input concentration showed a decrease in sorbed K in the upper layer (down to 15 cm depth), then the concentration showed the same trend of medium and high K concentration.

Fig.(13). Sorbed potassium distribution with depth for different K input concentrations (Experimental data)

The experimental data show higher values of sorbed and soluble K distribution through the soil profile than the modeling data. This result may be due to the soil heterogeneity as results of macropors and preferential flow of water and solute. Preferential flow, as opposed to uniform flow, results in irregular wetting of the soil profile as a direct consequence of water moving faster in certain parts of the soil profile than in others. Hendrickx and Flury (2001) defined preferential flow as ‘all phenomena where water and solutes move along certain pathways, while bypassing a fraction of the porous matrix. Thus, an important characteristic of preferential flow is that during wetting, part of the moisture front can propagate quickly to significant depths while bypassing a large part of the matrix pore-space. Water and solutes may move to far greater depths, and much faster, than would be predicted with the Richards equation using area-averaged moisture contents and pressure heads (Beven, 1981). The presence of macropores and other structural features, development of flow instabilities (i.e. fingering) caused by profile heterogeneities or water repellency (Hendrickxet al., 1993), and funneling ofFlow due to the presence of sloping soil layers that redirect downward water flow are probably the most important causes of preferential flow. While the latter two processes, i.e. flow instability and funneling, are usually caused by textural differences and other factors at scales significantly larger than the porescale, macropore flow and transport are usually generated at pore or slightly larger scales, including scales where soil structure first manifests itself (i.e. the pedon scale).

Accurate coupling of the fracture and matrix domains still represents the greatest challenge in terms of successfully describing non-equilibrium flow and transport in the vadose zone. Matrix–macropore (or matrix–fracture) interfaces can have very different properties than the bulk matrix due to the deposition of organic matter, various types of coatings, fine texture mineral particles, or various oxides and hydroxides on the aggregate exteriors or macropore walls; these coatings can markedly reduce rates of diffusion and mass flow between macropores and the soil matrix (Thomaet al., 1992).

Physical non-equilibrium occurs in the soil unsaturated zone when heterogeneity result in the generation of lateral differences (non-uniformity) either in water pressures or solute concentrations or both, during vertical flow and transport (Jarvis, 2007).

From the present study we need more basic research to cover the knowledge gaps according to the following points:

1-  More research is clearly needed on the role and importance of long-term leaching in macropores and kinetic non-equilibrium sorption effects.

2-  Further advances can be expected from researches that explain the combine geometric descriptions of soil macropores structure with continuous real-time measurements of solute transport.

3-  One more working concept is that a coarser, more heterogeneous structure (associated with a decrease in macro-porosity) promotes strong non-equilibrium macropore flow but only until critical limit is reached when macropore connectivity becomes limiting. These effects of pore volume and spectral dimension may be investigated for transport by diffusion, but not with gravity-driven convective processes such as macropore flow.

CONCLUSION

A general conclusion of the present study is that a model’s success or failure to represent the flow and transport processes internal to the transport domain should not be judged solely by the response prediction at a single outlet point. Therefore, we need more experimental data and more thorough studies of many soil types to evaluate how much and what type of information is required to fully parameterize selected non-equilibrium models. Until measurement techniques are fully developed, inverse parameter identification is an indispensable means for the application of complex non-equilibrium models. To facilitate such applications, the objective function for the inverse problem in HYDRUS-1D can be formulated in terms of a large number of variables involving not only boundary concentration fluxes, but also water and solute distributions within the soil profile and in different phases (Šimŭnekand van Genuchten, 2008).