Field-scale solute transport parameters derived from tracer tests in large undisturbed soil columns

Additional keywords: convection–dispersion model, mobile–immobile model, stochastic stream tube model, time-domain reflectometry, parameter optimisation, contaminant transport.

Protection of shallow groundwater systems will greatly benefit from our ability to accurately predict the field-scale transport of pollutants from the soil surface to the groundwater table. Remediation of shallow and deep vadose zones contaminated by organic and inorganic substances often depends on the success of delivery of surface-applied or injected amendments into the vadose zone. Field-scale solute transport modelling is indispensable for evaluating the feasibility of such clean-up methods (Brusseau et al. 2010; Zhang et al. 2012). However, field-scale solute transport is significantly influenced by the heterogeneous nature of flow and transport properties across the field (Jacques et al. 1998; Seuntjens et al. 2002; Mallants et al. 2011; Stauffer and Lu 2012). In the last three decades, research both theoretical and experimental has been carried out to investigate solute transport behaviour in field soils (Van Ommen et al. 1989; Ellsworth et al. 1996; Forrer et al. 1999; Ghodrati 1999; Dusek et al. 2008; Robinson et al. 2012). Theoretical investigations of solute transport in heterogeneous, stratified porous formations often rely on stochastic principles. For instance, the stream-tube model (STM) regards the field as a series of independent vertical stream tubes. The model parameters for each tube are viewed as realisations of a stochastic process. Bresler and Dagan (1981) investigated the combined effects of pore-scale dispersion and velocity variability on solute distributions in the field. Effects of local-scale transport parameters on field-scale transport for a variety of boundary and initial conditions were investigated by Toride and Leij (1996) assuming the convection–dispersion equation (CDE) in each stream tube. Ward et al. (1995) analysed the effects of variable layer thickness in a stochastic STM on solute spreading.

Testing of the predictive capacities of the above theoretical models is important and typically requires well-controlled flow conditions in the field, or alternatively, in the laboratory invoking flow conditions that mimic those in the field. Several studies have concentrated on the collection of field tracer data for this purpose (Jury 1982; Schulin et al. 1987a; Jaynes et al. 1988; Butters et al. 1989; Köhne and Gerke 2005). Solute transport at the field-scale, however, is known to be complicated by several processes, such as the presence of variable horizon thicknesses or tonguing (Van Wesenbeeck and Kackanoski 1991), and preferential flow (Roth et al. 1991; Zehe and Flühler 2001; Jarvis 2007). Furthermore, calibration and validation of theoretical models requires data collected from displacement experiments with well-controlled boundary conditions, which are often difficult to establish in field soils.

As a surrogate for actual field experiments, numerical experiments with hypothetical field soils have gained insight into the sensitivities of field-scale behaviour with respect to particular local-scale heterogeneities. For instance, Russo et al. (1989) investigated the influence of immobile water on solute transport subject to boundary conditions corresponding to cycles of infiltration, evaporation, and redistribution. Tseng and Jury (1994) compared the stochastic-convective log-normal transfer function model (CLT) and the CDE model to predict area-averaged transport in a hypothetical, two-dimensional heterogeneous field soil. Vanderborght et al. (1998) compared different approaches to calculate statistics of the solute travel time and displacement based on spatial covariance functions of the pore-water velocity. Based on simulated water flow in a two-dimensional, macroscopically homogeneous Miller-similar medium, Roth and Hammel (1996) demonstrated that the local structure of the velocity field (expressed by the form of the spatial autocorrelation function for the scaling factor) and the subscale hydrodynamic dispersion are of minor importance for the field-averaged transport process. Vogel and Roth (2003) proposed the ‘scaleway’ as a discrete hierarchical approach for modelling flow and transport, based on the explicit consideration of the full three-dimensional representation of the macro structure, while microscopic heterogeneities were averaged and replaced by effective descriptions.

In addition to field and numerical experiments, transport experiments for controlled laboratory conditions with undisturbed soil columns provide useful information on the local-scale solute behaviour if the natural layering is preserved (Schulin et al. 1987b; Seyfried and Rao 1987; Beven et al. 1993; Seuntjens et al. 2001). If the inherent spatial heterogeneity of the soil is represented by taking samples at different spatial locations across the field, and if the particular flow conditions at the local scale are also present at the field scale, e.g. steady-state vertical flow, then one should be able to evaluate the field-scale solute transport by appropriate averaging of the local-scale processes. Extrapolation of results obtained from laboratory experiments to field conditions may then yield satisfactory predictions.

