## Abstract

Groundwater has immense importance when it comes to sustaining agricultural and domestic water demand. However, improper use of groundwater has led to drastic increase in its contamination. For efficient management of groundwater it is important to have a thorough understanding of parameters that influence the groundwater flow and solute transport. In this study, HYDRUS-1D/2D software is used to analyse the effect of critical parameters that influence the behaviour of solute in groundwater flow and its transport under different soil conditions. Also, single porosity and dual porosity approaches are discussed with their respective sensitive parameters. Parameter variations are illustrated in the forms of concentration contours of the solute spread and the breakthrough curves obtained over the flow period. The results of this work can be useful in analysing flow and solute behaviour in the vadose zone.

## HIGHLIGHTS

Flow and contaminant transport using HYDRUS.

Variably saturated porous media with Richards equation.

Used single and dual porosity model.

Effect of critical parameters on spreading of solute.

Mass transfer coefficient affects the solute mass in porous media.

## INTRODUCTION

A major portion of surface water infiltrates into the subsurface and flows as groundwater through different layers of soils which behave as porous media. This infiltrated water collects solute particles from the soil surface and water bodies such as compounds found in pesticides and fertilizers applied in the crop fields. These dissolved solutes may behave as contaminants. Effective management of groundwater quality requires accurate predictions of contaminant loading to the groundwater system, which requires reliable estimates of lag times or solute velocity, dispersion, and transformations of contaminants as they traverse through the vadose zone (defined in this study as the zone of alluvial gravel between the base of the soil and the groundwater table) to the groundwater system (Dann *et al.* 2010). Also, fine-textured soils containing clay minerals and organic matter may act as sorption filters against solutes (Köhne *et al.* 2006).

A number of modelling approaches have been developed to simulate groundwater flow and solute transport, using both analytical and numerical techniques (Cunningham *et al.* 2005; Al-Jabri *et al.* 2006; Köhne *et al.* 2009; Wang *et al.* 2013; Batany *et al.* 2019). The most common method to study flow and solute transport is through the analysis of breakthrough curves. Rühle *et al.* (2013) conducted three tracer column experiments for over seven months to investigate the variations in the pathways of water flow and solute transport occurred due to the long experiment time periods. Breakthrough curves were analysed to evaluate water flow and transport parameters. It was found that a simple advection–dispersion model is appropriate to describe the initial homogenous transport; however, after some time the porous medium changes from uniform to non-uniform. Batany *et al.* (2019) studied the influence of flow rate and viscosity of the carrying liquid on non-reactive solute transport under saturated conditions in a macro-porous synthetic medium. Analysis of breakthrough curves (BTC) denoted early breakthrough, owing to the physical non-equilibrium induced by the preferential flow within the macropores which affects the flow rate and the coefficient of molecular diffusion of the solute.

These breakthrough curves plot concentration variation against the time period obtained by conducting experiments in the laboratory or by simulating the flow using numerical models such as MODFLOW, STANMOD, HYDRUS etc. HYDRUS software covers a vast range of approaches that can be selected for simulating equilibrium as well as non-equilibrium processes (Šimůnek & van Genuchten 2008). Hence, it is one of the most widely used tools to study groundwater flow and solute transport and in one, two and three dimensions (Vanderborght *et al.* 2005; Šimůnek *et al.* 2006; Dettmann *et al.* 2014; Gharabaghi *et al.* 2015). For example, Moradzadeh *et al.* (2014) studied the effect of potassium zeolite on nitrate and ammonium ion sorption and retention in a saturated sandy loam soil in laboratory conditions. Potassium zeolite was observed to enhance retention of ions in the soil and reduce the mobility of both nitrate and ammonium. The HYDRUS-1D model was applied to simulate leaching making use of the convection–dispersion equation (CDE) and mobile–immobile model (MIM). Lai & Ren (2016) used HYDRUS-1D coupled with PEST to simulate soil water distribution along the profiles and their dynamics during seven growing seasons. The overall results denoted fair agreement among the actual soil hydraulic parameters, such as hydraulic conductivity, saturated and residual water content, and that derived from the combination of HYDRUS-1D and PEST. Hence, this software can be used to efficiently determine effective hydraulic parameters of the equivalent soil profile. Nasta & Romano (2016) used a conventional unsteady drainage experiment leading to field capacity to identify effective soil hydraulic parameters of a layered soil profile. The SWAP model was used to numerically implement a synthetic drainage process in order to attain a flux-based field capacity criterion. Liang *et al.* (2017) adapted the HYDRUS-1D model to simulate overland flow at the soil surface using the diffusion wave equation. The results were in good agreement with the analytical solutions obtained by using a kinematic wave equation, which is generally used to simulate overland flow.

