Ubiquity of anomalous transport in porous media: Numerical evidence, continuous time random walk modelling, and hydrodynamic interpretation

Anomalous transport in porous media is commonly believed to be induced by the highly complex pore space geometry. However, this phenomenon is also observed in porous media with rather simple pore structure. In order to answer how ubiquitous can anomalous transport be in porous media, we in this work systematically investigate the solute transport process in a simple porous medium model with minimal structural randomness. The porosities we consider range widely from 0.30 up to 0.85, and we find by lattice Boltzmann simulations that the solute transport process can be anomalous in all cases at high Péclet numbers. We use the continuous time random walk theory to quantitatively explain the observed scaling relations of the process. A plausible hydrodynamic origin of anomalous transport in simple porous media is proposed as a complement to its widely accepted geometric origin in complex porous media. Our results, together with previous findings, provide evidence that anomalous transport is indeed ubiquitous in porous media. Consequently, attentions should be paid when modelling solute transport by the classical advection-diffusion equation, which could lead to systematic error.

with β ≠ 1). These two typical features of anomalous transport are in great contrast with the prediction of the classical advection-diffusion equation (ADE). Various approaches have been proposed to model the anomalous transport process by taking into account a widely distributed heterogeneous velocity field 2,20-37 . It is perhaps natural to ascribe the heterogeneity of the velocity field, and the resultant anomalous transport, to the complex geometry of pore space, and it is indeed so in some cases. For example, anomalous transport can emerge due to the non-trivial complex structure (scale-free connection patterns) of the underlying network in which the transport process takes place 20 . Also, experimental as well as numerical evidence has shown that the type of the transport process can undergo qualitative changes as the pore space structure becomes more complex 9 . However, on the other hand, we also notice there has been experimental and numerical evidence that even for some structurally "simple" porous media, the solute transport process can be anomalous at high Péclet numbers 4,14,38,39 . In such cases, the porous media are regular or very weakly complex in void space geometry; there can hardly be any non-trivial spatial structure that is complex enough to induce a highly heterogeneous velocity field. Hence, the structural complexity does not seem to be the major source of the anomaly. Concerning these facts, one may wonder how much geometric complexity is needed to induce anomalous transport, and whether anomalous transport can be persistently observed in porous media with simple pore space structure.
To address these questions, we in this work study the transport process for a toy model of porous media that is constructed to be structurally as simple as possible. By doing so, one can exclude, or at lease greatly suppress, the influence of structural complexity on solute transport. While, at the mean time, the model still contains some minimal degree of randomness to make the results generalizable. Also, this model of porous media is similar to the setting in some laboratory experiments 4 . We then use the lattice Boltzmann method (LBM) 14,[40][41][42][43][44][45][46] to numerically solve the Navier-Stokes equations (NSEs) and the ADE in turn for each realization of the porous medium with a varying porosity φ. Scaling exponents of the concentration spatial moments (M 1,2 , see below) are calculated to compare with those obtained theoretically by the continuous time random walk (CTRW) theory.
We find that anomalous transport is astonishingly prevalent in our simple porous medium model when advection plays a dominant role. The anomaly exponent β ≠ 1 persists for porosities ranging from φ ≈ 0.30 up to φ ≈ 0.85. The observed scalings of M 1 and M 2 can be explained by the CTRW theory, which takes advantage of the statistics of the reciprocal steady-state velocity field, i.e., the statistics of 1/u, where u is the magnitude of the flow velocity u. In fact, t = 1/u is interpreted as the waiting time at Pe = ∞ in the CTRW theory. We use a mixture model to describe the probability density function (PDF) of t, denoted w(t). We argue conceptually it is helpful to decompose the pore-scale flow field into two distinct parts with the first being a globally uniform velocity field and the second being a locally fluctuating velocity field. Physically, the first part is essentially a manifestation of Darcy's law and helps to explain the numerically observed linear scaling M t 1 ; the second part gives rise to the anomalous feature of the transport process and controls the scaling of M 2 . We find for the whole range of porosity considered in this work. The theoretical prediction that β = 3 − α is confirmed, which is an evident signature of anomalous transport.
The prevalence of anomalous transport in our simple model strongly indicates that such an anomaly is ubiquitous in general porous media. By showing an extreme case of anomalous transport in two dimensional Poiseuille flow, we further argue that the emergence of anomalous transport is a result of the joint action of hydrodynamics and pore space geometry. Qualitatively, the physical requirement of a no-slip boundary condition results in a quasi-parabolic velocity profile in throats as a solution to the NSEs. Combinations of throats constitute the main body of fluid pathways in porous media, which eventually lead to a heterogeneous flow field and hence anomalous transport 39 . For structurally simple porous media, hydrodynamics plays a crucial role in inducing non-uniform velocity profiles, while the influence of geometry manifests itself by enhancing the heterogeneity of the flow field as structural complexity grows. In our work, if we consider porosity as some "mean-field" measure of the pore space complexity, then we will see as the porosity φ is decreased, the exponent α follows, thus leading to a more heterogeneous flow field.
One implication of our results is that the conventional ADE on the Darcy scale may not be qualitatively adequate in describing transport dynamics in general porous media, due to the ubiquity of anomalous transport at high Péclet numbers. Thus, to avoid systematic error, approaches such as the CTRW theory may be adopted as more proper modelling tools 2 .

