# Numerical Analysis of the Advection-Diffusion of a Solute in Porous Media with Uncertainty

Julia Charrier

► To cite this version:

Julia Charrier. Numerical Analysis of the Advection-Diffusion of a Solute in Porous Media with Uncertainty. SIAM/ASA Journal on Uncertainty Quantification, 2015, 3 (1), 10.1137/130937457. hal-00862960

HAL Id: hal-00862960

https://hal.science/hal-00862960

Submitted on 17 Sep 2013

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Numerical analysis of the advection-diffusion of a solute in porous media with uncertainty

Julia Charrier*

September 17, 2013

## Abstract

We consider the problem of numerically approximating the solution of the coupling of the flow equation in a porous medium, with the advection-diffusion equation in the presence of uncertainty on the permeability of the medium. Random coefficients are classically used in the flow equation to modelize uncertainty. More precisely, we propose the numerical analysis of a method developed to compute the mean value of the spread of a solute introduced at the initial time, and the mean value of the macro-dispersion, defined as the temporal derivative of the spread. We consider a Monte-Carlo method to deal with the uncertainty, i.e. with the randomness of the permeability field. The flow equation is solved using a finite element method. The advection-diffusion equation is seen as a Fokker-Planck equation, and its solution is hence approximated thanks to a probabilistic particular method. The spread is indeed the expected value of a function of the solution of the corresponding stochastic differential equation, and is computed using an Euler scheme for the stochastic differential equation and a Monte-Carlo method. Error estimates on quantities generalizing the mean spread and the mean macro-dispersion are established, under some assumptions including the case of random fields of lognormal type (i.e. neither uniformly bounded from above nor below with respect to the random parameter) with low regularity, which is pertinent on an application point of view and rises several mathematical difficulties.

**Keywords:** uncertainty quantification, elliptic PDE with random coefficients, advection-diffusion equation, probabilistic interpretation of PDE, Monte-Carlo method, Euler scheme for SDE.

## 1 Introduction

Numerical modeling is an important key for the management and remediation of groundwater resources. The heterogeneity of natural geological formations has a major impact in the contamination of groundwater by migration of pollutants. In order to account for the limited knowledge of the geological characteristics and for the natural heterogeneity, stochastic models have been developed, see e.g. [7],[8]. The permeability of the porous media is then a random field. Our aim is then to study the migration of a contaminant in steady flow. The flow velocity is computed by solving an elliptic partial differential equation with random coefficients. The solute concentration is then the solution of an advection-diffusion equation, where the flow velocity, which is a random field, appears as a coefficient. The quantities we are interested in are finally the mean value of the spread of the solute, that is to say the mean value of the spatial variance of the solute, and the mean value of the dispersion, which is defined as the derivative of the spread with respect to the time. The determination of the large-scale dispersion coefficients has been widely debated in the last twenty five years, see e.g [3], [5], [9], [13], [25], [28] and [27].

Here we are interested in the case of a lognormal permeability field, which is a widely used model. Moreover we consider the case, physically pertinent, where the correlation length is small and the uncertainty

<footnote>*LATP, Aix-Marseille Université, CMI, 39 rue Joliot Curie, 13453 Marseille cedex 13, France. Email: julia.charrier@univ-amu.fr</footnote>
important. Therefore methods based on the approximation of the coefficients in a finite dimensional stochastic space, such as stochastic galerkin methods and stochastic collocation method would be highly expensive, and hence do not seem to be suitable to deal with such cases. Neither seem perturbation methods, since we suppose the uncertainty to be important. As regards the advection-diffusion equation, we focus on the advection-dominated model. Therefore it was chosen not to consider an Eulerian method, in order to avoid numerical diffusion. The method described and analyzed below is, up to a few modifications, the method proposed and implemented by A.Beaudoin, J.R. de Dreuzy and J.Erhel to compute the mean macro-dispersion in 2D, their numerical results can be found e.g. in [6]. A Monte-Carlo method is used to deal with the uncertainty. The solution of the steady flow equation is computed by using a finite element approximation. The solution of the advection-diffusion equation is approximated using a probabilistic method. We consider the stochastic differential equation associated to this Fokker Planck equation and its solution is approximated with an Euler scheme. A Monte-Carlo method provides finally an approximation of the solution of the Fokker Planck equation. All these steps together lead to an approximation of the mean spread. The mean macro-dispersion is then approximated by the same way (after applying Itô formula).

The aim of this paper is to provide the numerical analysis of the above described method, in particular
this works completes the paper [6] which proposes detailed numerical results. therefore we do not propose
numerical experiments in this work, which would be redundant with [6], but rather give a theoritical jus-
tification for the use of this method, and more important, study the speed of convergence with respect to
each discretization parameter in order to guide the choice of these parameters and hence to be able to use
efficiently this method.

More precisely we furnish a priori error estimates for the approximations of quantities generalizing the mean spread and the mean macro-dispersion, the spread and the macro-dispersion being the lead because of their physical interest. A specificity of this work is to address the coupling of the flow equation with the advection-diffusion equation, whereas most of the existing numerical analysis of methods for uncertainty quantification are limited to the flow equation, see e.g [1], [2]. A particularity of this work is also the use of numerical analysis tools from two different areas: finite element method and weak error analysis for SDEs. More importantly, we emphasize that in this work we deal with random permeability fields which are neither uniformly bounded from above nor below (with respect to the random parameter), moreover we suppose only low spacial regularity (Hölder regularity), whereas most of the works proposing analysis of numerical methods for flow equation in porous media with uncertainty suppose the permeability field to be smooth and uniformly bounded from above and below, which simplifies drastically the numerical analysis. Therefore the numerical analysis proposed here cannot be obtain by simply combining classical results. More precisely the first main difficulty remains in getting sharp explicit dependence on the random parameter $\omega$ (modelizing uncertainty) at each step, to take into account the fact that we cannot have uniform estimates with respect to the random parameter $\omega$. The second main difficulty is to deal with the discretization of SDE with non lipschitz drift and to include the error of spatial discretization in the wear error of time discretization.

After presenting the physical model in section 2, we describe in detail in section 3 the numerical method mentioned above. Section 4 is devoted to the numerical analysis of this method with slightly modified assumptions. More precisely we consider the case of a random permeability field, supposed to be almost surely periodic, with some integrability properties with respect to the random parameter and under regularity assumptions (Hölder continuous with respect to the spatial variable). We first give preliminary results. The first one is a bound of the finite element error in $W^{1,\infty}$ norm in the low regularity case. The bound has to be explicit with respect to $\omega$ in order to be integrated later with respect to $\omega$. The second one is a weak error result for the Euler scheme on a stochastic differential equation with additive noise and a $C^{0,\alpha}$ or $C^{1,\alpha}$ drift, which also takes into account the spatial discretization. Once again we need to track the dependence on $\omega$ sharply. After these preliminary results, we give the two main results of this paper, namely error results on the mean generalized spread and on its time derivative, the mean generalized macro-dispersion.

## 2 Physical model

We recall here, up to some minor modifications, the physical model studied in [6].
## 2.1 Steady flow equation

We consider an isotropic porous medium, we suppose the porosity to be constant, equal to 1. The domain $O$ is a box included in $\mathbb{R}^d$, with $d = 1, 2$, or $3$. The heterogeneity of the natural geological formations and the lack of data lead us to use a stochastic model, see e.g. [7],[8]. A classical case is to take an homogeneous lognormal field to modelize the permeability field:

$$a(\omega, x) = e^{g(\omega, x)}, \quad x \in O, \omega \in \Omega,$$

where $g$ is a gaussian field characterized by its mean $m$ and its covariance function, we suppose that the covariance function is of the form:

$$\operatorname{cov}[g](x, y) = \sigma^2 \exp \left( - \frac{\|x - y\|^\delta}{\ell} \right), \qquad (2.1)$$

for some $\delta > 0$. The random parameter is denoted by $\omega$. The case of an exponential covariance function corresponds to $\delta = 1$, and furnishes a model which is a reasonably good fit to some field data, see e.g [12] and [16].

The variance of the log hydraulic conductivity $\sigma^2$ is typically in the interval [1, 10], the correlation length $\ell$ typically ranges between 0.1m and 100m, whereas the size of the domain has to be at least hundred times the correlation length $l$. Classical laws governing the steady flow in porous media without source are mass conservation $\operatorname{div}(v) = 0$ and Darcy law $v = -a\nabla p$, where $v$ is the Darcy velocity and $p$ the hydraulic head. Finally, the hydraulic head is the solution of the following elliptic PDE with a random coefficient : for almost all $\omega$

$$\operatorname{div}(a(\omega, x)\nabla p(\omega, x)) = 0, \quad x \in O. \qquad (2.2)$$

This equation has to be subjected to boundary conditions (mixed boundary conditions for example), and the boundary condition is then imposed for almost all $\omega \in \Omega$.

Here, $\omega$ is then the parameter describing the randomness of the media. We recall that the Darcy velocity is then defined by

$$v(\omega, x) = -a(\omega, x)\nabla p(\omega, x).$$

## 2.2 Advection-diffusion equation

An inert solute is injected in the porous medium and transported by advection and diffusion. Here we consider only molecular diffusion, assumed to be homogeneous and isotropic. This type of solute migration is described by the advection-diffusion equation:

$$\frac{\partial c(\omega, x, t)}{\partial t} + v(\omega, x) \cdot \nabla_x c(\omega, x, t) - D \Delta_x c(\omega, x, t) = 0, \qquad (2.3)$$

where $D > 0$ is the molecular diffusion coefficient, $v$ the Darcy velocity defined previously and $c$ the solute concentration. We consider the case of advection-dominated model, i.e. the case where the Peclet number $Pe = \frac{\ell\|v\|_{\text{mean}}}{D}$ is large (typically $\ge 100$). The initial condition at $t = 0$ is the injection of the solute, i.e. for example $c(t = 0) = \frac{\ell}{|R|}$ where $R$ is a box included in $O$ (of volume $|R|$). Equation (2.3) has to be supplemented with boundary conditions on $\partial O$.

## 2.3 Spread and macro-dispersion

We now define the two quantities we want to compute.

First we introduce the center of mass of the solute distribution :

$$G(\omega, t) = \int_{O} c(\omega, x, t)xdx.$$
Our aim is then to compute $S(t)$ the mean spread of mass around $G$, and the mean macro-dispersion $\mathcal{D}(t)$, defined as its time derivative, i.e.

$$S(\omega, t) = \int_{O} c(\omega, x, t)(x - G(\omega, t))(x - G(\omega, t))^{t} dx, \quad S(t) = \mathbb{E}_{\omega}[S(\omega, t)]$$

and

$$\mathcal{D}(\omega, t) = \frac{1}{2} \frac{dS(\omega, t)}{dt}, \quad \mathcal{D}(t) = \mathbb{E}_{\omega}[\mathcal{D}(\omega, t)].$$

**Remark 2.1.** Note that the real quantities of interest are the mean spread or macro-dispersion in each direction, which are scalar quantities: for $i = 1, ..., d$

$$S_i(\omega, t) = \int_{O} c(\omega, x, t)(x_i - G_i(\omega, t))^2 dx, \quad S_i(t) = \mathbb{E}_{\omega}[S_i(\omega, t)]$$

and

$$\mathcal{D}_i(\omega, t) = \frac{1}{2} \frac{dS_i(\omega, t)}{dt}, \quad \mathcal{D}_i(t) = \mathbb{E}_{\omega}[\mathcal{D}_i(\omega, t)],$$

where

$$G_i(\omega, t) = \int_{O} c(\omega, x, t)x_i dx,$$

for example the transversal and longitudinal spread and macro-dispersion in the case $d = 2$. But these quantities correspond to the diagonal coefficients of the matrices $S$ and $D$ defined above, that is why for simplicity we will consider the matricial definitions of the spread and macro-dispersion $S$ and $D$ given above.

# 3 Description of the numerical method

## 3.1 A Monte-Carlo method to deal with uncertainty

As precised above, we suppose the uncertainty to be large, typically $\sigma^2 \in [1, 10]$, therefore perturbation type methods [3], [9], [13], [24], [25], [26], [27] do not seem to be suitable. Moreover, since we suppose $\ell$ to be small, $\sigma^2$ to be large and $cov[g]$ to be only lipschitz, stochastic galerkin and stochastic collocation methods (see e.g. [1], [2], [14], [10] and the references therein) do not seem to be adapted. Indeed, in such cases, the permeability field $a$ cannot be approximated correctly with a reasonable number of random variables. In particular, the eigenvalues of the Karhunen-Loève development are explicit in this case if we choose the 1-norm (see [11]), and we know that the number of term in the truncated Karhunen-Loève development should be much bigger than 100 (see [4] for example), which is not possible on a practical point of view. Therefore we choose to use a Monte-Carlo method to deal with uncertainty. More precisely, we consider $N$ independent realizations of the permeability field $a(x, \omega_1), ...a(x, \omega_N)$. For each $i$ from 1 to $N$, we compute approximations of the spread $S^i(t)$ and of the macro-dispersion $\mathcal{D}^i(t)$ corresponding to the permeability field $a^i$ as specified below, and we approximate the mean spread $S(t)$ by $\frac{1}{N} \sum_{i=1}^{N} S^i(t)$ and the mean macro-dispersion $\mathcal{D}(t)$ by $\frac{1}{N} \sum_{i=1}^{N} \mathcal{D}^i(t)$. For simplicity, the index $i$ as well as the random variable $\omega$ will be omitted in the remainder of this section, which is devoted to the description of the numerical method used to compute the solution of a deterministic problem: the computation of spread and dispersion.

## 3.2 Approximation of the flow velocity

The hydraulic head is defined as the solution of the following elliptic partial differential equation:

$$\mathrm{div}(a(x)\nabla p(x)) = 0, \quad x \in O,$$

submitted to boundary conditions.

We define then an approximation $p_h$ of $p$ in a finite elements space of continuous piecewise linear functions, with maximum space mesh $h$. The velocity $v$ is then approximated by $v_h(x) = -a(x)\nabla p_h(x)$.
### 3.3 A probabilistic particular method