The study aims to provide a detailed idea of parameters that govern flow and solute transport in porous media based on different approaches. The objectives can be stated as (1) identifying parameters that influence flow and solute behaviour in porous media, (2) analysing the behaviour of these parameters with respect to different initial and boundary conditions and (3) temporal and spatial variation of solute movement in terms of time and space with varying values of critical parameters. A comparison of these variations are also shown based on two modelling approaches viz. single and dual porosity models. The outcomes of the analysis are shown in the form of concentration contours and breakthrough curves obtained by HYDRUS 1D/2D software.

## MODELLING APPROACHES FOR ANALYSING WATER FLOW AND SOLUTE TRANSPORT

The fundamental equations that are used by HYDRUS to analyse groundwater flow and solute transport are the Richards equation and the advection–dispersion equation, respectively. This set of equations is governed by various parameters which are based on the type of conditions taken into consideration, such as equilibrium flow in the case of single porosity and non-equilibrium flow in the case of dual porosity.

### Single porosity flow and transport model

*et al.*2003):where

*z*represents the vertical co-ordinate positive upward [L],

*t*represents time [T],

*h*represents pressure head [L],

*θ*represents water content [L

^{3}L

^{−3}],

*S*is a sink term which represents the root water uptake or some other source or sink [T

^{−1}],

*K*(

*h*) is the unsaturated hydraulic conductivity function, and the saturated hydraulic conductivity

*K*

_{s}[LT

^{−1}]. Equation (1) can be expanded further for two dimensions and can be written as:where and are effective hydraulic conductivity [LT

^{−1}] in

*x*and

*z*directions, respectively, is matric potential [L], is volumetric moisture content [L

^{3}L

^{−3}], is specific storage [L

^{−1}], is differential water capacity = , [L

^{−1}],

*x*and are horizontal and vertical co-ordinates [L], respectively.

*c*represents solution phase concentration [ML

^{−3}],

*s*represents sorbed phase concentration [MM

^{−1}] = , represents the distribution coefficient [L

^{3}M

^{−1}],

*D*represents the hydrodynamic dispersion coefficient [L

^{2}T

^{−1}],

*q*represents volumetric fluid flux density [LT

^{−1}] evaluated using the Darcy–Buckingham law and

*ϕ*= sink–source term that accounts for various zero- and first-order or other reactions[ML

^{−3}T

^{−1}].

The solute particles get attached to the soil particles by the process of sorption and desorption. This process is generally explained by non-linear equations, and the Langmuir or Freundlich isotherms (Šimůnek *et al.* 2003) are the most widely used non-linear equations. Adsorption is one such process which may occur during transportation of certain solutes in porous media. Linear adsorption is considered for the sake of simplicity. Further, linear adsorption is found to deliver more symmetrical results and delay onset of instability in the flow (Rana *et al.* 2019).

*x*and

*z*directions respectively, and are hydrodynamic dispersion coefficients in

*x*and

*z*directions, respectively.

### Dual porosity flow and transport model

The single porosity flow model can be extended by considering the micro-porosity of the soil particles. According to an assumption, soil particles or aggregates have their own porosity; it means that pores are sufficiently small that water within the pores is considered immobile but available for plant extraction, which is known as micro-porosity. The concept of micro-porosity further leads to various cases of non-equilibrium flow and transport models such as the mobile–immobile water model, dual porosity model, dual permeability model etc. In the first case, these micro-pores retain water content which is immobile; however, solute present in the voids may move inside the micro-pores or out of them by the process of molecular diffusion. On the other hand, in the dual porosity concept the water and solute from the voids can move inside the micro-pore and vice versa is also true. Water flow into and out of the immobile zone is usually described using a first-order rate process. Solute can move into the immobile domain by molecular diffusion and advection with flowing water (Šimůnek *et al.* 2003).