Results
observation of anomalous transport in model porous media. The methods to generate the porous medium and to simulate fluid flow and solute transport are detailed in Methods. After generating the desired porous medium and obtaining the steady-state fluid velocity field, we investigate how the solute concentration profile is changed with time. In particular, to quantify the statistical features of the profile, by which the anomalous nature of the process can directly be recognized, we calculate the first and second central moments of the concentration field in the x direction, denoted M 1 and M 2 , respectively. Numerically, they are computed as Extensive numerical simulations are performed to investigate how M 1 and M 2 behave as t is increased, at various φ and Pe. We mainly change φ from 0.30 to 0.85. We adopt Pe = 50 and Pe = 0.5 in our simulations to qualitatively represent two limiting situations: In the first case, advection is dominating; in the second case, the effect of diffusion is more important. An astonishing finding is that anomalous transport is unexpectedly prevalent in the first case, characterized by an anomaly exponent β ≠ 1 and An illustrative example is presented in the following to show the emergence of anomalous transport in the advection dominant situation, in which the domain size is 500 × 100 and φ ≈ 0.8. Other parameters are: kinematic viscosity ν = 1.0, the density at the outlet ρ east = 1.0, and the density drop Δρ = 0.01. l = 9 is used to generate the solid block. We in Fig. 1(a) plot the steady-state flow field in this case, whose spacial resolution is five cells, i.e., the fluid velocity, as denoted by the arrow, is plotted every five cells. The seemingly narrow channels, however, almost always contain enough cells to ensure a reliable LBM simulation 40 . The mean velocity in the x direction is = . × − u 9 2 10 x 5 , hence the Mach number is = . × −  Ma 1 6 10 1 4 . We also analyze the statistics of the fluid velocity field and plot in Fig. 1(b) the PDF w(t) of t ≡ Δx/u = 1/u, where Δx is the spacial step in LBM. We note that w(t) has a peak at some t * . For t < t * , w(t) drops drastically, while for t > t * , w(t) displays a fat tail, scaling as α − − w t t ( ) 1 with α ≈ 1.36 roughly on a time interval t ∈ (10 4 , 10 5 ). If advection dominates the transport dynamics, then intuitively, this tail means solute particles may spend a much longer time at some cells than at others, and the solute particles are likely to move along the pressure drop direction at a non-uniform pace, demonstrating a spatially inhomogeneous concentration profile. However, when diffusion outweighs advection in affecting the transport process, a more homogeneous concentration field is expected to be observed. These are shown in www.nature.com/scientificreports www.nature.com/scientificreports/ In Fig. 2, we plot the concentration profiles, as well as the corresponding spatial distribution of accumulated concentration in the x direction ≡ ∑ c x i c i j ( ( )) ( , ) j , at times 0.1 × t est and 0.3 × t est , respectively. t est roughly denotes the time for a solute particle with average velocity to reach the boundary (see Methods). The non-uniform concentration profile is evident. As advection is dominant, solute particles that move along some streamline to reach a stagnant zone (  u u / 1 x locally) will get stuck there and only have very low probability to escape. This is because the diffusion effect is too weak to induce effective particle transitions from such stagnant zones to regions where u is large via the transport along some adjacent streamline. Such suppression of particle transition between nearby streamlines is supposed to be the main reason why the solute concentration is distributed heterogeneously at high Pe.
However, when the diffusion effect is strong on transport, the spatial distribution of the concentration field is qualitatively different. In this situation, molecular diffusion will help solute particles hop out of stagnant zones, Figure 1. (a) The steady-state velocity field u(x) for a porous medium with φ ≈ 0.80; (b) the associated waiting time PDF w(t) at Pe = ∞ is obtained by statistics of t = 1/u, which displays a fat tail that scales as ~− . t 2 36 when t ∈ (10 4 , 10 5 ). Throughout this work, the spacial resolution of the velocity field is five cells, i.e., the velocity at a position as represented by an arrow is plotted every five cells. The length of an arrow is proportional to the corresponding u.
t est roughly denotes the time for a solute particle with average velocity to reach the boundary (see Methods for its precise definition).
www.nature.com/scientificreports www.nature.com/scientificreports/ and overall the concentration profile will become spatially homogeneous, as shown in Fig. 3. Two snapshots of the concentration field for Pe = 0.5 at times 0.1 × t est and 0.3 × t est are also plotted, respectively, together with c(x). Compared with Fig. 2, it is even visually clear that solute particles are transported in a more uniform way. Note that M 2 reflects the "spread" of the distribution of solute concentration in space, hence at a given time t a larger M 2 (t) quantitatively corresponds to a more homogeneous transport process than a smaller M 2 (t) does. At 0.1 × t est or 0.3 × t est , M 2 under Pe = 0.5 is way greater than that under Pe = 50 (Actually, the former is always greater than the latter in our simulations). That is to say, diffusion enhances the transport process. Although random hop rates between nearby streamlines may be the same 47 , the velocity fields along nearby streamlines are by no means always similar. Streamlines undergo a sudden "compression" when entering a narrow throat, which is a typical structure in porous media. In the outside wider region, the velocity gradient between nearby streamlines is much smaller than that within the throat. If a particle hops from one streamline which has a low velocity in the throat (near the solid boundary) to another streamline which has a high velocity in the throat (near the center of the throat), then it has high probability to be transported quickly forward, without hopping back to the slower streamline in the throat. On average, as long as the particles hop to the "express way, " they leave the stagnant zone in a short time. This is why diffusion speeds up the transport through throats. When particles are trapped within some region almost enclosed by solid cells, diffusion is the only possible mechanism for particles to get out. So, as the diffusion effect is increased, the concentration profile becomes more homogeneous.
Next, we investigate the transport process more quantitatively by tracking the time evolution of M 1 and M 2 . In Figs 4 and 5, we plot how they evolve for Pe = 50 and Pe = 0.5, respectively. In the first case, we find that   with β = 1.64 for t ∈ (4 × 10 4 , 4 × 10 5 ), which is a signature of anomalous transport. (But β seems to begin decreasing after t = 5 × 10 5 due to the finiteness of Pe in simulation.) It is worth noting that there seems to be a quantitative relation between w(t) and M 2 that β ≈ 3 − α. In the latter case of Pe = 0.5, since the concentration profile is rather "normal, " then not surprisingly, we find both M 1 and M 2 scale linearly with time, and the transport is Fickian. We also investigate the scaling behaviors of M 1 and M 2 under various choices of Pe, respectively. It turns out that M t 1 is always present, whereas M 2 will undergo a gradual transition from M t 2 to . M t 2 1 64 as Pe is increased. We plot in Fig. 6 how M 2 changes with time t as Pe is increased from 0.5 to 50. From Fig. 6, we can also see that when Pe is higher than 20, there are some numerical errors in the early stage of transport (inherent in the LBM when τ c approaches 0.5, see Methods), however, the scaling of M 2 is clearly not influenced at later times.
Actually, we have also found the scaling relations for many other realizations of the porous media with various φ. Anomalous transport characterized by such relations is thus prevalent at high Pe in our model, at least for some time interval t ∈ (t lb , t ub ) with ≈t 10 10 ub 5 6 typically.
scaling relations at high péclet numbers. In the CTRW theory, a particle hops from one site to another based on the transition rule that the transition distance is drawn from a PDF λ(r) and the waiting time between successive transitions is drawn from a PDF w(t). In our setting of the initial condition (see Methods), the master equation of the transition process reads t) is the PDF of a particle just arrives at x at time t. The PDF of finding a particle at x at time t is being the probability that during a period of time t, a particle does not leave x. Since we are interested in the transport process in the direction of pressure drop, the vectorial x is simplified as x in the following, "intergarting out" the information in the y direction. By performing Fourier-Laplace transform to W(x, t) 2,31,32 , we obtain In this work, λ(k) is somewhat trivial, because a particle can only hop to adjacent cells, and is finite for any non-negative integer n. Then as k → 0 we have λ(k) ≈ 1 + ikl − k 2 σ 2 /2, where l = 〈x〉 λ and σ 2 = 〈x 2 〉 λ . On the contrary, w(z) is non-trivial, we notice that w(t) in our work can be effectively written as a sum of two distributions: where w 1 (t) represents a distribution (for example, approximately a truncated Gaussian) that mainly accounts for the small-t part of w(t), w 2 (t) features the long-time tail of w(t), and p, q are some positive weighting parameters with p + q = 1. We are here not aiming at an accurate decomposition of w(t); by writing equation (1), we emphasize the fact that within the time interval of numerical simulations, neither w 1 nor w 2 is overwhelmingly dominant. Heuristically, www.nature.com/scientificreports www.nature.com/scientificreports/ On the contrary, w 2 (t) is characterized by the t −1−α tail with 1 < α < 2, and for small z,  31 . Also note that k 0 According to equation (2), we , and as a result , whose inverse Laplace transform yields where Γ is the gamma function. In the long time limit, the term lAt dominates, hence 1 is expected to be observed. Similarly, 〈x 2 (z)〉 can be obtained as Then also by equation (4), it is straightforward to obtain 〈 〉≈ Combining the results for 〈x 2 (t)〉 and 〈x(t)〉, we arrive at 〈Δ , and finally the scaling relation of M 2 is obtained as Therefore, the scaling relations (5) and (6) are in accord with numerical results. It is worth noting that, however, such scalings are derived under the condition that z t 1/ ub , and they are only expected to be observed for time up to t t ub . For  t t ub , since Pe is not strictly infinity, the diffusion process begins to take effect and the scalings for normal transport processes will emerge that M t 1,2 . This may be the reason why there are deviations of β from 1.64 as time is increased in Fig. 4.
After establishing the relation between α and β, we can understand the prevalence of anomalous transport with the help of w(t). Keeping other major simulation parameters unchanged, we systematically investigate w(t) by varying porosities φ and solid block sizes l × l. In Fig. 7, typical w(t)'s are plotted for φ ∈ [0.30, 0.85] and l ∈ [9,17]. In each case, w(t) approximately displays a fat tail α − − t 1 . Also plotted are two lines with slope −2 and −3 on a log-log scale, respectively, to guide eyes. We can see 1 < α < 2 holds in all the cases we consider. In fact, according to the above theoretical work, as long as the decomposition of w(t) into two distinct parts is valid, then for any 1 < α < 2, the scaling relations (5) and (6) are expected to be observed. This implies that for a wide range of porosity, the transport dynamics at high Pe is anomalous, which originates from the fat tail of w(t). We also observe that α seems to decrease as φ is decreased. This reflects the influence of structure, since the porosity can be seen as a basic mean-field description of pore space geometry. In some sense, as φ is decreased, the pore space becomes more "complex, " and the transport process becomes more "anomalous. "

