Numerical approximation of head and ¯ux covariances in three dimensions using mixed ®nite elements Andrew I. James * & Wendy D. Graham Interdisciplinary Program in Hydrologic Sciences, University of Florida, Gainesville, FL 32611, USA (Received 14 April 1998; revised 27 September 1998; accepted 8 October 1998) A numerical method is developed for accurately approximating head and ¯ux covariances and cross-covariances in ®nite two- and three-dimensional domains using the mixed ®nite element method. The method is useful for determining head and ¯ux covariances for non-stationary ¯ow ®elds, for example those induced by injection or extraction wells, impermeable subsurface barriers, or non-stationary hydraulic conductivity ®elds. Because the numerical approximations to the ¯ux covariances are obtained directly from the solution to the coupled problem rather than having to dierentiate head covariances, the approximations are in general more accurate than those obtained from conventional ®nite dierence or ®nite element methods. Results for uniform ¯ow example problems are consistent with results from previously published ®nite domain analyses and demonstrate that head variances and covariances are quite sensitive to boundary conditions and the size of the bounded domain. Flux variances and covariances are less sensitive to boundary conditions and domain size. Results comparing approximations from lower-order Raviart±Thomas±Nedelec and higher order Brezzi±Douglas±Marini 9 ®nite element spaces indicate that higher order element space improve the esti- mate of the ¯ux covariances, but do not signi®cantly aect the estimate of the head covariances. Ó 1999 Elsevier Science Ltd. All rights reserved 1 INTRODUCTION Over the past two decades, stochastic methods have been widely used to describe the variability of ground- water systems. Of long standing interest is the determi- nation of head and velocity covariances, and head/ conductivity, head/velocity, and velocity/conductivity cross-covariances. These covariance functions have been determined analytically for two- and three-dimensional in®nite domains by Bakr et al., 3 Mizell et al., 20 Graham and Mclaughlin, 13 Rubin, 27 Rubin and Dagan, 30 and Zhang and Neumann. 37,38 Osnes, 24 Rubin and Da- gan 28,29 and Na and Vecchia 22 have used semi-analyt- ical methods to solve for head covariances in bounded domains. Osnes 25 used semi-analytical methods to solve for velocity covariances in two dimensional bounded domains. Numerical solutions for head covariances and head-conductivity cross-covariances in bounded do- mains have also been reported by McLaughlin and Wood, 18,19 Sun and Yeh, 35 Li and McLaughlin, 17 and Van Lent and Kitanidis. 36 Typically, when a numerical approach is taken, the covariances and cross-covari- ances involving velocity are obtained by numerically dierentiating the covariance functions of head. This can lead to loss of accuracy in the resulting functions if head covariances change rapidly over small distances within the domain, as may occur if conductivity is highly variable, if complex boundary conditions are present, and/or source and sink terms are present. Sun and Yeh 33 and Sun 32 point out that, in general, accurate numerical calculations of gradients of unknown functions in cou- pled inverse problems can be dicult. In this paper we examine a method for numerically evaluating covariances and cross-covariances between log hydraulic conductivity, head, and Darcy ¯ux using a Advances in Water Resources Vol. 22, No. 7, pp 729±740, 1999 Ó 1999 Elsevier Science Ltd Printed in Great Britain. All rights reserved 0309-1708/99/$ ± see front matter PII:S0309-1708(98)00044-X * Corresponding author. 729