© 2009 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim

In an industrial world there are unfortunately more and more polluted sites. These sites represent evidently a long-term danger for the groundwater, an consequently to our health. An active remediation of all these sites is practically impossible. In quantifying how dangerous a contaminated site is and if it shows a high natural attenuation potential, numerical simulations are playing an important role. The reliability of such simulations depends on many factors: the right mathematical model, suitable parameters, accurate numerical method, and so on. We present in this paper a mixed finite element method (MFEM) which uses stochastically generated parameters to simulate contaminant transport in heterogeneous aquifers.

The water flow through saturated aquifers (below the groundwater table) is mathematically modeled by

with $\psi (\mathbf{x})$$\psi (\mathbf{x})$psi(x)\psi(\mathbf{x})$\psi (\mathbf{x})$ denoting the pressure head, $\mathbf{Q}(\mathbf{x})$$\mathbf{Q}(\mathbf{x})$Q(x)\mathbf{Q}(\mathbf{x})$\mathbf{Q}(\mathbf{x})$ the water flux $K(\mathbf{x})$$K(\mathbf{x})$K(x)K(\mathbf{x})$K(\mathbf{x})$ the hydraulic conductivity, ${\mathbf{e}}_{\mathbf{z}}$${\mathbf{e}}_{\mathbf{z}}$e_(z)\mathbf{e}_{\mathbf{z}}${\mathbf{e}}_{\mathbf{z}}$ the vertical unit vector, $T$$T$TT$T$ a finite end time and $\mathrm{\Omega }\subset {\mathbb{R}}^d(d\ge 1)$$\mathrm{\Omega }\subset {\mathbb{R}}^d(d\ge 1)$Omega subR^(d)(d >= 1)\Omega \subset \mathbb{R}^{d}(d \geq 1)$\mathrm{\Omega }\subset {\mathbb{R}}^d(d\ge 1)$ the computational domain having a Lipschitz continuous boundary $\mathrm{\Gamma }$$\mathrm{\Gamma }$Gamma\Gamma$\mathrm{\Gamma }$. The equation (1) represents the mass conservation of water and (2) is the Darcy's law. Boundary conditions complete the model.

To describe the diffusive-advective-reactive transport of a one-component solute in saturated soil we consider the mathematical model:

with $c(t,\mathbf{x})$$c(t,\mathbf{x})$c(t,x)c(t, \mathbf{x})$c(t,\mathbf{x})$ denoting the concentration of the solute, $D$$D$DD$D$ the diffusion-dispersion coefficient and $r(\cdot )$$r(\cdot )$r(*)r(\cdot)$r(\cdot )$ a reaction term. Initial $c(t=0)=c_I$$c(t=0)=c_I$c(t=0)=c_(I)c(t=0)=c_{I}$c(t=0)=c_I$ and boundary conditions complete the model. We refer to [5,6] for a more extensive mathematical model including equilibrium/non-equilibrium sorption and saturated/unsaturated flow.

For the spatial discretization of the flow problem as well as of the transport we use the MFEM. The main advantages of the MFEM are the local mass conservation and, which is especially important for the flow problems, the water flux is an intrinsic part of the solution and it is continuous over the edges (sides) of the elements. The transport problem is discretized in time by the implicit Euler method. Unlike in stochastic finite element methods of estimating averages, this approach yields concentration fields and spatial concentration moments for samples of the random field. Ensembles of simulations can be used to estimate both mean values and fluctuations of observables of the transport process. Results of the stochastic simulations described here can be further used to assess the reliability for real cases of the ensemble average quantities provided by stochastic modeling of transport in groundwater.

The paper is structured as follows. In the next section the discretization method is presented. Section 3 introduces the stochastic approach. The fourth section shows first numerical results and we end with some conclusions.

