Published in Vadose Zone Journal 3:1309-1321 (2004)
© 2004 Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA
ORIGINAL RESEARCH
Inverse Mobile–Immobile Modeling of Transport During Transient Flow
Effects of Between-Domain Transfer and Initial Water Content
J. Maximilian Köhnea,b,*,
Sigrid Köhnea,b,
Binayak P. Mohantya and
Jirka 
im
nekc
a Department of Biological & Agricultural Engineering, Texas A&M University, Scoates Hall, College Station, TX 77843-2117
b 1912 Vinewood, Bryan, TX 77802
c University of California, Riverside, Dep. of Environmental Sciences, 900 University Avenue, A135 Bourns Hall, Riverside, CA 92521
* Corresponding author (mkoehne{at}cora.tamu.edu)
Received 16 September 2003.
 |
ABSTRACT
|
|---|
Mobile–immobile models (MIM) have rarely been used for inverse simulation of measurements of variably saturated water flow and contaminant transport. We evaluated two MIM approaches with water transfer across the mobile and immobile regions either based on relative saturation (Se) differences, MIM(Se), or pressure head (h) differences, MIM(h), for inverse simulation of transient water flow and Br– transport in aggregated soil. Six undisturbed Ap soil columns (14.7-cm length and diameter) at wet, medium, and dry initial water contents were subjected to application of 0.005 L Br– solution of 1000 mg L–1 and subsequent irrigation of 1 cm h–1 for 3 h. Two similar irrigations were applied after 7 and 14 d. Measurements comprised pressure heads at depths of 2.8 and 12.8 cm, average soil column water contents, outflow, and effluent solute concentrations. This experimental information was used for simultaneous optimization of van Genuchten soil hydraulic parameters for mobile and immobile regions, the dispersivity, and the first-order rate coefficients for water and solute transfer. In total, eight parameters were estimated for MIM(Se) and 10 parameters for MIM(h). The inverse MIM approaches described the various hydraulic and transport data adequately. Physical nonequilibrium was more pronounced for wet and dry than for intermediate initial moisture, while initial moisture had no obvious effect on the total Br– lost. For wet and dry initial conditions, parameter estimates seemed fairly reliable, with the exception of the highly correlated saturated water contents in mobile and immobile regions. MIM(h) yielded parameters that appeared physically more consistent with observations, but required smaller time steps than MIM(Se) to overcome oscillations of pressure heads. Both MIM approaches were found to be suitable for inverse simulation of physical nonequilibrium transport during variably saturated flow.
Abbreviations: BTC, breakthrough curve • MIM, mobile–immobile models
|
INTRODUCTION
|
|---|
DURING THE PAST several decades, increasing concentrations of surface-applied agrochemicals in ground and surface waters have given rise to environmental concerns and have intensified research efforts to understand underlying transport processes (Flury, 1996). Various studies showed that in fine-textured soils, chemicals may be rapidly channeled to groundwater through preferential pathways caused by aggregation, cracking, decaying plant roots, and earthworm burrows (e.g., Flury et al., 1994; Mohanty et al., 1997; Lennartz et al., 1999). One manifestation of preferential flow in soil is physical nonequilibrium, that is, distinct local (centimeter to decimeter scale) differences with respect to flow velocities, pressure heads, water contents, and solute concentrations, as was demonstrated using laboratory experiments with a macroporous soil column (Köhne and Mohanty, unpublished data, 2004).
The initial or antecedent water content (i.e., the moisture status of a soil before solute and rain application) is likely to affect the infiltration of rain and of surface-applied chemicals. However, investigations on the effect of the initial water content on solute transport have yielded ambiguous results. For example, White et al. (1986) and Shipitalo et al. (1990) found that lower initial water contents favor physical nonequilibrium during fast solute leaching, while other studies indicated that high initial water contents enhance leaching (Tillman et al., 1992; Clothier and Green, 1994; Granovsky et al., 1993). In this study we analyzed the effect of the antecedent water content on observations and simulation results of transient water and Br– movement in an aggregated Ap soil.
Two- and multiregion transport models can be used for the risk assessment of physical nonequilibrium transport of contaminants into groundwater. The classical MIM comprises a mobile region subject to steady-state flow and convective–dispersive solute transport, along with diffusive solute exchange with a stagnant (immobile) region (e.g., van Genuchten and Wierenga, 1976). The original MIM is not easily applicable to field transport situations involving transient flow. More complex dual-permeability models (e.g., Jarvis and Leeds-Harrison, 1987; Jarvis et al., 1991; Chen and Wagenet, 1992; Gerke and van Genuchten, 1993; Kohler et al., 2001) and multiple permeability models (e.g., Gwo et al., 1995; Hutson and Wagenet, 1995) may be used to describe transient, variably saturated flow and transport in two or more mobile regions, coupled with water and solute mass transfer terms. The application of dual- or multipermeability models is hampered by their large numbers of model parameters. Only few studies have been presented recently where model parameters were obtained by independent measurement (e.g., Köhne et al., 2002a, 2002b; Logsdon, 2002; Castiglione et al., 2003) or by inverse estimation (e.g., Schwartz et al., 2000; Kätterer et al., 2001; Dubus and Brown, 2002; Roulier and Jarvis, 2003).
The MIM approach in several studies has been extended to describe transport under variably saturated flow conditions (e.g., Russo et al., 1989; Zurmühl and Durner, 1998; Clothier et al., 1998; Larsson et al., 1999). However, advective transfer between regions was considered only in the study of Larsson et al. (1999). The most comprehensive MIM accounting for water and solute transfer between regions at variable saturations was recently presented by imnek et al. (2003). Here we will extend their approach by utilizing inverse parameter estimation to analyze experimental data for physical nonequilibrium solute transport under transient flow conditions.
An important aspect of this study is the effect of the advective component in the MIM solute transfer term on the inverse transport simulation results. Solute transfer has been identified as crucial for accurate prediction of nonequilibrium transport (e.g., Flühler et al., 1996; Wilson et al., 1998; Gerke and Köhne, 2004). Advective transfer directly depends on water transfer that in its turn can be either assumed to be proportional to the differences in saturation (e.g., Jarvis et al., 1991; Kohler et al., 2001) or the pressure head (e.g., Gerke and van Genuchten, 1993; Gwo et al., 1995) across regions. However, the effect of having two different advective terms on MIM simulation results for solute transport has not been evaluated to date.
The objectives of this work were (i) to study the effect of the initial soil moisture on physical nonequilibrium water flow and solute transport as observed and simulated for six aggregated Ap soil columns and (ii) to evaluate two MIM approaches in terms of their suitability for inverse simulation of measurements of variably saturated flow and physical nonequilibrium Br– transport. The approaches assume water and advective transfer either based on differences between the immobile and mobile regions in pressure heads, approach MIM(h), or saturations, approach MIM(Se).
|
MATERIALS AND METHODS
|
|---|
Six undisturbed soil columns (14.7-cm diameter and 14.7-cm length) were taken from the Ap horizon of a loamy soil in a field near Kiel (Northern Germany). The Ap soil consisted of 2.3, 43.8, and 53.9% (w/w) clay, silt, and sand, respectively, and showed subangular and crumbed peds of 0.5- to 1-cm diameter. Total C and CO–3–C of the soil materials were determined by dry-ashing at 1200°C and measuring the volume of CO2 produced after reaction with HCl, respectively. Organic C (1.83%) was calculated as the difference between total and CO–3–C. A soil pH of 6.7 was measured in 0.01 M CaCl2 solution at a soil/solution ratio of 1:2.5. The average bulk density, 
(1.37 g cm–3), and porosity, P (0.47), were determined from 20 100-cm3 samples. Assuming a mean particle density (p) of 2.6 g cm–3, porosity was calculated as P = 1 – p–1.
Tracer displacement experiments were conducted to study the effect of different initial water contents on transport behavior during transient flow. We used two replicates for each of three different initial water saturations, henceforth called "wet" (Columns Ap1, Ap2), "medium" (Ap3, Ap4), and "dry" (Ap5, Ap6). Conditions of selected experiments are summarized in Table 1. Rainfall at the top of the column was simulated by means of an automated sprinkler, while constant suctions in a range from 15 to 60 hPa were applied to the columns through the lower boundaries. Pressure heads at two depths (2.8 and 12.8 cm) were recorded using tensiometers. Pressure heads were assumed to represent the mobile region since the tensiometers' ceramic tip of 5-cm length was intersecting several interaggregate sections. Columns Ap1, Ap3, Ap4, and Ap6 were placed on automatically weight registering scales to monitor depth-averaged water contents. At the end of the transport experiment, the column was weighed and the soil was dried at 105°C to determine the gravimetric water content. The final volumetric water content was calculated by multiplying the gravimetric value with the bulk density, which was obtained by dividing the dry weight of the soil by the column volume. This final volumetric water content and the balance registrations were used to back-calculate depth-average volumetric water contents during the experiments. Since only four scales were available for automated measurements, for the Ap2 and Ap5 columns only the initial and the maximum water contents could be determined from the initial weights and the weights measured after the first irrigation.
The transport experiment was initiated by applying 5 cm3 (0.03 cm) of a KBr solution with a concentration of 1000 mg L–1 Br on the soil surface within 10 min by means of a chromatography syringe. Starting 20 min later, 3 cm of rain was sprinkled onto the columns at a rate of 1 cm h–1. Two more similar irrigation events were applied in 7-d intervals. The effluent of each column was collected in 16- to 18-cm3 volume fractions and analyzed for Br– (flux type) concentrations. For a complete mass balance, soil aliquots were taken from 3-cm-thick layers of the Ap columns to determine water content and the mass of Br– residing in these layers. The latter was determined by extracting 50 g dry soil with 50 cm3 H2O by overhead shaking for 1 h. Subsequently, the batch samples were rotated in a centrifuge and the supernatant solution was analyzed for Br– (resident type) concentration. Bromide concentrations were determined using an ion-chromatography system with a conductivity detector (GAT Wescan, Kontron Instruments, Exchingen, Germany). The solid phase was a 100-mm long anion cartridge (Methrom, Super-Sep 6.1009.000, Deutsche Methrom GmbH & Co. KG, Filderstadt, Germany) together with a 20-mm guard-cartridge (Bischoff chromatography, part 6302137, Metrohm AG, Herisau, Switzerland). The mobile phase was a 2.5-mM phthalic acid eluant (pH 4.2) that was degassed via an automatic degasser (Degasys, Chrom Tech, Apple Valley, MN). The detection limit of the system was 0.2 mg L–1 Br. Additional information about the experiments can be found in Meyer-Windel (1998). Bromide and herbicide transport in the wet Ap1 and Ap 2 soil columns have been discussed in detail by Meyer-Windel et al. (1999).
|
MODELS
|
|---|
The MIM approach assumes partitioning of the soil water content, 
, into a mobile region water content, m, and an immobile region water content, im (van Genuchten and Wierenga, 1976):
 | [1] |
