Simulation of water movement and its distribution in a soil column under a water source using pore-scale network modelling

The relationship between capillary pressure and saturation has a remarkable value in investigating interactions of two immiscible fluids in porous media. Pore network models, as predictors of fluids movement in porous media, can play undeniable role in determining the mentioned relationship. In the current study, the results of numerical pore network modelling, which represents pore structure of a real porous medium with networks of pore bodies connected with pore throats, are used for computing the macroscopic relationship between capillary pressure and saturation. A notable volume of soil is influenced by water source and according to the results of previous studies, it seems practically impossible to simulate water movement in macro scale dimensions by use of pore scale models. Meanwhile a suitable solution was found in this study, by considering a thin vertical soil column under the water source which was divided to smaller volumes as cells, by horizontal crossings. Each cell was considered as a pore network unit, so the soil column was consisted of series of pore network units which were vertically jointed to each other. The moisture distribution and also wetting front movement in the column were determined by application of pore network model, using the dynamic update saturation solution method. The solution was conducted between each pair of consequent cells in an alternative manner within a time step since arrival time of the water from an upper cell to the lower one. Moreover, for evaluation of the model ability, soil moisture profiles in a sandy soil of an experimental tank under the water source were studied. Comparison of the simulation and observation data confirmed the high ability of the column pore network model in prediction of the moisture distribution and wetting front movement in a soil column.


Introduction
The vadose zone of soil is a fundamental segment of the hydrological cycle, straightforwardly influencing infiltration, storm runoff, evapotranspiration and aquifer recharge.Understanding the way of water movement in the vadose zone and the measurement of its physical characteristics are the key factors to solving a variety of engineering issues.Examples of such issues are: forecasting runoff from given precipitation events for erosion control; sediment transport and flood control; estimation of the water accessibility for plant growth and the amount of water recharge to the aquifer.Therefore, the investigation of soil-water movement has intrigued researchers from different disciplines, e.g.soil science, hydrology, agriculture and environmental sciences for years.
In recent decades, different methods have been presented for prediction of water movement in soils which can be classified in three groups including a) numerical and analytical methods, b) empirical models and c) pore scale models.
Numerical and analytical methods are typically based on Richards' equation.Numerical modeling of Brandt et al. [1] and Taghavi et al. [2] can be categorized in this group.In numerical methods, hydraulic characteristics of soil e.g.residual soil water content, saturated soil water content, saturated hydraulic conductivity, etc have to be measured experimentally before starting the iterative numerical procedure.In other words, performing time consuming experiments are prerequisite of numerical methods.
Another group of methods are empirical models, which have been applied by a lot of researchers using experimental data in specific soil types.Empirical models are derived from data observed either in field or in laboratory.The models of Ben-Asher [3], Schwartzman and Zur [4] and Amin and Ekhmaj [5] etc. fall in the category of empirical models.
The third group consists of network models like pore network models.In order to properly represent the soil, they need to have the same void ratio and grain size distribution as the soil.Therefore, macro scale constitutive parameters and coefficients can be obtained by pore network models using basic soil properties.In the mentioned pore network models, soil moisture profile can be simulated by defining the pore distributions of specified soil type.Pore network models have been used for studying solute transport in porous media and petroleum engineering.The movements of fluids such as oil, water and gas in porous media have been the most interesting issues for researchers.Also, pore network models have been used for investigation of single and multi-phase flow processes [6,7], investigation of relationship between capillary pressure, saturation and interfacial area [8][9][10][11] and drainage and imbibition [12].So, applying pore network models for investigating infiltration phenomenon from water source improved the insight to the water movement in pore structure of soil.
In spite of increasing use of new and innovative pore scale models such as pore network models, lattice boltzmann models and smoothed particle hydrodynamic approach in simulating fluid flows in porous media, these models are still less known in water and soil sciences.Although pore network models have innovative characteristics and researchers have used them for investigation of fluid flow in porous media, these models have been applied on pore scale and till now they have not been used for simulation of soil moisture profile in macro scale directly.
A pore network model for representing a porous medium such as soil usually has a volume of few cubic centimeters or less.The network with these dimensions has dozens of pore bodies and pore throats and solving fluid flow governing hydraulic equations in this pore network model demands few days to be fulfilled by new modern computers.So, practically it seems impossible to simulate fluid flow in soil in macro scale using one single pore network model.In the current research, a notable idea is presented for joining a series of pore networks to each other, and for this purpose as a case study, the vertical movement of water in a sandy soil under a water source is investigated.

