\documentclass[10pt,a4paper]{article}
\usepackage[utf8]{inputenc}
\usepackage[english]{babel}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{graphicx}
\graphicspath{ {images/} }
\renewcommand{\baselinestretch}{1}
\usepackage[left=2.5cm,right=2cm,top=2.5cm,bottom=2cm]{geometry}
\usepackage[lofdepth,lotdepth]{subfig}
\usepackage{setspace}
\usepackage{float}
\usepackage{breqn}
\usepackage[colorlinks,urlcolor={blue}]{hyperref}
\usepackage{cite}
\usepackage{longtable}
\date{}

\begin{document}
\begin{singlespacing}
\title{\textbf{Adapting particle methods to model the dynamics of concentration gradients and chemical reactivity under advective diffusive transport conditions}}
\author{\normalsize A. Beaudoin$^1$, J.R. de Dreuzy and S. Huberson. \\
\scriptsize $^1$Institute Pprime, T\'el\'eport 2, Boulevard Marie et Pierre Curie, BP 30179, 86962 Futuroscope Chasseneuil Cedex France.}
\end{singlespacing}

\maketitle
\begin{singlespacing}
\begin{abstract}
Concentration organization and dynamics in heterogeneous porous media are
key physical factors driving chemical reactivity. At equilibrium, reactivity
depends not only on the concentration distribution but also on concentration
gradients. As high and low values of concentration and concentration
gradient do not superpose, we derive a transport equation for the
concentration gradients and set up an adapted particle method to approach 
them numerically. Particles dynamically optimize their organization to
provide highly accurate concentration gradients. The global strategy
combining separate particle methods for the concentration and its gradient
gives optimal predictions of reactivity.
\end{abstract}
\textbf{Keywords:} \textit{transport equation \sep concentration \sep gradient \sep particle methods \sep diffusion velocity method.}\\
\end{singlespacing}

\section{Introduction}
\label{introduction}
Particle methods are frequently used to model transport processes in porous and fractured media \cite{beaudoin1} \cite{cleary} \cite{eberts} 
\cite{herrera}. When transport by advection largely dominates dispersion and diffusion processes, particle methods offer a relevant alternative to Eulerian methods \cite{fernandez} \cite{zhang} \cite{zimmermann}. Even though recent developments in Eulerian methods have reduced numerical diffusion issues \cite{lautz} \cite{russell}, particle methods remain broadly appropriate to low dispersion conditions \cite{delay} \cite{salamon}. Particles efficiently adapt to flow structures and to local diffusive 
and dispersive conditions \cite{beaudoin2} \cite{noetinger} \cite{tyagi} \cite{wen}. 
They can integrate
exchanges between high and low flow zones with potentially sharp interfaces
(e.g. fracture/matrice) \cite{roubinet1} \cite{roubinet2} as well as chemical interactions with minerals and biofilms \cite{besnard} \cite{kapellos} 
\cite{sun} \cite{tartakovsky}. Because particles naturally resolve the multi-scale diversity of the transport conditions, they are frequently used to model flows with slight density variations
processes that emerge from their collective behavior. It is the case of the 
upscaled dispersion also called macro-dispersion that results from the 
differential influence of deterministic velocity correlations and stochastic 
dispersive/diffusive processes \cite{beaudoin3} \cite{dreuzy} \cite{tompson}. 

Particle methods become more involved when the process of interest is driven
by local interactions between particles as it is typically the case for 
reactions between solute species (homogeneous reactions). As reactivity is 
often  nonlinearly sensitive to solute concentrations \cite{appelo}, particle
methods have progressively evolved from random walk methods where independent particles are tracked following a Fokker-Planck equation \cite{kinzelbach} \cite{rahbaralam} 
\cite{tyagi} to meshfree methods with interacting particles
carrying concentration properties \cite{beaudoin1} \cite{herrera} \cite{tartakovsky}. Such methods provide continuous estimates of the
concentration field that can be coupled numerically to chemical reactivity
like any other Eulerian transport method \cite{dieuleveult} \cite{steefel}
\cite{younes}.
Particle methods are relevant to model the chemical control of reactivity expressed in terms of concentrations but not the physical control expressed 
in terms of concentration gradients \cite{fernandez} \cite{rubin}. In fact, reactivity is physically controlled by the diffusive mixing of solutes of different chemical concentrations, as mixing two waters with equilibrated solute concentrations generally results in an out-of-equilibrium solution \cite{leborgne1} \cite{simoni1} \cite{simoni2}. High concentration gradients thus promote diffusion, mixing and reactivity.

Particle methods are optimized to model concentration fields and not
concentration gradients \cite{fernandez} \cite{rubin}. Whether independent 
or interacting, particles become sparser in the medium because of spreading 
and dilution following the overall decrease of concentrations. Particles 
move apart and, after some time, can no longer resolve the spatial 
variations of the velocity field, challenging any approximation of 
concentration gradients. While particle reseeding is commonly used for 
concentrations \cite{beaudoin1} \cite{cleary} \cite{eberts} \cite{herrera}, 
we propose here a complementary strategy based on solving directly for 
concentration gradients instead of deriving concentration gradients from any 
approximation of the concentration field. Similar approaches have been studied
for convective flows of slightly varying density \cite{anderson}. 

We focus here on advection-diffusion processes in hetegeneous porous media. 
From the transport equation of concentrations, we derive the transport equation
of concentration gradients and show how both differ (section 2). We propose 
adapted particle methods to solve the concentration gradients (section 3).
we finally show how these methods, the particle method for the concentration (CPM)
and the particle method for the concentration gradient (GPM),
can be combined to compute reactivity rates when controlled both by concentrations 
and concentration gradients (section 5). 

\section{Transport equation of the concentration gradient}
\label{equation}
We focus here on advective and diffusive processes common in the transport of contaminants in porous media \cite{bear2} \cite{marsily}. In order to derive 
the transport equation of the concentration
gradient, the gradient operator $\nabla$ is applied to the advection - diffusion 
equation:
\begin{equation}
\label{eq1}
\nabla \left( \frac{\partial c}{\partial t} - \nabla \cdot (D {\bf g}) + \nabla \cdot ({\bf u}c) \right) = 0
\quad \mbox{with} \quad {\bf g} = \nabla c.
\end{equation}
$c$ is the concentration $[ML^{-3}]$, $D$ is the diffusion coefficient $[L^2T^{-1}]$, $\bf g$ is the concentration gradient $[ML^{-4}]$ and ${\bf u}$
is the flow velocity $[LT^{-1}]$. As $\nabla \wedge {\bf g}$ is zero and $D$ is
constant, the second term of equation (1) is given by:  
\begin{equation}
\label{eq2}
\nabla (\nabla \cdot (D {\bf g})) = \nabla \cdot (\nabla (D {\bf g})) + \nabla \wedge 
(\nabla \wedge (D {\bf g})) 
= \nabla \cdot (\nabla (D{\bf g})).
\end{equation}
Water flow being incompressible and irrotational in porous media, the third
term of equation (1) is given by: 
\begin{eqnarray}
\label{eq3}
\nabla ( \nabla \cdot ({\bf u}c) ) &=&  
\nabla ( (\nabla c) \cdot {\bf u}) + (\nabla \cdot {\bf u}) c) 
\nonumber \\
&=& ({\bf u} \cdot \nabla) {\bf g} + ({\bf g} \cdot \nabla) {\bf u} 
\\
&=& \nabla \cdot ({\bf u} \otimes {\bf g}) + ({\bf g} \cdot \nabla) {\bf u}  
\nonumber
\end{eqnarray}
where $\otimes$ is the tensorial product. The substitution of the two simplifications 
(Eq. \ref{eq2}) and (Eq. \ref{eq3}) in equation (\ref{eq1}) yields:
\begin{equation}
\label{eq4}
\frac{\partial {\bf g}}{\partial t} - \nabla \cdot (\nabla (D {\bf g})) + 
\nabla \cdot ({\bf u} \otimes {\bf g}) + ({\bf g} \cdot \nabla) {\bf u} = 0.
\end{equation}
The second term of the previous equation describes the diffusion of the
concentration gradient and the last two terms represent the transport of the
concentration gradient by the velocity field. Three distinct effects are 
included. The first is the displacement of the gradient application point by
the local velocity. This local velocity also induces a displacement and a
distorsion of concentration iso-contours. The variation of the distance between 
two such lines results in variations of the gradient intensity. These two
effects are accounted for by the last term of equation (\ref{eq4}). The last 
step consists in writing the diffusive term of equation (\ref{eq4}) as an
advective term by means of a diffusion velocity method \cite{beaudoin1}
\cite{beaudoin2} \cite{mycek1} \cite{mycek2} \cite{ogami}:
\begin{equation}
\label{eq5}
\frac{\partial {\bf g}}{\partial t} + \nabla \cdot (({\bf u_d}+{\bf u}) \otimes {\bf g}) + 
({\bf g} \cdot \nabla) {\bf u} = 0
\end{equation}
where ${\bf u_d}$ is the diffusion velocity $[LT^{-1}]$. The transport equation
of the concentration gradient (Eq. 4) is then transformed in a purely advection
equation with a source term (Eq. 5). Particle methods are well suited to solve 
the pure advection equation \cite{cottet} \cite{koumoutsakos1} \cite{chertock}. 
It must be noticed that the gradient vector $\bf g$ is not affected by this
transformation. Actually, the effect of the diffusion is to move the 
concentration iso-contours. The advection term is thus modified in order to
transport the gradient application point accordingly. This operation ensures 
the consistency between the concentration and its gradient at any point. The
expression of $\bf u_d$ is derived by means of a simple identification between
equations (\ref{eq4}) and (\ref{eq5}):
\begin{equation}
\label{eq6}
{\bf u_d} \otimes {\bf g} = - \nabla (D {\bf g}).
\end{equation}
It must be pointed out that the Particle Strength Exchange method can be
considered as a valuable alternative to the Velocity Diffusion Method
\cite{degond}. Both approaches have their advantages and drawbacks
\cite{zimmermann} \cite{beaudoin1}. In our case, the main advantages of the
velocity diffusion method is the straightforward inclusion of the
heterogeneity of porous media in the formulation and the ability of the
method to account for unbounded flows without any additional equation.
Respective avantages of both method should be further discussed for selecting
optimal strategies. Although the Diffusion Velocity Method was introduced more
than thirty years ago, there are still only few mathematical studies available
\cite{choquin}. In \cite{lions}, the convergence of the method was shown for
the case of the linear heat equation. Therefore, in the absence of any other
theoretical work, the problem arising when the Hessian is zero has been assumed 
to be a numerical difficulty for which different solutions can be applied (see
for example \cite{strickland}). Here, the simplest method was used by imposing 
a lower bound to the Hessian. 