This paper presents a solute-transport experiment carried out under controlled laboratory conditions using 14 undisturbed soil columns, 0.3 m in diameter and 1 m long, collected from a field transect. Local-scale solute breakthrough curves (BTCs) were measured at six different depths in the soil columns during steady saturated flow using time-domain reflectometry (TDR). The BTCs were analysed in terms of the one-dimensional CDE and the mobile–immobile non-equilibrium transport model (MIM). By assuming the field as a series of independent, vertical soil columns, field-scale mean transport was also analysed in terms of the STM with a stochastic pore-water velocity, v, and dispersion coefficient, D. The effect of the local-scale dispersion on field-scale transport was also discussed by considering the Peclet number, <v> x/D, in terms of the observation scale, x.

Thirty undisturbed soil columns, 0.3 m in diameter and 1 m in length (Fig. 1), were collected along a 31-m-long transect from the Bekkevoort experimental field, east of Leuven in central Belgium. The field plot was in an orchard, with the transect between two rows of trees. The surface between the trees was covered with permanent grass. The soil profile contained three different horizons: Ap, C1 and C2. The Ap and C1 horizon were pedogenetically identical, i.e. they had both developed on colluvial material. The C2 horizon was an old textural B horizon, which had been covered by erosion material. The average thickness of the Ap horizon was 0.25 m, whereas the mean thickness of the C1 and C2 horizons was 0.3 and 0.45 m, respectively, with a sharp boundary between all soil layers. Each layer had a fairly constant thickness along the transect. Physical and chemical properties of each soil horizon are listed in Table 1. The clay content increased from 12.7% in the Ap horizon to 16.6% in the Cl and 21.8% in the C2 horizon. The higher clay content in the C2 horizon is due to accumulation of migrating clay particles. The soil was classified as an Udifluvent (Eutric Regosol) (IUSS 2006) and contained a large number of macropores throughout the profile (Mallants et al. 1997). The macropores consisted mainly of decayed root channels and earthworm holes.

Fig. 1.  Diagram of the experimental setup for transport experiments on undisturbed soil columns. The inset (not to scale) shows the two-rod parallel transmission line: L is probe length, b is wire diameter, and 2d is the wire spacing.

Table 1.  Physical and chemical properties of the soil profile
CEC, Cation exchange capacity; OC, organic carbon

From the same field plot, additional 0.2-m-long, 0.2-m-diameter soil columns had previously been investigated to test the applicability of TDR for measuring transport processes at shallow depth (Mallants et al. 1994). The sampling locations were exactly between two wheel tracks so that the soil had not been compacted nor been in contact with any machinery. First, the soil was cleared of vegetation (grass) and the humus layer was removed. Sampling of the soil columns was done by means of polyvinyl chloride (PVC) cylinders. The PVC cylinders had an internal diameter of 0.3 m, and the inner wall was greased to have close contact between soil and the wall. The bottom end was provided with a sharp-edged steel ring to facilitate intrusion into the soil, and the cylinders were gently driven into the soil using a hydraulic jack. Once the cylinder was completely filled with soil, the soil around it was excavated gradually such that the column could be isolated and the bottom end cut off by means of a steel plate that was sharpened at one side. Next, a PVC end cap was placed on both ends of the cylinder. This unit was lifted hydraulically onto a truck and transported to the laboratory. Thirty undisturbed soil columns were sampled in this way with a separation distance of 1 m. At one location a sampling distance of 2 m was used because of difficulties in driving the cylinder into a locally compacted layer.

In the laboratory, a perforated 0.005-m-thick PVC plate was attached to the bottom end of the columns to prevent the soil from falling out. A nylon-type mesh was glued on top of the perforated plate to allow the application of a suction of up to 70 mbar (7 kPa). This assembly was attached to the bottom of the column together with an enclosed drainage system (Fig. 1). Drainage water from this column could be collected in small sampling bottles via a polyethylene (PET) tube. The set-up also allows a watertable to be maintained at the bottom of the column or a negative pressure to be imposed by means of a hanging water column. Solute-free water could be applied to the top of the columns by means of a constant-rate drip-irrigation system. A 0.01-m-thick fibre wool cover was put on top of the soil surface to reduce crust formation during irrigation. The top boundary condition was established by allowing water to pond for a maximum depth of 0.01 m. A spillway system diverted the excess water and secured a maximum water depth of 0.01 m. A second independent drip-irrigation line was constructed for the application of the tracer solution. In this way, it was easy to switch from solute-free water to tracer solution and vice versa. All columns were saturated from the bottom by stepwise increasing the water level in the small reservoir.