In what follows we make use of common notations in the functional analysis. By $⟨\cdot ,\cdot ⟩$$⟨\cdot ,\cdot ⟩$(:*,*:)\langle\cdot, \cdot\rangle$⟨\cdot ,\cdot ⟩$ we mean the inner product on $L^2(\mathrm{\Omega })$$L^2(\mathrm{\Omega })$L^(2)(Omega)L^{2}(\Omega)$L^2(\mathrm{\Omega })$ and $\Vert \cdot \Vert$$\Vert \cdot \Vert$||*||\|\cdot\|$\Vert \cdot \Vert$ the norm on it. The functions in $H(\mathrm{div};\mathrm{\Omega })$$H(\mathrm{div};\mathrm{\Omega })$H(div;Omega)H(\operatorname{div} ; \Omega)$H(\mathrm{div};\mathrm{\Omega })$ are vector valued, having a $L^2$$L^2$L^(2)L^{2}$L^2$ divergence. Furthermore, we let ${\mathcal{T}}_h$${\mathcal{T}}_h$T_(h)\mathcal{T}_{h}${\mathcal{T}}_h$ be a regular decomposition of $\mathrm{\Omega }\subset {\mathbb{R}}^d$$\mathrm{\Omega }\subset {\mathbb{R}}^d$Omega subR^(d)\Omega \subset \mathbb{R}^{d}$\mathrm{\Omega }\subset {\mathbb{R}}^d$ into closed $d$$d$dd$d$-simplices; $h$$h$hh$h$ stands for the mesh-size. For the discretization in time we let $N\in \mathbb{N}$$N\in \mathbb{N}$N inNN \in \mathbb{N}$N\in \mathbb{N}$ be strictly positive, and define the time step $\tau =T/N$$\tau =T/N$tau=T//N\tau=T / N$\tau =T/N$. We denote by $W_h$$W_h$W_(h)W_{h}$W_h$ the space of piecewise constant functions and by $V_h$$V_h$V_(h)V_{h}$V_h$ the lowest order Raviart - Thomas ( $R{T}_0$$R{T}_0$RT_(0)R T_{0}$R{T}_0$ ) space:

We can now formulate the fully discrete schemes for solving the flow and transport equations. As already mentioned before, the schemes are based on Euler implicit method and MFEM. The fully discrete scheme for the flow problem (1)- (2 reads as

Problem 2.1 Find $({\psi }_h,{\mathbf{Q}}_{\mathbf{h}})\in W_h\times V_h$$\left({\psi }_h,{\mathbf{Q}}_{\mathbf{h}}\right)\in W_h\times V_h$(psi_(h),Q_(h))inW_(h)xxV_(h)\left(\psi_{h}, \mathbf{Q}_{\mathbf{h}}\right) \in W_{h} \times V_{h}$({\psi }_h,{\mathbf{Q}}_{\mathbf{h}})\in W_h\times V_h$ such that there holds

Remark 2.2 The analysis of the scheme above is included as a particular case in [3,4], where saturated/unsaturated flow is treated.

We define now the fully discrete scheme for the equation (3) at level $n\in \{1,\dots ,N\}$$n\in \{1,\dots ,N\}$n in{1,dots,N}n \in\{1, \ldots, N\}$n\in \{1,\dots ,N\}$ :
Problem 2.3 Let $c_h^{n-1}\in W_h$$c_h^{n-1}\in W_h$c_(h)^(n-1)inW_(h)c_{h}^{n-1} \in W_{h}$c_h^{n-1}\in W_h$ and ${\mathbf{Q}}_h\in V_h$${\mathbf{Q}}_h\in V_h$Q_(h)inV_(h)\mathbf{Q}_{h} \in V_{h}${\mathbf{Q}}_h\in V_h$ be given. Find $(c_h^n,{\mathbf{q}}_h^n)\in W_h\times V_h$$\left(c_h^n,{\mathbf{q}}_h^n\right)\in W_h\times V_h$(c_(h)^(n),q_(h)^(n))inW_(h)xxV_(h)\left(c_{h}^{n}, \mathbf{q}_{h}^{n}\right) \in W_{h} \times V_{h}$(c_h^n,{\mathbf{q}}_h^n)\in W_h\times V_h$ such that for all $w_h\in W_h$$w_h\in W_h$w_(h)inW_(h)w_{h} \in W_{h}$w_h\in W_h$ there holds

and for all ${\mathbf{v}}_h\in V_h$${\mathbf{v}}_h\in V_h$v_(h)inV_(h)\mathbf{v}_{h} \in V_{h}${\mathbf{v}}_h\in V_h$ we have

Remark 2.4 For details regarding the implementation of the scheme above we refer to [5]. The convergence of the scheme was discussed in [6] in a very general framework including also the case of transport with equilibrium (nonlinear) sorption. The existence and uniqueness of a solution of Problem 2.3 is also treated there.

The algorithm based on Problems 2.1 and 2.3 was implemented in the software package $UG$$UG$UGU G$UG$ [1].

Heterogeneous aquifer systems are described by a statistically homogeneous normally distributed log-hydraulic conductivity $U(\mathbf{x})=\log (K(\mathbf{x})),\mathbf{x}\in {\mathbb{R}}^d$$U(\mathbf{x})=\log (K(\mathbf{x})),\mathbf{x}\in {\mathbb{R}}^d$U(x)=log(K(x)),xinR^(d)U(\mathbf{x})=\log (K(\mathbf{x})), \mathbf{x} \in \mathbb{R}^{d}$U(\mathbf{x})=\log (K(\mathbf{x})),\mathbf{x}\in {\mathbb{R}}^d$. The fluctuations random field $u(\mathbf{x})=U(\mathbf{x})-m$$u(\mathbf{x})=U(\mathbf{x})-m$u(x)=U(x)-mu(\mathbf{x})=U(\mathbf{x})-m$u(\mathbf{x})=U(\mathbf{x})-m$, where $m=E(U(\mathbf{x}))$$m=E(U(\mathbf{x}))$m=E(U(x))m=E(U(\mathbf{x}))$m=E(U(\mathbf{x}))$ is the constant expectation of $U$$U$UU$U$, is homogeneous in a wide sense (or second order homogeneous) with mean zero and a correlation function which depends only on the difference vector $\mathit{\rho }={\mathbf{x}}_1-{\mathbf{x}}_2$$\mathit{\rho }={\mathbf{x}}_1-{\mathbf{x}}_2$rho=x_(1)-x_(2)\boldsymbol{\rho}=\mathbf{x}_{1}-\mathbf{x}_{2}$\mathit{\rho }={\mathbf{x}}_1-{\mathbf{x}}_2$, i.e. $E\left(u\left({\mathbf{x}}_1\right)u\left({\mathbf{x}}_2\right)\right)=C({\mathbf{x}}_1,{\mathbf{x}}_2)=C(\mathit{\rho })$$E\left(u\left({\mathbf{x}}_1\right)u\left({\mathbf{x}}_2\right)\right)=C\left({\mathbf{x}}_1,{\mathbf{x}}_2\right)=C(\mathit{\rho })$E(u(x_(1))u(x_(2)))=C(x_(1),x_(2))=C(rho)E\left(u\left(\mathbf{x}_{1}\right) u\left(\mathbf{x}_{2}\right)\right)=C\left(\mathbf{x}_{1}, \mathbf{x}_{2}\right)=C(\boldsymbol{\rho})$E\left(u\left({\mathbf{x}}_1\right)u\left({\mathbf{x}}_2\right)\right)=C({\mathbf{x}}_1,{\mathbf{x}}_2)=C(\mathit{\rho })$. Homogeneous random fields have a spectral representation given by a stochastic Fourier integral

$Z(d\mathbf{k})$$Z(d\mathbf{k})$Z(dk)Z(d \mathbf{k})$Z(d\mathbf{k})$ is a complex random measure on ${\mathbb{R}}^d$${\mathbb{R}}^d$R^(d)\mathbb{R}^{d}${\mathbb{R}}^d$ with the property $E(|Z(d\mathbf{k}){|}^2)=F(d\mathbf{k})$$E\left(|Z(d\mathbf{k}){|}^2\right)=F(d\mathbf{k})$E(|Z(dk)|^(2))=F(dk)E\left(|Z(d \mathbf{k})|^{2}\right)=F(d \mathbf{k})$E(|Z(d\mathbf{k}){|}^2)=F(d\mathbf{k})$, where $F$$F$FF$F$ is the spectral distribution function of the homogeneous random field. $F$$F$FF$F$ is defined as the inverse Fourier transform of the correlation function

so that the spectral representation (9) preserves the correlation (10) [10]. In case of Gaussian homogeneous random fields $Z(d\mathbf{k})$$Z(d\mathbf{k})$Z(dk)Z(d \mathbf{k})$Z(d\mathbf{k})$ is a complex white noise random measure [2].

The standard Fourier method to generate homogeneous random fields consists of numerical evaluations of the spectral representation (9). The evaluation through a Monte Carlo integration approach is called "randomization method" [7]. The simplest variant of the randomization method has the form of a superposition of $N_p$$N_p$N_(p)N_{p}$N_p$ random sine modes,

where ${\zeta }_j,{\eta }_j,j=1,\dots ,N_p$${\zeta }_j,{\eta }_j,j=1,\dots ,N_p$zeta_(j),eta_(j),j=1,dots,N_(p)\zeta_{j}, \eta_{j}, j=1, \ldots, N_{p}${\zeta }_j,{\eta }_j,j=1,\dots ,N_p$ are independent Gaussian random variables of mean zero and unit variance, ${\sigma }^2={\int }_{{\mathbb{R}}^d}F(d\mathbf{k})$${\sigma }^2={\int }_{{\mathbb{R}}^d} F(d\mathbf{k})$sigma^(2)=int_(R^(d))F(dk)\sigma^{2}=\int_{\mathbb{R}^{d}} F(d \mathbf{k})${\sigma }^2={\int }_{{\mathbb{R}}^d}F(d\mathbf{k})$ is, according to (10), the variance of $u$$u$uu$u$, and ${\mathbf{k}}_j$${\mathbf{k}}_j$k_(j)\mathbf{k}_{j}${\mathbf{k}}_j$ are independent random vectors which are also independent of ${\zeta }_j$${\zeta }_j$zeta_(j)\zeta_{j}${\zeta }_j$ and ${\eta }_j$${\eta }_j$eta_(j)\eta_{j}${\eta }_j$, sampled according to a probability density $p(\mathbf{k})=2F(\mathbf{k})/{\sigma }^2$$p(\mathbf{k})=2F(\mathbf{k})/{\sigma }^2$p(k)=2F(k)//sigma^(2)p(\mathbf{k})=2 F(\mathbf{k}) / \sigma^{2}$p(\mathbf{k})=2F(\mathbf{k})/{\sigma }^2$. The accuracy of the method is influenced by the sampling of the wave numbers (components of ${\mathbf{k}}_j$${\mathbf{k}}_j$k_(j)\mathbf{k}_{j}${\mathbf{k}}_j$ ). Improvements are achieved if one divides the range of wave numbers into a finite number of bins, of uniform or logarithmic distributed dimensions, and one samples randomly wave numbers in each bin [2]. Since $U(\mathbf{x})=\log (K(\mathbf{x}))$$U(\mathbf{x})=\log (K(\mathbf{x}))$U(x)=log(K(x))U(\mathbf{x})=\log (K(\mathbf{x}))$U(\mathbf{x})=\log (K(\mathbf{x}))$ is normally distributed, the mean and the variance of the hydraulic conductivity can be finally computed from the relations $E(K)=e^{m+{\sigma }^2/2}$$E(K)=e^{m+{\sigma }^2/2}$E(K)=e^(m+sigma^(2)//2)E(K)=e^{m+\sigma^{2} / 2}$E(K)=e^{m+{\sigma }^2/2}$ and $\mathrm{var}(K)=e^{2m+{\sigma }^2}(e^{{\sigma }^2}-1)$$\mathrm{var}(K)=e^{2m+{\sigma }^2}\left(e^{{\sigma }^2}-1\right)$var(K)=e^(2m+sigma^(2))(e^(sigma^(2))-1)\operatorname{var}(K)=e^{2 m+\sigma^{2}}\left(e^{\sigma^{2}}-1\right)$\mathrm{var}(K)=e^{2m+{\sigma }^2}(e^{{\sigma }^2}-1)$.

A homogeneous random field is "ergodic" if spatial average unbiased estimators converge to ensemble averages. For instance, the random function given by the space average $\overline{u}=\frac{1}{\mathcal{V}(\mathrm{\Omega })}{\int }_{\mathrm{\Omega }}u(\mathbf{x})d\mathbf{x}$$\overline{u}=\frac{1}{\mathcal{V}(\mathrm{\Omega })}{\int }_{\mathrm{\Omega }} u(\mathbf{x})d\mathbf{x}$bar(u)=(1)/(V(Omega))int_(Omega)u(x)dx\bar{u}=\frac{1}{\mathcal{V}(\Omega)} \int_{\Omega} u(\mathbf{x}) d \mathbf{x}$\overline{u}=\frac{1}{\mathcal{V}(\mathrm{\Omega })}{\int }_{\mathrm{\Omega }}u(\mathbf{x})d\mathbf{x}$, where $\mathcal{V}(\mathrm{\Omega })$$\mathcal{V}(\mathrm{\Omega })$V(Omega)\mathcal{V}(\Omega)$\mathcal{V}(\mathrm{\Omega })$ is the volume of the domain $\mathrm{\Omega }\subset {\mathbb{R}}^d$$\mathrm{\Omega }\subset {\mathbb{R}}^d$Omega subR^(d)\Omega \subset \mathbb{R}^{d}$\mathrm{\Omega }\subset {\mathbb{R}}^d$, is an unbiased estimator of the expectation $E(\overline{u})=0$$E(\overline{u})=0$E( bar(u))=0E(\bar{u})=0$E(\overline{u})=0$. Then, the field $u$$u$uu$u$ is called ergodic in the mean if

The mean-square convergence of $\overline{u}$$\overline{u}$bar(u)\bar{u}$\overline{u}$ to the limit (12) is tantamount with the convergence to zero of the variance

Hence, a finite correlation range ${\lambda }_u={\int }_{\mathrm{\Omega }}C(\mathit{\rho })d\mathit{\rho }/C(\mathbf{0})$${\lambda }_u={\int }_{\mathrm{\Omega }} C(\mathit{\rho })d\mathit{\rho }/C(\mathbf{0})$lambda_(u)=int_(Omega)C(rho)d rho//C(0)\lambda_{u}=\int_{\Omega} C(\boldsymbol{\rho}) d \boldsymbol{\rho} / C(\mathbf{0})${\lambda }_u={\int }_{\mathrm{\Omega }}C(\mathit{\rho })d\mathit{\rho }/C(\mathbf{0})$ of the random field is a sufficient condition for ergodicity: $\mathrm{var}(\overline{u})={\lambda }_{u}C(\mathbf{0})/\mathcal{V}(\mathrm{\Omega })$$\mathrm{var}(\overline{u})={\lambda }_{u}C(\mathbf{0})/\mathcal{V}(\mathrm{\Omega })$var( bar(u))=lambda_(u)C(0)//V(Omega)\operatorname{var}(\bar{u})= \lambda_{u} C(\mathbf{0}) / \mathcal{V}(\Omega)$\mathrm{var}(\overline{u})={\lambda }_{u}C(\mathbf{0})/\mathcal{V}(\mathrm{\Omega })$ goes to zero for large $\mathcal{V}(\mathrm{\Omega })$$\mathcal{V}(\mathrm{\Omega })$V(Omega)\mathcal{V}(\Omega)$\mathcal{V}(\mathrm{\Omega })$ and the limit (12) is finite. Moreover, because $u(\mathbf{x})$$u(\mathbf{x})$u(x)u(\mathbf{x})$u(\mathbf{x})$ is a homogeneous Gaussian random function, a finite correlation range also ensures the ergodicity of the higher order moments and, in particular, of the correlation function $C(\rho )$$C(\rho )$C(rho)C(\rho)$C(\rho )$ [10].

Under routinely made assumptions, Darcy velocity fields and Lagrangian observables of the process inherit the ergodicity of the $u$$u$uu$u$ field [9]. It is therefore imperiously necessary that numerically generated log-hydraulic conductivity fields satisfactorily reproduce the ergodic behavior assumed for $u$$u$uu$u$.

In this section we present a numerical test of our method. We consider to solve the equations (1)-(2) and (3) in the domain $\mathrm{\Omega }=[0,300]\times [0,250]$$\mathrm{\Omega }=[0,300]\times [0,250]$Omega=[0,300]xx[0,250]\Omega=[0,300] \times[0,250]$\mathrm{\Omega }=[0,300]\times [0,250]$. The final time is $T=500$$T=500$T=500T=500$T=500$. The boundary conditions for the flow problem are: ${\psi }_{\mid x=0}=5,{\psi }_{\mid x=300}=0$${\psi }_{\mid x=0}=5,{\psi }_{\mid x=300}=0$psi_(∣x=0)=5,psi_(∣x=300)=0\psi_{\mid x=0}=5, \psi_{\mid x=300}=0${\psi }_{\mid x=0}=5,{\psi }_{\mid x=300}=0$ and everywhere else the flux is null. The initial position of the contaminant plume is $[30,80]\times [100,150]$$[30,80]\times [100,150]$[30,80]xx[100,150][30,80] \times[100,150]$[30,80]\times [100,150]$, where a constant concentration is imposed such that the total mass of solute is 1 . We assume everywhere homogeneous Neumann boundary conditions for the transport equation. The diffusion-dispersion coefficient is taken $D=0.01$$D=0.01$D=0.01D=0.01$D=0.01$ and we consider to have no reaction $(r(c)=0)$$(r(c)=0)$(r(c)=0)(r(c)=0)$(r(c)=0)$. All the dimensions should be understood in meters and days. The samples of the random conductivity field are generated with the method presented in Section 3. We used $N_p=6400$$N_p=6400$N_(p)=6400N_{p}=6400$N_p=6400$ periodic modes in (11) and uniform sampling of the wave numbers to generate a field $U$$U$UU$U$ with correlation length equal to 1 , mean $m=1.66$$m=1.66$m=1.66m=1.66$m=1.66$ and variance ${\sigma }^2=0.1$${\sigma }^2=0.1$sigma^(2)=0.1\sigma^{2}=0.1${\sigma }^2=0.1$ (which yield $E(K)=5$$E(K)=5$E(K)=5E(K)=5$E(K)=5$ and $\mathrm{var}(K)=3.22$$\mathrm{var}(K)=3.22$var(K)=3.22\operatorname{var}(K)=3.22$\mathrm{var}(K)=3.22$ ). The MFEM - Euler implicit scheme from Section 2 is used to compute numerical solutions of the flow and transport problems.

Figure 1 presents a sample of a Gaussian correlated hydraulic conductivity random field $K$$K$KK$K$ generated with (11). Figure 2 shows the correlation in the mean flow direction of the random field $u$$u$uu$u$, estimated by spatial averages over a rectangular domain $[0,4000]\times [0,100]$$[0,4000]\times [0,100]$[0,4000]xx[0,100][0,4000] \times[0,100]$[0,4000]\times [0,100]$. The agreement of spatial estimations with the analytical correlation functions indicate fairly good ergodic properties of the randomization method, which renders it suitable for large scale simulations of transport in heterogeneous

aquifers. The resulting concentration field is presented in Figure 3. The space average of the simulated velocity components shown in Figure 4 are close to the nominal values ( 0.083 and 0.0 respectively), which indicates the ergodicity in the mean of the simulated velocity field. Also shown in Figure 4 are the Lagrangian mean velocity $⟨\mathbf{V}⟩={\int }_{\mathrm{\Omega }}\mathbf{V}(\mathbf{x})c(\mathbf{x},t)d\mathbf{x}$$⟨\mathbf{V}⟩={\int }_{\mathrm{\Omega }} \mathbf{V}(\mathbf{x})c(\mathbf{x},t)d\mathbf{x}$(:V:)=int_(Omega)V(x)c(x,t)dx\langle\mathbf{V}\rangle=\int_{\Omega} \mathbf{V}(\mathbf{x}) c(\mathbf{x}, t) d \mathbf{x}$⟨\mathbf{V}⟩={\int }_{\mathrm{\Omega }}\mathbf{V}(\mathbf{x})c(\mathbf{x},t)d\mathbf{x}$ and the slope of the components of the first spatial concentration moment $\frac{1}{t}{\int }_{\mathrm{\Omega }}\mathbf{x}[c(\mathbf{x},t)-c(\mathbf{x},0)]d\mathbf{x}$$\frac{1}{t}{\int }_{\mathrm{\Omega }} \mathbf{x}[c(\mathbf{x},t)-c(\mathbf{x},0)]d\mathbf{x}$(1)/(t)int_(Omega)x[c(x,t)-c(x,0)]dx\frac{1}{t} \int_{\Omega} \mathbf{x}[c(\mathbf{x}, t)-c(\mathbf{x}, 0)] d \mathbf{x}$\frac{1}{t}{\int }_{\mathrm{\Omega }}\mathbf{x}[c(\mathbf{x},t)-c(\mathbf{x},0)]d\mathbf{x}$ which, in their turn, indicate the ergodicity of the Lagrangian observables [9].

These first numerical results show that a Monte Carlo approach based on MFEM and randomization method fulfil the basic ergodicity requirements for modeling transport in aquifer systems with spatial heterogeneity described by statistically homogeneous log-hydraulic conductivity fields. Results obtained with the method proposed here, which yields exact flow solutions, will be compared with Monte Carlo simulations for the same problem based on first order approximations of the velocity field [8]. The perspective to have accurate error estimations of single realizations MFEM solutions is expected to improve the assessment of the overall precision of the method. Further work will be done to investigate the probability distribution of the pressure head, velocity, and concentration, as well as the ergodicity with respect to reference stochastic models of the transport process [8]. Statistical analyses of the first two spatial moments of the concentration fields resulted from this MFEM Monte Carlo approach will also be done to get new insights on the memory effects due to the dependence on initial conditions which were previously investigated for approximated velocity fields [9].

[1] P. Bastian, K. Birken, K. Johanssen, S. Lang, N. Neuss, H. Rentz-Reichert and C. Wieners, UG-a flexible toolbox for solving partial differential equations, Comput. Visualiz. Sci. 1 (1997), pp. 27-40.
[2] P. R. Kramer, O. Kurbanmuradov and K. Sabelfeld, Comparative analysis of multiscale Gaussian random field simulation algorithms, J. Comp. Phys. 226, pp. 897-924 (2007).
[3] F. A. Radu, I. S. Pop and P. Knabner, Order of convergence estimates for an Euler implicit, mixed finite element discretization of Richards' equation, SIAM J. Numer. Anal. 42, pp. 1452-1478 (2004).
[4] F. A. Radu, I. S. Pop and P. Knabner, Error estimates for a mixed finite element discretization of some degenerate parabolic equations, Numer. Math. 109 (2), pp. 285-311 (2008).
[5] F. A. Radu, M. Bause, A. Prechtel and S. Attinger, A mixed hybrid finite element discretization scheme for reactive transport in porous media, Numerical Mathematics and Advanced Applications, K. Kunisch, G. Of and O. Steinbach (editors), Springer Verlag, Heidelberg (2008), pp. 513-520.
[6] F. A. Radu, I. S. PoP and S. ATTINGER, Analysis of an Euler implicit - mixed finite element scheme for reactive solute transport in porous media, Numerical Methods Partial Differential Equations, DOI:10.1002/num. 20436 (2009).
[7] K. Sabelfeld, Monte Carlo Methods in Boundary Value Problems, Springer, Berlin (1991).
[8] N. Suciu, C. Vamoş, J. Vanderborght, H. Hardelauf, and H. Vereecken, Numerical investigations on ergodicity of solute transport in heterogeneous aquifers, Water Resour. Res. 42, W04409, DOI:10.1029/2007WR004546 (2006).
[9] N. Suciu, C. Vamoş, H. Vereecken, K. Sabelfeld and P. Knabner, Memory effects induced by dependence on initial conditions and ergodicity of transport in heterogeneous media, Water Resour. Res. 44, W08501, DOI:10.1029/2007WR006740 (2008).
[10] A. M. Yaglom, Correlation Theory of Stationary and Related Random Functions, Volume I: Basic Results, Springer-Verlag, New York (1987).