- Email: [email protected]
The importance of observations on fluxes to constrain ground water model calibration
Available online at www.sciencedirect.com
Physics and Chemistry of the Earth 33 (2008) 1105–1110 www.elsevier.com/locate/pce
The importance of observations on fluxes to constrain ground water model calibration Chiara Vassena *, Cinzia Durante, Mauro Giudici, Giansilvio Ponzini Universita´ degli Studi di Milano, Dip. di Scienze della Terra, Milano, Italy Received 8 October 2007; received in revised form 18 January 2008; accepted 28 January 2008 Available online 5 February 2008
Abstract The aquifer system in the alluvial basin bordered by Adda, Po and Oglio rivers (Northern Italy) is characterised by a dual flow regime. In shallow sediments, which constitute a phreatic aquifer with high conductivity, great fluxes are driven by the interaction between ground water and the network of surface water, by the infiltration of rain and irrigation water, and by the fluxes drained from depression springs and river valley terraces. The underlying semiconfined aquifers are characterised by minor fluxes driven by water abstraction from wells of the public Water Works. Since most of the ground water flow occurs in the phreatic aquifer, an equivalent single layer 2D steady state flow model has been calibrated. The identification of the transmissivity field at the scale of the model has been obtained by solving an inverse problem with the comparison model method which requires an initial configuration, i.e., reference head, initial transmissivity field, source terms. Most of the head and source data are related to the phreatic aquifer, but most of the estimates of transmissivity are obtained with field tests conducted in deep wells pumping from the semiconfined aquifers, so that this kind of prior information cannot be used directly for model calibration. The inverse problem is underdetermined and a unique solution is not available. Furthermore information on surface hydrology is poor. Therefore many tests with different hypotheses about the initial configuration have been performed and some of them have been selected and used to initialise the automatic inversion procedure. Ó 2008 Elsevier Ltd. All rights reserved. Keywords: Mathematical modeling; Underdetermined inverse problem; Model calibration; Interaction between surface and ground water
1. Introduction The highly developed agricultural system of the Province of Cremona (Northern Italy) needs a great amount of water for irrigation which is presently derived from Adda, Po, Oglio and Serio rivers (Fig. 1). Limitations imposed by regional regulators and interannual climatic variations induce several farmers to ask for the permission to drill new water wells in the phreatic aquifer. Therefore the public administration has decided to develop a tool to support decisions and sustainable management of water resources. A research team of the University of Milan, merging the expertise of the Dipartimento di Scienze della Terra and of *
Corresponding author. E-mail address: [email protected] (C. Vassena).
1474-7065/$ - see front matter Ó 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.pce.2008.01.004
the Istituto di Idraulica Agraria (IIA), has developed a ground water flow model in the alluvial basin bordered by Adda, Po and Oglio rivers. The focus of this paper is on the model calibration. In fact, field tests provide values of the physical parameters that are relevant at scales different from those at which regional groundwater flow is modeled, so that the model parameters have to be computed as the solution of an inverse problem. Several review papers (e.g., Yeh, 1986; Carrera, 1988; Giudici, 2001; Carrera et al., 2005) provide comprehensive lists of the inverse methods proposed in groundwater hydrology and discuss the issues related to their application. For this case study the groundwater flow model has been calibrated with the comparison model method (CMM). This method was originally proposed by Scarascia and Ponzini (1972), successively developed by Ponzini and Lozej
1106
C. Vassena et al. / Physics and Chemistry of the Earth 33 (2008) 1105–1110
Fig. 1. Study area and contour plots of reference (solid line) and modelled (dashed line) piezometric head fields (years 2001–2003). Contour intervals: 5 m.
(1982), Ponzini and Crosta (1988), Ponzini et al. (1989), and applied to real aquifers (Associazione Irrigazione Est Sesia, 1979; Beatrizotti et al., 1983; Pasquier and Marcotte, 2005; Benoit et al., 2005). Recent advances about the application of the CMM to 3D transient flow have been proposed by Pasquier and Marcotte (2006) and Ponzini et al. (2007). The CMM requires the reference head field, the source terms and an initial tentative transmissivity field, which fix the water fluxes in the aquifer system. Data are of paramount importance for model development, especially for the model calibration and validation. Actually, data are often affected by errors, because it is difficult to perform experiments under controlled conditions in the field; furthermore, the modeler is forced to use already available data, often irregularly spaced and lacking in a check of the measurement procedures and accuracy. These remarks imply that the model predictions are not expected to be more accurate than the data used to develop and calibrate the model itself (Giudici, 2001). In this case study, data have been collected for purposes different from the development of a mathematical model so that some information is missing or uncertain. For instance, estimates of transmissivity are available from field tests performed in underlying semiconfined aquifers only, so that the initial transmissivity of the phreatic aquifer cannot be inferred from these data. The source terms, which include recharge from rain and irrigation, water extraction for drinking purposes and for agricultural and industrial use, drainage from depression springs and morphological terraces, recharge from artificial channels, and river-aquifer exchanges, are also estimated with uncertainty. In other words this case study is an example of an hydrological basin for which field data are not optimal, with reference
both to the space and time distribution and to the measurement quality. Therefore a sequence of tests coupling the CMM with different hypotheses about the source terms has been conducted. The choice of the best configuration (namely, parameters for the source terms) cannot be based on head data only, but other hydrological data are necessary, e.g. river discharge, spring fluxes, etc. In the case study this additional information is poor; as a consequence few configurations have been selected satisfying weak constraints and have been used to initialise the automatic inversion procedure.
2. The hydrogeological framework The study area (about 2500 km2) is located in the alluvial Po plain and is bordered by Adda, Po and Oglio rivers (Fig. 1); a three-year-long interval was considered for this study, from 2001 to 2003, because for this period the most complete data set was available. Current flow conditions of the aquifer system are still similar, even if the extraordinary repetition of drought periods could induce variation in the next years. The aquifer system is characterised by a dual flow regime. Shallow gravel and sandy-gravel sediments constitute a phreatic aquifer with high conductivity and an average thickness of about 50 m; the interaction between ground water and the network of surface water, the infiltration of rain and irrigation water, and the drainage from depression springs and river valley terraces induce large fluxes in this phreatic aquifer. The underlying semiconfined aquifers are discontinuous and characterised by minor
C. Vassena et al. / Physics and Chemistry of the Earth 33 (2008) 1105–1110
fluxes driven by leakage and by water extraction from wells of the public Water Works. The aquifer is fed by the inflow from north and by the recharge terms, subdivided into the recharge due to infiltration of the rain and of the water used for irrigation, F ðrechargeÞ , and the losses from irrigation channels, F ðchannelsÞ . Both F ðrechargeÞ and F ðchannelsÞ have been estimated by IIA with a model that computes a water mass balance at the ground surface and in the unsaturated zone, taking into account soil characteristics, the irrigation system, meteorological data, land use, crop production and growth (Facchi et al., 2004). The annual average estimates amount to 34 m3/s for F ðrechargeÞ and 29 m3/s for F ðchannelsÞ . Adda, Oglio and Po rivers are mainly acting as drains for the phreatic aquifer; Serio river is also draining the aquifer along most of its valley ðF ðriverÞ Þ. Drainage fluxes from depression springs, F ðspringsÞ , and from geomorphological terraces bounding the fluvial valleys, F ðterracesÞ , are sink terms for the mathematical model of the phreatic aquifer. They have to be added to the estimates of water abstraction fluxes, F ðabsÞ , for agricultural and industrial uses (about 11 m3/s from the phreatic aquifer) and for drinking water (about 3 m3/s from the semiconfined aquifer). We do not have reliable estimates of river-aquifer exchange fluxes. Ground water abstraction for irrigation purposes is estimated from IIA, whereas ground water uptake for industrial use is estimated from the database of the authorized water wells available from the local administrations. Abstraction for drinking water is known by measurements of the public Water Works. An interpolated head field hðintÞ has been obtained by kriging the piezometric head data (the time-averaged values measured in about 150 shallow wells from 2001 to 2003), the water levels of depression springs estimated at about 300 springs (values observed during a thorough survey in 2000), and the river levels estimated from about 300 transversal topographic sections measured along the main rivers (field observations support the assumption that these rivers are in direct hydraulic contact with the phreatic aquifer). A number of estimates of transmissivity are available from field tests in semiconfined aquifers, but we cannot profit from these data to infer the initial transmissivity of the phreatic aquifer, for which only four pumping tests are available. 3. The groundwater flow model The aquifer system of the Province of Cremona is modelled as an equivalent single-layer aquifer with 2D pseudostationary flow. The modelled piezometric head, h, is the solution of the forward problem:
r ðT rhÞ ¼ F h ¼ hðbcÞ
on X on @X
ð1Þ
where X is the aquifer domain and oX its border, T is the transmissivity, hðbcÞ is the prescribed head at the border, F is the source term.
1107
Dirichlet boundary conditions have been assigned, assuming hydraulic direct contact between Adda, Po, Oglio rivers and the aquifer. The source term, F, includes the following components. Interaction between the aquifer and the Serio river, F ðriverÞ , is modelled as ( h > bðriverÞ C ðriverÞ ðh hðriverÞ Þ ðriverÞ ¼ F C ðriverÞ ðhðriverÞ bðriverÞ þ lðriverÞ Þ h < bðriverÞ ð2Þ where h and hðriverÞ are the heads, respectively, of the aquifer and of the river; bðriverÞ and lðriverÞ are, respectively, the upper height (above m.s.l.) and the thickness of the river bed, C ðriverÞ is the conductance of the river bed. The drainage from depression springs, F ðspringsÞ , and the drainage from the geomorphological terraces, F ðterracesÞ , are modelled as F ðspringsÞ ¼ maxð0; C ðspringsÞ ðh bðspringsÞ ÞÞ F
ðterracesÞ
¼ maxð0; C
ðterracesÞ
ðh b
ðterracesÞ
ÞÞ
ð3Þ ð4Þ
where bðspringsÞ and bðterracesÞ are the drain base levels of, respectively, depression springs and geomorphological terraces and are estimated from field surveys (D’Auria and Zavagno, 2005) and from topographic maps; C ðspringsÞ and C ðterracesÞ are conductance coefficients of, respectively, springs at topographical depressions and terraces. The values of C ðriverÞ ; C ðspringsÞ ; C ðterracesÞ are assumed constant in space. Eq. (1) is solved with a conservative finite difference scheme applied to a square grid with spacing of 500 m, chosen on the basis of the data density and the goals of the model. The approximate solution is computed with an error on the global balance as small as 0.01% of the total inflow rate. 4. The model calibration The CMM, described in detail in Appendix A, is based on the forward solution of an auxiliary model, the comparison model, CM, which shares the same geometry, the same state equation, the same source terms, the same initial and boundary conditions, the same discretization as the model to be calibrated, but its transmissivity field is a tentative seed field. The transmissivity field is identified adapting iteratively the seed field with a procedure based on Darcy’s law and on the guess that the flow rate estimated from the solution to the CM is a good approximation of the real one. See Appendix A for a detailed explanation. The inverse problem is underdetermined and in principle admits infinite solutions. Furthermore in this case study no direct information on the values of C ðriverÞ ; C ðspringsÞ ; C ðterracesÞ is available and other source terms, for example recharge terms and losses from the irrigation channels, are estimated with some uncertainties. In order to explore a lot of possible solutions, several tests have been performed changing the values of C ðriverÞ ; C ðspringsÞ ; C ðterracesÞ and the source terms.
C. Vassena et al. / Physics and Chemistry of the Earth 33 (2008) 1105–1110
4.1. Search of the optimal source terms for the application of the CMM The source terms have been changed according to the following alternatives: values of C ðspringsÞ between 0.025 and 0.1 m2/s, values of C ðterracesÞ between 0.00625 and 0.05 m2/s and C ðriverÞ between 0.025 and 0.1 m2/s. These intervals were fixed after numerical tests and processing of field data. Values of the recharge terms and losses from channels correspondent to the 80% and to the 50% of the values estimated by IIA with a SVAT model have been also tested. For each test hðiniÞ is computed as the solution to (1) with a seed transmissivity field computed as described in section A.2, hðbcÞ ¼ hðintÞ and F ¼ F ðiniÞ , where F ðiniÞ ¼ F ðabsÞ F ðrechargeÞ þ F ðspringsÞ F ðchannelsÞ þ F ðterracesÞ þ F ðriverÞ
ð5Þ
Notice that F ðterracesÞ , F ðriverÞ and F ðspringsÞ depend on hðiniÞ . More than 2000 tests, fully mixing the alternatives described before, have been conducted and their performance has been evaluated with a simple statistics, MAE, the spatial average of the head residuals j hðiniÞ hðintÞ j. For more than 200 tests MAE < 1 m, which is a value consistent with the error on the head measurements and on the subsequent interpolation and regularisation; hence further criteria are needed to choose the best distribution of source terms. In order to reduce the uncertainty of the inverse solution due to the fact that this inverse problem is underdetermined, the following constraints, related to water fluxes, have been fixed: (1) total drainage from Serio river <7 m3/s, based on the observation that the average of Serio river discharge is about 20 m3/s and water uptake for irrigation in the area is about 5 m3/s; (2) total drainage from springs >10 m3/s, based on the observation that the average annual discharge measured in channels fed by springs is about 20 m3/s; however these channels are partly fed also by the irrigation network. Recall that all these values are time averages from 2001 and 2003, since the model is calibrated for the average annual flow. After the application of these constraints, 42 tests have been selected, and in order to extract the most promising ones, some further weak constraints have been added to the previous criteria:
According to these criteria five tests have been selected: for these tests MAE varies between 0.83 and 0.91 m. All tests share the same values of C ðspringsÞ and C ðriverÞ , respectively, 0.05 and 0.025 m2/s. Instead C ðterracesÞ ¼ 0:025 m2 =s for Tests 1 and 2 and C ðterracesÞ ¼ 0:05m2 =s for the remaining tests. For all tests the other terms of the balance, i.e., the exchanges through boundaries, the drainage terms and the recharge terms, differ from the values averaged among the five tests by less than 20%. 4.2. Results As an example of calibration results, the transmissivity field obtained using Test 1 as the initial configuration for the CMM is shown. After four iterations of the CMM, MAE < 0.45 m. The identified transmissivity field permits a good fit between observed and simulated heads (Fig. 2), shows the large scale heterogeneity of the aquifer system (Fig. 3) and permits a preliminary evaluation of the groundwater balance. A consistency check has been performed along Serio river. In fact a reference head field, independent from ðrefÞ river-aquifer interaction, hnoSerio , has been obtained by interpolating and filtering (see Section A.1) the data described in Section 2 excluding the hydrodynamic levels along the Serio river. Also F ðnoSerioÞ has been obtained filtering the difference F ðiniÞ F ðriverÞ , and a transmissivity field, T ðnoSerioÞ , has been identified. The balance error computed as BE ¼ r ðT ðnoSerioÞ rhðrefÞ Þ F ðnoSerioÞ is closed to F ðriverÞ and the similarity between the T field plotted in Fig. 3 and T ðnoSerioÞ can be considered as a validation of the calibration results and of the value of C ðriverÞ .
115 105
Modelled piezometric head (m)
1108
95 85 75 65 55 45 35
(1) more than 50% of the total number of springs should be active; this prevents from modeling few active springs with unrealistically high discharge rates; R (2) X ðF ðspringsÞ þ F ðterracesÞ þ F ðriverÞ Þ dX > 20 m3 =s, according to the fact that the average annual discharge measured in channels fed by springs and morphological terraces is about 20 m3/s.
25 15 15
25
35
45
55
65
75
85
95
105
115
Observed piezometric head (m) Fig. 2. Comparison between observed and modelled piezometric heads. For most of the points the residual error is less than ±2 m (black lines).
C. Vassena et al. / Physics and Chemistry of the Earth 33 (2008) 1105–1110
1109
Fig. 3. Identified LogT field (m2/s) for Test 1 (C ðriverÞ ¼ 0:025 m2 =s, C ðspringsÞ ¼ 0:05 m2 =s, C ðterracesÞ ¼ 0:025 m2 =s).
5. Conclusions Calibration of ground water flow models often does not provide a unique solution; this is the case for this aquifer, where a great number of configurations of source terms, and therefore of T fields obtained with the inverse procedure, are found. Moreover the results of the CMM depend also on the choice of the initial seed T field, T ð0Þ , and of the reference head field; the procedure to obtain T ð0Þ proposed in Appendix A proved to be the most robust among different alternatives and therefore can be applied with confidence to study poorly gauged hydrological basins. The lack of accurate data on surface hydrology prevented us from assigning strict constraints on the inversion results; nevertheless, constraints on exchange fluxes between surface and ground water were useful to limit the uncertainty. This shows, on one hand, the importance of the efforts aimed to monitor hydrological basins with the acquisition of field or remote, hard or soft data, and, on the other hand, the utility of inverse methods to check the consistency and possible deficiency of a data set. Moreover the application of inverse methods like the CMM to poorly gauged basins claims developments to handle the lack of information about the source terms in the inversion. As a general final remark, this application shows that merging surface and groundwater study is important not only for an integrated water balance at the catchment scale, but also for the characterization of the aquifer system. Acknowledgements This work was financially supported by the Province of Cremona. Dr. M. Pesaro and Dr. M.M. Cremonini are
kindly acknowledged for the continuous support and for releasing the permission to publish this paper. Prof. C. Gandolfi and his coworkers at the Istituto di Idraulica Agraria of the University of Milan are kindly acknowledged for the cooperation during the work, in particular for the modeling of the irrigation system and the estimate of the aquifer recharge. Appendix A. The comparison model method A.1. Basic properties of the method The CMM is explicitely based on Darcy’s law and on the forward solution of an auxiliary model, the comparison model, CM, which gives the solution hð0Þ , to the forward problem (1) for T ¼ T ð0Þ , F ¼ F ðrefÞ and hðbcÞ ¼ hðrefÞ . The reference head field hðrefÞ is obtained from the application of an averaging procedure (high-wavenumber-cut filter) to the piezometric head interpolated from the observations hðintÞ ; notice that hðrefÞ ¼ hðintÞ on oX. The reference source term F ðrefÞ is obtained from the application of a similar filter to F ðiniÞ . This filtering procedure has been introduced to limit the instability effects of high-wavenumber components of h on inversion (Giudici and Vassena, 2007) and is of paramount importance to obtain reliable results. Then rhð0Þ is compared with the observed hydraulic gradient rhðrefÞ , the basic assumption is that the tentative seed transmissivity field generates an hð0Þ field such that the auxiliary flow rate per unit length ~ Qð0Þ ¼ T ð0Þ rhð0Þ ~ ¼ T rhðrefÞ is a good approximation of the real one, Q where T is the unknown transmissivity field to be estimated.
1110
C. Vassena et al. / Physics and Chemistry of the Earth 33 (2008) 1105–1110
~ and Q ~ð0Þ have different directions. However In general Q ~ jffij ~ the assumptions that j Q Qð0Þ j and that the transmissivity is a scalar quantity lead to T ðx; yÞ ¼ T ð0Þ ðx; yÞ
j rhð0Þ ðx; yÞ j j rhðrefÞ ðx; yÞ j
ðA:1Þ
For 1D flow, Eq. (A.1) reduces to an exact solution of Eq. ~ jffij Q ~ð0Þ j (1), whereas for 2D flow the approximation j Q prevents the T field estimated from Eq. (A.1) to be an exact solution of Eq. (1). Therefore, for 2D flow, the CMM is applied with an iterative scheme, updating the T ð0Þ and hð0Þ fields at each iteration. In order to prevent the identified transmissivity from assuming an unjustifiably high value where the reference hydraulic gradient is small the transmissivity correction at each iteration k is given by DT ðkþ1Þ ¼ T ðkÞ
j rhðkÞ j j rhðrefÞ b j rhðrefÞ j
ðA:2Þ
where b is a weight which depends on the reference gradients. A.2. Determination of the seed field This subsection is focused on the determination of the T ð0Þ field, which is fundamental for a suitable initialization of the CMM. For a homogeneous aquifer with constant transmissivity s, source term F ð0Þ , and boundary conditions hðbcÞ the solution H to Eq. (1) is given by H ðsÞ ¼ hc þ s1 hs
ðA:3Þ
where hc and hs are the solutions to the following boundary value problems: ( r2 hs ¼ F ð0Þ on X ðA:4Þ hs ¼ 0 on @X ( r2 hc ¼ 0 on X ðA:5Þ hc ¼ hðbcÞ on @X The decomposition Eq. (A.3) permits to identify in an easy way a value of s for each cell of the grid such that $H is as close as possible to rhðrefÞ , or in other words the value of s 2 for which the function f ðsÞ ¼j rH ðsÞ rhðrefÞ j attains the minimum. The function f ðsÞ admits a stationary point at 2
smin ¼ j rhs j a1
ðA:6Þ ðrefÞ
where a ¼ ðrhc rh Þ rhs . If a < 0, then smin > 0, whereas if a > 0, then smin < 0 which is not physically consistent. Moreover, if a ! 0, e.g. when the source terms are small or null, smin can take unnatural high values. A straightforward method to overcome this problem is the assignment of an upper limit to s, say sup .
The source term F ð0Þ is the same as Eq. (5) but F ðriverÞ , and F ðterracesÞ are computed from Eqs. (2)–(4) with F ðrefÞ h ¼ h . Where a < 0, T ð0Þ is given by minðsmin ; sup Þ. Finally, T ð0Þ is extended where a > 0 by an interpolation with inverse distance weighting. ðspringsÞ
References Associazione Irrigazione Est Sesia, 1979. Le acque sotterranee della pianura irrigua Novarese-Lomellina. Studi e ricerche per la realizzazione di un modello gestionale, Technical report. Beatrizotti, G., Hansen, J.W., Spocci, R., 1983. Optimierung der beno¨tigten daten fu¨r ein numerisches Modell der Grundwasserbewirtschaftung im Lockergestein. Gas Wasser Abwasser 63, 469–476. Benoit, N., Pasquier, P., Marcotte, D., Nastev, M., 2005. Conditional stochastic inverse modelling of the Chaˆteauguay river aquifers. ModelCARE 2005, Pre-published proceedings, pp. 515–521. Carrera, J., 1988. State of art of the inverse problem applied to the flow and solute transport equations. In: Custodio, E. et al. (Eds.), Groundwater Flow and Quality Modeling. Kluwer Academic, Dordrech, pp. 549–583. Carrera, J., Alcolea, A., Medina, A., Hidalgo, J., Slooten, L.J., 2005. Inverse problems in hydrogeology. Hydrogeology Journal 13, 206– 222. D’Auria, G., Zavagno, F., 2005. I fontanili della provincia di Cremona, Provincia di Cremona. Facchi, A., Ortuani, B., Maggi, D., Gandolfi, C., 2004. Coupled SVATgroundwater model for water resources simulation in irrigated alluvial plains. Environmental Modelling and Software 19, 1053–1063. Giudici, M., 2001. Development, calibration and validation of physical models. In: Clarke, K.C., Parks, B.O., Krane, e M.C. (Eds.), Geographic Information Systems and Environmental Modeling. Prentice-Hall, Upper Saddle River (NJ), pp. 100–121. Giudici, M., Vassena, C., 2007. Spectral analysis of the balance equation of ground water hydrology. Transport in Porous Media. doi:10.1007/ s11242-007-9142-3. Pasquier, P., Marcotte, D., 2005. Solving the groundwater inverse problem by successive flux estimation. Geostatistics for Environmental Applications: Proceedings of the Fifth European Conference on Geostatistics for Environmental Applications, 297–308. Pasquier, P., Marcotte, D., 2006. Steady- and transient-state inversion in hydrogeology by successive flux estimation. Advances in Water Research 29, 1934–1952. Ponzini, G., Crosta, G., 1988. The comparison model method: a new arithmetic approach to the discrete inverse problem of groundwater hydrology, 1, One-dimensional flow. Transport in Porous Media 3, 415–436. Ponzini, G., Lozej, A., 1982. Identification of aquifer transmissivities: the comparison model method. Water Resources Research 18, 597– 622. Ponzini, G., Crosta, G., Giudici, M., 1989. Identification of thermal conductivities by temperature gradient profiles: one-dimensional steady state flow. Geophysics 54, 643–653. Ponzini, G., Giudici, M., Vassena, C., 2007. Comments on ‘‘Steady- and transient-state inversion in hydrogeology by successive flux estimation” by P. Pasquier and D. Marcotte. Advances in Water Research 30, 2051–2053. Scarascia, S., Ponzini, G., 1972. An approximate solution for the inverse problem in hydraulics. L’Energia Elettrica 49, 518–531. Yeh, W.-G.W., 1986. Review of parameter identification procedures in groundwater hydrology: the inverse problem. Water Resources Research 22, 95–108.
Physics and Chemistry of the Earth 33 (2008) 1105–1110 www.elsevier.com/locate/pce
The importance of observations on fluxes to constrain ground water model calibration Chiara Vassena *, Cinzia Durante, Mauro Giudici, Giansilvio Ponzini Universita´ degli Studi di Milano, Dip. di Scienze della Terra, Milano, Italy Received 8 October 2007; received in revised form 18 January 2008; accepted 28 January 2008 Available online 5 February 2008
Abstract The aquifer system in the alluvial basin bordered by Adda, Po and Oglio rivers (Northern Italy) is characterised by a dual flow regime. In shallow sediments, which constitute a phreatic aquifer with high conductivity, great fluxes are driven by the interaction between ground water and the network of surface water, by the infiltration of rain and irrigation water, and by the fluxes drained from depression springs and river valley terraces. The underlying semiconfined aquifers are characterised by minor fluxes driven by water abstraction from wells of the public Water Works. Since most of the ground water flow occurs in the phreatic aquifer, an equivalent single layer 2D steady state flow model has been calibrated. The identification of the transmissivity field at the scale of the model has been obtained by solving an inverse problem with the comparison model method which requires an initial configuration, i.e., reference head, initial transmissivity field, source terms. Most of the head and source data are related to the phreatic aquifer, but most of the estimates of transmissivity are obtained with field tests conducted in deep wells pumping from the semiconfined aquifers, so that this kind of prior information cannot be used directly for model calibration. The inverse problem is underdetermined and a unique solution is not available. Furthermore information on surface hydrology is poor. Therefore many tests with different hypotheses about the initial configuration have been performed and some of them have been selected and used to initialise the automatic inversion procedure. Ó 2008 Elsevier Ltd. All rights reserved. Keywords: Mathematical modeling; Underdetermined inverse problem; Model calibration; Interaction between surface and ground water
1. Introduction The highly developed agricultural system of the Province of Cremona (Northern Italy) needs a great amount of water for irrigation which is presently derived from Adda, Po, Oglio and Serio rivers (Fig. 1). Limitations imposed by regional regulators and interannual climatic variations induce several farmers to ask for the permission to drill new water wells in the phreatic aquifer. Therefore the public administration has decided to develop a tool to support decisions and sustainable management of water resources. A research team of the University of Milan, merging the expertise of the Dipartimento di Scienze della Terra and of *
Corresponding author. E-mail address: [email protected] (C. Vassena).
1474-7065/$ - see front matter Ó 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.pce.2008.01.004
the Istituto di Idraulica Agraria (IIA), has developed a ground water flow model in the alluvial basin bordered by Adda, Po and Oglio rivers. The focus of this paper is on the model calibration. In fact, field tests provide values of the physical parameters that are relevant at scales different from those at which regional groundwater flow is modeled, so that the model parameters have to be computed as the solution of an inverse problem. Several review papers (e.g., Yeh, 1986; Carrera, 1988; Giudici, 2001; Carrera et al., 2005) provide comprehensive lists of the inverse methods proposed in groundwater hydrology and discuss the issues related to their application. For this case study the groundwater flow model has been calibrated with the comparison model method (CMM). This method was originally proposed by Scarascia and Ponzini (1972), successively developed by Ponzini and Lozej
1106
C. Vassena et al. / Physics and Chemistry of the Earth 33 (2008) 1105–1110
Fig. 1. Study area and contour plots of reference (solid line) and modelled (dashed line) piezometric head fields (years 2001–2003). Contour intervals: 5 m.
(1982), Ponzini and Crosta (1988), Ponzini et al. (1989), and applied to real aquifers (Associazione Irrigazione Est Sesia, 1979; Beatrizotti et al., 1983; Pasquier and Marcotte, 2005; Benoit et al., 2005). Recent advances about the application of the CMM to 3D transient flow have been proposed by Pasquier and Marcotte (2006) and Ponzini et al. (2007). The CMM requires the reference head field, the source terms and an initial tentative transmissivity field, which fix the water fluxes in the aquifer system. Data are of paramount importance for model development, especially for the model calibration and validation. Actually, data are often affected by errors, because it is difficult to perform experiments under controlled conditions in the field; furthermore, the modeler is forced to use already available data, often irregularly spaced and lacking in a check of the measurement procedures and accuracy. These remarks imply that the model predictions are not expected to be more accurate than the data used to develop and calibrate the model itself (Giudici, 2001). In this case study, data have been collected for purposes different from the development of a mathematical model so that some information is missing or uncertain. For instance, estimates of transmissivity are available from field tests performed in underlying semiconfined aquifers only, so that the initial transmissivity of the phreatic aquifer cannot be inferred from these data. The source terms, which include recharge from rain and irrigation, water extraction for drinking purposes and for agricultural and industrial use, drainage from depression springs and morphological terraces, recharge from artificial channels, and river-aquifer exchanges, are also estimated with uncertainty. In other words this case study is an example of an hydrological basin for which field data are not optimal, with reference
both to the space and time distribution and to the measurement quality. Therefore a sequence of tests coupling the CMM with different hypotheses about the source terms has been conducted. The choice of the best configuration (namely, parameters for the source terms) cannot be based on head data only, but other hydrological data are necessary, e.g. river discharge, spring fluxes, etc. In the case study this additional information is poor; as a consequence few configurations have been selected satisfying weak constraints and have been used to initialise the automatic inversion procedure.
2. The hydrogeological framework The study area (about 2500 km2) is located in the alluvial Po plain and is bordered by Adda, Po and Oglio rivers (Fig. 1); a three-year-long interval was considered for this study, from 2001 to 2003, because for this period the most complete data set was available. Current flow conditions of the aquifer system are still similar, even if the extraordinary repetition of drought periods could induce variation in the next years. The aquifer system is characterised by a dual flow regime. Shallow gravel and sandy-gravel sediments constitute a phreatic aquifer with high conductivity and an average thickness of about 50 m; the interaction between ground water and the network of surface water, the infiltration of rain and irrigation water, and the drainage from depression springs and river valley terraces induce large fluxes in this phreatic aquifer. The underlying semiconfined aquifers are discontinuous and characterised by minor
C. Vassena et al. / Physics and Chemistry of the Earth 33 (2008) 1105–1110
fluxes driven by leakage and by water extraction from wells of the public Water Works. The aquifer is fed by the inflow from north and by the recharge terms, subdivided into the recharge due to infiltration of the rain and of the water used for irrigation, F ðrechargeÞ , and the losses from irrigation channels, F ðchannelsÞ . Both F ðrechargeÞ and F ðchannelsÞ have been estimated by IIA with a model that computes a water mass balance at the ground surface and in the unsaturated zone, taking into account soil characteristics, the irrigation system, meteorological data, land use, crop production and growth (Facchi et al., 2004). The annual average estimates amount to 34 m3/s for F ðrechargeÞ and 29 m3/s for F ðchannelsÞ . Adda, Oglio and Po rivers are mainly acting as drains for the phreatic aquifer; Serio river is also draining the aquifer along most of its valley ðF ðriverÞ Þ. Drainage fluxes from depression springs, F ðspringsÞ , and from geomorphological terraces bounding the fluvial valleys, F ðterracesÞ , are sink terms for the mathematical model of the phreatic aquifer. They have to be added to the estimates of water abstraction fluxes, F ðabsÞ , for agricultural and industrial uses (about 11 m3/s from the phreatic aquifer) and for drinking water (about 3 m3/s from the semiconfined aquifer). We do not have reliable estimates of river-aquifer exchange fluxes. Ground water abstraction for irrigation purposes is estimated from IIA, whereas ground water uptake for industrial use is estimated from the database of the authorized water wells available from the local administrations. Abstraction for drinking water is known by measurements of the public Water Works. An interpolated head field hðintÞ has been obtained by kriging the piezometric head data (the time-averaged values measured in about 150 shallow wells from 2001 to 2003), the water levels of depression springs estimated at about 300 springs (values observed during a thorough survey in 2000), and the river levels estimated from about 300 transversal topographic sections measured along the main rivers (field observations support the assumption that these rivers are in direct hydraulic contact with the phreatic aquifer). A number of estimates of transmissivity are available from field tests in semiconfined aquifers, but we cannot profit from these data to infer the initial transmissivity of the phreatic aquifer, for which only four pumping tests are available. 3. The groundwater flow model The aquifer system of the Province of Cremona is modelled as an equivalent single-layer aquifer with 2D pseudostationary flow. The modelled piezometric head, h, is the solution of the forward problem:
r ðT rhÞ ¼ F h ¼ hðbcÞ
on X on @X
ð1Þ
where X is the aquifer domain and oX its border, T is the transmissivity, hðbcÞ is the prescribed head at the border, F is the source term.
1107
Dirichlet boundary conditions have been assigned, assuming hydraulic direct contact between Adda, Po, Oglio rivers and the aquifer. The source term, F, includes the following components. Interaction between the aquifer and the Serio river, F ðriverÞ , is modelled as ( h > bðriverÞ C ðriverÞ ðh hðriverÞ Þ ðriverÞ ¼ F C ðriverÞ ðhðriverÞ bðriverÞ þ lðriverÞ Þ h < bðriverÞ ð2Þ where h and hðriverÞ are the heads, respectively, of the aquifer and of the river; bðriverÞ and lðriverÞ are, respectively, the upper height (above m.s.l.) and the thickness of the river bed, C ðriverÞ is the conductance of the river bed. The drainage from depression springs, F ðspringsÞ , and the drainage from the geomorphological terraces, F ðterracesÞ , are modelled as F ðspringsÞ ¼ maxð0; C ðspringsÞ ðh bðspringsÞ ÞÞ F
ðterracesÞ
¼ maxð0; C
ðterracesÞ
ðh b
ðterracesÞ
ÞÞ
ð3Þ ð4Þ
where bðspringsÞ and bðterracesÞ are the drain base levels of, respectively, depression springs and geomorphological terraces and are estimated from field surveys (D’Auria and Zavagno, 2005) and from topographic maps; C ðspringsÞ and C ðterracesÞ are conductance coefficients of, respectively, springs at topographical depressions and terraces. The values of C ðriverÞ ; C ðspringsÞ ; C ðterracesÞ are assumed constant in space. Eq. (1) is solved with a conservative finite difference scheme applied to a square grid with spacing of 500 m, chosen on the basis of the data density and the goals of the model. The approximate solution is computed with an error on the global balance as small as 0.01% of the total inflow rate. 4. The model calibration The CMM, described in detail in Appendix A, is based on the forward solution of an auxiliary model, the comparison model, CM, which shares the same geometry, the same state equation, the same source terms, the same initial and boundary conditions, the same discretization as the model to be calibrated, but its transmissivity field is a tentative seed field. The transmissivity field is identified adapting iteratively the seed field with a procedure based on Darcy’s law and on the guess that the flow rate estimated from the solution to the CM is a good approximation of the real one. See Appendix A for a detailed explanation. The inverse problem is underdetermined and in principle admits infinite solutions. Furthermore in this case study no direct information on the values of C ðriverÞ ; C ðspringsÞ ; C ðterracesÞ is available and other source terms, for example recharge terms and losses from the irrigation channels, are estimated with some uncertainties. In order to explore a lot of possible solutions, several tests have been performed changing the values of C ðriverÞ ; C ðspringsÞ ; C ðterracesÞ and the source terms.
C. Vassena et al. / Physics and Chemistry of the Earth 33 (2008) 1105–1110
4.1. Search of the optimal source terms for the application of the CMM The source terms have been changed according to the following alternatives: values of C ðspringsÞ between 0.025 and 0.1 m2/s, values of C ðterracesÞ between 0.00625 and 0.05 m2/s and C ðriverÞ between 0.025 and 0.1 m2/s. These intervals were fixed after numerical tests and processing of field data. Values of the recharge terms and losses from channels correspondent to the 80% and to the 50% of the values estimated by IIA with a SVAT model have been also tested. For each test hðiniÞ is computed as the solution to (1) with a seed transmissivity field computed as described in section A.2, hðbcÞ ¼ hðintÞ and F ¼ F ðiniÞ , where F ðiniÞ ¼ F ðabsÞ F ðrechargeÞ þ F ðspringsÞ F ðchannelsÞ þ F ðterracesÞ þ F ðriverÞ
ð5Þ
Notice that F ðterracesÞ , F ðriverÞ and F ðspringsÞ depend on hðiniÞ . More than 2000 tests, fully mixing the alternatives described before, have been conducted and their performance has been evaluated with a simple statistics, MAE, the spatial average of the head residuals j hðiniÞ hðintÞ j. For more than 200 tests MAE < 1 m, which is a value consistent with the error on the head measurements and on the subsequent interpolation and regularisation; hence further criteria are needed to choose the best distribution of source terms. In order to reduce the uncertainty of the inverse solution due to the fact that this inverse problem is underdetermined, the following constraints, related to water fluxes, have been fixed: (1) total drainage from Serio river <7 m3/s, based on the observation that the average of Serio river discharge is about 20 m3/s and water uptake for irrigation in the area is about 5 m3/s; (2) total drainage from springs >10 m3/s, based on the observation that the average annual discharge measured in channels fed by springs is about 20 m3/s; however these channels are partly fed also by the irrigation network. Recall that all these values are time averages from 2001 and 2003, since the model is calibrated for the average annual flow. After the application of these constraints, 42 tests have been selected, and in order to extract the most promising ones, some further weak constraints have been added to the previous criteria:
According to these criteria five tests have been selected: for these tests MAE varies between 0.83 and 0.91 m. All tests share the same values of C ðspringsÞ and C ðriverÞ , respectively, 0.05 and 0.025 m2/s. Instead C ðterracesÞ ¼ 0:025 m2 =s for Tests 1 and 2 and C ðterracesÞ ¼ 0:05m2 =s for the remaining tests. For all tests the other terms of the balance, i.e., the exchanges through boundaries, the drainage terms and the recharge terms, differ from the values averaged among the five tests by less than 20%. 4.2. Results As an example of calibration results, the transmissivity field obtained using Test 1 as the initial configuration for the CMM is shown. After four iterations of the CMM, MAE < 0.45 m. The identified transmissivity field permits a good fit between observed and simulated heads (Fig. 2), shows the large scale heterogeneity of the aquifer system (Fig. 3) and permits a preliminary evaluation of the groundwater balance. A consistency check has been performed along Serio river. In fact a reference head field, independent from ðrefÞ river-aquifer interaction, hnoSerio , has been obtained by interpolating and filtering (see Section A.1) the data described in Section 2 excluding the hydrodynamic levels along the Serio river. Also F ðnoSerioÞ has been obtained filtering the difference F ðiniÞ F ðriverÞ , and a transmissivity field, T ðnoSerioÞ , has been identified. The balance error computed as BE ¼ r ðT ðnoSerioÞ rhðrefÞ Þ F ðnoSerioÞ is closed to F ðriverÞ and the similarity between the T field plotted in Fig. 3 and T ðnoSerioÞ can be considered as a validation of the calibration results and of the value of C ðriverÞ .
115 105
Modelled piezometric head (m)
1108
95 85 75 65 55 45 35
(1) more than 50% of the total number of springs should be active; this prevents from modeling few active springs with unrealistically high discharge rates; R (2) X ðF ðspringsÞ þ F ðterracesÞ þ F ðriverÞ Þ dX > 20 m3 =s, according to the fact that the average annual discharge measured in channels fed by springs and morphological terraces is about 20 m3/s.
25 15 15
25
35
45
55
65
75
85
95
105
115
Observed piezometric head (m) Fig. 2. Comparison between observed and modelled piezometric heads. For most of the points the residual error is less than ±2 m (black lines).
C. Vassena et al. / Physics and Chemistry of the Earth 33 (2008) 1105–1110
1109
Fig. 3. Identified LogT field (m2/s) for Test 1 (C ðriverÞ ¼ 0:025 m2 =s, C ðspringsÞ ¼ 0:05 m2 =s, C ðterracesÞ ¼ 0:025 m2 =s).
5. Conclusions Calibration of ground water flow models often does not provide a unique solution; this is the case for this aquifer, where a great number of configurations of source terms, and therefore of T fields obtained with the inverse procedure, are found. Moreover the results of the CMM depend also on the choice of the initial seed T field, T ð0Þ , and of the reference head field; the procedure to obtain T ð0Þ proposed in Appendix A proved to be the most robust among different alternatives and therefore can be applied with confidence to study poorly gauged hydrological basins. The lack of accurate data on surface hydrology prevented us from assigning strict constraints on the inversion results; nevertheless, constraints on exchange fluxes between surface and ground water were useful to limit the uncertainty. This shows, on one hand, the importance of the efforts aimed to monitor hydrological basins with the acquisition of field or remote, hard or soft data, and, on the other hand, the utility of inverse methods to check the consistency and possible deficiency of a data set. Moreover the application of inverse methods like the CMM to poorly gauged basins claims developments to handle the lack of information about the source terms in the inversion. As a general final remark, this application shows that merging surface and groundwater study is important not only for an integrated water balance at the catchment scale, but also for the characterization of the aquifer system. Acknowledgements This work was financially supported by the Province of Cremona. Dr. M. Pesaro and Dr. M.M. Cremonini are
kindly acknowledged for the continuous support and for releasing the permission to publish this paper. Prof. C. Gandolfi and his coworkers at the Istituto di Idraulica Agraria of the University of Milan are kindly acknowledged for the cooperation during the work, in particular for the modeling of the irrigation system and the estimate of the aquifer recharge. Appendix A. The comparison model method A.1. Basic properties of the method The CMM is explicitely based on Darcy’s law and on the forward solution of an auxiliary model, the comparison model, CM, which gives the solution hð0Þ , to the forward problem (1) for T ¼ T ð0Þ , F ¼ F ðrefÞ and hðbcÞ ¼ hðrefÞ . The reference head field hðrefÞ is obtained from the application of an averaging procedure (high-wavenumber-cut filter) to the piezometric head interpolated from the observations hðintÞ ; notice that hðrefÞ ¼ hðintÞ on oX. The reference source term F ðrefÞ is obtained from the application of a similar filter to F ðiniÞ . This filtering procedure has been introduced to limit the instability effects of high-wavenumber components of h on inversion (Giudici and Vassena, 2007) and is of paramount importance to obtain reliable results. Then rhð0Þ is compared with the observed hydraulic gradient rhðrefÞ , the basic assumption is that the tentative seed transmissivity field generates an hð0Þ field such that the auxiliary flow rate per unit length ~ Qð0Þ ¼ T ð0Þ rhð0Þ ~ ¼ T rhðrefÞ is a good approximation of the real one, Q where T is the unknown transmissivity field to be estimated.
1110
C. Vassena et al. / Physics and Chemistry of the Earth 33 (2008) 1105–1110
~ and Q ~ð0Þ have different directions. However In general Q ~ jffij ~ the assumptions that j Q Qð0Þ j and that the transmissivity is a scalar quantity lead to T ðx; yÞ ¼ T ð0Þ ðx; yÞ
j rhð0Þ ðx; yÞ j j rhðrefÞ ðx; yÞ j
ðA:1Þ
For 1D flow, Eq. (A.1) reduces to an exact solution of Eq. ~ jffij Q ~ð0Þ j (1), whereas for 2D flow the approximation j Q prevents the T field estimated from Eq. (A.1) to be an exact solution of Eq. (1). Therefore, for 2D flow, the CMM is applied with an iterative scheme, updating the T ð0Þ and hð0Þ fields at each iteration. In order to prevent the identified transmissivity from assuming an unjustifiably high value where the reference hydraulic gradient is small the transmissivity correction at each iteration k is given by DT ðkþ1Þ ¼ T ðkÞ
j rhðkÞ j j rhðrefÞ b j rhðrefÞ j
ðA:2Þ
where b is a weight which depends on the reference gradients. A.2. Determination of the seed field This subsection is focused on the determination of the T ð0Þ field, which is fundamental for a suitable initialization of the CMM. For a homogeneous aquifer with constant transmissivity s, source term F ð0Þ , and boundary conditions hðbcÞ the solution H to Eq. (1) is given by H ðsÞ ¼ hc þ s1 hs
ðA:3Þ
where hc and hs are the solutions to the following boundary value problems: ( r2 hs ¼ F ð0Þ on X ðA:4Þ hs ¼ 0 on @X ( r2 hc ¼ 0 on X ðA:5Þ hc ¼ hðbcÞ on @X The decomposition Eq. (A.3) permits to identify in an easy way a value of s for each cell of the grid such that $H is as close as possible to rhðrefÞ , or in other words the value of s 2 for which the function f ðsÞ ¼j rH ðsÞ rhðrefÞ j attains the minimum. The function f ðsÞ admits a stationary point at 2
smin ¼ j rhs j a1
ðA:6Þ ðrefÞ
where a ¼ ðrhc rh Þ rhs . If a < 0, then smin > 0, whereas if a > 0, then smin < 0 which is not physically consistent. Moreover, if a ! 0, e.g. when the source terms are small or null, smin can take unnatural high values. A straightforward method to overcome this problem is the assignment of an upper limit to s, say sup .
The source term F ð0Þ is the same as Eq. (5) but F ðriverÞ , and F ðterracesÞ are computed from Eqs. (2)–(4) with F ðrefÞ h ¼ h . Where a < 0, T ð0Þ is given by minðsmin ; sup Þ. Finally, T ð0Þ is extended where a > 0 by an interpolation with inverse distance weighting. ðspringsÞ
References Associazione Irrigazione Est Sesia, 1979. Le acque sotterranee della pianura irrigua Novarese-Lomellina. Studi e ricerche per la realizzazione di un modello gestionale, Technical report. Beatrizotti, G., Hansen, J.W., Spocci, R., 1983. Optimierung der beno¨tigten daten fu¨r ein numerisches Modell der Grundwasserbewirtschaftung im Lockergestein. Gas Wasser Abwasser 63, 469–476. Benoit, N., Pasquier, P., Marcotte, D., Nastev, M., 2005. Conditional stochastic inverse modelling of the Chaˆteauguay river aquifers. ModelCARE 2005, Pre-published proceedings, pp. 515–521. Carrera, J., 1988. State of art of the inverse problem applied to the flow and solute transport equations. In: Custodio, E. et al. (Eds.), Groundwater Flow and Quality Modeling. Kluwer Academic, Dordrech, pp. 549–583. Carrera, J., Alcolea, A., Medina, A., Hidalgo, J., Slooten, L.J., 2005. Inverse problems in hydrogeology. Hydrogeology Journal 13, 206– 222. D’Auria, G., Zavagno, F., 2005. I fontanili della provincia di Cremona, Provincia di Cremona. Facchi, A., Ortuani, B., Maggi, D., Gandolfi, C., 2004. Coupled SVATgroundwater model for water resources simulation in irrigated alluvial plains. Environmental Modelling and Software 19, 1053–1063. Giudici, M., 2001. Development, calibration and validation of physical models. In: Clarke, K.C., Parks, B.O., Krane, e M.C. (Eds.), Geographic Information Systems and Environmental Modeling. Prentice-Hall, Upper Saddle River (NJ), pp. 100–121. Giudici, M., Vassena, C., 2007. Spectral analysis of the balance equation of ground water hydrology. Transport in Porous Media. doi:10.1007/ s11242-007-9142-3. Pasquier, P., Marcotte, D., 2005. Solving the groundwater inverse problem by successive flux estimation. Geostatistics for Environmental Applications: Proceedings of the Fifth European Conference on Geostatistics for Environmental Applications, 297–308. Pasquier, P., Marcotte, D., 2006. Steady- and transient-state inversion in hydrogeology by successive flux estimation. Advances in Water Research 29, 1934–1952. Ponzini, G., Crosta, G., 1988. The comparison model method: a new arithmetic approach to the discrete inverse problem of groundwater hydrology, 1, One-dimensional flow. Transport in Porous Media 3, 415–436. Ponzini, G., Lozej, A., 1982. Identification of aquifer transmissivities: the comparison model method. Water Resources Research 18, 597– 622. Ponzini, G., Crosta, G., Giudici, M., 1989. Identification of thermal conductivities by temperature gradient profiles: one-dimensional steady state flow. Geophysics 54, 643–653. Ponzini, G., Giudici, M., Vassena, C., 2007. Comments on ‘‘Steady- and transient-state inversion in hydrogeology by successive flux estimation” by P. Pasquier and D. Marcotte. Advances in Water Research 30, 2051–2053. Scarascia, S., Ponzini, G., 1972. An approximate solution for the inverse problem in hydraulics. L’Energia Elettrica 49, 518–531. Yeh, W.-G.W., 1986. Review of parameter identification procedures in groundwater hydrology: the inverse problem. Water Resources Research 22, 95–108.