To ensure steady-state saturated flow, soil water content, θ, and pressure head, ψ, values were monitored by TDR and tensiometers installed at six different depths (0.05, 0.15, 0.3, 0.45, 0.6, and 0.8 m), i.e. two observation depths for each horizon. The TDR probes consisted of two parallel waveguides with a distance of 0.025 m between the centre of the rods. The diameter of the stainless steel rods was 0.0005 m. The 0.25-m-long stainless-steel probes were directly connected to the inner and outer conductor of a coaxial cable (see inset Fig. 1) in a way similar to the probe design of Ledieu et al. (1986). The volume of soil that is effectively sampled for water content using this probe design corresponds to an approximate cylinder with a diameter of 0.066 m. Details for the calculation of the sampling volume for TDR probes can be found in Knight (1992). Volumetric water contents were inferred from the travel time of an electromagnetic wave through a volume of wet soil. Travel times were obtained from the screen of a Tektronix 1502B cable tester (Tektronix Inc., Beaverton, OR, USA) and were converted to soil water content using the calibration relationships established by measuring travel times for repacked soil at known but different water contents. The pressure head was measured by pushing a hypodermic needle into a rubber membrane mounted on an air-filled reservoir. The needle was connected to a pressure transducer, which converted the pressure in the needle to a digital pressure value on a LCD screen (Thies-Clima, Göttingen, Germany).

The steady-state water contents, θ, as observed with TDR show little horizontal and vertical variation. Generally, water contents at 0.45 m depth are slightly smaller than those at 0.15 and 0.80 m. The maximum difference in water content between soil horizons is only ~5% by volume. The bulk density, ρ_b, also changes little horizontally, except for some fluctuations occurring near the surface. Hence, the effects of variable θ and ρ_b on the transport process are not considered further.

When the soil was saturated, solute-free water was applied through the drip irrigation system at a rate sufficient to maintain ponding. After establishing steady-state flow conditions, the application of solute-free water was interrupted and the remaining water allowed to infiltrate. Next, a 7 × 10^–3 m CaC1₂ solution was applied continuously for a period of at least 79 h. This long application time guaranteed that, at least at the first observation depths (x = 0.05 m), the initial soil solution was completely displaced by the applied tracer solution (Mallants et al. 1996). For some extreme cases, solute had to be applied for as long as 664 h to saturate the soil at the first observation depth (x = 0.05 m). Such long application times resulted in an uncertain calibration, which would affect the derived transport parameter v and D (Mallants et al. 1996). Hence, only the columns with an application time, t₀, of 79 h will be discussed for field-scale mean transport in this paper (14 in total). This still presents a sufficiently large dataset for testing field-scale transport processes. Finally, solute-free water was subsequently added again to the soil until all of the solute had leached out of the columns. Tracer concentrations at the six observation depths were obtained by converting estimates of the bulk soil electrical conductivity, EC_a, to resident concentrations. The calibration methods for undisturbed soil were discussed in detail by Mallants et al. (1996). Values of EC_a were estimated with TDR using procedures as discussed by several authors (Kachanoski et al. 1992; Wraith et al. 1993; Vanclooster et al. 1993; Mallants et al. 1994; Ward et al. 1994).

Solute transport in each soil column was analysed in terms of the CDE and the MIM models. The CDE for non-reactive tracers is given by:

where c is the resident concentration (g/cm³), D is the local-scale dispersion coefficient (cm²/h), v is the pore-water velocity (cm/h), defined as q/θ with q the Darcian flux density, x is depth (cm), and t is time (h).

The MIM model, which is often invoked to describe asymmetrical concentration distributions displaying early breakthrough and tailing, is given by:

