MODELING WATER MOVEMENT IN HORIZONTAL COLUMNS USING FRACTAL THEORY ( 1 )

Fractal mathematics has been used to characterize water and solute transport in porous media and also to characterize and simulate porous media properties. The objective of this study was to evaluate the correlation between the soil infiltration parameters sorptivity (S) and time exponent (n) and the parameters dimension (D) and the Hurst exponent (H). For this purpose, ten horizontal columns with pure (either clay or loam) and heterogeneous porous media (clay and loam distributed in layers in the column) were simulated following the distribution of a deterministic Cantor Bar with fractal dimension ≈ 0.63. Horizontal water infiltration experiments were then simulated using Hydrus 2D software. The sorptivity (S) and time exponent (n) parameters of the Philip equation were estimated for each simulation, using the nonlinear regression procedure of the statistical software package SAS®. Sorptivity increased in the columns with the loam content, which was attributed to the relation of S with the capillary radius. The time exponent estimated by nonlinear regression was found to be less than the traditional value of 0.5. The fractal dimension estimated from the Hurst exponent was 17.5 % lower than the fractal dimension of the Cantor Bar used to generate the columns.


SUMMARY
Fractal mathematics has been used to characterize water and solute transport in porous media and also to characterize and simulate porous media properties.The objective of this study was to evaluate the correlation between the soil infiltration parameters sorptivity (S) and time exponent (n) and the parameters dimension (D) and the Hurst exponent (H).For this purpose, ten horizontal columns with pure (either clay or loam) and heterogeneous porous media (clay and loam distributed in layers in the column) were simulated following the distribution of a deterministic Cantor Bar with fractal dimension ≈ ≈ ≈ ≈ ≈ 0.63.Horizontal water infiltration experiments were then simulated using Hydrus 2D software.The sorptivity (S) and time exponent (n) parameters of the Philip equation were estimated for each simulation, using the nonlinear regression procedure of the statistical software package SAS ® .Sorptivity increased in the columns with the loam content, which was attributed to the relation of S with the capillary radius.The time exponent estimated by nonlinear regression was found to be less than the traditional value of 0.5.The fractal dimension estimated from the Hurst exponent was 17.5 % lower than the fractal dimension of the Cantor Bar used to generate the columns.Index terms: numerical modeling, Hydrus 2D, fractal model, Cantor bar, sorptivity, soil physics.
(1) Received for publication in May 2009 and approved in May 2010.

INTRODUCTION
Fractal mathematics has been used to characterize water and solute transport in porous media (Arya et al., 1999;Pachepsky et al., 2000) and also to characterize and simulate porous media properties (Anderson et al., 2000).From a non-fractal perspective, the position of any volumetric water content within a porous media column can be characterized by: x = λ(θ).Similarly, the process of one-dimensional infiltration of water into a uniform soil as a function of time can be characterized by the Philip (1957) equation: where I = Cumulative one-dimensional infiltration [L], and S = Sorptivity [L T -n ].
In both (Equations 1 and 2) the time exponent (n) was fixed at 0.5 (Clothier & Scotter, 2002).However in fractal characterization of water movement processes in the soil, the time empirical coefficient was estimated statistically, along with other hydraulic properties (Guerrini & Swartzendruber, 1994;Pachepsky & Timlin, 1998).
The diffusive process in a porous medium can be interpreted as a self-affine fractal.If this assumption holds, there is experimental evidence that the empirical coefficient n can be associated with fractal parameters such as the Hurst exponent (Guerrini & Swartzendruber, 1994).The objective of this study was to evaluate the correspondence between n and fractal parameters: fractal dimension (D) and the Hurst exponent (H).For this purpose, horizontal columns with heterogeneous porous media (i.e. two contrasting grain sizes) were generated following the horizontal distribution of the first four iterations of a deterministic Cantor Bar.Experiments of water transport were then simulated using Hydrus 2D software (Rassam et al., 2003) and n and S parameters were estimated for each simulation run using nonlinear regression.

