Numerical modeling of the effect of variation of boundary conditions on vadose zone hydraulic properties

Modelagem numérica do efeito da variação nas condições de contorno na estimativa de propriedades hidráulicas na zona vadosa

NOTA

Numerical modeling of the effect of variation of boundary conditions on vadose zone hydraulic properties

Modelagem numérica do efeito da variação nas condições de contorno na estimativa de propriedades hidráulicas na zona vadosa

Tairone Paiva Leão^I; Randall Gentry^II

^IUniversidade de Brasília – UnB. Faculdade de Agronomia e Veterinária. Campus Universitário Darcy Ribeiro - Asa Norte. Caixa Postal 4508, CEP 70910-970 Brasília (DF) Brasil. E-mail: [e-mail]

^IIUniversity of Tennessee, Department of Civil and Environmental Engineering, 223 Perkins Hall, Knoxville, TN 37996-2010, US. E-mail: [e-mail]

SUMMARY

An accurate estimation of hydraulic fluxes in the vadose zone is essential for the prediction of water, nutrient and contaminant transport in natural systems. The objective of this study was to simulate the effect of variation of boundary conditions on the estimation of hydraulic properties (i.e. water content, effective unsaturated hydraulic conductivity and hydraulic flux) in a one-dimensional unsaturated flow model domain. Unsaturated one-dimensional vertical water flow was simulated in a pure phase clay loam profile and in clay loam interlayered with silt loam distributed according to the third iteration of the Cantor Bar fractal object Simulations were performed using the numerical model Hydrus 1D. The upper and lower pressure heads were varied around average values of -55 cm for the near-saturation range. This resulted in combinations for the upper and lower constant head boundary conditions, respectively, of -50 and -60 cm, -40 and -70 cm, -30 and -80 cm, -20 and -90 cm, and -10 and -100 cm. For the drier range the average head between the upper and lower boundary conditions was set to -550 cm, resulting in the combinations -500 and -600 cm, -400 and -700 cm, -300 and -800 cm, -200 and -900 cm, and -100 and -1,000 cm, for upper and lower boundary conditions, respectively. There was an increase in water contents, fluxes and hydraulic conductivities with the increase in head difference between boundary conditions. Variation in boundary conditions in the pure phase and interlayered one-dimensional profiles caused significant deviations in fluxes, water contents and hydraulic conductivities compared to the simplest case (a head difference between the upper and lower constant head boundaries of 10 cm in the wetter range and 100 cm in the drier range).

Index terms: Numerical model, Hydrus 1D, uncertainty, boundary conditions, hydraulic conductivity, soil physics.

RESUMO

A estimativa precisa e acurada de fluxos hidráulicos na zona vadosa é um fator crítico para simular o transporte de água, nutrientes e contaminantes em sistemas naturais. O objetivo deste estudo foi simular o efeito de variações nas condições de contorno na estimativa de propriedades hidráulicas do solo (isto é, conteúdo de água, condutividade hidráulica efetiva em condições não saturadas e fluxo hidráulico) em um modelo de fluxo não saturado unidimensional. Fluxo vertical não saturado unidimensional foi simulado em um perfil de material franco-argiloso puro e em material franco-argiloso intercalado com camadas de material franco-siltoso distribuído seguindo a terceira iteração de um modelo do objeto fractal Conjunto de Cantor. As simulações foram realizadas no modelo numérico Hydrus 1D. Os potenciais totais superiores e inferiores na coluna unidimensional foram variados ao redor do valor médio de -55 cm para a faixa de valores próxima da saturação. Isso resultou em combinações para potenciais constantes nos limites superior e inferior da coluna de -50 e -60 cm; -40 e -70 cm; -30 e -80 cm; -20 e -90 cm; e -10 e -100 cm, para as condições de contorno superior e inferior, respectivamente. Para o intervalo de umidade definido como mais seco o valor médio de potencial de água no solo foi fixado em -550 cm, resultando em combinações de -500 e -600 cm; -400 e -700 cm; -300 e -800 cm; -200 e -900 cm; e -100 e -1.000 cm, para as condições de contorno superior e inferior, respectivamente. Houve aumento nos conteúdos de água, fluxos e condutividades hidráulicas com o aumento no gradiente entre os valores estipulados nos limites superior e inferior do domínio. Variação nas condições de contorno nos sistemas de fase pura e em camadas resultou em desvios significativos nos fluxos, conteúdos de água e condutividades hidráulicas com relação ao caso mais simples (uma diferença de potencial total entre os limites superior e inferior fixada em 10 cm para o intervalo mais úmido e 100 cm para o intervalo mais seco).