where θ is the volumetric water content, α is a first-order mass transfer coefficient (day^–1), and the subscripts m and im refer to the mobile and immobile liquid regions, respectively. The transport parameters in the CDE and MIM models were fitted to observed BTCs using an extended version of the CXTFIT non-linear parameter estimation program of Toride et al. (1995). Since the TDR only detected a relatively small region of the whole flow domain, this region may or may not include preferential flow channels. Following the approach of Mallants et al. (1994), independent estimates were first obtained of the pore-water velocities, v (for the MIM this is equivalent to v_m θ_m/θ), from the first normalised temporal moment of the BTC (i.e. calculating time moments is independent from assuming any transport process model), and D was determined for the CDE, and D_m, θ_m/θ and α for the MIM. In so doing, the number of fitting parameters is reduced by one, which increases the uniqueness of the remaining fitting parameters. It was further assumed that TDR measures total resident concentration, i.e. c_T = c for the CDE (Eqn 1) and c_T = (θ_mc_m + θ_imc_im)/θ for the MIM (Eqn 2).

When the field is viewed as a set of homogeneous vertical stream tubes without exchange between the stream tubes, field-scale transport may be described by averaging local-scale concentrations over all stream tubes (Bresler and Dagan 1981; Simmons 1982). In order to investigate the effect of local-scale velocity and dispersion fluctuations on field-scale transport, the stream-tube concept was applied to the observed column BTCs. It was also assumed that the local-scale transport is described by the CDE (Eqn 1). The field-scale mean concentration, <c>, in terms of the STM with stochastic v and D is given by (Toride and Leij 1996):

where f(v, D) is the joint probability density function (pdf) for v and D, which may be described with the bivariate lognormal distribution:

where:

in which µ is used for the mean and σ for the standard deviation of ln(v) or ln(D), and ρ_v, D is the correlation coefficient between Y_v and Y_D. It is assumed that all random parameters are statistically stationary and obey the ergodicity principle. An average concentration across the field, obtained from a sufficient number of samples, is then identical to the area-averaged concentration:

where A denotes the area of the field and n is the number of samples. The approach, hence, assumes that horizontal field-averaged transport behaviour can be estimated from the local-scale column BTCs. When v and D are perfectly correlated, i.e. ρ_v, D = ± 1, the local-scale D can be expressed in terms of v (Toride and Leij 1996):

where <v> and <D> are ensemble averages of v and D, given as <v> = exp(µ_v + σ²_v/2) and <D> = exp(µ_D + σ²_D/2). A special case for Eqn 9 occurs when ρ_v, D = 1 with σ_v = σ_D, in which case Eqn 9 reduces to:

where <D>/<v> = λ, with λ the field-scale dispersivity (m). The relationship between D and v following Eqn 10 has been extensively used to express the dependence of D on v (e.g. for a review of D–v relationships see Beven et al. 1993; Vanderborght and Vereecken 2007). The field-scale mean concentration <c> can be obtained by substitution of the appropriate analytical solutions for the CDE and the pdf into Eqn 4 (Toride et al. 1995).

A typical example of measured local-scale BTCs in one single column together with the fitted local-scale CDE and MIM transport models is given in Fig. 2a–f for six observation depths. Note that the fitted transport parameters reflect the average properties from the surface to the observed depth. Table 2 presents fitted parameter values for the two models, together with the measured pore-water velocity at depth x, v_o = q/θ, where θ is the steady-state water content at depth x. Since the pore-water velocity, v, was obtained from the first time moment of the observed BTC and subsequently fixed during parameter fitting for the CDE and MIM, the mean travel time τ to a given depth x is identical for the measured (by using time moments) and fitted curves. Comparison between v_o and v for all depths indicates consistently larger values for the former compared with the latter, except at x = 0.8 m where the two values are relatively close. Values for v differed significantly between the first and second observation depth (0.87 and 2.73, respectively) whereas flow became more uniform deeper in the profile. As discussed by Mallants et al. (1994), the soil was more heterogeneous near the surface. The TDR probe at x = 0.05 m presumably could not detect all of the preferential flow channels that were present in the soil. As a result, mainly the slower transport regions may have been measured, resulting in a relatively symmetrical BTC with a fairly small value for v.

Fig. 2.  Measured (Δ) and predicted BTCs using the CDE and MIM transport model at six depths in an undisturbed soil column.

