Nonlin. Processes Geophys., 16, 159–168, 2009 www.nonlin-processes-geophys.net/16/159/2009/ © Author(s) 2009. This work is distributed under the Creative Commons Attribution 3.0 License. Nonlinear Processes in Geophysics Characterizing water fingering phenomena in soils using magnetic resonance imaging and multifractal theory A. Posadas1,3 , R. Quiroz1 , A. Tannús2 , S. Crestana3 , and C. M. Vaz3 1 International Potato Center – CIP, P.O. Box 1558, Lima 12, Peru de Fı́sica e Informática, Instituto de Fı́sica de São Carlos, Universidade de São Paulo, Av. Dr. Carlos Botelho 1465, CEP 13560-250 São Carlos – SP, Brazil 3 Embrapa Instrumentação Agropecuária, Rua XV de Novembro, 1452, São Carlos, SP – CEP 13560-970, Brazil 2 Departamento Received: 1 September 2008 – Revised: 4 February 2009 – Accepted: 4 February 2009 – Published: 26 February 2009 Abstract. The study of water movement in soils is of fundamental importance in hydrologic science. It is generally accepted that in most soils, water and solutes flow through unsaturated zones via preferential paths or fingers. This paper combines magnetic resonance imaging (MRI) with both fractal and multifractal theory to characterize preferential flow in three dimensions. A cubic double-layer column filled with fine and coarse textured sand was placed into a 500 gauss MRI system. Water infiltration through the column (0.15×0.15×0.15 m3 ) was recorded in steady state conditions. Twelve sections with a voxel volume of 0.1×0.1×10 mm3 each were obtained and characterized using fractal and multifractal theory. The MRI system provided a detailed description of the preferential flow under steady state conditions and was also useful in understanding the dynamics of the formation of the fingers. The f (α) multifractal spectrum was very sensitive to the variation encountered at each horizontally-oriented slice of the column and provided a suitable characterization of the dynamics of the process identifying four spatial domains. In conclusion, MRI and fractal and multifractal analysis were able to characterize and describe the preferential flow process in soils. Used together, the two methods provide a good alternative to study flow transport phenomena in soils and in porous media. 1 Introduction The fluid flow through preferential paths or fingers is extremely important in hydrological and agricultural processes such as infiltration of water and the transport of agrochemiCorrespondence to: A. Posadas (a.posadas@cgiar.org) cals through the soil profile. Preferential paths increase the probability of underground water contamination – through a faster transportation of pesticides, heavy metals, radioactive waste and other contaminants – and thus constitute a phenomenon of particular interest. To understand this complex problem, some researchers concentrate their efforts in studying fluid transport in connection to the geometry of porous media (Lu et al., 1994). In a review of the principal theories of the fingering phenomena (Steenhuis et al., 1996), a model using a modification of the invasion percolation model (Glass and Yarrington, 1996), based on both laboratory and field studies, was introduced. Many experiments have shown that the fluid transport-porous medium coupling has auto-similarity or fractal characteristics in a range of defined scales (Katz and Thompson, 1985). The fingering phenomenon in soils, which is controlled by both capillarity and gravity force (Hill and Parlange, 1972; Glass et al., 1988, 1989; Crestana and Posadas, 1998), also presents fractal characteristics (Chang et al., 1994; Posadas and Crestana, 1993; Crestana and Posadas, 1998). These fractal characteristics have led to the creation of simulation models such as the Diffusion Limited Aggregation-DLA (Chen and Wilkinson, 1985) that simulates the viscous fingering phenomenon and the invasion percolation model (Wilkinson and Willemsen, 1983), which in turn simulates the capillary fingering phenomenon. Other researchers have confirmed the suitability of fractal and multifractal theory to describe and simulate preferential flow. Ogawa et al. (2002) applied fractal analysis to study preferential flow on field soils and obtained a good correlation between the surface fractal dimension and the exponent of a Van Genuchten expression applied to the particle size distribution of the soil. A multifractal analysis was successfully employed by Nittmann et al. (1987) and Mâløy et al. (1987) on viscous fingering structures observed in HeleShaw cells and in a mono-layer of glass beads, respectively. Published by Copernicus Publications on behalf of the European Geosciences Union and the American Geophysical Union. 160 A. Posadas et al.: Characterizing water fingering phenomena in soils Recent publications suggest that the multifractal formalism is applicable to three-dimensional systems. For example, it has been shown (Held and Illangasekare, 1995) that the width (internal energy or 1α) of the f (α) curve (multifractal spectrum), in the range of positive moments, quantifies displacement instability. On the other hand, a basic technical challenge in investigations of mass transport in soils is the need for a quantitative, visible and nondestructive monitoring of spatial and temporal water distribution, as a prelude to its quantitative analysis. Magnetic resonance imaging (MRI) is becoming an important tool for studies of patterns and mechanisms of water infiltration into soils (Amin et al., 1998; Posadas et al., 1996). This paper proposes an innovation in the characterization of preferential flow by combining MRI (Posadas et al., 1996; Crestana and Posadas, 1998) with multifractal theory (Chhabra et al., 1989b; Posadas et al., 2001, 2003) for a three dimensional description of the dynamics of fingers in sandy soils. 2 2.1 Materials and methods Experiment In order to visualize and characterize the fingering phenomenon in three dimensions, magnetic resonance imaging (MRI) and fractal and multifractal theory were employed in steady state conditions (Posadas, 1994; Onody et al., 1995; Posadas et al., 1996). The work was conducted at CNPDIA/EMBRAPA Laboratory, Brazil. A cubic 0.15×0.15×0.15 m3 double-layer quartz sand column was built, with the aim to generate preferential flow as described by Glass et al. (1989). The top layer was 2.0 cm high, consisting of fine texture quartz sand (particle size diameter in mm: 0.106<d<0.149), with a porosity (φ) of 0.448±0.030 m3 m−3 , saturated hydraulic conductivity (Ks ) of (6.3±0.08) 10−5 m s−1 and bulk density (ρb ) of 1500±80 kg m−3 and the lower 0.12 m layer was filled with coarse sand (particle size diameter in mm: 0.212<d<0.500), with a porosity (φ) of 0.312±0.030 m3 m−3 , saturated hydraulic conductivity (Ks ) of (26.7±0.15) 10−4 m s−1 and bulk density (ρb ) of 1800±100 kg m−3 . On the surface of the top layer, an acrylic plate of 0.15×0.15×0.03 m3 with 1.0 mm diameter holes was placed, in order to spread water uniformly over the surface. The top 2 mm section of the column remained free for water application. This column was placed into the head coils of the 500 gauss MRI systems (Fig. 1). The infiltration of water through the cubic column was studied under hydrodynamic steady state conditions. Porosity and saturated hydraulic conductivity for each soil type – fine and coarse – were previously estimated on other columns filled with either fine sand soil or coarse sand. The gravimetric method with water saturation and the constant head soil core (tank) method (Dane and Topp, 2002) with 10 replicates, were used respectively. All these parameters were Nonlin. Processes Geophys., 16, 159–168, 2009 determined before running the MRI experiments. The cubic column was first filled with the coarse material carefully placed at the bottom to achieve a homogeneous distribution. A single layer of filter paper was then placed to separate the soil types. A similar procedure was followed to place the fine layer on the top of the column. The process was repeated until a homogeneous substrate and the bulk density (ρb ), required for the desired saturated hydraulic conductivity were approximated. A second layer of filter paper was placed to separate the fine sand from the acrylic plate. A constant positive pressure of 2 mm of water was applied on top of the column by watering the system through an inlet placed above the acrylic plate, with 1 mm diameter holes. This system was coupled to a medical infusion pump hanged on a IV-stand equipment (Fig. 1a). This pump delivered a constant water flow and pressure and facilitated the halting of water flowing into the column. A water volume of 600 ml was applied during 3 min to the top of the cubic column at a constant rate of 200 ml min−1 . The steady-state flow of water was established when the first finger reached the bottom of the box. At that moment the water inflow to the column was turned off. The steady-state condition was demonstrated through NMR Spin-Echo experiments using the methodology described by Crestana and Posadas (1998); Posadas (1994); Posadas et al. (1996). Different 2-D and 3D experiments were conducted to verify the diffusion process for water redistribution, when the steady-state was reached. It has been observed previously that even after 24 h the fingers formed initially remained virtually unchanged (Crestana and Posadas, 1998). Horizontally-oriented, transverse and sagittal images of the column were recorded during 11.2 s just after the steady-state conditions were achieved as described by Posadas et al. (1996) and Crestana and Posadas (1998). All the recorded MRI images were obtained in 2-D slices and then processed using the software Image Processing and Analysis in Java-ImageJ (Rasband, 2007) and FracLab (Levy-Vehel and Mignot, 1994) developed by INRIA (2005). The contrast in each slice of grey images was enhanced with different tools from ImageJ and FracLab which permitted a good binarization using a thresholding filtering in the range 62–72 of the grey level (Fig. 2). These 2-D binary images were utilized for fractal and multifractal analysis using box-counting method following the gravitational force (horizontally-oriented slices). The binarization error was estimated by calculating the total area (160×160 pixels) of the first top slice and subtracting the binarized wet-area, which was estimated in 7.0% of error (see Table 1). 2.2 2.2.1 Theory MRI basics Subatomic particles such as protons have the quantum mechanical property of spin. Certain nuclei such as 1 H (protons) have a non-zero spin and therefore a magnetic moment. www.nonlin-processes-geophys.net/16/159/2009/ A. Posadas et al.: Characterizing water fingering phenomena in soils 161 Fig. 1. (a) 500 gauss MRI system (IFSC-USP Laboratory, Brazil) showing the cubic column used (b) Sketch of the cubic column. When these spins are placed in a strong external magnetic field they precess synchronously at a unique frequency, the Larmor frequency, which is characteristic of the species of atom and the strength of the magnetic field. Changes of this microscopic magnetization induced by a synchronous radiofrequency (RF) pulse result in a large proportion of the nuclei being excited to a higher energy state. When the RF pulses are discontinued the nuclei decay back to equilibrium. During that time, the energy dissipates from the excited protons into their environment (spin-lattice or spin-spin relaxations) and that energy is detected as an electrical signal by the RF receiver coil that has been previously tuned to detect the radiation at the Larmor frequency of protons associated with water. These signals are localized in space by the field gradient coils installed in the magnet. Therefore, by adjusting the strength of the magnetic field with gradient coils, the spatial differences in source of signals needed to produce an image are obtained. Images are reconstructed from the signals through the use of a Fourier transform (Mansfield and Morris, 1982; Shaw, 1984; Brown et al., 1998; Van As, 2007). 2.2.2 Fractal and multifractal analysis It is now widely accepted that physical systems that exhibit chaotic behavior are generic in nature. Since these systems lose information exponentially fast it is possible to follow and predict their motion in any detail only for short time scales (Chhabra et al., 1989b). To describe their longterm dynamic behavior, one must resort to suitable statistical descriptions. One such description is multifractal formalism (Chhabra et al., 1989b; Hentschel and Procaccia, 1983; Halsey et al., 1986; Chhabra and Jensen, 1989a). Multifractal theory permits the characterization of complex phenomena in a fully quantitative manner, for both temporal and spatial variations. Multifractal techniques and notions are increasingly widely recognized as the most appropriate and straightwww.nonlin-processes-geophys.net/16/159/2009/ Fig. 2. Image binarization processes. (a) Thresholding filtering in the range 62–72 of the grey level; (b) binarized image. forward framework within which to analyze and simulate not only the scale dependency of the geophysical observations, but also their variability over a wide range of scales (Mandelbrot, 1982; Schertzer and Lovejoy, 1996). The basic equation of the fractal theory expresses the relationship between the number and the size of the objects (Feder, 1988): N (ε) ∼ ε−D0 , (1) where N (ε) is the number of objects, ε is the scale and D0 is the fractal or capacity dimension. In this paper we will use the capacity dimension name, following the convention of Beck and Schlogl (1995). The box-counting technique is used to estimate the scaling properties of an image by covering it with boxes of size ε and counting the number of boxes containing at least one pixel representing the object under study: D0 = − lim ε→0 log N (ε) . log(ε) (2) Nonlin. Processes Geophys., 16, 159–168, 2009 162 A. Posadas et al.: Characterizing water fingering phenomena in soils Table 1. Selected nparameters of preferential flow infiltration and multifractal parameters. MRI slices numbers – 12 11 10 09 08 07 06 05 04 03 02 01 Depth (cm) ±error r 0.00 2.2±0.1 3.3±0.1 4.4±0.1 5.5±0.1 6.6±0.1 7.7±0.1 8.8±0.1 9.9±0.1 11.0±0.1 12.1±0.1 13.2±0.1 14.3±0.1 Wet-area percent ±error b Capacity dimension D0 (q=0)±MSE R 2 moments variation 1q – 100.00±7 77.40±7 53.42±7 44.64±7 36.55±7 35.20±7 30.60±7 25.20±7 22.70±7 21.30±7 19.30±7 14.80±7 – 2.00±0.01 1.94±0.02 1.81±0.05 1.79±0.04 1.76±0.05 1.73±0.04 1.71±0.04 1.59±0.05 1.60±0.02 1.44±0.03 1.39±0.02 1.34±0.02 – 1.0000 0.9995 0.9949 0.9944 0.9956 0.9933 0.9872 0.9828 0.9916 0.9748 0.9467 0.9985 – 14.20 11.20 11.76 6.56 4.44 3.26 2.92 2.22 1.22 0.74 0.72 0.28 1α±MSE domains Spatial 0.00 0.2098±0.04 1.3602±0.14 1.5187±0.16 1.4429±0.16 1.4784±0.13 1.1053±0.12 0.3437±0.03 0.4919±0.05 0.1767±0.05 0.0536±0.005 0.0774±0.008 0.0227±0.002 HZ LSH LSH-IP LSH-IP LSH-IP LSH-IP WSH WSH WSH EZ-Fractal EZ-Fractal EZ-Fractal Note: all the probability values (p-value) computed were, p<0.0001 error r is the reading error introduced for the measure rule; error b is the binarization error; MSE is defined as the mean square error; R 2 is the correlation coefficient; 1q is defined as qmax −qmin ; 1α is the internal energy variation defined as αmax −αmin ; HZ is the homogeneous zone; LSH is large spatial heterogeneity; IP is invasion percolation; WSH is weak spatial heterogeneity; EZ is the equilibrium zone. Provided the limit exists, the infinitum of N (ε) is approximated by varying the origin of the grid until the smallest number is found. Using Eq. (2), the capacity dimension D0 can be determined as the negative slope of log N (ε) versus log(ε), measured over a range of box widths. In a homogeneous system, the probability (P ) of a measured quantity (measure) varies with scale ε as (Chhabra et al., 1989b; Evertsz and Mandelbrot, 1992; Vicsek, 1992): P (ε) ∼ εD , (3) where D is a fractal dimension. For heterogeneous or nonuniform systems the probability within the ith region Pi varies as: αi Pi (ε) ∼ ε , (4) where f (α) can be considered as the generalized fractal dimension of the set of boxes with singularities α (Kohmoto, 1988). The exponent α can take on values from the interval [α−∞ , α+∞ ], and f (α) is usually a function with a single maximum at df (α(q))/dα(q)=0 (where q is the order moment of a statistic distribution). Thus, when q=0, fmax is equal to the capacity dimension, D0 (Gouyet, 1996; Vicsek, 1992). Multifractal sets can also be characterized on the basis of the generalized dimensions of the qth order moment of a distribution, Dq , defined as (Hentschel and Procaccia, 1983): Dq = lim ( ε→0 1 log µ(q, ε) ), q − 1 log(ε) (6) where µ(q, ε) is the partition function (Chhabra et al., 1989b): where αi is the Lipschitz-Hölder exponent or singularity strength, characterizing scaling in the ith region or spatial location (Feder, 1988). The parameter αi quantifies the degree of regularity in point xi . Loosely speaking, any measure µ of an interval [xi , xi +1x], behaves as (1x)αi (Halsey et al., 1986). For a uniform distribution one finds αi (x)=1 for all x. More generally, for any real value a>0 the distribution with density x a−1 on [0, 1] has αi (0)=a and αi (x)=1 for all xǫ[0, 1]. Values αi (x)<1 indicate, thus, a burst of the event around x on all levels, while αi (x)>1 is found in regions where events occur sparsely (Riedi, 1999). Similar αi values might be found at different positions in the space. The number of boxes N (α) where the probability Pi has singularity strengths between α and α+dα is found to scale as (Chhabra et al., 1989b; Halsey et al., 1986): where τ (q) is the correlation exponent of the qth order moment defined as (Halsey et al., 1986; Vicsek, 1992): N(α) ∼ ε−f (α) , τ (q) = (q − 1)Dq , Nonlin. Processes Geophys., 16, 159–168, 2009 (5) µ(q, ε) = N(ε) X q Pi (ε), (7) i=1 The generalized dimension Dq is a monotonic decreasing function for all real q’s within the interval [−∞+∞]. When q<0, µ emphasizes regions in the distribution with less concentration of a measure, whereas the opposite is true for q>0 (Chhabra and Jensen, 1989a). Also, the partition function scales as: µ(q, ε) ∼ ετ (q) , (8) (9) www.nonlin-processes-geophys.net/16/159/2009/ A. Posadas et al.: Characterizing water fingering phenomena in soils Fig. 3. Images obtained with the MRI system, showing three vertically-oriented sections of the fingering phenomenon in static conditions. Each section represents a slice, 2 cm thick, 15 cm wide and 15 cm high, of the 15×15×15 cm cubic soil column. The white areas represent wet zones and the grey areas represent dry zones. The connection between the power exponents f (α) (Eq. 5) and τ (q) (Eq. 9) is made via the Legendre transformation (Callen, 1985; Chhabra and Jensen, 1989a; Halsey et al., 1986): f (α(q)) = qα(q) − τ (q) and α(q) = dτ (q) dq (10) f (α) is a concave downward function with a maximum at q=0. When q takes the values of q=0, 1 or 2, (Eq. 6) is reduced to: log(N(ε)) D0 = lim , D1 = lim ε→0 log(ε) ε→0 N(ε) X µi (ε) log(µi (ε)) i=1 log(ε) log(C(ε)) , ε→0 log(ε) D2 = lim Fig. 4. Images obtained by MRI showing twelve horizontaloriented sections (slices 1.0 mm stick) of the cubic soil column following the direction of the infiltration. (a) Saturated section near the surface of the column (first layer); (l) section corresponding to the bottom of the soil column; and (b) through (k) represent intermediate situations. Light grey areas correspond to invading water within the pore networks and black areas correspond to both the non-invaded porous medium and quartz solid phase. Grey tones are function of the sand water content. (or average) information about a system (Voss, 1988). The entropy dimension is related to the information (or Shannon) entropy (Shannon and Weaver, 1949). The correlation dimension D2 is mathematically associated with the correlation function (Grassberger and Procaccia, 1983) and computes the correlation of measures contained in a box of size ε. The relationship between D0 ,D1 ,and D2 is, D2 ≤ D1 ≤ D0 , , (11) respectively, with C(ε) being the correlation function. The values D0 , D1 and D2 are known as the capacity dimension, the entropy dimension and the correlation dimension, respectively. The capacity dimension provides global www.nonlin-processes-geophys.net/16/159/2009/ 163 (12) The equality D0 =D1 =D2 occurs only if the fractal is statistically or exactly self-similar and homogeneous (Korvin, 1992). Following the methodology used by Posadas et al. (2001, 2003), multifractal theory was applied to the MR images of the fingering phenomena. The size of each 2-D image considered for the multifractal analysis was 160×160 pixels (or 15×15 cm2 ). Twelve binary-images slices of Nonlin. Processes Geophys., 16, 159–168, 2009 164 A. Posadas et al.: Characterizing water fingering phenomena in soils 2.5 y = -0.0503x + 2.0957 2 R = 0.966 f(q=0) 2 1.5 1 0.5 0 0 2 4 6 8 10 12 14 16 Depth (cm) Fig. 6. Capacity dimension, of the fingers, measured along the vertical depth. Fig. 5. Box counting plots for each binary image slices analyzed for q=0 with upper and lower boxes sizes. These plots show the scale range of fractal behavior of the systems. Also, two linear regressions estimating the capacity dimension for bottom and top layer slices are shown. horizontally-oriented sections, following the gravitational force, were analyzed using the multifractal algorithm (CIPMASS -downloadable after subscription at http://inrm.cip. cgiar.org/vlab). This algorithm was developed based on the method described in Chhabra and Jensen (1989a) and implemented in Matlab by Posadas et al. (2001). The spatial distribution of water concentration through each slice was partitioned in boxes size L, for L=8, 10, 16, 20, 32, 40 and 80 pixels (see Fig. 5). These upper and lower box sizes considered prevent the systematic biases in the small and large scales as mentioned by Vignes-Adler et al. (1991). The capacity dimension D0 was obtained from the maximum value of the multifractal spectrum when q=0 (Beck and Schlogl, 1995). The moments (q) in each slice considered ranged from 1q=14.20 to 0.28, as shown in Table 1, with steps variation of 0.02 which is fine enough to show the multifractal behavior in the very narrow range of q’s. 3 3.1 Results and discussion Magnetic resonance image analysis The results obtained with the MRI system are shown in Table 1 and depicted in Figs. 3 and 4. Figure 3 shows three acquisitions of the transverse plane of the cubic sand column at steady-state flow. In these three images it is possible to observe the three-dimensional character of the fingering pheNonlin. Processes Geophys., 16, 159–168, 2009 nomena and its spatial variability. After a few seconds of turning off water inflow, it was observed that the formation of fingers was completely halted. These findings were consistent across the five experiments conducted using different packed columns, but keeping the same physical parameters. Similar finger occupation behavior and steady-state conditions were observed. As expected the only thing that varied was the zoned occupied by the fingers. The accumulation of water at the bottom of the column can be seen in panel A. Since there is no evidence of water flow except through the preferential paths it can be inferred that, if the experiments replicate the process occurring in real systems, the movement of water and the solutes conveyed within it could reach shallow or deep water pools in the soil faster than in conditions where preferential flow is absent. The column under steady-state flow was “dissected” into 12 horizontally-oriented slices of 1.0 cm thick, as shown in Fig. 4. The light grey area depicts the concentration of water in the cross section. It can be seen that the water concentration throughout the profile follows a fixed path. Dark grey paths in the images represent dry sand areas. These findings support the existence of preferential flow under the conditions studied (Posadas et al., 1996). 3.2 Fractal and multifractal analysis Figures 5, 6, 7, 8 and Table 1 summarize the fractal and multifractal analysis. The capacity dimension (D0 , Fig. 6) evidences differences in the spatial distribution among the twelve layers of the columns depicted in Fig. 4. When D0 is near 2.0 (top layer), the system is more homogeneous i.e. most of the pores are filled with water (saturation of porous media), and corresponds to an apparent single water phase. The infiltration process throughout the entire top layer was slow and uniform with a constant vertical velocity like a plane wave (panel a in Fig. 4), following the gravitational force. As soon as the water crosses the interface between the fine and the coarse layers, a hydrodynamic instability seems to dominate de infiltration process and the fingers began to www.nonlin-processes-geophys.net/16/159/2009/ A. Posadas et al.: Characterizing water fingering phenomena in soils 165 MRI12 3 f(α) 2 1 0 0 1 2 3 α 4 MRI11 3 f(α) 2 1 0 0 1 2 3 α 4 Fig. 8. Variation of the capacity dimension and the internal energy of the system as a function of the depth along the cubic column. MRI05 3 f(α) 2 1 0 0 1 2 3 α 4 MRI01 3 f(α) 2 1 0 0 1 2 α 3 4 Fig. 7. Magnetic Resonance Images from a three-dimensional fingering phenomena and its multifractal spectrums following the gravitational direction (from top MRI12 to bottom MRI01). These sets represent binary images, where the black areas correspond to water and the white one correspond to the dry area. develop. This is depicted in the panel b in Fig. 4. This corresponds to a section with high heterogeneity in the distribution of water, which decreases towards the bottom of the column. This seems to be a good descriptor of the steady state conditions at which the images were taken. Figure 7 shows the multifractal analysis used to describe the heterogeneity of the spatial variability of the water in each cross-section throughout the profile of the column. A quick inspection along the column indicates how dynamic the system is, even though the analysis was made when the system reached a steady-state flow; a condition reached after 137 s of infiltration. It goes from the wetting instability – the first condition of the coarse texture substrate (Panel b in Fig. 4) to a hydrodynamic stability at the bottom (Panel l in Fig. 4), passing through chaotic sections in the intermediate portion (Panels c through k in Fig. 4). www.nonlin-processes-geophys.net/16/159/2009/ The spatial variations in all these conditions seem to be well characterized with the f (α) α spectrum as shown on the right panels in Fig. 7. For instance, the cross-section labeled MRI11 represents the wetting instability conditions in the column. This hydrodynamic instability is attributed to the change produced when the water flows from the fine to the coarse texture. The multifractal spectrum is characteristic of a heterogeneous system with variations on both sides of the maximum value. From the maximal capacity dimension (D0 ) to the left, the spectrum describes the behavior of the areas where water is present (positive q’s). The asymmetry toward the right from α=2 indicates domination of small or extremely small values of water. This is an indication of the existence of preference paths. The condition prevails, with small variations, down through the cross-section MRI07. Inspection of the variation in the internal energy (1α) shows that the 1α value for MRI11 is 1.36 (Table 1). Compared to the saturation condition presented in the layer with fine texture (1α=0.21) one can contrast how this multifractal parameter changes when the condition changes from a homogeneous (MRI12) to a heterogeneous wetting instability. The cross-sections MRI05 and MRI04 seem to be a transition section from instability to hydrodynamic stability. The corresponding 1α is around 0.4. The bottom three crosssections correspond to hydrodynamic stability, behaving as a pure fractal with a capacity dimension D0 changing from 1.44 to 1.34 and 1α<0.07 3.3 Spatial domains Different spatial domains can be identified following the depth of the soil column and the multifractal parameters. These domains are summarized in Table 1 and the spatial behavior of each slice can be seen in Fig. 4. The spatial domains were labeled as HZ (homogeneous zone), Nonlin. Processes Geophys., 16, 159–168, 2009 166 A. Posadas et al.: Characterizing water fingering phenomena in soils LSH (large spatial heterogeneity), WSH (weak spatial heterogeneity) and EZ (equilibrium zone). The equilibrium zone is attributed to the fact that multifractality does not occur and it behaves as a pure fractal. As an example let us analyze the slices MRI10 to MRI07, corresponding to the onset of observed fingers, with a rather constant capacity dimension, D0 of 1.73±0.10−1.81±0.10 (see also panels c through f in Fig. 4). The average D0 ∼ =1.82, might be associated with the invasion percolation capacity dimension without trapping (Gouyet, 1996). In the light of the invasion percolation theory (Wilkinson and Willemsen, 1983; Onody et al., 1995), this seems to be an unstable zone characteristic of a percolative pore network. This pore networks feature, from homogeneous to heterogeneous behavior, suggests that the gradient invasion percolation might be a good model to simulate preferential-infiltration processes in soils (Frette et al., 1992). An interesting feature of the multifractal analysis is that in spite of the fact that the image was taken under steadystate flow, the dynamic characteristics present in the profile is described by the f (α) α spectrum and its parameters. Figure 8 captures these characteristics by showing the dynamic of D0 and 1α across the column profile. Following the column depth, there is a reduction in the capacity dimension D0 , represented with the grey level. This can be seen as a first and gross approximation of the fingers’ dynamic. The inclusion of the internal energy 1α or the “ensemble average” of the total energy of the system (Reichl, 1998) – visualized through a proportional thickness of the line – enriches the description of the phenomena. It is evident that as soon as the water reaches the coarse layer there is a multifractal behavior associated to a spatial domain of high heterogeneity, which decays quasi exponentially as a function of the column depth until it dies down to an equilibrium zone with fractal behavior. The four spatial domains described above – HZ, LSH, WSH and EZ – are clearly differentiated when these two multifractal parameters are jointly interpreted. In terms of the magnitude, a stable zone can be inferred when D0 is around 2 and 1α∼0−0.2. The maximum instability zone from 4.4 to 7.7 cm of depth, corresponds to the invasion – percolation zone with D0 ∼1.73–1.81 and 1α∼1–1.5. The transition zone from 9 to 11 cm of depth presented D0 ∼1.59– 1.71 and 1α∼0.5–0.2. In turn, the equilibrium zone, from 12 to 14 cm of depth, showed D0 ∼1.34–1.44 and 1α∼0.02– 0.05; depth at which the hydrodynamic stability has been achieved. 4 Conclusions 1. Both the multifractal theory and the MRI system were good tools to characterize preferential water flow in soils. The MRI system is useful in visualizing the phenomenon for a better understanding of the process. The technique presented very innovative and encouraging results for the study Nonlin. Processes Geophys., 16, 159–168, 2009 of preferential flow in soils and porous media, allowing the observation, in a non-disturbing way, of the 3-D and random character of the phenomenon. Nonetheless, very few soil or geophysical research groups have access to this type of equipment. 2. The use of multifractal theory facilitates describing the dynamics of preferential flow and can be used to predict the outcomes under real conditions and to improve the accuracy of existing models, as the invasion percolation model, for example. The combination of these techniques opens a new set of options that must be further tested for different soil types and management conditions. 3. Also, using the multifractal parameters, capacity dimension D0 and the internal energy 1α it was possible to identify four spatial domains along the soil columns studied. These novel results open new research alternatives to study infiltration processes in soils and thus require further research. This is a challenge to be faced in the near future. Acknowledgements. The authors are indebted to several institutions and people. The final data analysis and the production of the manuscript were supported by the CIDA-Canada through the ALTAGRO Project and FAPESP (Brazilian Agency) Process No. 2007/58561-7. The authors would like to thank Victor Mares for the revision and improvement of the manuscripts and Ivonne Valdizan for patiently going through the process of learning Latex and converting the manuscript into that format, in addition to her normal contribution in guaranteeing the consistency and adequacy of the papers to the journal’s format. Edited by: A. Provenzale Reviewed by: an anonymous referee References Amin, M. H. G., Hall, L. D., Chorley, R. J., and Richards, K. S.: Infiltration into soils, with particular reference to its visualization and measurement by magnetic resonance imaging(MRI), Prog. Phys. Geog., 22(2), 135–165, 1998. Beck, C. and Schlogl, F.: Thermodynamics of Chaotic Systems: An introduction. Cambridge Nonlinear Science Series 4, Cambridge University Press, 306 pp., 1995. Brown, J. M., Kramer, P. J., Cofer, G. P., and Johnson, G. A.: Use of nuclear magnetic resonance microscopy for noninvasive observations of root-soil water relations, Theor. Appl. Climatol., 42(4), 229–236, 1998. Callen, H. B.: Thermodynamics and an introduction to thermostatics, John Wiley and Sons, New York, 2nd edn., 512 pp., 1985. Chang, W.-L., Biggar, J. W., and Nielsen, D. R.: Fractal description of wetting front instability in layered soils, Water Resour. Res., 30(1), 125–132, 1994. Chen, J. D. and Wilkinson, D.: Pore-scale viscous fingering in porous media, Phys. Rev. Lett., 58(18), 1892–1895, 1985. www.nonlin-processes-geophys.net/16/159/2009/ A. Posadas et al.: Characterizing water fingering phenomena in soils Chhabra, A. B. and Jensen, R. V.: Direct determination of the f (α) singularity spectrum, Phys. Rev. Lett., 62(12), 1327–1330, 1989a. Chhabra, A. B., Meneveu, C., Jensen, R. V., and Sreenivasan, K. R.: Direct determination of the f (α) singularity spectrum and its application to fully developed turbulence, Phys. Rev. A, 40(9), 5284–5294, 1989b. Crestana, S. and Posadas, D. A. N.: 2-D and 3-D fingering phenomenon in unsaturated soils investigated by fractal analysis, invasion percolation modeling and non-destructive image processing, in: Fractals in Soils Science, edited by: Baveye, P., Parlange, J.Y., and Stewart, B. A., Boca Raton, Florida, USA CRC Press, 293–332, 1998. Dane, J. H. and Topp, G. C. (Eds.): Methods of Soil Analysis: Part 4, Physical Methods, SSSA Book Series No 5, Soil Science Society of America, Madison, WI, 2002. Evertsz, C. J. G. and Mandelbrot, B. B.: Multifractal measures, in: Chaos and Fractals, edited by: Peitgen, H.-O., Jurgens, H., and Saupe, D., New Frontiers of Science, Springer-Verlag, New York, 921–953, 1992. Feder, J.: Fractals, Plenum Press, New York, 283 pp., 1988. Frette, V., Feder, J., Målfy, K. F., Jfssang, T., and Meaking., P.: Buoyancy-driven fluid migration in porous media, Phys. Rev. Lett., 68(21), 3164–3167, 1992. Glass, R. J. and Yarrington, L.: Simulation of gravity fingering in porous media using a modified invasion percolation model, Geoderma, 70(2–4), 231–252, 1996. Glass, R. J., Steenhuis, T. S., and Parlange, J.-Y.: Wetting front instability as a rapid and far-reaching hydrologic process in the vadoze zone, J. Contam. Hydrol, 3, 207–226, 1988. Glass, R. J., Steenhuis, T. S., and Parlange, J.-Y.: Wetting front instability: experimental determination of relationships between system parameters and two-dimensional unstable flow field behavior in initially dry porous media, Water Res. Res., 25(6), 1195–1207, 1989. Gouyet, J. F.: Physics and fractals structure, Springer, New York, 234 pp., 1996. Grassberger, P. and Procaccia, I.: Characterization of strange attractors, Phys. Rev. Lett., 50(5), 346–349, 1983. Halsey, T. C., Jensen, M. H., Kadanoff, L. P., Procaccia, I., and Shraiman, B. I.: Fractal measures and their singularities: The characterization of strange sets, Phys. Rev. A, 33(2), 1141–1151, 1986. Held, R. J. and Illangasekare, T. H.: Fingering of dense nonaqueous phase liquids in porous media 2. Analysis and classification, Water. Resour. Res., 31(5), 1223–1231, 1995. Hentschel, H. G. E. and Procaccia, I.: The infinite number of generalized dimensions of fractals and strange attractors, Physica D: Nonlinear Phenomena, 8(3), 435–444, 1983. Hill, D. E. and Parlange, J.-Y.: Wetting front instability in layered soils, Soil Sci. Soc. Am. J., 36(5), 697–702, 1972. INRIA: FracLab: APIS, http://complex.futurs.inria.fr/FracLab/ index.html, 2005. Katz, A. J. and Thompson, A. H.: Fractal sandstone pores: Implications for conductivity and pore formation, Phys. Rev. Lett., 54(12), 1325–1328, 1985. Kohmoto, M.: Entropy function for multifractals, Phys. Rev. A, 37(4), 1345–1350, 1988. www.nonlin-processes-geophys.net/16/159/2009/ 167 Korvin, G.: Fractals Models in the Earth Sciences, Elsevier, Amsterdam, The Netherlands, 396 pp., 1992. Levy-Vehel, J. and Mignot, P.: Multifractal segmentation of images, Fractals, 2, 371–377, 1994. Lu, T. X., Biggar, J. W., and Nielsen, D. R.: Water movement in glass bead porous media 2. Experiments of infiltration and finger flow, Water Resour. Res., 30(12), 3283–3290, 1994. Mâløy, K. J., Boger, F., Feder, J., and Jøssang, T.: Dynamics and structure of viscous fingering in porous media, in: timedependent effects in disordered materials, edited by: Pynn, R. and Riste, T., Plenum Press, NY, 111–118, 1987. Mandelbrot, B. B.: The fractal geometry of nature, W.H. Freeman and Company, NY, 2nd edn., 468 pp., 1982. Mansfield, P. and Morris, P. G.: NMR Imaging in biomedicine, Supplement 2, Advances in Magnetic Resonance, Academic Press, New York, 354 pp., 1982. Nittmann, J., Stanley, H. E., Toubul, E., and Daccord, G.: Experimental evidence of multifractality, Phys. Rev. Lett., 58(6), 619– 622, 1987. Ogawa, S., Baveye, P., Parlange, J.-Y., and Steenhuis, T.: Preferential flow in the field soils, Forma, 17(1), 31–53, 2002. Onody, R. N., Posadas, D. A. N., and Crestana, S.: Experimental studies of fingering phenomena in two dimensions and simulation using a modified invasion percolation model, J. Appl. Phys., 78(5), 2970–2976, 1995. Posadas, D. A. N. and Crestana, S.: Aplicação da teoria fractal na caracterização do fenômeno “fingering” em solos, Revista Brasileira de Ciências Sociais, 17(1), 1–8, 1993. Posadas, D. A. N.: Estudo do fenômeno “fingering” em um meio poroso através de imagens e teoria da percolação por invasão, Ph.D. thesis, IFSC/USP, São Carlos, SP-Brazil, 187 pp., 1994. Posadas, D. A. N., Tanns, A., Panepucci, C. H., and Crestana, S.: Magnetic resonance imaging as a non-invasive technique for investigating 3-D preferential flow occurring within stratified soil samples, Computers and Electronics in Agriculture, 14(4), 255– 267, 1996. Posadas, D. A. N., Gimenez, D., Bittelli, M., Vaz, C. M. P., and Flury, M.: Multifractal Characterization of soil particle-size distributions, Soil Sci. Soc. Am. J., 65(5), 1361–1367, 2001. Posadas, D. A. N., Gimenez, D., Quiroz, R., and Protz, R.: Multifractal characterization of soil pore systems, Soil Sci. Soc. Am. J., 67(5),1361–1369, 2003. Rasband, W.: ImageJ 1.38x, National Institute of Health, USA, available at: http://rsb.info.nih.gov/ij/, 2007. Reichl, L. E.: A modern course in statistical physics, John Wiley and Sons, Inc. NY, 2nd edn., 580 pp., 1998. Riedi, R. H.: An Introduction to Multifractals, Technical Report, Rice University, Unpublished Version, 5 September 1999. Shaw, D.: Fourier transform N.M.R. spectroscopy, Studies in physical and theoretical chemistry, NY Elsevier Publishing Company, 2nd edn., 304 pp., 1984. Schertzer, D. and Lovejoy, S.: EGS Richardson AGU Chapman NVAG3 Conference: Nonlinear variability in geophysics: scaling and multifractal processes, Nonlinear Proc. Geoph., 1(2–3), 77–79, 1996. Shannon, C. E. and Weaver, W.: The mathematical theory of communication, Urbana, University of Illinois Press, 125 pp., 1949. Nonlin. Processes Geophys., 16, 159–168, 2009 168 A. Posadas et al.: Characterizing water fingering phenomena in soils Steenhuis, T. S., Ritsema, C. J., and Dekker, L. W. (Eds): Fingered flow in unsaturated soil: from nature to model, Geoderma, Special Issue, 70(2–4), 83–326, 1996. Van As, H.: Intact plant MRI for the study of cell water relations, membrane permeability, cell-to-cell and long distance water transport, J. Exp. Bot., 58(4), 743–756, 2007. Vicsek, T.: Fractal growth phenomena, World Scientific Publishing Co., Singapore, 2nd edn., 380 pp., 1992. Nonlin. Processes Geophys., 16, 159–168, 2009 Vignes-Adler, M., Le Page, A., and Adler, P. M.: Fractal analysis of fracturing in two African regions, from satellite imagery to ground scale, Tectonophysics, 196(1–2), 69–86, 1991. Voss, R. F.: Fractals in nature: From characterization to simulation, in: The Science of Fractal Image, edited by: Peitgen, H.-O. and Saupe, D., Springer-Verlag, NY, 21–70, 1988. Wilkinson, D. and Willemsen, J. F.: Invasion percolation: A new form of percolation theory, J. Phys. A-Math. Gen., 16(14), 3365– 3376, 1983. www.nonlin-processes-geophys.net/16/159/2009/