Dual porosity considers macro- and micro-porosity, which can be represented by and , where the subscript denotes mobile region and immobile region, respectively. Total porosity can be expressed as .

*et al.*2003):where and are sink terms for both mobile and immobile regions, respectively [T

^{−1}], is the transfer rate for water between the inter- and intra-aggregate pore domains [T

^{−1}].

^{−3}], respectively, and are sorbed concentrations of the mobile and immobile regions [MM

^{−1}], respectively, is the dispersion coefficient in the mobile region [L

^{2}T

^{−1}], is volumetric fluid flux density in the mobile region [LT

^{−1}], and are sink–source terms that account for various zero- and first-order or other reactions in both regions [ML

^{−3}T

^{−1}], is the fraction of sorption sites in contact with the mobile water content (dimensionless), is the mass transfer term for solute between the mobile and immobile regions [ML

^{−3}T

^{−1}] and

*w*denotes the mass transfer coefficient [T

^{−1}]. In the above set of equations, Equation (10) governs solute transport in the macro-pore region and Equation (7) denotes mass balance for the micro-pore region. Equation (11) describes the rate of mass transfer between the micro-pore and macro-pore regions. In the concept of dual porosity, the second additive term in Equation (12) considers the advection taking place between micro- and macro-pore regions, such as

*c*is equal to for and for . The dual porosity flow and transport model can be explained using Equations (8)–(11). These equations make use of some additional parameters with respect to single porosity flow and transport. The most important of these parameters are , and

*w*. The mass transfer coefficient is a diffusion rate constant that relates the mass transfer rate, mass transfer area, and concentration change as driving force. In general, it is used to quantify the mass transfer between phases, immiscible and partially miscible fluid mixtures (or between a fluid and a porous medium). In respect to the dual porosity flow and transport model it helps in the evaluation of the mass transfer rate of water and the solute.

## METHODOLOGY

The study uses the HYDRUS model to set up various flow and solute transport conditions in the soil medium to understand the behaviour of water and solutes in a porous medium. The steps followed to prepare specific conditions in HYDRUS 1D/2D to simulate flow and solute transport in porous media are shown in Figure 1. As shown in the flowchart, in the beginning the geometric parameters are defined. This demarcates the area and plane of the flow domain. It also defines the origin of the plane and its dimensions. HYDRUS uses the finite element method to carry on the numerical calculations, hence, once the geometric parameters are set, it is important to introduce the required mesh details. After converting the domain into a proper mesh, the particulars of the main process need to be entered. This section requires a thorough knowledge of soil and solute properties that are being dealt with, the type of model approach that is required by a particular problem and other details such as root water uptake or solute uptake. The next step consists of time information which covers complete details of the time period of the flow simulation with initial and final timings, and time steps in the interval and its dimensions. The next step counts for the output information which implies the ways the output needs to be generated. One of the most critical steps of this process is introducing the boundary and initial conditions. These conditions can be specified in terms of pressure head, water content, concentration of solute etc. Darcy's law relates flow to gradient of head. Along a vertical, no-flow boundary, implies ; along a horizontal, no-flow boundary, implies (Salmasi & Azamathulla 2013). These conditions should be provided with utmost precision. Once the overall set-up is prepared, the model is ready to run. The outputs are given for different iterations in the forms of numerical values as well as plots. The outputs can be plotted for various parameters, at various depths and time intervals. These results are then analysed.

## IDENTIFICATION OF CRITICAL PARAMETERS

To identify critical parameters of single porosity and dual porosity transport models, a homogeneous soil domain of width 50 cm and depth 100 cm was considered. Initial head and concentration were specified as 100 cm and 1 gm/ml, respectively. To obtain breakthrough curves a node was introduced at a depth of 100 cm below the soil surface. Then a constant flux of 5 cm/day was introduced from the top surface and a seepage face was maintained at the bottom surface. At a seepage face boundary condition, water leaves the saturated flow domain and the pressure head of this face is taken to be uniformly equal to zero along the entire face. The seepage face boundary condition can be assigned to the nodes with negative calculated pressure, for example in the case of a free draining lysimeter or dike. It can be defined numerically as: if (*h* < 0) => *q* = 0 and if (*h* = 0) => *q* = q1.

