Evaluating the Impact of the Covariance of the Friction and Potential Gradient on Describing the Infiltration Process

During infiltration of water in soil, menisci form at the interface of water, grains, and air in the pores, inducing suction due to surface tension. Due to the random distribution of interconnected pores of different sizes, characteristic of porous media, differences in suction and friction inside pores give a diffusing infiltration front. The process of infiltration is often simulated by solving Richards’ equation in which the water flux is calculated with Darcy’s law. Underlying Darcy’s law is the assumption that the gradients in flow potential and the flow resistance due to viscous forces are independent from each other. This paper shows that these parameters are dependent and negatively correlated. A new method for calculating flows in unsaturated porous media has been developed to evaluate the impact of the covariance on infiltration predictions. The results show that the impact is significant and leads to a reduction in infiltration rate and mean friction experienced during infiltration. The method thereby provides a physical explanation for the subdiffusion observed during water infiltration in soil and is consequently expected to provide more insights into the processes of infiltration.


Introduction
Soil is characterized by a random distribution in pore sizes. During water infiltration, water menisci form at the interface of water, soil grains, and air in the pores, inducing suction due to surface tension . Differences in suction and friction inside pores give a diffusing infiltration front whereby the infiltration rate is higher in the smaller pores than in the larger pores. The process of infiltration is often simulated by solving Richards' equation (Richards 1931) which follows from substituting Darcy's law in the mass balance equation. Darcy (1856) discovered that for fully saturated porous media, the volume flux of water in soil is linearly dependent on the potential gradient over the porous media. Under fully saturated conditions, gradients in potential are deterministic and given by the imposed boundary  (Bonilla and Cushmann 2000) . The friction experienced by a laminar flow inside a pore thereby corresponds with the friction of a flow in an enclosed capillary tube. The friction is consequently inversely proportional to the pore radius R [m] squared (Pfitzner 1976;Sutera and Skalak 1993). The mean flow velocityw inside a water-filled pore with radius R is related to the deterministic spatial gradient in flow potential head dφ dη according tow Here, g [m/s 2 ] is the gravitational constant, and ν [m/s 2 ] is the kinematic viscosity of water. For low-velocity flows, the potential φ approximates the pressure head. From Eq. 1 follows that the expected value for the flow velocities inside a saturated porous medium E [w], which is characterized by the random distribution in pore radii R, follows from where n denotes the porosity of the soil and the term g 8nν E[R 2 ] resembles the expected value for the friction experienced by the flow in the water-filled capillary tubes (Dekking et al. 2005). For saturated soil conditions, the potential gradients dφ dη in the pores are deterministic and independent from each other (Bonilla and Cushmann 2000). A hydraulic conductivity parameter is now defined which is given by Multiplying the expected value for the flow discharge 1 E[π R 2 ] E[πR 2w ] with the porosity n and the expected value for the moisture content E(θ ) gives the volume flux of water Q [m/s]. It should thereby be noted that for fully saturated soils, the expected value for the moisture content is constant. Darcy's law is found when multiplying the hydraulic conductivity K [m/s] with the deterministic gradient in pressure head, which gives the specific discharge Q [m/s], or Unlike saturated soils, unsaturated soils are characterized by the presence of water menisci at the interface of air, water, and soil particles. Over water menisci, the surface tension of water induces a drop in pressure head resulting in suction. Young (1805) discovered that the drop in pressure head φ [m] over a water meniscus is inversely proportional to the pore radius R and follows from φ = − 2γ ρg R Here, γ denotes the surface tension of water, and ρ [kg/m 3 ] denotes the density of water. In soil, pore radii are randomly distributed. Consequently, the drop in pressure head at the boundary φ results in a stochastic Dirichlet boundary condition which gives a stochastic potential gradient field (Bonilla and Cushmann 2000). Over the years, different methods have been developed to account for the effects of water menisci on the infiltration process.
In 1911, Green and Ampt attempted to simulate the infiltration of water into unsaturated soil by representing soil as a series of independent capillary tubes between which no water exchange is possible. For a closed off capillary tube with a fixed diameter, the gradient in pressure head which is driving the flow is given by the difference in pressure head at either end of the tube φ| η 2 − φ| η 1 divided by the infiltration length L (Green and Ampt 1911). For a flow down an enclosed capillary tube, the cross-sectional averaged flow velocityw in a capillary tube is given byw Because Green and Ampt (1911) considered all tubes to be independent from each other, the infiltration profile follows from summing the individual solutions for the flow velocitiesw after multiplying them with the respective probability of occurrence. Equation 6 shows that the infiltration rate increases with an increase in pore radius R. A more rapid infiltration in the wider tubes than in the smaller tubes gives an increase in pressure gradients between tubes. For a porous medium like soil, this is unrealistic and conflicts observations. A criticism of this bundle-of-tubes description is therefore that it does not account for cross-flows between tubes (Bartley and Ruth 1999). Recently developed approaches therefore attempt to account for the exchange of water between pores (Dahle et al. 2005;Talbot and Ogden 2008), making the flow properties in pores dependent on the flow properties in other pores. These methods are, however, limited by the accuracy with which the exchange of water between pores is determined. In conclusion, bundle-of-tubes methods highlight two important aspects. One: An independent description of tubes gives an inaccurate description of the infiltration process, and two: Due to the dependence of the suction on the pore radius, the gradient in pressure head changes from being deterministic to stochastic. Richards (1931) took a continuum-based approach to modelling flows in partially saturated porous media. Richards' starting point was the mass balance equation which states that the change in moisture content inside a control volume equals the nett flux into that control volume, or Here, the first term refers to the change in overall moisture content, E[θ ] with time t, and the term on the right-hand side denotes the nett inflow of moisture. Richards (1931) then stated that the volume flux Q [m/s] could be described by Darcy's law (see Eq. 4) leading to Richards' equation which is given by where φ refers to the distribution in pressure heads. Equation 8 is closed by means of soilwater characteristic curves (Fredlund and Rahardio 1993;Fredlund and Xing 1994) which empirically relate the expected value in pressure head φ to the moisture content θ . By applying Darcy's law, Richards (1931) inherently assumed that the expected values of the viscous effects and the gradient in pressure head are independent of each other and that the pressure head is deterministic. The bundle-of-tubes description of soil has, however, shown that in unsaturated porous media, the viscous contributions to the friction and the gradients in pressure head are both stochastic parameters and dependent on the distribution in pore radii. In this case, the expected value for the flow follows from the product of the expected value of the hydraulic conductivity and the pressure gradient plus the covariance of the hydraulic conductivity and the pressure gradient (Dekking et al. 2005) or where E w denotes the expected value for the velocity inside a pore. Since the friction in a pore and the gradient in pressure head are negatively correlated, the covariance is negative.
Applying Eq. 2 to unsaturated porous media hence gives an overestimation of the mean velocity. Similarly, the water flux Q through a porous medium follows from The friction in a pore and the gradient in pressure head are negatively correlated. Darcy's law therefore overestimates the volume flux inside an unsaturated porous medium. The error introduced by applying Darcy's law thereby equals the covariance of the viscous effects and the potential gradients multiplied by the surface area.
Studies have been performed on the characteristics of Richards' equation in describing wetting fronts (Witelski 1996;Gilding 1991). According to these studies, Richards' equation describes a continuously expanding wet region whereby the speed of wetting is finite. Several discrepancies have been found between the solution and experimental data (Nielsen et al. 1962;Bell and Nur 1978;Roeloffs 1988). To improve the predictive capacity of Richards' equation, closed-form relationships have been developed that relate the hydraulic conductivity to the moisture content (Gardner 1958;van Genuchten 1980;Fredlund and Xing 1994) leading to a nonlinear partial differential equation. More recent improvements made to Richards' equation have a basis in the field of thermodynamics and aim to provide more insights into the physical processes of infiltrations like soil hysteresis (Hassanizadeh and Gray 1980;Niessner and Hassanizadeh 2008). Caputo (2000) noted that the permeability of a soil diminishes with time as if the porous medium has a memory and went on to produce diffusion models by including memory formalism. The memory effect is attributed to chemical reactions between the porous medium and the pore fluid and due to the transportation of particles through the channels (Caputo 2000;Sapora et al. 2017). Pachepsky et al. (2003) introduced a dependence of the diffusivity on time or distance to account for the effects of the memory by replacing the first-order time derivative by a fractional time derivative. Fractional diffusion equations are related to random walk approaches and have shown to be a synthetic and efficient tool to model memory effects and non-local interactions (Scalas et al. 2000;Evangelista et al. 2011). Mathematical models based on fractional calculus have led to improved predictions of the time-dependent decrease in diffusion observed during experiments. These models thereby aim to incorporate the past history of the infiltration process on the present state of the system (Freitas et al. 2017). Although a formal agreement with experiments is obtained, the physical proof that the effects of memory are the source of the a time-dependent decrease in diffusion is still lacking.
Although these models do give an explanation for the time-dependent decrease in diffusion, they all are still based on substituting Darcy's law in the mass balance equation and inherently do not yet account for the contribution of the covariance. This paper analyses the impact of the assumption underlying Darcy's law that the covariance of the gradient in pressure head and the friction approximates 0. To perform this analysis, a new modelling methodology has been developed. In the new method, one mass and two momentum balance equations are solved to describe a problem in one spatial dimension. One of the momentum balance equations describes the flow through a pore, and the other momentum balance equation describes the exchange of water between pores. The momentum balance equations are linked to each other via the pressure, for which a pressure relationship is solved. At the location of water menisci, the method accounts for the discontinuities in pressure field. The number of parameters needed to solve these equations has been kept equal to the number of parameters needed to solve Richards' equation. Soil-water characteristic curves are thereby used to derive a discretized pore space distribution (Fredlund and Xing 1994;Nimmo 2013) that serves as input for the model. The model output consists of a discretized distribution in flow velocities. The max flux then follows from the discrete sums of the discharges w j A j in the pores multiplied by their respective probability p j and the porosity n and divided by the expected value for the surface area inside a pore.
Consequently, the model is well suited for evaluating the impact of the expected values. Section 2 describes the development of the new method. The results are, respectively, presented and discussed in Sects. 3 and 4.