MATERIALS AND METHODS
Water infiltration was simulated using Hydrus 2D software (Rassam et al., 2003).The simulations were performed at the Earth and Planetary Sciences department of the University of Tennessee in Knoxville, US, from August 2004 to May 2005.Ten porous media scenarios were created in Hydrus 2D following the first four iterations of a deterministic Cantor Bar (CB) with fractal dimension: D = log(2)/ log(3) ≈ 0.6309 and b = 3 (Figure 1).Two materials with distinct particle sizes were used to generate the iterations of the CB, pure clay and pure loam.The criterion to define clay and loam materials was based on the clay, silt and sand content of each material (Figure 2).The van Genuchten (1980) parameters and saturated hydraulic conductivity (K s ) were higher in the loam material (Table 1).
(Figure 1).The columns were defined by a combination of two materials (clay and loam) and five iterations (i = 0 to 4) resulting in 10 treatments: two materials x five iterations = 10 columns.At the zero th iteration, infiltration was simulated in columns filled in two ways: pure clay (CLAY0) and pure loam (LOAM0) (Figure 3a).For the subsequent iterations, the two materials were layered horizontally within the column following the distribution of the CB.In the first iteration, one loam layer was inserted in the clay column (CLAY1) and one clay layer was inserted in the loam column (LOAM1) (Figure 3b).The filling pattern following the CB distribution was used to define each treatment until the fourth iteration (Figure 3e) as follows: CLAY0; CLAY1; CLAY2; CLAY3; CLAY4; and LOAM0; LOAM1; LOAM2; LOAM3; LOAM4.With increasing iterations, the content of the interlayered material (i.e.loam in the clay column and clay in the loam column) increased until exceeding the original material in the column (Figure 3).
Each simulation represented a 10-min water infiltration experiment for each column.The initial time step was 0.1 min and the minimum and maximum were 0.01 and 0.5 min, respectively.The "cumulative water boundary fluxes" output files of Hydrus 2D (Cum_Q.out)were used for the statistical analysis of the results.The files were imported into Microsoft Excel ® and converted to *.xls format.The Materials were mixed in 10 x 10 cm columns following the iterations of a deterministic CB   time and cumulative volume infiltration were fitted to equation 2 using nonlinear regression in the "nlin" procedure of SAS ® .The parameters sorptivity (S) and time exponent (n) were estimated by nonlinear regression procedures.Sorptivity was also estimated with the time exponent fixed at 0.5 (Clothier & Scotter, 2002).The estimated fractal dimension was calculated following Feder (1989): where D = Fractal dimension, d = Euclidean dimension (d = 1 for the Cantor Bar), H = Hurst exponent.
The Hurst exponent (H) was associated with the time exponent (n) following (Guerrini & Swartzendruber, 1994).The estimated fractal dimension could therefore be calculated as: where all definitions were as above.

Infiltration simulations
Simulations of infiltration versus time in pure clay and clay mixed with loam following the iterations of the Cantor Bar (CB) are illustrated in figure 4. The cumulative infiltration rate increased as the loam content increased with each iteration of the CB.An inverse behavior was observed in loam simulations where the cumulative infiltration rate decreased as the amount of clay increased within the columns (Figure 5).Soil sorptivity (S) was defined as the gravity-free, capillary-induced absorption (Clothier & Scotter, 2002).The greater infiltration in the columns with higher loam content is reflected in the Sorptivity (S) plots as a function of iteration level (Figure 6).Sorptivity decreased exponentially with each iteration level in loam simulations and increased exponentially for clay.White & Sully (1987) found that S decreased in soil materials following the order: fine Brindabella fine sand > Brindabella silty clay loam > Molonglo loam > Bungendore fine Sand (70 %) > Kaolinite (30 %).The overall sorptivity for increasing water potentials was also greater for coarser material (Pilliago sand) than in finer soil material (Murrumbateman silty clay loam) in a study conducted by White & Perroux (1989).Mapa et al. (1986) found that the overall S was greater in a silty clay loam (Typic Torrox) than clay soil (Vertic Haplustox).Talsma (1969) also found a decrease of S with particle size, and that S was one order of magnitude greater in sand than in clay loam soil.
The results of the simulations in Hydrus 2D and from the literature are evidence enough that sorptivity is greater in coarser materials.However this is somewhat contrary to the sorptivity definition and also to the intuitive knowledge that pure clays would have a greater capillary absorption potential than coarser materials.Some hypotheses on why the   sorptivity would be higher in a coarse than in finer media were therefore investigated.
Sorptivity is an integral soil-water property that contains information about the soil hydraulic conductivity, soil water potential and soil diffusivity (White & Perroux, 1987).In a first analysis, the dependence of S on saturated hydraulic conductivity (K s ) was investigated.S and K s are not independent variables; K s is roughly proportional to S 2 (Clothier & Scotter, 2002).Following this premise, and using the K s values 3.33 x 10 -3 and 1.72 x 10 -2 cm s -1 for clay and loam, respectively (Table 1), the estimated sorptivities values would be around 0.06 for clay and 0.13 for loam.These values, although roughly estimated, are of the same order of magnitude as those estimated based on simulations.Therefore the behavior observed in this research is once more confirmed.
Sorptivity is also defined as a measure of the capillary removal of water, being essentially a property of the medium.It is similar to the "capillary absorption coefficient" which is proportional to the square root of the capillary radius (Philip, 1957).Therefore, the smaller the capillary radius (i.e. the finer the soil particle size), the smaller the "capillary absorption coefficient" and thus, sorptivity.
Soil sorptivity is proportional to the pore diameter of the medium.The sorptivity of clay, despite a higher capillary rise potential, was not higher than in the loam material.The early-time infiltration rates, which are directly proportional to sorptivity, were also lower in clay.From these results and also from literature it is possible to conclude that the lower porosity of pure clays does not allow a higher infiltration rate, despite a higher capillary rise potential than in coarser materials.Also, the relationship between sorptivity and saturated hydraulic conductivity is not clear enough to define the real influence of K s and gravity on early-time infiltration.Clear quantitative and qualitative definitions on the significance of the term "early-time infiltration" are also needed to improve the interpretation of sorptivity experiments.

