Probability density function (PDF) models for particle transport in porous media

Mathematical models based on probability density functions (PDF) have been extensively used in hydrology and subsurface flow problems, to describe the uncertainty in porous media properties (e.g., permeability modelled as random field). Recently, closer to the spirit of PDF models for turbulent flows, some approaches have used this statistical viewpoint also in pore-scale transport processes (fully resolved porous media models). When a concentration field is transported, by advection and diffusion, in a heterogeneous media, in fact, spatial PDFs can be defined to characterise local fluctuations and improve or better understand the closures performed by classical upscaling methods. In the study of hydrodynamical dispersion, for example, PDE-based PDF approach can replace expensive and noisy Lagrangian simulations (e.g. trajectories of drift-diffusion stochastic processes). In this work we derive a joint position-velocity Fokker-Planck equation to model the motion of particles undergoing advection and diffusion in in deterministic or stochastic heterogeneous velocity fields. After appropriate closure assumptions, this description can help deriving rigorously stochastic models for the statistics of Lagrangian velocities. This is very important to be able to characterise the dispersion properties and can, for example, inform velocity evolution processes in Continuous Time Random Walk (CTRW) dispersion models. The closure problem that arises when averaging the Fokker Planck equation shows also interesting similarities with the mixing problem and can be used to propose alternative closures for anomalous dispersion.