Discussion
One issue we have not addressed is why w(t) typically displays a fat tail. The CTRW theory only takes advantage of this property to predict the scalings of anomalous transport, but itself does not provide any physical explanation why such a tail exists. At high Pe, w(t) is determined by the statistics of 1/u, and if we denote the PDF of u as f u (u), then w(t) = f u (1/t)|du/dt| = f u (1/t)t −2 . So, as long as f u (u) is a mildly increasing function of u on some interval that can be well approximated by α− u 1 , the t −1−α scaling will be observed. This might suggest that the fat tail in w(t) is indeed not that unusual. , where r is the distance to the center of the solid disk, and θ is the angle with respect to the direction of the unperturbed incoming flow with velocity u in . The waiting time t u 1/ , and the long-time tail of w(t) is mainly contributed by  u u / 1 in . Define ε = r − 1, then by expanding u(r) in the vicinity of r = 1, we have ε If a site is picked out at random and θ is uniformly distributed in (0, 2π), note that 1dxdy = rdrdθ = (1 + ε)dεdθ, then we have f ε (ε) ∝ (1 + ε), and consequently , and therefore | | w t f y dy dt ( ) ( ) / y , where f y (y) = 1 is the PDF of a randomly chosen y. We again find 3 . Quantitatively, we solve this problem using the LBM simulations. We also set Δx = Δt = 1 and the domain size n × m is 100 × 50. We choose ν = 1, ρ east = 1, and Δρ = 0.001. The no-slip condition is applied to the north and south boundaries, and the constant-pressure condition is applied to the west and east boundaries when simulating the NSEs. While, we apply the adsorption condition to the west and east boundaries when numerically solving the ADE. As the steady-state is achieved, = . × − u 6 38 10 x 4 . Then, we uniformly distribute solute particles in the y direction, and they are transported under the joint action of advection and diffusion. We in this case define the Péclet number as = u n D Pe / x m . Following similar theoretical analysis as shown above, we have at Pe = ∞ that While for M 2 , note that as 2 , then from the CTRW theory one expects M t t / log ( ) 2 2 4 32 when t is much more than 1/u 0 . In Fig. 8, Pe is set to be 100 and we clearly observe the good agreement between theory and simulations. In particular, M 2 becomes well represented by the theoretical scaling behavior when × Δ > t x u 10 / 10 x 4 . From the two examples presented above, we can see even for "normal" situations in fluid mechanics, the fat tail of w(t) exists at high Pe. In our porous medium model, the flow confronting a solid block and between two solid blocks can be approximated by these two cases. The above analyses might still be qualitatively reasonable, and at least partially justify that the ubiquitous fat tail of w(t) is more of a hydrodynamic than a geometric origin. In a recent work 39 , a similar analysis was made to show quantitatively that anomalous transport is present in weakly complex porous media and can be predicted by incorporating the Poiseuille flow field with a set of throat sizes that follow a power-law distribution, consistent with our results. We thus conclude that, under the joint influence of hydrodynamics and geometry, anomalous transport is indeed ubiquitous in porous media, especially at high Pe.
In this work, we mainly focus on transport dynamics in model geological systems where the heterogeneous flow field in pore space gives rise to anomalous transport. The setting of our model is relevant in such fields as petroleum engineering and groundwater science. It is worth noting that aside from geological systems,   www.nature.com/scientificreports www.nature.com/scientificreports/ anomalous diffusion is also observed in biological systems where, however, the main sources of the anomaly are the macromolecular complexity and spatiotemporal membrane heterogeneity 49 . Geometry-induced anomalous diffusion can be found in other cases as well 50,51 . In many biological systems, the spatial scale of interest is typically − − 10 10 m 9 6 , and experiments as well as numerical simulations can be performed using the single particle tracking technique to accurately capture each particle's dynamics. As a result, important statistics of a single-particle trajectory can be obtained for analyzing the time-averaged mean squared displacement. While in geological systems, although simulations can be done for rock samples of size . Hence, a macroscopic ensemble-averaged description of the process by CTRW modelling or a fractional ADE is perhaps more preferred in this sense. In this work, by using the LBM to directly simulate the evolution of the concentration field (equivalent to the single-particle PDF), we deal with the particle ensemble from the very beginning and neglect information such as the two-point time correlation. Consequently, our CTRW analysis only uses the most basic model and does not address issues like weak ergodicity breaking and ageing [52][53][54] , which are important in the analysis of anomalous diffusion in biological systems and surely worth investigating in geological systems as well. However, that will involve numerical algorithms totally different than the LBM, and we wish to leave these for future study.