Methodology
Porous media are here characterized as a random distribution of interconnected pores of different sizes. Solutions for flows in porous media are therefore given by the statistical ensemble of solutions for the flow in and between each of the pores. The ensemble of pore sizes has been discretized into k distinct sizes, each with a distinct radius R [m] and corresponding probability p such that i=k i=1 p i = 1. The shape of a pore is assumed to be cylindrical which is in line with the bundle-of-tubes-like description (Dahle et al. 2005;Talbot and Ogden 2008). The porous flow of a liquid with density ρ [kg/m 3 ] is consequently well described by the mass and momentum balance equations in cylindrical coordinates. In this section, the methodology is illustrated for the case of vertical infiltration of water into soil whereby the η coordinate direction has been replaced with Cartesianẑ direction.

Constitutive Equations
A laminar flow down the centre axis of a porous pipe is considered to have a parabolic velocity distribution given by where w [m/s] is the flow velocity as a function of r ∈ [0, R], where R [m] denotes the radius of the tube. The flow is maximum for r = 0. The fraction at the right-hand side scales the parabolic distribution in such a way that the cross-sectional averaged flow velocity in the pore equalsw.

Mass Balance
The exact flow field inside a cylindrical-shaped tube is of little interest as the assumption of a cylindrical shape is only a simplification of the pore structure. This warrants integration over the cross section of a pore. Due to the porous nature of soil, each individual pore in a porous medium is connected with surrounding pores of different radii. The circumference of a pore, when represented as a cylinder, consists of a fraction of pore spaces denoted by the porosity n and a fraction of solid tube walls given by the grains and denoted by (1 − n). The pore fraction in turn consists of j fractions, each of which resembles a connection with another pore with a distinct pore diameter, such that the sum of the individual fractions leads up to 1. Consequently, np i refers to that fraction of the circumference of the pore for which the pore is connected with a specific neighbouring pore denoted by index i. This has been Fig. 1 Representation of the wall fraction and neighbouring pore fractions of a porous medium divided into five discretized pore sizes (each pore has four neighbouring pores, or j = 4) illustrated in Fig. 1. Between pores, both mass and momentum are exchanged. The degree to which this exchange takes place depends on the probability of occurrence of the pores and their corresponding properties (See Fig. 1). Due to the presence of solid particles separating the pore spaces at the edge of the pore space, the flow in angular direction inside a pore is assumed negligible. The total exchange of water between pores is given by the discrete sum of the individual solutions for the flow between the pores. With R denoting the radius of the pore, A = π R 2 = i=k i=1 π p i R 2 , and withũ r referring to the cross-sectional averaged nett outflow from the tube in radial direction, the cross-sectional integrated mass balance equation follows from integrating the mass balance equation in cylindrical coordinates over the cross-sectional surface area of a cylinder and is given by Here, p i is the probability of occurrence of a surrounding pore, ρ [kg/m 3 ] is the density of water andw [m/s] is the cross-sectional averaged flow velocity down a pore.