Introduction
The evolution of solute and particle transport in heterogeneous porous media is determined by the heterogeneity of the medium, the consequent flow heterogeneity and small scale diffusion [33,20,4]. These processes determine the average transport behaviours and also the fluctuation dynamics. For example, small scale velocity fluctuations and mass transfer processes give rise to large scale transport dynamics characterised by hydrodynamic dispersion, but also non-Fickian transport characteristics such as long tails in breakthrough curves and non-linear evolution of solute dispersion [24,22,2,5,32]. The understanding of these behaviours plays a central role in a series of applications across different fields and applications ranging from geothermal energy to packed bead reactors.
Eulerian and Lagrangian velocity statistics have received increased attention over the last years because of their central role for the understanding of pore-scale dispersion and solute transport and their upscaling from the pore to the Darcy scale. The pioneering study of [29] used 3-dimensional particle tracking velocimetry to study hydrodynamic dispersion in porous media. [2] and [22] studied intermittency in purely advective Lagrangian velocity series in 2 and 3-dimensional synthetic porous media and modelled their evolution using a continuous time random walk approach. [34] analysed the scaling behaviours of Eulerian velocity statistics in 3-dimensional digitised rocks and find stretched exponential distributions for the stream-wise velocities. Similar statistics were found for the distribution of stream-wise Eulerian velocity for 2-dimensional synthetic fibrous media [25]. [19] analysed correlations in medium and flow properties as well as distributions of stream-wise velocities for different 3dimensional porous media. These authors find similar behaviours for flow and structural correlation functions. [15] and [28] used particle tracking velocimetry to study Lagrangian velocities, accelerations and dispersion in 3-dimensional bead packs, which are modelled using spatial velocity Markov models. [15] propose a model to relate pore size and velocity distributions. [3] analysed Eulerian velocity distributions in 2-dimensional synthetic porous media and revealed a relation between the pore aperture distribution and the distribution of the Eulerian velocity magnitude, which is explained by Poiseuille flow in individual pores combined with approximately constant pressure drops [33]. [26] analysed Lagrangian velocity time series in 3-dimensional porous media, based on which they propose a stochastic model for the evolution of particle velocities. [1] proposed a velocity model based on mass conservation and independence of the flow rates in neighboring pore throats to explain exponential tails for the stream-wise velocity distribution in simple 2-dimensional porous media. [32] conducted a thorough analysis of Eulerian and Lagrangian flow attributes for a digitised 3-dimensional Berea sandstone sample with the aim of predicting and characterising anomalous dispersion processes at the pore scale [32].
Most of these works focus on Eulerian and Lagrangian flow attributes for purely advective particle motion. The classical works by [20] and [33] con-sidered the stochastic motion of particles in porous media in the presence of diffusion to derive expressions for longitudinal and transverse hydrodynamic dispersion coefficients. [33] provides an approach for the estimation of the Eulerian velocity distribution based in 3-dimensional porous media based on Poiseuille flow in individual pores and the approximation of a constant pressure drop across different pores. It is argued that diffusion homogenises the particle velocities within pores. [5] used a similar approach to analyse the tailing of breakthrough curves in a 3-dimensional synthetic porous medium within a continuous time random approach parametrised by the Eulerian velocity distribution combined with diffusive mass transfer. [30] use a copula based method to investigate advective-diffusive particle motion in 3-dimensional porous media. Most of the works cited above are based on phenomenological understanding, rather than rigorously derived from the microscopic underlying mathematical equations. In fact, the standard upscaling techniques for Partial Differential Equations (PDEs) of porous media flow, such as homogenisation and volume averaging, cannot deal with non-Fickian dispersion and non-equilibrium effects [17]. We focus here on an alternative approach based on Probability Density Functions [35,27,31].
The derivation of macroscopic (averaged) transport models, both for particle positions or velocities, can also be seen as a coarse-graining or model reduction [12] of the high-dimensional Langevin governing equations for particle position and velocity. While this has been extensively studied for Hamiltoniantype systems [14,8], the coarse-graining of general advection-diffusion models has been only recently studied for the overdamped Langevin equation [10,23,16].
In this paper we introduce an equation for the joint distribution of particle position and velocity under advection and diffusion and apply these statistical mechanics approaches to study the evolution of the velocity PDF. In section 2, we derive a system of Langevin equations for the evolution of particle position and velocity and the equivalent Fokker-Planck equation for their joint probability density function (PDF). Section 3 focuses on the evolution of the marginal velocity PDFs and closure models. Section 4 provides some simplified models and exact solutions for the velocity PDF in deterministic linear shear flows.
2 Joint position-velocity PDF equation Scalar and particle (advective and diffusive) transport can be modelled with the Ito (also named overdamped Langevin) Stochastic Differential Equation eq :sd e_ po s eq :sd e_ po s with X(t = 0) = X 0 , where W is a n-dimensional Brownian motion, D a diffusion constant, and u(x) is a space-dependent constant velocity field. Equiv-alently this can be described by the Fokker-Planck equation eq :fp _p os eq :fp _p os with f (x; t = 0) = δ(x − X 0 ), and the probability density function (PDF) eq :fp _d ef eq :fp _d ef and the the operator · indicates the average with respect to the Brownian motion (stochastic average) 1 .
The problem with this formulation is that we loose completely the information on the trajectories and, particularly, on the Lagrangian velocities U = u(X) and their time evolution. One way to restore it is to introduce a joint position-velocity PDF. We can extend the Markovian variable to the joint variable (X, U) and write a system of Ito SDEs, by writing U using Ito's lemma and remembering that we have a steady flow field so that the time derivative disappears. Equivalently one can think of it as a Taylor expansion around the point X eq :du _e xp eq :du _e xp where u = u(X) is the velocity field evaluated in the position X and, and ∇ 2 indicating the Hessian operator. Substituting Eq. 1 in Eq. 4 we obtain an explicit SDE for the Lagrangian position and velocity. Considering only the first order terms and remembering that dW 2 = dt eq :sd e_ ve l eq :sd e_ ve l where G = G(X) = ∇u X is the deformation tensor evaluated at X and I n×n is the identity matrix of order n (spacial dimensions and length of vectors U and X).
The corresponding Fokker-Planck equation [21,11] for the joint PDF f = f (x, v; t) can be therefore written as 2 eq :fp _j oin t eq :fp _j oin t 1 In the following we will consider different averaging operators. To distinguish the molecular Brownian noise from random velocity fields, we will omit, for the random fields, the explicit dependence on a generic random event ω and denote with · the averaging with respect to the random field (either in an ensemble or spatial sense), while the average with respect with the noise will be denoted by · 2 We are dealing with a six-dimensional Fokker-Planck equation whose 2n × 2n diffusion matrix is eq :fp _j oin t2 eq :fp _j oin t2 where we have denoted the Frobenius double inner product with ":" notation, ∆ is the Laplacian and ∇ 2 is the Hessian. It is interesting to notice the presence of a mixed position-velocity diffusion term and of three different velocity-drift terms, two of them depending on diffusivity. The third drift term arises from the mixed position-velocity diffusion term. The other drift and diffusion coefficients can be instead taken out from the derivative operators unchanged, since they do not depend on v (assuming noninertial particles 3 , i.e., ∇ v u = 0 such that ∇ v G = 0) and we have assumed an incompressible flow field (∇ x · u = 0) 4 .
If we assume that u is the solution of a stationary Stokes flow, then the term ∆u = 1 µ ∇p, where p is the pressure.

