A New Approach to Generate Source Terms, Distributed in Space and Time for Transient Groundwater Flow Models

As presented in earlier papers (Fuchs and Fank, 1998; Fank, 1999) it is possible to generate initial and boundary conditions for transient ground water flow models using a set of measured groundwater hydrographs. In this paper we present • geostatistical models for the evaluation of spatial and time distributed recharge (RC) using measured ground water hydrographs, • limitations concerning the type of aquifer and the hydrological environment and • application to the western part of the “Leibnitzer Feld”, a shallow quaternary aquifer south of Austria. Zusammenfassung: Wie bereits in fr̈ uheren Arbeiten gezeigt (Fuchs and Fank, 1998; Fank, 1999) ist es mit Hilfe von Grundwasserstandsganglinien möglich, Anfangsund Randbedingungen f ür ein Grundwasserstr ömungsmodell zu generieren. In dieser Arbeit werden • auf Basis von gemessenen Grundwasserstandsganglinien geostatistische Modelle für die Berechnung r äumlich und zeitlich verteilter Grundwasserneubildung, • Grenzen bez̈ uglich des Aquifers und der hydrogeologischen Umgebung und • eine Anwendung f̈ ur den westlichen Teil des Leibnitzer Feldes, einen seichtliegenden Aquifer im S̈ udenÖsterreichs


Introduction
Quaternary sediments deposited from rivers in former glacial and peri-glacial regions play an important role for drinking water supply and distribution.The groundwater recharge from infiltrating precipitation leads to a high vulnerability to harmful chemicals.The increasing intensity of cultivation causes an increasing impact of nitrogen to the aquifer.This is especially important in the valleys of the Alps or their surrounding, which are used as settlement areas, areas for agricultural use, traffic areas and other economic purposes (gravel and sand are being used as building material).Numerical modelling of the movement of water and solids in the unsaturated and saturated zones of an aquifer enables to predict the quantity of impact of different measures on to the whole system.However, setting up a model and determining model parameters still calls for substantial interdisciplinary research efforts.Close ties between natural scientific and socio-economic components in combination with mathematical-technical approaches will lead to results which provide the basis for lawmakers for a better definition of the legal constraints of the various land use practices.Statistical and geostatistical evaluation of hydraulic head data may help to define initial and boundary conditions for transient groundwater models.The investigation area, the so called 'Leibnitzer Feld' is located in the lower Mur valley (Styria, Austria), which is filled up with quaternary gravels and sands with relatively high permeability and a thickness of 10 to 15 m.They represent an important aquifer which is intensively used for common water supply.The groundwater table monitoring network is described in Fuchs et al. (1995).For a part of the 'Leibnitzer Feld' with an extent of up to 15 km 2 (see Figure 1) a transient finite element groundwater flow model shall be developed.For the application of this groundwater flow model the definition of boundary values (hydraulic head) and initial conditions as well as space distributed system parameters (storage coefficient and permeability) for the solution of the partial differential flow equations are needed.The data set for the calibration of the model and the estimation of transient boundary conditions consist of 26 wells (see Figure 1), where the groundwater table was measured from 1992 to 1996 in weekly intervals.Using principal component analyses (PCA) former investigations showed, that-in homogeneous shallow phreatic aquifers-transient hydraulic head boundary conditions may be estimated using spatio-temporal models (Fuchs and Fank, 1998).The PCA (Flury and Riedwyl, 1983) is a mathematical transformation of the variables x ij and can be used to reduce highly correlated data sets to an uncorrelated form for the characterization and visualization of their real structure.The measured p hydrographs x i = (x 1j , . . ., x pj ) T , j = 1, . . ., n are calculated as linear functions of new vectors y i = (y 1j , . . ., y pj ) T : with u l = (u 1j , . . ., u pj ).The vectors y i , the PC (principal components), are calculated after the transformation solving the special Eigenvalue problem (S − λE)u = 0, having the real roots λ 1 ≥ . . .≥ λ p (S being the estimated covariance matrix from the n observations x ij and E identity matrix with rank p).The sum p of the eigenvalues λ describes the whole variability.The part of the variance explained by the first PC can be calculated by λ 1 /p.If the first principal component-as the result of a principal component transformation performed on the data set-explains more than 95 % of the whole variance, this can be used to analyses the spatial structure with respect to the temporal behavior.Kriging technique can be used to estimate the spatial distribution of the factor scores of the first principal component.Using the calculated kriging weights of unobserved grid points and the time series information of its explanatory observation wells, water table hydrographs (GW-HD) can be estimated.Using the provoked spatial-temporal models the definition of first kind boundary values for transient finite element flow models is possible (see Figure 2).The estimation of the storage coefficient (SC) may be done in combination of hydraulic head data and information about the unsaturated zone Fank (1999).In areas and time periods, where the rise of the groundwater table at shallow phreatic aquifers is mainly a reaction on the amount of percolating water through the unsaturated zone to groundwater it is possible to estimate the SC of the aquifer on infiltration events.In time periods where field capacity is reached in the unsaturated zone groundwater recharge (RC) is equal to the amount of precipitation minus evapotranspiration (calculated from meteorological parameters).SC in this case equals the amount of RC [mm] divided by the rise of the groundwater table [mm], taking into account the hydraulic drawdown of groundwater without recharge during the same period.The spatial distribution of SC is estimated using the kriging technique.The starting values of the hydraulic conductivity are taken from a steady state calibrated groundwater flow model at low water conditions, neglecting groundwater recharge from infiltrating precipitation.Assuming that the transient 1st order boundary conditions at the border of the model area and the estimation of time and space distributed RC are well estimated and therefor fixed values, hydraulic conductivity and SC are the calibration parameters for the transient groundwater flow model.
In this paper we want to present two solutions based on geostatistical models for the evaluation of spatial and time distributed RC using measured ground water hydrographs (see Figure 2) and their application to a shallow quaternary aquifer south of Austria (see Figure 1).The initial data set contains p=26 GW-HD measured weekly from 920701 to 951226 (i.e.u=208 different time steps).Figure 1 shows the geographical distribution of the GW-HD inside the test area.Both solutions depend on a method to calculate groundwater recharge from the fluctuation of the groundwater table in time presented in Fank (1999).Assuming that RC mainly depends on the infiltration of precipitation (other forms of RC must be negligible or constant in time), knowing the hydraulic drawdown of the groundwater table without recharge (estimated using the falling parts of groundwater hydrographs in dry periods) and the SC of the aquifer (estimation see before) RC for one time step (time period between two measurement dates of the groundwater table) equals the difference between the groundwater table at the end of the time step minus the groundwater table at the beginning of the time step minus the hydraulic drawdown during the time step multiplied with SC.RC for longer time periods is calculated by the summing up of RC for time steps.Time steps have to be set that the measured GW-HD explains the real fluctuation of groundwater exactly.For the investigation area the method is proved comparing the results to lysimeter measurements and soil water model calculations (Fank, 1999).