Termos de indexação: modelo numérico, Hydrus 1D, incerteza, condições de contorno, condutividade hidráulica, física do solo.

INTRODUCTION

Uncertainty is a factor of major influence in several areas in which scientists have to propose models and explanations of natural phenomena that cannot be comprehensively measured or characterized due to inherent spatial and/or temporal variability of processes and properties under evaluation. This is the case in unsaturated groundwater modeling, where hydraulic properties have to be estimated and characterized across areas/volumes that can range from a few centimeters to several kilometers. The modeling of natural systems is usually rather complex. The first major challenge is to develop an adequate conceptual model. If the conceptual model is inadequate, chances are that the whole modeling process will fail at some point, since the conceptual model is the foundation of the modeling process (Bredehoef, 2005). The second is to "feed" the mathematical model with representative estimations of the field parameters. Scaling and spatial variability should be considered. The measured values may represent an area of a few tens of square meters or a much larger area, depending on the characteristics of the subsurface material. Statistics and geostatistics are useful techniques in this context, since they can be used to define the minimum number of samples to adequately represent the population and spatial variability structure.

Several approaches have been used to evaluate the effect of uncertainty in groundwater models. Gaganis & Smith (2006) developed a stochastic per-datum approach to take uncertainty in groundwater predictions into account. Although not validated for real-world situations, the parameter-vector approach seems to be a better way of dealing with uncertainty than the usual single-value approach (Gaganis & Smith, 2006). Zhang et al. (2006) evaluated uncertainties in unsaturated zone flow and transport parameters at a Yucca Mountain site in variations of the fracture and matrix permeabilities and capillary strength, generating eight parameter sets, in addition to the baseline case. In general, uncertainties in matrix parameters caused larger uncertainties in simulated liquid fluxes than uncertainties in fracture properties.

The distribution of layers of materials with different physical and hence hydraulic properties is also a source of uncertainty in hydraulic models. Dillah & Protopapas (2000) evaluated uncertainty propagation in soils composed of two layers in various fine-over-coarse and coarse-over-fine configurations and using nonrandomized and randomized layering arrangements, by the Monte Carlo technique. They also evaluated the effect of varying initial and boundary conditions in the domain. The main conclusions of their study were: (a) randomization does not significantly affect the average head profiles, but does affect the variance profiles; (b) increase in uniform initial water contents causes a decrease in variance peaks, and (c) no characteristic trend for uncertainty behavior was observed in increases in supply rate (Dillah & Protopapas, 2000). Other alternative approaches to deal with uncertainty are maximum likelihood, Bayesian average (Neuman, 2003) and fuzzy set theory (Schultz & Huwe, 1997).

Hydrus is a family of numerical models used to simulate the flux of water, heat and solutes in the vadose zone. Currently three versions are available; the first, Hydrus 1D, was developed at the U.S. Salinity Lab, of the United States Department of Agriculture (Simunek et al., 2005) and may be downloaded free of charge, while 2D and 2D/3D are commercial versions (Simunek et al., 1999). Hydrus software has been extensively used for simulating fluxes in numerical domains, covering a diverse range of soils and scenarios (e.g. Simunek & van Genuchten, 2008; Simunek et al., 2009; Pontedeiro et al., 2010). Ideally, simulations should be matched and compared to field or laboratory data, but numerical models such as Hydrus can also be used for simulation only studies, as an approximate representation of real-world scenarios when such data are not available or difficult to obtain (e.g.: Pontedeiro et al., 2010). The accuracy and precision of the Hydrus simulations depend on the type of simulated problem (heat, water or solute transport), type of simulation (forward or inverse problem) and on the choice of values and discretization of statistical parameters, time parameters and initial boundary conditions, and on the choice of models available in the software for each simulation type (Simunek et al., 2005).