Velocity marginal PDF and closures
se c:m ar gi na l In certain cases, we can be interested solely on the evolution of the velocity PDF. We can therefore marginalise the general Eq. 6, by defining F = f dx = f and neglecting all the boundary terms arising from spatial derivatives, obtaining eq :fp _j oin t_ ma rg eq :fp _j oin t_ ma rg Now we can decompose f (x, v; t) = f (x|v; t)F (v; t) and, denoting conditional averaging over f (x|v; t) with superscript x|v, the equation can be written as eq :fp _j oin t_ ma rg _c on d_ clo s eq :fp _j oin t_ ma rg _c on d_ clo s In general, compared to the simplified case considered in Sec. 4.2, the flow map that links each position to a velocity might not be invertible, making the conditional expectations not trivial, depending on time too. We can, in fact, rewrite the conditional probability (see Appendix A), used to compute the closure terms, as where c(x; t) is the concentration field and we have highlighted now the intrinsic dependence of the conditional averaging on time. This means that the closure terms have to be evaluated on the domain Ω v = {x | u(x) = v} with the concentration field as a weight, i.e., for a generic function g = g(x), As we will see in Eq. 22, this closure gives an exact analytical formula for simple shear flows. More generally, in Eq. 9, one has to compute these closures as functions of velocity v and time t.
Constant shear: For the simple case of constant shear rate G, this can be explicitly closed and written as an advection diffusion equation eq :fp _j oin t_ ma rg _c on stG eq :fp _j oin t_ ma rg _c on stG

Equilibrium Eulerian velocity PDF
At equilibrium, the concentration c(x, t), for divergence-free velocity fields u, tends to a constant value 5 . Therefore, the conditional expectations over the conditional probability f (x|v) is now only a function of v, and does not depend on time. We identify the spatial averages against this distribution, for a generic function g, as is the equilibrium conditional probability. The stationary version of Eq. 9 needs to be satisfied by the equilibrium (Eulerian) velocity PDF F e (v): eq :F eq eq :F eq 5 Provided appropriate boundary conditions, such as local periodicity.

Perturbation near equilibrium
We can now decompose the velocity PDF near the equilibrium as: Applying Eq. 11 to Eq. 9, we obtain: eq :F pr im eta u eq :F pr im eta u and where we have neglected the higher order terms of they type g(x)

Interaction by Exchange with the Mean
We introduce here two simplifying assumptions that allows to close these terms. First, similarly to the IEM model (Interaction by Exchange with the Mean) in turbulent mixing [31], and to the Mori Projector in Statistical Mechanics [13], by assuming that velocity v and velocity gradient G are spatially distributed according to a joint Gaussian distribution, we can rewrite 6 the conditional expectations as inhomogeneous advection and diffusion terms in the velocity, as follows: eq :ie m eq :ie m eq :ie m2 eq :ie m2 eq :ie m3 eq :ie m3 where we have assumed that G = 0 and u = u − u. All averages are now unconditional averages and can be taken out of the derivatives. However they still depend on the concentration over time that acts as a weighting function in the spatial averaging. This means that the evolution equation for the velocity PDF is not not fully closed since it would require the solution of the concentration 6 using the formula for conditional expectations of joint Gaussian distributions field at all times. To overcome this limitation, one can here consider instead the stationary (equilibrium) distribution function to compute the averages. For incompressible flows and non-inertial non-reactive particles in infinite domain or finite domain with appropriate boundary conditions, this is generally equivalent to a uniform distribution in space. This makes the averages fully computable from the Eulerian velocity field, e.g., u u T is the spatial correlation matrix of the velocity field and GG T , that appears with an opposite sign due to the integration by parts, is a contraction of the full velocity gradient covariance G ⊗ G. The last term u ⊗ GG T , despite looking more complicated, can be equally computed as a simple space integration, given the velocity field and its derivatives.