Table 2.  Estimated transport parameters for the CDE (D, cm²/h) and MIM (D_m, cm²/h), θ_m/θ, α (h^–1), and α^–1 (day) transport model fitted to the data in Fig. 5
Pore-water velocity v was obtained independently from moment analysis and fixed during parameter fitting; v_o is measured pore-water velocity (cm/h). Values in parentheses are standard errors (SE). v_m = v/(θ/θ_m). where is the mean of y_obs

The dispersion coefficient D pertaining to the CDE varied considerably between the first three depths. From 0.45 m onwards, D became almost constant. The local-scale dispersivity, λ = D/v, was found to vary between 0.02 m at x = 0.05 m and 0.178 m for x = 0.8 m. Similar large values for λ have often been observed for undisturbed soils (Dyson and White 1987; Beven et al. 1993; Gelhar 1993; Vanderborght et al. 2001). Vanclooster et al. (1995), for example, reported a dispersivity λ = 0.31 m for their sandy soil in a 1-m-long lysimeter. Van Wesenbeeck and Kackanoski (1991) reported a maximum dispersivity of 4.81 m for a sampling depth of 1.83 m in a field test under ponding conditions.

The mobile-phase velocities, v_m = v (θ_m/θ), were less variable compared with v, although v_m was still somewhat smaller than v_o. Fractions of mobile water, θ_m/θ, showed an increasing trend with depth except at x = 0.15 m, where nearly all of the water was mobile, i.e. θ_m/θ = 0.99. The fitted first-order mass transfer coefficient, α, indicates that deeper in the soil profile, less time is required to equilibrate the immobile phase concentration with respect to the mobile phase. That is, as the solute moved deeper in the soil profile, more solute mixing occurred between preferential pathways and bypassed flow regions.

Based on the values for the coefficient of determination, r², the MIM model with three fitting parameters (D_m, θ_m/θ, and α) generally performed better than the CDE with a single fitting parameter (D). For observation depths greater than 15 cm, the CDE always underestimated the BTCs at small times, and overestimated the curves at large times. In order to compare the transport parameters between the CDE and the MIM, an effective dispersion coefficient, D_eff, was evaluated, obtained from the MIM model parameters (Valocchi 1985; van Genuchten and Dalton 1986):

in which the first term is the hydrodynamic dispersion in the mobile phase, and the second term reflects the source–sink effects of solute transfer into and out of immobile zones. Note that the hydrodynamic dispersion is generally proportional to v, whereas in the second term, D_eff is proportional to the square of v. Values of D_eff, calculated with Eqn 11, for x = 0.6 and x = 0.8 m were, respectively, 74 and 48 cm²/h, whereas the CDE estimates of D in Table 2 were 38.7 and 34.1 cm²/h. The larger values for D_eff are due to the solute exchange between mobile and immobile zones resulting in asymmetrical BTCs, which require a larger (effective) dispersion coefficient to describe it.

Parker and Valocchi (1986) indicated that transport in structured soils can be approached with the equilibrium CDE when the assumption of the local equilibrium (LE) is met. As the travel distance, x, increases, the mean residence time, xlv, becomes larger relative to the characteristic time for the solute exchange between mobile and immobile phases, 1/α. Validity of the LE assumption, therefore, improves as solute moves down. This LE effect could result in a good agreement between D_eff and D for the CDE at 0.8 m. These results illustrate the potentials of describing solute transport at greater depths with the CDE using an effective D given by Eqn 11.

In addition to the LE effect, this close agreement between CDE and MIM at x = 0.8 m could be caused in part by the presence of a compacted soil layer, which was observed in some columns. When horizontal solute mixing increased in the compacted layer, solute transport could be described by means of the convective–dispersive process (Jury et al. 1991). In a few columns this compaction resulted in much lower peak concentrations at x = 0.8 m (results not further presented here).

To evaluate the relationship between the dispersion coefficient, D, and the pore-water velocity, v, for the CDE, linear regression was used to estimate the coefficient, λ, and the exponent, n, in the following equation:

where n is a parameter generally ranging from 1 to 2 (Saffman 1959; Bear 1972; Beven et al. 1993), although values <1 and >2 have also been reported (Bowman and Rice 1986; Jaynes et al. 1988). As discussed for Eqn 11, when hydrodynamic dispersion (the first term in Eqn 11) is the main spreading mechanism, n will be close to 1, while values of n close to 2 represent processes governed by intra-aggregate diffusion (the second term in Eqn 11). Results of the regression analysis between ln(D) and 1n(v) are given in Fig. 3 using all data from three soil horizons. Combining all data at the 0.8 m depth leads to a ‘field-scale’ dispersivity <D>/<v> = λ of 0.656 m, and n = 1.92 with the correlation coefficient, ρ_v, D = 0.74

Fig. 3.  Correlation between transport parameters ln(D) and ln(v) from the CDE (complete data).

This relatively high value for n indicates that non-equilibrium effects on solute spreading were dominant compared with the hydrodynamic dispersion in the mobile phase. Under ponded conditions, Biggar and Nielsen (1976) found a similar correlation between v and D for their soil, i.e. ρ_v, D = 0.795 (by combining data from all plots and all depths). Padilla et al. (1999) determined an n value of 1.99 and a value of 0.085 cm for λ based on the CDE model for a transport distance of 0.25 m. Of further note in Fig. 3 is the greater heterogeneity in the ln(v) data for the Ap horizon compared with ln(v) values for the deeper soil layers (Cl and C2). This likely reflects again the more homogeneous nature of the flow regime deeper in the soil profile.

Values of v and D from all depths were combined and used to construct fractile diagrams. The raw and log-transformed v and D values are shown in Fig. 4. The regression models are only used to determine whether the parameters are normally or log-normally distributed. The data seem to be adequately described by both the normal and log-normal distributions. Based on the complete dataset, v and D had <v> = 0.82 cm/h with σ_v = 0.62 cm/h, and <D> = 52.8 cm²/h with σ_D = 1.6 cm²/h, where σ is the standard deviation of the log-transformed variable as used in Eqn 5. The coefficient of variation (CV) = (exp(σ²) – 1)^0.5 for a log-normally distributed v and D is, respectively, 68% and 350%. These results for λ lie in the range 0.15–1.56, which was compiled by Jury (1985) based on data in the literature. The 95% confidence intervals for v and D in terms of the log-normal distribution were 0.69–0.94 cm/h and 12.2–93.5 cm²/h, respectively. The larger confidence interval for the latter parameter indicates that D could be less accurately determined from the BTC data.

Fig. 4.  Fractile diagram (a) for v (●) and ln(v) (○), and (b) for D (●) and ln(D) (○); (x – µ)/σ represents probability units, with µ and σ the mean and standard deviation of x, where x is v, ln(v), D, or ln(D).

Values for the ensemble mean CDE parameters, together with the standard deviation for the log-transformed values for each depth, are given in Table 3. The relatively higher values for σ_v in the top layer demonstrate the reduced variability in the solute transport process at greater depths, with mean parameter values almost identical at 0.6 m and 0.8 m. Mean values for the MIM parameters are also listed in Table 4; parameter values at 0.6 m and 0.8 m are also identical as was found for the CDE. The fraction of mobile water increased with depth to a value of ~0.40. Judging from the results in Tables 3 and 4, a somewhat more homogeneous transport process occurred in the deeper layers compared with the near-surface layer.

Table 3.  Mean parameters for the CDE transport model

Table 4.  Mean parameters for the MIM transport model

The mean solute concentrations were evaluated for the 14 soil columns whose solute application time was 79 h. As shown in Eqn 8, it is assumed that the mean BTCs can be regarded as an area-averaged BTC across the horizontal plane. The procedure presented here considered flow and transport phenomena related to scale transition in undisturbed soil under identical boundary conditions (saturated and vertical flow). Measured field-scale BTCs and their 95% confidence limits (based on Student’s t distribution), together with the stream-tube predictions, given by Eqn 4, for all depths are presented in Fig. 5. Parameter values presented in Table 3 were used for the STM at each depth. The correlation coefficient, ρ_v, D = 0.74, as shown in Fig. 3, was used for all depths. Widest confidence limits are found (largest variances, see also Fig. 5) before and after the peak in the BTC. Although the STM parameters were obtained from the fitted local-scale parameters, the predicted and observed <c> distribution agree relatively well. However, the STM failed to predict the early arrival of solute and resulted in lower concentration peaks and longer tails. This result is caused primarily by the fact that the CDE does not describe accurately enough local-scale transport as was demonstrated in Fig. 2. Because the STM assumes a constant velocity with depth, it becomes clear from the above discussions that heterogeneities in the vertical direction cause predictions at greater depths based on parameters estimated at shallow depths to underestimate or overestimate the true field-scale solute distribution. Therefore, no attempt was made to evaluate the STM in this respect.