As mentioned above, one of the main sources of uncertainties in groundwater models is the formulation of the conceptual model. A wrong conceptualization of the problem can translate into wrong prescription of boundary conditions and parameter regions. The objective of this study was to simulate the effect of variation of boundary conditions (erroneous conceptualization of the model) on the estimation of hydraulic properties (i.e. water content, effective unsaturated hydraulic conductivity, and hydraulic flux) in a one-dimensional unsaturated flow model domain. The hypothesis of the study is that variations in boundary conditions will significantly affect flow predictions in an unsaturated system.

MATERIAL AND METHODS

Numerical simulations

Fluxes in a one-dimensional domain were simulated using Hydrus 1D software, version 3.01 (Simunek et al., 2005), from May to July 2008, at the Earth Sciences Department of the University of Tennessee, Knoxville. Two scenarios were created: in the first case, a vertical, steady-state downward water flow on a one-dimensional clay loam column (Figure 1a); in the second, a steady-state vertical flow on a clay loam column interlayered with silt loam (Figure 1b). The layers were distributed according to the 3^rd iteration of a Cantor Bar with fractal dimension D = log(2)/log(3) ≈ 0.6309 and b = 3 (Figure 1a and Leão & Perfect, 2010). Basically, in both domains, the upper and lower boundary conditions were fixed heads, while the left and right sides were no-flow boundaries (Figure 1). The flux was left to vary freely, to adjust to the fixed boundary conditions in the numerical model.

Five pressure head gradient scenarios were created for the wet and dry range simulations. For the wet range, the average head between the upper and lower boundary conditions was fixed at -55 cm. This resulted in combinations of upper and lower constant head boundary conditions of -50 and -60 cm, -40 and -70 cm, -30 and -80 cm, -20 and -90 cm, and -10 and -100 cm, respectively. For the drier range the average head between the upper and lower boundary conditions was set to -550 cm, resulting in the combinations of -500 and -600 cm, -400 and -700 cm, -300 and -800 cm, -200 and -900 cm, and -100 and -1,000 cm, for upper and lower boundary conditions, respectively. Both average pressure head scenarios were simulated in the layered and pure phase systems, resulting in 20 simulations. The ranges of pressure heads were chosen so as to generate an increase in the head gradient between the upper and lower boundary conditions around the fixed average values of -55 and -550 cm. The uncertainty was evaluated in terms of hydraulic flux deviation from the lowest gradient case: -50 to -60 cm in the wet range and -500 to -600 cm in the drier range. Deviations in hydraulic conductivity and water contents were also investigated.

The simulations were performed using the Brooks & Corey (1964) hydraulic models (Simunek et al., 2005):

and

and

where: S_e = Effective water saturation (-), α = Inverse of the air entry head value (cm^-1), h = Soil pressure head (cm), n = Pore-size distribution index (-), l = Pore-connectivity index (-), θ = Volumetric water content (cm³ cm^-3), θ_r = Residual water content (cm³ cm^-3), θ_s = Saturated water content (cm³ cm^-3), K(h) = Unsaturated hydraulic conductivity (cm h^-1), and K_s = Saturated hydraulic conductivity (cm h^-1).

The Brooks and Corey parameters for the clay loam and silt loam soils from Hydrus 1D software are presented in table 1. Hydrus 1D uses a modified form of Richard's equation to deal with the uniform (equilibrium) water movement in a partially saturated rigid porous medium (Simunek et al., 2005).

where: q is the hydraulic flux (cm h^-1), h = Pressure head (cm), θ = Volumetric water content (cm³ cm^-3), t = Time (h), x = Spatial coordinate (cm) (positive upward), S = Sink term (cm³ cm^-3 h), δ = Angle between flow direction and vertical axis (i.e., δ = 0 º for vertical flow, 90 º for horizontal flow, and 0 º < δ < 90 º for inclined flow), K(h, x) = Unsaturated hydraulic conductivity function (cm/h) given by:

where: K_r = Relative hydraulic conductivity, K_s(x) = Saturated hydraulic conductivity (cm h^-1). For downward vertical flow without water loss (e.g., by root absorption) cos δ = 0º = 1 and S is set to zero.

