- Email: [email protected]
Asymmetric diffusion in heterogeneous media
Physica A 395 (2014) 193–199
Contents lists available at ScienceDirect
Physica A journal homepage: www.elsevier.com/locate/physa
Asymmetric diffusion in heterogeneous media J. Alvarez-Ramirez ∗ , L. Dagdug, M. Meraz 1 División de Ciencias Básicas e Ingeniería, Universidad Autónoma Metropolitana-Iztapalapa, Apartado Postal 55-534, México D.F., 09340, Mexico
highlights • Diffusion across interfaces between materials with different diffusivity is studied. • Brownian particles run preferentially along the direction of decreasing diffusivity. • An uphill diffusion effect is found when particles move from high to low diffusivity layers.
article
info
Article history: Received 17 May 2013 Received in revised form 26 September 2013 Available online 22 October 2013 Keywords: Diffusion Heterogeneous systems Interface conditions Discontinuous interface density
abstract Diffusion in heterogeneous media (i.e., layers with discontinuous values of the diffusivity parameter) is found in different situations, including geology, catalyst, physiology and separation processes. This work considers an array of rectangular layers of different diffusivity as a simple model for studying diffusivity in heterogeneous media. Random walk simulations adapted for transport across interfaces between different materials are used for showing that heterogeneous media exhibit asymmetric transport in the sense that the diffusive flux depends of the transport direction. It is shown that Brownian particles run preferentially along the direction of decreasing diffusivity. An uphill diffusion effect is found when particles move from high to low diffusivity layers. In terms of Fick’s modeling, asymmetric transport is explained from the presence of a jump in the particle density at layer interfaces. © 2013 Elsevier B.V. All rights reserved.
1. Introduction Diffusion in domains with discontinuous diffusivity parameter (see Fig. 1) appears commonly in many scientific and technological fields, including geology [1–4], ecology [5–7], biomedical images [8], astrophysics [9,10] and molecular dynamics [11]. For instance, in ecology studies landscapes are considered as mosaics of patches of different habitat types with more or less well-defined patch boundaries [12]. In such case, the behavior of individuals in the vicinity of patch boundaries may be the determinant of overall movement behavior [13]. Also, interfaces between individual constituents play an important role in the effective transport properties of composite materials [14]. The modeling the diffusion process for systems involving discontinuous diffusivity requires the determination of matching conditions relating the limiting values on either side of the material interface of the particle density c (x, t ) and its spatial gradient. For the latter, imposing mass conservation across the interface fulfills the spatial gradient requirement. For the former, continuity of particle density c (0− , t ) = c (0+ , t ) is commonly assumed as motivated by intuition and physical experiments of diffusion through a near interface between two different media [2,3,15–18]. However, numerical
∗
Corresponding author. Tel.: +52 55 58044650; fax: +52 55 58044650. E-mail address: [email protected] (J. Alvarez-Ramirez).
1 Current address: Departamento de Biotecnología, Universidad Autónoma Metropolitana-Iztapalapa, Apartado Postal 55-534, 09340 México D.F., Mexico. 0378-4371/$ – see front matter © 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.physa.2013.10.027
194
J. Alvarez-Ramirez et al. / Physica A 395 (2014) 193–199
Fig. 1. Schematic description of a layered diffusion system with different diffusivity parameters.
stochastic simulations and experimental evidence indicated that the continuity of the particle density cannot be applied for heterogeneous systems with different value of the diffusivity parameter [19–25]. Theoretical work on random walk transport across interfaces with discontinuous diffusivity has indicated that particle density should meet a jump condition related to the ratio of diffusivities [26–28]. It was predicted that the probability density may have a point mass at the interface region, meaning that particles accumulate at the interface. Also, depending of the direction of flux, particles can move against density gradients in the interface vicinity, exhibiting a type of uphill diffusion induced by diffusivity discontinuity. A drift was detected of the particles’ individual positions in the direction of the diffusion coefficient gradient, in the absence of any external force or concentration gradient [29–32]. Some numerical results have pointed out that drift velocity and driving force are not directly proportional in the case of inhomogeneous suspensions, where a space dependent mobility induces an additional contribution to the drift velocity [33]. It has also been argued that the classical Fick law is not applicable in the vicinity of sharp inhomogeneities and that the Fokker–Planck equation should be conceded preference [34,35]. The present work considers the diffusion of particles in heterogeneous (layered) systems as described in Fig. 1. In a first step, in order to gain insight into the effects of the discontinuous diffusion parameter, random walk simulations adapted for discontinuous media [36,37] are performed. For D− ̸= D+ , it is found that the resident particle density can be considered as discontinuous at the interface, and the form of such discontinuity depends of the particle direction. In this way, the commonly used matching condition c (0− , t ) = c (0+ , t ) should be discarded from the modeling of the diffusion system. In order to close the problem, an interface jump condition is imposed [26–28], which is shown to be equivalent to the continuity of the chemical potential spatial gradient. This jump condition allows the recovery of the salient features near the interface as observed from the random walk simulations. In particular, in agreement with random walk simulations, asymmetric transport in the sense that the diffusive flux is dependent of the transport direction is shown to be exhibited by heterogeneous media. Also, an uphill diffusion effect is found when particles move from high to low diffusivity layers. The paper is organized as follows. Section 2 describes the numerical algorithm for performing random walk simulations across boundaries with discontinuous diffusivity. Section 3 discusses the numerical results in terms of the particle density jump condition at the system interface. Section 4 closes the work with some conclusions. 2. Random walk across the interface The motion of Brownian particles in a rectangular layered system like that in Fig. 1 is considered. Given that only the effects of the discontinuous interface are studied, the simulation of Brownian random walkers is reduced to only the horizontal direction, so the simulation problem is one-dimensional. In each layer, corresponding to high and low diffusivities, particles undergo Brownian motion governed by the Langevin equations dx dt
√ =
2Dξ (t )
with the diffusivity D given by D = D− , D = D+ ,
−L < x < 0 0 < x < +L.
Here, ξ (t ) is a zero-mean white Gaussian noise with autocorrelation function ⟨ξ (t ), ξ (t ′ )⟩ = δ(t − t ′ ). The physical domain is constrained to the interval x ∈ (−L, L). For each layer, the Langevin equation was numerically integrated by using
J. Alvarez-Ramirez et al. / Physica A 395 (2014) 193–199
195
the Euler algorithm. A particular situation arises when the particle crosses the interface located at x = 0. The simulation problem is that particle can a start at, e.g., the left of the interface (x < 0) where D = D− and end at the right of the interface (x > 0) with D = D+ . The question is, what diffusivity value should be used in√ the Langevin equation for accurate simulation of the physical phenomenon? The integration of the Langevin equation dx = 2Dξ (t ) in the vicinity of the interface follows dt the ideas of skewed diffusion [36]. For a given position x(t ) < 0, the procedure can be described as follows:
√ • Step 1. Compute the predicted position x(t + 1t ) = x(t ) + 2D− 1t ξ (t ). • Step 2. If x(t + 1t ) < 0, repeat Step 1. Otherwise, compute the time fraction that the particle spends moving in the (x(t )/ξ (t ))2 interface left-hand side by solving the equation 0 = x(t ) + 2D− 1tf ξ (t ). This gives 1tf = . 2D− • Step 3. Update the particle position with the remaining time after crossing the interface; namely, x(t + 1t ) = 2D+ (1t − 1tf )ξ (t ).
A similar procedure is implemented for particle position at the right-hand side of the interface, x(t ) > 0. The mean passage time was obtained as ensemble averages over NT = 105 trajectories starting in a domain boundary. Particle trajectory starts at the left boundary x = −L and finishes when the particle achieves the right boundary x = +L. The passage time is taken as the time the particle spends to move from the left boundary to the right boundary. Trajectories starting in the left boundary and returning to it are discarded for avoiding bias in the estimation of the mean passage time. To monitor the transport process, the variable of interest is the flux of particles, J. For a given trajectory, let TFPT be the first passage time (i.e., the time required by a trajectory to pass from the left to the right boundary). If TMFPT = ⟨TFPT ⟩ is the mean value of the first passage time obtained by averaging over NT trajectories, the flux can be expressed by J =
1 TMFPT
.
The flux J is an index of the ease with which the particles are transported along the layered system. Physically, the quantity J indicates the mean number of particles leaving the system domain. The flux J is a quantity of interest in diverse fields because it can be measured and, for instance in chemical reactors and membranes, reflect the ability of the physical system for passive mass transport. The effect of the transport direction (from low to high diffusivity and vice versa) in the resident particle density profile and the particle flux was evaluated by means of numerical simulations. To this end, let us introduce a diffusivity ratio as rD = Dhigh /Dlow , where Dhigh = max{D− , D+ } and Dlow = min{D− , D+ }, so that rD ≥ 1. The numerical simulations described below were performed for Dhigh = 1, L = 1 and 1t = 10−5 . The resident particle density profile is computed as c (x, t ) = N (x, t )/NT , where N (x, t ) is the number of particles at the position x for given time t. Recall that NT = 105 . For two different values of the diffusivity ratio rD and for the time t = 0.5, Fig. 2(a) and (b) present the resident particle density profiles when Brownian particles undergo low-to-high (LH) and high-to-low (HL) diffusivity transport, respectively. In the former case, D− < D+ while in the latter case D+ < D− . To model the physical transport of particles across the inhomogeneous system, all trajectories started in the left (resp., right) boundary for particles undergoing LH (resp., HG) diffusivity transport. The following can be noticed: (a) The resident particle density exhibits a sharp transition at the interface, and this transition can be considered as a discontinuity for space scales of the order of 0.025L. This result is in agreement with previous findings pointing out that the continuity of the particle density cannot be applied for macroscopic modeling of diffusion across interfaces [21–28]. (b) The geometry of the discontinuity depends of the order of the layers. That is, while the LH diffusion transport induces a monotonous discontinuity where the resident particle density decreases abruptly at the interface, the HL diffusion transport yields a non-monotonous discontinuity where the resident particle density increases abruptly at the interface. It is expected that the differences in the geometry of the particle density discontinuity has an important effect in the transport of particles across the system domain. Also for t = 0.5, Fig. 3(a) shows the particle flux J (rD ) normalized with respect to the particle flux for rD = 1. It is noted that the flux depends of the layer ordering. In fact, the HL mass flux is higher than the LH mass flux, indicating an asymmetry in the diffusion transport when interfaces are present. That is, the results in Fig. 3(a) indicate that mass flux runs preferentially along the direction of decreasing diffusivity. Fig. 3(b) displays the mass flux ratio γ (rD ) = JHL (rD )/JLH (rD ), where JLH (rD ) and JHL (rD ) are the LH and the HL diffusivity fluxes, respectively. It is shown that the diffusion transport asymmetry, quantified in terms of γ (rD ), increases with the diffusivity ratio rD , and such asymmetry is a consequence of the interface effects in the diffusion of Brownian particles. 3. Interface matching conditions For individual layers in Fig. 1, the average mass transport can be described by means of the mass balance in combination with the first Fick law. However, given the heterogeneous nature of the transport system, and motivated by the results obtained from random walk simulation, the following question can be posed: What interface (i.e., matching) conditions should be imposed in order to recover the behavior described in Figs. 2 and 3? Undoubtedly, the mass conservation across the interface should be imposed. However, given the particle density profile shown in Fig. 2, it is apparent that continuity of resident particle density cannot be imposed as an interface matching condition. To address this issue, consider the simple
196
J. Alvarez-Ramirez et al. / Physica A 395 (2014) 193–199
Fig. 2. Resident particle density profiles for (a), D− < D+ and (b) D− > D+ . Note the abrupt resident particle density change in a vicinity of the interface.
layered system in Fig. 1. The description of the transport phenomenon along the horizontal direction for individual domains can be made with the standard diffusion equation:
∂ c (x, t ) ∂ 2 c (x, t ) = D− , −L < x < 0 ∂t ∂ x2 (1) 2 ∂ c (x, t ) ∂ c (x, t ) = D+ , 0 < x < +L ∂t ∂ x2 where c (x, t ) is the mass density and the diffusivity parameters satisfy D− ̸= D+ . A complete modeling involves imposing transport matching conditions at the interface located at x = 0. Mass conservation is a mandatory interface matching condition that is established according to the Fick law of the form
− D−
∂ c (0+ , t ) ∂ c (0− , t ) = −D+ . ∂x ∂x
(2) ∂ c (0− ,t )
∂ c (0+ ,t )
It is noted that D− ̸= D+ implies gradient discontinuity at the interface. For instance, > ∂ x if D− < D+ . ∂x Given the parabolic nature of the system Eq. (1), a second interface matching condition is required for closing the model. Mass density continuity c (0− , t ) = c (0+ , t ) is commonly imposed as motivated by intuition and physical experiments of diffusion through a near interface between two different media [2,3,15–18]. However, Figs. 2 and 3 and reported numerical stochastic simulations and experimental evidence indicated that the continuity of the particle density cannot be applied for cases with different value of the diffusivity parameter [19–25]. As a consequence, the particle density should present a neat jump at the interface. Ovaskainen and Cornell [26] provided a rigorous proof that normal diffusion across interfaces should meet the following jump condition: D− c (0− , t ) = D+ c (0+ , t ).
(3)
In this way, since D− ̸= D+ one has that a mass density discontinuity is present at the interface (i.e., c (0− , t ) ̸= c (0+ , t )). Subsequently, Marseguerra and Zoia [27] showed that this jump condition is equivalent to v− c (0− , t ) = v+ c (0+ , t ), where v− and v+ are the velocities of the random walkers at the left- and right-hand side of the interface, respectively. It has been also shown [28] that the interface jump condition can be extended for anomalous diffusion transport. The interface matching ∂ ln c (0− ,t )
∂ ln c (0+ ,t )
= . Physically, the above relationship can be interpreted conditions (2) and (3) can be combined to give ∂x ∂x as the equality of the chemical potential gradient across the discontinuous interface [31,38].
J. Alvarez-Ramirez et al. / Physica A 395 (2014) 193–199
197
Fig. 3. (a) Diffusion flux as a function of the diffusivity ratio rD = Dhigh /Dlow when the transport goes from a low diffusivity layer to a high diffusivity layer, and vice versa. It is noted that the flux exhibits asymmetries depending on the ordering of the diffusion layers. (b) Diffusion flux ratio γ = JLH /JHL as a function of the diffusivity ratio rD .
It will be shown that the interface condition Eq. (3) recovers the transport asymmetry obtained by means of Brownian particle simulations, as already described by Figs. 2 and 3. To this end, for the sake of simplicity, the case of non-transient diffusion transport with fixed boundary particle density is considered. That is, the diffusion transport is governed by 0 = D− 0 = D+
d2 c (x, t ) dx2 2 d c (x, t )
,
−L < x < 0 (4)
,
dx2 with boundary conditions
0 < x < +L
c (−L) = cl
(5)
c (+L) = cr
and interface conditions given by Eqs. (2) and (3). Physically, Eqs. (1)–(3) can be seen as a model describing the diffusion transport across a membrane composed of two different materials. The solution to Eq. (1) is given by c (x) = a− x + b− , c (x) = a+ x + b+ ,
−L < x < 0 0 < x < +L .
(6)
The set of constants {a− , a+ , b− , b+ } is obtained after the application of the boundary and interface conditions. It can be shown that these constants are given by a− = a+ = b− = b+ =
r± cr − cl 2L r± cr − cl 2Lr± r± cr + cl 2 r± cr + cl 2r±
(7)
198
J. Alvarez-Ramirez et al. / Physica A 395 (2014) 193–199
where r± = D+ /D− . Note that r± = 1 implies that a− = a+ and b− = b+ . The following results can be established according to Eqs. (6) and (7):
• The interface mass density is discontinuous. Eq. (3) imposes a jump in the mass density at the interface. In fact, c (0− ) = b− and c (0+ ) = b+ . In addition, c (0− ) = r± c (0+ ). In this way, in agreement with the Brownian particle simulations exhibited in Fig. 2, c (0− ) < c (0+ ) for r± < 1, and c (0− ) > c (0+ ) for r± > 1 regardless of the left and right boundary particle densities cl and cr . This is in agreement with experimental findings that used two liquids having different gelatin concentrations and thus different diffusivity [34]. Such experiments showed the spontaneous appearance of a relatively sharp transition centered on the interface between the two liquids. The numerical results discussed in the above section indicate that the sharp transition is a geometric effect induced by diffusivity discontinuity. On the other hand, the interface matching condition given by Eq. (3) establishes that the mass density discontinuity arises for balancing chemical potential across the interface. D + cr • The mass flux is asymmetric. The flux across the layered system is given by J = D− cl − . Assume that cl > cr , so that the 2L particles are transported from the left boundary (i.e., x = −L) to the right boundary (i.e., x = +L). The mass flux when particles go from the high diffusivity layer to low diffusivity is given by JHL =
Dhigh cl −Dlow cr 2L
. Similarly, the mass flux when
particles are transported from the low diffusivity layer to high diffusivity is given by JLH = mass flux ratio γ = JLH /JHL is given by Dlow cl − Dhigh cr
γ (r D ) =
Dhigh cl − Dlow cr
=
cl − rD cr rD cl − cr
Dlow cl −Dhigh cr 2L
. One has that the
.
(8)
Recall that rD = Dhigh /Dlow . Let β = cl /cr be the boundary particle density ratio. The derivative of the flux ratio is
γ ′ (rD ) =
1 − β2
(β rD − 1)2
.
(9)
As in the Brownian particle simulations, it was assumed that β > 1. Then, one has that γ ′ (rD ) < 1, which is in agreement with the numerical findings shown in Fig. 3(b). In this way, the value of the mass flux depends of how layers are ordered with respect to the flux direction. Diffusion transport is faster for the HL transport direction than for LH transport direction. It is noted that this mass flux asymmetry is not found when continuous interface particle density is assumed. Analogously to electric rectifiers in electronic circuits and to the recently considered problem of heat transport rectifiers [39,40], the fact that mass is transported preferentially in one direction to the other suggests that layered material with discontinuous diffusivity parameter can be used as mass transport rectifiers. Here, particles move easily in one direction while finding high resistance in the reverse direction. Potential applications of mass transport rectification can be for partially impregnated catalyst where the catalyst pellet is composed of a mesoporous inert support containing islands or eggshells of active microporous support [41,42]. • Interface resistance can be negative. The mass flux across the layered system can be expressed as J =
1ceff
(10)
Reff
where 1ceff is an effective particle density gradient and Reff is the corresponding effective mass transport resistance. These functions are given by
1ceff = cl − r± cr
(11)
Reff = R− + Ri + R+
(12)
and
where the standard resistance for individual layers is given as R− = L/D− and R+ = L/D+ . Analogous to Kapitza interface resistance for heat transfer [43], Ri (r± ) represents the interface resistance for mass transfer and is given by
Ri = L
D+ − D− D+ D−
.
(13)
In this regard, the following can be pointed out: (a) With respect to the physical particle density gradient 1c = cl − cr , the effective particle density gradient 1ceff = cl − r± cr is reduced (resp., magnified) for r± > 1 (resp., r± < 1). Also, for non-layered (i.e., homogeneous) systems, the physical particle density gradient 1ceff (1) = cl − cr is recovered. (b) Eq. (13) reflects the particle resistance to jumping from one layer to the other. As expected, one has that Ri = 0 when D+ = D− . Interestingly, negative (resp., positive) interface resistance is obtained for D+ < D− (resp., D+ > D− ). Since cl > cr was assumed, negative (resp., positive) interface resistance is obtained for the HL (resp., LH) transport direction. In this way, increased mass flux in the HL transport direction is obtained by the combined effect of magnified particle density gradient and negative interface resistance. It should be pointed out that negative transport resistances have been also observed in heat transport through shaped graded materials [44]. Analogously to the results exhibited in Fig. 2(b),
J. Alvarez-Ramirez et al. / Physica A 395 (2014) 193–199
199
thermal energy is apparently transferred from a region of lower temperature to a region of higher temperature across composed materials of graded conductivity. Finally, an apparent paradoxical situation is found when particles move in the HL direction. Fig. 2(b) shows that the resident density exhibits a sharp increase in a vicinity of the interface. This indicates that particles diffuse up against density gradient (i.e., uphill diffusion transport). The interface is a geometric anomaly where the conditions for equilibrium transport are not met. The sharp increase of the resident particle density in the interface vicinity appears because particles are quickly transported in the high-diffusivity layer. In this way, particles are accumulated in the interface vicinity where particle transport faces a higher resistance. Uphill diffusion induced by discontinuous interfaces has been already observed in experiments for, e.g., boron transport in ultrashallow junctions [45,46] and carbon diffusion in steel layers [47]. It should ∂ ln c (0− ,t )
∂ ln c (0+ ,t )
= at the interface, particles are still transported along negative gradients of be remarked that, since ∂x ∂x chemical potential (i.e., partial molar free energy), preserving the thermodynamic constraints for diffusion processes [47,48]. 4. Conclusions Numerical simulations of Brownian trajectories in heterogeneous media (i.e., layers with discontinuous values of the diffusivity parameter) indicated that the commonly assumed condition of continuous interface particle density for diffusion transport cannot be applied since the particle density in a vicinity of the interface exhibits sharp positive and negative changes. Analogously to electric rectifiers, the mass flux is not symmetric in the sense that particles flow more easily in one direction than in the reverse direction. In fact, particle transport is faster when particles move from high diffusivity to low diffusivity layers than in the reverse direction. By imposing a jump condition for particle density at the interface, simple analytical results showed that particles are transported along chemical potential gradients, and this jump condition is able to explain the asymmetric mass flux of layered systems. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38] [39] [40] [41] [42] [43] [44] [45] [46] [47] [48]
F. Delay, Ph. Ackerer, C. Danquigny, Vadose Zone J. 4 (2005) 360. H. Hoteit, R. Mose, A. Younes, F. Lehmann, Ph. Ackerer, Math. Geol. 34 (2002) 435. E.M. LaBolle, J. Quastel, G.E. Fogg, J. Gravner, Water Resour. Res. 36 (2000) 651. H. Zhan, Zh. Wen, G. Huang, D. Sun, J. Contam. Hydrol. 107 (2009) 162. R.S. Cantrell, C. Cosner, Theor. Popul. Biol. 55 (1999) 189. W.F. Fagan, R.S. Cantrell, Ch. Cosner, Amer. Nat. 153 (1999) 165. J.T. Cronin, Ecology 84 (2003) 1506. E. Fieremans, D.S. Novikov, J.H. Jensen, J.A. Helpern, NMR Biomed. 23 (2010) 711. M. Zhang, Astrophys. J. 541 (2000) 428. A. Marcowith, F. Casse, Astronom. Astrophys. 515 (2010) A90. M. Mascagni, N.A. Simonov, SIAM J. Sci. Comput. 26 (2004) 339. ¯ A. Okubo, Diffusion and Ecological Problems: Mathematical Models, Springer Science + Business Media, 1980. H. Kaiser, Oecologia 56 (1983) 249. K.K. Chawla, Composite Materials: Science and Engineering, Springer, 2009. H.S. Carslaw, J.C. Jaeger, Conduction of Heat in Solids, Clarendon, Oxford, 1959. E.M. LaBolle, G.E. Fogg, A.F.B. Tompson, Water Resour. Res. 32 (1996) 583. E.M. LaBolle, J. Quastel, G.E. Fogg, Water Resour. Res. 34 (1998) 1685. J.Y. Parlange, F. Stagnitti, J.L. Starr, R.D. Braddock, J. Hydrol. 70 (1984) 251. O. Ovaskainen, S.J. Cornell, Proc. Natl. Acad. Sci. 103 (2006) 12781. M.Th. Van Genuchten, J.C. Parker, Soil Sci. Soc. Am. J. 48 (1984) 703. R.C. Schwartz, K.J. McInnes, A.S.R. Juo, L.P. Wilding, D.L. Reddel, Water Resour. Res. 35 (1999) 671. K.S. Novakowski, Water Resour. Res. 28 (1992) 2399. K.S. Novakowski, Water Resour. Res. 28 (1992) 2423. F.J. Leij, J.H. Dane, M.Th. van Genuchten, Soil Sci. Soc. Am. J. 55 (1991) 124. J.-L. Barrat, F. Chiaruttini, Mol. Phys. 101 (2003) 1605. O. Ovaskainen, S.J. Cornell, J. Appl. Probab. 40 (2003) 557. M. Marseguerra, A. Zoia, Ann. Nucl. Energy 33 (2006) 1396. N. Korabel, E. Barkai, Phys. Rev. E 83 (2011) 051113. P. Lançon, G. Batrouni, L. Lobry, N. Ostrowsky, Europhys. Lett. EPL 54 (2001) 28. P. Lançon, G. Batrouni, L. Lobry, N. Ostrowsky, Physica A 304 (2002) 65. M. Smerlak, A. Youssef, ArXiv Preprint. arXiv:1206.3441. F. Marchesoni, Materials 6 (2013) 3598. M. Yang, M. Ripoll, J. Chem. Phys. 136 (2012) 204508. B.P. Van Milligen, P.D. Bons, B.A. Carreras, R. Sánchez, Eur. J. Phys. 26 (2005) 913. F. Sattin, Phys. Lett. A 372 (2008) 3941. A. Lejay, G. Pichot, J. Comput. Phys. 231 (2012) 7299. L.E. Reichl, A Modern Course in Statistical Physics, Wiley & Sons, New York, 1997. S.R. De Groot, P. Mazur, Non-Equilibrium Thermodynamics, Dover Publications, 1984. W. Kobayashi, Y. Teraoka, I. Terasaki, Appl. Phys. Lett. 95 (2009) 171905. S. Narayana, Y. Sato, Phys. Rev. Lett. 108 (2012) 214303. S. Minhas, J.J. Carberry, J. Catal. 14 (1969) 270. A. Graviilidis, A. Varma, M. Morbidelli, Catal. Rev. Sci. Eng. 35 (1993) 399. E.T. Swartz, R.O. Pohl, Rev. Modern Phys. 61 (1989) 605. C.Z. Fan, Y. Gao, J.P. Huang, Appl. Phys. Lett. 92 (2008) 251907. H.C.N. Wang, C.C. Wang, C.S. Chang, T. Wang, P.B. Griffin, C.H. Diaz, IEEE Electron Device Lett. 22 (2001) 65. M. Ferri, S. Solmi, D. Giubertoni, M. Bersani, J.J. Hamilton, M. Kah, N.E.B. Cowern, J. Appl. Phys. 102 (2007) 103707. L.S. Darken, Trans. AIME 180 (1949) 430. D. Kondepudi, Introduction to Modern Thermodynamics, John Wiley & Sons, Ltd., Chichester, UK, 1998.
Contents lists available at ScienceDirect
Physica A journal homepage: www.elsevier.com/locate/physa
Asymmetric diffusion in heterogeneous media J. Alvarez-Ramirez ∗ , L. Dagdug, M. Meraz 1 División de Ciencias Básicas e Ingeniería, Universidad Autónoma Metropolitana-Iztapalapa, Apartado Postal 55-534, México D.F., 09340, Mexico
highlights • Diffusion across interfaces between materials with different diffusivity is studied. • Brownian particles run preferentially along the direction of decreasing diffusivity. • An uphill diffusion effect is found when particles move from high to low diffusivity layers.
article
info
Article history: Received 17 May 2013 Received in revised form 26 September 2013 Available online 22 October 2013 Keywords: Diffusion Heterogeneous systems Interface conditions Discontinuous interface density
abstract Diffusion in heterogeneous media (i.e., layers with discontinuous values of the diffusivity parameter) is found in different situations, including geology, catalyst, physiology and separation processes. This work considers an array of rectangular layers of different diffusivity as a simple model for studying diffusivity in heterogeneous media. Random walk simulations adapted for transport across interfaces between different materials are used for showing that heterogeneous media exhibit asymmetric transport in the sense that the diffusive flux depends of the transport direction. It is shown that Brownian particles run preferentially along the direction of decreasing diffusivity. An uphill diffusion effect is found when particles move from high to low diffusivity layers. In terms of Fick’s modeling, asymmetric transport is explained from the presence of a jump in the particle density at layer interfaces. © 2013 Elsevier B.V. All rights reserved.
1. Introduction Diffusion in domains with discontinuous diffusivity parameter (see Fig. 1) appears commonly in many scientific and technological fields, including geology [1–4], ecology [5–7], biomedical images [8], astrophysics [9,10] and molecular dynamics [11]. For instance, in ecology studies landscapes are considered as mosaics of patches of different habitat types with more or less well-defined patch boundaries [12]. In such case, the behavior of individuals in the vicinity of patch boundaries may be the determinant of overall movement behavior [13]. Also, interfaces between individual constituents play an important role in the effective transport properties of composite materials [14]. The modeling the diffusion process for systems involving discontinuous diffusivity requires the determination of matching conditions relating the limiting values on either side of the material interface of the particle density c (x, t ) and its spatial gradient. For the latter, imposing mass conservation across the interface fulfills the spatial gradient requirement. For the former, continuity of particle density c (0− , t ) = c (0+ , t ) is commonly assumed as motivated by intuition and physical experiments of diffusion through a near interface between two different media [2,3,15–18]. However, numerical
∗
Corresponding author. Tel.: +52 55 58044650; fax: +52 55 58044650. E-mail address: [email protected] (J. Alvarez-Ramirez).
1 Current address: Departamento de Biotecnología, Universidad Autónoma Metropolitana-Iztapalapa, Apartado Postal 55-534, 09340 México D.F., Mexico. 0378-4371/$ – see front matter © 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.physa.2013.10.027
194
J. Alvarez-Ramirez et al. / Physica A 395 (2014) 193–199
Fig. 1. Schematic description of a layered diffusion system with different diffusivity parameters.
stochastic simulations and experimental evidence indicated that the continuity of the particle density cannot be applied for heterogeneous systems with different value of the diffusivity parameter [19–25]. Theoretical work on random walk transport across interfaces with discontinuous diffusivity has indicated that particle density should meet a jump condition related to the ratio of diffusivities [26–28]. It was predicted that the probability density may have a point mass at the interface region, meaning that particles accumulate at the interface. Also, depending of the direction of flux, particles can move against density gradients in the interface vicinity, exhibiting a type of uphill diffusion induced by diffusivity discontinuity. A drift was detected of the particles’ individual positions in the direction of the diffusion coefficient gradient, in the absence of any external force or concentration gradient [29–32]. Some numerical results have pointed out that drift velocity and driving force are not directly proportional in the case of inhomogeneous suspensions, where a space dependent mobility induces an additional contribution to the drift velocity [33]. It has also been argued that the classical Fick law is not applicable in the vicinity of sharp inhomogeneities and that the Fokker–Planck equation should be conceded preference [34,35]. The present work considers the diffusion of particles in heterogeneous (layered) systems as described in Fig. 1. In a first step, in order to gain insight into the effects of the discontinuous diffusion parameter, random walk simulations adapted for discontinuous media [36,37] are performed. For D− ̸= D+ , it is found that the resident particle density can be considered as discontinuous at the interface, and the form of such discontinuity depends of the particle direction. In this way, the commonly used matching condition c (0− , t ) = c (0+ , t ) should be discarded from the modeling of the diffusion system. In order to close the problem, an interface jump condition is imposed [26–28], which is shown to be equivalent to the continuity of the chemical potential spatial gradient. This jump condition allows the recovery of the salient features near the interface as observed from the random walk simulations. In particular, in agreement with random walk simulations, asymmetric transport in the sense that the diffusive flux is dependent of the transport direction is shown to be exhibited by heterogeneous media. Also, an uphill diffusion effect is found when particles move from high to low diffusivity layers. The paper is organized as follows. Section 2 describes the numerical algorithm for performing random walk simulations across boundaries with discontinuous diffusivity. Section 3 discusses the numerical results in terms of the particle density jump condition at the system interface. Section 4 closes the work with some conclusions. 2. Random walk across the interface The motion of Brownian particles in a rectangular layered system like that in Fig. 1 is considered. Given that only the effects of the discontinuous interface are studied, the simulation of Brownian random walkers is reduced to only the horizontal direction, so the simulation problem is one-dimensional. In each layer, corresponding to high and low diffusivities, particles undergo Brownian motion governed by the Langevin equations dx dt
√ =
2Dξ (t )
with the diffusivity D given by D = D− , D = D+ ,
−L < x < 0 0 < x < +L.
Here, ξ (t ) is a zero-mean white Gaussian noise with autocorrelation function ⟨ξ (t ), ξ (t ′ )⟩ = δ(t − t ′ ). The physical domain is constrained to the interval x ∈ (−L, L). For each layer, the Langevin equation was numerically integrated by using
J. Alvarez-Ramirez et al. / Physica A 395 (2014) 193–199
195
the Euler algorithm. A particular situation arises when the particle crosses the interface located at x = 0. The simulation problem is that particle can a start at, e.g., the left of the interface (x < 0) where D = D− and end at the right of the interface (x > 0) with D = D+ . The question is, what diffusivity value should be used in√ the Langevin equation for accurate simulation of the physical phenomenon? The integration of the Langevin equation dx = 2Dξ (t ) in the vicinity of the interface follows dt the ideas of skewed diffusion [36]. For a given position x(t ) < 0, the procedure can be described as follows:
√ • Step 1. Compute the predicted position x(t + 1t ) = x(t ) + 2D− 1t ξ (t ). • Step 2. If x(t + 1t ) < 0, repeat Step 1. Otherwise, compute the time fraction that the particle spends moving in the (x(t )/ξ (t ))2 interface left-hand side by solving the equation 0 = x(t ) + 2D− 1tf ξ (t ). This gives 1tf = . 2D− • Step 3. Update the particle position with the remaining time after crossing the interface; namely, x(t + 1t ) = 2D+ (1t − 1tf )ξ (t ).
A similar procedure is implemented for particle position at the right-hand side of the interface, x(t ) > 0. The mean passage time was obtained as ensemble averages over NT = 105 trajectories starting in a domain boundary. Particle trajectory starts at the left boundary x = −L and finishes when the particle achieves the right boundary x = +L. The passage time is taken as the time the particle spends to move from the left boundary to the right boundary. Trajectories starting in the left boundary and returning to it are discarded for avoiding bias in the estimation of the mean passage time. To monitor the transport process, the variable of interest is the flux of particles, J. For a given trajectory, let TFPT be the first passage time (i.e., the time required by a trajectory to pass from the left to the right boundary). If TMFPT = ⟨TFPT ⟩ is the mean value of the first passage time obtained by averaging over NT trajectories, the flux can be expressed by J =
1 TMFPT
.
The flux J is an index of the ease with which the particles are transported along the layered system. Physically, the quantity J indicates the mean number of particles leaving the system domain. The flux J is a quantity of interest in diverse fields because it can be measured and, for instance in chemical reactors and membranes, reflect the ability of the physical system for passive mass transport. The effect of the transport direction (from low to high diffusivity and vice versa) in the resident particle density profile and the particle flux was evaluated by means of numerical simulations. To this end, let us introduce a diffusivity ratio as rD = Dhigh /Dlow , where Dhigh = max{D− , D+ } and Dlow = min{D− , D+ }, so that rD ≥ 1. The numerical simulations described below were performed for Dhigh = 1, L = 1 and 1t = 10−5 . The resident particle density profile is computed as c (x, t ) = N (x, t )/NT , where N (x, t ) is the number of particles at the position x for given time t. Recall that NT = 105 . For two different values of the diffusivity ratio rD and for the time t = 0.5, Fig. 2(a) and (b) present the resident particle density profiles when Brownian particles undergo low-to-high (LH) and high-to-low (HL) diffusivity transport, respectively. In the former case, D− < D+ while in the latter case D+ < D− . To model the physical transport of particles across the inhomogeneous system, all trajectories started in the left (resp., right) boundary for particles undergoing LH (resp., HG) diffusivity transport. The following can be noticed: (a) The resident particle density exhibits a sharp transition at the interface, and this transition can be considered as a discontinuity for space scales of the order of 0.025L. This result is in agreement with previous findings pointing out that the continuity of the particle density cannot be applied for macroscopic modeling of diffusion across interfaces [21–28]. (b) The geometry of the discontinuity depends of the order of the layers. That is, while the LH diffusion transport induces a monotonous discontinuity where the resident particle density decreases abruptly at the interface, the HL diffusion transport yields a non-monotonous discontinuity where the resident particle density increases abruptly at the interface. It is expected that the differences in the geometry of the particle density discontinuity has an important effect in the transport of particles across the system domain. Also for t = 0.5, Fig. 3(a) shows the particle flux J (rD ) normalized with respect to the particle flux for rD = 1. It is noted that the flux depends of the layer ordering. In fact, the HL mass flux is higher than the LH mass flux, indicating an asymmetry in the diffusion transport when interfaces are present. That is, the results in Fig. 3(a) indicate that mass flux runs preferentially along the direction of decreasing diffusivity. Fig. 3(b) displays the mass flux ratio γ (rD ) = JHL (rD )/JLH (rD ), where JLH (rD ) and JHL (rD ) are the LH and the HL diffusivity fluxes, respectively. It is shown that the diffusion transport asymmetry, quantified in terms of γ (rD ), increases with the diffusivity ratio rD , and such asymmetry is a consequence of the interface effects in the diffusion of Brownian particles. 3. Interface matching conditions For individual layers in Fig. 1, the average mass transport can be described by means of the mass balance in combination with the first Fick law. However, given the heterogeneous nature of the transport system, and motivated by the results obtained from random walk simulation, the following question can be posed: What interface (i.e., matching) conditions should be imposed in order to recover the behavior described in Figs. 2 and 3? Undoubtedly, the mass conservation across the interface should be imposed. However, given the particle density profile shown in Fig. 2, it is apparent that continuity of resident particle density cannot be imposed as an interface matching condition. To address this issue, consider the simple
196
J. Alvarez-Ramirez et al. / Physica A 395 (2014) 193–199
Fig. 2. Resident particle density profiles for (a), D− < D+ and (b) D− > D+ . Note the abrupt resident particle density change in a vicinity of the interface.
layered system in Fig. 1. The description of the transport phenomenon along the horizontal direction for individual domains can be made with the standard diffusion equation:
∂ c (x, t ) ∂ 2 c (x, t ) = D− , −L < x < 0 ∂t ∂ x2 (1) 2 ∂ c (x, t ) ∂ c (x, t ) = D+ , 0 < x < +L ∂t ∂ x2 where c (x, t ) is the mass density and the diffusivity parameters satisfy D− ̸= D+ . A complete modeling involves imposing transport matching conditions at the interface located at x = 0. Mass conservation is a mandatory interface matching condition that is established according to the Fick law of the form
− D−
∂ c (0+ , t ) ∂ c (0− , t ) = −D+ . ∂x ∂x
(2) ∂ c (0− ,t )
∂ c (0+ ,t )
It is noted that D− ̸= D+ implies gradient discontinuity at the interface. For instance, > ∂ x if D− < D+ . ∂x Given the parabolic nature of the system Eq. (1), a second interface matching condition is required for closing the model. Mass density continuity c (0− , t ) = c (0+ , t ) is commonly imposed as motivated by intuition and physical experiments of diffusion through a near interface between two different media [2,3,15–18]. However, Figs. 2 and 3 and reported numerical stochastic simulations and experimental evidence indicated that the continuity of the particle density cannot be applied for cases with different value of the diffusivity parameter [19–25]. As a consequence, the particle density should present a neat jump at the interface. Ovaskainen and Cornell [26] provided a rigorous proof that normal diffusion across interfaces should meet the following jump condition: D− c (0− , t ) = D+ c (0+ , t ).
(3)
In this way, since D− ̸= D+ one has that a mass density discontinuity is present at the interface (i.e., c (0− , t ) ̸= c (0+ , t )). Subsequently, Marseguerra and Zoia [27] showed that this jump condition is equivalent to v− c (0− , t ) = v+ c (0+ , t ), where v− and v+ are the velocities of the random walkers at the left- and right-hand side of the interface, respectively. It has been also shown [28] that the interface jump condition can be extended for anomalous diffusion transport. The interface matching ∂ ln c (0− ,t )
∂ ln c (0+ ,t )
= . Physically, the above relationship can be interpreted conditions (2) and (3) can be combined to give ∂x ∂x as the equality of the chemical potential gradient across the discontinuous interface [31,38].
J. Alvarez-Ramirez et al. / Physica A 395 (2014) 193–199
197
Fig. 3. (a) Diffusion flux as a function of the diffusivity ratio rD = Dhigh /Dlow when the transport goes from a low diffusivity layer to a high diffusivity layer, and vice versa. It is noted that the flux exhibits asymmetries depending on the ordering of the diffusion layers. (b) Diffusion flux ratio γ = JLH /JHL as a function of the diffusivity ratio rD .
It will be shown that the interface condition Eq. (3) recovers the transport asymmetry obtained by means of Brownian particle simulations, as already described by Figs. 2 and 3. To this end, for the sake of simplicity, the case of non-transient diffusion transport with fixed boundary particle density is considered. That is, the diffusion transport is governed by 0 = D− 0 = D+
d2 c (x, t ) dx2 2 d c (x, t )
,
−L < x < 0 (4)
,
dx2 with boundary conditions
0 < x < +L
c (−L) = cl
(5)
c (+L) = cr
and interface conditions given by Eqs. (2) and (3). Physically, Eqs. (1)–(3) can be seen as a model describing the diffusion transport across a membrane composed of two different materials. The solution to Eq. (1) is given by c (x) = a− x + b− , c (x) = a+ x + b+ ,
−L < x < 0 0 < x < +L .
(6)
The set of constants {a− , a+ , b− , b+ } is obtained after the application of the boundary and interface conditions. It can be shown that these constants are given by a− = a+ = b− = b+ =
r± cr − cl 2L r± cr − cl 2Lr± r± cr + cl 2 r± cr + cl 2r±
(7)
198
J. Alvarez-Ramirez et al. / Physica A 395 (2014) 193–199
where r± = D+ /D− . Note that r± = 1 implies that a− = a+ and b− = b+ . The following results can be established according to Eqs. (6) and (7):
• The interface mass density is discontinuous. Eq. (3) imposes a jump in the mass density at the interface. In fact, c (0− ) = b− and c (0+ ) = b+ . In addition, c (0− ) = r± c (0+ ). In this way, in agreement with the Brownian particle simulations exhibited in Fig. 2, c (0− ) < c (0+ ) for r± < 1, and c (0− ) > c (0+ ) for r± > 1 regardless of the left and right boundary particle densities cl and cr . This is in agreement with experimental findings that used two liquids having different gelatin concentrations and thus different diffusivity [34]. Such experiments showed the spontaneous appearance of a relatively sharp transition centered on the interface between the two liquids. The numerical results discussed in the above section indicate that the sharp transition is a geometric effect induced by diffusivity discontinuity. On the other hand, the interface matching condition given by Eq. (3) establishes that the mass density discontinuity arises for balancing chemical potential across the interface. D + cr • The mass flux is asymmetric. The flux across the layered system is given by J = D− cl − . Assume that cl > cr , so that the 2L particles are transported from the left boundary (i.e., x = −L) to the right boundary (i.e., x = +L). The mass flux when particles go from the high diffusivity layer to low diffusivity is given by JHL =
Dhigh cl −Dlow cr 2L
. Similarly, the mass flux when
particles are transported from the low diffusivity layer to high diffusivity is given by JLH = mass flux ratio γ = JLH /JHL is given by Dlow cl − Dhigh cr
γ (r D ) =
Dhigh cl − Dlow cr
=
cl − rD cr rD cl − cr
Dlow cl −Dhigh cr 2L
. One has that the
.
(8)
Recall that rD = Dhigh /Dlow . Let β = cl /cr be the boundary particle density ratio. The derivative of the flux ratio is
γ ′ (rD ) =
1 − β2
(β rD − 1)2
.
(9)
As in the Brownian particle simulations, it was assumed that β > 1. Then, one has that γ ′ (rD ) < 1, which is in agreement with the numerical findings shown in Fig. 3(b). In this way, the value of the mass flux depends of how layers are ordered with respect to the flux direction. Diffusion transport is faster for the HL transport direction than for LH transport direction. It is noted that this mass flux asymmetry is not found when continuous interface particle density is assumed. Analogously to electric rectifiers in electronic circuits and to the recently considered problem of heat transport rectifiers [39,40], the fact that mass is transported preferentially in one direction to the other suggests that layered material with discontinuous diffusivity parameter can be used as mass transport rectifiers. Here, particles move easily in one direction while finding high resistance in the reverse direction. Potential applications of mass transport rectification can be for partially impregnated catalyst where the catalyst pellet is composed of a mesoporous inert support containing islands or eggshells of active microporous support [41,42]. • Interface resistance can be negative. The mass flux across the layered system can be expressed as J =
1ceff
(10)
Reff
where 1ceff is an effective particle density gradient and Reff is the corresponding effective mass transport resistance. These functions are given by
1ceff = cl − r± cr
(11)
Reff = R− + Ri + R+
(12)
and
where the standard resistance for individual layers is given as R− = L/D− and R+ = L/D+ . Analogous to Kapitza interface resistance for heat transfer [43], Ri (r± ) represents the interface resistance for mass transfer and is given by
Ri = L
D+ − D− D+ D−
.
(13)
In this regard, the following can be pointed out: (a) With respect to the physical particle density gradient 1c = cl − cr , the effective particle density gradient 1ceff = cl − r± cr is reduced (resp., magnified) for r± > 1 (resp., r± < 1). Also, for non-layered (i.e., homogeneous) systems, the physical particle density gradient 1ceff (1) = cl − cr is recovered. (b) Eq. (13) reflects the particle resistance to jumping from one layer to the other. As expected, one has that Ri = 0 when D+ = D− . Interestingly, negative (resp., positive) interface resistance is obtained for D+ < D− (resp., D+ > D− ). Since cl > cr was assumed, negative (resp., positive) interface resistance is obtained for the HL (resp., LH) transport direction. In this way, increased mass flux in the HL transport direction is obtained by the combined effect of magnified particle density gradient and negative interface resistance. It should be pointed out that negative transport resistances have been also observed in heat transport through shaped graded materials [44]. Analogously to the results exhibited in Fig. 2(b),
J. Alvarez-Ramirez et al. / Physica A 395 (2014) 193–199
199
thermal energy is apparently transferred from a region of lower temperature to a region of higher temperature across composed materials of graded conductivity. Finally, an apparent paradoxical situation is found when particles move in the HL direction. Fig. 2(b) shows that the resident density exhibits a sharp increase in a vicinity of the interface. This indicates that particles diffuse up against density gradient (i.e., uphill diffusion transport). The interface is a geometric anomaly where the conditions for equilibrium transport are not met. The sharp increase of the resident particle density in the interface vicinity appears because particles are quickly transported in the high-diffusivity layer. In this way, particles are accumulated in the interface vicinity where particle transport faces a higher resistance. Uphill diffusion induced by discontinuous interfaces has been already observed in experiments for, e.g., boron transport in ultrashallow junctions [45,46] and carbon diffusion in steel layers [47]. It should ∂ ln c (0− ,t )
∂ ln c (0+ ,t )
= at the interface, particles are still transported along negative gradients of be remarked that, since ∂x ∂x chemical potential (i.e., partial molar free energy), preserving the thermodynamic constraints for diffusion processes [47,48]. 4. Conclusions Numerical simulations of Brownian trajectories in heterogeneous media (i.e., layers with discontinuous values of the diffusivity parameter) indicated that the commonly assumed condition of continuous interface particle density for diffusion transport cannot be applied since the particle density in a vicinity of the interface exhibits sharp positive and negative changes. Analogously to electric rectifiers, the mass flux is not symmetric in the sense that particles flow more easily in one direction than in the reverse direction. In fact, particle transport is faster when particles move from high diffusivity to low diffusivity layers than in the reverse direction. By imposing a jump condition for particle density at the interface, simple analytical results showed that particles are transported along chemical potential gradients, and this jump condition is able to explain the asymmetric mass flux of layered systems. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38] [39] [40] [41] [42] [43] [44] [45] [46] [47] [48]
F. Delay, Ph. Ackerer, C. Danquigny, Vadose Zone J. 4 (2005) 360. H. Hoteit, R. Mose, A. Younes, F. Lehmann, Ph. Ackerer, Math. Geol. 34 (2002) 435. E.M. LaBolle, J. Quastel, G.E. Fogg, J. Gravner, Water Resour. Res. 36 (2000) 651. H. Zhan, Zh. Wen, G. Huang, D. Sun, J. Contam. Hydrol. 107 (2009) 162. R.S. Cantrell, C. Cosner, Theor. Popul. Biol. 55 (1999) 189. W.F. Fagan, R.S. Cantrell, Ch. Cosner, Amer. Nat. 153 (1999) 165. J.T. Cronin, Ecology 84 (2003) 1506. E. Fieremans, D.S. Novikov, J.H. Jensen, J.A. Helpern, NMR Biomed. 23 (2010) 711. M. Zhang, Astrophys. J. 541 (2000) 428. A. Marcowith, F. Casse, Astronom. Astrophys. 515 (2010) A90. M. Mascagni, N.A. Simonov, SIAM J. Sci. Comput. 26 (2004) 339. ¯ A. Okubo, Diffusion and Ecological Problems: Mathematical Models, Springer Science + Business Media, 1980. H. Kaiser, Oecologia 56 (1983) 249. K.K. Chawla, Composite Materials: Science and Engineering, Springer, 2009. H.S. Carslaw, J.C. Jaeger, Conduction of Heat in Solids, Clarendon, Oxford, 1959. E.M. LaBolle, G.E. Fogg, A.F.B. Tompson, Water Resour. Res. 32 (1996) 583. E.M. LaBolle, J. Quastel, G.E. Fogg, Water Resour. Res. 34 (1998) 1685. J.Y. Parlange, F. Stagnitti, J.L. Starr, R.D. Braddock, J. Hydrol. 70 (1984) 251. O. Ovaskainen, S.J. Cornell, Proc. Natl. Acad. Sci. 103 (2006) 12781. M.Th. Van Genuchten, J.C. Parker, Soil Sci. Soc. Am. J. 48 (1984) 703. R.C. Schwartz, K.J. McInnes, A.S.R. Juo, L.P. Wilding, D.L. Reddel, Water Resour. Res. 35 (1999) 671. K.S. Novakowski, Water Resour. Res. 28 (1992) 2399. K.S. Novakowski, Water Resour. Res. 28 (1992) 2423. F.J. Leij, J.H. Dane, M.Th. van Genuchten, Soil Sci. Soc. Am. J. 55 (1991) 124. J.-L. Barrat, F. Chiaruttini, Mol. Phys. 101 (2003) 1605. O. Ovaskainen, S.J. Cornell, J. Appl. Probab. 40 (2003) 557. M. Marseguerra, A. Zoia, Ann. Nucl. Energy 33 (2006) 1396. N. Korabel, E. Barkai, Phys. Rev. E 83 (2011) 051113. P. Lançon, G. Batrouni, L. Lobry, N. Ostrowsky, Europhys. Lett. EPL 54 (2001) 28. P. Lançon, G. Batrouni, L. Lobry, N. Ostrowsky, Physica A 304 (2002) 65. M. Smerlak, A. Youssef, ArXiv Preprint. arXiv:1206.3441. F. Marchesoni, Materials 6 (2013) 3598. M. Yang, M. Ripoll, J. Chem. Phys. 136 (2012) 204508. B.P. Van Milligen, P.D. Bons, B.A. Carreras, R. Sánchez, Eur. J. Phys. 26 (2005) 913. F. Sattin, Phys. Lett. A 372 (2008) 3941. A. Lejay, G. Pichot, J. Comput. Phys. 231 (2012) 7299. L.E. Reichl, A Modern Course in Statistical Physics, Wiley & Sons, New York, 1997. S.R. De Groot, P. Mazur, Non-Equilibrium Thermodynamics, Dover Publications, 1984. W. Kobayashi, Y. Teraoka, I. Terasaki, Appl. Phys. Lett. 95 (2009) 171905. S. Narayana, Y. Sato, Phys. Rev. Lett. 108 (2012) 214303. S. Minhas, J.J. Carberry, J. Catal. 14 (1969) 270. A. Graviilidis, A. Varma, M. Morbidelli, Catal. Rev. Sci. Eng. 35 (1993) 399. E.T. Swartz, R.O. Pohl, Rev. Modern Phys. 61 (1989) 605. C.Z. Fan, Y. Gao, J.P. Huang, Appl. Phys. Lett. 92 (2008) 251907. H.C.N. Wang, C.C. Wang, C.S. Chang, T. Wang, P.B. Griffin, C.H. Diaz, IEEE Electron Device Lett. 22 (2001) 65. M. Ferri, S. Solmi, D. Giubertoni, M. Bersani, J.J. Hamilton, M. Kah, N.E.B. Cowern, J. Appl. Phys. 102 (2007) 103707. L.S. Darken, Trans. AIME 180 (1949) 430. D. Kondepudi, Introduction to Modern Thermodynamics, John Wiley & Sons, Ltd., Chichester, UK, 1998.