Eddy dispersion closure
An alternative closure for the first term, is based on the perturbation expansion In this case, the purely advective (first) term in Eq. 8, with an eddy dispersion hypothesis, can be written as is expected to depend on v, and possibly be a tensor. In the case of G being a matrix formed by Gaussian random fields with identical radial covariance, this reduces to D G ef f (v) = |v|D G,0 ef f , with constant D G,0 ef f and Euclidean norm |v|. Since, in most cases (e.g., periodic flow fields) the average shear G is zero, this approximation confirms the intuition that the change in particle velocities purely due to advection homogenises quicker for high velocities while slow velocities are mostly explored by diffusion.

Deterministic shear flows
Consider the following 2D shear-flow dynamics eq :2d sh ea r eq :2d sh ea r where u is a shear flow depending on y only (and possibly in a random manner), v is a generic function (specified later), and W is a two-dimensional Wiener process. In the following we will consider a few special cases of this dynamics. Let us assume D x = 0, i.e., molecular diffusion is acting only on y. The joint position-velocity PDF can be written for f = f (x, y, u, v; t) where the dependence on the velocity field is only through the velocity gradient 7 σ = ∂ u ∂y and, from now on, u represents only the internal velocity coordinate.
Under this assumptions, and a few simplifications can be made on the general Fokker-Planck Eq. 6, leading to the following equation: (18) eq :2d sh ea r4 _p df eq :2d sh ea r4 _p df where σ = ∂ σ ∂y . Since no significant dynamics happens in the y− component of the velocity, v, if we start from an initial condition f 0 = δ(v)g(u, x, y), the term σv ∂ f ∂u and the mixed term can be disregarded.
As done for the general case (see Sec. 3), we now perform the marginalisation to obtain an equation for the velocity PDF F (u, t). In an infinite domain, the boundary terms in x and y disappear. Integrating along the channel height y we obtain eq :2d sh ea r_ ma rg ina l eq :2d sh ea r_ ma rg ina l where M σ (u; t) = σ f (y|u; t) dy and M σ 2 (u; t) = σ 2 f (y|u; t) dy. As shown in Sec. 3, this depends, in general on time, as f (y|u; t) ∝ δ(y − y(u))c(y; t). However, as expected in a stratified flow, the evolution of velocities is due solely to the diffusion. Defining a Péclet number as Pe = u0L D , andt = Pe t, we obtain: eq :2d sh ea r_ ma rg ina l_ ad im eq :2d sh ea r_ ma rg ina l_ ad im

Constant shear (Couette flow)
Couette flow is characterised by σ = σ 0 . The velocity marginal PDF does not need any closure. This reduces to a simple heat equation for F (u).
eq :2d sh ea r_ co ns ta nt _p df eq :2d sh ea r_ co ns ta nt _p df 7 We assume the velocity field is smooth enough, i.e. with a Gaussian correlation. In this way the derivative exists and it is Gaussian itself with Gaussian correlation function.

Linear shear (Hagen-Poiseuille flow)
se c:s he ar This is the case of a two-dimensional channel flow with σ = −2y u0 L 2 and σ = −2 u0 L 2 . In this case, there exists a deterministic relation between y and u, i.e., u(y) = u0 L 2 (L − y)(L + y) = u 0 (1 − y 2 L 2 ), such that y 2 = L 2 1 − u u0 and the solution for long time limits, near equilibrium, can be computed analytically.
Equilibrium closure: With the equilibrium closure, considering a constant concentration c(x, y), M σ and M σ 2 are constant in time. Eq. 19 can be therefore be closed as an inhomogeneous diffusion in the velocity space: eq :2d sh ea r_ lin ea r_ pd f2 eq :2d sh ea r_ lin ea r_ pd f2 The stationary solution of this equation is given by: eq :2d sh ea r_ lin ea r_ lon gt im e eq :2d sh ea r_ lin ea r_ lon gt im e This solution is, in fact, consistent with the fact that, when particles are uniformly spread throughout the channel, their velocity is simply given by the inverse function of the velocity profile. In Fig. 1, we test the equilibrium closure and its convergence to the Eulerian velocity PDF. Eq. 22 is solved with the Matlab library Chebfun [9], for L = 1, u 0 = 1, D = 1, on the domain [0, 1], with no-flux boundary conditions. The initial conditions (crosses) are such that particles (or solutes) are injected only in 20% highest (red dot-dashed, wrt 20% lowest, blue dashed) velocity regions, and the evolution towards the Eulerian velocity distribution (black continuous curve with asterisks) is depicted for different times (t = 0.005; 0.1; 0.2; 0.4; 0.8; 1.6). As it can be seen, both the initial conditions converges relatively fast towards the equilibrium distribution. In this equation only diffusion is the driving force for this relaxation, therefore this is completely defined by the Péclet number.