Fractal characterization of infiltration process
Sorptivity estimated with the time exponent parameter (n) fixed at 0.5 was not statistically different from sorptivity estimated with n also estimated in the nonlinear regression procedures (Figure 7).The variation of S estimated with varying n at each iteration level for clay and loam is illustrated in figure 8.There are no visual differences between figures 6 and 8.The use of n fixed at 0.5 or n estimated by regression procedures had no influence on the estimated sorptivity value.The estimated n values were therefore used in the fractal characterization of the infiltration process (Guerrini & Swartzendruber, 1994;Pachepsky & Timlin, 1998).
As described earlier, the estimated time exponent (n) was associated with the fractal parameter Hurst exponent (H) (Guerrini & Swartzendruber, 1994).The n = H values ranged from 0.43 to 0.49 ( 0.48; σ = 0.028) in clay and from 0.47 to 0.49 ( 0.49; σ = 0.008) in loam, where average and σ = standard deviation.The variation of the time exponent at each iteration level for clay and loam is shown in figure 9.The mean fractal dimension (D c ) value calculated from n = H values of clay and loam, using equation 4, was 0.52.The fractal dimension of the Cantor Bar used to generate the soil columns (D) was approximately 0.63.Therefore, the variation between D c and D was relatively small (around 17.5 %).
The Hurst exponent associated with the time exponent was less than the traditional 0.5 for all materials and iteration levels.The character of a Brownian motion is expressed and classified as either regular or fractional, according to the value of the Hurst exponent.The conventional diffusion is associated with a regular Brownian motion if H ≡ 0.5 and as a fractional Brownian motion if 0 < H ≠ 0.5 < 1  ( Guerrini & Swartzendruber, 1994).According to the results of this study, the Brownian motion process in "early-time" water infiltration in soil can be classified as fractional, since H < 0.5.

CONCLUSIONS
Hydrus 2D adequately simulated the "early-time" infiltration process for pure clay and pure loam and clay and loam materials mixed in different proportions following the horizontal distribution of a Cantor Bar with fractal dimension D ≈ 0.63.Soil sorptivity increased as the loam content increased in the columns.This behavior was attributed to the dependence of the sorptivity parameter on capillary radius.The time exponent parameter estimated by nonlinear regression was found to be less than the traditional value of 0.5.The fractal dimension estimated from the Hurst exponent was 17.5 % lower than the fractal dimension of the Cantor Bar used to generate the columns.Since the value of H was different from 0.5, the water diffusion process within the columns was characterized as fractional Brownian motion.

Figure
Figure 2. Textural triangle for classification of soil materials according to clay, silt and sand contents.

Figure
Figure 1.Cantor Bar with an initiator of 10 cm in length and fractal parameters D ≈ ≈ ≈ ≈ ≈ 0.6309 and b = 3.

Figure 3 .
Figure 3. Distribution of the materials within the columns according to the iterations of the Cantor Bar.In (a) CLAY0 (left) and LOAM0; in (b) CLAY1 (left) and LOAM1; in (c) CLAY2 (left) and LOAM2; in (d) CLAY3 (left) and LOAM3; in (e) CLAY4 (left) and LOAM4.

Figure 6 .
Figure 6.Sorptivity parameter estimated with time exponent fixed at 0.5 at each iteration level for clay and loam materials.

Figure 5 .
Figure 5. Infiltration rate as a function of time for loam material at each iteration level (clay content is zero at iteration zero and increases with each iteration).

Figure 4 .
Figure 4. Infiltration rate as a function of time in clay material at each iteration level (loam content is zero at iteration zero and increases with each iteration).

Figure 8 .
Figure 8. Sorptivity parameter estimated with the time exponent also estimated by regression at each iteration level for clay and loam materials.

Figure 1 :
Figure 1:1 plot of Sorptivity (S) estimated with time exponent (n) fixed at 0.5 and with time exponent also estimated by regression.

Figure 9 .
Figure 9. Variation of estimated time exponent (n) at each iteration level for clay and loam materials.