Network topology
Pore network models can be classified in two groups: structured and unstructured.Up to now, most pore network models had a two-dimensional or threedimensional lattice with uniform dividing.In such a network, pore bodies are placed at the lattice nodes, which are all similarly divided, and pore throats are lined up along the grid coordinates.The quantity of pore throats associated with a pore body is called coordination number.

Network geometry
Network geometry is identified with the geometrical state of pore bodies and cross sections of pore throats.Most researchers used pore bodies of cubic or spherical shapes; Joekar-Niasar and Hassanizadeh [13] introduced the pore bodies as truncated octahedrons and Ryazanov et al. [14] utilized n-cornered star-shaped pores in their researches.

Computational algorithms
Continuum-scale formulations of two-phase flow in porous media are typically solved for the pressure and saturation of wetting or nonwetting phases.Therefore, the accompanying system of equations for an inflexible porous medium and incompressible immiscible nonreactive fluids ought to be solved: where, superscripts n and w denote nonwetting and wetting phases, respectively, φ is the porosity, μ is the viscosity of phase α, S is saturation of phase α, u is the velocity of phase α, K is the intrinsic permeability tensor, k is the relative permeability of phase α, P is the pressure of phase α and P c is the capillary pressure.
In this study, we have used a pore network model which has been proposed by Joekar-Niasar et al. [11].The pore network model has a three-dimensional regular lattice structure with constant coordination number of six.Pore bodies have cubic shape and pore throats have square cross sections.Figure 1 shows a schematic presentation of two pore bodies which are connected with a pore throat.

Numerical solution of Richards' equation using finite difference scheme
The Richards' equation is the most general method to compute soil moistures and hydrological fluxes, such as infiltration in porous media [15].Consider two-and/or three-dimensional isothermal uniform Darcian flow of water in a variably saturated rigid porous medium and assume that the air phase plays an insignificant role in the liquid flow process.The governing flow equation for these conditions is given by the following modified form of the Richards' equation [16]: where C is the specific capacity of water, h is the pressure head, r, z are radial and vertical directions, t is time, S is sink term and K denotes the unsaturated hydraulic conductivity function.
The general form of Richards' equation after discretization, linearization and simplification can be expressed as [17]: In Eq. 6, superscript n refers to the current time step and superscript n+1/2 denotes the arithmetic mean of a parameter at time steps n and n+1.Furthermore, unsaturated hydraulic conductivity function (K) is defined by: where K s is the saturated hydraulic conductivity, S e is the effective saturation, m and l are shape parameters.For K, we take the geometrical mean as proposed by Vauclin et al. [18].
It should be noted that a third-type (Cauchy type) boundary condition is used to prescribe the water flux from water source in the soil surface and constant water content has been used as initial condition in numerical scheme [19].

Laboratory experiments
For evaluating the accuracy of the proposed model, experiments of water infiltration under surface water source were conducted on a sandy soil (90% sand, 5% silt and 5% clay).Laboratory experiments were carried out using a 120 cm × 120 cm × 120 cm transparent Plexiglas box (as shown in Fig. 2).Air dried sand with a mean particle size d 50 = 0.4 mm was compacted at predetermined dry bulk density of 1.5 g.cm -3 .A polyethylene pipeline connected to water reservoir was laid on the soil surface which had a 16 mm outside diameter, a wall thickness of 2 mm and used for supplying water on top surface of the soil column.During operation, wetted soil dimensions were measured using high performance photography and analyzed by Digimizer Software.The observed soil moisture profile had a high degree of the horizontal symmetry and the flow in the centerline of this profile was considered as a pipe with only one dimensional vertical flow.So, the main aim of the experiments was to measure the depth of the wetted soil in the centerline in accordance with the location of the proposed column pore network model for simulation of vertical flow.

Evaluation parameters
Several statistical criteria can be considered for evaluating the estimated depth of wetted soil.In this study the following criteria were used: correlation coefficient (R), mean absolute error (MAE), root mean squared error (RMSE) and index of agreement (IA).
where x i is predicted distance from water source (cm), y i is observed distance from water source (cm) and n denotes the number of values.