Momentum Balance
For a laminar flow down a cylindrical pore, the shear stress is related to the flow velocity according to (Schlichting 1979) where μ [Pa-s −1 ] is the dynamic viscosity of water. Substituting Eq. 12 in Eq. 14 shows that the viscous effects experienced by the flow in a tube are linearly dependent on the cross-sectional averaged flow velocity in a tubew, which is in line with Darcy's law, and Hagen-Poiseuille (Sutera and Skalak 1993). Consequently, integrating the momentum balance equation in cylindrical coordinates for flows in theẑ coordinate direction over the cross-sectional area of a cylindrical-shaped pore results in The friction experienced by the flow due to viscous forces is assumed to be due to the flow interaction with the solid particles, denoted by (1 − n) in the last term of Eq. 15. The horizontal flow between tubes follows from integrating the radial momentum balance equation over the cross section of the tube. For the acceleration terms, this gives Here, B denotes the body force in radial direction. The last term on the right-hand side denotes the friction experienced by the radial flow. It should be noted that due to the unknown size of the connections between the pores in radial direction, the term K r [m/s], which resembles the local conductivity in a pore, is at this stage yet unknown.

Pressure Relationship
Due to the exchange of water in and out of the pores and the effects of this exchange on the pressure gradient, the cross-sectional integrated pressure in a pore is a function of the pore diameter and z. An expression for the pressure follows from the Poisson equation which is found by taking the divergence of the generic momentum balance equations. Integrating the resulting expression over the angular direction from 0 to 2π p i and over the radial direction It has thereby been assumed thatũ r | r =0 = 0 and that A and p are constant in time.

