Access the full text.

Sign up today, get DeepDyve free for 14 days.

Journal of Hydrology and Hydromechanics
, Volume 60 (3) – Sep 1, 2012

/lp/de-gruyter/uncertainty-analysis-of-a-dual-continuum-model-used-to-simulate-0Fpazkw1Rn

- Publisher
- de Gruyter
- Copyright
- Copyright © 2012 by the
- ISSN
- 0042-790X
- DOI
- 10.2478/v10098-012-0017-0
- Publisher site
- See Article on Publisher Site

J. Hydrol. Hydromech., 60, 2012, 3, 194205 DOI: 10.2478/v10098-012-0017-0 MICHAL DOHNAL, TOMÁS VOGEL, MARTIN SANDA, VLADIMÍRA JELÍNKOVÁ Czech Technical University in Prague, Faculty of Civil Engineering, Thákurova 7, 166 29 Prague, Czech Republic; Mailto: dohnalm@mat.fsv.cvut.cz A one-dimensional dual-continuum model (also known as dual-permeability model) was used to simulate the lateral component of subsurface runoff and variations in the natural 18O content in hillslope discharge. Model predictions were analyzed using the GLUE generalized likelihood uncertainty estimation procedure. Model sensitivity was evaluated by varying two separate triplets of parameters. The first triplet consisted of key parameters determining the preferential flow regime, i.e., the volumetric proportion of the preferential flow domain, a first-order transfer coefficient characterizing soil water exchange between the two flow domains of the dual-continuum system, and the saturated hydraulic conductivity of the preferential flow domain. The second triplet involved parameters controlling exclusively the soil hydraulic properties of the preferential flow domain, i.e., its retention curve and hydraulic conductivity function. Results of the analysis suggest high sensitivity to all parameters of the first triplet, and large differences in sensitivity to the parameters of the second triplet. The sensitivity analysis also confirmed a significant improvement in the identifiability of preferential flow parameters when 18O content was added to the objective function. KEY WORDS: Hillslope Discharge, Preferential Flow, Dual-Permeability Model, Sensitivity Analysis, GLUE, 18O content. Michal Dohnal, Tomás Vogel, Martin Sanda, Vladimíra Jelínková: ANALÝZA NEJISTOT PI MODELOVÁNÍ PODPOVRCHOVÉHO ODTOKU ZE SVAHU METODOU DUÁLNÍHO KONTINUA S VYUZITÍM IZOTOPU KYSLÍKU 18O JAKO PIROZENÉHO STOPOVACE. J. Hydrol. Hydromech., 60, 2012, 3; 27 lit., 10 obr., 2 tab. K simulacím laterální slozky podpovrchového proudní a zmn koncentrace izotopu kyslíku 18O ve vod vytékající ze svahu byl pouzit jednorozmrný model vyuzívající pístupu duálního kontinua. Nejistota modelových pedpovdí byla odhadnuta s vyuzitím metody zobecnné vrohodnosti (GLUE). Citlivost modelu byla zjisována pomocí variací dvou samostatných trojic parametr. První trojice sestávala z klícových parametr pro urcení rezimu preferencního proudní, tj. objemového podílu preferencní domény proudní, penosového koeficientu charakterizujícího výmnu vody mezi obma doménami duálního systému a nasycené hydraulické vodivosti preferencní domény. Druhá trojice zahrnovala výhradn parametry urcující hydraulické charakteristiky preferencní domény proudní, tj. retencní kivku a funkci hydraulické vodivosti. Z výsledk analýzy vyplývá vysoká citlivost modelu na vsechny parametry z první trojice a velké rozdíly v citlivostech parametr druhé trojice. Analýza dále potvrdila významné zlepsení zjistitelnosti parametr preferencní domény v pípad, kdy je do cílové funkce pidána koncentrace izotopu kyslíku 18O. KLÍCOVÁ SLOVA: odtok ze svahu, preferencní proudní, model duální permeability, analýza citlivosti, GLUE, koncentrace 18O. Introduction Predictions obtained with hydrological models are always subject to considerable uncertainty. Several sources of uncertainty can be distinguished: inherent model uncertainty, uncertainty in the model parameters and uncertainty associated with the formulation of initial and boundary conditions. The 194 estimation of prediction uncertainty is an important topic in the past and recent hydrologic literature (e.g., Beven and Binley, 1992; Abbaspour and Yang, 2006). From the point of view of parameter estimation and model calibration, uncertainties may lead to a plateau in the response surface of the objective function or the presence of several distinct optimal parameter sets in different regions of parameter space. The situation when many parameter sets produce the same goodness-of-fit value, and thus lead to similarly acceptable representations of the system, has been called "equifinality" by Beven and Binley, (1992). To enable an assessment of prediction uncertainty, they proposed a methodology based on generalized likelihood uncertainty estimation (GLUE). The GLUE methodology has been used in various hydrological studies, including those dealing with the TOPMODEL (e.g. Blazková et al., 2002; Lamb et al., 1998), the MIKE-SHE model (Christiaens and Feyen, 2002), and the HYDRUS-1D model (Hansson and Lundin, 2006). To our knowledge no comparable numerical study has been conducted to analyze prediction uncertainties associated with the parameterization of preferential flow in dual-continuum modeling of subsurface runoff. Dual-continuum models (often also termed dualpermeability models) of soil water flow and solute transport account for preferential flow effects by dividing the flow domain into a preferential flow domain and a soil matrix domain (Gerke and van Genuchten, 1993). These models have been successfully used in a wide range of applications, including runoff generation (Pavelková et al., 2012), cadmium and chlorotoluron transport in macroporous soils (Dusek et al., 2006; Kodesová et al., 2005) and bromide transport in tile-drained fields (Gerke et al., 2007). More recently, Vogel et al. (2010b) successfully used a dual-continuum model to explain observed variations of 18O in hillslope discharge. Simnek et al. (2003) noted that very few studies have looked at the possibility of optimizing preferential flow models. The low sensitivity of parameters associated with the retention capacity of the preferential flow domain was reported by Dohnal (2008) to cause inverse modeling problems. Laloy et al. (2010) experienced convergence problems with the HYDRUS-1D code, including complex parameter correlations, for larger values of soil water transfer coefficient in dual-permeability simulations during inversion. On the other hand, Dolezal et al. (2007) obtained reasonably good agreement between observed and simulated soil water pressures based on inversely estimated soil hydraulic parameters of a dual-continuum model. Our present study complements a previous study of Vogel et al. (2010b) in which 18O was used as a natural tracer to study subsurface hillslope runoff. The present study adopted the same modeling tools, except that they are placed within the GLUE framework in which the parameter values from Vogel et al. (2010b) are treated as reference values, and whereby multiple model responses are calculated for randomly generated parameter combinations. Specific objectives of the present study are: (i) to analyze the performance of a dual-continuum model using the GLUE methodology, (ii) to elucidate the model sensitivity to a variety of combinations of preferential flow parameter values, including those representing weak preferential flow conditions associated with low interfacial resistances, and (iii) to evaluate the uncertainty of predicted discharges and isotopic signatures. Material and methods Uhlíská catchment Uhlíská is a small (1.78 km2) experimental catchment with long-term monitoring of hydrological and climate variables, situated in the western part of Jizera Mts., Czech Republic. Continuous hydrological and micro-meteorological observations are performed at the Tomsovka experimental hillslope site. Measurements are situated in young forest gradually succeeding grass dominated vegetation. Uhlíská was subject to large-scale deforestation caused by acid rains at the end of 1980s. The hillslope is currently covered with patches of grass (Calamagrostis villosa) interspersed with spruce trees (Picea abies). The shallow sandy loam soil at Tomsovka is classified as Dystric Cambisol. Runoff formation at the hillslope site is affected by the presence of preferential pathways in the shallow layered soil profile. Preferential flow effects have been attributed to highly conductive pathways along decayed tree roots, as well as to soil structure. The micro-structural component of preferential flow was studied under similar pedological conditions by Císlerová and Votrubová (2002) and Snhota et al. (2010). During major rainfallrunoff events at Tomsovka, a considerable amount of water moves laterally in the subsurface along the shallow soilbedrock interface. The lateral component of subsurface runoff was measured in an experimental trench. The trench consisted of two sections: A and B (see Fig. 1). The contributing area of each section was estimated to be about 100 m2 (for details see Sanda and Císlerová, 2009). 195 Soil water as well as effluent water collected at the trench during major rainfall-runoff events were sampled for natural 18O concentrations. Soil water was sampled at approximately monthly intervals by means of suction cups installed at depths of 30 and 60 cm below the soil surface. Hillslope discharge (A and B combined) was sampled at 6-h intervals. Soil water flow model The modeling approach is based on the assumption that soil water flow takes place in a system of two parallel flow domains: the soil matrix domain (SM domain), and the preferential flow domain (PF domain). Vertical movement of water is described by a dual set of one-dimensional Richards equations (e.g., Gerke and van Genuchten, 1993). These equations are coupled through a first-order soil water transfer term, which allows for the dynamic exchange of water between the PF and SM domains. The dual system of governing equation can be written as w f f t h f + 1 - w f S f - ws K ar (h f - hm ) , wf K f z z (1) h wm m = wmK m m + 1 - wmS m + ws K ar (h f - hm ) , t z z where the subscript m denotes the SM domain, the subscript f denotes the PF domain, h the soil water pressure head [m], K the unsaturated hydraulic conductivity [m s-1], the volumetric soil water content [m3 m-3], S the rate of the local root water uptake [s-1], ws the first-order interdomain soil water transfer coefficient at saturation [m-1 s-1], Kar the relative unsaturated hydraulic conductivity of the SMPF interface, and wm and wf the volumetric fractions of soil occupied by the respective flow domains, such that wm + wf = 1. (2) In this expression, is the specific interfacial area [m-1], a characteristic length [m] associated with the soil water transfer between the PF and SM domains, Ksa the saturated interfacial hydraulic conductivity [m s-1] and Ksm the saturated conductivity of the soil matrix [m s-1]. The transport of natural tracer 18O in the dualcontinuum system is described by a similar dual set of advectiondispersion equations (Vogel et al., 2010b). These equations are not restated here. The governing equations for soil water flow and 18O transport are solved using the numerical model S1D developed at the Czech Technical University in Prague. The current version of the S1D code is described in Vogel et al. (2010a). Scaling of soil hydraulic characteristics The S1D model employs a scaling procedure designed to simplify the description of spatial variability in the soil hydraulic properties of heterogeneous soil profiles. The scaling procedure consists of a set of linear scaling transformations which relate the individual soil hydraulic characteristics (h) and K(h) to reference characteristics *(h*) and K*(h*) (Vogel et al., 1991). The scaling procedure, as applied in the present study, makes use of three spatially invariant scaling factors representing the PFdomain properties. Their use leads to the following scaling transformations: Fig. 1. Schematic of the experimental trench for collecting subsurface hillslope discharge. Details about how governing equations are coupled physically and numerically can be found in Vogel et al. (2010a), where the transfer coefficient ws is further parameterized as ws = K sa K sa K sm . (3) f (h f ) = + f [ * (h* ) - ] , fr f f fr (4) h f = hf h* , f K f (h f ) = Kf K * (h* ) , f f (5) (6) where f, hf and Kf are the PF-domain scaling factors for soil water content, pressure head and hydraulic conductivity, respectively. GLUE methodology The GLUE methodology (Beven and Binley, 1992) utilizes Monte Carlo simulations to determine a quantitative measure of model performance (the likelihood measure) for each combination of model parameters drawn from a selected parametric space. The resulting model performance distributions are then used to assess the model sensitivity to individual parameters and to estimate the model prediction uncertainty. The simulations with high likelihood measure, which simulate the behavior of a real system reasonably well, were called behavioral by Beven and Binley (1992). Less successful simulations were called non-behavioral. Model performance measure: In this study, the Nash-Sutcliffe efficiency (NSE; Nash and Sutcliffe, 1970) was used to evaluate model performance: ing the effect of inclusion of different types of observations in the objective function. In our case, uniform prior distributions were assumed for all parameters in Monte Carlo sampling. The posterior parameter distributions were constructed as normalized cumulative efficiency distributions (NCED). NCEDs were determined for each parameter by taking into account the ratio between the number of successful model runs associated with a particular parameter value and the number of successful runs executed for the entire parametric space. A model run was rated as successful (for the purpose of calculating NCED) when the corresponding model efficiency value counted among the best 1000 values in a given parametric space. Prediction uncertainty: To estimate prediction uncertainty, first the model responses are labeled by a performance measure and ranked to determine a cumulative probability distribution. Then, uncertainty quantiles - representing prediction limits - are derived from the distribution for a selected level of significance. Prediction limits are estimated separately for each simulation time and each type of observation (Beven and Binley, 1992). Modeling hillslope runoff responses at Tomsovka The S1D software was used to explain the observed soil water dynamics, subsurface hillslope discharge rate and variations in natural isotope 18O in resident and effluent soil water. The numerical simulations were performed for two consecutive growing seasons (2007 and 2008). The soil water pressure profile measured by tensiometers was used as initial condition at the beginning of each simulated period. S1D further accounted for soil atmosphere interactions involving natural rainfall and plant transpiration. Direct evaporation from the soil surface was considered negligible due to the presence of dense grass cover and natural surface mulch. The lower boundary condition at a depth of 75 cm below the soil surface was formulated as a unit hydraulic gradient condition for both domains of the dual-continuum system. Episodic outflow generated at the lower boundary of the PF domain was assumed to feed the saturated subsurface flow process during major rainfallrunoff events. To obtain the subsurface discharge rates to the trench, the vertical PF-domain fluxes were simply multiplied by the upslope contributing area. The soil water flux computed at the lower boundary of the SM domain was assumed to percolate vertically to 197 NSE = 1- n i=1(xs,i n i=1(xs,i - xo,i )2 - xo ) 2 (7) where n is the number of observation times, xo,i the observed value at time ti, xs,i the corresponding simulated value and xo is the mean of xo,i. NSE values may range from - to 1. As shown by Gupta et al. (2009) model predictions with NSE values less than zero are generally not better than simply using the observed mean as a predictor. Scatter diagrams: The relationship between model efficiency and position in parametric space can be visualized by means of scattergrams (Beven and Binley, 1992). Each dot in a scattergram represents a single model run (one process simulation). The dot scatter corresponding to a particular value of a parameter reflects the effect of variations in the remaining parameters on model efficiency. Sensitivity analysis: Useful information about model sensitivity can be obtained by comparing prior and posterior parameter distributions (e.g. Hansson and Lundin, 2006). Significant difference between the two distributions for a parameter indicates high model sensitivity to that parameter. The posterior distributions are also useful when study- deeper horizons. The two upslope microcatchments, connecting to trenches A and B were modeled as identical, i.e., no differences in geometric and material properties, and initial and boundary conditions were considered. Daily potential transpiration was calculated using the PenmanMonteith formula (Monteith, 1981) using micrometeorological data observed directly at the Tomsovka site. Root water uptake was modeled using the approach of Feddes et al. (1978). The vertical distribution of the local uptake rate was assumed to be constant to the depth of 20 cm and then to decrease linearly to the depth of 70 cm. The soil profile was divided into four layers. The unsaturated soil hydraulic properties of each layer were characterized using modified van Genuchten- Mualem expressions as discussed by Vogel et al. (2000). The reference values of soil hydraulic parameters were adopted from Vogel et al (2010b), Tab. 1. The volumetric fraction of the PF domain, wf, was set equal to 0.07 at the soil surface and to 0.05 at the lower boundary of the flow domain, varying linearly between the two points. The value of ws was also estimated to decrease linearly between the soil surface and the lower boundary, from 1 to 0.01 cm-1 d-1, consistently with the decreasing hydraulic conductivity of the soil matrix. More details about the applied modeling approach and model assumptions can be found in previous studies of Vogel et al. (2010b) and Dusek et al. (2012). T a b l e 1. Reference soil hydraulic parameters (adopted from Vogel et al., 2010b). Domain SM Layer 1 2 3 4 Depth [cm] 0-8 8-20 20-70 70-75 0-75 r [-] 0.20 0.20 0.20 0.20 0.01 s [-] 0.55 0.54 0.49 0.41 0.60 VG [cm-1] 0.050 0.050 0.020 0.020 0.050 nVG [-] 2.00 1.50 1.20 1.20 3.00 Ks [cm d-1] 567 67 17 1.3 5000 hs [cm] 0.00 -0.69 -1.48 -1.88 0.00 PF r and s are the residual and saturated soil water contents, respectively, Ks the saturated hydraulic conductivity, hs the air-entry value of Vogel et al. (2000), and VG and nVG are empirical parameters Application of GLUE procedure In the following section we describe procedures that were used to select model parameters and their variation ranges for the model sensitivity/uncertainty analysis. Furthermore, we describe the composition of the NSE objective function, used to evaluate the model performance. Monte Carlo sampling of preferential flow parameters The number of parameters necessary for a full characterization of the dual-continuum system is quite large. Some reduction of this number, based on reasonable simplifying assumptions, is needed to make the sensitivity/uncertainty analysis effective. The choice of parameters and determination of their feasible variation ranges are therefore crucial steps in a successful GLUE application. Two independent triplets of parameters were selected to assess the sensitivity of the simulated hillslope responses to preferential flow parameters. The first triplet included key parameters controlling the volume and intensity of preferential flow in the dual-continuum 198 model, i.e. the interdomain soil water transfer coefficient, ws, the volumetric fraction of the PF domain, wf, and the PF-domain saturated hydraulic conductivity, Ksf. The second triplet included the PF-domain scaling factors: f, hf and Kf. The Monte Carlo sampling procedure was used to generate a large number of parameter combinations (in the order of 104) from the selected ranges of parameter triplets (Tab. 2), with each parameter value being drawn randomly from a uniform distribution. For the purpose of Monte Carlo sampling, the parameters wf and ws were assumed to vary according to the following expressions (as illustrated in Fig. 2): ws (z) = C * (z) ws w f (z) = C w w* (z), f (8) where C and Cw are depth invariant multiplicative coefficients [-], representing the parameters wf and ws in random sampling. The sampled ranges of these coefficients are given in Tab. 2. The range of ws was determined so as to examine a broad spectrum of flow situations involving different degrees of preferential flow and inter- domain communication. This was achieved by transforming the estimated ranges of the parameters , and Ksa using Eq. (3), as described below. Unlike the range of ws, the feasible ranges of these parameters could be more easily determined based on their physical interpretation. The range of the characteristic length was set according to our hypothesis that preferential flow at Tomsovka is mainly caused by the presence of decayed roots, which are about 1 to 100 cm apart. Based on the same hypothesis, the specific interfacial area was identified with the specific root surface area, which was estimated from the root length density and rhizosphere radius (de Jong van Lier, 2008). Hence, and were assigned variation ranges of (1,100) cm and (0.001, 1) cm-1. The saturated interfacial conductivity, Ksa, was assumed to be equal to the saturated conductivity of the soil matrix Ksm, which is a relatively easily measurable quantity. It is treated here as invariant within each soil layer and assigned the reference values in Tab. 1. Using Eq. (3), the sampling range of ws was determined from the estimated ranges of , and Ksa as ws (Ksm×10-5, Ksm) cm-1 d-1, as illustrated in Fig. 3. The figure shows how the estimated ranges of and are projected on the ws variation range. The parameter ws also depends on the position below the soil surface (via Ksm). To avoid the use of depth dependent parameters, the soil water transfer coefficient, ws, was represented in the Monte Carlo sampling scheme by the spatially invariant parameter C (Eq. (8), Tab. 2 and Fig. 2). The variation range of the volumetric fraction of the PF domain, wf, was estimated as wf (0.01, 0.21) at the soil surface and wf (0.01, 0.15) at the bottom of the soil profile. The volumetric fraction of the PF domain, wf, was (similarly to ws) represented in Monte Carlo sampling by the spatially invariant parameter Cw (Eq. (8), Tab. 2 and Fig. 2). The saturated hydraulic conductivity of the PF domain, Ksf, was assumed to be constant with depth. The range used in Monte Carlo sampling of Ksf was chosen so as to retain the character of preferential flow near the lower bound and to preserve a laminar flow regime near the upper bound of the Ksf range. For the second triplet, the PF-domain scaling factors of soil water content and pressure head, f and hf, were allowed to vary within one order of magnitude about their reference value of 1. The sampling interval of the scaling factor Kf was set consistently with the range of the parameter Ksf in the first triplet. The effect of varying scaling factors f and hf on the shape of the PF-domain retention curve is shown in Fig. 4. T a b l e 2. S1D parameter ranges used in the Monte Carlo sampling scheme. Sampled triplet 1st 2nd Parameter C [-] Cw [-] Ksf [cm d-1] f [-] hf [-] Kf [-] Lower bound 0.006 0.15 170 0.2 0.1 0.015 Upper bound 65 3.0 17 000 1.4 10 2.0 ws (cm-1 d-1) 0.001 0.1 Cw = 0.15 wf ( - ) 10 15 Depth (cm) C = 6 ×10-3 C = 1.0 C = 65 Cw = 1.0 Cw = 3.0 Fig. 2. Prescribed depth-dependent variations of the interdomain soil water transfer coefficient, ws, and the PF-domain volumetric fraction, wf, in the GLUE procedure. The solid lines represent the reference case for which the multipliers C and Cw are equal to 1 (see Eq. (8)). Dashed lines correspond to the lower and upper bounds of the respective multipliers. = 10 = 0.1 10 (cm) = 1 = 0.01 Fig. 3. Sampling range of ws [cm-1 d-1] derived using Eq. (3) based on variations of and . - 100000 Reference hf = 0.1 pared with the observed data. Four different types of observations were used to calculate NSE: P1. Resident concentrations of 18O in soil water; P2. Flux concentrations of 18O in effluent water (A + B combined); P3. Subsurface hillslope discharge measured in trench section A; P4. Subsurface hillslope discharge measured in trench section B. Two alternative NSE objective functions, corresponding to the two trench sections (A and B), were defined. One objective function was associated with trench section A and consisted of the weighted contributions of P1, P2 and P3 observations. The second objective function was associated with trench section B and consisted of the weighted contributions of P1, P2 and P4 observations. Results and discussion (cm-1) Pressure head (cm) - 10000 - 1000 hf = 10 f = 0.2 f = 1.4 - 100 - 10 -1 0.0 0.2 0.4 0.6 Water content ( - ) 0.8 Fig. 4. Prescribed variation of the PF-domain retention curve in the GLUE procedure. The reference retention curve is defined by f = 1 and hf = 1. According to Vogel et al. (2010a), the flow regime in a dual-continuum system can be characterized by four dimensionless parameters: , wf, Ksf/Ksm and Ksa/Ksm. These parameters represent an attractive parameter set for Monte Carlo sampling. Note that using the assumption of Ksa/Ksm = 1, which seems to be a valid simplification for our experimental site, this combination of parameters reduces to a set which is equivalent to our first triplet (ws, wf and Ksf). Composition of objective function The subsurface flow and transport simulations were performed for each parameter set and com- As a result of the GLUE analysis, we obtained 13 111 behavioral simulations (NSE > 0) for the first parameter triplet and 1 107 for the second triplet. 6 889 simulations for the 1st triplet and 18 892 for the 2nd triplet were marked as non-behavioral (NSE 0). As shown in the scattergrams (Fig. 5), the parameters C , Cw and Ksf are quite well identified (i.e., there is a well-defined NSE maximum in each of the three respective scattergrams). The agreement between model responses and observations is better for lower C , and higher Ksf values. For Cw, the best model efficiency was obtained in close proximity to the reference value. As far as the 2nd parameter triplet is concerned, the scaling factor of the pressure head, hf, has a very well-defined maximum, while the other two scaling factors, f and Kf, have poorly resolved maxima. Although Ksf affects model efficiency quite significantly for the 1st triplet, the functionally equivalent parameter Kf (see Eq. (6)) seems to play an unimportant role for the 2nd triplet. This is caused by a strong physical link between f and Kf, which leads to substantial correlation between the two parameters. The significantly lower number of behavioral simulations obtained for the 2nd triplet was caused by the exceedingly wide range of hf. This could have been amended by narrowing the sampling range of hf and repeating the whole Monte Carlo procedure for the 2nd triplet. However, since the number of behavioral simulations (1107) obtained with the present hf range was still reasonably large, the repetition of the Monte Carlo procedure was deemed unnecessary. Of the 1 107 behavioral simulations, 90% fell inside the interval 0.8 < hf < 1.2 (which corresponds to 0.042 < VG < 0.063). This is in a good agreement with our perception that the preferential pathways at Tomsovka are mainly biopores characterized by relatively flat retention curves (see Fig. 4). A simple comparison of model efficiencies obtained for trench sections A and B is shown in Fig. 6. The figure depicts the shape of the upper envelope of the respective scattergrams. The envelope curves were constructed alternatively with the P3 and P4 observations in the objective function, i.e., by taking into account the discharge observations from trench sections A or B. The envelope curves are quite similar for both sections. However, different NSE maxima for parameters Cw and Ksf can be distinguished. In general, the measured subsurface discharge of section B can be explained better by the model than the discharge of section A. We further note that the model run based on the reference parameters (Tab. 1) belongs to the best 2% behavioral simulations (in terms of model effi- ciency) for the 1st triplet and to the best 10% for the 2nd triplet of parameters. Model sensitivity The posterior distributions of model parameters are shown in Fig. 7, where the individual distributions were conditioned on different types of observations (P1, P2 and P4). The distributions were constructed from the best 1000 model runs. Most posterior distributions differ significantly from the uniform prior distributions, except for Kf and partly also for f. The distributions, which were conditioned on the subsurface hillslope discharge (P4) only, indicate that the model is most sensitive to the fraction of the PF domain wf (represented by the parameter Cw) and to the parameter hf. The posterior distributions conditioned on the 18 O signature in the subsurface discharge (P2) suggest large sensitivity to the soil water transfer coefficient ws (represented by C) and fair sensitivity to all other parameters except Kf. Also, the parameter f was found to influence the 18O content in hillslope discharge, while having little effect on the magnitude of discharge and 18O content in soil water sampled by suction cups. 20 40 C ( - ) 2 Cw ( - ) 0.8 0.6 0.4 0.2 0 0.2 0.8 f ( - ) 1.4 5 hf ( - ) Fig. 5. Scattergrams for the selected model parameters (1st row = 1st triplet, 2nd row = 2nd triplet). Only scattergrams related to trench section B (conditioned on P1, P2 and P4 observations) are shown here. 20 40 C ( - ) 2 Cw ( - ) 0.9 0.2 0.8 f ( - ) 1.4 0 5 hf ( - ) 10 Section A Section B Fig. 6. Upper envelopes of the scattergrams for trench sections A and B. For this comparison, the objective function was composed of a single observation type: P3 for section A, and P4 for section B. NCED ( - ) 0 0 1 20 40 C ( - ) 60 1 2 Cw ( - ) NCED ( - ) 0 0.2 0.8 f ( - ) 1.4 hf ( - ) Resident (P1) Effluent (P2) Discharge (P4) Prior Fig. 7. Normalized cumulative posterior distributions of parameters used in the sensitivity analysis. The distributions were conditioned on three different types of observations relevant to trench section B: P1 - 18O content of soil water; P2 - 18O content in hillslope discharge (A+B); P4 - discharge measured in trench section B. Cumulative prior distributions of parameters appear as diagonal lines. Prediction uncertainty Fig. 8 shows prediction limits for the simulated subsurface hillslope discharge observed in trench section B. The limits were determined as 5% and 95% prediction quantiles from the best 1000 behavioral parameter sets. This implies that they do not follow any particular simulation. The observed cumulative discharge somewhat departs from the predicted limits in the middle of the season, but then returns to the predicted range in later stages of the season. In Fig. 9, the suction cup 18O data are compared with the prediction median calculated from the simulated SM-domain resident 18O concentrations. The prediction median, instead of uncertainty limits, is shown because of very small differences between the estimated prediction quantiles. The narrow prediction limits are caused by the fact that the resident 18O concentrations, determined by suction cups, represent water contained in the soil matrix rather than in preferential flow pathways. Thus, the variations of the PF-domain parameters do not affect the simulation results significantly. The good agreement between the observed and predicted values is due to a reasonable choice of SM-domain reference parameters (Tab.1). measured 18O contents were inside the prediction limits. Some of the short-term fluctuations, however, fell outside. -6 Observation Prediction median -8 18O () -10 -12 16-Jun-07 10-Nov-07 5-Apr-08 30-Aug-08 Fig. 9. Prediction median for the simulated resident concentration of 18O in the SM-domain (conditioned on P1, P2 and P4 observations). The observed concentrations were measured in soil water extracted by suction cup at depth of 60 cm. -4 Observation -6 Cumulative specific discharge (cm) Observation 18O () -8 -10 -12 15-Jun-07 29-Jul-07 11-Sep-07 25-Oct-07 0 14-May-08 7-Jul-08 30-Aug-08 23-Oct-08 Fig. 10. Prediction limits for 18O concentrations in subsurface discharge (conditioned on P1, P2 and P4 observations). The observed concentrations were measured in effluent water collected from both trench sections (A and B sections combined). Fig. 8. Prediction limits for the simulated subsurface hillslope discharge (conditioned on P1, P2 and P4 observations). The observed discharges were measured in trench section B. Summary and conclusions Uncertainty of dual-continuum modeling of subsurface flow and 18O transport was studied by means of Monte Carlo simulations. Key model parameters were sampled randomly from their estimated prior distributions. Due to the fully coupled numerical solution scheme of the governing equations used in the S1D software code, the conver203 The isotope signatures of the subsurface discharge, observed at trench sections A and B combined, are compared with the corresponding model responses in Fig. 10. The prediction limits were calculated from simulated flux concentrations at the lower boundary of the PF domain. The majority of gence problems reported by Laloy et al. (2010) were significantly diminished. Consequently, substantially broader parameter ranges could be sampled. The GLUE methodology was used to analyze the sensitivity of the model to variations of selected parameters and to determine prediction uncertainty for the simulated hillslope discharge and 18O concentration (both resident and effluent). The sensitivity analysis was performed for two separate triplets of parameters. The first triplet consisted of key parameters determining the preferential flow regime, i.e., the proportion of the preferential flow domain, the resistance of the interface separating the PF domain from the soil matrix domain and the hydraulic conductivity of the PF domain. The second triplet involved parameters controlling the soil hydraulic properties of the PF domain, i.e. its retention curve and hydraulic conductivity function. Results of the sensitivity analysis suggest high model sensitivity to all parameters of the 1st triplet. Among the parameters of the 2nd triplet, the scaling factor for pressure head, hf, showed high sensitivity over a very narrow posterior feasibility range, which indicates that the analysis could have been simplified by fixing hf at its reference value. The scaling factor for the soil water content, f, was found to be insensitive from the point of view of hillslope discharge responses, but moderately sensitive in terms of 18O fluctuations in effluent water. The inclusion of 18O data in the analysis significantly improved the identifiability of the parameter C, which controls the transfer between the preferential and soil matrix domains. Application of the GLUE procedure confirmed the importance of considering the preferential nature of subsurface flow when modeling runoff generation at the hillslope scale. Acknowledgments. The research was supported by the Czech Science Foundation-projects 205/09/ 0831 and 205/08/1174.

Journal of Hydrology and Hydromechanics – de Gruyter

**Published: ** Sep 1, 2012

Loading...

You can share this free article with as many people as you like with the url below! We hope you enjoy this feature!

Read and print from thousands of top scholarly journals.

System error. Please try again!

Already have an account? Log in

Bookmark this article. You can see your Bookmarks on your DeepDyve Library.

To save an article, **log in** first, or **sign up** for a DeepDyve account if you don’t already have one.

Copy and paste the desired citation format or use the link below to download a file formatted for EndNote

Access the full text.

Sign up today, get DeepDyve free for 14 days.

All DeepDyve websites use cookies to improve your online experience. They were placed on your computer when you launched this website. You can change your cookie settings through your browser.