The Richards' equation is used to describe one-dimensional transient water flow in the mobile region, while a source–sink term accounts for advective water exchange with the variably saturated immobile region (slightly modified after imnek et al., 2003):
 | [2] |
where t (T) is time, z (L) is vertical distance positive downward, K (L T–1) is the hydraulic conductivity as a function of pressure head in the mobile region, hm (L), and 
w (T–1) is the water transfer rate between the mobile and the immobile region.
We used two alternative approaches for calculating the water transfer term w. First, w was assumed to be proportional to the difference in effective saturation, Se (0 
Se 1), between the mobile and immobile regions, MIM(Se)
 | [3] |
where 
(T–1) is a first-order water transfer rate coefficient and where Sem (Seim) is the effective saturation of the mobile (immobile) region defined as
 | [4a] |
 | [4b] |
where m (im) (L3 L–3) denotes the mobile (immobile) water content with residual and saturated constraints of m,r (im,r) and m,s (im,s), respectively. To describe m in Eq. [4a], the following form of van Genuchten's (1980) equation was assumed
 | [5] |
where m and nm are parameters defining the shape of the retention function for the mobile region.
Alternatively, w was assumed to be proportional to the difference in pressure heads between the mobile and immobile regions, MIM(h).
 | [6] |
where * (L–1 T–1) is a first-order water transfer rate coefficient. The MIM(h) approach requires a retention function of the immobile region analogous to Eq. [5]. Water transfer terms based on pressure head differences (Eq. [6]) are commonly perceived as being more physically realistic than water transfer terms based on saturation differences (Eq. [3]), since the difference in pressure heads is the actual driving force (imnek et al., 2003). However, note that the first-order term (Eq. [6]) approximates the real, unknown pressure head profile as a difference between only two values (i.e., hm and him). Due to this first-order simplification, no substantial conceptual differences exist between first-order terms based on saturation and pressured head. In fact, Eq. [6] and [3] would become identical if in Eq. [3] the effective saturations in both mobile and immobile regions were calculated using van Genuchten (1980) relations (Eq. [5]) with the same parameter values for 
and n. Hence, the main difference between these two approaches is that while MIM(Se) assumes similarity in the retention curves of the mobile and immobile regions, the MIM(h) approach allows different retention curves for each region.
To assess the potential effect of initial moisture on the mobile water content, the fraction of mobile water at saturation, ßs, was calculated as
 | [7] |