Integrating over the Infiltration Length
At this stage, each pore in a porous medium is assumed to be connected with few surrounding pores. Due to the random arrangement of the soil matrix and corresponding random distribution in pore sizes, for each pore it is unknown what the sizes of the directly adjacent pores are. Consequently, no local information is available on the pressure in each of the directly adjacent pores, or on the friction experienced by the flow towards or from the adjacent pores. This inhibits quantifying the flow exchange to and from each pore individually. Joekar-Niasar et al. (2008) addressed this problem by resolving to a network model in which the pore structure is randomized. Here, this problem is addressed by integrating each specific pore diameter over the infiltration length. Each individual pore size is thereby evaluated on the sum of the inflow and outflow from and to other pore sizes. It is thereby assumed that, over the length of infiltration, each specific pore size is connected with the full distribution of pore sizes larger and smaller than the pore size under evaluation. After removing the pore size under evaluation, the distribution of surrounding pore sizes is obtained. The cumulative sum of the probabilities of occurrences of the surrounding pore sizes per definition adds up to one, which has been accounted for by normalizing the probability of each surrounding pore p i by the discrete sum of the probabilities of the surrounding pores.
Each time, the flow properties that correspond to a single pore radius are evaluated. Pore radii are thereby not assumed to vary in time or space.