The solute concentration is defined as the solution of (2.3). The domain $O$ is chosen such that a very small amount of the solute reaches the boundary. Therefore, in practice, it is harmless to replace (2.3) by:

$$ \left\{ \begin{array}{ll} \frac{\partial c}{\partial t}(x,t) + v(x) \cdot \nabla c(x,t) - D \Delta c(x,t) & = 0, & x \in \mathbb{R}^d \text{ and } t \in [0, T] \\ c(x,0) & = c_0(x), & x \in \mathbb{R}^d, \end{array} \right. $$

where $v$ is extended to $\mathbb{R}^d$ in some way (see Section 4 for more details). Since $\mathrm{div}(v) = 0$, this is a Fokker-Planck equation. A probabilistic particular method was chosen to approximate the solution of this Fokker-Planck equation. This is motivated among others to avoid numerical diffusion, since we focus of the advection dominated case. We then define the associated stochastic differential equation:

$$ \left\{ \begin{array}{ll} dX(t) & = v(X(t))dt + \sqrt{2D}dW(t) \\ X(0) & = X_0, \end{array} \right. $$

where $X_0$ is a random variable with density $c_0$ with respect to the Lebesgue measure. It is classical that $X(t)$ admits then $c(x,t)dx$ as density. The law of $X$ can be approximated by a Monte-Carlo method. We take $M$ independent realizations of approximations of $X$ using an Euler scheme and the approximated flow velocity $v_h: X_{n,h}^1, \dots, X_{n,h}^M$.

$$ \left\{ \begin{array}{ll} X_{n,h}^j(t_{k+1}) & = X_{n,h}^j(t_{k+1}) + v_h(X_{n,h}^j(t_k))\Delta t + \sqrt{2D\Delta t}N_k^j \text{ for } t \in [t_k, t_{k+1}], \\ X_{n,h}^j(0) & = X_0, \end{array} \right. $$

where $T = n\Delta t$, $t_k = k\Delta t$ and the $N_k^j$ are independent $d$-dimensional mean-free gaussian random vector with identity as covariance. Finally we approximate the mass center $G(t)$ by $G_{n,h}^M(t) = \frac{1}{M} \sum_{j=1}^M X_{n,h}^j(t)$, the spread $S(t)$ by

$$ S_{n,h}^M(t) = \frac{1}{M} \sum_{j=1}^{M} (X_{n,h}^j(t) - G_{n,h}^M(t))(X_{n,h}^j(t) - G_{n,h}^M(t))^t. $$

Indeed we have $G(t) = \mathbb{E}[X(t)]$, $S(t) = \mathbb{E}[(X(t) - G(t))(X(t) - G(t))^t]$. Moreover the macro-dispersion is defined by $\mathcal{D}(t) = \frac{1}{2}\frac{d}{dt}S(t)$, which thanks to Itô formula, is equal to

$$ \frac{1}{2}\mathbb{E}[(X(t) - G(t))(v(X(t)) - V(t))^t + (v(X(t)) - V(t))(X(t) - G(t))^t] + DId, $$

where we have used the notation $V(t) = \mathbb{E}[v(X(t))]$. We approximate hence the macro-dispersion $\mathcal{D}(t)$ by

$$ DId + \frac{1}{2M} \sum_{j=1}^{M} (X_{n,h}^j(t) - G_{n,h}^M(t))(v(X_{n,h}^j(t)) - V_{n,h}^M(t))^t + (v(X_{n,h}^j(t)) - V_{n,h}^M(t))(X_{n,h}^j(t) - G_{n,h}^M(t))^t, $$

where we have defined $V_{n,h}^M(t) = \frac{1}{M} \sum_{j=1}^M v(X_{n,h}^j(t))$.

For more details on a possible numerical implementation, see [6]. Note however that in this case, the numerical method used to compute the mean value of the macro-dispersion is slightly different. The derivative is computed using the increase corresponding to a small time step. See section 4 for more details concerning the difference of the convergence rates obtained for each method.

# 4 Numerical analysis of the method

## 4.1 Notations and assumptions

We consider $O$ a box of $\mathbb{R}^d$, and $(\Omega, \mathcal{F}, \mathbb{P})$ a probability space. For $k \in \mathbb{N}$ and $0 < \alpha \le 1$, we denote by $C_b^k$ the space of functions which are $k$ times differentiable with bounded derivatives and $C_b^{k, \alpha}$ the space of $C_b^k$
functions such that any k-th derivative is α-hölder continuous. For $f \in C_b^{k,\alpha}$ we introduce the associated norm:

$$ \|f\|_{C_b^{k,\alpha}} = \max\{\|f\|_{\infty}, \|f'\|_{\infty}, \dots, \|f^{(k)}\|_{\infty}, |f^{(k)}|_{C^{0,\alpha}}\}, \text{ where the semi-norm of an } \alpha\text{-Hölder continuous function } g \text{ is defined by} $$

$$ |g|_{C^{0,\alpha}} = \max_{x \neq y} \frac{|g(x) - g(y)|}{|x - y|^{\alpha}}. $$

The numerical analysis of the above algorithm requires the solution of (2.2) to be sufficiently regular with respect to $x \in O$. Unfortunately, this is not the case with the above described model. Indeed we are dealing with an elliptic equation on a rectangular domain with mixed boundary conditions. This limitates the smoothness of the solution. Also, note that the advection-diffusion equation is set on the full space $\mathbb{R}^d$ and it is not clear how to extend the velocity field on $O$ to $\mathbb{R}^d$. We consider that this is a technical problem and avoid it by replacing the mixed boundary conditions in (2.2) by periodic boundary conditions, so that the solutions have the smoothness naturally associated to the smoothness of the permeability field, and the extension to $\mathbb{R}^d$ is trivial. Another way could be to truncate the velocity field close to the boundary of $O$ and to then extend it smoothly by zero on $\mathbb{R}^d$. The final solution would not be very different, since in practice the domain $O$ is chosen very large with respect to the box $R$ and a very small amount of the solute reaches the boundary. The solution of (2.2) with mixed boundary conditions being smooth inside the domain, the same analysis as below would give a similar result. We chose the periodic boundary conditions to simplify the presentation. We consider then the flow equation: for almost all $\omega$

$$ \left\{ \begin{array}{ll} \text{div}(a(\omega, x)\nabla p(\omega, x)) = f(x) & \text{on } \mathbb{R}^d, \\ \\ \int_O p(\omega, x)dx = 0, & \end{array} \right. \qquad (4.1) $$

and $p$ is $O$-periodic ($O$ being a box). The right hand side $f$ takes into account non homogeneous boundary conditions and is supposed to be $O$-periodic and such that $\int_O f(x)dx = 0$ (compatibility condition).

**Remark 4.1.** We could also consider a random second member $f$, which would provide similar results through a straightforward extension of the proofs below.

We now introduce two different sets of assumptions on the random coefficient $a$ and the second member $f$. The second assumption requires more spatial regularity than the first one. These two assumptions lead to different orders of convergence in the error estimate for the generalized spread. The error estimate for the generalized macro-dispersion will only be obtained under the strongest assumption 4.3.

**Assumption 4.2.** The permeability field is such that for some $0 < \alpha < 1$, we have for any finite $q \ge 1$, $a \in L^q(\Omega, C_b^{0,\alpha}(\mathbb{R}^d))$ and $\frac{1}{a_{min}} \in L^q(\Omega)$, where we have defined for almost all $\omega$, $a_{min}(\omega) = \min_{x \in O} a(\omega, x)$. We also suppose that for almost all $\omega$, $a(\omega, \cdot)$ is $O$-periodic. Moreover we suppose that $f \in L^r_{loc}(\mathbb{R}^d)$ for some $r > d$ such that $\alpha < 1 - \frac{d}{r}$ and is also $O$-periodic.

**Assumption 4.3.** The permeability field is such that for some $0 < \alpha < 1$, we have for any finite $q \ge 1$, $a \in L^q(\Omega, C_b^{1,\alpha}(\mathbb{R}^d))$ and $\frac{1}{a_{min}} \in L^q(\Omega)$, where we have defined for almost all $\omega$, $a_{min}(\omega) = \min_{x \in O} a(\omega, x)$. We also suppose that for almost all $\omega$, $a(\omega, \cdot)$ is $O$-periodic. Moreover we suppose that $f \in W_{loc}^{1,r}(\mathbb{R}^d)$ for some $r > d$ such that $\alpha < 1 - \frac{d}{r}$ and is also $O$-periodic.

Let us comment shortly the spatial regularity obtained for the realizations of the permeability field $a$ in the important example of an homogeneous lognormal field, that is to say $a = e^g$ with $\text{cov}[g](x, y) = k(x - y)$, where $k$ only depend of the norm of its argument. If the function $k$ is lipschitz continuous, then Assumption 4.2 is fullfiled for any $\alpha < 1/2$ (except of course the periodicity property, which is on a pratical point of view quite artificial). And assumption 4.3 is fullfiled (except the periodicity property once again) if $k$ belong to $C^{2,2\alpha}$ with $\alpha < 1/2$ or belong to $C^{3,2\alpha-1}$ with $\alpha \ge 1/2$. For a proof in the case where $k$ is supposed to be Lipschitz, see [18]. The extension to the case where $k$ belongs to any Hölder space is straightforward.
## 4.2 Solution of the flow equation and its approximation using finite elements

**Proposition 4.4. Equation (4.1) admits a unique solution p.**

1. If Assumption 4.2 holds for some $0 < \alpha < 1$, then $p \in L^q(\Omega, C_b^{1,\alpha}(\mathbb{R}^d))$ for any finite $q \ge 1$.

2. If Assumption 4.3 holds for some $0 < \alpha < 1$, then $p \in L^q(\Omega, C_b^{2,\alpha}(\mathbb{R}^d))$ for any finite $q \ge 1$.

*Proof.* 1. The result is a consequence of an extension to the case of periodic boundary conditions of Theorem 3.1 of [19]: if $f, g, a$ are O-periodic such that $f$ belongs to $L_{loc}^r(\mathbb{R}^d)$ for some $r > d$ such that $\alpha < 1 - \frac{d}{r}$, $g$ belongs to $C_b^{0,\alpha}(\mathbb{R}^d, \mathbb{R}^d)$ and $a$ belongs to $C_b^{0,\alpha}(\mathbb{R}^d)$ for some $0 < \alpha < 1$, with for any $x$ in $\mathbb{R}^d$ $a(x) \ge a_{min} > 0$, we have then that the O-periodic solution $u$ of:

$$ \left\{ \begin{array}{l} \mathrm{div}(a\nabla u) = f + \mathrm{div} g \quad \text{on } \mathbb{R}^d, \\ \\ \int_O u(x)dx = 0, \end{array} \right. \qquad (4.2) $$

which is classically unique in the space $H_{per}^1$ of locally $H^1$ functions which are O-periodic, also belongs to $C_b^{1,\alpha}(\mathbb{R}^d)$ with

$$ \|u\|_{C_b^{1,\alpha}(\mathbb{R}^d)} \le P_1(a_{min}, \|a\|_{C_b^{0,\alpha}(\mathbb{R}^d)})(\|f\|_{L^r(O)} + \|g\|_{C_b^{0,\alpha}(\mathbb{R}^d)}), $$

where $P_1$ is a polynomial function, whose coefficients do not depend on $a, f$ and $g$. This result can be obtained by adapting the proof of Theorem 3.1 of [19] to the case where the spatial domain is the torus associated to the box $O$. We then apply this inequality (with $g = 0$) for almost all $\omega$ to $p(\omega, \cdot)$, where $p$ is the solution of (4.1), which yields that for almost all $\omega$ we have

$$ \|p(\omega, \cdot)\|_{C_b^{1,\alpha}(\mathbb{R}^d)} \le P_1(a_{min}(\omega), \|a(\omega)\|_{C_b^{0,\alpha}(\mathbb{R}^d)}) \|f\|_{L^r(O)}. $$

The first part of the Proposition follows then from this inequality, Assumption 4.2 and Hölder inequality.

2. From the previous regularity result on the solution of (4.2), we also deduce a similar $C^{2,\alpha}$ regularity result under additional regularity assumptions on $a, f$ and $g$. More precisely, if we suppose moreover that $f$ belongs to $W_{loc}^{1,r}(\mathbb{R}^d)$ and $a$ belongs to $C_b^{1,\alpha}(\mathbb{R}^d)$, then for any $1 \le i \le d$, if $u$ is the solution of (4.2) (with $g = 0$) then $\frac{\partial u}{\partial x_i}$ solves the equation (4.2) with $\frac{\partial f}{\partial x_i}$ instead of $f$ and $g = -\mathrm{div}(\frac{\partial a}{\partial x_i}\nabla u)$. We can therefore deduce that the solution $u$ belongs to $C_b^{2,\alpha}(\mathbb{R}^d)$ with

$$ \begin{aligned} \|u\|_{C_b^{2,\alpha}(\mathbb{R}^d)} &\le P_1(a_{min}, \|a\|_{C_b^{1,\alpha}(\mathbb{R}^d)})(\|f\|_{W^{1,r}(O)} + \|u\|_{C_b^{1,\alpha}(\mathbb{R}^d)}\|a\|_{C_b^{0,\alpha}(\mathbb{R}^d)}) \\ &\le P_2(a_{min}, \|a\|_{C_b^{1,\alpha}(\mathbb{R}^d)})(1 + \|f\|_{W^{1,r}(O)}), \end{aligned} $$

where $P_2$ is a polynomial function, whose coefficients do not depend on $a$ and $f$. It remains to apply this result to $p(\omega, \cdot)$, where $p$ is the solution of (4.1) and where we have fixed $\omega$, which yields that for almost all $\omega$ we have

$$ \|p(\omega, \cdot)\|_{C_b^{2,\alpha}(\mathbb{R}^d)} \le P_2(a_{min}(\omega), \|a(\omega)\|_{C_b^{1,\alpha}(\mathbb{R}^d)})(1 + \|f\|_{W^{1,r}(O)}). $$

We conclude thanks to Hölder inequality and Assumption 4.3. $\square$

Let $V_h$ be a finite element space of O-periodic, continuous, piecewise linear functions, whose integral on the domain $O$ is equal to 0, associated to a shape-regular family of simplicial triangulations of $O$, parametrized by its mesh width $h$. We consider for almost any $\omega$ the finite element approximation $p_h(\omega, \cdot)$ of $p(\omega, \cdot)$ in the finite element space $V_h$.

We define then for almost any $\omega$ the Darcy velocity: $v(\omega, x) = -a(\omega, x)\nabla p(\omega, x)$ and its approximation $v_h(\omega, x) = -a(\omega, x)\nabla p_h(\omega, x)$. Note that the obtained velocity $v_h$ is discontinuous.
**Proposition 4.5.** 1. Let Assumption 4.2 hold for some $0 < \alpha < 1$, then $v \in L^q(\Omega, C_b^{0,\alpha}(\mathbb{R}^d))$ for any finite $q \ge 1$.

2. Let Assumption 4.3 hold for some $0 < \alpha < 1$, then $v \in L^q(\Omega, C_b^{1,\alpha}(\mathbb{R}^d))$ for any finite $q \ge 1$.

3. In both cases (i.e. if Assumption 4.2 holds) we have $v_h \in L^q(\Omega, L^\infty(\mathbb{R}^d))$ for any finite $q \ge 1$.

*Proof.* The first two results follow easily from Hölder inequality and Proposition 4.4 together with Assumptions 4.2 and 4.3. In order to prove the third result, we first notice that for any $h$ and almost every $\omega$ we have

$$ \begin{aligned} \|\nabla p_h(\omega)\|_{L^2(O)} &\le \sqrt{\frac{a_{max}(\omega)}{a_{min}(\omega)}} \|\nabla p\|_{L^2(O)} \\ &\le \sqrt{\frac{a_{max}(\omega)}{a_{min}(\omega)}} \sqrt{|O|} \|p\|_{C_b^1(\mathbb{R}^d)}. \end{aligned} $$

Moreover $V_h$ is a finite dimensional space, and therefore the norm $\|\nabla \cdot \|_{L^\infty(O)}$ is equivalent to the $\|\nabla \cdot \|_{L^2(O)}$ norm on $V_h$ which implies, using the previous bound, that there exists a constant $C_h$ such that for almost all $\omega$ we have

$$ \|\nabla p_h(\omega)\|_{L^\infty(\mathbb{R}^d)} \le C_h \|p(\omega)\|_{C_b^1(\mathbb{R}^d)}. $$

It remains then to apply Proposition 4.4, Assumption 4.2 and Hölder inequality. □

We have then the following error bound for the computation of the Darcy velocity.

**Proposition 4.6.** 1. Let Assumption 4.2 hold for some $0 < \alpha < 1$, then for any $1 \le q < +\infty$, there exists a constant $C_1(q, \alpha)$ such that for any $h > 0$ we have

$$ \|v - v_h\|_{L^q(\Omega, L^\infty(\mathbb{R}^d))} \le C_1(q, \alpha) h^\alpha |\ln(h)|. $$

2. Let Assumption 4.3 hold for some $0 < \alpha < 1$, then for any $1 \le q < +\infty$, there exists a constant $\tilde{C}_1(q)$ such that for any $h > 0$ we have

$$ \|v - v_h\|_{L^q(\Omega, L^\infty(\mathbb{R}^d))} \le \tilde{C}_1(q) h |\ln(h)|. $$

*Proof.* We notice that Theorem 4 of [29] can be adapted to the case of periodic boundary conditions. It implies then that if assumption 4.2 holds, then for almost all $\omega$ we have

$$ \|(p - p_h)(\omega)\|_{W^{1,\infty}(\mathbb{R}^d)} \le \frac{\|a(\omega)\|_{C_b^0(\mathbb{R}^d)}}{a_{min}(\omega)} h^\alpha |\ln h| \|p(\omega)\|_{C_b^{1,\alpha}(\mathbb{R}^d)}. $$

And if Assumption 4.3 holds, then for almost all $\omega$ we have

$$ \|(p - p_h)(\omega)\|_{W^{1,\infty}(\mathbb{R}^d)} \le \frac{\|a(\omega)\|_{C_b^0(\mathbb{R}^d)}}{a_{min}(\omega)} h |\ln h| \|p(\omega)\|_{C_b^2(\mathbb{R}^d)}. $$

The results follows then from these bounds, Hölder inequality and assumptions 4.2 and 4.3 respectively. □

In the next two subsections, the variable $\omega$ modelizing uncertainty will be fixed. Therefore, for the sake of readability the variable $\omega$ will be omitted. However, we will give explicit bounds, which will enable us to track the dependence on $\omega$ and to integrate with respect to $\omega$ in the final step.
## 4.3 The advection-diffusion equation

We consider an initial condition $c_0 \in L^1(\mathbb{R}^d)$, with $\int_{\mathbb{R}^d} c_0(x)dx = 1$, $c_0(x) \ge 0$ for any $x \in \mathbb{R}^d$ and take $v \in C_b^{0,\alpha}(\mathbb{R}^d, \mathbb{R}^d)$ for some $0 < \alpha < 1$. We consider then the following advection-diffusion equation:

$$ \left\{ \begin{array}{l} \frac{\partial c}{\partial t}(x,t) + v(x) \cdot \nabla c(x,t) - D \Delta c(x,t) = 0, \quad x \in \mathbb{R}^d \text{ and } t \in [0, T] \\ c(x,0) = c_0(x), \quad x \in \mathbb{R}^d. \end{array} \right. \qquad (4.3) $$

We consider also the stochastic differential equation associated to this Fokker-Planck equation:

$$ \left\{ \begin{array}{l} dX(t) = v(X(t))dt + \sqrt{2D}dW(t), \\ X(0) = X_0, \end{array} \right. \qquad (4.4) $$

We suppose that $X_0$ admits $c_0(x)dx$ as density. The SDE (4.4) admits a unique strong solution (it is very classical if $v$ is lipschitz continuous, and it follows for example from [21] otherwise). We recall the following well known result (see [20] e.g. for the link with the SDE and [22] e.g. for the regularity result).

**Proposition 4.7.** The equation (4.3) admits a unique solution $c \in C^0([0, T], C^2(\mathbb{R}^d)) \cap C^0([0, T], L^2(\mathbb{R}^d))$ and $X(t)$ admits $c(x, t)dx$ as density.

## 4.4 Time discretization

We consider $(\Omega', \mathcal{F}', \mathcal{P}')$ an another probability space, whose generic variable is denoted by $\xi$. Here we give bounds of the weak error resulting both from the time discretization of a stochastic differential equation (with an additive noise and a $C^{0,\alpha}$ or $C^{1,\alpha}$ drift) and from the spatial approximation of the drift. Let $v \in C_b^{0,\alpha}(\mathbb{R}^d, \mathbb{R}^d)$ for some $0 < \alpha < 1$. We denote by $X^x$ the solution of the following stochastic differential equation:

$$ \left\{ \begin{array}{l} dX^x(t) = v(X^x(t))dt + \sqrt{2D}dW(t), \\ X^x(0) = x. \end{array} \right. \qquad (4.5) $$

For the same reasons as seen in the case of equation (4.4), the SDE (4.5) admits a unique strong solution. We denote by $X_n^x$ the numerical approximation of $X^x$ using an Euler scheme (as in section 3), where the mesh of the time discretization is $\Delta t = \frac{T}{n}$, and $t_k = k\Delta t$ for $0 \le k \le n$. We extend $X_n^x$ to a function defined for all $t \ge 0$ by:

$$ \left\{ \begin{array}{l} dX_n^x(t) = v(X_n^x(t_k))dt + \sqrt{2D}dW(t), \quad \text{for } t_k \le t \le t_{k+1}, \\ X_n^x(0) = x. \end{array} \right. \qquad (4.6) $$

We also define the Euler scheme with an approximated velocity $\tilde{v}$, where $\tilde{v} \in L^\infty(\mathbb{R}^d)$:

$$ \left\{ \begin{array}{l} d\tilde{X}_n^x(t) = \tilde{v}(\tilde{X}_n^x(t_k))dt + \sqrt{2D}dW(t), \quad \text{for } t_k \le t \le t_{k+1}, \\ \tilde{X}_n(0) = x. \end{array} \right. \qquad (4.7) $$

The next subsection will be devoted to the estimation of the error committed when we approximate the law of $X$ by the law of $\tilde{X}_n$.

## 4.5 Weak error for both time and space discretizations

We denote by $C_b^{1;2}([0, T] \times \mathbb{R}^d)$ the space of functions of $(t, x)$ which admit one derivative with respect to $t$ and two derivatives with respect to $x$, all these derivatives being continuous and bounded on $[0, T] \times \mathbb{R}^d$. For $u \in C_b^{1;2}([0, T] \times \mathbb{R}^d)$, we introduce the natural norm

$$ \|u\|_{C_b^{1;2}([0,T] \times \mathbb{R}^d)} = \max \left\{ \|u\|_{C_b^0([0,T] \times \mathbb{R}^d)}, \left\| \frac{\partial u}{\partial t} \right\|_{C_b^0([0,T] \times \mathbb{R}^d)}, \left\| \frac{\partial u}{\partial x} \right\|_{C_b^0([0,T] \times \mathbb{R}^d)}, \left\| \frac{\partial^2 u}{\partial x^2} \right\|_{C_b^0([0,T] \times \mathbb{R}^d)} \right\}. $$
We introduce the following Kolmogorov equation associated to the previous SDE (4.4):

$$ egin{cases} \frac{\partial u}{\partial t}(t, x) = D\Delta u(t, x) + v(x) \cdot \nabla u(t, x) \\ u(0, x) = \varphi(x). \end{cases} \tag{4.8} $$

**Proposition 4.8.** Let $0 < \alpha < 1$, $\varphi \in C_b^{1,\alpha}(\mathbb{R}^d)$ and $v \in C_b^{0,\alpha}(\mathbb{R}^d)$, then the Kolmogorov equation (4.8) admits a unique solution $u$ and for any $T > 0$ there exists a constant $C_2(\alpha, T)$ such that we get

$$ \sup_{t \in [0, T]} \|u(t)\|_{C_b^{1,\alpha}(\mathbb{R}^d)} \le C_2(\alpha, T) (\|\varphi\|_{C_b^{1,\alpha}(\mathbb{R}^d)} + \|\varphi\|_{C_b^0(\mathbb{R}^d)} \|v\|_{C_b^0(\mathbb{R}^d)}^{1+\alpha}). \tag{4.9} $$

*Proof.* To begin with, we recall a classical result whose proof can be found in [22] page 184 : for $0 < \alpha < 1$, if $v \in C_b^{0,\alpha}(\mathbb{R}^d)$ and $\varphi \in C_b^{2,\alpha}(\mathbb{R}^d)$ the equation (4.8) admits a unique solution $u \in C_b^{1;2+\alpha}([0, T] \times \mathbb{R}^d)$. However this result does not enables us to conclude : we have a weaker regularity assumption on $\varphi$ and, more importantly, we need an explicit bound for $u$ (for a weaker norm). In a first step we suppose that $\varphi \in C_b^{2,\alpha}(\mathbb{R}^d)$ and will weaken this assumption later by regularization. We can then consider the unique solution $u$ of the equation (4.8). It follows from [22] that $u \in C_b^{1;2+\alpha}([0, T] \times \mathbb{R}^d)$. We denote by $S(t)$ the heat semi-group on $\mathbb{R}^d$, and using it we get that

$$ u(t) = S(t)\varphi + \int_0^t S(t-s)(v \cdot \nabla u)(s)ds. $$

It is classical that there exists a constant $C$ such that for any $t \in [0, T]$, $S(t)$ is both a continuous operator from $C_b^{1,\alpha}(\mathbb{R}^d)$ to itself with norm $C$, and a continuous operator from $C_b^0(\mathbb{R}^d)$ to $C_b^{1,\alpha}(\mathbb{R}^d)$ with norm $Ct^{-\frac{1+\alpha}{2}}$. We deduce that for any $t \in [0, T]$

$$ \begin{aligned} \|u(t)\|_{C_b^{1,\alpha}(\mathbb{R}^d)} &\le C \|\varphi\|_{C_b^{1,\alpha}(\mathbb{R}^d)} + C \|v\|_{C_b^0(\mathbb{R}^d)} \sup_{t \in [0,T]} \|\nabla u(t)\|_{C_b^0(\mathbb{R}^d)} \int_0^t (t-s)^{-\frac{1+\alpha}{2}} ds \\ &\le C \|\varphi\|_{C_b^{1,\alpha}(\mathbb{R}^d)} + C \|v\|_{C_b^0(\mathbb{R}^d)} \sup_{t \in [0,T]} \|\nabla u(t)\|_{C_b^0(\mathbb{R}^d)} T^{\frac{1-\alpha}{2}}. \end{aligned} $$

We now use a classical interpolation inequality : for any $h \in C_b^{1,\alpha}(\mathbb{R}^d)$ we have

$$ \|\nabla h\|_{C_b^0(\mathbb{R}^d)} \le 3 \|h\|_{C_b^{0,\frac{\alpha}{1+\alpha}}(\mathbb{R}^d)} \|\nabla h\|_{C_b^{0,\alpha}(\mathbb{R}^d)}^{\frac{1}{1+\alpha}} $$

We deduce then from this interpolation inequality combined with Young inequality that there exists a constant $c$ such that for any $h \in C_b^{1,\alpha}(\mathbb{R}^d)$ and $\varepsilon > 0$ we have

$$ \|\nabla h\|_{C_b^0(\mathbb{R}^d)} \le \varepsilon \|\nabla h\|_{C_b^{0,\alpha}(\mathbb{R}^d)} + \frac{c}{\varepsilon^{\alpha}} \|h\|_{C_b^0(\mathbb{R}^d)}. \tag{4.10} $$

Using this inequality we get :

$$ \sup_{t \in [0, T]} \|u(t)\|_{C_b^{1,\alpha}(\mathbb{R}^d)} \le C \|\varphi\|_{C_b^{1,\alpha}(\mathbb{R}^d)} + C T^{\frac{1-\alpha}{2}} \|v\|_{C_b^0(\mathbb{R}^d)} \left( \varepsilon \sup_{t \in [0, T]} \|\nabla u(t)\|_{C_b^{0,\alpha}(\mathbb{R}^d)} + \frac{c}{\varepsilon^{\alpha}} \sup_{t \in [0, T]} \|u(t)\|_{C_b^0(\mathbb{R}^d)} \right). $$

It remains to take $\varepsilon = \frac{1}{2CT^{\frac{1-\alpha}{2}} \|v\|_{C_b^0(\mathbb{R}^d)}}$ and we finally get

$$ \sup_{t \in [0, T]} \|u(t)\|_{C_b^{1,\alpha}(\mathbb{R}^d)} \le 2C \|\varphi\|_{C_b^{1,\alpha}(\mathbb{R}^d)} + c(2CT^{\frac{1-\alpha}{2}} \|v\|_{C_b^0(\mathbb{R}^d)})^{1+\alpha} \sup_{t \in [0, T]} \|u(t)\|_{C_b^0(\mathbb{R}^d)} \tag{4.11} $$

Moreover we have classically $u(t, x) = \mathbb{E}[\varphi(X^x(t))]$ for any $t \in [0, T]$ and any $x \in \mathbb{R}^d$ (it follows from Itô formula applied for $t \in [\varepsilon, T]$ for any $\varepsilon > 0$) therefore

$$ \sup_{t \in [0, T]} \|u(t)\|_{C_b^0(\mathbb{R}^d)} \le \|\varphi\|_{C_b^0(\mathbb{R}^d)}. $$
Note that more generally, we deduce from Itô formula that for any $t \le s \in [0, T]$ and $x \in \mathbb{R}^d$,

$$ \mathbb{E}[u(T-t, x)] = \mathbb{E}[u(T-s, X^x(s-t)] \quad (4.12) $$

which will be used below. We deduce from (4.11) together with this inequality that there exists a constant $C_2(\alpha, T)$ such that

$$ \sup_{t \in [0, T]} \|u(t)\|_{C_b^{1, \alpha}(\mathbb{R}^d)} \le C_2(\alpha, T)(\|\varphi\|_{C_b^{1, \alpha}(\mathbb{R}^d)} + \|\varphi\|_{C_b^0(\mathbb{R}^d)} \|v\|_{C_b^0(\mathbb{R}^d)}^{1+\alpha}). $$

Using a regularization argument on $\varphi$ and the previous inequality, it can easily be seen that if we suppose $\varphi$ to belong only to $C^{1, \alpha}(\mathbb{R}^d)$, we have then $u(t) \in C_b^{1, \alpha}(\mathbb{R}^d)$ for any $t \in [0, T]$ and the inequality (4.9) still holds. $\square$

**Proposition 4.9.** 1. Let $0 < \alpha < 1$, then for any $x \in \mathbb{R}^d$, $\varphi \in C_b^{1, \alpha}(\mathbb{R}^d)$, $v \in C_b^{0, \alpha}(\mathbb{R}^d)$, and $\tilde{v} \in L^\infty(\mathbb{R}^d)$, then recalling that $X^x$ and $\tilde{X}_n^x$ are respectively the solutions of (4.5) and (4.7), we have

$$ |\mathbb{E}[\varphi(X^x(T)) - \varphi(\tilde{X}_n^x(T))]| \le P_3(\|v\|_{C_b^{0, \alpha}(\mathbb{R}^d)}, \|\varphi\|_{C_b^{1, \alpha}(\mathbb{R}^d)}((\Delta t)^{\frac{\alpha}{2}} + \|v - \tilde{v}\|_{L^\infty(\mathbb{R}^d)})), $$

where $P_3$ is a polynomial function of $\|\varphi\|_{C^{1, \alpha}(\mathbb{R}^d)}$ and $\|v\|_{C_b^{0, \alpha}(\mathbb{R}^d)}$, whose coefficients only depend on $\alpha$ and $T$.

2. Let $0 < \alpha < 1$, then for any $x \in \mathbb{R}^d$, $\varphi \in C_b^{1, \alpha}(\mathbb{R}^d)$, $v \in C_b^{1, \alpha}(\mathbb{R}^d)$ and $\tilde{v} \in L^\infty(\mathbb{R}^d)$, recalling that $X^x$ and $\tilde{X}_n^x$ are respectively the solutions of (4.5) and (4.7), we have

$$ |\mathbb{E}[\varphi(X^x(T)) - \varphi(\tilde{X}_n^x(T))]| \le \tilde{P}_3(\|v\|_{C_b^{1, \alpha}(\mathbb{R}^d)}, \|\varphi\|_{C_b^{1, \alpha}(\mathbb{R}^d)}((\Delta t)^{\frac{1+\alpha}{2}} + \|v - \tilde{v}\|_{L^\infty(\mathbb{R}^d)} + \|v - \tilde{v}\|_{L^\infty(\mathbb{R}^d)}^{1+\alpha}), $$

where $\tilde{P}_3$ is a polynomial function of $\|\varphi\|_{C^{1, \alpha}(\mathbb{R}^d)}$ and $\|v\|_{C_b^{1, \alpha}(\mathbb{R}^d)}$, whose coefficients only depend on $\alpha$ and $T$.

**Remark 4.10.** We note that the particular case $v = \tilde{v}$ corresponds to the classical case of the weak error of the Euler scheme on a SDE with a drift having less than $C^2$ regularity This has been studied in [23], nevertheless, we give a simple proof in the case of an additive noise, which moreover yields for a $C^{1, \alpha}$ drift a better weak convergence order in our case (namely $\frac{1+\alpha}{2}$) than the general result of [23] (namely $\frac{1}{2-\alpha}$). More importantly, we need to get an explicit dependence of the constant on $v$ and $\varphi$, which cannot be found in classical literature such as [23]. Moreover we also include the spatial discretization in the weak error, which is not classical up to our knowledge.

*Proof.* 1. Let $0 < \alpha < 1$, $x \in \mathbb{R}^d$, $\varphi \in C_b^{1, \alpha}(\mathbb{R}^d)$, $v \in C_b^{0, \alpha}(\mathbb{R}^d)$ and $\tilde{v} \in L^\infty(\mathbb{R}^d)$, then using (4.12), the total weak error can be classically (see [15] for example) expressed as

$$ \begin{aligned} E &= \mathbb{E}[\varphi(X^x(T))] - \mathbb{E}[\varphi(\tilde{X}_n^x(T))] \\ &= u(T, x) - \mathbb{E}[u(0, \tilde{X}_n^x(T))], \end{aligned} $$

and can hence be split into $E = \sum_{i=0}^{n-1} E_i$, where:

$$ \begin{aligned} E_i &= \mathbb{E}[u(T - t_i, \tilde{X}_n^x(t_i))] - \mathbb{E}[u(T - t_{i+1}, \tilde{X}_n^x(t_{i+1}))] \\ &= \mathbb{E}[u(T - t_{i+1}, X^{\tilde{X}_n^x}(t_i)(\Delta t))] - \mathbb{E}[u(T - t_{i+1}, \tilde{X}_n^{\tilde{X}_n^x}(t_i)(\Delta t))] \\ &= \mathbb{E}[\mathbb{E}[u(T - t_{i+1}, X^{\tilde{X}_n^x}(t_i)(\Delta t)) - u(T - t_{i+1}, \tilde{X}_n^{\tilde{X}_n^x}(t_i)(\Delta t))|\tilde{X}_n^x(t_i)]] \\ &= \mathbb{E}[e_i(\tilde{X}_n^x(t_i))], \end{aligned} $$

where $e_i(y) = \mathbb{E}[u(T - t_{i+1}, X^y(\Delta t)) - u(T - t_{i+1}, \tilde{X}_n^y(\Delta t))]$, by using the Markov property of the solution $X$ and of the discretized solution $\tilde{X}_n$ as well as (4.12).
In order to bound the term $e_i$, we first notice that we have for any $s \in ]0, \Delta t]$:

$$ \begin{aligned} \|X^y(s) - y\| &= \left\| \int_0^s v(X^y(t))dt + \sqrt{2DW}(s) \right\| \\ &\le \|v\|_{C_b^0(\mathbb{R}^d)} s + \sqrt{2D}|W(s)|, \end{aligned} $$

which implies that

$$ \begin{aligned} \mathbb{E}[\|X^y(s) - y\|^\alpha] &\le \mathbb{E} \left[ \left( \|v\|_{C_b^0(\mathbb{R}^d)} \Delta t + \sqrt{2D\Delta t} \frac{|W(s)|}{\sqrt{s}} \right)^\alpha \right] \\ &\le (\Delta t)^{\frac{\alpha}{2}} (\|v\|_{C_b^0(\mathbb{R}^d)} T^{\frac{\alpha}{2}} + (2D)^{\frac{\alpha}{2}} \mathbb{E}[|Y|^\alpha]), \\ &\le (\Delta t)^{\frac{\alpha}{2}} (\|v\|_{C_b^0(\mathbb{R}^d)} T^{\frac{\alpha}{2}} + 2(2D)^{\frac{\alpha}{2}}), \end{aligned} \quad (4.13) $$

where $Y$ is a standard normal deviate. Using this bound together with the equality

$$ X^y(\Delta t) - \tilde{X}_n^y(\Delta t) = \int_0^{\Delta t} v(X^y(s)) - v(y)ds + \Delta t(v(y) - \tilde{v}(y)), $$

we deduce the following bound for $e_i$:

$$ \begin{aligned} |e_i(y)| &\le \|u(T - t_{i+1})\|_{C_b^1(\mathbb{R}^d)} \mathbb{E} \left[ \|X^y(\Delta t) - \tilde{X}_n^y(\Delta t)\| \right] \\ &\le \sup_{t \in [0,T]} \|u(t)\|_{C_b^1(\mathbb{R}^d)} \left( \|v\|_{C_b^{0,\alpha}(\mathbb{R}^d)} \int_0^{\Delta t} \mathbb{E}[\|X^y(s) - y\|^\alpha]ds + \Delta t \|v - \tilde{v}\|_{L^\infty(\mathbb{R}^d)} \right) \\ &\le \Delta t C_2(\alpha, T) (\|\varphi\|_{C_b^{1,\alpha}(\mathbb{R}^d)} + \|\varphi\|_{C_b^0(\mathbb{R}^d)} \|v\|_{C_b^0(\mathbb{R}^d)}^{1+\alpha}) \\ &\quad \times (\|v\|_{C_b^{0,\alpha}(\mathbb{R}^d)} (\Delta t)^{\frac{\alpha}{2}} (\|v\|_{C_b^0(\mathbb{R}^d)} T^{\frac{\alpha}{2}} + 2(2D)^{\frac{\alpha}{2}}) + \|v - \tilde{v}\|_{L^\infty(\mathbb{R}^d)}) \\ &\le \Delta t P_3 (\|v\|_{C_b^{0,\alpha}(\mathbb{R}^d)}, \|\varphi\|_{C_b^{1,\alpha}(\mathbb{R}^d)})(\Delta t)^{\frac{\alpha}{2}} + \|v - \tilde{v}\|_{L^\infty(\mathbb{R}^d)}). \end{aligned} \quad (4.14) $$

where we have used Proposition 4.8 and where $P_3$ is a polynomial function of $\|v\|_{C_b^{0,\alpha}(\mathbb{R}^d)}$ and $\|\varphi\|_{C_b^{1,\alpha}(\mathbb{R}^d)}$ whose coefficients only depend on $\alpha$ and $T$. It remains to take the sum over $i$ to get the bound for the total error $E$.

2. Let now make the additional assumption that $v \in C_b^{1,\alpha}(\mathbb{R}^d)$.

Using a Taylor expansion with integral remainder of $u$ at order one with respect to $x$, we get:

$$ \begin{aligned} e_i(y) &= \mathbb{E} \left[ \int_0^1 D_x u(T - t_{i+1}, X^y(\Delta t) + \theta(\tilde{X}_n^y(\Delta t) - X^y(\Delta t))). (X^y(\Delta t) - \tilde{X}_n^y(\Delta t)) d\theta \right] \\ &= \mathbb{E} \left[ D_x u(T - t_{i+1}, y) (X^y - \tilde{X}_n^y) \right] \\ &\quad + \mathbb{E} \left[ \int_0^1 (D_x u(T - t_{i+1}, X^y + \theta(\tilde{X}_n^y - X^y)) - D_x u(T - t_{i+1}, y)). (X^y - \tilde{X}_n^y) d\theta \right] \\ &= D_x u(T - t_{i+1}, y) \mathbb{E} \left[ X^y - \tilde{X}_n^y \right] \\ &\quad + \mathbb{E} \left[ \int_0^1 (D_x u(T - t_{i+1}, X^y + \theta(\tilde{X}_n^y - X^y)) - D_x u(T - t_{i+1}, y)). (X^y - \tilde{X}_n^y) d\theta \right], \end{aligned} $$

where in the last equality we have denoted $X^y(\Delta t)$ by $X^y$ and $\tilde{X}_n^y(\Delta t)$ by $\tilde{X}_n^y$ for the sake of readability. These shorter notations will be also used in the remainder of the proof, when there is no ambiguity.
We have then

$$
\begin{aligned}
|e_i(y)| &\le \sup_{t \in [0,T]} \|u(t)\|_{C_b^1(\mathbb{R}^d)} \|\mathbb{E}[X^y(\Delta t) - \tilde{X}_n^y(\Delta t)]\| \\
&\quad + \sup_{t \in [0,T]} \|D_x u(t)\|_{C_b^{0,\alpha}(\mathbb{R}^d)} \mathbb{E} \left[ \|X^y(\Delta t) - y\|^\alpha \|\tilde{X}_n^y(\Delta t) - X^y(\Delta t)\| \right] \\
&\quad + \sup_{t \in [0,T]} \|D_x u(t)\|_{C_b^{0,\alpha}(\mathbb{R}^d)} \mathbb{E} \left[ \|\tilde{X}_n^y(\Delta t) - X^y(\Delta t)\|^{1+\alpha} \right]
\end{aligned}
$$

In order to bound these terms, we need to bound $X^y(\Delta t) - \tilde{X}_n^y(\Delta t)$:

$$
\begin{aligned}
X^y(\Delta t) - \tilde{X}_n^y(\Delta t) &= \int_0^{\Delta t} v(X^y(s)) - v(y)ds + \Delta t(v(y) - \tilde{v}(y)) \\
&= \int_0^{\Delta t} Dv(y)(X^y(s) - y)ds + \Delta t(v(y) - \tilde{v}(y)) \\
&\quad + \int_0^{\Delta t} \int_0^1 (Dv(y + \theta(X^y(s) - y)) - Dv(y))(X^y(s) - y)dsd\theta.
\end{aligned}
\tag{4.15}
$$

We deduce first that

$$
\begin{aligned}
\|\mathbb{E}[X^y(\Delta t) - \tilde{X}_n^y(\Delta t)]\| &\le \|v\|_{C_b^1(\mathbb{R}^d)} \int_0^{\Delta t} \|E[X^y(s) - y]\| ds + \Delta t \|v - \tilde{v}\|_{L^\infty(\mathbb{R}^d)} \\
&\quad + \|Dv\|_{C_b^{0,\alpha}(\mathbb{R}^d)} \int_0^{\Delta t} \mathbb{E}[\|X^y(s) - y\|^{1+\alpha}] ds \\
&\le \|v\|_{C_b^1(\mathbb{R}^d)} \Delta t^2 \|v\|_{C_b^0(\mathbb{R}^d)} + \Delta t \|v - \tilde{v}\|_{L^\infty(\mathbb{R}^d)} \\
&\quad + \|Dv\|_{C_b^{0,\alpha}(\mathbb{R}^d)} \int_0^{\Delta t} \mathbb{E}[\|X^y(s) - y\|^{1+\alpha}] ds,
\end{aligned}
\tag{4.16}
$$

where we have used that

$$
X^y(s) - y = \int_0^s v(X^y(t))dt + \sqrt{2DW}(s),
$$

which implies that

$$
\begin{aligned}
\|\mathbb{E}[X^y(s) - y]\| &= \left\| \int_0^s \mathbb{E}[v(X^y(t))]dt \right\| \\
&\le \Delta t \|v\|_{C_b^0(\mathbb{R}^d)}.
\end{aligned}
$$

In order to bound the second term of (4.16), we first recall that

$$
\begin{aligned}
\|X^y(s) - y\| &= \left| \int_0^s v(X^y(t))dt + \sqrt{2DW}(t) \right| \\
&\le \|v\|_{C_b^0(\mathbb{R}^d)}s + \sqrt{2D}|W(s)|,
\end{aligned}
$$

which implies that for any $s \in ]0, \Delta t]$ we have

$$
\begin{aligned}
\mathbb{E}[\|X^y(s) - y\|^{1+\alpha}] &\le \mathbb{E} \left[ \left( \|v\|_{C_b^0(\mathbb{R}^d)} \Delta t + \sqrt{2D\Delta t} \frac{|W(s)|}{\sqrt{s}} \right)^{1+\alpha} \right] \\
&\le 2(\Delta t)^{\frac{1+\alpha}{2}} (\|v\|_{C_b^0(\mathbb{R}^d)} T^{\frac{1+\alpha}{2}} + 2(2D)^{\frac{1+\alpha}{2}}).
\end{aligned}
\tag{4.17}
$$
This inequality finally enables us to bound $\|\mathbb{E}[X^y(\Delta t) - \tilde{X}_n^y(\Delta t)]\|$:

$$ \begin{aligned} \|\mathbb{E}[X^y(\Delta t) - \tilde{X}_n^y(\Delta t)]\| & \le \|v\|_{C_b^1(\mathbb{R}^d)} \Delta t^2 \|v\|_{C_b^0(\mathbb{R}^d)} + \Delta t \|v - \tilde{v}\|_{L^\infty(\mathbb{R}^d)} \\ & \quad + 2 \|Dv\|_{C_b^{0,\alpha}(\mathbb{R}^d)} (\Delta t)^{1+\frac{1+\alpha}{2}} (\|v\|_{C_b^0(\mathbb{R}^d)}^{1+\alpha} T^{\frac{1+\alpha}{2}} + 2(2D)^{\frac{1+\alpha}{2}}) \\ & \le \Delta t \|v - \tilde{v}\|_{L^\infty(\mathbb{R}^d)} + \|v\|_{C_b^{1,\alpha}(\mathbb{R}^d)} (\Delta t)^{1+\frac{1+\alpha}{2}} \\ & \quad \times (\|v\|_{C_b^{1,\alpha}(\mathbb{R}^d)} T^{\frac{1-\alpha}{2}} + 2T^{\frac{1+\alpha}{2}} \|v\|_{C_b^{1,\alpha}(\mathbb{R}^d)}^{1+\alpha} + 4(2D)^{\frac{1+\alpha}{2}}). \end{aligned} \quad (4.18) $$

We now use (4.15) together with (4.17) and Hölder inequality to get a bound for $\mathbb{E}[\|X^y(\Delta t) - \tilde{X}_n^y(\Delta t)\|^{1+\alpha}]$:

$$ \begin{aligned} \mathbb{E}[\|X^y(\Delta t) - \tilde{X}_n^y(\Delta t)\|^{1+\alpha}] & \le 2(\Delta t)^\alpha \mathbb{E} \left[ \int_0^{\Delta t} \|v(X^y(s)) - v(y)\|^{1+\alpha} ds \right] + 2(\Delta t)^{1+\alpha} \|v - \tilde{v}\|_{L^\infty(\mathbb{R}^d)}^{1+\alpha} \\ & \le 2 \|Dv\|_{C_b^1(\mathbb{R}^d)}^{1+\alpha} (\Delta t)^\alpha \mathbb{E} \left[ \int_0^{\Delta t} \|X^y(s) - y\|^{1+\alpha} ds \right] \\ & \quad + 2(\Delta t)^{1+\alpha} \|v - \tilde{v}\|_{L^\infty(\mathbb{R}^d)}^{1+\alpha} \\ & \le 4 \|v\|_{C_b^1(\mathbb{R}^d)}^{1+\alpha} (\Delta t)^{3(\frac{1+\alpha}{2})} (\|v\|_{C_b^0(\mathbb{R}^d)}^{1+\alpha} T^{\frac{1+\alpha}{2}} + 2(2D)^{\frac{1+\alpha}{2}}) \\ & \quad + 2(\Delta t)^{1+\alpha} \|v - \tilde{v}\|_{L^\infty(\mathbb{R}^d)}^{1+\alpha}. \end{aligned} \quad (4.19) $$

It remains to bound $\mathbb{E}[\|X^y(\Delta t) - y\|^\alpha \|\tilde{X}_n^y(\Delta t) - X^y(\Delta t)\|]$ to get a bound for $e_i$. To get such a result, we use bounds similar to the one used to get inequalities (4.19) and (4.13),

$$ \begin{aligned} \mathbb{E} [\|X^y(\Delta t) - y\|^\alpha \|\tilde{X}_n^y(\Delta t) - X^y(\Delta t)\|] & \\ & \le \mathbb{E} \left[ \int_0^{\Delta t} \|X^y(\Delta t) - y\|^\alpha \|\tilde{v}(y) - v(X^y(s))\| ds \right] \\ & \le \mathbb{E} \left[ \|X^y(\Delta t) - y\|^\alpha \left( (\|Dv\|_{C_b^0(\mathbb{R}^d)} \int_0^{\Delta t} \|X^y(s) - y\| ds + \Delta t \|v - \tilde{v}\|_{L^\infty(\mathbb{R}^d)}) \right) \right] \\ & \le 2 \|v\|_{C_b^1(\mathbb{R}^d)} (\Delta t)^{1+\frac{1+\alpha}{2}} (\|v\|_{C_b^0(\mathbb{R}^d)}^{1+\alpha} T^{\frac{1+\alpha}{2}} + 2(2D)^{\frac{1+\alpha}{2}}) \\ & \quad + \Delta t T^{\frac{\alpha}{2}} (\|v\|_{C_b^0(\mathbb{R}^d)}^{\alpha} T^{\frac{\alpha}{2}} + 2(2D)^{\frac{\alpha}{2}}) \|v - \tilde{v}\|_{L^\infty(\mathbb{R}^d)} \end{aligned} \quad (4.20) $$

The estimates (4.18),(4.19),(4.20) lead to the following bound for $e_i$:

$$ \begin{aligned} |e_i(y)| & \le \Delta t \sup_{t \in [0,T]} \|u(t)\|_{C_b^1(\mathbb{R}^d)} [\|v - \tilde{v}\|_{L^\infty(\mathbb{R}^d)} \\ & \quad + \|v\|_{C_b^{1,\alpha}(\mathbb{R}^d)} (\Delta t)^{\frac{1+\alpha}{2}} (\|v\|_{C_b^{1,\alpha}(\mathbb{R}^d)} T^{\frac{1-\alpha}{2}} + 2T^{\frac{1+\alpha}{2}} \|v\|_{C_b^{1,\alpha}(\mathbb{R}^d)}^{1+\alpha} + 4(2D)^{\frac{1+\alpha}{2}})] \\ & \quad + \Delta t \sup_{t \in [0,T]} \|D_x u(t)\|_{C_b^{0,\alpha}(\mathbb{R}^d)} [2\|v\|_{C_b^1(\mathbb{R}^d)} (\Delta t)^{\frac{1+\alpha}{2}} (\|v\|_{C_b^0(\mathbb{R}^d)}^{1+\alpha} T^{\frac{1+\alpha}{2}} + (2D)^{\frac{1+\alpha}{2}}) \\ & \quad + T^{\frac{\alpha}{2}} (\|v\|_{C_b^0(\mathbb{R}^d)}^{\alpha} T^{\frac{\alpha}{2}} + 2(2D)^{\frac{\alpha}{2}}) \|v - \tilde{v}\|_{L^\infty(\mathbb{R}^d)}] \\ & \quad + \Delta t \sup_{t \in [0,T]} \|D_x u(t)\|_{C_b^{0,\alpha}(\mathbb{R}^d)} [4\|v\|_{C_b^1(\mathbb{R}^d)}^{1+\alpha} (\Delta t)^{\frac{1+\alpha}{2}} T^{\alpha} (\|v\|_{C_b^0(\mathbb{R}^d)}^{1+\alpha} T^{\frac{1+\alpha}{2}} + 2(2D)^{\frac{1+\alpha}{2}}) \\ & \quad + 2T^{\alpha} \|v - \tilde{v}\|_{L^\infty(\mathbb{R}^d)}^{1+\alpha}] \end{aligned} $$
The final result on the weak error follows by taking the sum over $i$ of these inequalities, recalling that $n\Delta t = T$ and using Proposition 4.8:

$$|E| \le \tilde{P}_3(\|v\|_{C_b^{1,\alpha}(\mathbb{R}^d)}, \|\varphi\|_{C_b^{1,\alpha}(\mathbb{R}^d)}[(\Delta t)^{\frac{1+\alpha}{2}} + \|v-\tilde{v}\|_{L^\infty(\mathbb{R}^d)} + \|v-\tilde{v}\|_{L^\infty(\mathbb{R}^d)}^{1+\alpha}]),$$

where $\tilde{P}_3$ is a polynomial function of $\|\varphi\|_{C_b^{1,\alpha}(\mathbb{R}^d)}$ and $\|v\|_{C_b^{1,\alpha}(\mathbb{R}^d)}$ whose coefficients only depend on $\alpha$ and $T$. $\square$

<figure><img src="image_2.png" /></figure>

## 4.6 Total error on the generalized spread

We recall that $(\Omega, \mathcal{F}, \mathbb{P})$ and $(\Omega', \mathcal{F}', \mathbb{P}')$ are two probability spaces, with generic variables $\omega \in \Omega$ and $\xi \in \Omega'$. We define the stochastic process $X(\omega, \xi, t)$ as the solution for almost all $\omega \in \Omega$ of the following stochastic differential equation:

$$\begin{cases} dX(\omega, \xi, t) = v(\omega, X(\omega, \xi, t))dt + \sqrt{2D}dW(\xi, t), \quad x \in \mathbb{R}^d, t \ge 0, \\ X(\omega, \xi, 0) = X_0(\xi), \end{cases} \qquad (4.21)$$

where $v$ is defined as in subsection 4.2, $W$ is a $d$-dimensionnal brownian motion on $(\Omega', \mathcal{F}', \mathbb{P}')$ and $X_0$ admits $c_0$ as density, as defined in section 4.3. Then we define for any $1 \le i \le N$, $1 \le j \le M$ and almost all $\omega$ the approximations $X_{n,h}^{i,j}(\omega, \xi, t)$ by:

$$\begin{cases} dX_{n,h}^{i,j}(\omega, \xi, t) = v_h^i(\omega, X_{n,h}^{i,j}(\omega, \xi, t_k))dt + \sqrt{2D}dW^j(\xi, t), \quad \text{for } t \in [t_k, t_{k+1}] \\ X_{n,h}^{i,j}(\omega, \xi, 0) = X_0^{i,j}(\xi), \end{cases} \qquad (4.22)$$

where $v_h^i$ is the finite element approximation of $v^i$ as defined in subsection 4.2, the $W^j$ are independent $d$-dimensionnal brownian motion with unit covariance and the $X_0^{i,j}$ are independent random variables of density $c_0$. We define the following quantity of interest, which is a generalization of the spread defined is section 1 by: $\mathbb{E}_\omega[\psi(\mathbb{E}_\xi[\varphi(X(\omega, \xi, T))])]$, for some vector-valued functions $\varphi$ and $\psi$.

**Definition 4.11.** We define the total error on the generalized spread by:

$$Er(\omega, \xi) = \mathbb{E}_{\omega}[\psi(\mathbb{E}_{\xi}[\varphi(X(\omega, \xi, T))])] - \frac{1}{N} \sum_{i=1}^{N} \psi \left( \frac{1}{M} \sum_{j=1}^{M} \varphi(X_{n,h}^{i,j}(\omega, \xi, T)) \right).$$

**Theorem 4.12.** Let $\varphi \in C_b^{1,\alpha}(\mathbb{R}^d, \mathbb{R}^{d'})$ and $\psi \in C_b^1(\mathbb{R}^{d'}, \mathbb{R}^{d''})$ for some $d', d'' \in \mathbb{N}$ and $0 < \alpha < 1$.

1. Let assumption 4.2 hold with the same $\alpha$ as above, then there exists a constant $C$ independent of $h, M, N$, and $\Delta t$ such that

$$\|Er\|_{L^2_{\Omega \times \Omega'}} \le C \left( (\Delta t)^{\frac{\alpha}{2}} + h^{\alpha} |\ln(h)| + \frac{1}{\sqrt{M}} + \frac{1}{\sqrt{N}} \right).$$

2. Let assumption 4.3 hold with the same $\alpha$ as above, then there exists a constant $C$ independent of $h, M, N$, and $\Delta t$ such that

$$\|Er\|_{L^2_{\Omega \times \Omega'}} \le C \left( (\Delta t)^{\frac{1+\alpha}{2}} + h |\ln(h)| + \frac{1}{\sqrt{M}} + \frac{1}{\sqrt{N}} \right).$$

**Remark 4.13.** An estimate of the error on the spread as defined in Section 1 follows from the cases where $\varphi(x) = xx^t$, $\psi(x) = x$ and $\varphi(x) = x$, $\psi(x) = xx^t$. For simplicity, we treat only the case where $\varphi$ and $\psi$ are bounded with bounded derivatives. The extension to the case where $\psi$ has polynomial growth is straightforward.
Proof. 1. Let Assumption 4.2 hold. We split the error into three terms:

$$Er(\omega, \xi) = Er1 + Er2(\omega) + Er3(\omega, \xi),$$

where we define:

$$ \begin{aligned} Er1 &= \mathbb{E}_{\omega}[\psi(\mathbb{E}_{\xi}[\varphi(X(\omega, \xi, T))])] - \mathbb{E}_{\omega}[\psi(\mathbb{E}_{\xi}[\varphi(X_{n,h}^{i,j}(\omega, \xi, T))])] \\ Er2(\omega) &= \mathbb{E}_{\omega}[\psi(\mathbb{E}_{\xi}[\varphi(X_{n,h}^{i,j}(\omega, \xi, T))])] - \frac{1}{N} \sum_{i=1}^{N} \psi(\mathbb{E}_{\xi}[\varphi(X_{n,h}^{i,j}(\omega, \xi, T))]) \\ Er3(\omega, \xi) &= \frac{1}{N} \sum_{i=1}^{N} \left( \psi(\mathbb{E}_{\xi}[\varphi(X_{n,h}^{i,j}(\omega, \xi, T))]) - \psi\left(\frac{1}{M} \sum_{j=1}^{M} \varphi(X_{n,h}^{i,j}(\xi, T))\right) \right). \end{aligned} $$

The first error term $Er1$ takes account for both the space discretization and the time discretization. For almost all $\omega$ and for any $1 \ge i \ge N$, we have,

$$ \|\mathbb{E}_{\xi}[\varphi(X^i(\omega, \xi, T))] - \mathbb{E}_{\xi}[\varphi(X_{n,h}^i(\omega, \xi, T))]\| \le P_3(\|v^i(\omega)\|_{C_b^{0,\alpha}(\mathbb{R}^d)}, \|\varphi\|_{C_b^{1,\alpha}(\mathbb{R}^d)}) \\ \times [(\Delta t)^{\frac{\alpha}{2}} + \|(v^i - v_h^i)(\omega)\|_{L^\infty(\mathbb{R}^d)}] $$

where we have used a straightforward extension of Proposition 4.9 to the case of a test function $\varphi$ with vectorial values. This inequality holds for almost all $\omega$, then by taking the expected value of the image by $\psi$ and by using Proposition 4.6 with $\tilde{v} = v_h$, we obtain thanks to Hölder inequality the existence of a constant $C_4$ such that:

$$ \begin{aligned} \|Er1\| &\le \mathbb{E}_{\omega}[\|D\psi\|_{C_b^0(\mathbb{R}^{d'}, \mathcal{L}(\mathbb{R}^{d'}, \mathbb{R}^{d''}))} P_3(\|v^i(\omega)\|_{C_b^{0,\alpha}(\mathbb{R}^d)}, \|\varphi\|_{C_b^{1,\alpha}(\mathbb{R}^d)})((\Delta t)^{\frac{\alpha}{2}} + \|(v^i - v_h^i)(\omega)\|_{L^\infty(\mathbb{R}^d)})] \\ &\le \|D\psi\|_{C_b^0(\mathbb{R}^{d'}, \mathcal{L}(\mathbb{R}^{d'}, \mathbb{R}^{d''}))} \|P_3(\|v^i(\omega)\|_{C_b^{0,\alpha}(\mathbb{R}^d)}, \|\varphi\|_{C_b^{1,\alpha}(\mathbb{R}^d)})\|_{L^2}(\|\Delta t\|^{\frac{\alpha}{2}} + C_1(2)h^\alpha |\ln h|) \\ &\le C_4((\Delta t)^{\frac{\alpha}{2}} + C_1(2)h^\alpha |\ln h|), \end{aligned} $$

where we have used the fact that $P_3$ is a polynomial function, together with the fact that $\|v\|_{C_b^{0,\alpha}(\mathbb{R}^d)}$ belongs to $L^q(\Omega)$ for any $1 \le q < +\infty$ (see Assumption 4.2 and Proposition 4.5). The random variables $(Y_i)_{1 \le i \le N}$ defined by $Y_i(\omega) = \psi(\mathbb{E}_{\xi}[\varphi(X_{n,h}^{i,j}(\omega, \xi, T))])$ being independent, identically distributed and belonging to $L^2(\Omega)$, we have:

$$ \begin{aligned} \|Er2(\omega)\|_{L^2_{\omega}} &\le \frac{\|Y_i - \mathbb{E}[Y_i]\|_{L^2_{\omega}}}{\sqrt{N}} \\ &\le \frac{2\|Y_i\|_{L^2_{\omega}}}{\sqrt{N}} \\ &\le \frac{\|\psi\|_{C_b^0(\mathbb{R}^{d'}, \mathbb{R}^{d''})}}{\sqrt{N}}. \end{aligned} $$

Indeed, for almost all $\omega$, we have $|Y_i(\omega)| \le \|\psi\|_{C_b^0(\mathbb{R}^{d'}, \mathbb{R}^{d''})}$.

Analogously, for any $1 \le i \le N$ and almost all $\omega$, the random variables $(Z_j)_{1 \le j \le M}$ defined by $Z_j(\xi) =$
$\mathbb{E}_{\xi}[\varphi(X_{n,h}^{i,j}(\omega, \xi, T))]$ are independent, identically distributed $L^2(\Omega')$ random variables, therefore we get:

$$
\begin{aligned}
\left\| \mathbb{E}[Z_j] - \frac{1}{M} \sum_{j=1}^{M} Z_j(\xi) \right\|_{L_{\xi}^2} &\le \frac{\|Z_j - \mathbb{E}[Z_j]\|_{L_{\xi}^2}}{\sqrt{M}} \\
&\le \frac{2\|Z_j\|_{L_{\xi}^2}}{\sqrt{M}} \\
&\le \frac{2\|\varphi\|_{C_b^0(\mathbb{R}^d, \mathbb{R}^{d'})}}{\sqrt{M}}.
\end{aligned}
$$

For all $1 \le i \le N$ and almost all $\omega$,

$$
\begin{aligned}
& \left\| \psi(\mathbb{E}_{\xi}[\varphi(X_{n,h}^{i,j}(\omega, \xi, T))]) - \psi\left(\frac{1}{M} \sum_{j=1}^{M} \varphi(X_{n,h}^{i,j}(\omega, \xi, T))\right) \right\| \\
& \le \|D\psi\|_{C_b^0(\mathbb{R}^{d'}, \mathcal{L}(\mathbb{R}^{d'}, \mathbb{R}^{d''}))} \left\| \mathbb{E}_{\xi}[\varphi(X_{n,h}^{i,j}(\omega, \xi, T))] - \frac{1}{M} \sum_{j=1}^{M} \varphi(X_{n,h}^{i,j}(\omega, \xi, T)) \right\|,
\end{aligned}
$$

thus

$$
\left\| \psi(\mathbb{E}_{\xi}[\varphi(X_{n,h}^{i,j}(\omega, \xi, T))]) - \psi\left(\frac{1}{M} \sum_{j=1}^{M} \varphi(X_{n,h}^{i,j}(\omega, \xi, T))\right) \right\|_{L_{\xi}^2} \le \frac{2\|\psi\|_{C_b^1(\mathbb{R}^{d'}, \mathbb{R}^{d''})} \|\varphi\|_{C_b^0(\mathbb{R}^d, \mathbb{R}^p)}}{\sqrt{M}}.
$$

This bound holds for any $1 \le i \le N$ and almost all $\omega$, therefore taking the sum over $i$ and the $L_{\omega}^2$ norm yields finally the following bound for Er3:

$$
\| \mathrm{Er3}(\omega, \xi) \|_{L_{\omega}^2 L_{\xi}^2} \le \frac{2 \| \psi \|_{C_b^1(\mathbb{R}^{d'}, \mathbb{R}^{d''})} \| \varphi \|_{C_b^0(\mathbb{R}^d, \mathbb{R}^p)}}{\sqrt{M}}.
$$

2. The case where Assumption 4.3 holds is totally similar, except the fact that we use the second parts of Propositions 4.9 and 4.6 instead of their first parts.

<figure><img src="image_3.png" /></figure>

4.7 Total error on the generalized macro-dispersion

The stochastic process $X(\omega, \xi, t)$ and its approximations $X_{n,h}^{i,j}(\omega, \xi, t)$ are defined as previously by respectively (4.21) and (4.22). Here we first define for $\bar{\varphi} : \mathbb{R}^d \mapsto \mathbb{R}^{d'}$ and $\bar{\psi} : \mathbb{R}^{d'} \mapsto \mathbb{R}^{d''}$ a quantity of interest which generalizes the macro-dispersion, namely the quantity

$$
\frac{d}{dt} \mathbb{E}_{\omega} [\bar{\psi}(\mathbb{E}_{\xi}[\bar{\varphi}(X(\omega, \xi, T))])],
$$

which is the time derivative of the quantity of interest considered in the previous section (which we called the
generalized spread). If we suppose that the test functions $\varphi$ and $\psi$ are smooth with bounded derivatives, an
application of Itô lemma yields the following equality : for almost all $\omega$ we have that $t \mapsto \mathbb{E}_{\xi}[\bar{\varphi}(X(\omega, \xi, T))]$
is continuously differentiable, its differential being

$$
t \mapsto \mathbb{E}_{\xi}[D\bar{\varphi}(X(\omega, \xi, T)).v(X(\omega, \xi, T)) + D\Delta\bar{\varphi}(X(\omega, \xi, T))].
$$

Thus for almost all $\omega$,

$t \mapsto \bar{\psi}(\mathbb{E}_{\xi}[\bar{\varphi}(X(\omega, \xi, T))])$
is continuously differentiable and its differential

$$t \mapsto D\bar{\psi}(E_{\xi}[\bar{\varphi}(X(\omega, \xi, T)))].\mathbb{E}_{\xi}[D\bar{\varphi}(X(\omega, \xi, T)).v(X(\omega, \xi, T)) + D\Delta\bar{\varphi}(X(\omega, \xi, T))]$$

can be bounded, uniformly with respect to $t$, by

$$\|\bar{\psi}\|_{C_b^1(\mathbb{R}^{d'}, \mathcal{L}(\mathbb{R}^{d'}, \mathbb{R}^{d''}))} (\|\bar{\varphi}\|_{C_b^1(\mathbb{R}^d, \mathcal{L}(\mathbb{R}^d, \mathbb{R}^{d'}))} \|v(\omega)\|_{C_b^0(\mathbb{R}^d)} + D \|\bar{\varphi}\|_{C_b^2(\mathbb{R}^d), \mathcal{B}(\mathbb{R}^d \times \mathbb{R}^d, \mathbb{R}^{d'})}),$$

which belongs to $L^1(\Omega)$. Therefore we deduce that

$$\frac{d}{dt} \mathbb{E}_{\omega}[\bar{\psi}(\mathbb{E}_{\xi}[\bar{\varphi}(X(\omega, \xi, T))])] = \mathbb{E}_{\omega}[D\bar{\psi}(E_{\xi}[\bar{\varphi}(X(\omega, \xi, T)))].\mathbb{E}_{\xi}[D\bar{\varphi}(X(\omega, \xi, T)).v(X(\omega, \xi, T)) + D\Delta\bar{\varphi}(X(\omega, \xi, T))]],$$

which leads naturally to the following approximation, using the same ideas as in the approximation of the generalized spread (see the previous subsection) :

$$\frac{1}{N} \sum_{i=1}^{N} \left[ D\bar{\psi} \left( \frac{1}{M} \sum_{j=1}^{M} \bar{\varphi}(X_{n,h}^{i,j}(T)) \right) \cdot \left( \frac{1}{M} \sum_{j=1}^{M} (D\bar{\varphi}(\tilde{X}_{n}^{i,j}(T))).v_{h}^{i}(X_{n,h}^{i,j}(T)) + D\Delta\bar{\varphi}(X_{n,h}^{i,j}(T)) \right) \right],$$

where here and below we have omitted the dependence on $\omega$ and $\xi$ for the sake of readability. Note that $\Delta\bar{\varphi}$ denotes the vector whose coordinate $(\Delta\bar{\varphi})_i$ is the Laplacian of the coordinate $\bar{\varphi}_i$ of $\bar{\varphi}$. In this subsection we give a bound for the error between the quantity and its approximation defined above.

**Definition 4.14.** The error is then defined by

$$\begin{aligned} \bar{E}r(\omega, \xi) &= \frac{d}{dt} \mathbb{E}_{\omega}[\bar{\psi}(\mathbb{E}_{\xi}[\bar{\varphi}(X(T))])] \\ &\quad - \frac{1}{N} \sum_{i=1}^{N} \left[ D\bar{\psi} \left( \frac{1}{M} \sum_{j=1}^{M} \bar{\varphi}(X_{n,h}^{i,j}(T)) \right) \cdot \frac{1}{M} \sum_{j=1}^{M} (D\bar{\varphi}(X_{n,h}^{i,j}(T))).v_{h}^{i}(X_{n,h}^{i,j}(T)) + D\Delta\bar{\varphi}(X_{n,h}^{i,j}(T)) \right] \\ &= \mathbb{E}_{\omega}[D\bar{\psi}(E_{\xi}[\bar{\varphi}(X(T))]) \cdot \mathbb{E}_{\xi}[D\bar{\varphi}(X(T)).v(X(T)) + D\Delta\bar{\varphi}(X(T))] \\ &\quad - \frac{1}{N} \sum_{i=1}^{N} \left[ D\bar{\psi} \left( \frac{1}{M} \sum_{j=1}^{M} \bar{\varphi}(X_{n,h}^{i,j}(T)) \right) \cdot \frac{1}{M} \sum_{j=1}^{M} (D\bar{\varphi}(X_{n,h}^{i,j}(T))).v_{h}^{i}(X_{n,h}^{i,j}(T)) + D\Delta\bar{\varphi}(X_{n,h}^{i,j}(T)) \right]. \end{aligned}$$

We have then the following bound for this error on the generalized macro-dispersion.

**Theorem 4.15.** Let $\bar{\varphi} \in C_b^{3,\alpha}(\mathbb{R}^d, \mathbb{R}^{d'})$ and $\bar{\psi} \in C_b^2(\mathbb{R}^{d'}, \mathbb{R}^{d''})$ for some $d', d'' \in \mathbb{N}$ and $0 < \alpha < 1$ such that Assumption 4.3 holds. There exists a constant $c$ independent of $h, M, N$ and $\Delta t$ such that

$$\|\bar{E}r(\omega, \xi)\|_{L_{\omega}^2} \le C \left( (\Delta t)^{\frac{1+\alpha}{2}} + h|\ln(h)| + \frac{1}{\sqrt{M}} + \frac{1}{\sqrt{N}} \right)$$
Proof. We split the error $\bar{E}r$ into two terms.

$$
\begin{aligned}
\bar{E}r(\omega, \xi) &= \mathbb{E}_{\omega}[D\bar{\psi}(E_{\xi}[\bar{\varphi}(X(T))]) \cdot \mathbb{E}_{\xi}[D\bar{\varphi}(X(T)) \cdot v(X(T)) + D\Delta\bar{\varphi}(X(T))]] \\
&\quad - \frac{1}{N} \sum_{i=1}^{N} \left[ D\bar{\psi} \left( \frac{1}{M} \sum_{j=1}^{M} \bar{\varphi}(X_{n,h}^{i,j}(T)) \right) \cdot \frac{1}{M} \sum_{j=1}^{M} (D\bar{\varphi}(X_{n,h}^{i,j}(T))) \cdot v_h^i(X_{n,h}^{i,j}(T)) + D\Delta\bar{\varphi}(X_{n,h}^{i,j}(T)) \right] \\
&= \mathbb{E}_{\omega}[D\bar{\psi}(E_{\xi}[\bar{\varphi}(X(T))]) \cdot \mathbb{E}_{\xi}[D\bar{\varphi}(X(T)) \cdot v(X(T)) + D\Delta\bar{\varphi}(X(T))]] \\
&\quad - \frac{1}{N} \sum_{i=1}^{N} \left[ D\bar{\psi} \left( \frac{1}{M} \sum_{j=1}^{M} \bar{\varphi}(X_{n,h}^{i,j}(T)) \right) \cdot \frac{1}{M} \sum_{j=1}^{M} (D\bar{\varphi}(X_{n,h}^{i,j}(T))) \cdot v^i(X_{n,h}^{i,j}(T)) + D\Delta\bar{\varphi}(X_{n,h}^{i,j}(T)) \right] \\
&\quad + \frac{1}{N} \sum_{i=1}^{N} \left[ D\bar{\psi} \left( \frac{1}{M} \sum_{j=1}^{M} \bar{\varphi}(X_{n,h}^{i,j}(T)) \right) \cdot \frac{1}{M} \sum_{j=1}^{M} (D\bar{\varphi}(X_{n,h}^{i,j}(T))) \cdot v^i(X_{n,h}^{i,j}(T)) + D\Delta\bar{\varphi}(X_{n,h}^{i,j}(T)) \right] \\
&\quad - \frac{1}{N} \sum_{i=1}^{N} \left[ D\bar{\psi} \left( \frac{1}{M} \sum_{j=1}^{M} \bar{\varphi}(X_{n,h}^{i,j}(T)) \right) \cdot \frac{1}{M} \sum_{j=1}^{M} (D\bar{\varphi}(X_{n,h}^{i,j}(T))) \cdot v_h^i(X_{n,h}^{i,j}(T)) + D\Delta\bar{\varphi}(X_{n,h}^{i,j}(T)) \right]
\end{aligned}
$$

The difference between the first two terms can be bounded thanks to a variant of the second part of Theorem 4.12 by taking $\varphi$ and $\psi$ defined respectively by $\varphi(v, x) = (\varphi_1(x) := \bar{\varphi}(x), \varphi_2(v, x) := D\bar{\varphi}(x) \cdot v(x) + D\Delta\bar{\varphi}(x))$ and $\psi(x, y) = D\bar{\psi}(x) \cdot y$. However we have to adapt the proof of the second point of Theorem 4.12 to bound the difference between the first two terms since $\varphi$ also depend on $v$ (and moreover $\psi$ and $D\psi$ are not bounded).

$$
\begin{aligned}
& \left\| \mathbb{E}_{\omega}[D\bar{\psi}(E_{\xi}[\bar{\varphi}(X(T))]) \cdot \mathbb{E}_{\xi}[D\bar{\varphi}(X(T)) \cdot v(X(T)) + D\Delta\bar{\varphi}(X(T))]] \right. \\
& \qquad \left. - \frac{1}{N} \sum_{i=1}^{N} \left[ D\bar{\psi} \left( \frac{1}{M} \sum_{j=1}^{M} \bar{\varphi}(X_{n,h}^{i,j}(T)) \right) \cdot \frac{1}{M} \sum_{j=1}^{M} (D\bar{\varphi}(X_{n,h}^{i,j}(T))) \cdot v^i(X_{n,h}^{i,j}(T)) + D\Delta\bar{\varphi}(X_{n,h}^{i,j}(T)) \right] \right\|_{L_{\omega, \xi}^2} \\
& = \left\| \mathbb{E}_{\omega}[\psi(\mathbb{E}_{\xi}[\varphi(v^i, X^i(T))))] - \frac{1}{N} \sum_{i=1}^{N} \psi \left( \frac{1}{M} \sum_{j=1}^{M} \varphi(v^i, X_{n,h}^{i,j}(T)) \right) \right\|_{L_{\omega, \xi}^2}
\end{aligned}
\quad (4.23)
$$

In order to bound this term, we follow the proof of Theorem 4.12. First we split the error (4.23) into three terms :

$$
\left\| \mathbb{E}_{\omega}[\psi(\mathbb{E}_{\xi}[\varphi(v^i, X^i(T))))] - \frac{1}{N} \sum_{i=1}^{N} \psi \left( \frac{1}{M} \sum_{j=1}^{M} \varphi(v^i, X_{n,h}^{i,j}(T)) \right) \right\|_{L_{\omega, \xi}^2} = \bar{E}r1 + \bar{E}r2(\omega) + \bar{E}r3(\omega, \xi),
$$

where we define

$$
\begin{aligned}
\bar{E}r1 &= \mathbb{E}_{\omega}[\psi(\mathbb{E}_{\xi}[\varphi(v(\omega), X(\omega, \xi, T))])] - \mathbb{E}_{\omega}[\psi(\mathbb{E}_{\xi}[\varphi(v^i(\omega), X_{n,h}^{i,j}(\omega, \xi, T))])] \\
\bar{E}r2(\omega) &= \mathbb{E}_{\omega}[\psi(\mathbb{E}_{\xi}[\varphi(v^i(\omega), X_{n,h}^{i,j}(\omega, \xi, T))])] - \frac{1}{N} \sum_{i=1}^{N} \psi(\mathbb{E}_{\xi}[\varphi(v^i(\omega), X_{n,h}^{i,j}(\omega, \xi, T))]) \\
\bar{E}r3(\omega, \xi) &= \frac{1}{N} \sum_{i=1}^{N} \left( \psi(\mathbb{E}_{\xi}[\varphi(v^i(\omega), X_{n,h}^{i,j}(\omega, \xi, T))]) - \psi\left(\frac{1}{M} \sum_{j=1}^{M} \varphi(v^i(\omega), X_{n,h}^{i,j}(\omega, T))\right) \right).
\end{aligned}
$$
In order to bound $\bar{E}r1$ we first notice that for any $1 \le i \le N$ and almost all $\omega$ we have

$$ \|\mathbb{E}_{\xi}[\varphi(v^i(\omega), X^i(\omega, \xi, T))] - \mathbb{E}_{\xi}[\varphi(v^i(\omega), X^i_{n,h}(\omega, \xi, T))]\| \le \tilde{P}_3(\|v^i(\omega)\|_{C_b^{1,\alpha}(\mathbb{R}^d)}, \|\varphi(v^i(\omega), \cdot)\|_{C_b^{1,\alpha}(\mathbb{R}^d)}^2) \quad (4.24) \\ \times [(\Delta t)^{\frac{1+\alpha}{2}} + \|(v^i - v_h^i)(\omega)\|_{L^\infty(\mathbb{R}^d)} + \|(v^i - v_h^i)(\omega)\|_{L^\infty(\mathbb{R}^d)}^{1+\alpha}], $$

where we have used Theorem 4.9. We note that there exists a constant $C_5$ such that we have for almost all $\omega$

$$ \|\varphi(v^i(\omega), \cdot)\|_{C_b^{1,\alpha}(\mathbb{R}^d)} \le C_5(1 + \|v^i(\omega)\|_{C_b^{1,\alpha}(\mathbb{R}^d)}\|\bar{\varphi}\|_{C_b^{3,\alpha}(\mathbb{R}^d)}), \quad (4.25) $$

which implies that $\tilde{P}_3(\|v^i(\omega)\|_{C_b^{1,\alpha}(\mathbb{R}^d)}, \|\varphi(v^i(\omega), \cdot)\|_{C_b^{1,\alpha}(\mathbb{R}^d)}^2)$ belongs to $L^q(\Omega)$ for any finite $q \ge 1$. Besides this, we note that for any $(x, y) \in \mathbb{R}^{d'} \times \mathbb{R}^{d''}$ we have

$$ \|D\psi(x, y)\|_{\mathcal{L}(\mathbb{R}^{d'} \times \mathbb{R}^{d''}, \mathbb{R}^{d''})} \le \|\bar{\psi}\|_{C_b^2(\mathbb{R}^{d'}, \mathcal{B}(\mathbb{R}^{d'} \times \mathbb{R}^{d'}, \mathbb{R}^{d''}))} (1 + \|y\|). \quad (4.26) $$

We deduce from this and (4.24), (4.26) and (4.24) and Proposition 4.5 that

$$ \begin{aligned} \|\bar{E}r1\| &\le \mathbb{E}_{\omega}[\sup_{t \in [0,1]} \|D\psi\|_{\mathcal{L}(\mathbb{R}^{d'}, \mathbb{R}^{d''})}(t\varphi(v(\omega), X(\omega, \xi, T)) + (1-t)\varphi(v^i(\omega), X_{n,h}^{i,j}(\omega, \xi, T))) \\ &\quad \times \tilde{P}_3(\|v^i(\omega)\|_{C_b^{1,\alpha}(\mathbb{R}^d)}, \|\varphi(v^i(\omega), \cdot)\|_{C_b^{1,\alpha}(\mathbb{R}^d)}^2) \\ &\quad \times [(\Delta t)^{\frac{1+\alpha}{2}} + \|(v^i - v_h^i)(\omega)\|_{L^\infty(\mathbb{R}^d)} + \|(v^i - v_h^i)(\omega)\|_{L^\infty(\mathbb{R}^d)}^{1+\alpha}]] \\ &\le \mathbb{E}_{\omega}[\|\bar{\psi}\|_{C_b^2(\mathbb{R}^{d'}, \mathcal{B}(\mathbb{R}^{d'} \times \mathbb{R}^{d'}, \mathbb{R}^{d''}))}(1 + \|\varphi_2(v(\omega), X(\omega, \xi, T))\| + \|\varphi_2(v^i(\omega), X_{n,h}^{i,j}(\omega, \xi, T))\|) \\ &\quad \times \tilde{P}_3(\|v^i(\omega)\|_{C_b^{1,\alpha}(\mathbb{R}^d)}, \|\varphi(v^i(\omega), \cdot)\|_{C_b^{1,\alpha}(\mathbb{R}^d)}^2) \\ &\quad \times [(\Delta t)^{\frac{1+\alpha}{2}} + \|(v^i - v_h^i)(\omega)\|_{L^\infty(\mathbb{R}^d)} + \|(v^i - v_h^i)(\omega)\|_{L^\infty(\mathbb{R}^d)}^{1+\alpha}]] \\ &\le \mathbb{E}_{\omega}[\|\bar{\psi}\|_{C_b^2(\mathbb{R}^{d'}, \mathcal{B}(\mathbb{R}^{d'} \times \mathbb{R}^{d'}, \mathbb{R}^{d''}))}(1 + 2\|\bar{\varphi}\|_{C_b^2(\mathbb{R}^d, \mathbb{R}^{d'})}(\|v^i(\omega)\|_{C_b^0(\mathbb{R}^d)} + D)) \\ &\quad \times \tilde{P}_3(\|v^i(\omega)\|_{C_b^{1,\alpha}(\mathbb{R}^d)}, \|\varphi(v^i(\omega), \cdot)\|_{C_b^{1,\alpha}(\mathbb{R}^d)}^2) \\ &\quad \times [(\Delta t)^{\frac{1+\alpha}{2}} + \|(v^i - v_h^i)(\omega)\|_{L^\infty(\mathbb{R}^d)} + \|(v^i - v_h^i)(\omega)\|_{L^\infty(\mathbb{R}^d)}^{1+\alpha}]] \\ &\le C_6((\Delta t)^{\frac{1+\alpha}{2}} + h|\ln h|), \end{aligned} $$

where $C_6$ is a constant.

In order to bound $\bar{E}r2$ we introduce for $1 \le i \le N$ the random variables $Y_i(\omega) = \psi(\mathbb{E}_{\xi}[\varphi(v^i(\omega), X_{n,h}^{i,j}(\omega, \xi, T))])$. They are independent, identically distributed random variable in $L^2(\Omega)$, and hence

$$ \begin{aligned} \|\bar{E}r2(\omega)\|_{L^2_{\omega}} &\le \frac{\|Y_i - \mathbb{E}[Y_i]\|_{L^2_{\omega}}}{\sqrt{N}} \\ &\le \frac{2\|Y_i\|_{L^2_{\omega}}}{\sqrt{N}} \\ &\le \frac{2\|\bar{\psi}\|_{C_b^1(\mathbb{R}^{d'})} \|\mathbb{E}_{\xi}[D\bar{\varphi}(X_{n,h}^{i,j}(\omega, \xi, T)).v^i(\omega) + D\Delta\bar{\varphi}(X_{n,h}^{i,j}(\omega, \xi, T))]\|_{L^2_{\omega}}}{\sqrt{N}} \\ &\le \frac{2\|\bar{\psi}\|_{C_b^1(\mathbb{R}^{d'})} \|(\|\bar{\varphi}\|_{C_b^1(\mathbb{R}^d)} \|v^i(\omega)\|_{C_b^0(\mathbb{R}^d)} + D)\|\bar{\varphi}\|_{C_b^2(\mathbb{R}^d)}\|_{L^2_{\omega}}}{\sqrt{N}} \\ &\le \frac{C_7}{\sqrt{N}}, \end{aligned} $$
for some constant $C_7$. And finally we bound $\bar{E}r3$ by adapting the proof of Theorem 4.12 similar to what precedes : for any $1 \le i \le N$ and almost all $\omega$, the random variables $(Z_j)_{1\le j\le M}$ defined by $Z_j(\xi) = \mathbb{E}_{\xi}[\varphi(X_{n,h}^{i,j}(\omega, \xi, T))]$ are independent, identically distributed $L^2(\Omega')$ random variables, therefore we get:

$$ \left\| \mathbb{E}[Z_j] - \frac{1}{M} \sum_{j=1}^{M} Z_j(\xi) \right\|_{L_{\xi}^2} \le \frac{2 \|\bar{\varphi}\|_{C_b^2(\mathbb{R}^d, \mathbb{R}^{d'})} (1 + \|v^i(\omega)\|_{C_b^0(\mathbb{R}^d)} + D)}{\sqrt{M}} $$

For all $1 \le i \le N$ and almost all $\omega$,

$$ \begin{aligned} & \left\| \psi(\mathbb{E}_{\xi}[\varphi(X_{n,h}^{i,j}(\omega, \xi, T))]) - \psi\left(\frac{1}{M} \sum_{j=1}^{M} \varphi(X_{n,h}^{i,j}(\omega, \xi, T))\right) \right\| \\ & \le \|\bar{\psi}\|_{C_b^1(\mathbb{R}^{d'})} 2 \|\bar{\varphi}\|_{C_b^2(\mathbb{R}^d, \mathbb{R}^{d'})} (1 + \|v^i(\omega)\|_{C_b^0(\mathbb{R}^d)} + D) \left\| \mathbb{E}_{\xi}[\varphi(X_{n,h}^{i,j}(\omega, \xi, T))] - \frac{1}{M} \sum_{j=1}^{M} \varphi(X_{n,h}^{i,j}(\omega, \xi, T)) \right\|, \end{aligned} $$

thus

$$ \begin{aligned} & \left\| \psi(\mathbb{E}_{\xi}[\varphi(X_{n,h}^{i,j}(\omega, \xi, T))]) - \psi\left(\frac{1}{M} \sum_{j=1}^{M} \varphi(X_{n,h}^{i,j}(\omega, \xi, T))\right) \right\|_{L_{\xi}^2} \\ & \le \frac{4 \|\bar{\psi}\|_{C_b^1(\mathbb{R}^{d'}, \mathbb{R}^{d''})} \|\bar{\varphi}\|_{C_b^2(\mathbb{R}^d, \mathbb{R}^{d'})}^2 (1 + \|v^i(\omega)\|_{C_b^0(\mathbb{R}^d)} + D)^2}{\sqrt{M}}. \end{aligned} $$

This bound holds for any $1 \le i \le N$ and almost all $\omega$, therefore taking the sum over $i$ and the $L_{\omega}^2$ norm and using Proposition 4.5 we get finally the following bound for $\bar{E}r3$:

$$ \|\bar{E}r3(\omega, \xi)\|_{L_{\omega}^2 L_{\xi}^2} \le \frac{C_8}{\sqrt{M}}, $$

for some constant $C_8$

Moreover we can easily bound the difference between the last two terms appearing in $\bar{E}r$ as follows.

$$ \begin{aligned} & \frac{1}{N} \sum_{i=1}^{N} \left[ D \bar{\psi} \left( \frac{1}{M} \sum_{j=1}^{M} \bar{\varphi}(X_{n,h}^{i,j}(T)) \right) \cdot \frac{1}{M} \sum_{j=1}^{M} (D \bar{\varphi}(X_{n,h}^{i,j}(T))) \cdot v(X_{n,h}^{i,j}(T)) + D \Delta \bar{\varphi}(X_{n,h}^{i,j}(T)) \right] \\ & - \frac{1}{N} \sum_{i=1}^{N} \left[ D \bar{\psi} \left( \frac{1}{M} \sum_{j=1}^{M} \bar{\varphi}(X_{n,h}^{i,j}(T)) \right) \cdot \frac{1}{M} \sum_{j=1}^{M} (D \bar{\varphi}(X_{n,h}^{i,j}(T))) \cdot v_h(X_{n,h}^{i,j}(T)) + D \Delta \bar{\varphi}(X_{n,h}^{i,j}(T)) \right] \\ & = \frac{1}{N} \sum_{i=1}^{N} \left[ D \bar{\psi} \left( \frac{1}{M} \sum_{j=1}^{M} \bar{\varphi}(X_{n,h}^{i,j}(T)) \right) \cdot \frac{1}{M} \sum_{j=1}^{M} (D \bar{\varphi}(X_{n,h}^{i,j}(T))) \cdot (v(X_{n,h}^{i,j}(T)) - v_h(X_{n,h}^{i,j}(T))) \right]. \end{aligned} $$

And we have, using again Proposition 4.6

$$ \left\| \frac{1}{N} \sum_{i=1}^{N} \left[ D \bar{\psi} \left( \frac{1}{M} \sum_{j=1}^{M} \bar{\varphi}(X_{n,h}^{i,j}(T)) \right) \cdot \frac{1}{M} \sum_{j=1}^{M} (D \bar{\varphi}(X_{n,h}^{i,j}(T))) \cdot (v(X_{n,h}^{i,j}(T)) - v_h(X_{n,h}^{i,j}(T))) \right] \right\|_{L_{\omega, \xi}^2} \\ \le \|D \bar{\psi}\|_{C_b^0(\mathbb{R}^d, \mathcal{L}(\mathbb{R}^d, \mathbb{R}^{d'}))} \|D \bar{\varphi}\|_{C_b^0(\mathbb{R}^{d'}, \mathcal{L}(\mathbb{R}^{d'}, \mathbb{R}^{d''}))} \tilde{C}_1(2) h |\ln h|. $$
**Remark 4.16.** We note that a slightly different numerical method was used in [6] to compute the mean dispersion : the time derivative was computed by using the increase on a small time step $\Delta s$, leading under additional assumptions to an error bound of the form

$$ \|\bar{E}r\|_{L^2(\Omega \times \Omega')} \le C \left( \Delta t + \Delta s + h |\ln h| + \frac{1}{\sqrt{N}} + \frac{1}{\sqrt{M \Delta s}} \right). $$

For more details on this error estimate, see [17]. Note that this numerical method requires a condition of type CFL to get convergence and seems to be less efficient than the numerical studied in this paper. The numerical comparaison of these two methods will be the subject of a future work.

<figure><img src="image_4.png" /></figure>

## References

[1] I. Babuška, F. Nobile, and R. Tempone. A stochastic collocation method for elliptic partial differential equations with random input data. *SIAM J. Numer. Anal.*, 45(3):1005-1034, 2007.

[2] I. Babuška, R. Tempone, and G. E. Zouraris. Galerkin finite element approximations of stochastic elliptic partial differential equations. *SIAM J. Numer. Anal.*, 42(2):800-825, 2004.

[3] A. Bellin, P. Salandin, and A. Rinaldo. Simulation of dispersion in heterogeneous porous formations: Statistics, first-order theories, convergence of computations. *Water Resour. Res.*, 28:2211-2227, 1992.

[4] K. A. Cliffe, M. B. Giles, R. Scheichl, and A. L. Teckentrup. Multilevel Monte Carlo methods and applications to elliptic PDEs with random coefficients. *Comput. Vis. Sci.*, 14(1):3-15, 2011.

[5] G. Dagan. *Flow and transport in porous formations*. Springer Verlag. 1989.

[6] J.-R. de Dreuzy, A. Beaudoin, and J. Erhel. Asymptotic dispersion in 2d heterogeneous porous media determined by parallel numerical simulation. *Water Resourc. Res.*, 43, 2007.

[7] G. de Marsily, F. Delay, J. Goncalves, P. Renard, V. Teles, and S. Violette. Dealing with spatial heterogeneity. *Hydrogeol. J*, 13:161-183, 2005.

[8] P. Delhomme. Spatial variability and uncertainty in groundwater flow parameters, a geostatistical approach. *Water Resourc. Res.*, pages 269-280, 1979.

[9] M. Dentz, H. Kinzelbach, S. Attinger, and W. Kinzelbach. Temporal behavior of a solute cloud in a heterogeneous porous medium: 3. numerical simulations. *Water Resour. Res.*, 7:1118, 2002.

[10] X. Dongbin. Fast numerical methods for stochastic computations: a review. *Commun. Comput. Phys.*, 5(2-4):242-272, 2009.

[11] D. Zhang and L. Zhiming. An efficient, high-order perturbation approach for flow in random porous media via karhunen-loève and polynomial expansions. *J. Comput. Phys.*, 194(2):773-794, 2004.

[12] R. Freeze. A stochastic conceptual analysis of one-dimensional groundwater flow in nonuniform homogeneous media. *Water Resour. Res.*, 11:725-741, 1975.

[13] L. W. Gelhar. *Stochastic Subsurface Hydrology*. Engelwood Cliffs. 1993.

[14] R. G. Ghanem and P. D. Spanos. *Stochastic finite elements: a spectral approach*. Springer-Verlag, New York, 1991.
[15] C. Graham, Th. G. Kurtz, S. Méléard, Ph. E. Protter, M. Pulvirenti, and D. Talay. *Probabilistic models for nonlinear partial differential equations*, volume 1627 of *Lecture Notes in Mathematics*. Springer-Verlag, Berlin, 1996. Lectures given at the 1st Session and Summer School held in Montecatini Terme, May 22–30, 1995, Edited by Talay and L. Tubaro, Fondazione C.I.M.E.. [C.I.M.E. Foundation].

[16] R.J. Hoeksema and P.K. Kitanidis. Analysis of the spatial structure of properties of selected aquifers. *Water Resour. Res.*, 21:536–572, 1985.

[17] J.Charrier. Numerical analysis of the advection diffusion of a solute in random media. *INRIA Research rapport, available at http://hal.inria.fr/inria-00581244*, 2011.

[18] J.Charrier. Strong and weak error estimates for elliptic partial differential equations with random coefficients. *SIAM J. Numer. Anal.*, 50(1):216–246, 2012.

[19] J.Charrier and A.Debussche. Weak truncation error estimates for elliptic pdes with lognormal coefficients. *Stochastic Partial Differential Equations: Analysis and Computations*, 2013.

[20] I. Karatzas and S. Shreve. *Brownian motion and stochastic calculus*, volume 113 of *Graduate Texts in Mathematics*. Springer-Verlag, New York, second edition, 1991.

[21] N. V. Krylov and M. Röckner. Strong solutions of stochastic equations with singular time dependent drift. *Probab. Theory Related Fields*, 131(2):154–196, 2005.

[22] A. Lunardi. *Analytic semigroups and optimal regularity in parabolic problems*. Progress in Nonlinear Differential Equations and their Applications, 16. Birkhäuser Verlag, Basel, 1995.

[23] R. Mikulevičius and E. Platen. Rate of convergence of the Euler approximation for diffusion processes. *Math. Nachr.*, 151:233–239, 1991.

[24] S. P. Neumann. Eulerian-lagrangian theory of transport in space-time nonstationary velocity fields: Exact nonlocal formalism by conditional moments and weak approximations. *Water Resour. Res.*, 29:633–645, 1993.

[25] Y. Rubin. Stochastic modeling of macrodispersion in heterogeneous porous media. *Water Resour. Res.*, 26:133–141, 1990.

[26] P. Salandin and V. Fiorotto. Dispersion tensor evaluation in heterogeneous media for finite peclet values. *Water Resour. Res.*, 36:1449–1455, 2000.

[27] H. Schwarze, U. Jaekel, and H. Vereecken. Estimation of macrodispersion by different approximation methods for flow and transport in randomly heterogeneous media. *Transp. Porous Media*, 43:265–287, 2001.

[28] N. Suciu, C. Vamos, and K. Sabelfeld. Ergodic simulations for diffusion in random velocity fields. *Monte Carlo and Quasi-Monte Carlo methods*, pages 659–668, 2007.

[29] A. Teckentrup. Multilevel monte carlo methods for highly heterogeneous media. *Proceedings of the Winter Simulation Conference 2012*, 2012.