The MIM formulation for nonreactive solute transport comprises a convection–dispersion equation for the mobile region coupled via a sink–source term to the immobile region (imnek et al., 2003):
 | [8a] |
 | [8b] |
 | [8c] |
where cm and cim (M L–3) are solute concentrations in the mobile and immobile domains, respectively, qm (L T–1) is the flux density in the mobile region, s (M L–3 T–1) is the solute transfer rate between mobile and immobile region, with w calculated either using Eq. [3] for MIM(Se) or Eq. [6] for MIM(h), s (T–1) is the first-order solute transfer rate coefficient, and Dm (L2 T–1) is the dispersion coefficient in the mobile region defined as
 | [9a] |
in which
 | [9b] |
where 
m (L) is the dispersivity of the mobile region, D0 (L2 T–1) is the molecular diffusion coefficient, 
m is the tortuosity coefficient of the mobile domain evaluated according to Millington and Quirk (1961), and m,s is the saturated water content of the mobile domain.
The above MIM equations cannot be solved analytically. They were implemented in the numerical HYDRUS-1D computer model (imnek et al., 2003). To ensure an accurate numerical solution, the numerical grid consisted of 200 elements of only 0.07-cm length, while the maximum time-step was restricted to 0.001 d in a time-step adaptive Picard iteration solution scheme. The effect of using different spatial and temporal resolutions on numerical simulation results was analyzed (see Results section).
The initial condition for water flow was represented by a pressure head profile linearly interpolated between measurements of tensiometers installed at the 2.8- and 12.8-cm depth. We assumed initial equilibrium with identical pressure heads in the mobile and immobile regions. The upper flow boundary condition was a variable flux-type condition equal to the time-variable rates of irrigation (22 cm d–1) and evaporation (0.3 cm d–1) set according to the experimental conditions. Upon ponding this boundary condition may evolve into a program-controlled head-type condition. For the lower boundary, an unsaturated seepage face at a user-specified pressure head (constant head imposed at the lower boundary of soil column) was implemented in HYDRUS-1D. The initial condition for solute transport was set to zero. The upper and lower boundaries for solute transport were represented by specified flux and (zero) gradient conditions, respectively.
Parameter Estimation
For describing mobile–immobile water flow using Eq. [2], [3], [4], and [5], MIM(Se) requires eight parameters (Ks, m,r, m,s, m, nm, im,r, im,s, and ). The MIM(h) approach (Eq. [2], [4], [6]) additionally requires im and nim. To reduce the number of unknowns, m,r and im,r were fixed to values of zero consistent with r values of zero fitted to water retention and hydraulic conductivity data for all Ap soil columns (Meyer-Windel, 1998). This limited the total number of unknown water flow parameters to six [MIM(Se)] and eight [MIM(h)], respectively. To calculate solute transport using Eq. [8] and [9], both inverse MIM approaches additionally require estimation of s and , whereas D0 (1.797 cm2 d–1) was assumed to be a known value for Br (Atkins, 1990). Hence, full descriptions of mobile–immobile water flow and solute transport assuming zero residual water contents require 8 and 10 parameters for MIM(Se) and MIM(h), respectively.
The inverse parameter estimation was performed by minimization of the objective function 
(imnek et al., 1998):
 | [10] |