The numerical model domain was a one-dimensional column (height 10 cm). Hydrus 1D software was set to simulate a direct solution (forward problem) of water flow in response to a fixed hydraulic gradient, given by the upper and lower boundary conditions. The time discretization scenario was a 100 h experiment, with initial, minimum and maximum time steps of 0.01 h, 0.0001 h and 1 h, respectively. The hydraulic properties of the simulated soils were obtained from Hydrus 1D soil catalog. In Hydrus the soil profile is discretized into N-1 adjoining elements, where the ends of the elements are located at the nodal points, and N is the number of nodes. A mass-lumped linear finite element scheme was used to discretize the mixed form of Richards' equation, which is equivalent to a standard finite elements scheme (Simunek et al., 2005). Hydrus assumes the pressure head boundary conditions as first-type (Diriletch) boundary conditions. As Hydrus 1D assumes that the vertical coordinate z is positive upward, unsaturated hydraulic conductivity for unidirectional vertical hydraulic flux can be calculated as:

where K(θ) is the unsaturated hydraulic conductivity (cm h^-1).

Statistical analyses

Statistical analyses were performed using R software (Fox, 2002). The procedure lm was used to perform linear regression analysis, to compare the hydraulic properties estimated from pure phase with layered numerical simulation schemes:

where: H_P = Hydraulic property estimated from pure phase system, H_L = Hydraulic property estimated from layered system, β₁ = Slope of the regression line, β₀ = Intercept of the regression line. The validity of the linear regression analyses was evaluated using the F test from regression and the validity of the coefficients by the t test for regression (Davis, 2002).

RESULTS AND DISCUSSION

Pure phase system

The pure phase simulated columns was composed of clay loam material (Table 1). The pressure head distribution in the 1D column is shown in figure 2a for the near-saturation range (NS) and figure 2b for the drier range (DR). Although the two cases were one order of magnitude apart, the behavior of the pressure head profiles was very similar. As the difference in the upper and lower constant head boundaries increases, from 10 to 90 cm in NS and 100 to 900 cm in DR, the nonlinearity of the head distribution profiles increases (Figure 2). The volumetric water content in the columns increases with the increase in head difference (Figure 3a,b). The increase in water content with the head difference is mainly due to the increase in nonlinearity in the pressure head distribution. As the head difference increases, the head along the column tends to move asymptotically closer to the upper (higher) boundary condition, with a sharp head decrease close to the bottom of the column (Figure 2).

The Hydrus 1D forward model begins with a linear distribution of heads along the column and then solves the optimization problem using a hydraulic model (e.g. van Genuchten-Mualem, modified van Genuchten, Kosugi, Brooks & Corey) (Simunek et al., 2005). In this study the Brooks and Corey model was used for being the simplest and therefore less prone to cause numerical instabilities. The increase in gradient and water content along the column with the increase in head difference leads to an increase in hydraulic flux (Figure 3c,d). The flux increase is closer to linear in NS and more exponential-like in DR. This can be explained by the fact that the increase in the differences between constant head boundary conditions is one order of magnitude higher in DR than in NS. For the same reason, the estimated fluxes in the DR are one order of magnitude lower than in NS. One immediate conclusion from the data is that perturbations in boundary conditions in unsaturated media will be more problematic under near-saturation than under drier conditions. This is critical when modeling contaminant transport in the vadose zone, where fluxes are easily under- or overestimated, resulting in inadequate assessment of subsurface contamination risks.

Calculated unsaturated hydraulic conductivities increase with the head differences in both NS and DR (Figure 3e,f). The hydraulic conductivity values in DR are two orders of magnitude lower than in NS. This confirms that perturbations in boundary conditions can cause larger errors in hydraulic conductivities estimations at higher water contents (near-saturation) than in dry soils.

Layered system

The layered column contained clay loam interlayered with silt loam soil (Figure 1a). The horizontal layers were distributed according to the third iteration of a deterministic Cantor Bar (Leão & Perfect, 2010). The inclusion of layers of coarse grain material introduced changes in the behavior of hydraulic properties in the system. A previous study also found variations in another hydraulic parameter, soil sorptivity, as the content of coarse-grained material increased in columns simulated in 2D environment (Leão & Perfect, 2010). The pressure head profiles were not as smooth as in the case of the pure system. The instabilities in the head distribution are due to the silt loam layers with distinct hydraulic parameters (Figure 4, and also Table 1). Aside from small-scale instabilities, the overall shape of the pressure head profiles in the layered media was similar to the pure phase system. Nonlinearity increases as the pressure head difference between the upper and lower boundary conditions increases. The effective water contents in the layered media were generally higher than in the pure phase profiles (Figure 5a,b). This is because at high to intermediate pressure heads, as adopted here, the water-holding capacity of silt loam was greater than of clay loam soil.

