HAL Id: in2p3-00623484
http://hal.in2p3.fr/in2p3-00623484
Submitted on 20 May 2020
HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entic 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 diusion de documents
scientiques 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.
Three-dimensional Lagrangian Turbulent Diusion of
Dust Grains in a Protoplanetary Disk: Method and
First Applications
S. Charnoz, L. Fouchet, J. Aléon, Manuel Moreira
To cite this version:
S. Charnoz, L. Fouchet, J. Aléon, Manuel Moreira. Three-dimensional Lagrangian Turbulent Diusion
of Dust Grains in a Protoplanetary Disk: Method and First Applications. The Astrophysical Journal,
American Astronomical Society, 2011, 737, pp.33. �10.1088/0004-637X/737/1/33�. �in2p3-00623484�
The Astrophysical Journal, 737:33 (17pp), 2011 August 10 doi:10.1088/0004-637X/737/1/33
C
2011. The American Astronomical Society. All rights reserved. Printed in the U.S.A.
THREE-DIMENSIONAL LAGRANGIAN TURBULENT DIFFUSION OF DUST GRAINS IN A
PROTOPLANETARY DISK: METHOD AND FIRST APPLICATIONS
S
´
ebastien Charnoz
1,5
, Laure Fouchet
2
,J
´
er
ˆ
ome Aleon
3
, and Manuel Moreira
4
1
Laboratoire AIM, Universit
´
e Paris Diderot/CEA/CNRS, UMR 7158, 91191 Gif sur Yvette cedex, France; charnoz@cea.fr
2
Physikalisches Institut, Universit
¨
at Bern, CH-3012 Bern, Switzerland
3
Centre de Spectrom
´
etrie Nucl
´
eaire et de Spectrom
´
etrie de Masse, CNRS/IN2P3, Universit
´
e Paris Sud 11, B
ˆ
atiment 104, 91405 Orsay Campus, France
4
Institut de Physique du Globe de Paris, Universit
´
e Paris-Diderot, UMR CNRS 7154, 1 rue Jussieu, 75238 Paris cedex 05, France
Received 2010 December 20; accepted 2011 May 17; published 2011 July 26
ABSTRACT
In order to understand how the chemical and isotopic compositions of dust grains in a gaseous turbulent
protoplanetary disk are altered during their journey in the disk, it is important to determine their individual
trajectories. We study here the dust-diffusive transport using Lagrangian numerical simulations using the popular
“turbulent diffusion” formalism. However, it is naturally expressed in an Eulerian form, which does not allow the
trajectories of individual particles to be studied. We present a simple stochastic and physically justified procedure
for modeling turbulent diffusion in a Lagrangian form that overcomes these difficulties. We show that a net diffusive
flux F of the dust appears and that it is proportional to the gas density (ρ) gradient and the dust diffusion coefficient
D
d
: (F = D
d
/ρ × grad(ρ)). It induces an inward transport of dust in the disk’s midplane, while favoring outward
transport in the disk’s upper layers. We present tests and applications comparing dust diffusion in the midplane and
upper layers as well as sample trajectories of particles with different sizes. We also discuss potential applications
for cosmochemistry and smoothed particle hydrodynamic codes.
Key words: methods: numerical – protoplanetary disks
Online-only material: color figures
1. INTRODUCTION
The transport of solids in turbulent, gaseous protoplanetary
disks involves highly diverse physical processes such as gas
drag, turbulence, photophoretic gas pressure and radiation pres-
sure. Modeling these processes is important and necessary as
an increasing amount of observational data testifies that, both
in our own solar system and in protoplanetary disks orbit-
ing distant stars, large-scale transport of solids occurs over
tenths of astronomical units (AU). Samples from the comet
81P/Wild2 (Brownlee et al. 2006; Zolensky et al. 2006;
Westphal et al. 2009) and observations of comet 9P/Tempel
1 (Lisse et al. 2006) have revealed the presence of crystalline
silicates, which may have originated within 2 AU of the Sun and
then been transported outward by some mechanism. In this re-
spect, perhaps the most spectacular result of the Stardust mission
to comet Wild2 is the discovery of refractory inclusions formed
at high temperatures (1500 K) in the innermost regions of
the solar protoplanetary disk (1 AU) in a comet originating in
the Kuiper Belt (e.g., Zolensky et al. 2006; Simon et al. 2008;
Matzel et al. 2010). Similar observations made with the Spitzer
Space Telescope have revealed that crystalline silicates are also
found in T Tauri disks with crystalline-to-amorphous silicate
ratios varying greatly from one disk to another (Van Boekel
et al. 2005; Olofsson et al. 2009;Watsonetal.2009). These
results suggest that global dust transport processes are indeed
active, but that their nature and efficiency may differ signifi-
cantly from one disk to another. Since dust settling is expected
to modify the photometric appearance of a disk (Dullemond &
Dominik 2004), observations can reveal the disk structure and
signs of transport. For example, in the GG Tauri circumbinary
disk (Duch
ˆ
ene et al. 2004; Pinte et al. 2007), multi-wavelength
5
Author to whom any correspondence should be addressed.
observations have revealed a radial size sorting of dust particles
at least qualitatively consistent with the radial migration induced
by gas drag.
Many studies have addressed different physical aspects of
dust transport (for radiation pressure, see, e.g., Vinkovi
´
c 2009;
photophoresis, see, e.g., Krauss & Wurm 2005; stellar wind,
see, e.g., Shu et al. 2001; turbulent diffusion, see, e.g., Gail
2001;Ciesla2009; photoevaporation, see, e.g., Alexander &
Armitage 2009). These studies use various tools and formalisms,
either Eulerian or Lagrangian, which are most often designed to
study one specific aspect of transport. For example, multi-fluid
simulations are appropriate for describing the turbulent transport
of fine dust (see, e.g., Fromang & Nelson 2009), but they are very
computationally demanding. At the opposite end, populations
of large particles (which are less collisional) are not accurately
described by a fluid approach and Lagrangian approaches seem
to be more adapted given that they track individual trajectories
(see, e.g., Johansen & Youdin 2007; Johansen et al. 2007).
It would thus be useful to build a three-dimensional (3D)
modular code that allows for a simple coupling of the different
physical processes of dust transport. As a first step, the present
paper describes a 3D Lagrangian dust transport code in which
the equation of motion is numerically solved with as few
approximations as possible. While the motion of small dust
particles that are tightly coupled to gas is analytically well
known, the motion of particles loosely coupled to gas requires
direct integration. Here we focus on the motion of dust particles
under gas drag and turbulent diffusion. Other processes will
be incorporated into future work. Because the effect of gas
drag in a laminar disk using a Lagrangian description is largely
documented (see, e.g., Barri
`
ere-Fouchet et al. 2005), this paper
focuses mainly on the most difficult aspect: the inclusion of
turbulent diffusion in a Lagrangian code. We describe below the
various features of this code.
1
The Astrophysical Journal, 737:33 (17pp), 2011 August 10 Charnoz et al.
1.1. A Lagrangian Approach
A Lagrangian code treats each particle individually, which
means that it is not well suited for representing a fluid system
where the mean free path is short. It is, however, well designed
for tracking pointlike particles in a gas disk that have very few
pairwise interactions. One very attractive advantage of this ap-
proach is that it affords the possibility of tracking individual
particle trajectories, and thus, of reconstructing the thermody-
namical history of particles in a protoplanetary disk environ-
ment. This should have many applications for cosmochemistry.
For instance, chondrules and refractory inclusions, which are
the main millimeter- to-centimeter-sized components of prim-
itive chondritic meteorites, are thought to have undergone a
complex multi-stage evolution during the first 2–3 Myr of so-
lar system history, before being incorporated into asteroidal or
cometary planetesimals at various heliocentric distances. Chon-
drules (Mg–Fe-rich) and refractory inclusions (also known as
calcium–aluminum inclusions, CAIs hereafter) were formed in
the inner solar system and experienced multiple heating and ir-
radiation events before, or during, their transport to the region
where they were finally incorporated into a meteorite. As they
are the major building blocks of rocky planetesimals, track-
ing their thermodynamical history provides key information for
deciphering the physics and chemistry of planetary formation.
Similarly, the history of frozen volatile-rich dust grains from
outer solar system regions is essential to understanding the so-
lar protoplanetary disk and to unraveling the origin of planetary
volatiles, such as water or organic molecules with prebiotic po-
tential.
1.2. Turbulent Diffusion
The disk is expected to be magnetorotational instability (MRI)
turbulent, inducing an efficient mixing of the dust component
and the gas component. A direct approach would be to couple
a particle-based dust transport code with a 3D magnetohydro-
dynamic (MHD) simulation of a turbulent gas disk (as in Fro-
mang & Nelson 2005; Johansen & Youdin 2007; Johansen et al.
2007). However, while conceptually simple, MHD simulations
remain very computationally demanding, and consequently are
currently limited to a few thousand orbits at most. For this rea-
son, we turn our attention to a simplified turbulence model,
the popular turbulent diffusion model, which mimics turbulent
transport as a diffusive process through a Brownian motion with
an efficiency parameter, the dust diffusion coefficient D
d
.This
diffusion coefficient is thought to be comparable in magnitude
to the turbulent viscosity coefficient. More recently, Fromang
&Nelson(2009) have shown that D
d
may increase with the
distance from the midplane. The description of turbulence as a
diffusive process, though not fully accurate, is widely used and
underpinned by vast literature. We have therefore built on earlier
studies to develop an efficient Lagrangian code. For example,
Lagrangian diffusion has been extensively used in environmen-
tal studies for the transport of air and oceanborne pollutants
(see Wilson & Sawford 1996 for a review). In the planetary
science literature, several models couple gas drag and turbu-
lent diffusion within a Lagrangian framework, but this is either
in a form not adapted to large particles (which decouple from
the gas and undergo oscillations) or limited to some specific
prescription of the gas density field. For example, in Ciesla
(2010) and in Hughes & Armitage (2010), a one-dimensional
(1D) stochastic diffusion model (partly similar to ours; see
Section 2) is applied, but the treatment of gas density variations
is different. Our method is readily applicable to any 3D gas disk
and thus much more general than that used in Ciesla (2010). It is
also more physically justified than in Hughes & Armitage’s
(2010) model, which is mainly empirical. Another major differ-
ence is our use of an implicit solver to integrate the dust motion,
which allows us to accurately follow any range of particle size,
from the most tightly gas-coupled particles, such as polycyclic
aromatic hydrocarbons (PAHs), to those totally decoupled from
the gas, such as kilometer-sized planetesimals (see the exam-
ples in Section 4.2). We will see, in particular, that proper im-
plementation of turbulent diffusion is not straightforward, with
one frequently encountered problem being that of satisfying the
“good mixing condition,” i.e., reaching an asymptotic state in
which dust is well mixed with the gas for any gas spatial density
distribution, in accordance with the Second Law of Thermo-
dynamics. This requires a special treatment of diffusion in a
Lagrangian approach and constitutes the core of this paper. The
case of a non-constant diffusion is also addressed.
1.3. Three-dimensional System
Another useful physical aspect of our code is that we consider
dust motion in 3D. Transport of dust to altitudes high above
the midplane is expected due to the efficient vertical diffusion
induced by turbulence. Most studies treat radial mixing in the
protoplanetary disk through a one dimensional (1D) approach
(see, e.g., Gail 2001; Brauer et al. 2008; Hughes & Armitage
2010). Yet, two-dimensional (2D) models that explicitly treat
the vertical motion of dust particles (see, e.g., Takeuchi & Lin
2002; Dullemond & Dominik 2004;Ciesla2007; Tscharnuter
&Gail2007) show that the disk is stratified, which may have an
impact on global transport. For example, Takeuchi & Lin (2002)
show that the radial velocity of dust depends quadratically on
Z inducing an outward gas drag in the disk’s upper layers (for
above ∼1.5 pressure scale heights, see Takeuchi & Lin 2002,
Equations (11) and (17)). This depends on the gas density profile
only, and is thus a robust result. We shall also see that radial
diffusion is more effective at altitudes far from the midplane,
due to the shallower slope of the radial density gradient for
Z > 0 (see Section 3.2). The importance of including the
vertical dimension is also emphasized by Bai & Stone (2010)in
the context of dust transport in a dead zone. Thus, in order
to ensure that the approach remains as general as possible,
it is important to incorporate the vertical dimension into the
system and integrate the motion of each particle with the fewest
approximations possible. In the present paper, the azimuthal
direction is also explicitly included. However, as there is no
planet for the moment, the disk remains azimuthally symmetric,
while the system is intrinsically evolved in 3D.
For practical use, this code will be referred to as LIDT3D (for
Lagrangian implicit dust transport in 3D).
The paper is organized as follows: in Section 2, we describe
the procedure used to introduce turbulent diffusion into a
Lagrangian form. It should be noted that the gas disk considered
in this paper is a simple and non-evolving parameterized gaseous
disk (as in Takeuchi & Lin 2002) in order to test the dust
transport algorithm. This algorithm, however, is independent
of the choice of disk and can be readily extended to any disk
sampled on a 3D grid. In Section 3, we present several tests
aimed at reproducing known results about turbulent diffusion in
a gaseous disk, and in Section 4 we discuss some applications
to compare dust diffusion in 2D (at the disk midplane only) and
in 3D and to show individual trajectories of different-sized dust
grains extracted from the simulations.
2
The Astrophysical Journal, 737:33 (17pp), 2011 August 10 Charnoz et al.
2. NUMERICAL IMPLEMENTATION OF
TURBULENT DIFFUSION
2.1. Basic Concepts and Definitions
In a Lagrangian code, the motion of an individual dust particle
is described using the classical Newtonian formalism:
d v
dt
=
−→
F
∗
m
−
v −v
g
τ
, (1)
where F
∗
is the gravitational force of the central star, the second
term is the gas drag force, v is the particle’s velocity, v
g
is the
gas velocity, and m is the particle’s mass. The dust stopping time
τ is in the Epstein regime:
τ =
aρ
s
ρC
s
, (2)
where a is the dust grain radius, ρ
s
is the dust material density,
ρ is the gas density, and C
s
is the local sound velocity. Particles
with sizes larger than the mean free path of gas molecules may
be in the “Stokes regime,” where the exact expression for τ
depends on the gas Reynolds number (see Equation (10) of
Birnstiel et al. 2010). However, our discussion and the method
described here do not closely depend on the expression for the
stopping time τ. We introduce the Stokes number St, which
is the particle coupling time τ divided by the eddy turnover
time τ
ed
(St = τ/τ
ed
). τ
ed
is about 1/Ω
k
(see, e.g., Fromang
& Papaloizou 2006), with Ω
k
representing the local Keplerian
frequency (Ω
k
= (GM
∗
/r
3
)
1/2
), where M
∗
, G, and r are the star’s
mass, the universal gravitational constant, and the distance to the
star projected along the disk midplane, respectively, such that St
∼ τ Ω
k.
In the laminar case, once the gas velocity field is known
by applying, for example, a numerical method like that used by
Tscharnuter & Gail (2007), or an analytical model like that in
Takeuchi & Lin (2002; as is used here), Equation (1) is easy to
solve numerically in the absence of turbulence. However, since
gas is turbulent, the gas velocity field is highly complex, then
it is more convenient to introduce a turbulent diffusion model
into this equation to take the stochastic motion of the gas into
account (see e.g., Hinze 1959; Youdin & Lithwick 2007). We
first transform Equation (1) by decomposing v
g
into a mean-
field contribution v
g
plus a turbulent fluctuation δv
g
, such that
v
g
=v
g
+ δv
g,
so that Equation (1)isrewritten:
d v
dt
=
−→
F
∗
m
−
v −v
g
τ
+
δv
g
τ
, (3)
where δv
g
/τ is the acceleration of dust induced by the turbulent
gas velocity field that induces a random walk (or turbulent dif-
fusion) of the particle, consistent with the “turbulent diffusion”
model. We assume that v
g
can be identified with the unper-
turbed laminar flow. This gas flow can either be numerically
computed as in Ciesla (2009) or Tscharnuter & Gail (2007)
or taken from analytical models as in Takeuchi & Lin (2002).
Whereas the former method is more versatile, we choose to
adopt the latter in the present paper for the sake of simplicity
and also to enable us to test our dust transport model against
analytical results. Yet, the method for including turbulent diffu-
sion that we present below is not dependent on this choice and a
gas disk sampled on a 2D or 3D grid can be easily implemented
to provide the v
g
field.
We introduce δv
T
and δr
T
(kicks in velocity and position)
defined as δr
T
= dt × δv
Tv
, where dt is the time step. They
are decomposed into the following elements by integrating
Equation (3):
δv = δ v
M
+ δ v
T
(4)
δr = δr
M
+ δr
T
, (5)
where δv
M
and δv
T
are velocity increments due to the central
star gravity and gas drag with the disk mean velocity field (δv
M
)
and gas drag with the turbulent velocity field (δv
T
). As regards
positions, δr
M
is position increments due to initial velocity plus
the star’s gravity field and gas drag with the mean flow. δr
T
corresponds to the position increment due to turbulent diffusion
that induces a particle’s random walk. The turbulent random
walk of dust is entirely contained in the δr
T
and δv
T
terms, as
discussed in Section 2.2.Thetermδr
M
is obtained by direct
numerical integration of the equation:
−→
δr
M
=
t+dt
t
v(t) dt (6)
with
d v
dt
=
−→
F
∗
m
−
v −v
g
τ
. (7)
Equations (6) and (7) can be solved using a variety of implicit
or explicit variable-order methods (see, e.g., Press et al. 1992).
An implicit method is highly recommended here given that
Equation (1) is a stiff equation due to the two very different
timescales at play: the stopping timescale τ and the Keplerian
orbital timescale T
k
. The smallest particles have a very short
value of τ (meaning that they are tightly coupled to the gas)
that may be much smaller than T
k
: for example, 0.1 μm and
1 mm particles at 1 AU in the minimum-mass solar nebula have
τ/T
k
about 10
−7
and 10
−3
, respectively. If an explicit integrator
such as the popular explicit fourth-order Runge–Kutta scheme
is used, the integration time step will be limited to a fraction
of the gas-coupling timescale τ . It will then be impossible to
include the smallest particles (with a small Stokes number) or
to integrate them over long timescales.
In the present work. we use a Bulirsch–Stoer scheme with a
semi-implicit solver (Bader & Deuflhard 1983) as described in
Press et al. (1992, chapter 16.6). This yields excellent results in
terms of accuracy and rapidity, at least down to a Stokes number
of about 10
−8
when we use the Bulirsch–Stoer extrapolation
method up to the eighth order in the case of a laminar flow.
However, the taking into account of turbulent diffusion (see the
following section) is only first-order accurate due to operator
split. In this case, the Bulirsch–Stoer scheme is used with a
second-order integrator. Due to the adaptive time step scheme,
it yields excellent results in terms of stability whereas the
intrinsic accuracy is only first order. However, extensive tests
(see Section 4) have shown that our results precisely match
analytical expectations even when turbulence is included.
2.2. The “Good Mixing” Problem
A diffusion process is usually described by Fick’s law, which
relates the diffusive flux of material F
T
to the density gradient:
−→
F
T
=−D ·
−−→
grad(ρ), (8)
where D is the diffusion coefficient and ρ is the local material
density. For the specific case of dust in the solar nebula, D will
be written hereafter as D
d
and ρ is written as ρ
d
.Tomimica
3
The Astrophysical Journal, 737:33 (17pp), 2011 August 10 Charnoz et al.
random walk of the particles and obtain a diffusion flux obeying
Equation (8), a well-known method is to express δr
T
or δv
T
as Gaussian random variables with 0 mean and a variance
dependant on the diffusion coefficient and the time step (see,
e.g., Wilson & Sawford 1996; Hughes & Armitage 2010; and
Appendix A of this paper). We now go on to consider a 1D
variable only, but the procedure is readily generalized to 3D.
We express the variance of δr
T
as σ
2
r
: the variance of δv
T
as σ
2
v
,
δr
T
and δr
T
being their respective mean. In the “position”
representation, the kick on position is a random Gaussian with
δr
T
=
δr
T
=0
σ
2
r
= 2Ddt
. (9)
In the “velocity” representation, the kick on velocity is a random
Gaussian variable δv
T
constructed so that, after time integration,
it induces the same average mean dispersion on positions as δr
T
(i.e., (σ
2
r
)
1/2
= dt × (σ
2
V
)
1/2
):
δv
T
=
⎧
⎨
⎩
δv
T
=0
σ
2
v
=
2D
dt
. (10)
Both methods are possible: one of the two representations must
be chosen.
These kinds of procedures are known to yield good results
for the random walk of a solute particle (i.e., dust) with
constant D
d
inside a solvent (i.e., gas) of uniform density. It
has been used several times in environmental studies, notably to
study the diffusion of pollutants in the atmosphere or ocean
(see, e.g., Wilson & Sawford 1996). For transport of dust
in the protoplanetary disk, Youdin & Lithwick (2007) and
Hughes & Armitage (2010) use similar procedures in which
δv
T
is a discrete random variable allowed to take two values,
+(2D/dt)
1/2
or −(2D/dt)
1/2
.
However, problems arise when the solvent (i.e., gas) has a
non-uniform density in space. In this case, the simple procedure
described above is no longer valid. It can be easily verified that
this method does not satisfy the “good mixing condition”: at
steady state, a solute (i.e., dust) has a uniform concentration
throughout the solvent (i.e., gas) in order to reach maximum
entropy. In other words, it must have the same spatial density as
the solvent times a constant multiplicative factor. The procedure
described above (Equation (9) or Equation (10)) will lead
inexorably to a uniform density distribution of the solute
(i.e., dust) throughout the whole space, even though the solvent
(i.e., gas) has a non-uniform density. For our astrophysical
case, this means that dust would diffuse everywhere in space,
out of the gaseous disk itself, in the absence of central star
gravity. This problem has long been identified in environmental
studies and several solutions exist with different derivations
(see, e.g., Sothl & Thomson 1999 or Ermak & Nasstrom 2000).
Yet most of them are either expressed in terms of gas velocity
fluctuations (which are not known or not “well documented”
in the protoplanetary disk literature) or deal only with the self-
diffusion of non-inertial material rather than being expressed
in terms of diffusion coefficient, as is required here. For the
case of protoplanetary disks, we wish to use only the diffusion
coefficient and gas macroscopic properties in order to compare
our results with previous analytical studies. To cure this problem
of “good mixing,” it could be tempting to introduce a spatial
dependence in the diffusion coefficient D to enforce a uniform
concentration at steady state. But this would run counter to
the usual evaluations of the dust diffusion coefficient D
d
in a
turbulent disk: in an α-turbulent isothermal disk, D
d
is assumed
to be close to the turbulent viscosity D
d
∼ αC
s
H. Since C
s
and H
depend only on r in an isothermal disk, D
d
is a constant function
of Z. Hughes & Armitage (2010) identified this problem of dust-
to-gas ratio and propose an empirical method in which δv
T
is
a non-Gaussian discrete random variable. It is designed so that
there is higher probability of a dust particle migrating toward
higher density regions with weights depending on the local
density profile. However, although this procedure is successful
in the framework of their study, it is performed in the radial
direction only and its extension to 3D systems is somewhat
unclear since the number of possible directions becomes infinite.
Ciesla (2010) presents a physically motivated derivation of
the dust displacement algorithm, but treats the problem in the
vertical direction only, for a disk density in the form ρ(z) =
ρ
0
exp(−z
2
/2H
2
), and within the limit of very small particles
that are strongly coupled to the gas and thus follow the same
path as the gas molecules in the absence of turbulence. We
present below a heuristic derivation to properly account for gas
density variations. While sharing some similarities with Ciesla
(2010), this derivation is more general and not constrained by
all the previously mentioned limitations or assumptions.
2.3. Good Diffusion with a Constant Diffusion Coefficient
For the sake of clarity, we restrict ourselves to the case of a
constant diffusion coefficient in the current section. The more
general case of Lagrangian diffusion with a varying diffusion
coefficient is treated in Section 2.4. We first recall some basic
principles of a 1D Gaussian random walk of a solute particle in a
uniform and static solvent. We call X the position of a Lagrangian
particle. For a particle with Lagrangian velocity V and constant
diffusion coefficient D, the spatial position increment dX due to
random walk as described in Equation (9)is
dX = Vdt +(2Ddt)
1/2
W, (11)
where W is a Gaussian random variable with mean 0 and
variance σ
2
= 1 (such that (2Ddt)
1/2
× W is a random
variable with mean 0 and variance 2Ddt, as in Equation (9)). A
basic result of Einstein–Brownian diffusion is that the resulting
motion has an average position X so that dX/dt = V, and
that the standard deviation (X−X)
2
evolves according to
d/dt((X−X)
2
) = 2D. In a Lagrangian simulation involving
a large number of test particles with positions that evolve
according to Equation (11), the local ensemble average is
a solution to the Eulerian advection diffusion equation of a
solute in a solvent with uniform volume density, which reads in
Cartesian coordinates:
∂ρ
∂t
+
∂ρV
∂x
=
∂
∂x
D
∂ρ
d
∂x
, (12)
also written as
∂ρ
∂t
+
∂
∂x
ρV −D
∂ρ
∂x
= 0. (13)
Bearing in mind that the mass conservation equation reads
∂ρ/∂t + ∂F/∂x = 0, where F is the flux, the term ρV in the
parentheses of Equation (13) is the advective flux and −D∂ρ/
∂x is the diffusive flux. Unfortunately, since the solvent (gas)
density is not uniform, the transport equation of the solute (dust)
does not have exactly the same form as Equation (12). Thus, it
4