Solution 1
Based on the initial spatio-temporal data set GW − HD i , i = 1, . . ., p for all observation points time series RC − HD i ,i = 1, . . ., p with u = 208 different values (i.e.weekly sums) are estimated using the method presented in Fank (1999).Then structural analysis (i.e.variography) on these RC-HD is performed.As RC exhibits nearly the same behavior in every direction at any time the experimental semivariograms γ t (h), t = 1, . . ., u for each time step t are calculated only for the isotropic overall case (see Figure 3) according to where • n h number of pairs separated by a distance h To ensure that kriging estimation works proper the spatial correlation 2γ(h) must be based on some "positive definite" theoretical models (Deutsch and Journel, 1997).Such models are fitted to the experimental variograms, and include spherical, exponential, Gaussian, etc. models.Combinations of such models are called nested variograms and are also valid.The experimental semivariograms γ i (h) show different shapes so nearly the whole range of theoretical models are used for variogram modelling.Figure 3 demonstrates for example empirical variograms (dotted lines) and their variogram models (bold solid lines) for 4 different time steps (930420,931012,940705,940830) of weekly sum of groundwater recharge.For 940830 a spherical model is used, 931012 represents a nested model, whereas 930420 and 940705 are linear models.Kriging is expressed in a socalled "kriging system", which relates the variogram between the samples, the variogram between each sample to the location to be estimated, and the unknown kriging weights.The ordinary kriging system can be expressed in terms of the semivariogram function γ(h) (see e.g.Journel and Huijbregts, 1978): Repeating this procedure 208 times we get the geostatistical recharge grids RC − M OD t , t = 1, . . ., u.At least for each grid element m, m = 1, . . ., n the appropriate 208 RC-values are assembled and so RC − HD m , m = 1, . . ., n as source terms for the transient groundwater flow model are obtained.