Simulation of time dependent vertical water movement using column pore network model (CPNM)
A considerable volume of the soil is influenced by a water source and practically it seems impossible to simulate water movement in macro scale using pore scale models.So, a notable solution is proposed for simulating vertical flow in the current research and it is based on dividing a considered thin column in the center part of the volume using cells to smaller volumes, in which pore network model can be applied for each cell.In other words, the volume of the mentioned soil column consists of severed pore network models which connect to each other in vertical direction.The upper surface of each pore network model is considered as inflow boundary condition and the lower surface as drainage boundary condition.
In this research, the size distribution of pore bodies in a porous medium is computed using grain distribution of the experimental soil by considering a truncated lognormal distribution.Then, these pore bodies are arranged randomly in a porous medium.After that, the space between the centers of two adjacent pore bodies in a generated porous medium is considered twice the radial size of the maximum pore body (see Fig. 1).This established pore network size is selected in an iterative manner, so that to be as a representative volume of the studied soil.The selection procedure of this optimumsized pore network is defined in the following section of the article.To facilitate the explanations in this discussion, the optimum size of pore network is named cell.For better imagination of the column pore network model, the connections of three typical cells, in which each cell has 4 pore bodies in each direction, are illustrated in Fig. 3.

Simulation of time-dependent soil moisture's vertical profile using pore network model coupled with numerical solution of Richards' equation (PNMRE)
It should be noted that the computation of PNMRE by finite difference scheme were done by an algorithm which has been developed in Wolfram Mathematica 10.0.To attain an optimum pore network size, arrangement of computations was carried out utilizing the same algorithm.For this reason, 6 different networks having 12, 16, 20, 24, 28 and 32 pore bodies in each bearing were considered.These simulations took 4 to 170 hours long with a high speed processing computer.
It was resulted that the capillary pressure -saturation curve changed with the network smaller than the network size of 28×28×28 pore bodies and afterward these changes were negligible for greater networks.Subsequently, the resulted curve from the selected network with 32 pore bodies in each direction was utilized for simulation of soil moisture profile in macro scale.In other words, the resulted pressure-saturation curve was used in solving the Richards' equation by finite difference scheme.The mentioned procedure started with initial pressures of all nodes of the network and the next pressures were computed by equation 6.Then, the correspondent saturations of all nodes were computed by the pressure-saturation curve resulted from pore network modeling.This procedure was continued until two pressure values of all nodes were less than a constant prescribed value.Finally, by using the prepared algorithms, simulations of the pore network models were done numerically in two different scenarios including CPNM and PNMRE.In the CPNM, the simulations using developed algorithms took 15-20 days on the high speed processing computer.Consequently, distances of the wetting front from soil surface for different time durations from the beginning of the experiment were computed, which are presented in Fig. 4. Also, the corresponding statistical evaluations are given in Table 1.As it is clear from Fig. 4 and Table 1, the CPNM by having less root mean squared error and higher index of agreement in comparison with PNMRE has been simulated the depth of soil moisture profile precisely.
Figure 5 shows the simulated water content distributions using PNMRE and CPNM models in different time durations from the beginning of the experiment.Each figure contains comparisons of the water contents simulated by the both mentioned methods in the vertical direction under the water source.few centimeters in the distance) between two simulated curves are not clearly visible.However, according to table 1, the CPNM gives precise predictions in comparison with the PNMRE.

Conclusion
In the current study, capability of the pore network model in simulating time-dependent moisture vertical movement under the water source has been investigated in two separate parts including the pore network model coupled with numerical solving of Richards' equation and column pore network model.For achieving the optimum size of pore network model, series of computations have been done using the developed algorithms and by considering grain size distribution of the experimental soil.As a result, the network with 32 pore bodies in each direction has been selected as the optimum network for simulation of soil moisture profile in macro scale.In column pore network modeling, the simulations using the developed algorithms took 15-20 days on the high speed processing computer.The results proved that the CPNM in comparison with PNMRE had been simulated the depth of soil moisture profile with higher precision.As a conclusion, it can be stated that the connecting pore network cells and generating column pore network have increased the precision of pore network modeling and proved the efficiency of the presented idea about connecting the pore network models.

Figure 1 .
Figure 1.Schematic presentation of two pore bodies and the connected pore throats.

Figure 2 .
Figure 2. A sample of sandy soil moisture profile

Figure 3 .
Figure 3. Connection of three typical cells in a vertical column, in which each cell has four pore bodies in each direction

Figure 4 .
Figure 4. Observed and predicted distances of wetting front from soil surface for different time durations

Figure 5 .
Figure 5. simulated water content distributions using PNMRE and CPNM models for different time durations It is clear from the plots in Fig. 5 that, in general, the spatial distributions of the water content under the water source for the both models are in excellent agreement. )

Table 1 .
Statistical results of PNMRE and CPNM models in experiments