Mass Balance Integrated
Integrating the mass balance equation in each discrete pore size over the infiltration length inẑ direction from z 1 to z 2 , where z 1 denotes the soil surface and z 2 the infiltration front, gives after applying Leibniz integration rule By integrating each pore size over the infiltration length, the density no longer depends on z. The overbar . in Eq. 19 denotes the infiltration length averaged value. A new outflow parameter is now defined asq i =ũ r L. Applying this parameter to the over the infiltration length integrated mass balance gives

Momentum Balance Integrated
A positive consequence of integrating over the infiltration length is that the exchange of water between the pores occurs over the full distribution of pore sizes. The pressure gradient between two distinct pores averaged over the infiltration length is considered deterministic as the pores are fully saturated over the infiltration length. Hence, the volume flux between the pores could be determined using Darcy's law whereby the average friction experienced by the flow is denoted by the saturated hydraulic conductivity K s [m/s]. Integrating the radial momentum balance equation over z from z 1 to z 2 now gives The integral of the pressure gradient in radial direction given in Eq. 21 yet needs to be determined. Applying integration by parts gives Here, E denotes the expected value. Hence, when integrating the pressure gradients over the domain [0, R], the sum of the integral of the pressure gradients becomes equal to the expected value for the nett driving pressure multiplied by the area of each pore. When ∂ p i ∂r = 0, as would be the case for a uniform-distributed pore sizes, then no pressure differences are present between the pores and the flow in horizontal direction is 0.

Pressure Relationship Integrated
Equation 17 is now integrated over theẑ coordinate direction. Hereby, ∂ L ∂t = w| z 2 , and u r | z 1 = 0. Water is thereby assumed incompressible leading to Spatial gradients in flow accelerations in groundwater flows are small allowing for gradients in advective acceleration to be omitted. Furthermore, substituting Z = −ρg and substituting Eq. 23 in Eq. 17 give an expression for the over the infiltration length averaged pressure inside a poreP. This expression also guarantees that mass is overall conserved during the exchange of water between the pores.
For each pore size, the distribution of the surrounding pores is unique as the pore size under investigation is removed from the pore size distribution (see Fig. 1). This leads to k unique equations for k pore sizes, and it becomes possible to resolve the flow for each individual pore.

Body Forces
The final step needed to solve the mass and momentum balance equations is an expression for the body force B. For a situation at rest, the capillary rise is inversely proportional to the pore radius. Pressure profiles thereby linearly decrease over the capillary rise height and the mean pressure isP =P | z 1 +P| z 2 2 (see Fig. 2). With no radial body forces (B) in place, differences in boundary pressure P| z 2 cause a net horizontal force which results in a flow between tubes of different diameters despite an equal pressure gradient in the pores, and an equal bottom boundary pressure (see Fig. 2). Consequently, a body force must be in place to correct this. A body force in horizontal direction could only originate from the capillary forces on that fraction of the sides of the pores in contact with air over the height equal to the difference in infiltration length between pore radii. Hence, the body force is a function of P(L i − L) where L is the pore under investigation and L i denotes the infiltration length in a pore with index i. Between two pore sizes, the pressure acting over the height L − L i would equal the average in pressure boundaries given by 1 2 P| z 2 + P i | z 2 . Hence, the integrated body forces are given by The gradient in body forces with respect to R, as given in the integrated pressure Poisson equation (see Eq. 24), is given by Substituting Eqs. 25 and 26 in, respectively, Eqs. 21 and 24 gives, respectively, the radial momentum balance equation and pressure relationship used for resolving the volume flux between pore sizes. The equations have been solved numerically according to the method prescribed in "Appendix A."

