Quantification of pore-size spectrums by solute breakthrough curves

were conducted with sand columns (5 cm id and 5 cm length), repacked with washed sand with a particle size of 2–1, 1–0.45, 0.45–0.325 and <0.325 mm in diameter. The resulting breakthrough curves were divided into approximately 20 segments, and for each segment, a concentration of Cl in an out-flowing effluent was used with corresponding effluent volume and travel time to calculate corresponding pore water velocity 10


Introduction
Quantitative description of effective pore-size spectrum in soils is challenging work due to complex nature of pore-geometry. As reported by Bouma (1991) and Deeks et 20 al. (1999), several methods have been developed and studied to quantitatively analyze pore characteristics in soils and similar porous media. Methods such as computerized classification and recognition, image analysis, and digital binary images provide a 2-D interpretation of soil morphology. A serial sectioning technique can provide a 3-D interpretation of soil pore space morphology, however, this method requires sophisticated 25 equipment that is not readily available and is also time consuming to perform (Deeks 8374 et al., 1999). Water desorption and mercury intrusion methods are also used widely to determine spectrum of soil pore-size. Since water desoprtion methods are generally conducted on small destructed soil samples, this method addresses only the matrix pore-size spectrum (Kung et al., 2005). In addition, as soil is dried and then liquid is forced into it under pressure, a mercury intrusion method may overestimate the ef- 5 fective pore radii (Radulovich et al., 1989). The details on the principles of these and other methods can be found elsewhere (Bouma, 1991;Lal and Shukla, 2004). Besides the aforementioned methods of soil pore spectrum analysis, water (Radulovich et al., 1989) and solute breakthrough curves (Deeks et al., 1999(Deeks et al., , 2008Kung et al., 2000Kung et al., , 2005Kung et al., , 2006 were utilized to analyze effective pore-size spectrum and corresponding 10 water flux distribution in soils. Assessment of pore-size distribution is a major target in the evaluation of water, solute, and air transport in soils (Lal and Shukla, 2004). Kung et al. (2006) stressed that quantitative methodologies for evaluating effective pore-size spectrum are urgently needed.
15 Radulovich et al. (1989) conducted a water breakthrough study, an extension of a concept of the solute breakthrough curves developed by Nielsen and Biggar (1961), to characterize pore-size spectrum and corresponding water flow characteristics of soil macropores in a well-aggregated Oxic Dystropept. They measured water breakthrough outflow on undisturbed soil cores from a drained state to a saturated state, and calcu-20 lated macropore-size and -frequency by travel time and the amount of water collected at consecutive steps. They calculated the micropore-size specrum with a pressure plate apparatus, and then they combined macropore-size and micropore-size spectrums to obtain an entire pore-size spectrum of the soils. Kung et al. (2005) used solute breakthrough curves to analyze pore-size distribution in two intact isolated blocks from a 25 grassland, using the functional relationship between time needed for a tracer to travel a known length of the soil block and radii of the pores through which displacing solution flows with a known velocity. Deeks et al. (1999) related observed changes in solute concentration following a miscible displacement of a 250 mg l −1 Cl solution to 8375 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | characterize size-distribution of large soil conduits. In that study, they calculated volume of discharging traced fluid in a specified time together with the observed rate of movement of tracer within the soils to calculate pore-radius by Poiseuille's equation. Later, Deeks et al. (2008) characterized flow paths and saturated conductivity distribution in 2.4-x 3.4-x 1-m lysimeters. They applied the tracer onto the soil surface and 5 collected soil solution in different locations of the lysimeters, and then they classified the breakthrough curves into statistically distinct pathway types. Kung et al. (2005) proposed field trace mass flux breakthrough patterns to quantify equivalent pore-size spectrum. They used an improved tile drain monitoring protocol, applying a tracer to a narrow strip to a tile-line to measure breakthrough pattern of a tracer mass flux, 10 and used the measured breakthrough curves as surrogates to drive pore spectrum of preferential pathways.
Many miscible displacement studies have shown that one of the most important physical features of the porous media was a magnitude of the volume of water not readily displaced, often termed as immobile or less mobile water. In this study, it was 15 hypothesized that the volume and tracer concentration of discharging traced fluid in a miscible displacement of a conservative tracer can be used to calculate pore-size spectrum. The objectives of this study were (i) to develop a procedure to utilize breakthrough curves of conservative tracers to calculate effective pore-size spectrum for entire range of effective pore-size and (ii) to calculate corresponding pore water veloc-20 ity spectrum in a porous medium. The theory and principles developed were applied to data of BTCs of chloride from five-cm long and five-cm wide sand columns repacked with different-sized homogenous sand. These small columns were preferred to keep the media as uniform as possible and different-sized homogenous sand was repacked to evaluate model's response to gradually decreased particle-size. 25 Hydrodynamic dispersion and convection are the main processes controlling the transport of a conservative tracer (i.e. chloride) in soils and similar porous media (i.e. columns of glass beads). In a miscible displacement, volume of displacing effluent (∆Q), needed to transport a known amount of the conservative tracer by convective 5 flow (assuming no lateral mass exchange and hydrodynamic dispersion for now) from the inlet to outlet of a column may be calculated as follows: where, M is the concentration (in molar) of tracer at the inlet of the column (stock solution), ∆M is the difference in the concentration of the tracer at the outlet between 10 times t i and t i +1 , and ∆V is the volume of effluent collected between t i and t i +1 . Applied to the BTC, Eq.
(2) results in: where, ∆ (C/C 0 ) is the difference in dimensionless concentration (it is equal to ∆M/M in Eq. 1), ∆ (P/P 0 ) is difference in dimensionless pore volume, and P 0 (L 3 ) is the pore vol-15 ume of porous media (amount of effluent held in the column). The property ∆ (P/P 0 )P 0 is equal to ∆V in Eq. (1). It is assumed that the porous system is made up bundles of different sized capillaries. Each size-class (bundle) has a mean equivalent diameter, which controls the rate of effluent flow and transport characteristics. The Eq.
(2) assumes intrinsically that original solute entering the capillaries arrives in the outlet of 20 the capillaries; that is no solutes are lost or gained by lateral diffusive mass exchange and/or no solutes are transported by hydrodynamic dispersion. This means that, the Eq.
(2) simply applies to a piston flow. However, the piston flow never occurs in practice, as lateral mass exchange between effluent in the effective capillaries and in matrix capillaries (pores) takes place, and part of the solute is transported by hydrodynamic 25

8377
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | dispersion. Therefore, the combined effect of lateral mass exchange and hydrodynamic dispersion results in a "dilution of original solution (solution applied from the inlet of the column or simply stock solution)" during the transport in the capillary. In a piston flow, when the dimensionless pore number is one (P/P 0 = 1.0), relative solute concentration at the outlet of the column is unity (C/C 0 = 1.0), meaning that all the effective 5 capillaries in the system are saturated by displacing solution. Thus, to account for lateral mass exchange and hydrodynamic dispersion effect, the right hand side of Eq.
(2) may be multiplied by P/P 0(C/C0=1.0) where all the effective capillaries are saturated by displacing effluent. This results in Eq. (3): where (P/P 0 ) (C/C0=1.0) is (P/P 0 ) at C/C 0 = 1.0 and ∆Q a is the hydrodynamic dispersionadjusted value of ∆Q. Derivation of adjustment factor (P/P 0 ) (C/C0=1.0) ) will be discussed in more detail in Sect. 5. Mean pore-water velocity v i (LT −1 ) for a capillary class i can be calculated by: where, τ is mean tourtuousity of the pores (−), l is the length of the column (cm), and t ci is the cumulative of the time (s) elapsed from the moment at which displacing effluent is first applied to the column. The value for τ was assumed as 1.1 for small columns repacked with well-sorted sands, the lowest value for soils (Radulavich et al., 1989). Radii of capillaries in a size-class can be calculated by (Jury et al., 1991); where, γ is the kinematic viscosity for water (10 −4 cm s −1 ) (Freeze and Cherry, 1979). If volume flow rate per unit area is known, number of pores (n i ) with equivalent radius in a size-class (r i ) can be calculated by Eq. (6) (Jury et al., 1991): 8378 where, In Eq. (7), t i +1 − t i is the time elapsed (s), during which volume of traced effluent is collected to calculate ∆Q i a (see Eq. 3).

5
The average value of pore water velocity (v b ) for all the effective capillaries in the column can be calculated by taking geometric mean of values calculated for all the capillary bundles (size-classes corresponding to segments of BTC) corresponding to entire range of BTC. Therefore, v b of the column can be approximated by: where v 1 is the v i corresponding to the first segment and v n is the v i corresponding to the last segment of the breakthrough curve. Total amount of mobile water in the column can be calculates as: Finally, mobile water fractionation coefficient, β, can be calculated by Eq. (10): where, P 0 is the amount of water held in the column (pore volume; cm 3 ) during the test.

8379
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

Materials and methods
The theory and methodology proposed in Sect. 2 were applied to experimental breakthrough curves (BTCs) of chloride, obtained on small sand columns (5 cm length and 5 cm id) repacked with washed sand.

5
Triplicate sand columns were repacked with sand, screened through 2, 1, 0.425 and 0.325 mm openings. To achieve a uniform packing, the bottoms of the columns were gently tapped on the laboratory bench 20 times during packing.

Miscible displacement tests
Before conducting a miscible displacement test, both ends of the soil column were 10 supported with a fabric. The column then was gradually saturated with a 0.01 M KBr solution from the bottom (van Genuchten and Wierenga, 1977). Upon saturation, the column was positioned on an upright stand. The inlet at the top of the column was connected to a mini-disc infiltrometer with a bottom of 5 cm diameter that just fitted to the top of the column with 5 cm of id. After steady state flow was established under zero 15 tension, approximately 6.0 pore volumes of tracer solution of 0.05 CaCl 2 in 0.01 M KBr solution was introduced into the column with the disc infiltrometer under zero tension. The effluent was collected manually, and analyzed for chloride by the potentiometric titration method (Adriano and Doner, 1982) Following the miscible displacement tests, the column was taken out and placed in 20 an oven with a constant temperature of 105 • C to measure bulk density and calculate total porosity f . The total porosity was calculated by; where D b is bulk density (M L −3 ) of the column and D p is the density of particles (M L −3 ), assumed as 2.65 Mg m −3 . The saturated water content of the column was deemed equal to the volume of pores in sand column. Dimensionless concentrations of chloride were calculated, dividing the concentration of chloride in the effluent measured by potentiometric titration method to the concen-5 tration in the stock solution measured by the same method.

Calculation of v , r, and n for capillary classes (bundles)
Each breakthrough curve was divided into approximately 20 segments, the first segment being with the smallest C/C 0 > 0 and the last one with greatest C/C 0 . For each segment, v (cm s −1 ) was calculated by Eq. (4) using l = 5.0 cm (length of the column), 10 τ = 1.1 (Radulavıich et al., 1989), and t i = time elapsed from the moment at which the displacing effluent (Cl) applied to the column. Radius of each capillary class was calculated using the calculated v i values with Eq. (5), and number pores (n) in each capillary class was calculated using the calculated r i values and Q i values with Eq. (6). The values for Q i was calculated by Eq. (7), dividing the value calculated by Eq.
where t i beginning and t i +1 is the ending time for step (segment) i . Finally, a mobile water partitioning coefficient β was calculated by Eq. (9)

Verification of calculations
Mean pore water velocity of the saturated column (v b ) was calculated using the vvalues calculated for all the segments (pore-size classes) with Eq. (8). To verify the 20 modeling results, the measured values of v at saturated column (v bm ) were plotted (Hscatterplot) against v b values calculated by Eq. (8). Coalescence of the data around a 45 • -diagonal was interpreted as a qualitative measure of the similarity of measured and calculated values of mean pore water velocity (and also mean saturated hydraulic conductivity). To evaluate the similarity quantitatively, the measured values were correlated 25 to modeled values. The extent of correlation coefficient was considered in assessing 8381 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | the success of the technique developed. Since it takes a considerably long time to complete a miscible displacement test in unsaturated conditions, the validation was constrained to saturated conditions.

Results
The properties of the sand columns used in miscible displacement of chloride are given 5 in Table 1. As expected, lowest bulk density D b occurred with the smallest particle size. The triplicates were quite similar in D b , total porosity f , saturated water content θ s , and pore volume at saturation P 0 . Pore water velocity ν (cm s −1 ) decreased gradually versus the decreasing particle-size of the repacked sand. The breakthrough curves of chloride and corresponding v and pore number (n) in 10 each capillary size class are presented in Figs. 1-4. In general, the BTCs gradually shifted to right with decreasing particle-size as indicated by trajectories of BTCs at P/P 0 against C/C 0 = 0.5. Nielsen and Biggar (1961) attributed the right-shift in their study to amount of water not displaced during the miscible displacement. Figures 1-4 show that considerable amount water was not displaced in the columns and that 15 calculated values of mobile water fraction (β) in the columns were not correlated to particle-size. Angulo-Jaramillo et al. (1996) measured mobile water content with a disc infiltrometer, adding a Cl tracer after pre-wetting the soil with water. Their results showed that the mobile water content was a function of mean pore-size, described by an s-shaped curve showing the continuous change in the mobile water content ratio 20 from the microporus to the macroporus system. That β in this study was not related to particle-size may be attributed to the difference in calculation of β between this study and study by Angulo-Jaramillo et al. (1996). The variable β represents fraction of mobile water -mobile water/(mobile water plus immobile or less mobile water) -in the column. The parameter β is somehow similar replicates, the effective pore-size spectrum differs considerably in both number of pores and range of pore-size. For example, most of the transport was controlled by the pores with radii between 0.01 and 0.006 cm in the first replicate (the upper graph in Fig. 1), however, the transport was controlled by pores of completely different sizes in the third replicate (the bottom graph in Fig. 1). In addition, the peak location and height shift of 5 n considerably differs in replicates given in Figs. 1-4. Similar conclusions can be made for other tests conducted with lower particle sizes. For example, Fig. 2 shows that the third replicate has a highly different pore-size spectrum than other two replicates. In addition, although the breakthrough curves of the first and the second replicates are relatively similar; their corresponding pore-size spectrum 10 differs considerably. Second and third replicates in 0.424 mm treatment (Fig. 3) are highly similar in β, however, pore size distribution and v distribution are quite dissimilar. In 0.325 mm treatment (Fig. 4), although first and second replicates resembled in β, their corresponding pore spectrums were exceedingly different.
Mean pore water velocity on saturated columns (v b ) was calculated by Eq. (8) were highly associated as indicated by high correlation coefficient (r = 0.89, P < 0.01) calculated between measured and predicted values.

Discussion
Hydraulic conductivity can be viewed as contributions from all individual pathways to be lumped into a single value (Kung et al., 2005(Kung et al., , 2006. However, this lumped variable may 25 be inadequate when matrix flow is not dominant (Gish et al., 2004). Kung et al. (1995) 8383 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | stressed that it is intrinsically wrong to use a lumped hydraulic conductivity value to predict convective contaminant transport through preferential flow pathways. Results of this study highly agreed to their conclusions. The lumped hydraulic conductivity value measured for the first replicate of 2-1 mm size class treatment was 0.044 cm s −1 . However, the Fig. 1 shows that conductivity of water saturated domains ranges from 5 0.026 to 0.124 cm s −1 and that tremendous of flux is controlled by pores corresponding to water flux rate greater than 0.044 cm s −1 , suggesting that use of a single lumped parameter may lead erroneous conclusions as suggested by Kung et al. (2000Kung et al. ( , 2005Kung et al. ( , and 2006). Similar conclusions can be made from the other graphs in Figs. 1-4. In our calculations, the measured K s for the 2-1 mm size treatment were 0.044, 10 0.046, and 0.043 for replications 1, 2, and 3, respectively. By these values, it may be concluded that the three replicates are highly similar in K s . However, Fig. 1 shows tremendous differences in range, number, and distribution of pores in these replicates. In the first replicate, the transport is mainly controlled by pores with radii between 0.008 and 0.006 cm ( Fig. 1 top graph), while in second and especially in third repli-15 cate, completely different pore-sizes controlled the water and chemical fluxes. Similarly, for the r < 0.325 mm particle-size treatment (Fig. 4), the water and solute flux was controlled completely different pore-sizes despite their highly similar K s values. These pore-spectrums, measured in highly homogenized conditions (small columns repacked with sieved washed sands) were considerably different. In structured media, far more 20 differences are expected, and these show that a single lumped parameter is inadequate to represent all these differences in pore characteristics that ultimately have an important impact on water and chemical fluxes in soils. Shape of BTCs is generally interpreted to qualitatively assess physical and chemical characteristics of transmitting media. The method developed in this study made 25 this interpretation more fruitful, translating the details in a BTC to better interpretable pore-size spectrum and corresponding pore-water velocity distribution information. For example, in Fig. 3, it is highly difficult to compare the BTCs of the second and third replicates in terms of pore-size spectrum and pore water velocity distributions without graphs given in middle and right columns since two graphs are similar. However, the pore-size distribution and corresponding pore-water velocity distribution graphs reveal many details in pore characteristics.
The capacity of the soils to store and transmit water and nutrients is hinged on the porous nature of the soils (Kung et al., 2005). Attempts have been made to better rep-5 resent the pore characterisitics in flow domains in water and solute transport studies. The two-domain approach (van Genuchten and Wierenga, 1976) and multiple domain approach (Gwo et al., 1995) were proposed to compensate the drawback form use of a single lumped variable, mean hydraulic conductivity (Kung et al., 2005). Durner and Fluhler (1996) suggested that a continuous pore spectrum is necessary to simulate 10 chemical transport and proposed that multiple domain approach should be used as an alternative to the lumped parameter, hydraulic conductivity. The stream tube model from Toride et al. (1995) was fundementally similar to multiple domain approach. The stream tube model assumes that the each stream tube can be treated theoretically as an independent domain. The multiple domain approach is highly daunting in mea-15 suring hydraulic conductivity for a soil that is divided further into more domains. Also, defining the boundaries among domains and tubes is very difficult and the success of using multiple domains depends firmly on the accurate measurement of hydraulic conductivity of these domains (Kung et al., 2005). The procedure developed in this study allows us to divide the system with as many domains (segments) as desired and to 20 define boundaries among domains.
Methodlogy, partly similar to one developed in this study, was used by others (Kung et al., 2005) to quantify the spectrum of macropore-type preferential pathways using BTCs of conservative tracers. They used an improved tile drain monitoring protocol and they applied tracers to a narrow strip near the tile line and measured breakthrough 25 patterns of tracer mass flux. They used BTCs as surrogates to find the parameters of the proposed function to derive the pore-spectrum of preferential pathways. However, their model did not account for hydrodynamic dispersion and chemical adsorption effects. In a similar fashion, Deeks et al. (2008) characterized flow paths and saturated 8385 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | conductivity in a large soil block in relation to chloride BTCs. In another study, Deeks et al. (1999) found that solute movement was predominantly confined to pores with radii between 0.3 and 1.0 mm. The method developed by Radulovich et al. (1989) predicts pore-size spectrum as combination of water breakthrough curves and moisture release curves. However, use of the moisture release curves is criticized since these curves 5 are typically measured using homogenized soils where the soil structural pores were destroyed.
The model developed in this study can be used to measure the size spectrum of all hydraulically active pores in a porous medium. The model accounts for hydrodynamic dispersion as well as lateral mass exchange between mobile domains and between 10 mobile and stagnant (or less mobile) domains. The following mental experiment leads us to intrinsically introduce an "adjustment factor (P/P 0(C/C0=1.0) )" used in (Eq. 3) to account for effect of hydrodynamic dispersion and lateral mass exchange: In "piston flow" (Fig. 6), whole pores in the system are in same radius and saturated by the traced effluent where P/P 0 = 1.0. The entire effluent in the system contributes transport and, 15 therefore, there is no immobile water in the column. Furthermore, all the chemical is transported by convective flow and no hydrodynamic dispersion presents in the system as indicated by P/P 0 = 1.0 against C/C 0 = 1.0. Therefore, the conditions for the piston flow may be used as a reference point (P/P 0 = 1.0 against C/C 0 = 1.0) to account for "dilution effect" caused by hydrodynamic dispersion and lateral exchange between mo-20 bile/less mobile and immobile regions. Any discrepancy (early appearance and delay) from this reference point can be accounted as effect of hydrodynamic dispersion and lateral solute exchange among the flow domains. Therefore, as an account of hydrodynamic dispersion, the BTC diverse from the piston flow, taking different values of P/P 0 against C/C 0 as a response to dilution effect caused by hydrodynamic dispersion 25 and interactions among flow domains. Apart from this, we may multiply the Eq. (2) by P/P 0(C/C0=1.0) as an adjustment factor, resulting in Eq.
(3). The model developed in this study differs from the those by others in the following aspects. Firstly, the developed model calculates the amount of mobile water for any desired segment in a BTC. This allows to calculate corresponding number for pores by well known Pouiseille equation, which in turn, leads calculate size spectrum of all effective pores in the system. And secondly, as mentioned above, the model accounts for effects of hydrodynamic dispersion and lateral interactions among domains in a porous system. 5 The spectrum of hydraulically active pores should be determined to better model water and chemical flux in soils. The pore-size is the key property controlling the convective transport through soils since the flux rate is proportional to the forth power of equivalent pore radius according to well known Poiseuille equation (Jury et al., 1991). In this study, it is assumed that the adjacent flow domains (domains made up similar 10 sized capillary pores) were interacting.
Water flow in soils depends mainly on pore geometry and pore-size distribution, which are controlled by soil structure and soil texture. It is considerably difficult to quantify the relation between soil structure and water and solute transport in soils (Deeks et al., 1999). The theory and applications adapted here suggests there are several 15 directions for future research efforts. Development and application of laboratory techniques for measurement of effective pore-size distribution together with convective flow of water in soils with diverse structure may result in a significant breakthrough in understanding of water and chemical transport in structured soils.
The model developed in this study has several drawbacks. First, it takes quite long 20 time to complete a miscible disparagement test in structured clay soils. In this case, the P/P 0 of the BTC may be extrapolated to C/C 0 = 1.0, and this resulting extrapolated value of P/P 0 can be used as an adjustment factor in Eq. (3). Second, in highly structured or macroporous soils, it may be very difficult to catch large macropores, which is very important in preferential transport of water and solutes. This problem can be tack- 25 led, conducting a more frequent sampling at the early stage of test. Finally, the BTC should be an increasing function to make calculations, otherwise one may come up with unreasonable values of v, n, and r. Therefore, the measurement device should be very precise that it can measure slight changes in solute concentrations in consecutive 8387 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | samples. The same is also true for measuring first concentration breakthrough in very early stage of test where largest pores start to discharge the displacing effluent.

Conclusions
A theory was proposed to calculate effective pore-size spectrum by breakthrough curves (BTCs) of conservative tracers in soil and similar porous media. The theory 5 was tested using breakthrough curves of chloride from small columns repacked with clean uniform sand material. Miscible displacement tests of CL were conducted on sand columns repacked with sieved sands of 1-2, 0.425-1, 0.3. 25-0.25 and, <0.325 mm, and the resulting BTCs were used to calculate distributions of pore-water velocity and corresponding pore-size 10 spectrums. The rationale behind use of these small columns was to evaluate model in simplest and uniform conditions to decrease number of uncontrolled variables.
Decreased particle-size resulted in increased pore number, slightly increased water content, and decreased mean pore water velocity. A considerable difference occurred among treatments (particle-sizes) as well as replicates of treatments shown by range 15 and height shift of pore-size spectrums. Replicates similar in water content, mean pore water velocity, bulk density, and mobile water content yielded different pore-size spectrums as shown by range and height shift of graphs. Pore water velocity distribution in the columns showed that use of a lumped parameter of hydraulic conductivity may not be adequate in chemical transport modeling studies. The developed model ac-20 counts for hydrodynamic dispersion and lateral mass exchange between mobile and immobile/less mobile domains. Another study is currently underway to evaluate the proposed model on disturbed and undisturbed soil columns (30 cm long and 8.5 cm id) to extend its use to more complicated conditions.