\section{Particle method for solving the concentration gradient equation}
\label{model} 

In this section, we develop a particle method to solve the transport equation
for the concentration gradient (Eq. 5). We show how the diffusion velocity is
estimated, and how concentrations are derived from concentration gradients.

\subsection{Particles carrying the concentration gradient}
\label{particles}
We build up an initial set of particles carrying the concentration gradient
on the basis of a regular grid to define the initial values of the components of ${\bf g}$ \cite{cottet} \cite{koumoutsakos1}. Each particle
${\cal P}_i$ is defined by its location ${\bf x}_i$ $[L]$ corresponding to
the center of one grid cell and its gradient vector ${\bf g}_i$ $[ML^{-2}]$
in which each component corresponds to its quantity contained in the grid 
cell:
\begin{equation}
{\bf g_i} = \int_{|{\cal P}_i|} {\bf g} ({\bf x'}) d {\bf x'}
\approx
{\bf g} ({\bf x_i}) |{\cal P}_i|
\end{equation}
where $|{\cal P}_i|$ is the support of the particle ${\cal P}_i$ $[L^{2}]$.
When a regularly spaced cartesian grid is used, the result is a set of 
equidistant particles with different weights according to their location on the grid. The concentration gradient can be approximated by:
\begin{equation}
{\bf g}_h = \sum_{i} {\bf g_i} \zeta_{\epsilon} \left( {\bf x} - {\bf x}_{i} \right)
\quad \mbox{with} \quad \zeta_{\epsilon} ({\bf x}) = (1/{\epsilon}^{d}) \zeta 
({\bf x}/\epsilon)
\end{equation} 
where ${\bf g}_h$ is the approximated concentration gradient, $\zeta_{\epsilon}$ is a smoothing function and $d$ is the Euclidean
dimension of the problem. In this work, a 2D Gaussian smoothing function 
was selected for the  two-dimensional examples later shown:
\begin{equation}
\zeta_{\epsilon} ({\bf x}) = \frac{1}{\pi {\epsilon}^{2}} \exp 
\left(
- \frac{\mid {\bf x} \mid ^2}{\epsilon ^2}
\right).
\end{equation}
The smoothing parameter $\epsilon$ $[L]$ is proportional to the grid cell
size $h$ $[L]$ used to generate the particles, $\epsilon / h$ is held constant
\cite{beaudoin1} \cite{cottet} \cite{zimmermann}. To solve the transport equation for the
concentration gradient (Eq. \ref{eq5}), it is written in a Lagrangian framework
yielding the discrete form:
\begin{eqnarray}
\frac{d {\bf x}_i}{dt} = {\bf u} ({\bf x}_i) + {\bf u_d} ({\bf x}_i)  
\quad \mbox{and} \quad
\frac{d {\bf g}_{i}}{dt} = - {\bf g}_{i} \nabla {\bf u} ({\bf x}_i).
\end{eqnarray}
This set of differential equations is numerically solved by using a 4th order accurate Runge - Kutta scheme \cite{beaudoin1} \cite{mayer} \cite{zimmermann}. The time step
$\delta t$ is fixed according to the CFL conditions $\delta t = \mbox{min} (h/|{\bf u}|, \sqrt{h^2 / D})$ $[T]$ \cite{beaudoin1} \cite{moura} 
\cite{zimmermann}. 
This condition is probably not optimal since it has been shown, at least for
the case of a linear convection equation, that a so-called "Lagrangian CFL
condition" should be a better choice \cite{cottet2}. However, the former
condition, although somewhat conservative, was prefered in our work because
of the presence of an additional diffusion term. Subjected to the effects of
flow and diffusion velocities, the set of particles becomes progressively non
uniform. In order to maintain the overlap relation $\epsilon / h$ = constant,
a classical regridding process is used here all the $N_{\delta t}$ time steps
\cite{beaudoin1} \cite{cottet} \cite{zimmermann}. 

\subsection{Estimation of the diffusion velocity}
\label{diffusion}
To express the diffusion velocity ${\bf u_d}$, the divergence operator $\nabla
\cdot$ is applied to equation (\ref{eq6}):
\begin{equation}
\label{eq10}
\nabla \cdot ({\bf u_d} \otimes {\bf g}) = - \nabla \cdot (\nabla (D{\bf g})). 
\end{equation} 
As the diffusion coefficient $D$ is constant, the right hand side reads:
\begin{equation}
\nabla \cdot ({\bf u_d} \otimes {\bf g}) = - D \Delta {\bf g}.
\end{equation}
From the definition of the tensorial product, the previous equation can be rewritten by:
\begin{equation}
({\bf u_d} \cdot \nabla) {\bf g} = - D \Delta {\bf g}.
\end{equation}
The components of the diffusion velocity ${\bf u_d}$ are determined by solving
the system of equations (13) by means of Cramer's rule \cite{brunetti}:
\begin{equation}
\left \{
\begin{array}{l}
u_{dx} = 
- D \left( 
\Delta g_x \frac{\partial g_y}{\partial y}
-
\Delta g_y \frac{\partial g_x}{\partial y}
\right) 
/
\left( 
\frac{\partial g_x}{\partial x} \frac{\partial g_y}{\partial y} - 
\frac{\partial g_x}{\partial y} \frac{\partial g_y}{\partial x} 
\right)
\\
u_{dy} = 
-D \left(
\Delta g_y
\frac{\partial g_x}{\partial x}
-
\Delta g_x
\frac{\partial g_y}{\partial x}
\right)
/
\left(
\frac{\partial g_x}{\partial x} \frac{\partial g_y}{\partial y} - 
\frac{\partial g_x}{\partial y} \frac{\partial g_y}{\partial x}
\right)
\end{array}
\right .
\end{equation}
where ($u_{dx}$,$u_{dy}$) and ($g_x$,$g_y$) are respectively the components
of the diffusion velocity ${\bf u_d}$ and of the concentration gradient 
${\bf g}$.  

\subsection{Derivation of the concentration from its gradient}
\label{concentration}
The computation of the concentration from particles carrying the concentration
gradient ${\bf g}$ can not be performed by means of a direct convolution product
because the relation between $c$ and $\bf g$ is a first order partial 
differential equation. We derive the concentration by solving the Poisson equation giving the concentration $c$ 
from its gradient $\bf g$ obtained from equation (1 - right):
\begin{equation}
\label{eq7}
\Delta c = \nabla \cdot {\bf g}.
\end{equation}
An integral solution of which, for an unbounded domain with zero external conditions, reads:
\begin{equation}
\label{eq8}
c ({\bf x}) = \int_{{\cal R}^2} G({\bf x},{\bf x'}) \nabla \cdot {\bf g} ({\bf x'}) 
d {\bf x'}
\quad \mbox{with} \quad 
G ({\bf x}) = \frac{1}{2 \pi} \mbox{log} \mid {\bf x} \mid
\quad \mbox{in 2D}
\end{equation}
where $G$ is the Green function of Poisson equation \cite{hancock} \cite{hunt}. In this paper, the boundary conditions for the concentration were either
unbounded domain conditions or periodic conditions. The flow was actually
bounded in two directions because of the use of a finite volume solver. The
absence of boundaries was simulated by using a flow domain much larger than 
the domain occupied by the particles. The velocity diffusion method only
requires the knowledge of the flow velocity and concentration within this 
last domain.     

\section{Concentration versus gradient particle methods to simulate reactive transport in heterogeneous porous media}
We show how the two particle methods can be used to simulate reactive transport
in heterogeneous porous media. We first describe the reactive transport model
and the medium properties. We then compare the accuracy and convergence of
the concentration (CPM) and gradient (GPM) particle methods. We eventually
propose a more optimal combination of both particle methods to improve the
estimate of reactivity.

\subsection{Chemical reactivity rate}
The simple case of a dissolution/precipitation reaction at equilibrium is considered here. The precipitate component $P$ is in equilibrium with the
solutes $A$ and $B$ of concentrations $c_A$ and $c_B$ $(A+B \Leftrightarrow P)$
\cite{appelo}. Assuming that $P$ is in a pure phase, its activity is equal to
one, while activities of solute species are equal to their concentrations. 
$c_A$ and $c_B$ are then related by:
\begin{equation}
c_{A} \cdot c_{B} = K
\end{equation}
with $K$ the equilibrium constant $[ML^{-3}]$. The reactive transport equation of
each specy can be expressed generically for $c$ (either $c_A$ or $c_B$) as:
\begin{equation}
\frac{\partial c}{\partial t}  - \nabla \cdot (D \nabla c) + 
\nabla \cdot ({\bf u}c) = r.
\end{equation}
The reactivity rate $r$ $[ML^{-3}T^{-1}]$ is derived by combining the
equilibrium equation (17) and the reactive transport equations (18) for 
$c_A$ and $c_B$:
\begin{equation}
r = \frac{2K}{((\delta c)^2 + K)^{3/2}} 
\nabla (\delta c) D \nabla (\delta c)
\end{equation}
where $\delta c = c_A - c_B$ is the difference of concentrations between
the species $A$ and $B$ \cite{simoni1}. $\delta c$ is not modified by reactivity but only 
by transport processes. In fact, as $\delta c$ is simply obtained by subtracting equations (18) expressed for each specy, it follows the classical transport equation: 
\begin{equation}
\frac{\partial (\delta c)}{\partial t} - \nabla \cdot (D \nabla (\delta c))
+ \nabla \cdot ({\bf u} (\delta c)) = 0.
\end{equation}
The first term of the left-hand side of equation (19) comes from the chemical
system. The second term shows that reactivity is physically induced by the
local gradient of solute concentration imbalances \cite{simoni2}. This simple
example shows how mixing and reactivity are activated by interactions between
advective and dispersive processes traduced in the local-scale dynamics of
concentration gradients \cite{tanguy1} \cite{leborgne4}. 

\subsection{Benchmark description}
The dissolution/precipitation reaction is modeled in a 2D heterogeneous
porous medium defined by a rectangular domain of dimensions $L_x$ = 2048 $m$
and $L_y$ = 512 $m$ \cite{beaudoin3} \cite{gelhar} \cite{tompson}. The porous
medium is characterized by a random hydraulic conductivity field $K$ 
$[LT^{-1}]$ following a stationary log-normal probability distribution, 
parametrized by a mean $m$ and a Gaussian covariance function $C$ given by \cite{dagan} \cite{axness}:
\begin{equation}
C ({\bf r}) = \sigma^2 \exp \left( - \left( \frac{|{\bf r}|}{\lambda} \right)^2 \right)
\end{equation}
where $|{\bf r}|$ represents the separation distance between two points $[L]$.
The mean $m$, the variance $\sigma^2$ and the correlation length $\lambda$
are respectively equal to 0, 4 and 10 $m$. The generation of a correlated
log-normal field is performed via a Fourier transform by means of the
software FTTW \cite{frigo} \cite{gutjahr}. Steady state flows are described 
by mass conservation and Darcy's law \cite{bear2} \cite{marsily}:
\begin{equation}
\nabla \cdot {\bf u} = 0
\quad \mbox{and} \quad
{\bf u} = - K \nabla \phi
\end{equation}
where ${\bf u}$ and $\phi$ are the Darcean velocity $[LT^{-1}]$ and the 
hydraulic head $[L]$. The boundary conditions are permeameter-like with 
homogeneous Neumann on the upper and lower sides and Dirichlet $\phi = 0$ 
on the left side and $\phi = L_x$ on the right side. The flow equations are 
discretized on a regular grid using a Finite Volume scheme with the mesh resolution equal to 1 $m$ in the two directions $x$ and $y$ \cite{chavent}
\cite{roberts}. The resulting linear system is solved with the parallel
algebraic multigrid method of HYPRE \cite{falgout}. At initial time $t$ = 0 
$h$, the difference of concentrations between the species $A$ and $B$,
$\delta c_0$, is injected at the point source ${\bf x_o}$ = ($x_o$ = 256 $m$,
$y_o$ = 256 $m$). The set of particles is initialized by using a 
gaussian function. The 
re-griding frequency of particles $N_{\delta t}$ is fixed to $10$ time steps for 
the two particle methods. The Peclet number is defined by $P_e = \lambda |{\bf
u}_{mean}| / D$. The norm of the mean flow velocity $|{\bf u}_{mean}|$ is equal to 
1 $mh^{-1}$. Thus the Peclet number $P_e$ is equal to 100 for $D$ = 0.1 $m^2 h^{-1}$.

\subsection{Comparison of the two particle methods}

On Fig. 2 (top), concentrations $\delta c$ between the species $A$ and
$B$ obtained by the concentration (CPM) and gradient (GPM) particle
methods are broadly different. While $\delta c$ obtained by CPM remains
channelled along the velocity field as it should be \cite{goc} 
\cite{kerrou} \cite{knudby} \cite{tsang}, GPM smoothes it critically
because of the necessary integration between the concentration gradient
and the concentration itself (Eq. 16). Smoothing effects accumulate
with time eventually and remove the multiple concentration peaks at
time $t$ = 150 $h$. This effect is especially marked in the red square
zone indicated on Fig. 2.
 
Derivation of the two components $g_x$ and $g_y$ of the concentration
gradient $\bf g$ from the concentration $\delta c$ obtained with CPM
induces oscillations. These oscillations occur everywhere in the solute
plume because heterogeneity enhances concentration gradients. 

Because each particle method gives accurate numerical results for each
term of the reactivity rate $r$ (Eq. 19), we propose a strategy
combining CPM and GPM. The chemical term depending on $\delta c$ is
estimated by CPM. The physical term depending on $\nabla \delta c = 
\bf g$ is estimated by GPM. On Fig. 3, the reactivity rate $r$, 
obtained by this combination of particle methods, is presented at two
times, $t$ = 75 and 150 $h$, with a grid cell size $h$ = 0.64 $m$. No
oscillation appears on the numerical results. Spatial variability 
of the velocity field broadly induces channeling and localised 
reactivity.

The convergence of the two particle methods proposed is assessed on 
the basis of the integral $R$ of the reactivity rate $r$ over the computational domain:  
\begin{equation}
R (t) = \int_{L_x \times L_y} r ({\bf x}, t) dx dy.
\end{equation}
As shown on Fig. 4, the evolution of $R$ with CPM is not monotonous 
and thus does not converge. The two other particle methods give results 
close and display monotonous trends. The combined particle method 
especially gives very consistent results for the three particle 
resolutions showing good convergence properties. It also shows that a 
limited number of particles gives acceptable approximations of 
reactivity. 

In terms of performances, the particle number $N_p$ increases with the
same trend for the two particle methods, CPM and GPM, as the 
concentration plume spreads in the porous medium (see Fig. 5). As the 
spreading eventually scales in both directions with the square root of 
time, $N_p$ scales linearly in the intermediary asymptotic regime. 
Table 1 gives the iteration number $N_{iter}$, the time step $\delta t$
and the CPU time obtained by CPM and GPM for the three values of the
grid cell size $h$ tested here. As CPM and GPM give the same time
evolution of $N_p$ (see Fig. 5), the CPU times obtained by the two
particle methods are very close. Then an optimal solution is the
combination of the two particle methods with a grid cell size $h$ = 
1.28 $m$. This solution gives accurate predictions of reactivity with a 
reasonable CPU time. 

\section{Conclusions}
\label{conclusions}
Concentration gradients are generally estimated from concentrations on the 
basis of finite difference methods. This approach depends not only on the order
of accuracy of finite differences but also on the order of the method used for 
estimating the concentration. In this work, solving directly the transport
equation of concentration gradient has improved the description of 
concentration gradients. Obtaining this transport equation has revealed a
stretching term that takes into account a vectorial adaptation of concentration
gradients to the topology of iso-concentration lines. Meshfree particle methods
inherently conserve this characteristic. While estimates of concentration
gradients are accurate, reconstructed concentrations are too smooth as a result
of the necessary integration. As reactivity eventually depends on the
concentration and its gradient, we propose a global strategy combining separate
particle methods for the concentration and its gradient. Prospective results
show that this solution gives accurate predictions of reactivity. As the
CPU time stays reasonable, the futur work is to use this solution for studying
reactivity in higher heterogeneous porous media. The extension to 3-D cases 
will be also considered. It is straigthfoward and does not require any new
development exepted drastic improvements in the computational efficiency of 
the code.

\section*{Acknowledgments} The French National Research Agency ANR is
acknowledged for its financial founding through the H2MNO4 project 
(ANR - 12 - MONU - 0012). 

\bibliographystyle{elsarticle-num}
\begin{thebibliography}{00}
\bibitem{anderson} C.R. Anderson, A vortex method for flows with slight density
variations, J. Comput. Phys., 61 (1985), 417-444.
\bibitem{tanguy1} P. de Anna, M. Dentz, A. Tartakovsky and T. Le Borgne, The filamentary structure of mixing fronts and its control on reaction kinetics in 
porous media flows, Geophys. Res. Lett., 41 (2014) 4586-4593.
\bibitem{appelo} C. Appelo and D. Postma, Geochemistry, Groundwater 
and Pollution, AA Balkema, Brookfield, 1994.
\bibitem{bear2} J. Bear, Hydraulics of groundwater, McGraw-Hill, New
York, 1979.
\bibitem{besnard} K. Besnard, J.R. de Dreuzy, P. Davy and L. Aquilina, A modified Lagrangian-volumes method to simulate nonlinearly and kinetically sorbing solute transport in heterogeneous porous media, J. Contam. Hydrol.,
120-121 (2011) 89-98.
\bibitem{beaudoin1} A. Beaudoin, S. Huberson and E. Rivoalen, Simulation of anisotropic diffusion by means of a diffusion velocity method, J. Comput. 
Phys., 186 (2003) 122-135.
\bibitem{beaudoin3} A. Beaudoin, J. Erhel and J.R. de Dreuzy, A comparison between a direct and a multigrid sparse linear solvers for highly heterogeneous flux computations, Eccomas CFD 2006, volume CD, 2006.
\bibitem{beaudoin2} A. Beaudoin, S. Huberson and E. Rivoalen, A particle
method for solving Richard's equation, C.R. Acad. Sci., Ser. IIb: Mec., 
339 (2011) 257-261.
\bibitem{brunetti} M. Brunetti, Old and new proofs of cramer’s rule, Applied Mathematical
Sciences, 8 (2014) 6689-6697.
\bibitem{chavent} G. Chavent and J. E. Roberts, A unified physical presentation of mixed,
mixed-hybrid finite elements and standard finite difference approximations for the
determination of velocities in waterflow problems, Adv. Water Resour., 14 (1991) 329-348.
\bibitem{chertock} A. Chertock and D. Levy, Particle methods for dispersive
equations, J. Comput. Phys., 171 (2001) 708-730.
\bibitem{choquin} J.P. Choquin and S. Huberson, Particle simulation of viscous flow, Comp. and fluids,
17 (1989), p 397.
\bibitem{cleary} P.W. Cleary and J.J. Monaghan, Conduction modelling using smoothed particle hydrodynamics, J. Comput. Phys., 148 (1999) 227-264.
\bibitem{cottet} G.H. Cottet and P.D. Koumoutsakos, Vortex Methods 
Theory and Practice, Cambridge University Press, Mars 2000.
\bibitem{cottet2} G.H. Cottet, J.M. Etancelin, F. Perignon and C. Picard, High
order semi-Lagrangian particle methods for transport equations: numerical analysis
and implementation issues, ESAIM: M2AN, 48 (2014), 1029-1060.
\bibitem{dagan} G. Dagan, Solute transport in heterogeneous porous 
formation, J. Fluid Mech., 145 (1984) 151-177.
\bibitem{degond} P. Degond, S. Mas-Gallic, The weighted particle method for
convection–diffusion equations, Part 2: The anisotropic case, Math.
Comput., 53 (1989) 509.
\bibitem{delay} F. Delay, G. Porel and P. Sardini, Modelling diffusion in a
heterogeneous rock matrix with a time-domain Lagrangian method and an
inversion procedure, C.R. Acad. Sci., Ser. IIa: Sci. Terre Planets, 
334 (2002) 967-973.
%\bibitem{dentz} M. Dentz, T. Le Borgne, A. Englert and B. Bijeljic, Mixing, spreading and %reaction in heterogeneous media : A brief review, J. Contam.
%Hydrol., 120-121 (2011) 1-17.
\bibitem{dieuleveult} C. de Dieuleveult, J. Erhel and M. Kern, A global 
strategy for solving reactive transport equations, J. Comput. Phys., 228
(2009) 6395-6410.
\bibitem{dreuzy} J.R. de Dreuzy, A. Beaudoin and J. Erhel, Asymptotic dispersion in 2D heterogeneous porous media determined by parallel numerical simulations, Water Resour. Res., 43 (2007).
\bibitem{eberts} S.M. Eberts, J.K. Bohlke, L.J. Kauffman and B.C. Jurgens, Comparison of particle-tracking and lumped-parameter age-distribution models for evaluating vulnerability of production wells to contamination, Hydrogeol.
J., 20 (2012) 263-282.
\bibitem{falgout} R.D. Falgout, J.E. Jones and U.M. Yang, Pursuing
scalability for hypre’s conceptual interfaces, ACM Trans. Math. Softw., 31 (2005) 326-350.
\bibitem{fernandez} D. Fernandez-Garcia and X. Sanchez-Vila, Optimal 
reconstruction of concentrations, gradients and reaction rates from particle
distributions, J. Comput. Phys., 120-121 (2011) 99-114.
\bibitem{frigo} M. Frigo and S.G. Johnson, The design and implementation of FFTW3, paper
presented at the IEEE, 93 (Special Issue on Program Generation, Optimization, and Platform
Adaptation), 2005.
\bibitem{axness} L. Gelhar and C. Axness, Three dimensional stochastic analysis of macrodispersion in aquifers, Water Resour. Res., 19 (1983) 161-180.
\bibitem{gelhar} L. Gelhar, Stochastic subsurface hydrology, Engelwood Cliffs, New
Jersey, 1993.
\bibitem{goc} R. Le Goc, J.R. de Dreuzy and P. Davy, Statistical 
characteristics of flow as indicators of channeling in heterogeneous
porous and fractured media, Advances in Water Resources, 33 (2010) 257-269.
\bibitem{gutjahr} A.L. Gutjahr, Fast fourier transforms for random field
generation, Project Report for Los Alamos Grant to New Mexico Tech, Contract number 4-R58-2690R, Department of Mathematics, New Mexico Tech, Socorro, New Mexico, 1989.
\bibitem{hancock} M.J. Hancock, Method of Green’s functions, Linear Partial Differential
Equations, 2006
\bibitem{herrera} P.A. Herrera, M. Massabo and R.D. Beckie, A meshless method to simulate solute transport in heterogeneous porous media, Adv. Water Resour., 32 (2009) 413-429.
\bibitem{hunt} R. Hunt, Chapter 2 - Poissons Equation, Lecture Notes, University of
Cambridge, 2007.
\bibitem{huysmans} M. Huysmans and A. Dassargues, Review of the use of
Peclet numbers to determine the relative importance of advection and diffusion in low permeability environments, Hydrogeol. J., 2004.
\bibitem{kapellos} G.E. Kapellos, A.S. Terpsichori and A.C. Payatakes,
Hierarchical simulator of biofilm growth and dynamics in granular porous
materials, Adv. Water Resour., 30 (2007) 1648-1667.
\bibitem{kerrou} J. Kerrou, P. Renard, H.J. Hendricks Franssen and I. Lunati, Issues in characterizing heterogeneity and connvectivity in 
non-multiGaussian media, Adv. Water Resour., 31 (2008) 147-159.
\bibitem{kinzelbach} W. Kinzelbach, The random-walk method in pollutant transport simulation, Groundwater flow and quality modelling 1988, Springer, New York, NATO ASI Serie C, 227-245.
\bibitem{konikow1} L.F. Konikow, Applying dispersive changes to 
lagrangian particles in groundwater transport models, Transport Porous
Med., 85 (2010) 437-449.
\bibitem{konikow2} L.F. Konikow, D.J. Goode and G.Z. Hornberger, A 
Three-dimensional method of characteristics solute-transport model
(MOC3D), U.S. Geological Survey Water Resources Investigation Report, 96-4267, 1996.
\bibitem{koumoutsakos1} P. Koumoutsakos, Multiscale flow simulations
using particles, Annu. Rev. Fluid Mech., 37 (2005) 457-487.
\bibitem{knudby} C. Knudby and J. Carrera, On the relationship between
indicators of geostatistical flow and transport connectivity, Adv. Water
Resour., 28 (2005) 405-421.
\bibitem{lautz} L.K. Lautz and D. I. Siegel, Modeling surface and ground water mixing in the hyporheic zone using MODFLOW and MT3D, Adv. Water Resour., 29 (2006) 1618-1633.
\bibitem{leborgne1} T. Le Borgne, M. Dentz, D. Bolster, J. Carrera, J.R. de Dreuzy and P. Davy, Non-Fickian mixing : Temporal evolution of the scalar dissipation rate in heterogeneous porous media, Adv. Water Resour., 3 
(2010) 1468-1475.
%\bibitem{leborgne2} T. Le Borgne, M. Dentz, P. Davy, D. Bolster, J. Carrera, J.R. de Dreuzy 
%and O. Bour, Persistence of incomplete mixing : A key to anomalous transport, Phys. Rev. E, 
% 84 (2011).
\bibitem{leborgne4} T. Le Borgne, M. Dentz and E. Villermaux, Stretching, 
coalescence, and mixing in porous media, Phys. Rev. Lett., 110 (2013) 204501.
%\bibitem{leborgne3} T. Le Borgne, M. Dentz and E. Villermaux, The lamellar
%description of mixing in porous media, J. Fluid Mech., 770 (2015) 458-498.
\bibitem{lions} P.L. Lions and S. Mas-Gallic, Une m\'ethode 
particulaire d\'eterministe pour des \'equations diffusives non lin\'eaires, CRAS Sciences - S\'erie I - Mathematica, 332 (2001), 
369-376. 
\bibitem{marsily} G. de Marsily, Quantitative hydrogeology : 
groundwater hydrology for engineers, Academic press, Inc., San Diego,
California, 1986.
\bibitem{mayer} A. Mayer, Numerical approach of physical problems, 
Chapter 16 : Integration of ordinary differential equations, Lecture
Notes, University of Namur, 2014.
\bibitem{moura} C.A. de Moura and C.S. Kubrusly, The 
Courant-Friedrichs-Lewy (CFL) Condition, Birkhauser, 2013. 
\bibitem{mycek1} P. Mycek, G. Pinon, G. Germain and E. Rivoalen,
Formulation and analysis of a diffusion-velocity particle model for
transport-dispersion equations, Comput. Appl. Math., (2014) 1-27.
\bibitem{mycek2} P. Mycek, G. Pinon, G. Germain and E. Rivoalen, A 
self-regularising DVM-PSE method for the modelling of diffusion in
particle methods, C.R. Acad. Sci., Ser. IIb: Mec., 341 (2013) 9-10.
\bibitem{noetinger} B. Noetinger, D. Roubinet, A. Russian, T. Leborgne, F.
Delay, M. Dentz, J.R. de Dreuzy and P. Gouze, Random walk methods for
modeling hydrodynamic transport in porous and fractured media from pore to
reservoir scale, Transport Porous Med., (2016) 1-41.
\bibitem{ogami} Y. Ogami and T. Akamatsu, Viscous flow simulation 
using the discrete vortex model - the diffusion velocity method,
Comput. Fluids, 19 (1991) 433-444.
\bibitem{pokrajac} D. Pokrajac and R. Lazic, An efficient algorithm for 
high accuracy particle trakcing in finite elements, Adv. Water Resour., 
25 (2002) 353-369.
\bibitem{rahbaralam} M. Rahbaralam, D. Fernandez-Garcia and X. Sanchez-Vila,
Do we really need a large number of particles to simulate bimolecular 
reactive transport with random walk methods ? A kernel density estimation
approach, J. Comput. Phys., 303 (2015) 95-104.
\bibitem{roberts} J.E. Roberts and J.M. Thomas, Mixed and hybrid methods, in Handbook of Numerical Analysis 2, Finite Element Methods part 1, edited by P.G. Ciarlet and J.L. Lions, 1991, 523-639, Elsevier Sci.
\bibitem{roubinet1} D. Roubinet, H.H. Liu and J.R. de Dreuzy, A new particle-tracking approach to simulating transport in heterogeneous fractured porous media, Water Resour. Res., 46 (2010).
\bibitem{roubinet2} D. Roubinet, J.R. de Dreuzy and D.M. Tartakovsky, Particle-tracking simulations of anomalous transport in hierarchically fractured rocks, Comput. Geosci., 50 (2013) 52-58.
\bibitem{roy} C.J. Roy and M.M. Hopkins, Discretization error estimates 
using exact solutions to nearby problems, IAA Paper 2003-0629.
\bibitem{rubin} J. Rubin, Transport of reacting solutes in porous media - relation between mathematical nature of problem formulation and chemical nature of reactions, Water Resour. Res., 19 (1983) 1231-1252.
\bibitem{russell} T.F. Russell and M.A. Celia, (2002), An overview of
research on Eulerian-Lagrangian localized adjoint methods (ELLAM),
Adv. Water Resour., 25 (8-12) 1215-1231.
\bibitem{salamon} P. Salamon, D. Fernandez-Garcia and J.J Gomez-Hernandez, A review and numerical assessment of the random walk particle tracking method, J. Contam. Hydrol., 87 (2006) 277-305.
\bibitem{simoni1} M. de Simoni, J. Carrera, X. Sanchez-Villa and A. Guadagnini, A procedure for the solution of multicomponent reactive transport problems, Water Resour. Res., 41 (2005).
\bibitem{simoni2} M. de Simoni, X. Sanchez-Vila, J. Carrera, M.W. Saaltink, A mixing ratios-based formulation for multicomponent reactive transport, Water Resour. Res., 43 (2007).
\bibitem{steefel} C.I. Steefel, D.J. DePaolo and P.C. Luchtner, Reactive transport modeling : An essential tool and a new research approach for the Earth sciences, Earth Planet. Sci. Lett., 240 (2005) 539-558.
\bibitem{strickland} J.H. Strickland, S.N. Kempka and W.P. Wolf, 
Viscous diffusion of vorticity in unsteady wall layers using the
diffusion velocity concept, ESAIM Proceedings, 1 (1996), 13-151. 
\bibitem{sun} N.Z. Sun, A finite cell method for simulating the mass transport process in porous media, Water Resour. Res., 35 (1999) 3649-3662.
\bibitem{tartakovsky} A.M. Tartakovsky, P. Meakin, T.D. Scheibe and R.M. Eichler West, Simulations of reactive transport and precipitation with smoothed particle hydrodynamics, J. Comput. Phys., 222 (2007) 654-672.
\bibitem{tompson} A. Tompson and L. Gelhar, Numerical simulation of 
solute transport in three dimensional, randomly heterogeneous porous
media, Water Resour. Res., 26 (1990) 2541-2562.
\bibitem{tsang} C.F. Tsang and I. Neretnieks, Flow channeling in 
heterogeneous fractured rocks, Rev. Geophys., 36 (1998) 275-298.
\bibitem{tyagi} M. Tyagi, P. Jenny, I. Lunati and H-A. Tchelepi,
A lagrangian, stochastic modeling framework for multi-phase flow in
porous media, J. Comp. Phys., 227 (2008) 6696-6714.
\bibitem{wen} X.H. Wen and J.J. Gomez-Hernandez (1996), The constant
displacement scheme for tracking particles in heterogeneous aquifers,
Ground Water, 34 (1), 135-142.
\bibitem{wexler} E.J. Wexler, Analytical solutions for one-, two-, and 
three-dimensional solute transport in ground-water systems with uniform
flow : U.S. Geological Survey Techniques of Water-Resources Investigations, Book 3, Chapter B7, 1992.
\bibitem{younes} A. Younes and P. Ackerer, Solving the advection-diffusion
equation with the Eulerian-Lagrangian localized adjoint method on
unstructured meshes and non uniform time stepping, J. Comput. Phys.,
208 (2005) 384-402.
\bibitem{zhang} Z. Zhang and Q. Chen, Comparison of the Eulerian and Lagrangian methods for predicting particle transport in enclosed spaces, Atmos. Environ., 41 (2007) 5236-5248.
\bibitem{zheng} C. Zheng nd G.D. Bennet, Applied contaminant transport 
modeling : second edition, John Wiley and Sons, New York, 2002.
\bibitem{zimmermann} S. Zimmermann, P. Koumoutsakos and W. Zinzelbach, 
Simulation of pollutant transport using a particle method, J. Comput. Phys., 173 (2001) 322-347. 

\end{thebibliography}

\clearpage

\begin{figure}
\begin{center}
\includegraphics[width=12cm]{figure111.eps}
\includegraphics[width=12cm]{figure112.png}
\end{center}
\caption{Top : sketch of benchmark, bottom : horizontal component of flow velocity.}
\label{figure8}
\end{figure}

\clearpage

\begin{figure}
\begin{center}
\includegraphics[width=6.5cm]{concentration-concentration-nx=200-t=150.png}
\includegraphics[width=6.5cm]{concentration-gradient-nx=200-t=150.png}
\includegraphics[width=6.5cm]{dcdx-concentration-nx=200-t=150.png}
\includegraphics[width=6.5cm]{dcdx-gradient-nx=200-t=150.png}
\includegraphics[width=6.5cm]{dcdy-concentration-nx=200-t=150.png}
\includegraphics[width=6.5cm]{dcdy-gradient-nx=200-t=150.png}
\end{center}
\caption{Difference of concentration $\delta c$ (top), horizontal (middle)
and vertical (bottom) components, $g_x$ and $g_y$, of gradient $\bf g$ of
concentration difference $\delta c$ between the species $A$ and $B$ (colored
contour lines) obtained by the concentration (left) and gradient (right)
particle methods at time $t$ = 150 $h$ for a value of the grid cell size $h$
= 0.64 $m$ with a Peclet number $P_e$ = 100. Arrows show the flow field.}
\label{figure9}
\end{figure}

\clearpage

\begin{figure}
\begin{center}
\includegraphics[width=10cm]{reaction-mixed-nx=200-t=75.png}
\includegraphics[width=10cm]{reaction-mixed-nx=200-t=150.png}
\end{center}
\caption{Reactive rate $r$ (colored contour lines) obtained by the combined
particle method at $t$ = 75 $h$ (top) and 150 $h$ (bottom) for a grid cell
size $h$ = 0.64 $m$ with a Peclet number $P_e$ = 100.}
\label{figure9}
\end{figure}

\clearpage

\begin{figure}
\begin{center}
\includegraphics[width=13cm]{figure10.eps}
\end{center}
\caption{Integral $R$ of the reactivity rate $r$ over the computational domain as a function
of the grid cell size $h$ obtained by CPM (blue), GPM (green) and combined particle method
(black) at time $t$ = 150 $h$ with a Peclet number $P_e$ = 100.}
\label{figure9}
\end{figure}

\clearpage

\begin{figure}
\begin{center}
\includegraphics[width=13cm]{particle-number.eps}
\end{center}
\caption{Time evolution of particle number $N_p$ obtained by the concentration (blue) and
gradient (green) particle methods for three values of the grid cell size $h$ = 0.64 $m$ 
(dotted line), 1.28 $m$ (solid line) and 2.56 $m$ (dashed  line) with a Peclet number $P_e$ 
= 100.}
\label{figure14}
\end{figure}

\clearpage

\begin{table}
\begin{center}
\noindent
\begin{tabular}{|c|c|c|c|c|}
\hline
$h$ $(m)$ & $N_{iter}$ & $\delta t$ $(h)$ & CPU time $(h)$ & CPU time $(h)$ \\
 & & & with CPM & with GPM \\
\hline 
0.64 & 4688 & 0.032 & 11.56 & 12.87 \\
1.28 & 2344 & 0.064 & 1.55 & 1.73 \\
2.56 & 1172 & 0.128 & 0.24 & 0.27 \\
\hline
\end{tabular}
\caption{Iteration number $N_{iter}$, time step $\delta t$ and CPU time for 
a physical time of 150 $h$ obtained by the concentration (CPM) and gradient
(GPM) particle methods with the three values of the grid cell size $h$ =
0.64, 1.28 and 2.56 $m$ for a Peclet number $P_e$ = 100. The numerical
simulations were performed on a single processor Intel(R) Core(TM) i7-4800MQ
CPU 2.70GHz.}
\end{center}
\end{table}

\end{document}