The flow was simulated for five days and various parameters were checked for sensitivity analysis, by varying the value of one parameter at a time, thereby keeping values of all other parameters constant. Hence, the parameters whose variations influenced the flow and solute transport to the maximum extent were identified as the critical parameters.

### Single porosity transport model

The single porosity transport model represents the equilibrium flow and in the steady state depends on four basic parameters, namely water content (*θ*), flux density (*q*), dispersivity (*λ*) and distribution coefficient (), as shown in Table 1. The first two parameters govern the water flow whereas the other two parameters are directly responsible for the solute transport. Figure 2(a) represents the variation in breakthrough curves with changing values of dispersivity (*λ*). It is a characteristic property of the porous medium which differs with the spatial components of the medium and influences the dispersion of solute.

Model name . | Critical parameters . |
---|---|

Single porosity flow and transport model | θ, , λ, |

Dual porosity flow and transport model | , q_{mo}, λ, , , w |

Model name . | Critical parameters . |
---|---|

Single porosity flow and transport model | θ, , λ, |

Dual porosity flow and transport model | , q_{mo}, λ, , , w |

Dispersivity in the direction of flow is known as longitudinal dispersivity and perpendicular to the direction of flow is known as transverse dispersivity. Generally, longitudinal dispersivity is ten times the transverse dispersivity. A small value of dispersivity leads to a decrease in the value of solute dispersion and increase in the rate of adsorption in the solute in soil. Henceforth, the time required to reach adsorption equilibrium reduces with the decrease in the value of dispersivity. As evident from the figure, as the value of dispersivity (longitudinal) increases, the removal time of solute from the soil decreases and the tailing of the solute becomes less distinct. The tailing part of the BTC in the figures represents the time taken for the removal of solute from the porous medium as the concentration finally reaches zero from the initial concentration. Figure 2(b) represents the variations in breakthrough curve on varying the values of distribution coefficient, which represents the proportion of solute mass adsorbed onto the surface of particles of the porous medium. Here, the linear form of adsorption is considered for the sake of simplicity. The rate of adsorption is directly proportional to the distribution coefficient. Hence, a contrast to the variation is shown by *λ*, and higher values of increase the removal time of the solute. However, the tailing curves of the solute concentration are almost unaffected.

### Dual porosity model