where m is the number of different sets of measurements, n represents the number of observations in a particular measurement set, Oj(z,ti) are observations such as pressure heads in the mobile region at the 2.8- and 12.8-cm depths, average soil column water contents (averaged over the column depth and over mobile and immobile regions), cumulative water fluxes across the lower boundary, and effluent concentrations at time ti at location z, Ej(z,ti,b) are the corresponding estimated space-time variables for the vector b of optimized parameters (im,r, im,s, m,r, m,s, m, nm, Ks, , , s, and additionally im, nim for MIM(h), and vj and wi,j are weighting factors associated with a particular measurement set or point, respectively. We assumed in this study that the weighting coefficients wi,j in Eq. [10] are equal to one; that is, the variances of the errors within a particular measurement set are the same. Only data that were measured at larger time intervals (and hence were underrepresented as compared with the more frequent measurements) were assigned larger weights wi,j, calculated as the ratio of their measurement time interval to the average time interval. The weighting factors vj were obtained in two steps. An initial value for vj was calculated for each inverse simulation as (Clausnitzer and Hopmans, 1995)
 | [11] |
which assumes that vj is inversely related to the variance 
2j within the jth measurement set and to the number of nj measurements of the set. However, using Eq. [11] was not sufficient to yield equal contributions of the four measurement sets (i.e., water contents, pressure heads, water outflow, and Br– concentrations in the outflow) to the optimized objective function (Eq. [10]) and resulted in poor simulations of those measurement sets that were underrepresented in the objective function. Therefore, the weight vj of every measurement set was fine-tuned by repeated manual calibration in consecutive inverse simulation runs to achieve equal individual contributions of the four measurement sets to the overall sum of squared residuals.
The Levenberg–Marquardt algorithm was used to minimize the objective function (Eq. [10]) and to calculate confidence intervals and correlation matrices for the parameters (imnek et al., 1998). Our inverse modeling strategy relied on simultaneous estimation of hydraulic and transport parameters from combined water flow and solute transport information. This approach takes advantage of all information contained in the data, since concentrations are a function of water flow also (Medina and Carrera, 1996). For coupled single-domain flow and transport, simultaneous estimation of hydraulic and transport parameters was found to yield smaller estimation errors and confidence intervals for model parameters than sequential inversion of hydraulic properties from water content and pressure head data followed by inversion of transport properties from concentration data (imnek et al., 2002).
To facilitate direct comparisons of the goodness-of-fit for different data sets, the RMSE between observations (Oi) and model estimates (Ei) for a particular data set was normalized with the arithmetic mean of observations, 
, to yield the coefficient of variation (CV, %) as
 | [12] |
An average coefficient of variation, 
, over the four measurement sets of a particular Ap column was calculated as a statistical measure of the overall agreement between model and observations as follows:
 | [13] |
where CVk is the coefficient of variation (Eq. [12]) for measurement set k, and m denotes the total number of CV values or data sets (in this study m = 4).
To enable a more physically based interpretation of the simulated MIM transport processes, we analyzed the degree of physical (non-)equilibrium for different initial water contents. For steady-state flow, an index frequently used for defining physical nonequilibrium is the Damköhler number (Damköhler, 1936; Bahr and Rubin, 1987; Vanderborght et al., 1997; Michalak and Kitanidis, 2000). However, the Damköhler number cannot be used for transient flow conditions. In this study, we introduce an empirical physical equilibrium (PE) index to quantify the degree of physical equilibrium during transient flow by describing the combined effects of advective and diffusive first-order mass transfers. The PE index (%) is defined as
 | [14] |
where (T–1) (Eq. [3], [6]) is a first-order water transfer rate coefficient, s (T–1) (Eq. [8]) is a first-order solute transfer rate coefficient (Table 1), and max and s,max are maximum parameter constraints for and s, respectively, indicative of water and solute transfer equilibrium. Theoretically, instantaneous equilibrium exists only if and s are infinite. Practically, parameter constraints max and s,max could be operationally defined, for example, by using values that force water contents and solute concentrations in mobile and immobile regions to reach 99% equilibrium within the experimental time scale. However, such a definition would require individual simulations for every set of parameters, and initial and boundary conditions to identify max and s,max. Using a more practical approach, max and s,max in this study were both simply set equal to a large value of 10 d–1; hence, max+ s,max equaled 20 d–1. This was found to effectively ensure equilibration of water contents or pressure heads and concentrations within the experimental time scale. For the parameter combinations tested in this study, increasing values for or s to above 10 d–1 did not substantially change the results (not further shown here). By definition, the PE index is inversely related to the degree of physical nonequilibrium and varies between 0% (perfect nonequilibrium, no mass transfer) to 100% (inter-region equilibrium due to very fast mass transfer).
|
RESULTS AND DISCUSSION
|
|---|
Water Flow Results
Initially Wet Ap Soil Columns (Ap1, Ap2)
Experimental water flow data and model simulations for the initially wet soil columns are shown in Fig. 1
. Water contents, pressure heads in the mobile region at depths of 2.8 and 12.8 cm, and the cumulative outflow exhibited a recurrent pattern of increasing values during infiltration and decreasing or constant (cumulative outflow) values during redistribution (Fig. 1a–1h). For the Ap2 column, only the initial water content and a value after the first irrigation event were known (Fig. 1b). Observations were comparable for Ap1 and Ap2, reflecting similar water flow dynamics in the initially wet replicate soil columns. According to Fig. 1, the MIM(Se) and MIM(h) approaches approximated observed water contents, pressure heads, and cumulative outflow equally well. Since infiltration occurs on a smaller time scale than drainage, Fig. 2
shows a detailed view of the short-term responses to one irrigation event. The third irrigation event is presented here to illustrate the maximum deviations between model results and data, as opposed to the first event, where known initial conditions resulted in a closer match with the data (not shown). The detailed view still shows reasonable agreement between hydraulic data measured during the third irrigation event and both MIM approaches (Fig. 2).
View larger version (38K):
[in this window]
[in a new window]
|
Fig. 1. Water flow data and simulation results obtained with MIM(Se) and MIM(h) for the initially wet columns Ap1 (left) and Ap2 (right). (a, b) Water contents, (c, d) pressure heads in the mobile region at depths of 2.8 cm and (e, f) 12.8 cm, and (g, h) the cumulative outflow.
|
|
Medium Initial Water Content in Ap Soil Columns (Ap3, Ap4)
For the replicate Ap3 and Ap4 soil columns of medium initial water content, hydraulic dynamics as reflected by water contents, mobile region pressure heads at the 2.8-cm depth, and cumulative outflow were again comparable (Fig. 3)
. All observations were matched satisfactorily by MIM(Se) and MIM(h), both on average during the experimental duration of 20 d (Fig. 3) and during single irrigation events (not shown).
View larger version (36K):
[in this window]
[in a new window]
|
Fig. 3. Water flow data and simulation results obtained with MIM(Se) and MIM(h) for the initially medium wet columns Ap3 (left) and Ap4 (right). (a, b) Water contents, (c, d) pressure heads in the mobile region at the 2.8-cm depth, and (e, f) the cumulative outflow.
|
|
Initially Dry Ap Soil Columns (Ap5, Ap6)
The Ap5 and Ap6 soil columns were characterized by low initial water contents (see Table 1). Hydraulic dynamics were comparable for Ap5 and Ap6 (Fig. 4)
, except for differences in water content variations during the first 7 d of the experiment (no water contents could be measured for Ap5 after 7 d) (Fig. 4a). However, it is likely that these apparent differences were caused by less reliable measurements of Ap5 water contents. Except for the water contents for Ap5, hydraulic data were reasonably well described with both MIM(Se) and MIM(h) during the 20-d experiment (Fig. 4a–4f) and during the irrigation events (detailed view not shown).
View larger version (33K):
[in this window]
[in a new window]
|
Fig. 4. Water flow data and simulation results obtained with MIM(Se) and MIM(h) for the initially dry columns Ap5 (left) and Ap6 (right). (a, b) Water contents, (c, d) pressure heads in the mobile region at the 2.8-cm depth, and (e, f) the cumulative outflow.
|
|
Comparison of Hydraulic Data, Simulations, and Parameters
The maximum water contents that were reached after the irrigations declined in the order: Ap1,2 (wet); Ap3,4 (medium); Ap5,6 (dry) (see Fig. 1a, 1b, 3a, 3b, 4a, and 4b), indicating incomplete saturation during irrigation of the initially drier soil columns. The MIM(h) simulation results for the saturated water content, s, summed for the mobile and immobile regions (s,m + s,im), reflected the decreasing order of maximum water contents with decreasing initial water content, whereas respective s values obtained with MIM(Se) did not consistently decrease with drier initial conditions (Table 2). The s values obtained with MIM(h) were hence more consistent with the observations, even if statistical CV values indicated that MIM(Se) and MIM(h) on average described hydraulic observations of water contents, pressure heads, and outflow equally well (Table 3).
View this table:
[in this window]
[in a new window]
|
Table 3. Coefficient of variation between measurement sets and model estimates for cumulative outflow, qcum(t), Br– concentration in outflow, cout(t), pressure heads in the mobile region at depths of 2.8 and 12.8 cm, hm(z,t), and water contents averaged over the column depth, (t).
|
|
For the initially wet columns Ap1 and Ap2, independently measured and inverse estimates of retention functions, (h), and hydraulic conductivity functions, K(h), are displayed in Fig. 5
. Measurements of (h) and K(h) were obtained at static and steady-state conditions, respectively (Meyer-Windel, 1998). The estimated bimodal functions represent the summation (Durner, 1993) of van Genuchten subcurves based on parameters obtained by the inverse numerical solution of MIM(h) and MIM(Se) using all available transient flow and transport data. Good agreement between independent measurements and predictions of (h) and K(h) was achieved particularly for MIM(h), while the inverse MIM(Se) approach overestimated s for Ap1 (Fig. 5).
View larger version (18K):
[in this window]
[in a new window]
|
Fig. 5. (a) Measured water retention and (b) hydraulic conductivity data and bimodal Durner (1993) functions estimated inversely using MIM(Se) and MIM(h) for the Ap1 and Ap2 columns.
|
|
Incomplete saturation during infiltration may indicate the presence of immobile regions. Hence, one may assume that the initially dry columns that reached smaller maximum saturations also had a smaller fraction of mobile water at saturation, ßs (Eq. [7]). However, ßs assumed values between 11 and 41%, independent of the initial water content (Table 4). Moreover, the actual fraction of mobile water at partial saturation may have been variable during transient flow. Overall, the saturation differences were apparently too small (compared with other causes for variation) to exert a notable effect on ßs.
View this table:
[in this window]
[in a new window]
|
Table 4. Estimated fraction of mobile water at saturation, ßs, dispersivity of the mobile region, transfer rate coefficients, and physical equilibrium index (PE).
|
|
Among the fitted hydraulic mobile–immobile model parameters, only the MIM(h) parameters im, nim, and s (s,m + s,im) showed a tendency to decrease for drier initial conditions (Table 2), thus indicating an effect of initial water content on the water flow process. The remaining parameters (m, nm, Ks, s,m, s,im) neither depended on the invoked model [MIM(h) vs. MIM(Se)] nor on the initial water content. The water transfer coefficients were found to strongly affect the solute transport simulation, and therefore will be included in the discussion below.
Solute Transport Results
For all Ap columns, the breakthrough curves (BTC) in terms of the fraction of applied Br– concentration vs. cumulative outflow are shown in Fig. 6
. The BTCs exhibited an asymmetric shape with a long tail typical of physical nonequilibrium. The low peaks of only 0.004 to 0.008 (4–8 mg L–1) of the applied concentration (
1000 mg L–1, Table 1) indicate that only mild physical nonequilibrium existed for all columns, as opposed to preferential flow, where much higher concentration peaks would be expected. The initially wet column Ap1 and the initially dry Ap5 and Ap6 columns exhibited an early rise, whereas the Ap2 column and the initially medium wet Ap3 and Ap4 columns displayed a somewhat delayed initial rise. Note that the Ap2 column had in fact an initial water content somewhere between those of the initially wet Ap1 and medium wet Ap3 and Ap4 columns (Table 1). Overall, the MIM(Se) and MIM(h) described the BTC reasonably well. According to CV values (Table 3), Br– concentrations of Ap1, Ap2, Ap3, and Ap4 were better matched with MIM(Se), whereas the MIM(h) approach better described the BTC of the initially dry Ap5 and Ap6 columns.
View larger version (40K):
[in this window]
[in a new window]
|
Fig. 6. Relative Br– concentration vs. cumulative outflow: data for all Ap columns and simulation results obtained with MIM(Se) and MIM(h).
|
|
The fraction of mobile water (discussed above), the dispersivity for the mobile region of the fitted two-region models, m, the rate coefficients for transfer of water, (*), and solute, s, and the PE index are given in Table 4. The m values varied between 1 mm and 3.75 cm, which is within the range of m values reported by others for column scale studies (e.g., Schwartz et al., 2000). Estimated values for , *, and s for MIM(h) were smaller for Ap1, Ap5, and Ap6 than for Ap2, Ap3, and Ap4 (Table 4). This indicates restricted advective and diffusive transfer (nonequilibrium) for initially wet and dry as compared with initially intermediate moisture conditions (equilibrium). For MIM(h), the s value actually reached its maximum constraint (s = 10) for Ap2, Ap3, and Ap4 (Table 4), indicating enhanced diffusive solute transfer for concentration equilibration across regions. For Ap2, caution must be taken to interpret s with respect to solute transport processes since the BTC of Ap2 was only described adequately for the first 2 cm of outflow, whereas the low tailing concentrations could not be described (Fig. 6). For MIM(Se), optimized values for water and solute transfer rate coefficients were not as consistent as for MIM(h). For example, MIM(Se) gave a large value for Ap6 and a small for Ap5 for similar initial water contents (Table 4).
The PE values obtained with MIM(h) were small for the high (Ap1) and low (Ap5 and Ap6) initial moisture conditions, while PE values were large for medium (Ap3, Ap4) initial moisture, indicating more physical nonequilibrium for initially wet and dry than for intermediate conditions. For MIM(Se), PE showed same tendencies, but differences were not as pronounced among PE values for the various water contents (Table 4). The simulations with low PE values (Ap1, Ap5, Ap6) displayed secondary concentration peaks typical of physical nonequilibrium; these were caused by diffusive back transfer from the immobile into the mobile region during no-flow periods. Secondary peaks were also observed for the experimental BTC of Ap2. Considered together, the experimental data, the model results, and the PE values indicate a higher degree of physical nonequilibrium for initially wet and dry conditions than for medium initial soil moisture. Note, however, that different initial moisture contents did not systematically affect the total mass loss of Br for the Ap columns under study, which amounted to 62, 47, 51, 73, 41, and 45% of the applied Br mass for Column Ap1, 2, 3, 4, 5, and 6, respectively.
We now discuss differences between MIM(h) and MIM(Se) on the basis of simulation results for Ap1, Ap4, and Ap6 after 70 min (i.e., 60 min after solution application and 50 min after initiating the first rain event), and we suggest a mechanistic interpretation for the transport at different initial water contents. Effective saturations (Se) as simulated with MIM(Se) (Fig. 7a)
indicate that the wetting front (Se = 1) in the mobile region at 70 min almost reached the bottom of the Ap1 column (14.7 cm), followed by Ap4 and Ap6. Behind the respective wetting fronts, saturation differences, assumed as the driving force for water transfer in MIM(Se), were smaller for the initially wet Ap1 than for initially intermediate wet Ap4 and initially dry Ap6 (Fig. 7a). Additionally, for MIM(Se) the water transfer coefficient, , increased in the order: Ap1, Ap4, and Ap6 (Table 4). Accordingly, water transfer into the immobile region was smallest for Ap1, larger for Ap4, and largest for Ap6 (Fig. 7b). For the initially dry Ap6, water laterally imbibed into the immobile region, thereby slowing down vertical water flow in the mobile region. By contrast, mobile water in the wet Ap1 could percolate faster due to less horizontal transfer into the relatively wet immobile region. The differences in Se values between the mobile and immobile regions of the Ap6 column (Fig. 7a) were two orders of magnitude smaller than similar differences in h values (Fig. 7c). However, water transfer calculated using both MIM approaches was comparable for the Ap6 column (Fig. 7b, 7e), since fitted with MIM(Se) was 50 times that of * fitted with MIM(h) (Table 4). For Ap1 and Ap4, estimated values for and * were of the same order of magnitude (Table 4), but water transfer in the lower profile of Ap1 and Ap4 was larger for MIM(h) (Fig. 7e) than for MIM(Se) (Fig. 7b).
View larger version (29K):
[in this window]
[in a new window]
|
Fig. 7. Depth profiles at 70 min characterizing water flow for the Ap1, Ap4, and Ap6 columns as simulated with MIM(Se) (top): (a) water saturation and (b) water transfer rate (w), and as simulated with MIM(h) (bottom): (c) pressure heads for Ap6 and (d) for Ap1 and Ap4, and (e) water transfer rate. Symbols: mo, mobile region; im, immobile region.
|
|
Figure 8
shows simulated concentration profiles and solute transfer rates at 70 min. For MIM(Se), differences between Br– in the mobile and immobile regions were somewhat more pronounced for Ap1 than for Ap4 (Fig. 8a, 8b), corresponding to a larger overall (advective–diffusive) Br– transfer rate s for Ap4 (Fig. 8c). For MIM(h), the Ap4 concentration profiles (Fig. 8d, 8e) and Br– transfer rate (Fig. 8e) indicate an advanced state of equilibration throughout the profile. The above results support the interpretation of having physical nonequilibrium for high and low initial moisture and transport close to equilibrium at medium initial moisture contents.
The following explanation can be given for solute transfer at 70 min being directed into the immobile region, except in the upper 2 cm (Fig. 8c, 8f). The advective transfer component was always directed into the immobile region. At the end of the tracer application (at 10 min), diffusive mass transfer in the uppermost 2 cm of the soil profile was also directed into the immobile region, thus leading to a large overall solute transfer rate, s (Fig. 9a, 9b)
. The first irrigation subsequently flushed the Br– downward through the mobile region, thus leaving higher concentrations in the immobile region, which reversed the concentration gradient in the upper 2 cm and directed diffusive transfer back into the mobile region after 70 min (Fig. 8c, 8f).
View larger version (15K):
[in this window]
[in a new window]
|
Fig. 9. Bromide transfer rate (s) into the immobile (im) region in the 0- to 2-cm depth layer at the end of the Br– application (10 min) for the Ap1, Ap4, and Ap6 columns. Simulation results obtained with (a) MIM(Se) and (b) MIM(h).
|
|
Parameter Uniqueness and Overall Model Evaluation
Having a large number of flow and transport parameters—8 MIM(Se) and 10 MIM(h) parameters—that needed to be identified at the same time evokes the question of whether or not the estimated parameter values were unique. Unfortunately, the typical analysis of parameter uniqueness using two-dimensional response surfaces of the objective function as a function of parameter pairs (e.g., imnek and van Genuchten, 1996) is not readily possible for such a large number of parameters. The number of parameter pairs to be analyzed using response surfaces in our study would be 28 [MIM(Se)] and 45 [MIM(h)], respectively. Moreover, these response surfaces as two-dimensional cross sections through the 8- or 10-dimensional parameter space may not accurately reflect the true behavior of the objective function. To assess the reliability of the presented parametric values, we (i) evaluated confidence intervals for the estimated parameter values, (ii) analyzed correlation matrices, and (iii) studied how different initial parameter values affected final parameter results.
Confidence intervals of the parameters were in general relatively small. Typically, an approach with small CV values like MIM(Se) for Ap1 (Table 3) yielded smaller confidence intervals (Table 2). For both MIM(Se) and MIM(h), smaller confidence intervals for the water transfer rate coefficients were obtained for initially wet (Ap1) and dry (Ap5,6) conditions than for initially medium (Ap3,4) and medium wet (Ap2) conditions (Table 4). This indicates that for initially wet and dry conditions, where physical nonequilibrium was more pronounced, parameter estimation using the MIM approaches was more reliable than for medium initial moisture conditions.
Correlation matrices for Ap1, Ap4, and Ap6 are presented in Table 5 for the MIM(Se) approach and in Table 6 for the MIM(h) approach. In Tables 5 and 6, moderately large correlation coefficients (italics denote |r| > 0.5) and large coefficients are identified (parentheses denote |r| > 0.75). Correlations between parameters obtained with MIM(Se) in general appear relatively small: for Ap1, 96% are below 0.75 and 82% are smaller than 0.5 (Ap4: 89% < 0.75, 82% < 0.5; Ap6: 93% < 0.75, 71% < 0.5) (Table 5). Particularly for Ap1 and Ap6, small correlations for water and solute transfer coefficients indicate reliable estimates of and s (Tables 5 and 6). By contrast, for Ap4 two large correlations (|r| > 0.75) obtained for and s indicate that for intermediate initial moisture, and s may have been less well defined. For both MIM(Se) (Table 5) and MIM(h) (Table 6), correlations between m,s and im,s were always large, indicating that their simultaneous estimation did not yield reliable and unique parameter values. This in turn explains that no consistent effect of initial water content on m,s and im,s could be found.
View this table:
[in this window]
[in a new window]
|
Table 5. Correlation, r, between MIM(Se) parameters for inverse simulation of the Ap1, Ap4, and Ap6 transport studies.
|
|
View this table:
[in this window]
[in a new window]
|
Table 6. Correlation, r, between MIM(h) parameters for inverse simulation of the Ap1, Ap4, and Ap6 transport studies.
|
|
Experimental methods for independent, approximate experimental determination of m and im have previously been proposed (e.g., Clothier et al., 1995). However, we did not utilize these methods here since we intended to elucidate the full capacity of the inverse parameter identification method by analyzing the extent to which parameters can be identified at the same time without resorting to additional, time-consuming experimental methods. Apart from the correlation between m and im, correlations between parameters obtained with MIM(h) were generally moderate for Ap1 and Ap6: 96% were below 0.75, and 71% (Ap1) and 82% (Ap6) were less than 0.5. For Ap4, high correlation values for nim, im, nm, and m indicate that these parameters may not have been uniquely identified for intermediate initial water contents. Correlations for * and s displayed mostly moderate values, indicating proper parameter identification as required for calculations of the PE index.
Simulations using MIM(Se) and MIM(h) were repeated for Ap1, Ap4, and Ap6 by increasing each initial parameter by 30% of its previous value. In general, final optimized parameters were within ±10% of their prior optimized values, thus indicating reliable parameter identification for Ap1 and Ap6. As an example, the simulation for Ap1 using MIM(Se) with 30% higher initial parameter values gave the following final estimates (previous optimization result are given in parentheses): m,s = 0.128 (0.115); m = 0.09 (0.07) cm–1; nm = 1.22 (1.25); Ks = 25.8 (22.0) cm d–1; im,s = 0.312 (0.310); = 0.41 (0.40) d–1; m = 1.17 (1.08) cm; s = 2.29 (2.20) d–1, while for MIM(h) the following final estimates were obtained for Ap1: m,s = 0.150 (0.147); m = 0.059 (0.034) cm–1; nm = 1.25 (1.25); Ks = 22.5 (22.5) cm d–1; im,s = 0.242 (0.214); im = 0.021 (0.016) cm–1; nim = 1.68 (2.00); = 0.06 (0.19) d–1; m = 2.04 (2.26) cm; s = 1.68 (1.34) d–1. For Ap4, parameter identification using MIM(Se) with 30% higher initial parameter values was somewhat less accurate than for Ap1 and Ap6: m,s = 0.178 (0.128); m = 0.049 (0.048) cm–1; nm = 1.29 (1.26); Ks = 34.9 (21.6) cm d–1; im,s = 0.216 (0.248); = 0.50 (1.32) d–1; m = 2.35 (1.91) cm; s = 10 (10) d–1. The inverse solution of MIM(h) for input parameters increased by 30% was not completed within 12 d, at which time we aborted the calculation.
The inverse MIM approaches overall gave acceptable representations of the data. Most of the 8 [MIM(Se)] and 10 [MIM(h)] flow and transport parameters could be reasonably well identified at the same time for initially wet and dry conditions. The parameters m,s and im,s could not be reliably identified at the same time and hence should be determined independently. For medium initial water contents, large correlations between several parameters and dependence of the final on the initial parameter values demonstrated that simultaneous fitting did not represent a well-posed problem for both MIM approaches. This is because medium initial moisture conditions produced near-equilibrium transport where MIM parameters are less well identifiable due to overparameterization of the model. It has previously been found that for soils with spherical aggregates, MIM parameters are not always identifiable from miscible displacement experiments (Vanderborght et al., 1997). The problem of parameter nonuniqueness is often caused not only by the large number of optimized parameters, but also by the lack of sufficient experimental information. Parameter identifiability can therefore be increased either by reducing the number of optimized parameters, or by including appropriate additional experimental information (Hopmans and imnek, 1999), such by using simultaneously information about multiple variables. Note that we defined the objective function using four different sets of variables (i.e., pressure heads, water contents, cumulative water flow, and solute concentrations). This combination of variables yielded a relatively good identification of most parameters in our study.
Numerical Discretization Effects on Simulation Results
For the above comparison of inverse MIM(h) and MIM(Se) approaches, we always used the same small values for the maximum time step dtmax (0.001 d) and the space increment dz (0.074 cm). Numerical pressure head oscillations for MIM(h) occurred when the time step adaptive iteration scheme was not constrained by a small value for dtmax and calculated larger time steps during drainage. A forward simulation for the Ap5 column shows that numerical oscillations in the pressure head, h, were more effectively reduced by decreasing dtmax than by decreasing dz (Fig. 10)
. Keeping dz at 0.5 cm (29 finite elements), the oscillations in h could be eliminated by reducing dtmax from 0.03 to 0.01 d (Fig. 10a, 10b, 10c). Constraining dtmax to 0.03 d and reducing dz to 0.037 cm (400 elements) did not eliminate oscillations (Fig. 10d, 10e, 10f). Even for the largest selected dz (0.5 cm) and dtmax (0.3 d) the MIM(Se) approach did not show oscillations in h (Fig. 10a). The need for smaller time steps to eliminate oscillations makes the MIM(h) approach more time-consuming than the MIM(Se) approach. Numerical oscillations for MIM(h) were clearly caused by the mass transfer term, which for larger time steps could lead to overshooting of the equilibrium state between the two pore regions. Numerical oscillations could have been eliminated if the time stepping scheme had been adjusted to prevent this overshooting. This was done only for the MIM(Se) approach for which the adjustment was trivial due to the similarity assumption between the two pore regions, and for which accordingly no oscillations were registered.
View larger version (37K):
[in this window]
[in a new window]
|
Fig. 10. Pressure heads for the Ap5 column as simulated with MIM(Se) and MIM(h) assuming different maximum time steps (dtmax) while keeping the space increment (dz) constant at (a, b, c) 0.5 cm and (d, e, f) assuming different space increments while keeping time steps constant at 0.03 d.
|
|
|
CONCLUSIONS
|
|---|
Including experimental data of the pressure head, water content, cumulative outflow and effluent concentration into the objective functions of the MIM(Se) and MIM(h) approaches allowed us to identify model parameters, most of which appeared to be reliable as reflected by relatively small confidence intervals, generally moderate linear correlation between parameters, and only minor dependence of the final optimized parameter values on initial parameter estimates. As an exception, strong parameter correlations were found between m,s and im,s, indicating that these parameters should be determined independently. In spite of having two more adjustable parameters, MIM(h) (10 parameters) produced results that were superior to those of MIM(Se) (8 parameters) only for dry initial conditions, which indicates that at least for other than dry initial conditions, the inverse problem for MIM(h) was not well-posed. Numerical oscillations in the pressure heads for MIM(h) could be eliminated by reducing the maximum time step. The results further show that for aggregated soils, physical nonequilibrium can be more pronounced for wet and dry than for intermediate initial water contents. Even if the initial water content influences the transport process, no distinct effect was found on total nonreactive solute loss. For MIM(h), physical nonequilibrium was reflected by smaller values for the water and solute transfer coefficients and relatively large values for the related PE index. The respective parameters for MIM(Se) displayed similar tendencies but were less consistent. Overall, the MIM(Se) inverse approach was found to be well suited for simulating physical nonequilibrium solute transport during variably saturated flow. Results indicate that the MIM(h) approach can yield higher accuracy along with more consistent parameter values if the observations contain sufficient information for successful fitting. A future theoretical analysis using synthetic datasets could complement this study to further clarify such inverse optimization issues as parameter sensitivity and uniqueness for both MIM(Se) and MIM(h).
|
REFERENCES
|
|---|
- Atkins, P.W. 1990. Physical chemistry. 4th ed. W.H. Freeman, New York.
- Bahr, J.M., and J. Rubin. 1987. Direct comparison of kinetic and local equilibrium formulations for solute transport affected by surface reactions. Water Resour. Res. 23:438–452.
- Castiglione, P., B.P. Mohanty, P.J. Shouse, J. imnek, M.Th. van Genuchten, and A. Santini. 2003. Lateral water diffusion in an artificial macroporous system: Modeling and experimental evidence. Available at www.vadosezonejournal.org. Vadose Zone J. 2:212–221.[Abstract/Free Full Text]
- Chen, C., and R.J. Wagenet. 1992. Simulation of water and chemicals in macropore soils. Part I. Representation of the equivalent macropore influence and its effect on soilwater flow. J. Hydrol. (Amsterdam) 130:105–126.
- Clausnitzer, V., and J.W. Hopmans. 1995. Non-linear parameter estimation: LM
TOP
ABSTRACT
INTRODUCTION