Fig. 5.  Measured and predicted field-scale mean resident concentration distributions at six depths.

In addition to the mean concentration, the variance can be used to characterise fluctuations in the local-scale concentration at a particular depth. The variance in terms of the STM with stochastic v and D may be given by (Bresler and Dagan 1981):

Field-scale variances were evaluated in a way similar to field-scale mean concentrations (Toride et al. 1995). The variance of the concentration characterises the degree of variation over the horizontal plane at a given depth x and time t. Figure 6 presents observed variances for 14 columns against time as well as predicted variances described with Eqn 13 using the same parameter values as in Fig. 5. Predicted and observed variances agree well at the first observation depth. Deeper in the profile, the predicted variance underestimates the observed variance and this difference increases with depth. The observed disparities at greater depths are probably the result of the relatively small number of observations used to calculate the transport parameter used in Eqn 5. Another reason is the observation that at greater depths the MIM is systematically better than the CDE model. In the calculation of the variances, a model was used that assumes the CDE is the governing local transport model. Because the CDE is not the best local-scale model, the variance estimations based on the CDE may explain in part the observed discrepancies.

Fig. 6.  Observed and predicted variances for resident concentrations at six depths.

Notice that the variance profile has a double peak with a relative minimum at the time peak concentration was observed. This bimodal behaviour for a pulse input was also demonstrated in a numerical study by Burr et al. (1994). At the time of the peak, the CV ranges from 8.9% at the first observation depth to 27% at x = 0.8 m. That is, relative variability in the local-scale peak concentrations are relatively small, and largest at the greatest observation depth. These small values of CV are the result of the relatively large solute application time. Continuous solute injection counteracted the horizontal concentration fluctuation caused by random solute transport and resulted in smaller variances (Toride and Leij 1996). The magnitude of the concentration variance itself only minimally decreases with increasing travel depth due to an increased mixing as the solute traverses through the medium.

Bresler and Dagan (1981) demonstrated that the effect of local-scale dispersion on field-scale transport is negligible compared with the effect of heterogeneous flow; most previous studies of stochastic transport have, therefore, neglected local-scale dispersion (Jury 1982; Destouni and Cvetkovic 1991). Toride and Leij (1996) discussed the STM using travel time moments. They concluded that the Peclet number, <v> x/<D> (= x/λ), in terms of the observation scale, x, is an index to determine whether or not local-scale dispersion can be neglected, while the effect of the standard deviation, σ_D, and the correlation coefficient, ρ_v, D, on solute spreading is relatively minor compared with the effect of <D>.

Bresler and Dagan (1981) assumed a local-scale dispersivity, λ, of 0.03 m. Their corresponding Peclet number was 26.7 at x = 0.8 m, whereas, as shown in Table 3, values for x/λ in our experiment are considerably smaller because of the large dispersion coefficient due to the mobile–immobile non-equilibrium effect. The following discussion demonstrates the effect of local-scale dispersion on field-scale mean transport using parameter values obtained at 0.8 m depth. At this depth, all heterogeneities are included in the BTCs.

First, comparison was made of the fitted BTC with the observed mean concentrations with two hypothetical cases: (1) a deterministic dispersion assuming σ_D = 0; (2) a constant dispersivity λ = <D>/<v> across the field, as shown in Eqn 10 assuming σ_D = σ_v, with a perfect positive correlation, ρ_v, D = 1. Figure 7 shows the mean BTCs at 0.8 m described with the STM using fitted and hypothetical parameter values. Fitted (case C) and hypothetical (case A and B) parameter values used in Fig. 7 are presented in Table 5. The observation that the predicted peak values are smaller than the mean observations indicates that the STM with parameters derived from a local-scale CDE model may not fully capture the extreme non-equilibrium effects observed in this soil.

Fig. 7.  Predicted mean BTCs using the STM for hypothetical cases A (deterministic dispersion) and B (constant dispersivity), and using fitted parameters (case C).

Table 5.  Hypothetical parameter values for the stream tube model
Case C represents fitted values