Methods
Generation of porous Media. A porous medium in this work is composed of n × m cells; each is of unit size. The state of a cell indexed by (i, j) with i = 1, …, n and j = 1, …, m is denoted s(i, j) with s(i, j) = 1 for a solid cell and s(i, j) = 0 for a fluid cell. To achieve a statistically homogeneous porous medium, we divide the n × m region into = × n n m S 20 20 subregions, each of the size 20 × 20. In every subregion, a solid block consisting of l × l solid cells is generated randomly in position with probability φ/(1 − l 2 /n S ). For a given l, the minimum achievable porosity is 1 − l 2 /n S . By adjusting l from 9 to 17, we generate porous media with porosities ranging from roughly . . 0 30 0 85 in this work. Adopting l in this range also ensures the Knudsen number Kn to be small so that the LBM simulation results are reliable 40 . After generating a porous medium from some given porosity φ, we calculate the realized porosity is the number of solid cells. If |φ real − φ| < Δφ th , where Δφ th is a prescribed threshold and is chosen to be less than 0.03 in our simulation, we will proceed with the generated porous medium; otherwise, we will re-sample until the criterion is met.
By the generating rule above, the porous media thus constructed are statistically homogeneous. We use such a simplest model to suppress the structural effect on transport as much as possible, while at the mean time a minimal level of randomness is maintained to make the results more generalizable. simulation of Fluid Flow. Once the porous medium is generated, we numerically solve the NSEs by the standard single relaxation time D2Q9 scheme of the LBM 14,40,41 : ξ , where the subscript i = 0, …, 8, denoting the i th direction, f i is the single-particle distribution function, f i eq is the local equilibrium distribution function, τ is the dimensionless relaxation time, Δt = 1 is the time step used in this work, and ξ i is the discrete velocity with ξ 0 = (0, 0) T , ξ i = (cos((i − 1)π/2), sin((i − 1)π/2)) T for i = 1, 2, 3, 4, and ξ π π π π = + − + − i i 2 (cos( /4 ( 5) /2), sin( /4 ( 5) /2)) i T for i = 5, 6, 7, 8. We also choose 43  www.nature.com/scientificreports www.nature.com/scientificreports/ ν = (τ − 0.5)/3 is the kinematic viscosity, and the pressure p = ρ/3. Up to second order accuracy, this set of equations is equivalent to the classical NSEs: ∇ ⋅ = u 0, and ν + ⋅ ∇ = − + ∇ In the two-dimensional domain, boundary conditions are set as in a previous work 14 . We keep a constant pressure difference Δp between west (inlet) and east (outlet) boundaries; this is numerically realized by keeping a constant density difference Δρ 44 , due to the linear dependence of p on ρ. North and south boundary conditions are periodic. The no-slip condition is applied to the interface between fluid and solid cells with a second order accurate bounce-back method 45 . Initially, the fluid is set to be at rest. After performing the LBM simulation for a transient period of time, we obtain a steady-state velocity field; the steady state is considered to be reached when the criterion < simulation of solute transport. Having obtained the steady-state velocity field, we proceed to simulate the solute transport process, which is also done by the LBM with a single relaxation time D2Q5 scheme 14 : where i = 0, …, 4, denoting the i th direction, g i is the single particle distribution function, g i eq is the local equilibrium distribution function, τ c is another dimensionless relaxation time. g i eq is chosen as 14 ω ξ = + . ⋅ . This numerical scheme in the small Knudsen number limit leads to 14 + ∇ ⋅ = ∇