Results
In this section, the numerical solutions of the equations from Sect. 2.2 are given for two cases. For both cases, the pore size distribution is based on the Van Genuchten equation. When the steady-state potential is expressed in terms of the pore radius by means of φ = 2γ R , then this equation becomes (Fredlund and Xing 1994)  Here, m, w, and s are curve fitting parameters whereby m = 1− 1 w . The resulting distribution is depicted in Fig. 3 for w = 3, and s = 1. The values for P(R < R) range from 1 no+2 to no no+2 , where no is the discrete number of pore sizes indicative of the distribution. This results shown correspond with no = 200. A further increase in the number of tubes only provided a small increase in accuracy. The values corresponding with the lowest and highest cumulative probability have been used to, respectively, determine the upper and lower limit values for φ, which have been converted into radii using Eq. 5. The radii are thereby assumed to increase at equal distant steps from R min to R max . For various deterministic constant positive pressures described at a boundary location below a horizontal phreatic surface, the simulated capillary rise results in a rest saturation profile that corresponds with a moisture distribution given by Van Genuchten equation. The rate of capillary rise is thereby faster in the smaller pores than in the larger pores, indicating that the exchange of mass and momentum is accurately predicted by the model. Mass and momentum are thereby exchanged from the wider pores towards the narrower pores. Figure 4 demonstrates the model output for a case of horizontal infiltration in a soil with a pore space distribution given by Fig. 3. Here, the infiltration profiles are given at 80 s intervals whereby the outside water pressure is given by 0.50 m water column. For a higher positive pressure head, the degree of advection of the infiltration front was found to increase. The corresponding internal exchange of water between pores of different sizes is depicted for each of the time steps in Fig. 5. As can be seen, the larger pores show a net outflow and the smaller pores a net inflow. Water is hence rerouted from the larger pores towards the smaller pores. The exchange of water between the pores thereby decreases with time leading to the characteristic decrease in diffusion of the infiltration profile with time as illustrated in Fig. 4.
At the boundary where the soil is saturated, a water flux conforming to Darcy's law is found by multiplying the expected value of the potential gradient by the expected value of the friction, both of which follow from the model. To determine the expected value for the friction which conforms to Darcy, the effects of the exchange of radial momentum were set at 0. The prediction of the specific discharge according to Darcy as depicted on the horizontal axis in Fig. 6

now represents
These results have been set out against the expected value for the mass flux that follows directly from the model, given on the vertical axis, which follows from Figure 6 shows both predictions set out against each other for t > 2 s. During the first 2 s, the model is initializing. The initial conditions of the model were chosen such that no pressure gradients were present between the pores. After 2 s, the effects of this assumption were negligible. The difference in prediction indicates the effects of the covariance. The solid line denotes the perfect fit line. As time progresses, the nett effects of the covariance become smaller as indicated by the decrease in deviation from the perfect fit line when following the predictions in Fig. 6 from the top right to the bottom left. The relative effect, indicated by the ratio between predictions, however appears to remain constant. Substituting Darcy's law in the mass balance equation gives Richards' equations (see Eq. 8) (Richards 1931). By substituting the Boltzmann variable λ = η t 0.5 , Eq. 8 becomes where the hydraulic conductivity K is a function of the moisture content. For Eq. 30 to be valid, the relation between the expected value of the location of the infiltration front L and time Here, lg refers to the natural logarithm. In Fig. 7, the values for α have been set out against time for t > 2.8 s. As indicated by the values for α < 0.5, subdiffusion takes place during the initial stages of infiltration. Over time, the behaviour of the system converges to that of Fickian diffusion. Pachepsky et al. (2003) noted that Richards' equation is not general enough to simulate water transport in various soils. One of the improvements suggested by Pachepsky et al. (2003) was to introduce a dependence of the diffusivity on time or distance via a fractional derivative of the water content. This approach has been the focus of more recent studies which highlight the effectiveness of such an approach on capturing the subdiffusion of the infiltration front (Freitas et al. 2017;Evangelista et al. 2011;Sapora et al. 2017). Figure 7 shows that during the initial stages of infiltration, the change in expected location of the infiltration front with time is described by a subdiffusion process as α < 0.5, whereby L :: t α . With time, the process asymptotically approaches a Fickian diffusion problem indicated by α = 0.5. The values for α following from the model are in line with the values given by Pachepsky et al. (2003) who found values between α = 0.344 and α = 0.479 for infiltration depths up to 50 cm. Consequently, the model appears to capture the time-dependent diffusivity of the infiltration front well by simulating the exchange of mass and momentum between the pores. Richards' equation follows from substituting Darcy's law in the mass balance equation. An error is introduced in Richards' equation by ignoring the impact of the covariance of the potential gradients and hydraulic conductivity in calculating the volume flux. Deviations from the Fickian diffusion indicated by the values for α < 0.5 are therefore indicative of the impact of the divergence of the covariance. The degree in deviation from Fickian diffusion during the initial infiltration (see Fig. 7) indicates that the effects of the divergence of the covariance are significant. It is important to note that the degree of subdiffusion decreases as the infiltration front progresses. Consequently, it is related to the degree of exchange of water between the different-sized pores and pressure differences between the pores. A more extensive study to quantify the effects of subdiffusion is recommended.