Solution 2
Like solution 1 solution 2 starts with the estimation of RC − HD i , i = 1, . . ., p from groundwater hydrographs using the method presented in Fank (1999).Parallel a principal component analyses (PCA) with GW − HD i , i = 1, . . ., p is performed (Fuchs and Fank, 1998).The PCA of the 26 groundwater hydrographs in the investigation area (see Figure 1) for the period 1992 to 1995 shows that the first PC explains more than 95 % of the whole variance.Therefore the factor scores of the first PC can be used for the analyses of the spatial structure with respect to the temporal behavior of the groundwater table.An experimental variogram of the factor scores is calculated and a linear variogram model is fitted.Based on this variogram model and on the factor scores of the first PC the ordinary kriging algorithm is performed to get the kriging weights KW m (F S i ), [m = 1, . . ., n; i = 1, . . ., p] of influencing observation points for each grid element at n (n p) predefined co-ordinates.The first PC is an overall variogram over observation time, therefore this model can be used for all observed time points with an error resulting from

Results
Table 1 and Figure 4 show a comparison of the results of estimation using solution 1 and solution 2 for the selected points EP1 ... EP6 (location 1).At the research station in Wagna the method to estimate RC from the fluctuation of the groundwater table in time was proofed in comparing the results to lysimeter measurements and soil water modelling (Fank, 1999).The residuals to the mean value of different methods for the period 1992 to 1995 are nearby 1 %, for the single years the differences are between -10 % and +8 % of RC.Taking into account that the input data to estimate RC are of very high quality at research station Wagna the significance of estimating yearly sum RC from the fluctuation of the groundwater in the investigation area may be of 85 %.The spatial distribution of RC at EP1 ... EP6 (see Table 1) represent the influence of different soil conditions and land use (agriculture, forestry, settlement) on RC as described in Fank (1999).Due to the hydro-meteorological conditions and the vegetation cycle of the culture crops in the investigation area the main time period for RC was in November and December with more than 30 % of yearly mean RC during the period 1992 to 1995.During the summer months (July and August) RC is less than 5 % of yearly mean RC.The method to estimate RC from the fluctuation of the groundwater table leads to nearly the same results (33 % and 6 %) for groundwater measurement point "Wagna".Mean monthly RC in Figure 4 for estimation points EP1 ... EP6 show the same fluctuation in time, the differences between them are the consequence of different soil conditions and land use practices.

Conclusions
The presented methods to estimate the distribution of RC in space and time using statistical and geostatistical algorithm may help to generate source terms for every node of groundwater flow models in areas where the application of soil water models is not possible due to the lack of input data.Both methods show very well comparable results in the yearly sum of RC as well as in the mean monthly sum of RC for the period 1992 to 1995.The overhead to the generation of the experimental variograms as well as the calibration of variogram models for longer time periods is very time expensive in comparison with the second solution.This process can be supported by automatic software procedures only in a low extent.Objective test procedures like cross validation are seldom used because for variogram calibration other relevant information (without a clear mathematical context) are needed.

Figure 1 :
Figure 1: Topographic situation with groundwater observation points and elected points for recharge estimation

Figure 2 :
Figure 2: Flow chart to visualize the generation of initial and boundary conditions for transient groundwater flow models using measured groundwater hydrographs.GW-HD = groundwater hydrograph, SC = storage coefficient, RC-HD = recharge hydrograph, RC-MOD = geostatistical recharge model, KW = kriging weights, FS = factor scores of the 1st principal component, i = number of groundwater hydrographs Lagrange parameter• λ j kriging weights • γ(s i , s j ) mean variogram value between samples • γ(s i , V ) mean variogram value between samples and location to be estimated Based on the theoretical variogram models the ordinary kriging algorithm of the RCvalues for time step t for every grid element m m = 1, . . ., n at n (n p) predefined co-ordinates gives the kriging weights having the contribution of every observation point to the resulting estimation value.

Figure 3 :
Figure 3: Empirical variogram structures, number of pairs at lag distances, variogram models for 4 different data sets of weekly sum of groundwater recharge and elected search radius for kriging interpolation

Figure 4 :
Figure 4: Comparison of mean monthly groundwater recharge (mm) from the period 1992 to 1995 at different locations (EP 1 to EP 6) evaluated using different solution methods

Table 1 :
Yearly sum of groundwater recharge (mm) at estimation points EP 1 to EP 6 comparing different methods from 1992 to 1995