In the dual porosity approach the total water content is divided into two forms viz. mobile water content by ) and immobile water content (. Hence, the parameters on which the dual porosity model depends (apart from those which also influence the uniform transport model) are immobile water content (, the fraction of sorption sites (*f*), and the mass transfer coefficient (*w*), as shown in Table 1.

Figure 3(a) illustrates the breakthrough curve with varying values of *w* which represents the first-order rate coefficient for non-equilibrium adsorption. Increment in the values of mass transfer coefficient accelerates the rate of adsorption, thereby helping the solute reach equilibrium faster. Decreasing the value of *w* takes the dual porosity transport condition towards the uniform transport condition, such that if *w* is equal to zero, Figures 2(a) and 3(a) would show the same curves. Therefore, the curve of mass transfer coefficient (0.1 per day) most likely resembles the case of uniform transport in the mobile region (macro-pores). Also, this breakthrough curve shows the most tailing, resulting from the slow release of solute from the mobile zone into immobile liquid (Šimůnek & van Genuchten 2008). Further, as the value of *w* is increased, the concentration of solute in the mobile and immobile regions starts to reach equilibrium. Hence, the greater the value of *w*, the faster is the rate of reaching equilibrium. For Figure 3(b), the fraction of sorption sites (*f*) is taken to be equal to the ratio of water content in the mobile region to the total water content (and the variation of breakthrough curves is shown altering the values of *f*. As evident from the figure, a smaller value of *f* represents a smaller value of mobile water, thereby resulting in early removal of solute and more pronounced tailing.

### Analysing of critical parameters in various soil and flow conditions

In order to precisely model groundwater flow and solute transport it is important to have precise knowledge of the behaviour of critical parameters under different soil and solute flow conditions. Here, various examples are taken with reference to previously conducted studies using various modelling techniques to analyse the behaviour of selected critical parameters (Haverkamp *et al.* 1977; Kirkland *et al.* 1992; Li 1996; Khaleel *et al.* 2002; Soraganvi & Mohan Kumar 2009). Breakthrough curves for various curves are plotted to study the pattern and extent of the variations against time with respect to various changing properties. Also, plots depicting the spread of the solute along with the depth and concentration variations are discussed.

#### Test problem 1: transient infiltration through a sand column

To simulate infiltration in a dry homogeneous soil profile, a rectangular sand column was considered of height 0.7 m and 0.05 m width. Infiltration was considered as a constant flux of 3.29 m/day from the top of the soil profile and a constant suction head of 0.61 m was maintained at the bottom. The values of saturated and residual water contents were taken as *θ*_{s} = 0.287 and *θ*_{r} = 0.075 in accordance with Haverkamp *et al.* (1977). The values of van Genuchten parameters were *α* = 3.6 per m and *n**=* 1.56 (no units). Further, the effect of *K*_{s} on the moisture content was analysed by varying values of *K*_{s} as 0.8, 0.34 and 0.05 cm/hr, as depicted in Figure 4. The simulated results were found to be in good agreement with the experimental results.

As evident from Figure 4, as the value of saturated hydraulic conductivity increases, the value of saturated water content decreases. The lower the value of saturated moisture content, the more quickly it reaches the maximum value and then remains constant for the entire soil column.

#### Test problem 2: flow and transport through sand column

*K*

_{s}= 18 cm/hr,

*θ*

_{s}= 0.27,

*θ*

_{r}= 0.06,

*α*= 0.036 per cm and

*n*

*=*1.56 (no units). Similarly, the dispersion coefficients were adopted as mentioned in the aforementioned literature, as:

With reference to Li (1996), the transport process was considered to be advection-dominant and the Peclet number to be on the order of 50. The pressure head and concentration profiles taken after a period of 1 hr are shown in Figure 5(a) and 5(b), respectively.

By keeping other parameters constant, the pressure head was plotted against depth of the soil column for varying values of as shown in Figure 5(a). A greater value of implies easier transmission of water through the soil in saturated conditions, thereby the pressure gradient is gentler for greater values of , and as decreases the pressure gradient along the depth of soil column becomes steeper. Figure 5(b) denotes the variation of concentration along the depth of the soil column for increasing values of . Greater values of indicate greater adsorption of the solute onto the soil particles, hence, the peak concentration is reached at smaller depth resulting in steep concentration profiles. The solute needs to cover greater depth in order to reach the peak concentration and vice versa.

#### Test problem 3: flow and transport through very dry layered soil

For this study, the case of a dry layered soil has been taken for reference as reported in the paper by Kirkland *et al.* (1992). To attain a perched water table, a rectangular flow domain of 500 cm × 300 cm is considered, in which a 300 cm × 100 cm layer of sand is bounded by clay layers at both the sides and bottom. And there is another layer of sand below the clay, as shown in Figure 6. The central 100 cm length on the top of the sand layer is used for introducing the solute flux at a concentration of 0.01 gm/ml at a rate of 50 cm/day. To introduce a very dry soil condition the initial head is set at −500 m. The study is aimed at analysing the variations in critical parameters of groundwater flow and solute transport when subjected to very dry conditions based on the approach of single and dual porosity.

The values for various parameters are taken with reference to Kirkland *et al.* (1992), as specified below.

The simulations using HYDRUS software are made after 1.5 days from the time that flux is introduced in the soil.

##### Single porosity model

The spread of solute below the earth surface is compared by varying the critical parameters, one at a time. The values of *λ* and are fixed at 50 cm and 0, respectively, as the default case. As the flux enters the soil domain the distribution of the solute starts spreading downwards in symmetrical patterns. Figure 7 illustrates the concentration contours after 1.5 days from the introduction of flux at the soil surface. The yellow elliptical portion formed in the perched aquifer region denotes a region of greater concentration owing to the faster movement of solute in the sand region due to its high hydraulic conductivity. Similarly, the surrounding region made up of clay shows lower concentration of the solute due to slower movement of solute. A higher value of *λ* results in higher dispersion of solute through the porous media, however, as the value of *λ* decreases from 50 to 5 cm, the extent of dispersion also decreases thereby, and the red parts show the region of high concentration in the centre. Hence, it can be said higher conductivity and low values of *λ* resulted in accumulation of higher concentration. In the next comparison, the value of *λ* is kept constant and distribution coefficient is increased gradually. It is observed that as the value of increases from 0 to 1, the spread of solute decreases. Higher values of distribution coefficient denote higher rates of adsorption, hence the solute gets adsorbed in the porous medium thereby restricting the spread of the solute. If the default values of all the parameters are kept constant and the flux is decreased the spread of solute after 1.5 days is restricted to a limited depth, however, establishing similar patterns as in the default case.

Breakthrough curves are plotted to examine the variations in the concentration against time when critical parameters are altered as mentioned in the above section. Five nodes are marked in the domain at depths of 10, 70, 125, 200 and 275 cm from the top surface for observation. In the default case it can be observed, as shown in Figure 8, that the arrival of solute is quick thereby reaching a peak value and then the removal of solute is slower and smoother. With increasing the depth from the surface the time of arrival of solute is observed to increase and the concentration peak decreases. However, the removal of solute for all the nodes takes place more or less simultaneously, making a gradual slope, unlike in the case of arrival. Now, as the value of *λ* is decreased, the arrival of solute starts getting delayed and a slight increase in the peak concentration for each node can be observed. The time of removal of solute for each node decreases, and the slope of the plot follows similar patterns as the time of arrival. However, lower values of *λ* result in early removal of solute for each node, hence, the removal of solute at each node is no longer simultaneous. As the value of *λ* further decreases, these variations become more pronounced. In the next case, the value of *λ* is kept constant at 50 and is increased from 0 to 1. It is observed that on increasing the values of , the time of arrival for the peak concentration increases and the peak concentration drastically reduces at all the nodes. Also, the time for removal of the solute from the respective depths increases. Figure 8 shows that the slope of the variations gets smoother as the value of increases. Further, reducing the value of flux leads to delay in arrival of the solute. Here, decreasing the flux to 5 cm/day, the solute reaches only up to 70 cm from the top surface in the duration of five days.

##### Dual porosity model

For the dual porosity model the most sensitive parameters are found to be *f* and *w* and 0.5 is taken as the default value for both the parameters. The value of *f* implies the fraction of water content in the mobile region to total water content. At very low value of *f* the spread of solute increases as a much lower amount of solute is present in the mobile region of the porous medium, thereby, allowing free movement of solute. However, as we increase the value of *f* a greater volume of water and solute gets collected in the immobile region of the porous medium by which the concentration of solute increases but the spread of the solute becomes restricted, as shown in Figure 9. The variation of *w*, keeping the values of *f* and other parameters constant, is also illustrated in Figure 9. Variation of *w* affects the spread of solute very slightly, however, a small value of *w* decreases the rate of adsorption, and hence, accumulation of higher concentration takes place. Similarly, taking higher values of *w* enhances the process of adsorption, thereby the concentration of solute present in the region is not very high, however, the spread of the solute remains more or less the same. If the default values of all the parameters are kept constant and the flux is decreased, the spread of solute after 1.5 days is restricted to a very limited depth, however, establishing similar patterns as in the default case.

Similar cases were evaluated considering dual porosity altering the parameters found critical for this conceptual model as illustrated in Figure 10. The concentration peak can be observed to be very sharp, especially near the top surface, and as the depth from the top surface increases the sharpness of the peak starts reducing. For very deep layers of the order of 200 and 275 cm, the BTC very gently increases to the peak value and then declines very smoothly showing very gradual arrival and removal of solute in the soil domain at greater depths. When the value of *f* decreases, the peak concentration tends to increase; as the water content in the mobile zone is less the solute gets concentrated. If we decrease the value of *f*, implying reduction in the value of mobile water content, the solute concentration gets stronger and the peak concentration increases, as is evident from Figure 10. This variation becomes more visible in the BTC plots as the depth increases. Similar variation is also observed in the case of the value of *w* is decreasing. The lower the value of *w*, the higher is the concentration peak and vice versa. If the flux is drastically decreased while keeping the other parameters constant, the arrival of solute is delayed to a large extent.

#### Test problem 4: flow and transport through homogeneous media

To study contaminant flow in a homogeneous soil medium, a field of 5 m × 5 m was considered with uniform hydraulic conductivity of ln = 0.752. The other reactive soil parameters were taken from Khaleel *et al.*(2002), such as *θ*_{s} = 0.397, *θ*_{r} = 0.027, *α*_{v} = 4.30596/m, *n*_{v} = 1.82212 and *D* = 0.000012 m^{2}/day (in *x*- and *y*- directions). The solute was made to enter at a flux rate of 1 m/day applied at the top surface, in the centre 3 m width, continuously for a period of one day. Apart from the flux inflow section, the top boundary and the side boundaries are considered no-flow boundaries. An initial suction head of 1.7 m was specified uniformly for the whole domain. The results were simulated after 1.5 days from the initial introduction of flux.

HYDRUS was used to examine the spread of solute in the soil domain and breakthrough curves were plotted for different values of critical parameters for single and dual porosity models. In case of single porosity models the effects of three parameters on the solute transport are studied viz. *λ*, *K*_{d} and *q*, which are illustrated in Figures 11 and 12. The outcomes for the homogeneous soil domain were similar to those of the perched water table domain, as explained in the above section, however, the shape of the spread in both the cases is quite different owing to the difference in soil properties. The soil properties for a homogeneous soil condition are uniform throughout the domain; hence, the spread of solute inhibits a uniform curve. On decreasing *λ*, the solute gets more concentrated in some areas under the soil as denoted by the red colour in Figure 11. Breakthrough curves were simulated after 1.5 days of initiation of flux and illustrated alike variations of various parameters that were as in the case of the perched water table. However, due to perfect homogeneity the concentration curves had gentler slope in comparison with the previous case.

As the results for the perched water table and homogeneous soil medium for the single porosity transport model were highly relatable, for the dual porosity model the spread of solute was examined only on the basis of mass transfer coefficient. The effect of variation of *w* on the solute spread in the homogeneous soil medium, keeping all other parameters constant, is shown in Figure 13. The mass transfer coefficient (*w*) influences the extent of adsorption of the solute taking place in the soil. A lower value of mass transfer suggests the model is behaving closer to the single porosity concept, thereby, the concentration is higher in a larger part of the coverage area. As the value of *w* increases, the adsorption of solute onto the soil particles increases thereby reducing the concentration in the macro-pore region of the soil, as is evident from Figure 13. Further, the variation of solute mass transfer in the domain is also illustrated in Figure 13, and it was observed that as the value of *w* increases, transfer of solute mass decreases as the dual porosity concept becomes more profound. In case of dual porosity, a reasonable proportion of solute is present in the immobile or the micro-pore region and does not take part in the mass transfer. Hence, increasing the mass transfer coefficient implies a decrement in mass transfer in the macro-pore region.

## CONCLUSION

This study provides a brief description of most influential parameters controlling flow and solute transport in porous media. As aforementioned in the objectives, critical parameters are identified after running several simulations under different soil and flow conditions. Then patterns of these critical parameters are presented. After identification of critical parameters various test problems are discussed presenting various initial and boundary condition scenarios. All these results are evaluated using two approaches to model groundwater flow and solute transport in porous media, namely the uniform flow and transport model and dual porosity model. The single porosity transport model considers the movement of water and solute through macro-pores whereas the dual porosity transport model considers mobile and immobile regions for water and solute movement. The uniform flow and transport model uses the Richards equation and advection–dispersion equation to analyse water flow and solute transport, respectively, whereas certain modified versions of these equations are used to analyse dual porosity flow and solute transport in the unsaturated porous medium. Both models are sensitive to their respective sets of parameters whose effects are more prominent on the modelled outcomes than the rest of the parameters, as specified in Table 1. The behaviour of various parameters was observed with respect to space in the form of solute spread diagrams and time in the form of breakthrough curves. HYDRUS software package was used to obtain the spread of solute in order to illustrate the concentration contours and solute mass transfer in the soil domain. The study of these parameters is significant for better understanding of the application, execution and responses of flow and solute transport models.

## CONFLICT OF INTEREST

The authors declared that there is no conflict of interest.

## DATA AVAILABILITY STATEMENT

All relevant data are included in the paper or its Supplementary Information.