Discussion
In the model, the water exchange between the different-sized pores follows from the mass and momentum balance equations. Pressure differences between the different pore sizes decrease as the infiltration front diffuses. This causes a reduction in the exchange of water from the larger towards the smaller pores and a decrease in subdiffusion of the infiltration front. Bundle-of-tube models (Talbot and Ogden 2008) and dual-permeability models (Gerke and Van Genuchten 1993) account for a redistribution of water between the different-sized pores. Characteristic of these types of models is that mass is redistributed, but momentum is not. In this type of models, the expected value of the pressure gradient in the individual pores is assumed to be independent on those in the surrounding pores. In the presented model, the pressures between the different pore sizes are made dependent on the other pores by solving the momentum balance equations for the flow between the pores. The presented model thereby shows similarities with the model of Sapora et al. (2017) who made the flow rate in a point dependent on the gradient in pressures in all nonadjacent points. Sapora et al. (2017), however, assumed that the sum of horizontal pressure differences is 0, indicating that the sum of the exchange of water between the pores is 0. Here is assumed that the sum of the exchange of water between the pores is equal to the sum of the spatial gradients in the exchange of water. Diffusion of the pressures between the different pore sizes and the time-dependent reduction in mass exchange between the different pore sizes thereby explain the subdiffusive behaviour of the infiltration front. The results may be used to extend the method of Sapora et al. (2017) by adding fractional time derivatives. Figures 4 and 5 show the effect of the pressure differences between the different pore sizes and how diffusion of the pressure differences between the pores results in a continuous decrease in the rate of exchange of water between the different pore sizes. Diffusion of the rate of exchange implies that there is one pore size for which the nett inflow and outflow are the same. The rate of diffusion is thereby limited by the balance between the pressure differences and the flow resistance. A larger positive pressure head on the soil surface gives an increase in the rate of displacement of the infiltration front which corresponds with findings by Freyberg et al. (1980). During desorption, pressure differences between the pores change. Water exchange tends to be dominant in the direction of the larger pores towards the smaller pores due to the higher capillary action in the smaller pores. The lack of water exchange during desorption therefore results in a higher average flow resistance and a higher matric suction. This corresponds with the soil hysteresis phenomena described by Tami et al. (2004). To improve the understanding of the hysteresis phenomena, further studies to the dependence of pore pressure gradients and the exchange of water between the pores, indicated by the effects of the covariance, are therefore recommended.
A requirement for any model to capture the infiltration process correctly is the need to predict the volume flux correctly. Due to differences in pore sizes, the average infiltration rate is a function of the rate of exchange of water between the different pore sizes. In Richards' equation, the average volume flux is determined by Darcy's law which is based on the assumption that the friction in the pores and the potential gradient are independent on each other. The developed model was used to evaluate this assumption. Figure 6 shows that the impact of the dependence between the potential gradient and friction is significant. Ignoring the effect of the covariance gives an overestimation of the specific discharge predictions of approximately 48% (See Fig. 6). It should be noted that the quantification of this impact is influenced by the method of choice whereby pore sizes are averaged over the infiltration length and soil is schematized as a distribution of interconnected tubes. However, considering that the model physically captures the subdiffusion of the infiltration front indicates that in the case of unsaturated soils, the average pressures in a pore are dependent on those in the neighbouring pores. The effects of the dependence between the friction and capillary pressures on predicting the volume flux are thereby significant. The significant impact of the covariance also introduces questions about the impact of the continuum assumption underlying Richards' equation (Gilding 1991). The problem with the continuum approach is that the divergence of the covariance cannot be accurately determined and the effects of the covariance on the inflow boundary conditions are a priori unknown. An advantage of the method presented herein over the Richards' equation solvers is that the pressure boundary conditions and body forces are well defined in the model. The constitutive equations and solution method given here for the case of a 1D infiltration problem could thereby be extended to 2D or 3D problems. The model is thereby able to reproduce the subdiffusion of the infiltration front. Spatial and temporal fractional approaches may also be able to capture the impact more accurately. Further research on how to capture the effects of the covariance of the friction and potential gradient in these types of models is therefore recommended.

