Simple Hydraulic Conductivity Estimation by the Kalman Filtered Double Constraint Method by M.A. El-Rawy1,2 , O. Batelaan1,3 , and W. Zijl1
Abstract This paper presents the Kalman Filtered Double Constraint Method (DCM-KF) as a technique to estimate the hydraulic conductivities in the grid blocks of a groundwater flow model. The DCM is based on two forward runs with the same initial grid block conductivities, but with alternating flux-head conditions specified on parts of the boundary and the wells. These two runs are defined as: (1) the flux run, with specified fluxes (recharge and well abstractions), and (2) the head run, with specified heads (measured in piezometers). Conductivities are then estimated as the initial conductivities multiplied by the fluxes obtained from the flux run and divided by the fluxes obtained from the head run. The DCM is easy to implement in combination with existing models (e.g., MODFLOW). Sufficiently accurate conductivities are obtained after a few iterations. Because of errors in the specified head-flux couples, repeated estimation under varying hydrological conditions results in different conductivities. A time-independent estimate of the conductivities and their inaccuracy can be obtained by a simple linear KF with modest computational requirements. For the Kleine Nete catchment, Belgium, the DCM-KF yields sufficiently accurate calibrated conductivities. The method also results in distinguishing regions where the head-flux observations influence the calibration from areas where it is not able to influence the hydraulic conductivity.
Introduction Do you calibrate the hydraulic conductivity in your groundwater flow model based on head observations, while you specify the recharge? Then you can estimate the hydraulic conductivity in the neighbourhood of the observations by performing alternating model runs specifying head and recharge. In this paper we propose, demonstrate and discuss the Kalman Filtered Double Constraint Method (DCM-KF). The combined approach of the DCM with the KF makes a simple but effective form of direct inverse calibration possible. To model a groundwater flow system, two types of problems must be solved: (1) the forward problem (simulation) and (2) the inverse problem (calibration). Inverse problems arise in various fields such as medicine, meteorology and hydrology whenever mathematical models are used to determine unknown properties of a system from the observations of a response of that system. Usually the 1 Department of Hydrology and Hydraulic Engineering, Vrije Universiteit Brussel, Pleinlaan 2, 1050 Brussels, Belgium. 2 Corresponding author: Department of Civil Engineering, Faculty of Engineering, Minia University, Minia 61111, Egypt; fax: +2-086-2346674; [email protected] 3 National Centre for Groundwater Research and Training, School of the Environment, Flinders University, GPO Box 2100, Adelaide, SA 5001, Australia. Received August 2013, accepted March 2014. © 2014, National Ground Water Association. doi: 10.1111/gwat.12217
NGWA.org
unknown properties are represented as physical parameters characterizing the model and they are systematically adjusted in order to match the model outputs with the observed values. Inverse problems in groundwater modeling have been studied widely during the last four decades (de Marsily et al. 2000) and have been reviewed in detail by Yeh (1986), Kuiper (1986), Carrera (1988), Keidser and Rosbjerg (1991), Sun (2004), Carrera et al. (2005), and Nilsson et al. (2007). Calibration is the most challenging and time-consuming task of groundwater flow modeling (Hill and Tideman 2007; Cao et al. 2006) and involves optimizing the model parameters to honor groundwater heads, fluxes or concentration observations (McLaughlin and Townley 1996; Kitanidis 1997; Carrera et al. 2005; Pinault and Schomburgk 2006). Calibration of a groundwater model consists of the determination of parameters such as hydraulic conductivity, specific yield, storage coefficient, dispersivity and system geometry. Sometimes also recharge and other boundary conditions, as well as initial conditions are considered as parameters. The main benefit of calibrating hydraulic conductivities is that they often contain the dominant source of uncertainty. If measured heads and fluxes are the only input data, direct inversion may become mathematically unstable as small errors in head data can result in large errors in hydraulic conductivity (Sun 2004). However, for the more standard indirect parameter estimation Groundwater
1
approach to inverse problems the issue is not primarily with classic definitions of stability. When both conductivities and fluxes are estimated while only measured heads are specified, only ratios between conductivities and fluxes can be determined. Only if flux data are specified conductivities and fluxes can be determined separately (Haitjema 2006; Poeter and Hill 1997). More importantly, when many parameters (e.g., many grid block conductivities) have to be estimated while only a few measured heads and nonzero fluxes can be specified, the parameter estimation problem is mathematically underdetermined, thus ill-posed. Therefore, the main issue of calibration is parameter nonuniqueness, which can be reduced by putting additional constraints on the problem. Calibration is performed by trial and error, or automatically by minimizing the difference between the simulated and measured hydraulic heads and hence these are indirect inverse methods. Indirect methods are methods that can handle any type of parameterized process, independent from the mathematical equations and boundary conditions that govern the process (Sun 2004; Hill and Tideman 2007; Doherty and Hunt 2009, 2010a). Indirect calibration can also use any type of data as input, among which soft data (Tikhonov 1963a, 1963b; Doherty 2003; Hunt et al. 2007). A gradient method is a popular choice for automatic calibration. It has been used in, for example, the pilot point and sequential self-calibration methods (Fasanino et al. 1986; RamaRao et al. 1995; LaVenue et al. 1995; Wen et al. 1998, 1999; Doherty 2003). Gradient methods have good convergence rates and can, therefore, find the parameters that match the measured data very accurately (Sun 2004; Hager and Zhang 2006; Ewald 2006). Several reviews and comparisons of indirect inverse methods investigated the precision and robustness of these methods. Zimmerman et al. (1998) compared seven geostatistical methods and how suited they are for making probabilistic forecasts of solute transport through aquifers with spatial variable hydraulic conductivities. Hendricks Franssen et al. (2009) review and compare seven recent inverse approaches. These reviews confirm that some methods perform well for a given type of heterogeneity, while they perform poor for another partly due to subjective decisions in translating the conceptual to numerical calibration processes (Rajanayaka and Kulasiri 2001). It is also apparent that there is no agreement on how to select the most suitable method for a given problem. In contrast to indirect inverse methods, direct inverse methods determine hydraulic conductivities directly from the groundwater flow equations, comparable to solving the forward problem for the hydraulic head directly from the flow equations. In contrast to other disciplines like medical imaging and geophysics, there is a widespread believe among hydrogeologists that direct calibration methods “do not work.” Although direct methods to determine conductivity fields have already been proposed in the 1950s, Nelson (1962, 1968) was the first to present a practical 2
M.A. El-Rawy et al. Groundwater
application, followed by Emsellem and de Marsily (1971), Neuman (1973), and Sagar et al. (1975). Neglecting specific and phreatic storage, Darcy’s law was combined with the mass conservation equation to obtain an equation relating the log-conductivity field to the head field. It was shown that only if the head is known everywhere in the modeling domain, and if the specific flux is known at the in- or outflow boundaries, integration along the streamlines results in a successful estimation of the conductivity field. Calder´on (1980) and Borcea (2002) showed that for two- and three-dimensional flow a sufficiently smooth conductivity field is uniquely determined by both the head and the specific flux on the boundary. The DCM has been proposed as one of the methods to solve the Calderon problem (Borcea 2002). This method was originally proposed by Wexler et al. (1985) and has subsequently been applied in the context of electrical impedance tomography for medical and geophysical applications (Wexler 1988; Yorkey and Webster 1987; Kohn and Vogelius 1987; Kohn and McKenney 1990). Brouwer et al. (2008), Trykozko et al. (2008, 2009), El-Rawy et al. (2010), and Zijl et al. (2010) introduced DCM as a calibration alternative in groundwater flow modeling problems. A direct approach based on Darcy’s law has also been used for upscaling hydraulic conductivities (Warren and Price 1961; Durlofsky 1991; Zijl and Trykozko 2001). In the remainder of this paper, the Methodology section presents the essentials of the DCM, as well as the KF, which will be used to sequentially update the DCM determined hydraulic conductivity estimates. This section closes with explaining how the “measured” recharge fluxes are obtained. The Application and Results section presents an example application of the DCMKF for the Kleine Nete catchment, Belgium. It discusses the role of the uncertainty in the initial conductivities and how the methodology allows to determine areas where the hydraulic conductivity can be and cannot be updated. The last section presents conclusions as well as recommendations for future research.
Methodology In contrast to the direct inverse approach of Nelson (1962, 1968) we do not attempt to determine the complete conductivity field from measured hydraulic data. Instead, our aim is to calibrate an initially specified conductivity field (i.e., soft data). DCM is a relatively simple and instructive method to assist hydrogeologists in the calibration of grid block conductivities in a groundwater flow model (e.g., MODFLOW). Keeping extension to a three-dimensional MODFLOW model in mind, the basic steps of DCM for one-dimensional flow are 1. Assign initial (old) hydraulic conductivity K old to each grid block. 2a. Run the model with all known flow rates through the boundaries (e.g., recharge rates on the top boundary) NGWA.org
K old = K new
Specify old conductivities Kold
Specify all known fluxes (recharge rates)
Specify all known heads (piezometers) MODFLOW Flux run
MODFLOW Head run
F
H
Fluxes Q from flux run
Fluxes Q from head run
Update rule: K
new
old
F
H
= K (Q ) /(Q )
Correct anisotropy ratios No
Stop check criterion Yes Stop iterations
Figure 1. Calibration of grid block conductivities by the DCM.
and into the wells (the internal boundaries). This socalled flux run results in the specific flux q F (L/T) through the grid blocks. 2b. Run the model again, now with all known heads on the boundaries and in the wells including the piezometers. This so-called head run results in the head gradient ∂h H /∂x in the grid blocks and specific flux q H (L/T) through the grid blocks. 3. Apply Darcy’s law to replace the old grid block conductivities by the new ones K new = − q F /(∂h H /∂x ). A groundwater flow model based on these new grid block conductivities honors both the flux and the head conditions on the boundary and in the piezometers; therefore, this model may be considered as calibrated (Brouwer et al. 2008). Substitution of q H = − K old ∂h H /∂x and Q = qA, where A (L2 ) is the grid block’s cross section and Q (L3 /T) is the flux, yields a result that characterizes and explains the DCM in the update rule K new = K old
QF QH
(1)
where Q F and Q H are calculated at the block level. Basic equations; assumptions; avoidance of negative K new values; dealing with almost zero flow; anisotropy; straightforward generalization to 3-dimensional flow; and proof of convergence is described in the online Supporting Information. Iteration means that new conductivities are considered as old ones, after which another flux model and another head model are run, until a minimum head difference is found, usually 2 to 5 runs (Figure 1). Taking into account the inaccuracy of the specified heads (e.g., NGWA.org
measured in piezometers) and the specified fluxes (e.g., the “measured” recharge rates), this solution is generally sufficiently accurate. Repeating DCM under different flow conditions will result in a time series of estimated conductivities. If the specified (“measured”) fluxes and heads were error-free, the conductivities would be the same for all repeated DCM estimations. For each grid block n = 1, . . . , N with conductivities Kxn , Kyn , Kzn we considera time series of M calibrated logconductivities zmn,i = ln Kin m , m =1, . . . , M , with mean n,i log-conductivities z n,i = (1/M ) M m=1 zm and observation errors M
2 1 n,i zmn,i − z n,i (2) s = M −1 m=1
Following a suggestion by Brouwer et al. (2008) we apply a KF to obtain estimates xmn,i of the grid block logconductivities with estimation uncertainties σmn,i . Because the KF incorporates Bayesian inference rules for conditional probabilities, the estimation uncertainties will become appreciably smaller than the observation errors s n,i when the number of DCM repetitions increases. In other words, although the observation error remains constant, the reliability of the estimates increases when the estimates are reproduced more and more times by the KF. In the online Supporting Information we describe in more detail how the Linear KF can be applied. In the DCM the reliability of the calibration is determined by the accuracy of the measured head-flux couples. As fluxes have to be determined (“measured”) in a sufficiently accurate way, recharge was estimated M.A. El-Rawy et al. Groundwater
3
Figure 2. Location, topography, and rivers of the Kleine Nete basin.
using WetSpa, a spatially distributed rainfall-runoff model (Liu et al. 2003; Dams et al. 2008, 2012). In the online Supporting Information we describe in more detail the WetSpa recharge methodology.
Application and Results The Kleine Nete Catchment The Kleine Nete catchment is situated about 65 km northeast of Brussels, Belgium, and it comprises 581 km2 (Figure 2). The topography of the basin ranges from 3 to 48 m, while the mean slope is 0.36%. The dominant soil texture is sand, though in the valleys some loamy sand, sandy loam and sandy clay is present. The annual precipitation in the area is about 840 mm. The distributed recharge was for 1992 half-monthly estimated with WetSpa, resulting in a total recharge of about 500,000 m3 /d. Heads were measured half-monthly (24 times for 1992) in the piezometers, and there are 565 pumping wells abstracting a total of about 50,000 m3 /d. Groundwater flow is simulated with a MODFLOW model (Harbaugh et al. 2000) similar to Dams et al. (2008, 2012). The model has two layers, horizontal grid block dimensions of 50 m and a modeling domain of 511 rows and 722 columns. WetSpa estimated recharge is 4
M.A. El-Rawy et al. Groundwater
assigned by the RCH package while the pumping wells are handled by the WELL package. The watershed boundaries of the model are assumed to represent no-flow conditions. The distributed hydraulic conductivity for the aquifers is based on the hydrogeological classification system for Flanders (HCOV) (Table S1) (Cools et al. 2006). Figure 3 shows the occurrence of the Tertiary aquifers. Model layer 1 includes the Quaternary, Pleistocene, Pliocene aquifers, while model layer 2 represents the productive Miocene aquifer (HCOV 025, Figure 3b). Rivers, canals and lakes are simulated with MODFLOW’s RIVER package, while groundwater drainage to the small streams and wetlands is simulated by the DRAIN package. The drainage level is set to the highest elevation in the soil profiles where oxidation appears, which is assumed to be equal to the level that the groundwater will reach before it is discharged (Batelaan and De Smedt 2004). In the DCM-estimation the anisotropy ratios were kept within the interval 1.1 ≤ K x /K z ≈ K y /K z ≤ 10, which is reasonable for the Kleine Nete catchment (El-Rawy et al. 2010; El-Rawy 2013). Case 1: Soft Data Without Uncertainty In this case the DCM is applied with the same soft data (same initial conductivities) for each repetition. Only the measured head-flux couples are different at each NGWA.org
(a)
(b)
Figure 3. (a) Presence of the Tertiary aquifers, geologic faults, and piezometers. (b) Cross-section along B-B -B showing the different Tertiary formations (adapted from Dams et al. 2012). The HCOV codes are given in Table S1.
(a)
(b)
Figure 4. Mean values of log-conductivities for 24 observation maps for case 1 with no soft data uncertainty: (a) top layer, (b) bottom layer.
observation time. For each of the 737,884 grid blocks a time series of 24 different log-conductivities has been “observed” by the DCM (a log-conductivity is defined as x = ln(K /K ref ), with K ref = 1 m/d). The maps of mean values of the log-conductivities are presented in Figure 4a and 4b for respectively the top and bottom layer and range respectively from 3.9 to 0.01 with an average conductivity of 15.03 m/d (Figure 4a), and from 2.65 to 0.002 with an average conductivity of 4.85 m/d (Figure 4b). The standard deviations, or observation uncertainties, have been calculated from Equation 2 and are presented in Figure 5. The observation uncertainties range from 1.73 to 0.0 in the top layer and from 1.15 to 0.0 in the bottom layer. Grid blocks with (almost) zero variance indicate uncertainties in the “terra incognita” regions, that is, locations so far away from the measurement locations that the hydraulic conductivities are hardly updated by DCM. However, although not taken into account, the uncertainty of the soft data (the initial conductivities) will generally be much larger than the uncertainty in the conductivities that are updated by the DCM calibration. Hence, only the uncertainties in the “measurement ranges,” that is, in the grid blocks close to measured head-flux couples are realistic. NGWA.org
From Figure 5a and 5b we can conclude that the bottom layer shows slightly lower uncertainty than the top layer, but the difference is not pronounced. The lower uncertainty is due to the fact that the pumping wells are located in the bottom layer, where flow induced by the relatively well known abstraction rates dominates the flow induced by the relatively uncertain recharge. The above presented case illustrates that DCM is attractive because it is a conceptually simple method that can easily be implemented and used in combination with a standard groundwater flow model. The DCM iterations are terminated after having reached a minimum in difference between measured and calculated heads (or a minimum of the Darcy residual; see the Supporting Information). Such a minimum is reached after a few iterations. Looking for smaller minimums or the global minimum is not required in view of the observation inaccuracy obtained from the time series of estimated conductivities. Convergence to high accuracy is not required in realistic hydrogeological problems; it does not make sense to obtain an accuracy that is much higher than the accuracy introduced by the uncertain head and the generally even more uncertain flux estimates. Therefore, in our applications we used an initial calculation followed by two iterations. M.A. El-Rawy et al. Groundwater
5
(a)
(b)
Figure 5. Observation uncertainties (standard deviations) of hydraulic conductivities for case 1 with no soft data uncertainty: (a) top layer, (b) bottom layer.
The log-conductivities estimated by the KF range from 3.90 to 0.01 with an average conductivity of 15.18 m/d in the top layer, and from 2.65 to 0.002 with an average conductivity of 4.85 m/d in the bottom layer. Comparing these results with the results of the mean values of the log-conductivities of the 24 observations (Figure 4), we find that the mean value and the KF estimate are almost the same. In Figure 6 the initial conductivities, the timeaveraged DCM-estimated conductivities and the DCM-KF estimated hydraulic conductivities are shown at the locations of the piezometers in the top (27 piezometers; Figure 6a) and the bottom layer (34 piezometers; Figure 6b). The figure shows that in the measurement ranges, that is, the regions around the piezometers the conductivities determined by the DCM and by the DCM-KF do hardly differ from each other, but differ from the initial conductivities. The standard deviations of the estimation uncertainty range between 0.4 and 0.0 in the top layer (Figure 7a) and between 0.26 and 0.0 in the bottom layer (Figure 7b). As a result of the KF these standard deviations are about 77% smaller than the standard deviations of the observation uncertainty of the log-conductivities obtained by the DCM. Figure 6 also shows a comparison of the standard deviations of the log-conductivities obtained from the time-averaged DCM with the standard deviations obtained from the DCM-KF estimate at locations of piezometers. At piezometer 1-0081 and 1-0098 the hydraulic conductivity has much less uncertainty than at other piezometers, even before applying the KF. On the other hand, piezometer 1-0402 has the highest uncertainty in the log-conductivities. The average standard deviation for all piezometers in the top layer is 0.92 after time-averaged DCM and 0.22 after the KF. For the bottom layer these standard deviations are 0.71 and 0.16, respectively. It is clear that the standard deviation averaged over the piezometers in the top layer is higher than the standard deviations averaged over the piezometers in the bottom layer. This is consistent with the deeper flux argument presented above. 6
M.A. El-Rawy et al. Groundwater
In this case, we have assumed that there is no uncertainty in the soft data, that is, we have kept the initial conductivities the same for each DCM conductivity estimation. As a consequence, we found the greatest variances in the grid blocks near the piezometers. This is caused by the fact that conductivities in the “terra incognita” will hardly be updated by DCM. Hence, KF calculations make their variances very small or zero. Case 2: Soft Data with Uncertainty In practice, the uncertainty in the soft data should be significantly greater than the observation error found near the piezometers. A constructive way to avoid misinterpretation of the variances in the “terra incognita” is to determine realistic variances by specifying at the different observation times different possible distributions of the soft conductivity data. In this case the 24 applications of DCM are not only based on 24 different measured head-flux couples, but also on 24 different soft data sets. To demonstrate the principle we have chosen an extremely simple series of different initial conductivity patterns. For each grid block i the different initial patterns have a mean logconductivity pattern x i (i = 1, 2, . . . , 737,884) that is equal to the pattern used in case 1 without soft data uncertainty. One of the simplest ways to obtain soft data set m is xmi = x i (1 + c (m − (M + 1) /2)), where m = 1, 2, . . . , M is the distribution number and c is a constant. i . We observe that in our case, where M = 24, x i = xm=12.5 i i For c = 0.045 we obtain xm = x (1 + 0.045 (m − 12.5)). Admittedly, this is not a realistic random distribution, but we consider this simple approach as sufficient to show the effect of uncertainty in the initial conductivities, that is, the uncertainty in the soft data. Figure 8a and 8b presents the mean values of the log-conductivities calculated for the top and the bottom layer respectively. They range from 3.83 to 0.005 with an average conductivity of 14.73 m/d in the top layer and from 2.65 to 0.002 with an average conductivity of 4.76 m/d in the bottom layer. Comparison of these log-conductivities with those of the case with zero soft NGWA.org
(a)
(b)
Figure 6. Horizontal hydraulic conductivities at locations of piezometers in the top (a) and bottom (b) layer: initial conductivities K0 (blue); after the time-averaged DCM (red) with their standard deviations; and after the KF (black) with their standard deviations. The error bars represent standard deviation.
(a)
(b)
Figure 7. Estimation uncertainties (standard deviations) of hydraulic conductivities for case 1 with no soft data uncertainty: (a) top layer; (b) bottom layer.
data error (Figure 4a and 4b) shows that there are small differences. The observation uncertainties (Figure 9a and 9b) range from 1.76 to 0.004 and from 1.16 to 0.0 for the top and bottom layer respectively. We can conclude that, like in case 1, the top layer has more uncertainty than the bottom layer, due to the uncertainty in recharge, which in this case is added to the uncertainty in the soft data. NGWA.org
Comparison with Figure 5a and 5b shows that the red regions are the same for the case with no soft data uncertainty (c = 0) and the case with appreciable soft data uncertainty (c = 0.045). As a consequence, in these regions the variances represent a realistic estimation. In the regions further away from the red regions the presented observation uncertainties, blue in case of zero soft data uncertainty and yellow/green in case of M.A. El-Rawy et al. Groundwater
7
(a)
(b)
Figure 8. Mean values of log-conductivities for 24 observations, for case 2 with soft data uncertainty c= 0.045: (a) top layer, (b) bottom layer.
(a)
(b)
Figure 9. Observation uncertainties (standard deviation) of hydraulic conductivities for case 2 with soft data uncertainty c = 0.045: (a) top layer; (b) bottom layer.
appreciable soft data uncertainty, differ considerably. The regions with relatively large uncertainty difference between case 1 (low uncertainty) and case 2 (relatively high uncertainty) represent the “terra incognita,” that is, regions where, on basis of the available piezometers, the inverse model cannot improve the perceived hydraulic conductivities. To assign conductivities to these regions the hydrogeologist has to use additional knowledge, for instance knowledge about the subsurface’s facies distributions, or information about continuity of layers obtained from boreholes, etc., often combined with a technique to determine conductivity from grain size data obtained from samples (Rogiers et al. 2012). There are several geostatistical approaches that incorporate already available geological knowledge for solving the inverse problem in groundwater (Doherty 2003; Doherty et al. 2010). In this respect multiple-point geostatistical models based on geological analogs might be especially helpful (Huysmans et al. 2008; Huysmans and Dassargues 2009, 2011). The log-conductivities estimated by the KF for the 24 observations range from 3.89 to 0.004 with an average conductivity of 15.18 m/d for the top layer, and from 8
M.A. El-Rawy et al. Groundwater
2.65 to 0.002 with an average conductivity of 4.83 m/d for the bottom layer. If we compare these results with the results of the mean values of the log-conductivities of the 24 observations (Figure 8a and 8b), we do not observe big differences. In order to clarify this point, we have calculated the layer-averages of the horizontal hydraulic conductivities at the locations of the piezometers for case 1 and case 2. Table 1 shows clearly that in case 2 the layer-averages of the conductivities at the locations of the piezometers have smaller values than those averages in case 1, both for the time-averaged DCM and for the DCM-KF estimates. It is clear that in the “measurement ranges,” that is, at the locations of the piezometers, the time-averaged DCM conductivities, as well as the estimated DCM-KF conductivities have been changed with respect to the initial hydraulic conductivities. It is also shown that the differences between the time-averaged DCM conductivities and the estimated DCM-KF conductivities are relatively small. Comparing the hydraulic conductivities of case 2 (c = 0.045) with case 1 (c = 0), we can conclude that at the locations of the piezometers both the DCM and the DCM-KF conductivities in case 2 have a lower value NGWA.org
Table 1 Layer-Average of the Horizontal Hydraulic Conductivity at the Locations of the Piezometers for the Top and the Bottom Layer: Initial Conductivities, DCM-Estimated Conductivities and DCM-KF-Estimated Conductivities for Cases 1 and 2 Kh avg (m/d) Top Layer
Before DCM After DCM After KF
Bottom Layer
Case 1
Case 2
Case 1
Case 2
16.09 6.63 6.99
16.09 6.09 6.27
7.24 4.04 4.16
7.24 3.67 3.79
Table 2 Mean Absolute Error (MAE) and the Root Mean Square Error (RMSE) for Initial, DCM-Estimated and DCM-KF-Estimated Conductivities with Zero and Nonzero Soft Data Uncertainty Initial
DCM
DCM-KF
Zero Nonzero Zero Nonzero Zero Nonzero MAE 0.8 RMSE 1.13
0.8 1.13
0.6 0.9
0.74 1.06
0.6 0.89
0.73 1.05
with respect to the initial conductivities than those of case 1. The estimation uncertainties calculated by the KF, that is, the standard deviations of the estimates, range from 0.38 to 0.001 for the top layer (Figure 10a) and from 0.27 to 0.0 for the bottom layer (Figure 10b). The estimation uncertainty calculated by the KF is about 77% smaller than the observation uncertainty. Differences Between the Two Cases In order to calculate the differences between the observed and simulated heads we have used a forward MODFLOW model with flux boundary conditions and three different conductivity fields: (1) the initial grid block conductivities (the soft data); (2) the time-averaged conductivities obtained from DCM-calibration; (3) the conductivities estimated by the KF from the 24 DCMobservations. MODFLOW is run for the hypothetical case where the recharge flux is equal to the average of the 24 time-dependent recharge rates. The calculated heads are then compared with the time-averaged measured heads. The mean absolute error (MAE) and the root mean square error (RMSE) for the cases 1 and 2 are presented in Table 2. For the case 1 with zero soft data uncertainty we found that the MAE and RMSE obtained from the KF estimates are relatively small, 0.60 m and 0.89 m, respectively. These errors deviate considerably from the errors obtained from the model with initial grid block NGWA.org
conductivities (MAE = 0.8 m and RMSE = 1.13 m). This shows that the estimated conductivities do represent an improvement for this time-averaged case. This is not surprising as the KF estimate is a time-averaged value. In case 2 with non-zero uncertainty of the soft data (Table 2), the MAE and RMSE values before and after application of the KF are slightly higher than the values for the case with zero soft data uncertainty. This is expected, as uncertainties in the case with nonzero soft data uncertainty are higher and more realistic than uncertainties in the case with zero soft data uncertainty. Measurement Ranges and Terra Incognita Using both DCM hydraulic conductivity estimates resulting from the initial conductivity pattern of case 1 and the series of different initial conductivity patterns of case 2, it is possible to assess in which regions conductivities can be determined by the hydraulic data and in which regions this is impossible. To delineate these regions, we have determined for each grid block the relative difference between the standard deviations of case 1 (without soft data uncertainty) and the standard deviations of case 2 (with soft data uncertainty) (Figure 11a and 11b). In the “measurement ranges” near the piezometers the observed conductivities and their uncertainties (both observation and estimation uncertainties) are approximately the same for cases 1 and 2 (the grey regions in Figure 11). For these regions, in which the standard deviations are mainly determined by the uncertain head-flux measurements, we may assume that the results represent a realistic conductivity estimate. However, in the regions further away from the piezometers (from the orange to blue regions in Figure 11) the standard deviations differ considerably between case 1 and case 2, representing the “terra incognita.” In these regions additional hydraulic measurements (piezometers) are needed to find a more reliable conductivity pattern.
Conclusions 1. The DCM consists of alternating flux and head specified (MODFLOW) runs using observations of head-flux couples; the difference between the runs provides insight in areas where the hydraulic conductivity can and cannot be adjusted from its initial value. The DCM is based on two (MODFLOW) model runs based on the same initially chosen grid block conductivities: •
The flux run has as many specified flux data as available (abstraction well rates, indirectly “measured” recharge rates derived from measured rainfall, runoff and evapotranspiration). • The head run has as many specified head data as available (including the heads measured in piezometers), like in T´oth’s (2009) flow systems approach, where instead M.A. El-Rawy et al. Groundwater
9
(a)
(b)
Figure 10. Estimation uncertainties (standard deviation) σ of hydraulic conductivities for case 2 with soft data uncertainty c = 0.045: (a) top layer, (b) bottom layer.
(a)
(b)
Figure 11. Relative differences between the standard deviations of case 1 and case 2: (a) top layer, (b) bottom layer.
of the recharge the head is specified on the phreatic surface. The differences in solution of the flux and the head model give insight in the so-called “terra incognita,” that is, the regions where the “measured” head-flux couples cannot improve the initial hydraulic conductivity, and in the so-called “measurement ranges,” that is, the regions where the “measured” head-flux couples do adjust the initial conductivity. In addition, the differences in the solution of the two runs can be used to update the grid block conductivities by application of Darcy’s law. As a consequence, the DCM is essentially based on the principle that calibration is meaningful only if we can specify more head and flux data than required for a forward problem. That is, we need to know head-flux couples in a number of wells or at a number of boundary points. These hydraulic data have to be acquired from sources independent from the model to be calibrated. 2. DCM is cheap and can support trial and error or PEST calibration. 10
M.A. El-Rawy et al. Groundwater
The DCM is relatively cheap in computation time and memory requirements. Also if we apply the “trialand-error” or a more comprehensive indirect calibration method like PEST (Doherty and Hunt 2010b), running both a flux model and a head model may enhance insight in the calibration problem and can assist in delineating a meaningful zonation. In addition, DCMestimated conductivities in the measurement ranges could be used as pilot points (Certes and de Marsily 1991), for instance in a subsequent PEST calibration (Hunt et al. 2007). In this sense the DCM and its underlying idea of a “flux run” and a “head run” may be considered as a calibration assisting tool. 3. Repeated DCM estimations with KF postprocessing, results in conductivities with lower uncertainty. The calibration assisting methodology can be enhanced by repeated DCM estimations under different hydrological conditions. For instance, estimation in the winter season with relatively high rainfall and high water levels, and estimation in the summer season with NGWA.org
relatively little rainfall and low water levels. If the “headflux” measurements were without measurement error, we would find the same conductivities for all hydrological conditions. In reality we find a series of time-dependent conductivities, from which we can determine the time-averaged mean values and standard deviation or observation error. Using this observation error as input of a KF, an estimation of the time-independent grid block conductivities can be obtained from the time series of observed (DCM calibrated) conductivities. The KF then determines for the estimated conductivities an uncertainty that is appreciably smaller than the observation error. 4. DCM-KF is a valuable calibration tool and enhances hydrological understanding. In conclusion we may say that the calibration strategy underlying the DCM and subsequent application of the KF works well while enhancing hydrological understanding. Therefore this approach to calibration may be considered as a valuable tool to assist hydrogeologists in their calibration of the hydraulic conductivities.
Acknowledgments The authors gratefully acknowledge Randall J. Hunt and two anonymous reviewers for their constructive review, which considerably improved this paper. Also James McCallum is acknowledged for suggestions for improvement. The first author acknowledges Erasmus Mundus External Cooperation Window and Vrije Universiteit Brussel for the funding of his scholarship.
Supporting Information Additional Supporting Information may be found in the online version of this article: Figure S1. Flow through one-dimensional domain partitioned in sections with lengths L1 and L2 , conductivities K 1 and K 2 . Table S1. Overview of the Hydrostratigraphy of the Study Area
References Batelaan, O., and F. De Smedt. 2004. SEEPAGE, a new MODFLOW DRAIN package. Ground Water 42, no. 4: 576–588. Borcea, L. 2002. Electrical impedance tomography. Inverse Problems 18, no. 6: R99–R136. Brouwer, G.K., P.A. Fokker, F. Wilschut, and W. Zijl. 2008. A direct inverse model to determine permeability fields from pressure and flow rate measurement. Mathematical Geosciences 40, no. 8: 907–920. Calder´on, A.P. 1980. On an Inverse Boundary Value Problem, Seminar on Numerical Analysis and its Applications to Continuum Physics (Rio de Janeiro, 1980), 65–73. Rio de Janeiro: Soc. Brasil Mat. Cao, W., W.B. Bowden, T. Davie, and A. Fenemor. 2006. Multi-variable and multi-site calibration and validation of
NGWA.org
SWAT in a large mountainous catchment with high spatial variability. Hydrological Processes 20, no. 5: 1057–1073. Carrera, J. 1988. State of the art of the inverse problem applied to the flow and solute transport equations. In Groundwater Flow and Quality Modelling, ed. E. Custodio, A. Guirgui, and J.P.L. Ferreira, 549–583. Netherlands: D. Reidel Publishing Company/Springer. Carrera, J., A. Alcolea, A. Medina, J. Hidalgo, and L.J. Slooten. 2005. Inverse problem in hydrogeology. Hydrogeology Journal 13, no. 1: 206–222. Certes, C., and G. de Marsily. 1991. Application of the pilot points method to the identification of aquifer transmissivities. Advances in Water Resources 14, no. 5: 285–300. Cools, J., Y. Meyus, S.T. Woldeamlak, O. Batelaan, and F. De Smedt. 2006. Large-scale GIS-based hydrogeological modelling of Flanders: A tool for groundwater management. Environmental Geology 50, no. 8: 1201–1209. Dams, J., E. Salvadore, T. Van Daele, V. Ntegeka, P. Willems, and O. Batelaan. 2012. Spatio-temporal impact of climate change on the groundwater system. Hydrology Earth System Sciences 16: 1517–1531. DOI:10.5194/hess-16-1517-2012. Dams, J., S.T. Woldeamlak, and O. Batelaan. 2008. Predicting land-use change and its impact on the groundwater system of the Kleine Nete catchment, Belgium. Hydrology Earth System Sciences 12: 1369–1385. DOI:10.5194/hess-121369-2008. de Marsily, G., J.P. Delhomme, A. Coudrain-Ribstein, and A.M. Lavenue. 2000. Four decades of inverse problems in hydrogeology. Geological Society of America. Special Paper 348. Doherty, J. 2003. Ground water model calibration using pilot points and regularization. Ground Water 41, no. 2: 170–177. Doherty, J., and R.J. Hunt. 2010a. Response to comment on: Two statistics for evaluating parameter identifiability and error reduction. Journal of Hydrology 380: 489–496. DOI:10.1016/j.jhydrol.2009.10.012. Doherty, J.E., and R.J. Hunt. 2010b. Approaches to highly parameterized inversion: A guide to using PEST for groundwater-model calibration. U.S. Geological Survey Scientific Investigations Report 2010-5169, 60. http://pubs.usgs.gov/sir/2010/5169/. Doherty, J.E., R.J. Hunt, and M.J. Tonkin. 2010. Approaches to highly parameterized inversion: A guide to using PEST for model-parameter and predictive-uncertainty analysis. U.S. Geological Survey Scientific Investigations Report 20105211, 71. Doherty, J., and R.J. Hunt. 2009. Two statistics for evaluating parameter identifiability and error reduction. Journal of Hydrology 366: 119–127. DOI:10.1016/j.jhydrol.2008.12.018. Durlofsky, L.J. 1991. Numerical calculation of equivalent grid block permeability tensors for heterogeneous porous media. Water Resources Research 27: 699–708. El-Rawy, M.A. 2013. Calibration of hydraulic conductivities in groundwater flow models using the double constraint method and the Kalman filter. PhD thesis, Vrije Universiteit Brussel, Brussels, Belgium. El-Rawy, M.A., W. Zijl, O. Batelaan, and G.A. Mohammed. 2010. Application of the double constraint method combined with MODFLOW. In Proceedings of the Valencia IAHR International Groundwater Symposium, September 22–24. Emsellem, Y., and G. de Marsily. 1971. An automatic solution for the inverse problem. Water Resources Research 7, no. 5: 1264–1283. Ewald, C.-O. 2006. The Malliavin Gradient Method for the calibration of stochastic dynamical models. Applied Mathematics and Computer 175, no. 2: 1332–1352. Fasanino, G., J.E. Molinard, G. de Marsily, and V. Pelc´e. 1986. Inverse modeling in gas reservoirs. SPE 15592, 61st SPE
M.A. El-Rawy et al. Groundwater
11
Annual Technical Conference and Exhibition, October 5–8, New Orleans, Louisiana. Hager, W.W., and H. Zhang. 2006. A survey of nonlinear conjugate gradient methods. Pacific Journal of Optimization 2, no. 1: 35–58. Haitjema, H. 2006. The role of hand calculations in ground water flow modeling. Ground Water 44, no. 6: 786–791. DOI:10.1111/j.1745-6584.2006.00189.x. Harbaugh, A.W., E.R. Banta, M.C. Hill, and M.G. McDonald. 2000. MODFLOW-2000. The U.S. Geological Survey modular ground-water model—User guide to modularization concepts and the ground-water flow process. U.S. Geological Survey Open-File Report 00-92, 121. Reston, Virginia: USGS. Hendricks Franssen, H.J., A. Alcolea, M. Riva, M. Bakr, N. Van der Wiel, F. Stauffer, and A. Guadagnini. 2009. A comparison of seven methods for the inverse modelling of groundwater flow. Application to the characterisation of well catchments. Advances in Water Resources 32, no. 6: 851–872. Hill, M.C., and C.R. Tideman. 2007. Effective Groundwater Model Calibration: With Analysis of Data, Sensitivities, Predictions, and Uncertainty. Hoboken, New Jersey: John Wiley and Sons, Inc. Hunt, R.J., J. Doherty, and M.J. Tonkin. 2007. Are models too simple? Arguments for increased parameterization. Ground Water 45, no. 3: 254–261. DOI:10.1111/j.1745-6584.2007. 00316.x. Huysmans, M., and A. Dassargues. 2011. Direct multiple-point geostatistical simulation of edge properties for modeling thin irregularly shaped surfaces. Mathematics Geosciences 43, no. 5: 521–536. Huysmans, M., and A. Dassargues. 2009. Application of multiple-point geostatistics on modelling groundwater flow and transport in a cross-bedded aquifer (Belgium). Hydrogeology Journal 17, no. 8: 1901–1911. Huysmans, M., L. Peeters, G. Moermans, and A. Dassargues. 2008. Relating small-scale sedimentary structures and permeability in a cross-bedded aquifer. Journal of Hydrology 361, no. 1–2: 41–51. Keidser, A., and D. Rosbjerg. 1991. A comparison of four inverse approaches to groundwater flow and transport parameter Identification. Water Resources Research 27, no. 9: 2219–2232. Kitanidis, P.K. 1997. The minimum structure solution to the inverse problem. Water Resources Research 33, no. 10: 2263–2272. Kohn, R.V., and A. McKenney. 1990. Numerical implementation of a variational method for electrical impedance tomography. Inverse Problems 6, no. 3: 389–414. Kohn, R.V., and M. Vogelius. 1987. Relaxation of a variational method for impedance computed tomography. Communications on Pure Applied Mathematics 40, no. 6: 745–777. Kuiper, L.K. 1986. A comparison of several methods for the solution of the inverse problem in two-dimensional steady-state groundwater flow modeling. Water Resources Research 22, no. 5: 705–714. LaVenue, A.M., B.S. RamaRao, G. de Marsily, and M.G. Marietta. 1995. Pilot point methodology for automated calibration of an ensemble of conditionally simulated transmissivity fields, 2, Application. Water Resources Research 31, no. 3: 495–516. Liu, Y.B., S. Gebremeskel, F. De Smedt, L. Hoffmann, and L. Pfister. 2003. A diffusive transport approach for flow routing in GIS-based flood modeling. Journal of Hydrology 283, no. 1–4: 91–106. McLaughlin, D., and L.R. Townley. 1996. A reassessment of the groundwater inverse problem. Water Resources Research 32, no. 5: 1131–1161.
12
M.A. El-Rawy et al. Groundwater
Nelson, R.W. 1968. In-place determination of permeability distribution for heterogeneous porous media through analysis of energy dissipation. Society of Petroleum Engineers Journal 8, no. 1: 33–42. http://www.onepetro.org/mslib/ servlet/onepetropreview?id=00001554. Nelson, R.W. 1962. Conditions for determining areal permeability distribution by calculation. Society of Petroleum Engineers Journal 2, no. 3: 223–224. http://www.onepetro.org/ mslib/servlet/onepetropreview?id=00000371. Neuman, S.P. 1973. Calibration of distributed parameter groundwater flow models viewed as a multiple-objective decision process under uncertainty. Water Resources Research 9, no. 4: 1006–1021. Nilsson, B., A.L. Højberg, J.C. Refsgaard, and L. Troldborg. 2007. Uncertainty in geological and hydrogeological data. Hydrology Earth System Sciences 11: 1551–1561. Pinault, J.L., and S. Schomburgk. 2006. Inverse modeling for characterizing surface water/groundwater exchanges. Water Resources Research 42, no. 8. DOI:10.1029/2005WR 004587. Poeter, E.P., and M.C. Hill. 1997. Inverse models-A necessary next step in ground-water modeling. Ground Water 35, no. 2: 250–260. Rajanayaka, C., and D. Kulasiri. 2001. Investigation of a parameter estimation method for contaminant transport in aquifers. Journal of Hydroinformatics 3: 203–213. RamaRao, B.S., A.M. LaVenue, G. de Marsily, and M.G. Marietta. 1995. Pilot point methodology for automated calibration of an ensemble of conditionally simulated transmissivity fields: 1. Theory and computational experiments. Water Resources Research 31, no. 3: 475–493. Rogiers, B., D. Mallants, O. Batelaan, M. Gedeon, M. Huysmans, and A. Dassargues. 2012. Estimation of hydraulic conductivity and its uncertainty from grain-size data using GLUE and artificial neural networks. Mathematical Geosciences 44, no. 6: 739–763. Sagar, B., S. Yakowitz, and L. Duckstein. 1975. A direct method for the identification of the parameters of dynamic nonhomogeneous aquifers. Water Resources Research 11, no. 4: 563–570. Sun, N.-Z. 2004. Inverse Problems in Groundwater Modeling. Dordrecht: Kluwer Academic Publishers. Tikhonov, A.N. 1963a. Solution of incorrectly formulated problems and the regularization method. Soviet Mathematics Doklady 4: 1035–1038. Tikhonov, A.N. 1963b. Regularization of incorrectly posed problems. Soviet Mathematics Doklady 4: 1624–1637. T´oth, J. 2009. Gravitational Systems of Groundwater Flow . Cambridge: Cambridge University Press. Trykozko, A., G.A. Mohammed, and W. Zijl. 2009. Downscaling: The inverse of upscaling. In Conference on Mathematical and Computational Issues in the Geosciences. SIAM GS 2009 , June 15–18, Leipzig. Trykozko, A., G.K. Brouwer, and W. Zijl. 2008. Downscaling: A complement to homogenization. International Journal of Numerical Analysis and Modeling 5, no. supp: 157–170. Warren, J.E., and H.S. Price. 1961. Flow in heterogeneous porous media. Society of Petroleum Engineers Journal 1, no. 3: 153–169. DOI:10.2118/1579-G. Wen, X.-H., C.V. Deutsch, J.J. Gomez-Hernandez, and A.S. Cullick. 1999. A program to create permeability fields that honor single-phase flow rate and pressure data. Computers and Geosciences 25, no. 3: 217–230. Wen, X.-H., C.V. Deutsch, and A.S. Cullick. 1998. Highresolution reservoir models integrating multiple-well production data. Society of Petroleum Engineers Journal 3, no. 4: 344–355. Wexler, A. 1988. Electrical impedivity imaging in two and three dimensions. Clinical Physics and Physiological Measurement 9, no. Suppl A: 29–33.
NGWA.org
Wexler, A., B. Fry, and M.R. Neuman. 1985. Impedivitycomputed tomography algorithm and system. Applied Optics 24, no. 23: 3985–3992. Yeh, W.W.-G. 1986. Review of parameter identification procedures in groundwater hydrology: The inverse problem. Water Resources Research 22, no. 2: 95–108. Yorkey, T.J., and J.G. Webster. 1987. A comparison of impedivity tomographic reconstruction algorithms. Clinical Physics and Physiological Measurement 8, no. Suppl A: 55–62. Zijl, W., G.A. Mohammed, O. Batelaan, and F. De Smedt. 2010. Constraining methods for direct inverse modeling. XVIII International Conference on Water Resources, CMWR,
NGWA.org
CIMNE, Barcelona. http://congress.cimne.com/cmwr2010/ Proceedings/docs/p62.pdf. Zijl, W., and A. Trykozko. 2001. Numerical homogenization of the absolute permeability using the conformal-nodal and mixed-hybrid finite element method. Transport in Porous Media 44: 33–62. Zimmerman, D.A., G.D. de Marsily, C.A. Gotway, M.G. Marietta, C.L. Axness, R.I. Beauheim, R.I. Bras, J. Carrera, G. Dagan, and P.B. Davies. 1998. A comparison of seven geostatistically based inverse approaches to estimate transmissivities for modeling advective transport by groundwater flow. Water Resources Research 34, no. 6: 1373–1413.
M.A. El-Rawy et al. Groundwater
13
Abstract This paper presents the Kalman Filtered Double Constraint Method (DCM-KF) as a technique to estimate the hydraulic conductivities in the grid blocks of a groundwater flow model. The DCM is based on two forward runs with the same initial grid block conductivities, but with alternating flux-head conditions specified on parts of the boundary and the wells. These two runs are defined as: (1) the flux run, with specified fluxes (recharge and well abstractions), and (2) the head run, with specified heads (measured in piezometers). Conductivities are then estimated as the initial conductivities multiplied by the fluxes obtained from the flux run and divided by the fluxes obtained from the head run. The DCM is easy to implement in combination with existing models (e.g., MODFLOW). Sufficiently accurate conductivities are obtained after a few iterations. Because of errors in the specified head-flux couples, repeated estimation under varying hydrological conditions results in different conductivities. A time-independent estimate of the conductivities and their inaccuracy can be obtained by a simple linear KF with modest computational requirements. For the Kleine Nete catchment, Belgium, the DCM-KF yields sufficiently accurate calibrated conductivities. The method also results in distinguishing regions where the head-flux observations influence the calibration from areas where it is not able to influence the hydraulic conductivity.
Introduction Do you calibrate the hydraulic conductivity in your groundwater flow model based on head observations, while you specify the recharge? Then you can estimate the hydraulic conductivity in the neighbourhood of the observations by performing alternating model runs specifying head and recharge. In this paper we propose, demonstrate and discuss the Kalman Filtered Double Constraint Method (DCM-KF). The combined approach of the DCM with the KF makes a simple but effective form of direct inverse calibration possible. To model a groundwater flow system, two types of problems must be solved: (1) the forward problem (simulation) and (2) the inverse problem (calibration). Inverse problems arise in various fields such as medicine, meteorology and hydrology whenever mathematical models are used to determine unknown properties of a system from the observations of a response of that system. Usually the 1 Department of Hydrology and Hydraulic Engineering, Vrije Universiteit Brussel, Pleinlaan 2, 1050 Brussels, Belgium. 2 Corresponding author: Department of Civil Engineering, Faculty of Engineering, Minia University, Minia 61111, Egypt; fax: +2-086-2346674; [email protected] 3 National Centre for Groundwater Research and Training, School of the Environment, Flinders University, GPO Box 2100, Adelaide, SA 5001, Australia. Received August 2013, accepted March 2014. © 2014, National Ground Water Association. doi: 10.1111/gwat.12217
NGWA.org
unknown properties are represented as physical parameters characterizing the model and they are systematically adjusted in order to match the model outputs with the observed values. Inverse problems in groundwater modeling have been studied widely during the last four decades (de Marsily et al. 2000) and have been reviewed in detail by Yeh (1986), Kuiper (1986), Carrera (1988), Keidser and Rosbjerg (1991), Sun (2004), Carrera et al. (2005), and Nilsson et al. (2007). Calibration is the most challenging and time-consuming task of groundwater flow modeling (Hill and Tideman 2007; Cao et al. 2006) and involves optimizing the model parameters to honor groundwater heads, fluxes or concentration observations (McLaughlin and Townley 1996; Kitanidis 1997; Carrera et al. 2005; Pinault and Schomburgk 2006). Calibration of a groundwater model consists of the determination of parameters such as hydraulic conductivity, specific yield, storage coefficient, dispersivity and system geometry. Sometimes also recharge and other boundary conditions, as well as initial conditions are considered as parameters. The main benefit of calibrating hydraulic conductivities is that they often contain the dominant source of uncertainty. If measured heads and fluxes are the only input data, direct inversion may become mathematically unstable as small errors in head data can result in large errors in hydraulic conductivity (Sun 2004). However, for the more standard indirect parameter estimation Groundwater
1
approach to inverse problems the issue is not primarily with classic definitions of stability. When both conductivities and fluxes are estimated while only measured heads are specified, only ratios between conductivities and fluxes can be determined. Only if flux data are specified conductivities and fluxes can be determined separately (Haitjema 2006; Poeter and Hill 1997). More importantly, when many parameters (e.g., many grid block conductivities) have to be estimated while only a few measured heads and nonzero fluxes can be specified, the parameter estimation problem is mathematically underdetermined, thus ill-posed. Therefore, the main issue of calibration is parameter nonuniqueness, which can be reduced by putting additional constraints on the problem. Calibration is performed by trial and error, or automatically by minimizing the difference between the simulated and measured hydraulic heads and hence these are indirect inverse methods. Indirect methods are methods that can handle any type of parameterized process, independent from the mathematical equations and boundary conditions that govern the process (Sun 2004; Hill and Tideman 2007; Doherty and Hunt 2009, 2010a). Indirect calibration can also use any type of data as input, among which soft data (Tikhonov 1963a, 1963b; Doherty 2003; Hunt et al. 2007). A gradient method is a popular choice for automatic calibration. It has been used in, for example, the pilot point and sequential self-calibration methods (Fasanino et al. 1986; RamaRao et al. 1995; LaVenue et al. 1995; Wen et al. 1998, 1999; Doherty 2003). Gradient methods have good convergence rates and can, therefore, find the parameters that match the measured data very accurately (Sun 2004; Hager and Zhang 2006; Ewald 2006). Several reviews and comparisons of indirect inverse methods investigated the precision and robustness of these methods. Zimmerman et al. (1998) compared seven geostatistical methods and how suited they are for making probabilistic forecasts of solute transport through aquifers with spatial variable hydraulic conductivities. Hendricks Franssen et al. (2009) review and compare seven recent inverse approaches. These reviews confirm that some methods perform well for a given type of heterogeneity, while they perform poor for another partly due to subjective decisions in translating the conceptual to numerical calibration processes (Rajanayaka and Kulasiri 2001). It is also apparent that there is no agreement on how to select the most suitable method for a given problem. In contrast to indirect inverse methods, direct inverse methods determine hydraulic conductivities directly from the groundwater flow equations, comparable to solving the forward problem for the hydraulic head directly from the flow equations. In contrast to other disciplines like medical imaging and geophysics, there is a widespread believe among hydrogeologists that direct calibration methods “do not work.” Although direct methods to determine conductivity fields have already been proposed in the 1950s, Nelson (1962, 1968) was the first to present a practical 2
M.A. El-Rawy et al. Groundwater
application, followed by Emsellem and de Marsily (1971), Neuman (1973), and Sagar et al. (1975). Neglecting specific and phreatic storage, Darcy’s law was combined with the mass conservation equation to obtain an equation relating the log-conductivity field to the head field. It was shown that only if the head is known everywhere in the modeling domain, and if the specific flux is known at the in- or outflow boundaries, integration along the streamlines results in a successful estimation of the conductivity field. Calder´on (1980) and Borcea (2002) showed that for two- and three-dimensional flow a sufficiently smooth conductivity field is uniquely determined by both the head and the specific flux on the boundary. The DCM has been proposed as one of the methods to solve the Calderon problem (Borcea 2002). This method was originally proposed by Wexler et al. (1985) and has subsequently been applied in the context of electrical impedance tomography for medical and geophysical applications (Wexler 1988; Yorkey and Webster 1987; Kohn and Vogelius 1987; Kohn and McKenney 1990). Brouwer et al. (2008), Trykozko et al. (2008, 2009), El-Rawy et al. (2010), and Zijl et al. (2010) introduced DCM as a calibration alternative in groundwater flow modeling problems. A direct approach based on Darcy’s law has also been used for upscaling hydraulic conductivities (Warren and Price 1961; Durlofsky 1991; Zijl and Trykozko 2001). In the remainder of this paper, the Methodology section presents the essentials of the DCM, as well as the KF, which will be used to sequentially update the DCM determined hydraulic conductivity estimates. This section closes with explaining how the “measured” recharge fluxes are obtained. The Application and Results section presents an example application of the DCMKF for the Kleine Nete catchment, Belgium. It discusses the role of the uncertainty in the initial conductivities and how the methodology allows to determine areas where the hydraulic conductivity can be and cannot be updated. The last section presents conclusions as well as recommendations for future research.
Methodology In contrast to the direct inverse approach of Nelson (1962, 1968) we do not attempt to determine the complete conductivity field from measured hydraulic data. Instead, our aim is to calibrate an initially specified conductivity field (i.e., soft data). DCM is a relatively simple and instructive method to assist hydrogeologists in the calibration of grid block conductivities in a groundwater flow model (e.g., MODFLOW). Keeping extension to a three-dimensional MODFLOW model in mind, the basic steps of DCM for one-dimensional flow are 1. Assign initial (old) hydraulic conductivity K old to each grid block. 2a. Run the model with all known flow rates through the boundaries (e.g., recharge rates on the top boundary) NGWA.org
K old = K new
Specify old conductivities Kold
Specify all known fluxes (recharge rates)
Specify all known heads (piezometers) MODFLOW Flux run
MODFLOW Head run
F
H
Fluxes Q from flux run
Fluxes Q from head run
Update rule: K
new
old
F
H
= K (Q ) /(Q )
Correct anisotropy ratios No
Stop check criterion Yes Stop iterations
Figure 1. Calibration of grid block conductivities by the DCM.
and into the wells (the internal boundaries). This socalled flux run results in the specific flux q F (L/T) through the grid blocks. 2b. Run the model again, now with all known heads on the boundaries and in the wells including the piezometers. This so-called head run results in the head gradient ∂h H /∂x in the grid blocks and specific flux q H (L/T) through the grid blocks. 3. Apply Darcy’s law to replace the old grid block conductivities by the new ones K new = − q F /(∂h H /∂x ). A groundwater flow model based on these new grid block conductivities honors both the flux and the head conditions on the boundary and in the piezometers; therefore, this model may be considered as calibrated (Brouwer et al. 2008). Substitution of q H = − K old ∂h H /∂x and Q = qA, where A (L2 ) is the grid block’s cross section and Q (L3 /T) is the flux, yields a result that characterizes and explains the DCM in the update rule K new = K old
QF QH
(1)
where Q F and Q H are calculated at the block level. Basic equations; assumptions; avoidance of negative K new values; dealing with almost zero flow; anisotropy; straightforward generalization to 3-dimensional flow; and proof of convergence is described in the online Supporting Information. Iteration means that new conductivities are considered as old ones, after which another flux model and another head model are run, until a minimum head difference is found, usually 2 to 5 runs (Figure 1). Taking into account the inaccuracy of the specified heads (e.g., NGWA.org
measured in piezometers) and the specified fluxes (e.g., the “measured” recharge rates), this solution is generally sufficiently accurate. Repeating DCM under different flow conditions will result in a time series of estimated conductivities. If the specified (“measured”) fluxes and heads were error-free, the conductivities would be the same for all repeated DCM estimations. For each grid block n = 1, . . . , N with conductivities Kxn , Kyn , Kzn we considera time series of M calibrated logconductivities zmn,i = ln Kin m , m =1, . . . , M , with mean n,i log-conductivities z n,i = (1/M ) M m=1 zm and observation errors M
2 1 n,i zmn,i − z n,i (2) s = M −1 m=1
Following a suggestion by Brouwer et al. (2008) we apply a KF to obtain estimates xmn,i of the grid block logconductivities with estimation uncertainties σmn,i . Because the KF incorporates Bayesian inference rules for conditional probabilities, the estimation uncertainties will become appreciably smaller than the observation errors s n,i when the number of DCM repetitions increases. In other words, although the observation error remains constant, the reliability of the estimates increases when the estimates are reproduced more and more times by the KF. In the online Supporting Information we describe in more detail how the Linear KF can be applied. In the DCM the reliability of the calibration is determined by the accuracy of the measured head-flux couples. As fluxes have to be determined (“measured”) in a sufficiently accurate way, recharge was estimated M.A. El-Rawy et al. Groundwater
3
Figure 2. Location, topography, and rivers of the Kleine Nete basin.
using WetSpa, a spatially distributed rainfall-runoff model (Liu et al. 2003; Dams et al. 2008, 2012). In the online Supporting Information we describe in more detail the WetSpa recharge methodology.
Application and Results The Kleine Nete Catchment The Kleine Nete catchment is situated about 65 km northeast of Brussels, Belgium, and it comprises 581 km2 (Figure 2). The topography of the basin ranges from 3 to 48 m, while the mean slope is 0.36%. The dominant soil texture is sand, though in the valleys some loamy sand, sandy loam and sandy clay is present. The annual precipitation in the area is about 840 mm. The distributed recharge was for 1992 half-monthly estimated with WetSpa, resulting in a total recharge of about 500,000 m3 /d. Heads were measured half-monthly (24 times for 1992) in the piezometers, and there are 565 pumping wells abstracting a total of about 50,000 m3 /d. Groundwater flow is simulated with a MODFLOW model (Harbaugh et al. 2000) similar to Dams et al. (2008, 2012). The model has two layers, horizontal grid block dimensions of 50 m and a modeling domain of 511 rows and 722 columns. WetSpa estimated recharge is 4
M.A. El-Rawy et al. Groundwater
assigned by the RCH package while the pumping wells are handled by the WELL package. The watershed boundaries of the model are assumed to represent no-flow conditions. The distributed hydraulic conductivity for the aquifers is based on the hydrogeological classification system for Flanders (HCOV) (Table S1) (Cools et al. 2006). Figure 3 shows the occurrence of the Tertiary aquifers. Model layer 1 includes the Quaternary, Pleistocene, Pliocene aquifers, while model layer 2 represents the productive Miocene aquifer (HCOV 025, Figure 3b). Rivers, canals and lakes are simulated with MODFLOW’s RIVER package, while groundwater drainage to the small streams and wetlands is simulated by the DRAIN package. The drainage level is set to the highest elevation in the soil profiles where oxidation appears, which is assumed to be equal to the level that the groundwater will reach before it is discharged (Batelaan and De Smedt 2004). In the DCM-estimation the anisotropy ratios were kept within the interval 1.1 ≤ K x /K z ≈ K y /K z ≤ 10, which is reasonable for the Kleine Nete catchment (El-Rawy et al. 2010; El-Rawy 2013). Case 1: Soft Data Without Uncertainty In this case the DCM is applied with the same soft data (same initial conductivities) for each repetition. Only the measured head-flux couples are different at each NGWA.org
(a)
(b)
Figure 3. (a) Presence of the Tertiary aquifers, geologic faults, and piezometers. (b) Cross-section along B-B -B showing the different Tertiary formations (adapted from Dams et al. 2012). The HCOV codes are given in Table S1.
(a)
(b)
Figure 4. Mean values of log-conductivities for 24 observation maps for case 1 with no soft data uncertainty: (a) top layer, (b) bottom layer.
observation time. For each of the 737,884 grid blocks a time series of 24 different log-conductivities has been “observed” by the DCM (a log-conductivity is defined as x = ln(K /K ref ), with K ref = 1 m/d). The maps of mean values of the log-conductivities are presented in Figure 4a and 4b for respectively the top and bottom layer and range respectively from 3.9 to 0.01 with an average conductivity of 15.03 m/d (Figure 4a), and from 2.65 to 0.002 with an average conductivity of 4.85 m/d (Figure 4b). The standard deviations, or observation uncertainties, have been calculated from Equation 2 and are presented in Figure 5. The observation uncertainties range from 1.73 to 0.0 in the top layer and from 1.15 to 0.0 in the bottom layer. Grid blocks with (almost) zero variance indicate uncertainties in the “terra incognita” regions, that is, locations so far away from the measurement locations that the hydraulic conductivities are hardly updated by DCM. However, although not taken into account, the uncertainty of the soft data (the initial conductivities) will generally be much larger than the uncertainty in the conductivities that are updated by the DCM calibration. Hence, only the uncertainties in the “measurement ranges,” that is, in the grid blocks close to measured head-flux couples are realistic. NGWA.org
From Figure 5a and 5b we can conclude that the bottom layer shows slightly lower uncertainty than the top layer, but the difference is not pronounced. The lower uncertainty is due to the fact that the pumping wells are located in the bottom layer, where flow induced by the relatively well known abstraction rates dominates the flow induced by the relatively uncertain recharge. The above presented case illustrates that DCM is attractive because it is a conceptually simple method that can easily be implemented and used in combination with a standard groundwater flow model. The DCM iterations are terminated after having reached a minimum in difference between measured and calculated heads (or a minimum of the Darcy residual; see the Supporting Information). Such a minimum is reached after a few iterations. Looking for smaller minimums or the global minimum is not required in view of the observation inaccuracy obtained from the time series of estimated conductivities. Convergence to high accuracy is not required in realistic hydrogeological problems; it does not make sense to obtain an accuracy that is much higher than the accuracy introduced by the uncertain head and the generally even more uncertain flux estimates. Therefore, in our applications we used an initial calculation followed by two iterations. M.A. El-Rawy et al. Groundwater
5
(a)
(b)
Figure 5. Observation uncertainties (standard deviations) of hydraulic conductivities for case 1 with no soft data uncertainty: (a) top layer, (b) bottom layer.
The log-conductivities estimated by the KF range from 3.90 to 0.01 with an average conductivity of 15.18 m/d in the top layer, and from 2.65 to 0.002 with an average conductivity of 4.85 m/d in the bottom layer. Comparing these results with the results of the mean values of the log-conductivities of the 24 observations (Figure 4), we find that the mean value and the KF estimate are almost the same. In Figure 6 the initial conductivities, the timeaveraged DCM-estimated conductivities and the DCM-KF estimated hydraulic conductivities are shown at the locations of the piezometers in the top (27 piezometers; Figure 6a) and the bottom layer (34 piezometers; Figure 6b). The figure shows that in the measurement ranges, that is, the regions around the piezometers the conductivities determined by the DCM and by the DCM-KF do hardly differ from each other, but differ from the initial conductivities. The standard deviations of the estimation uncertainty range between 0.4 and 0.0 in the top layer (Figure 7a) and between 0.26 and 0.0 in the bottom layer (Figure 7b). As a result of the KF these standard deviations are about 77% smaller than the standard deviations of the observation uncertainty of the log-conductivities obtained by the DCM. Figure 6 also shows a comparison of the standard deviations of the log-conductivities obtained from the time-averaged DCM with the standard deviations obtained from the DCM-KF estimate at locations of piezometers. At piezometer 1-0081 and 1-0098 the hydraulic conductivity has much less uncertainty than at other piezometers, even before applying the KF. On the other hand, piezometer 1-0402 has the highest uncertainty in the log-conductivities. The average standard deviation for all piezometers in the top layer is 0.92 after time-averaged DCM and 0.22 after the KF. For the bottom layer these standard deviations are 0.71 and 0.16, respectively. It is clear that the standard deviation averaged over the piezometers in the top layer is higher than the standard deviations averaged over the piezometers in the bottom layer. This is consistent with the deeper flux argument presented above. 6
M.A. El-Rawy et al. Groundwater
In this case, we have assumed that there is no uncertainty in the soft data, that is, we have kept the initial conductivities the same for each DCM conductivity estimation. As a consequence, we found the greatest variances in the grid blocks near the piezometers. This is caused by the fact that conductivities in the “terra incognita” will hardly be updated by DCM. Hence, KF calculations make their variances very small or zero. Case 2: Soft Data with Uncertainty In practice, the uncertainty in the soft data should be significantly greater than the observation error found near the piezometers. A constructive way to avoid misinterpretation of the variances in the “terra incognita” is to determine realistic variances by specifying at the different observation times different possible distributions of the soft conductivity data. In this case the 24 applications of DCM are not only based on 24 different measured head-flux couples, but also on 24 different soft data sets. To demonstrate the principle we have chosen an extremely simple series of different initial conductivity patterns. For each grid block i the different initial patterns have a mean logconductivity pattern x i (i = 1, 2, . . . , 737,884) that is equal to the pattern used in case 1 without soft data uncertainty. One of the simplest ways to obtain soft data set m is xmi = x i (1 + c (m − (M + 1) /2)), where m = 1, 2, . . . , M is the distribution number and c is a constant. i . We observe that in our case, where M = 24, x i = xm=12.5 i i For c = 0.045 we obtain xm = x (1 + 0.045 (m − 12.5)). Admittedly, this is not a realistic random distribution, but we consider this simple approach as sufficient to show the effect of uncertainty in the initial conductivities, that is, the uncertainty in the soft data. Figure 8a and 8b presents the mean values of the log-conductivities calculated for the top and the bottom layer respectively. They range from 3.83 to 0.005 with an average conductivity of 14.73 m/d in the top layer and from 2.65 to 0.002 with an average conductivity of 4.76 m/d in the bottom layer. Comparison of these log-conductivities with those of the case with zero soft NGWA.org
(a)
(b)
Figure 6. Horizontal hydraulic conductivities at locations of piezometers in the top (a) and bottom (b) layer: initial conductivities K0 (blue); after the time-averaged DCM (red) with their standard deviations; and after the KF (black) with their standard deviations. The error bars represent standard deviation.
(a)
(b)
Figure 7. Estimation uncertainties (standard deviations) of hydraulic conductivities for case 1 with no soft data uncertainty: (a) top layer; (b) bottom layer.
data error (Figure 4a and 4b) shows that there are small differences. The observation uncertainties (Figure 9a and 9b) range from 1.76 to 0.004 and from 1.16 to 0.0 for the top and bottom layer respectively. We can conclude that, like in case 1, the top layer has more uncertainty than the bottom layer, due to the uncertainty in recharge, which in this case is added to the uncertainty in the soft data. NGWA.org
Comparison with Figure 5a and 5b shows that the red regions are the same for the case with no soft data uncertainty (c = 0) and the case with appreciable soft data uncertainty (c = 0.045). As a consequence, in these regions the variances represent a realistic estimation. In the regions further away from the red regions the presented observation uncertainties, blue in case of zero soft data uncertainty and yellow/green in case of M.A. El-Rawy et al. Groundwater
7
(a)
(b)
Figure 8. Mean values of log-conductivities for 24 observations, for case 2 with soft data uncertainty c= 0.045: (a) top layer, (b) bottom layer.
(a)
(b)
Figure 9. Observation uncertainties (standard deviation) of hydraulic conductivities for case 2 with soft data uncertainty c = 0.045: (a) top layer; (b) bottom layer.
appreciable soft data uncertainty, differ considerably. The regions with relatively large uncertainty difference between case 1 (low uncertainty) and case 2 (relatively high uncertainty) represent the “terra incognita,” that is, regions where, on basis of the available piezometers, the inverse model cannot improve the perceived hydraulic conductivities. To assign conductivities to these regions the hydrogeologist has to use additional knowledge, for instance knowledge about the subsurface’s facies distributions, or information about continuity of layers obtained from boreholes, etc., often combined with a technique to determine conductivity from grain size data obtained from samples (Rogiers et al. 2012). There are several geostatistical approaches that incorporate already available geological knowledge for solving the inverse problem in groundwater (Doherty 2003; Doherty et al. 2010). In this respect multiple-point geostatistical models based on geological analogs might be especially helpful (Huysmans et al. 2008; Huysmans and Dassargues 2009, 2011). The log-conductivities estimated by the KF for the 24 observations range from 3.89 to 0.004 with an average conductivity of 15.18 m/d for the top layer, and from 8
M.A. El-Rawy et al. Groundwater
2.65 to 0.002 with an average conductivity of 4.83 m/d for the bottom layer. If we compare these results with the results of the mean values of the log-conductivities of the 24 observations (Figure 8a and 8b), we do not observe big differences. In order to clarify this point, we have calculated the layer-averages of the horizontal hydraulic conductivities at the locations of the piezometers for case 1 and case 2. Table 1 shows clearly that in case 2 the layer-averages of the conductivities at the locations of the piezometers have smaller values than those averages in case 1, both for the time-averaged DCM and for the DCM-KF estimates. It is clear that in the “measurement ranges,” that is, at the locations of the piezometers, the time-averaged DCM conductivities, as well as the estimated DCM-KF conductivities have been changed with respect to the initial hydraulic conductivities. It is also shown that the differences between the time-averaged DCM conductivities and the estimated DCM-KF conductivities are relatively small. Comparing the hydraulic conductivities of case 2 (c = 0.045) with case 1 (c = 0), we can conclude that at the locations of the piezometers both the DCM and the DCM-KF conductivities in case 2 have a lower value NGWA.org
Table 1 Layer-Average of the Horizontal Hydraulic Conductivity at the Locations of the Piezometers for the Top and the Bottom Layer: Initial Conductivities, DCM-Estimated Conductivities and DCM-KF-Estimated Conductivities for Cases 1 and 2 Kh avg (m/d) Top Layer
Before DCM After DCM After KF
Bottom Layer
Case 1
Case 2
Case 1
Case 2
16.09 6.63 6.99
16.09 6.09 6.27
7.24 4.04 4.16
7.24 3.67 3.79
Table 2 Mean Absolute Error (MAE) and the Root Mean Square Error (RMSE) for Initial, DCM-Estimated and DCM-KF-Estimated Conductivities with Zero and Nonzero Soft Data Uncertainty Initial
DCM
DCM-KF
Zero Nonzero Zero Nonzero Zero Nonzero MAE 0.8 RMSE 1.13
0.8 1.13
0.6 0.9
0.74 1.06
0.6 0.89
0.73 1.05
with respect to the initial conductivities than those of case 1. The estimation uncertainties calculated by the KF, that is, the standard deviations of the estimates, range from 0.38 to 0.001 for the top layer (Figure 10a) and from 0.27 to 0.0 for the bottom layer (Figure 10b). The estimation uncertainty calculated by the KF is about 77% smaller than the observation uncertainty. Differences Between the Two Cases In order to calculate the differences between the observed and simulated heads we have used a forward MODFLOW model with flux boundary conditions and three different conductivity fields: (1) the initial grid block conductivities (the soft data); (2) the time-averaged conductivities obtained from DCM-calibration; (3) the conductivities estimated by the KF from the 24 DCMobservations. MODFLOW is run for the hypothetical case where the recharge flux is equal to the average of the 24 time-dependent recharge rates. The calculated heads are then compared with the time-averaged measured heads. The mean absolute error (MAE) and the root mean square error (RMSE) for the cases 1 and 2 are presented in Table 2. For the case 1 with zero soft data uncertainty we found that the MAE and RMSE obtained from the KF estimates are relatively small, 0.60 m and 0.89 m, respectively. These errors deviate considerably from the errors obtained from the model with initial grid block NGWA.org
conductivities (MAE = 0.8 m and RMSE = 1.13 m). This shows that the estimated conductivities do represent an improvement for this time-averaged case. This is not surprising as the KF estimate is a time-averaged value. In case 2 with non-zero uncertainty of the soft data (Table 2), the MAE and RMSE values before and after application of the KF are slightly higher than the values for the case with zero soft data uncertainty. This is expected, as uncertainties in the case with nonzero soft data uncertainty are higher and more realistic than uncertainties in the case with zero soft data uncertainty. Measurement Ranges and Terra Incognita Using both DCM hydraulic conductivity estimates resulting from the initial conductivity pattern of case 1 and the series of different initial conductivity patterns of case 2, it is possible to assess in which regions conductivities can be determined by the hydraulic data and in which regions this is impossible. To delineate these regions, we have determined for each grid block the relative difference between the standard deviations of case 1 (without soft data uncertainty) and the standard deviations of case 2 (with soft data uncertainty) (Figure 11a and 11b). In the “measurement ranges” near the piezometers the observed conductivities and their uncertainties (both observation and estimation uncertainties) are approximately the same for cases 1 and 2 (the grey regions in Figure 11). For these regions, in which the standard deviations are mainly determined by the uncertain head-flux measurements, we may assume that the results represent a realistic conductivity estimate. However, in the regions further away from the piezometers (from the orange to blue regions in Figure 11) the standard deviations differ considerably between case 1 and case 2, representing the “terra incognita.” In these regions additional hydraulic measurements (piezometers) are needed to find a more reliable conductivity pattern.
Conclusions 1. The DCM consists of alternating flux and head specified (MODFLOW) runs using observations of head-flux couples; the difference between the runs provides insight in areas where the hydraulic conductivity can and cannot be adjusted from its initial value. The DCM is based on two (MODFLOW) model runs based on the same initially chosen grid block conductivities: •
The flux run has as many specified flux data as available (abstraction well rates, indirectly “measured” recharge rates derived from measured rainfall, runoff and evapotranspiration). • The head run has as many specified head data as available (including the heads measured in piezometers), like in T´oth’s (2009) flow systems approach, where instead M.A. El-Rawy et al. Groundwater
9
(a)
(b)
Figure 10. Estimation uncertainties (standard deviation) σ of hydraulic conductivities for case 2 with soft data uncertainty c = 0.045: (a) top layer, (b) bottom layer.
(a)
(b)
Figure 11. Relative differences between the standard deviations of case 1 and case 2: (a) top layer, (b) bottom layer.
of the recharge the head is specified on the phreatic surface. The differences in solution of the flux and the head model give insight in the so-called “terra incognita,” that is, the regions where the “measured” head-flux couples cannot improve the initial hydraulic conductivity, and in the so-called “measurement ranges,” that is, the regions where the “measured” head-flux couples do adjust the initial conductivity. In addition, the differences in the solution of the two runs can be used to update the grid block conductivities by application of Darcy’s law. As a consequence, the DCM is essentially based on the principle that calibration is meaningful only if we can specify more head and flux data than required for a forward problem. That is, we need to know head-flux couples in a number of wells or at a number of boundary points. These hydraulic data have to be acquired from sources independent from the model to be calibrated. 2. DCM is cheap and can support trial and error or PEST calibration. 10
M.A. El-Rawy et al. Groundwater
The DCM is relatively cheap in computation time and memory requirements. Also if we apply the “trialand-error” or a more comprehensive indirect calibration method like PEST (Doherty and Hunt 2010b), running both a flux model and a head model may enhance insight in the calibration problem and can assist in delineating a meaningful zonation. In addition, DCMestimated conductivities in the measurement ranges could be used as pilot points (Certes and de Marsily 1991), for instance in a subsequent PEST calibration (Hunt et al. 2007). In this sense the DCM and its underlying idea of a “flux run” and a “head run” may be considered as a calibration assisting tool. 3. Repeated DCM estimations with KF postprocessing, results in conductivities with lower uncertainty. The calibration assisting methodology can be enhanced by repeated DCM estimations under different hydrological conditions. For instance, estimation in the winter season with relatively high rainfall and high water levels, and estimation in the summer season with NGWA.org
relatively little rainfall and low water levels. If the “headflux” measurements were without measurement error, we would find the same conductivities for all hydrological conditions. In reality we find a series of time-dependent conductivities, from which we can determine the time-averaged mean values and standard deviation or observation error. Using this observation error as input of a KF, an estimation of the time-independent grid block conductivities can be obtained from the time series of observed (DCM calibrated) conductivities. The KF then determines for the estimated conductivities an uncertainty that is appreciably smaller than the observation error. 4. DCM-KF is a valuable calibration tool and enhances hydrological understanding. In conclusion we may say that the calibration strategy underlying the DCM and subsequent application of the KF works well while enhancing hydrological understanding. Therefore this approach to calibration may be considered as a valuable tool to assist hydrogeologists in their calibration of the hydraulic conductivities.
Acknowledgments The authors gratefully acknowledge Randall J. Hunt and two anonymous reviewers for their constructive review, which considerably improved this paper. Also James McCallum is acknowledged for suggestions for improvement. The first author acknowledges Erasmus Mundus External Cooperation Window and Vrije Universiteit Brussel for the funding of his scholarship.
Supporting Information Additional Supporting Information may be found in the online version of this article: Figure S1. Flow through one-dimensional domain partitioned in sections with lengths L1 and L2 , conductivities K 1 and K 2 . Table S1. Overview of the Hydrostratigraphy of the Study Area
References Batelaan, O., and F. De Smedt. 2004. SEEPAGE, a new MODFLOW DRAIN package. Ground Water 42, no. 4: 576–588. Borcea, L. 2002. Electrical impedance tomography. Inverse Problems 18, no. 6: R99–R136. Brouwer, G.K., P.A. Fokker, F. Wilschut, and W. Zijl. 2008. A direct inverse model to determine permeability fields from pressure and flow rate measurement. Mathematical Geosciences 40, no. 8: 907–920. Calder´on, A.P. 1980. On an Inverse Boundary Value Problem, Seminar on Numerical Analysis and its Applications to Continuum Physics (Rio de Janeiro, 1980), 65–73. Rio de Janeiro: Soc. Brasil Mat. Cao, W., W.B. Bowden, T. Davie, and A. Fenemor. 2006. Multi-variable and multi-site calibration and validation of
NGWA.org
SWAT in a large mountainous catchment with high spatial variability. Hydrological Processes 20, no. 5: 1057–1073. Carrera, J. 1988. State of the art of the inverse problem applied to the flow and solute transport equations. In Groundwater Flow and Quality Modelling, ed. E. Custodio, A. Guirgui, and J.P.L. Ferreira, 549–583. Netherlands: D. Reidel Publishing Company/Springer. Carrera, J., A. Alcolea, A. Medina, J. Hidalgo, and L.J. Slooten. 2005. Inverse problem in hydrogeology. Hydrogeology Journal 13, no. 1: 206–222. Certes, C., and G. de Marsily. 1991. Application of the pilot points method to the identification of aquifer transmissivities. Advances in Water Resources 14, no. 5: 285–300. Cools, J., Y. Meyus, S.T. Woldeamlak, O. Batelaan, and F. De Smedt. 2006. Large-scale GIS-based hydrogeological modelling of Flanders: A tool for groundwater management. Environmental Geology 50, no. 8: 1201–1209. Dams, J., E. Salvadore, T. Van Daele, V. Ntegeka, P. Willems, and O. Batelaan. 2012. Spatio-temporal impact of climate change on the groundwater system. Hydrology Earth System Sciences 16: 1517–1531. DOI:10.5194/hess-16-1517-2012. Dams, J., S.T. Woldeamlak, and O. Batelaan. 2008. Predicting land-use change and its impact on the groundwater system of the Kleine Nete catchment, Belgium. Hydrology Earth System Sciences 12: 1369–1385. DOI:10.5194/hess-121369-2008. de Marsily, G., J.P. Delhomme, A. Coudrain-Ribstein, and A.M. Lavenue. 2000. Four decades of inverse problems in hydrogeology. Geological Society of America. Special Paper 348. Doherty, J. 2003. Ground water model calibration using pilot points and regularization. Ground Water 41, no. 2: 170–177. Doherty, J., and R.J. Hunt. 2010a. Response to comment on: Two statistics for evaluating parameter identifiability and error reduction. Journal of Hydrology 380: 489–496. DOI:10.1016/j.jhydrol.2009.10.012. Doherty, J.E., and R.J. Hunt. 2010b. Approaches to highly parameterized inversion: A guide to using PEST for groundwater-model calibration. U.S. Geological Survey Scientific Investigations Report 2010-5169, 60. http://pubs.usgs.gov/sir/2010/5169/. Doherty, J.E., R.J. Hunt, and M.J. Tonkin. 2010. Approaches to highly parameterized inversion: A guide to using PEST for model-parameter and predictive-uncertainty analysis. U.S. Geological Survey Scientific Investigations Report 20105211, 71. Doherty, J., and R.J. Hunt. 2009. Two statistics for evaluating parameter identifiability and error reduction. Journal of Hydrology 366: 119–127. DOI:10.1016/j.jhydrol.2008.12.018. Durlofsky, L.J. 1991. Numerical calculation of equivalent grid block permeability tensors for heterogeneous porous media. Water Resources Research 27: 699–708. El-Rawy, M.A. 2013. Calibration of hydraulic conductivities in groundwater flow models using the double constraint method and the Kalman filter. PhD thesis, Vrije Universiteit Brussel, Brussels, Belgium. El-Rawy, M.A., W. Zijl, O. Batelaan, and G.A. Mohammed. 2010. Application of the double constraint method combined with MODFLOW. In Proceedings of the Valencia IAHR International Groundwater Symposium, September 22–24. Emsellem, Y., and G. de Marsily. 1971. An automatic solution for the inverse problem. Water Resources Research 7, no. 5: 1264–1283. Ewald, C.-O. 2006. The Malliavin Gradient Method for the calibration of stochastic dynamical models. Applied Mathematics and Computer 175, no. 2: 1332–1352. Fasanino, G., J.E. Molinard, G. de Marsily, and V. Pelc´e. 1986. Inverse modeling in gas reservoirs. SPE 15592, 61st SPE
M.A. El-Rawy et al. Groundwater
11
Annual Technical Conference and Exhibition, October 5–8, New Orleans, Louisiana. Hager, W.W., and H. Zhang. 2006. A survey of nonlinear conjugate gradient methods. Pacific Journal of Optimization 2, no. 1: 35–58. Haitjema, H. 2006. The role of hand calculations in ground water flow modeling. Ground Water 44, no. 6: 786–791. DOI:10.1111/j.1745-6584.2006.00189.x. Harbaugh, A.W., E.R. Banta, M.C. Hill, and M.G. McDonald. 2000. MODFLOW-2000. The U.S. Geological Survey modular ground-water model—User guide to modularization concepts and the ground-water flow process. U.S. Geological Survey Open-File Report 00-92, 121. Reston, Virginia: USGS. Hendricks Franssen, H.J., A. Alcolea, M. Riva, M. Bakr, N. Van der Wiel, F. Stauffer, and A. Guadagnini. 2009. A comparison of seven methods for the inverse modelling of groundwater flow. Application to the characterisation of well catchments. Advances in Water Resources 32, no. 6: 851–872. Hill, M.C., and C.R. Tideman. 2007. Effective Groundwater Model Calibration: With Analysis of Data, Sensitivities, Predictions, and Uncertainty. Hoboken, New Jersey: John Wiley and Sons, Inc. Hunt, R.J., J. Doherty, and M.J. Tonkin. 2007. Are models too simple? Arguments for increased parameterization. Ground Water 45, no. 3: 254–261. DOI:10.1111/j.1745-6584.2007. 00316.x. Huysmans, M., and A. Dassargues. 2011. Direct multiple-point geostatistical simulation of edge properties for modeling thin irregularly shaped surfaces. Mathematics Geosciences 43, no. 5: 521–536. Huysmans, M., and A. Dassargues. 2009. Application of multiple-point geostatistics on modelling groundwater flow and transport in a cross-bedded aquifer (Belgium). Hydrogeology Journal 17, no. 8: 1901–1911. Huysmans, M., L. Peeters, G. Moermans, and A. Dassargues. 2008. Relating small-scale sedimentary structures and permeability in a cross-bedded aquifer. Journal of Hydrology 361, no. 1–2: 41–51. Keidser, A., and D. Rosbjerg. 1991. A comparison of four inverse approaches to groundwater flow and transport parameter Identification. Water Resources Research 27, no. 9: 2219–2232. Kitanidis, P.K. 1997. The minimum structure solution to the inverse problem. Water Resources Research 33, no. 10: 2263–2272. Kohn, R.V., and A. McKenney. 1990. Numerical implementation of a variational method for electrical impedance tomography. Inverse Problems 6, no. 3: 389–414. Kohn, R.V., and M. Vogelius. 1987. Relaxation of a variational method for impedance computed tomography. Communications on Pure Applied Mathematics 40, no. 6: 745–777. Kuiper, L.K. 1986. A comparison of several methods for the solution of the inverse problem in two-dimensional steady-state groundwater flow modeling. Water Resources Research 22, no. 5: 705–714. LaVenue, A.M., B.S. RamaRao, G. de Marsily, and M.G. Marietta. 1995. Pilot point methodology for automated calibration of an ensemble of conditionally simulated transmissivity fields, 2, Application. Water Resources Research 31, no. 3: 495–516. Liu, Y.B., S. Gebremeskel, F. De Smedt, L. Hoffmann, and L. Pfister. 2003. A diffusive transport approach for flow routing in GIS-based flood modeling. Journal of Hydrology 283, no. 1–4: 91–106. McLaughlin, D., and L.R. Townley. 1996. A reassessment of the groundwater inverse problem. Water Resources Research 32, no. 5: 1131–1161.
12
M.A. El-Rawy et al. Groundwater
Nelson, R.W. 1968. In-place determination of permeability distribution for heterogeneous porous media through analysis of energy dissipation. Society of Petroleum Engineers Journal 8, no. 1: 33–42. http://www.onepetro.org/mslib/ servlet/onepetropreview?id=00001554. Nelson, R.W. 1962. Conditions for determining areal permeability distribution by calculation. Society of Petroleum Engineers Journal 2, no. 3: 223–224. http://www.onepetro.org/ mslib/servlet/onepetropreview?id=00000371. Neuman, S.P. 1973. Calibration of distributed parameter groundwater flow models viewed as a multiple-objective decision process under uncertainty. Water Resources Research 9, no. 4: 1006–1021. Nilsson, B., A.L. Højberg, J.C. Refsgaard, and L. Troldborg. 2007. Uncertainty in geological and hydrogeological data. Hydrology Earth System Sciences 11: 1551–1561. Pinault, J.L., and S. Schomburgk. 2006. Inverse modeling for characterizing surface water/groundwater exchanges. Water Resources Research 42, no. 8. DOI:10.1029/2005WR 004587. Poeter, E.P., and M.C. Hill. 1997. Inverse models-A necessary next step in ground-water modeling. Ground Water 35, no. 2: 250–260. Rajanayaka, C., and D. Kulasiri. 2001. Investigation of a parameter estimation method for contaminant transport in aquifers. Journal of Hydroinformatics 3: 203–213. RamaRao, B.S., A.M. LaVenue, G. de Marsily, and M.G. Marietta. 1995. Pilot point methodology for automated calibration of an ensemble of conditionally simulated transmissivity fields: 1. Theory and computational experiments. Water Resources Research 31, no. 3: 475–493. Rogiers, B., D. Mallants, O. Batelaan, M. Gedeon, M. Huysmans, and A. Dassargues. 2012. Estimation of hydraulic conductivity and its uncertainty from grain-size data using GLUE and artificial neural networks. Mathematical Geosciences 44, no. 6: 739–763. Sagar, B., S. Yakowitz, and L. Duckstein. 1975. A direct method for the identification of the parameters of dynamic nonhomogeneous aquifers. Water Resources Research 11, no. 4: 563–570. Sun, N.-Z. 2004. Inverse Problems in Groundwater Modeling. Dordrecht: Kluwer Academic Publishers. Tikhonov, A.N. 1963a. Solution of incorrectly formulated problems and the regularization method. Soviet Mathematics Doklady 4: 1035–1038. Tikhonov, A.N. 1963b. Regularization of incorrectly posed problems. Soviet Mathematics Doklady 4: 1624–1637. T´oth, J. 2009. Gravitational Systems of Groundwater Flow . Cambridge: Cambridge University Press. Trykozko, A., G.A. Mohammed, and W. Zijl. 2009. Downscaling: The inverse of upscaling. In Conference on Mathematical and Computational Issues in the Geosciences. SIAM GS 2009 , June 15–18, Leipzig. Trykozko, A., G.K. Brouwer, and W. Zijl. 2008. Downscaling: A complement to homogenization. International Journal of Numerical Analysis and Modeling 5, no. supp: 157–170. Warren, J.E., and H.S. Price. 1961. Flow in heterogeneous porous media. Society of Petroleum Engineers Journal 1, no. 3: 153–169. DOI:10.2118/1579-G. Wen, X.-H., C.V. Deutsch, J.J. Gomez-Hernandez, and A.S. Cullick. 1999. A program to create permeability fields that honor single-phase flow rate and pressure data. Computers and Geosciences 25, no. 3: 217–230. Wen, X.-H., C.V. Deutsch, and A.S. Cullick. 1998. Highresolution reservoir models integrating multiple-well production data. Society of Petroleum Engineers Journal 3, no. 4: 344–355. Wexler, A. 1988. Electrical impedivity imaging in two and three dimensions. Clinical Physics and Physiological Measurement 9, no. Suppl A: 29–33.
NGWA.org
Wexler, A., B. Fry, and M.R. Neuman. 1985. Impedivitycomputed tomography algorithm and system. Applied Optics 24, no. 23: 3985–3992. Yeh, W.W.-G. 1986. Review of parameter identification procedures in groundwater hydrology: The inverse problem. Water Resources Research 22, no. 2: 95–108. Yorkey, T.J., and J.G. Webster. 1987. A comparison of impedivity tomographic reconstruction algorithms. Clinical Physics and Physiological Measurement 8, no. Suppl A: 55–62. Zijl, W., G.A. Mohammed, O. Batelaan, and F. De Smedt. 2010. Constraining methods for direct inverse modeling. XVIII International Conference on Water Resources, CMWR,
NGWA.org
CIMNE, Barcelona. http://congress.cimne.com/cmwr2010/ Proceedings/docs/p62.pdf. Zijl, W., and A. Trykozko. 2001. Numerical homogenization of the absolute permeability using the conformal-nodal and mixed-hybrid finite element method. Transport in Porous Media 44: 33–62. Zimmerman, D.A., G.D. de Marsily, C.A. Gotway, M.G. Marietta, C.L. Axness, R.I. Beauheim, R.I. Bras, J. Carrera, G. Dagan, and P.B. Davies. 1998. A comparison of seven geostatistically based inverse approaches to estimate transmissivities for modeling advective transport by groundwater flow. Water Resources Research 34, no. 6: 1373–1413.
M.A. El-Rawy et al. Groundwater
13