The effect of σ_D and ρ_v, D on the mean solute concentration is almost negligible. In other words, the deterministic (σ_D = 0) and stochastic (σ_D = 0) dispersion process gave identical results. Similar results were found by Bresler and Dagan (1981), Amoozegar-Fard et al. (1982) and others, which also showed how little was the influence of fluctuations of local-scale dispersion on field-scale solute distribution. Hence, a constant dispersivity given by Eqn 10 is an acceptable simplification for practical applications of the STM model when little information is available on σ_D (Toride and Leij 1996).

In addition to σ_D and ρ_v, D, local-scale dispersion is often neglected for field-scale solute transport (e.g. Jury 1982; Destouni and Cvetkovic 1991). However, the relatively small Peclet numbers, <v>x/<D> = x/λ are around one for x > 0.05 m, as shown in Table 3, suggest that local-scale dispersion is important. We therefore evaluated the effect of <D> on the predicted field-scale BTC at 0.8 m depth assuming two additional hypothetical cases (case D and E in Table 5), whose local-scale dispersion is given by <D>/10 and <D>/100, respectively, where <D> is the fitted value. Figure 8 demonstrates that reduced dispersion leads to higher peaks and less spreading for relatively large <D>. Ignoring local-scale dispersion thus results in quite different field-scale solute distributions in cases where solute spreading is caused by physical non-equilibrium effects.

Fig. 8.  Predicted mean BTCs using the STM for hypothetical cases D (<D>/10) and E (<D>/l00) and using fitted parameters (case C).

A transport experiment on 1-m-long undisturbed soil columns collected along a transect in the field was used to evaluate the stochastic representation of the CDE in terms of predicting field-scale transport. Local-scale transport parameters were obtained for the classical CDE model and alternatively for the non-equilibrium MIM transport model. The latter gave a better description of the observed early breakthrough and the tailing. The effective dispersion coefficient, D_eff, in terms of the MIM model agreed relatively well at 0.8 m with D for the CDE. The assumption of the local equilibrium (LE) was nearly met as the travel distance approached the bottom of the column. Because of lateral solute mixing between mobile and immobile zones, estimated values for the CDE dispersion coefficient D were quite large.

Estimated transport parameters v and D varied significantly across the field with the coefficient of variation of 68% and 350%, respectively. A linear regression analysis between ln(v) and ln(D) resulted in D = 31v^1.92 with a positive correlation of ρ_v, D = 0.74; the overall field-scale dispersivity value is 0.656 m. The exponent close to 2 indicates that intra-aggregate diffusion is the dominant mechanism compared with the hydrodynamic dispersion for the solute spreading. Values for σ_v decreased from 0.55–0.71 in the Ap horizon to 0.36–0.38 in the C2, indicating that flow and transport processes were more heterogeneous in the top horizon than the two deeper layers.

Field-scale transport was assumed to be described by averaging of the local scale BTCs. When field-scale mean transport was predicted by means of STM using transport parameters obtained from the local-scale BTC, the predicted BTC agreed relatively well with the observed BTC as a whole. However, peak concentrations were underestimated and the tails were overestimated.

A comparison between the measured and predicted concentration variances revealed that for the first observation depth, both agreed well. Deeper in the soil, the STM gave consistently smaller values for the standard deviation.

In order to demonstrate the effect of local-scale dispersion on field-scale transport, mean observed BTCs were compared with predicted BTCs assuming (1) deterministic dispersion (σ_D = 0), and (2) a constant dispersivity across the field (σ_D = σ_v and ρ_v, D = 1). The results showed that the contribution from local-scale dispersion to field-scale solute spreading was insignificant and suggest that solute spreading at the field-scale may be approximated using a heterogeneous flow field (stochastic v) and a constant dispersion coefficient <D>. However, values of <D> could not be ignored, as evidenced by the considerably larger concentration peaks and reduced spreading when 10 or 100 times smaller values of <D> were used in comparison with the mean fitted <D>. The large values of <D> correspond to small values of x/<λ> = <v> x/<D>. Using the latter index we demonstrated local-scale dispersion cannot be neglected where physical non-equilibrium dominates the transport process; this is often true for structured or macroporous soil such as these, where large values of the local-scale dispersion can be expected to have an important effect on field-scale solute spreading.