Conclusion
This paper presents the results of a study to the impact of the dependence between the friction and potential gradient, which is ignored by Darcy's law underlying Richards' equation. Both the friction and potential gradient are a function of the pore radius and are randomly distributed. Darcy's law must therefore be extended to include the covariance of the potential gradient and the friction. A new modelling method was presented which accounts for this dependence in simulating the infiltration process in unsaturated porous media. In this model, pore pressures have been made dependent on the pressures in the surrounding pores. Pressure differences between pores thereby balance the friction losses due to the exchange of water between the pores. The model output shows that ignoring the effects of the covariance gives a significant over-prediction of the volume flux and prevents the time-dependent diffusion process from being captured. The impact of the covariance, however, decreases as the infiltration depth increases. The degree of subdiffusion thereby approaches Fickian diffusion. Further studies to the impact of the divergence of the covariance are recommended to further improve the accuracy with which the infiltration process can be modelled. superscript n and the new time step by superscript n+1 . The new time step is reached after the iterations have converged at which moment the variable at the new time step equals the last guess from the iterations. The length L has implicitly been accounted for in solving the radial momentum balance equation to improve model stability. The length at a new time step follows from the mass balance equation which has been discretized as L n+1 = L n − t 2n Rq n+1 r + tw n | z 1 ; please note thatw| z 1 [m/s] is defined at the old time step and is considered as a quasi-steady inflow boundary condition. The value forq n+1 r is given by the final iteration step forq r . In the radial momentum balance equation (see Eq. 31), the local and advective acceleration term and the friction term have been discretized as ρ Aq * * −q n t + 2π Rp iq * L * + ρgn K s q * * A semi-implicit expected value for the driving force in the momentum balance equation in radial direction (See Eq. 31) follows from replacing the infiltration lengths by the expression for L n+1 given by Eq. 35. This gives Rq * * i + tw n i | z 1 p i − 2π R P * L n − t 2n Rq * * + tw n | z 1 The body forces in Eq. 31 are represented by the capillary action at the side walls of the tube and are hence given by the difference in infiltration length between the tubes. The same procedure as followed for deriving the implicit discretized pressure relationship has been used to find an implicit discretized relationship for the body force. This results in whereby i=k i=1 p i = 1. Hence, for tubes with a larger infiltration length than the tube under regard, the body forces are negative reducing the outflow from a tube, and for tube sizes smaller than the tube under regard, the body forces are positive giving a stronger outflow.
In the Poisson equation, the pressure terms consisting of the right-hand side of Eq. 24 have been discretized as Equation 39 has been solved by rewriting it as a matrix vector multiplication given by AP * * = V . Matrix A thereby contains the known parameters and vectorP * * the unknown pressures. Vector V follows from the input values given by the previous iteration step of the momentum balance equations and contains the pressure boundary conditions given by the last two terms in Eq. 39 and the parameters given by the left-hand side of Eq. 33 which have been discretized as (40) the total volume flux now follows from the sum of the flow down all tubes scaled by the area of the individual tubes.
An advantage of the proposed method is that the exchange of water between the pores is based on the mass and momentum balance equations instead of an a priori assumed redistribution profile. The exchange of mass and momentum is thereby also accounted for when determining the overall redistribution of water in soil.