Discussion and conclusions
In this paper, we have derived a joint position-velocity Probability Density Function (PDF) model for passive particles flowing in a porous medium or a heterogeneous flow field. This is achieved by extending the phase space and writing down the Stochastic Differential Equation (SDE) and the corresponding Fokker-Planck Partial Differential Equation (PDE). While the (Smoluchowki) position-only PDF equation contains classical advection terms u∇ x c, this new formulation requires only first and second-order derivatives of the flow field that often have higher-frequency oscillations and therefore, better mixing/averaging properties. This aspect, possibly leading to more accurate macroscopic models, will be explored in our future works.
Another advantage of this joint position-velocity formulation is that it can be easily extended to full second-order model with velocity increments given by a given drag force: eq :sd e_ ve l2 eq :sd e_ ve l2 where St is the Stokes number, that depends on particle size and density and the slip velocity (difference between particle velocity and local fluid velocity) is U−u(X)). In this alternative model, it is more appropriate to put the random fluctuating term in the momentum, as in the classical Langevin dynamics. This problem is more complex as the particle velocities are no longer divergence free. This case will be studied in future works.
Since the joint PDF model is defined on a high-dimensional space (6+1 dimensions). We have therefore study the reduced equations obtained by integrating (marginalising) over space (assuming therefore periodicity or ergodicity) to focus on the evolution of Lagrangian particle velocities. Connections with the theory of coarse-graining are established and several possible closures are proposed. This approach, first derived for a general flow field, is then applied for two-dimensional shear flows. Here, the effect of the equilibrium, and short-time approximations are analysed.
The study of the evolution of particle Lagrangian velocities is motivated by recent studies that highlight the importance of modelling the evolution of Lagrangian velocities [7,6] for predicting anomalous transport in heterogeneous media. For the first time, in this work, we propose a rigorous approach to derive an evolution equation for the Lagrangian velocity PDF.
In Fig. 2, as a further motivation for further research in this area, we report the evolution of the velocity (magnitude) PDFs in a three-dimensional random pore-scale packing of non-spherical grains studied in [18,6]. The initial conditions (crosses) are such that particles (or solutes) are injected only in 20% highest (red dot-dashed, wrt 20% lowest, blue dashed) velocity regions, and the evolution towards the Eulerian velocity distribution (black continuous curve with asterisks) is depicted for different times and different Péclet numbers. While, in [6], we have demonstrated that the the overall transport in this medium can be well approximated by representing the porous medium as a collection of cylindrical pores, these results cannot be approximated by the theoretical models in Sec. 4.2. In fact, a simple rescaling of time based on the diffusion coefficient (or equivalently the inverse Péclet number) cannot explain the differences shown in Fig. 2. This can be due to the terms we have neglected in Eq. 18, or to the pores heterogeneity that should be represented with an infinite sum of representative pores (see [6]), with different local Péclet numbers. Future works based on the application of theoretical study this type of three-dimensional random porous media will be performed, as well as on the applicability of closures for particle velocities for improved transport models.
eq :ve l_ pd f_ de f eq :ve l_ pd f_ de f where the last step is well defined only locally or when flow field is invertible and x 0 (v) is the inverse of the function u(x) (i.e. the point in space whose velocity is v. This is possibly defined only locally). This highlights the relation between the joint PDF f and the usual spatial concentration PDF c(x, t) = δ(x − X) and the role of the velocity gradient G. In particular, in an homogeneous state, when c(x) = const, the velocity PDF is equivalent to the inverse of the velocity gradient determinant composed with the inverse velocity function