The instabilities in the pressure head distributions and increased water contents in the layered systems affected the hydraulic flux (shown in Figure 5c,d. The silt loam layers caused an increase in the fluxes in the 1D profile. In the DR the difference was not very significant compared to the pure phase system (Figures 3d and 5d). However, in the near-saturation range, the fluxes in the layered system were about twice as high as those observed in the pure phase system (Figures 3c and 5c). In the same gradients, the inclusion of coarser-grained layers increased the fluxes. The effective hydraulic conductivities were also higher in the layered systems (Figure 5e,f). Heterogeneity plays an important role in the water and contaminant transport in unsaturated media. Incorrect estimations of physical and hydraulic properties in soil layers/horizons can significantly affect the estimations of mass transport in unsaturated soil/sediment profiles.

Comparison of results

The results of water contents, fluxes and hydraulic conductivities between the pure and layered systems were quantitatively compared (Figure 6). In general, the agreement of the estimates of the systems was very good. The coefficient of determination was 0.99 for water contents and fluxes and 0.98 for hydraulic conductivities. The intercepts (β₀) were not significantly different from zero for fluxes (Pr > |t| = 0.26077) and hydraulic conductivities (Pr > |t| = 0.4781). The slopes of the regression lines were all below one (β₁ = 0.7145 for water contents; β₁ = 0.6305 for fluxes; and β₁ = 0.6539 for hydraulic conductivity). Clay loam predicted fluxes were 0.63 times lower than those predicted for the system with silt loam layers, showing that coarse material layers can increase fluxes in the vadose zone.

Hydraulic conductivity

Nonlinear regression analysis was used to compare hydraulic conductivity between layered and pure phase systems. In both the wet and dry ranges, the relationship between hydraulic conductivity and pressure head was clearly exponential. Linearization by logarithm transformation did not make the data suitable for linear regression. Therefore, a simple exponential model was fitted to the curves using nonlinear regression procedures with the R software (Fox, 2002):

where: K(h) = Unsaturated hydraulic conductivity (cm h^-1), h = Pressure head difference between boundary conditions (cm), γ and β = Empirical coefficients.

Since the coefficients from the regression in Excel were different from those of R, only the R coefficients were reported (Table 2). Based on standard errors and estimates of the parameters g and b, the functions for the layered and pure phase systems were significantly different in the wet and not significantly different in the dry range, indicating that the presence of layers did not affect hydraulic conductivity in drier soil. However unsaturated hydraulic conductivities differed significantly between water content ranges (i.e. wet and dry), as expected, since unsaturated hydraulic conductivity is a function of the water content.

Error analysis

To calculate the error in estimated fluxes that would be caused by setting erroneous boundary conditions for a steady state flow system in porous materials we assumed that the lowest gradients (10 cm in a wetter and 100 cm in a drier range), with upper and lower pressure heads BC of -50 and -60 cm for the wet and -500 and -600 cm for the dry range, would be the true values, while the increasing simulated gradients would represent an erroneous conceptualization of the system. The error was calculated as the percent increase in hydraulic flux with the increase in pressure head gradient. For the pure phase and layered systems, both under wetter and drier conditions, the error increased exponentially with the increase in pressure head gradient (Table 3). For pure phase system the maximum error reached 1131.7 % for wet range and 4670.7 % for dry range, corresponding to gradients of 90 and 900 cm, respectively. Similar results were found for layered systems, where the maximum error reached 1406.9 % for the wet and 4852.2 % for the dry range. In other words, a relatively small error in the conceptualization of a flow system in the vadose zone can induce rather large errors in the predicted fluxes. At the end of a 24 h period the cumulative flux through the pure phase system simulated in this study would amount to 1.74 cm in the baseline case (10 cm gradient) increasing to 21.10 cm for the 90 cm gradient. This could lead to drastic changes in the calculated recharge rates for an aquifer or induce considerable errors in the estimated amount of contaminants or plant nutrients transported to groundwater and/or available to root absorption.

CONCLUSIONS

1. Variations in boundary conditions induced larger errors in calculated hydraulic conductivities at near-saturation water content values.

2. Predicted unsaturated hydraulic conductivities did not differ between the simulated interlayered and pure phase systems under drier conditions but were significantly different in the wetter range.

3. Coarser-grained silt loam layers interlayered with finer clay loam layers increased predicted fluxes in unsaturated media.

4. Errors in prescribed boundary conditions can induce serious errors in estimated fluxes in unsaturated media, which were greater under drier conditions in the layered system.

ACKNOWLEDGEMENTS

The authors wish to thank the University of Brasília, Decanato de Pesquisa e Pós-Graduação for assuming the cost of the publication fees.

LITERATURE CITED

Received for publication in November 2009 and approved in September 2010.

• BREDEHOEF, J. The conceptualization problem-surprise. Hydrogeol. J., 13:37-46, 2006.
• BROOKS, R.H. & COREY, A.T. Hydraulic properties of porous media. Fort Collins, Colorado State University, 1964. (Hydrol. Paper 3)
• DAVIS, J.C. Statistics and data analysis in geology. New York, John Wiley and Sons, 2002. 638p.
• DILLAH, D.D. & PROTOPAPAS, A.L. Uncertainty propagation in layered unsaturated soils. Transport Porous Med., 38:273-290, 2000.
• FOX, J. An R and S-Plus companion to applied regression. Sage Publications, 2002. 328p.
• GAGANIS, P. & SMITH, L. Evaluation of the uncertainty of groundwater model predictions associated with conceptual errors: A per-datum approach to model calibration. Adv. Water Res., 29:503-514, 2006.
• LEAO, T.P. & PERFECT, E. Modeling water movement in horizontal columns using fractal theory. R. Bras. Ci. Solo, 34:1463-1468, 2010.
• NEUMAN, S.P. Maximum likelihood Bayesian averaging of uncertain model predictions. Stoch. Environ. Res. Risk A., 17:291-305, 2003.
• PONTEDEIRO, E.M.; van GENUCHTEN, M.T.; COTTA, R.M. & SIMUNEK, J. The effects of preferential flow and soil texture on risk assessments of a NORM waste disposal site. J. Haz. Mat., 174:648-655, 2010.
• SCHULTZ, K. & HUWE, B. Water flow modeling in the unsaturated zone with imprecise parameters using a fuzzy approach. J. Hydrol., 201:211-229, 1997.
• SIMUNEK, J. & HOPMANS, J.W. Parameter optimization and nonlinear fitting. In: DANE, J.H. & TOPP, G.C., eds. Methods of soil analysis: Physical methods. 3.ed. Madison, American Society of Agronomy, 2002. p.139-157.
• SIMUNEK, J. & van GENUCHTEN, M.T. Modeling Nonequilibrium Flow and Transport Processes Using HYDRUS. Vadoze Zone J., 7:782-797, 2008.
• SIMUNEK, J.; SEJNA, M. & van GENUCHTEN, M.T. The HYDRUS-2D software package for simulating two-dimensional movement of water, heat, and multiple solutes in variably saturated media. Version 2.0, IGWMC - TPS - 53. Golden, International Ground Water Modeling Center, Colorado School of Mines, 1999. 251p.
• SIMUNEK, J.; van GENUCHTEN, M.T. & SEJNA, M. The HYDRUS-1D Software package for simulating the movement of water, heat, and multiple solutes in variably saturated media, Version 3.0. Riverside, Department of Environmental Sciences, University of California Riverside, 2005. 270p. (HYDRUS Software Series 1)
• SIMUNEK, J.; JACQUES, D.; TAWARAKAVI, N.K.C. & van GENUCHTEN, M.T. Selected HYDRUS modules for modeling subsurface flow and contaminant transport as influenced by biological processes at various scales. Biologia, 64:465-469, 2009.
• ZHANG, K.; WU, Y-S. & HOUSEWORTH, J.E. Sensitivity analysis of hydrogeological parameters in modeling flow and transport in the unsaturated zone of Yucca Mountain, Nevada, USA. Hydrogeol. J., 14:1599-1619, 2005.

Publication Dates

• Publication in this collection
  25 Apr 2011
• Date of issue
  Feb 2011

History

• Received
  Nov 2009
• Accepted
  Sept 2010