Abstract
We study first passage behaviors in the flow through three-dimensional random fracture networks. Network and flow heterogeneity lead to the emergence of heavy-tailed first passage time distributions that evolve with increasing distance between the start and target planes, and transition toward stable laws. Analysis of the spatial memory of the first passage process shows that particle motion can be quantified stochastically by a time domain random walk conditioned on the initial velocity data. This approach identifies advective tortuosity, the velocity point distribution and the average fracture link length as key quantities for the prediction of first passage times. Using this approach, we develop a theory for the evolution of first passage times with increasing distance between the start and target planes and the convergence towards stable laws.
| Original language | English (US) |
|---|---|
| Article number | 248501 |
| Journal | Physical review letters |
| Volume | 123 |
| Issue number | 24 |
| DOIs | |
| State | Published - Dec 9 2019 |
Bibliographical note
Funding Information:https://orcid.org/0000-0002-4224-2847 Hyman Jeffrey D. 1 ,* https://orcid.org/0000-0002-3940-282X Dentz Marco 2 https://orcid.org/0000-0003-2327-3725 Hagberg Aric 3 https://orcid.org/0000-0002-4961-6899 Kang Peter K. 4 Computational Earth Science Group(EES-16), Earth and Environmental Sciences Division, 1 Los Alamos National Laboratory , Los Alamos, New Mexico 87545, USA 2 Spanish National Research Council (IDAEA-CSIC), Barcelona 08034, Spain Computer, Computational, and Statistical Sciences Division, 3 Los Alamos National Laboratory , Los Alamos New Mexico 87545, USA Department of Earth and Environmental Sciences, 4 University of Minnesota—Twin Cities , Minneapolis, Minnesota 55455, USA * Corresponding author. jhyman@lanl.gov 9 December 2019 13 December 2019 123 24 248501 16 April 2019 © 2019 American Physical Society 2019 American Physical Society We study first passage behaviors in the flow through three-dimensional random fracture networks. Network and flow heterogeneity lead to the emergence of heavy-tailed first passage time distributions that evolve with increasing distance between the start and target planes, and transition toward stable laws. Analysis of the spatial memory of the first passage process shows that particle motion can be quantified stochastically by a time domain random walk conditioned on the initial velocity data. This approach identifies advective tortuosity, the velocity point distribution and the average fracture link length as key quantities for the prediction of first passage times. Using this approach, we develop a theory for the evolution of first passage times with increasing distance between the start and target planes and the convergence towards stable laws. U.S. Department of Energy 10.13039/100000015 89233218CNA000001 Los Alamos National Laboratory 10.13039/100008902 20180621ECR 20170103DR Triad National Security Basic Energy Sciences 10.13039/100006151 H2020 European Research Council 10.13039/100010663 Korea Environmental Industry and Technology Institute 10.13039/501100003654 Ministry of Environment 10.13039/501100003562 2018002440003 Seventh Framework Programme 10.13039/100011102 617511 First passage times play a central role in a wide variety of processes that depend on or are conditioned by the notion of first encounters including chemical transformations in fluctuating environments [1] , solute migration in geological repositories [2] , transport in brain microvascular networks [3] , and optimal search strategies [4,5] . However, the identification of measurable media properties and quantification of their control on first passage times, which is critical for transport predictions in complex media, is often problem dependent and illusive. For hydrodynamic transport in subsurface fracture networks, and natural and engineered media in general, spatial heterogeneity is the determining factor. Natural fractured media exhibit multiscale features spanning several orders of magnitude [6] . It has been ubiquitously found in field and numerical experiments that first passage times in such media are characterized by heavy tails, and that transport is typically non-Fickian or anomalous [7–14] . These behaviors have been modeled and interpreted through the use of continuous time and time-domain random walks [15–19] , which account for spatial heterogeneity through a distribution of mass transfer timescales, as well as multirate mass transfer approaches [20,21] , which model transport under mass transfer between fast mobile and slow immobile spatial continua. The complexity of random fracture networks makes it difficult to link structural and hydrodynamic properties to large scale transport behaviors [22–24] . Observed power-law tails in first passage time distributions indicates there may be universal behavior [7,8,10] , while observed dependencies on the injection conditions may indicate otherwise [25,26] . In this Letter, we address these questions for hydrodynamic transport in three-dimensional random fracture networks by studying the evolution of first passage time distributions with increasing distance between the start and target planes. We investigate spatial memory in terms of streamwise particle velocities, and use this information to construct a stochastic particle model that accounts for both universal and nonuniversal aspects of first passage times as a function of the distance between the start and target planes. The model is used to develop a theory that provides conditions on the Eulerian velocity distribution for first passage times to converge to a stable density. We verify these predictions through comparison with the model and discrete fracture network (DFN) particle tracking simulations. We consider a stochastically generated DFN using the dfnworks software [27,28] . Fracture radii are sampled from the truncated power-law distribution p r ( r ) ∝ ( r / r 0 ) - 1 - γ for r 0 ≤ r < r u , where r 0 is the minimum and r u the maximum radius. Geological fractured media exhibit hierarchical structures due to the interaction of different fracture growth processes that gives rise to length distributions that can be characterized by power laws [6] . We set here γ = 2.6 , r u = 10 2 r 0 , and consider a cubic domain of size equal to r u . Apertures b are constant within each fracture but vary between fractures being positively correlated to the fracture size [29] . The resulting network shown in Fig. 1 is semigeneric in that it is not meant to represent a particular field site, but values are selected to mimic the characteristics of large scale crystalline rock [6,30] . The dimensionless density of the network is ≈ 14 times larger than the critical percolation density ensuring that there are multiple paths between the flow boundaries [31,32] . 1 10.1103/PhysRevLett.123.248501.f1 FIG. 1. Illustration of the stochastically generated three-dimensional discrete fracture network. Fractures are colored according to their radii, darker colors indicating larger fractures. Pressure gradients in flows through low-permeability fractured rock are typically small enough that the flow is laminar and the Stokes equation is an appropriate model [33] . Furthermore, the large contrast between fracture length and aperture, typically several orders of magnitude, allows for the Stokes equation to be reduced to the Reynolds equation [34] . Flow through the DFN is driven by a constant mean pressure gradient ∇ P across the domain aligned with the 1 axis of the coordinate system. In the following, flow velocities are scaled by the characteristic velocity v 0 = b 0 2 | ∇ P | / 12 μ , where b 0 is the minimum fracture aperture and μ dynamic viscosity. Note that in the Stokes regime velocities scale linearly with the pressure difference and thus with the characteristic velocity. Details of network generation, numerical flow simulations, and specific flow conditions are provided in the Supplemental Material [35] , which includes Refs. [36–43] . Advective transport through the DFN is simulated using particle tracking. The trajectory x ( t ; a ) of a particle starting at a at time t = 0 is given by the kinematic relationship d x ( t , a ) d t = u [ x ( t , a ) ] . (1) where u ( x ) is the Eulerian velocity field. Although apertures are uniform within each fracture, u ( x ) varies within each fracture plane due to the location of intersections, boundary conditions, and fracture position within the entire network. One million particles are inserted and tracked through the system. The initial particle distribution is flux weighted, this means the particle density at a position in the injection plane ρ ( a ) is weighted by the magnitude v e ( a ) = | u ( a ) | of the local flow velocity. At fracture intersections, particles are distributed proportional to the outgoing fluxes. This means a complete mixing rule is applied. The first passage time to cross a target plane at a linear distance x 1 from the inlet is t a ( x 1 ; a ) = min [ t | x 1 ( t , a ) ≥ x 1 ] . The probability density function (PDF) f ( t ; x 1 ) of first passage times sampled over the ensemble of particles is f ( t ; x 1 ) = ∫ d a ρ ( a ) δ [ t - t a ( x 1 ; a ) ] , (2) where δ ( t ) is the Dirac delta. Figure 2 shows f ( t ; x 1 ) from the DFN simulations at three target planes within the domain (solid lines). Time is rescaled by the characteristic advective time τ v = r 0 / v 0 . We observe relatively abrupt early arrivals at distances close to the start plane, which become smoother as distance increases. The late time behavior is characterized by power-law tails as t - 1 - β with β = 3 / 2 at all distances. 2 10.1103/PhysRevLett.123.248501.f2 FIG. 2. First passage time distributions at (black) x 1 = 25 r 0 , (green) 50 r 0 , and (orange) 1 0 2 r 0 from the start plane. Time is rescaled by the advection travel scale τ v = r 0 / v 0 . Solid lines denote the direct numerical simulations, symbols are the predictions from the stochastic model, the dashed lines shows the t - 1 - β long time scaling with β = 3 / 2 . In order to understand these behaviors, we examine the individual particle motion, and construct a model that accounts for the physical processes at the origin of the observed power-law scaling. Figure 3 shows a particle’s velocity magnitude and fractures visited along the trajectory. This particle passes through 16 fractures that are denoted with a different shaded band. On the subfracture scale, the velocities fluctuate moderately. In contrast, major changes in velocity typically occur as a particle transitions between fractures. Thus, the number of major velocity transitions a particle makes depends on the number of fractures through which it passes, a topological and geometric property of the network. 3 10.1103/PhysRevLett.123.248501.f3 FIG. 3. A particle’s velocity magnitude v s sampled along a trajectory. The particle passes through 16 fractures. Each fracture, with different radii and permeability, is denoted with a different shaded band. Particle velocities fluctuate moderately on the subfracture scale, but there are major velocity changes as particles transition between fractures. These observations suggest a stochastic approach that models particle transitions over the characteristic spatial fluctuation scale of particle velocities. Velocities decorrelate as particles transition between fractures. Thus, the velocity correlation along pathlines can be estimated by computing the distance between centers of intersections within fracture planes. For this network, the mean distance between intersections is 153 m. This hydrodynamically independent value sets the scale of velocity transitions denoted by l c . Thus particle motion can be represented by a time-domain random walk (TDRW) type approach that advances the particle position by the constant distance l c along the trajectory, for which it takes the transition time τ = l c / v s , where v s is the velocity magnitude sampled equidistantly along a particle path [44] . This approach assumes Lagrangian ergodicity, which means that the velocity statistics sampled between particles and along streamlines are the same. The steady state distribution p s ( v ) of v s is characterized by two power-law regimes, which can be attributed to primary flow channels where the highest velocities exist and secondary structures where low velocities are present [45] . Such a distribution is well modeled by the Burr distribution p s ( v ) = α v c ( v / v c ) β - 1 [ 1 + ( v / v c ) β ] α / β - 1 , (3) with α = 1.7 , β = 1.5 and v c = 2.48 v 0 , see Fig. 4 . Under ergodic conditions, p s ( v ) is related to the PDF p e ( v ) of the Eulerian velocity magnitude v e ( x ) through flux weighting as p s ( v ) = v p e ( v ) / ⟨ v e ⟩ [26,44,46] . This relationship connects flow and transport attributes. Furthermore, under ergodic conditions the PDF of initial velocities p 0 ( v ) is equal to the steady velocity PDF p s ( v ) [44] . The differences illustrated in Fig. 4 are due to the finite size of the start plane. The particle velocity statistics evolve from the initial to the steady state distribution. Nonetheless, both p s ( v ) and p 0 ( v ) show the power-law behavior ∝ ( v / v c ) β - 1 at low velocities, which determines the late time scaling of the particle travel time distributions, as shown in the following. 4 10.1103/PhysRevLett.123.248501.f4 FIG. 4. (top panel) Illustration of (empty circles) steady Lagrangian velocity PDF p s ( v ) , (full circles) initial velocity PDF p 0 ( v ) , and (solid line) Burr distribution (3) . (bottom panel) Tortuosity as a function of streamwise distance. The full circles correspond to the tortuosity values at x 1 = 25 , 50, 10 2 r 0 . Note that this approach models the stochastic particle motion along a trajectory, which is tortuous as a result of the complex network and flow structure. This means the distance traveled is not equal to the streamwise distance between the start and the target planes. We account for this feature in terms of the advective tortuosity χ ( x 1 ) , the ratio between the average trajectory length at a given linear distance, which is illustrated in Fig. 4 . In the limit of x 1 → ∞ , χ ( x 1 ) converges to χ ∞ = ⟨ v e ⟩ / ⟨ u 1 ( x ) ⟩ [47] , which here is χ ∞ = 2.7 . The average trajectory length from the start to the target plane at x 1 = 10 3 m is 2.4 × 10 3 m . This means, tortuosity is not stationary at this distance. Furthermore, the average number of transitions for a particle to cross the domain is 2.4 × 10 3 / l c ≈ 16 . The first passage time t a ( x 1 ) derived from the TDRW approach is t a ( x 1 ) = ∑ i = 0 ⌈ ℓ ( x 1 ) / l c ⌉ - 1 τ i , τ i = l c v i , (4) where ℓ ( x 1 ) = χ ( x 1 ) x 1 is the average trajectory length. The ceiling function ⌈ ℓ ( x 1 ) / l c ⌉ = n v ( x 1 ) counts the number of different velocities experienced by the particle until first passage at x 1 . The first time increment τ 0 is distributed according to ψ 0 ( t ) = l c p 0 ( l c / t ) / t 2 , where p 0 ( v ) is the velocity PDF at the inlet plane. The random values τ i for i ≥ 1 are identical and independently distributed according to ψ s ( t ) = l c p s ( l c / t ) / t 2 , where p s ( v ) is the steady state velocity distribution. Figure 2 compares the predictions of the model with the DFN simulations. The stochastic approach predicts both the tailing behavior and peak times, while the early arrivals at short distances between the start and target plane are not captured. In fact, the first passage time PDF shows relatively abrupt early arrivals at short distances. This behavior can be traced back to the fact that, at short distances, the early arrivals are due to only a few long and fast fractures that connect start to target planes. In this sense, early arrivals are not ergodic as particles only have access to a very limited part of the velocity spectrum, while the TDRW approach samples velocities at each step from the full velocity distribution. The long time behavior on the other hand is due to particle motion along more tortuous paths that consist of larger numbers of slow and short fractures. The particles that contribute to the tails, in this sense form an ergodic subset, which is captured by the TDRW model. With increasing distance from the start plane, x c ≫ l c , more fractures contribute also to early arrivals, which then become ergodic. Thus, particle motion and first passage times are well represented by the stochastic particle motion. The construction of this stochastic model and its validation against the DFN simulations provides the starting point to develop a theory of long term properties of first passage times in random fracture networks. The model states first passage times as the sums of independent random time increments, which is equivalent to construction of the first passage time distribution f ( t ; x 1 ) as the n v -fold convolution of the transition time distributions ψ 0 ( t ) and ψ s ( t ) . In Laplace space, this can be written as f * ( λ ; x 1 ) = f n v * ( λ ) = ψ 0 * ( λ ) ψ s * ( λ ) n v - 1 , (5) where the asterisk denotes Laplace transformed quantities and λ the Laplace variable. Also note that the characteristic tails of the velocity distributions (cf. small velocities in Fig. 4 ) show that the transition time distributions ψ 0 ( t ) and ψ s ( t ) behave asymptotically as t - 1 - β . In fact, the Burr distribution (3) for the steady state velocity PDF p s ( v ) implies that ψ s ( t ) ∼ α τ c ( t / τ c ) - 1 - β , (6) for t ≫ τ c = l c / v c . Thus, for λ τ c ≪ 1 , the Laplace transforms of ψ 0 ( t ) and ψ s ( t ) can be expanded as ψ k * ( λ ) ≈ exp [ - λ ⟨ τ k ⟩ + a k ( λ τ c ) β + ⋯ ] , (7) where k = 0 , s , ⟨ τ k ⟩ is the mean transition time and a k a constant. The dots indicate contributions of order λ 2 . Using the approximation (7) for λ τ c ≪ 1 , we obtain f n * ( λ ) = exp ( - λ ⟨ t n ⟩ ) g n * ( λ ⟨ t n ⟩ 1 / β ) , (8) where [35] g n * ( λ ) = exp ( a 0 + ( n - 1 ) a s ⟨ t n ⟩ ( λ τ c ) β + ⋯ ) . (9) In the limit n → ∞ , the g n * ( λ ) converges toward g β ( λ ) = exp [ a s ( λ τ c ) β / ⟨ τ s ⟩ ] because higher order contributions go to zero at least as fast as n ( β - 2 ) / β [35] . This means, for n ≫ 1 , g n * ( λ ) in Eq. (8) can be replaced by g β * ( λ ) . By inverse Laplace transform of the resulting expression, we then obtain the following scaling form for f n ( t ) : f n ( t ) = g β [ ( t - ⟨ t n ⟩ ) / ⟨ t n ⟩ 1 / β ] ⟨ t n ⟩ 1 / β , (10) where g β ( y ) is a stable density that is totally skewed to the right [35,42,43] . This implies that ϕ ( y , x 1 ) = f [ t ( y ) , x 1 ] τ m 1 / β with y = ( t - τ m ) / τ m 1 / β and τ m = ⟨ t a ( x 1 ) ⟩ converges toward g β ( y ) for x 1 ≫ ℓ c . Figure 5 (top panel) shows the DFN data and TDRW prediction rescaled according to ϕ ( y , x 1 ) . The data in the tails collapse, which indicates the validity of the proposed reasoning. Note the TDRW approach identifies the first passage times as sums of independent random variables. Thus, the observed convergence to a stable density is the result of the generalized central limit theorem [43,48] because the distributions of time increments τ are heavy tailed. This is due to the behavior of the velocity distribution as p s ( v ) ∝ ( v / v c ) β - 1 with β < 2 for small velocities v < v c . For β > 2 , the first passage time distributions would converge to a stable distribution with a stability parameter 2, this means a Gaussian distribution. Thus, the theory identifies the features of the velocity point distribution that lead to the emergence of heavy-tailed first passage time distributions. Examples for networks with different exponents β that emphasize the robustness of these results are given in the Supplemental Material [35] . 5 10.1103/PhysRevLett.123.248501.f5 FIG. 5. (top panel) The predicted Lévy stable density (solid blue) and rescaled first passage time PDFs from the (symbols) DFN simulations and (lines) TDRW at distances (black) 25 r 0 , (green) 50 r 0 , and (orange) 1 0 2 r 0 . The dashed and dash-dotted blue lines show TDRW predictions at distances 1 0 3 r 0 and 1 0 4 r 0 from the inlet. The TDRW uses 1 0 7 particles. (bottom panel) Kullback-Leibler divergence between the rescaled first passage time distributions and the Lévy stable density plotted as a function of number of convolutions. The full circles correspond to the distributions shown in the top panel. The convergence towards the stable law depends on the number of independent transition times and thus velocities sampled along the path between start and target planes as illustrated in Fig. 5 . The convergence behavior is further investigated in the TDRW model, which is used to extrapolate the first passage time distributions to distances beyond the size of the study domain. The tails converge to the stable limit law for distances larger or equal than 1 0 4 m , which corresponds to around 200 convolutions, while the early time behavior converges at slower rate towards the Lévy stable distribution g β ( y ) . The rate of convergence is quantified in terms of the Kullback-Leibler divergence d k l ( n ) [49] shown in the bottom panel of Fig. 5 . It is fit by a decaying power-law d k l ( n ) ∼ n - 2 / 3 , which confirms convergence of f ( t ; x 1 ) toward the stable limit law. The presented analysis and theory identify the network properties that control first passage times in random three-dimensional fracture networks. Network topology and heterogeneity control the emergence of heavy-tailed first passage time distributions, which are constructed as convolutions of the distributions of transition times over hydrodynamically independent network features. The velocity point distribution controls the tailing behavior and thus the asymptotic stable density to which the first passage times converge. The convergence toward the universal stable behavior is controlled by the number of transitions between fractures along pathways and thus by the average fracture link length, and the tortuosity of trajectories, which are measures of the network topology. Despite the apparent complexity of three-dimensional random fracture networks, the first passage behavior and the emerging heavy tailed densities, can be understood and predicted in terms of only a few network and flow characteristics, which can be determined from transport-independent network attributes. The integration of these features in a TDRW modeling approach allows for the prediction of first passage time distributions at asymptotic and preasymptotic distances between start and target planes. J. D. H. thanks the Department of Energy (DOE) Basic Energy Sciences program (LANLE3W1) for support. J. D. H. and A. H. are thankful for support from the U.S. Department of Energy through the Los Alamos National Laboratory. Specifically, support through the Laboratory-Directed Research and Development Program Grants No. 20180621ECR and No. 20170103DR. Los Alamos National Laboratory is operated by Triad National Security, LLC, for the National Nuclear Security Administration of U.S. Department of Energy (Contract No. 89233218CNA000001). M. D. gratefully acknowledges funding from the European Research Council under the European Union’s Seventh Framework Programme (FP7/2007-2013)/ERC Grant Agreement No. 617511 (MHetScale). P. K. K. acknowledges a grant from the Korea Environment Industry & Technology Institute (KEITI) through Subsurface Environmental Management (SEM) Project, funded by the Korea Ministry of Environment (MOE) (2018002440003). [1] 1 O. Bénichou , C. Chevalier , J. Klafter , B. Meyer , and R. Voituriez , Nat. Chem. 2 , 472 ( 2010 ). NCAHBB 1755-4330 10.1038/nchem.622 [2] 2 Radionuclide Behaviour in the Natural Environment , edited by C. Poinssot and H. Geckeis ( Woodhead Publishing , Sawston, Cambridge, 2012 ). [3] 3 M. Peyrounette , Y. Davit , M. Quintard , and S. Lorthois , PLoS One 13 , e0189474 ( 2018 ). POLNCL 1932-6203 10.1371/journal.pone.0189474 [4] 4 O. Bénichou , C. Loverdo , M. Moreau , and R. Voituriez , Rev. Mod. Phys. 83 , 81 ( 2011 ). RMPHAT 0034-6861 10.1103/RevModPhys.83.81 [5] 5 S. Hwang , D.-S. Lee , and B. Kahng , Phys. Rev. Lett. 109 , 088701 ( 2012 ). PRLTAO 0031-9007 10.1103/PhysRevLett.109.088701 [6] 6 E. Bonnet , O. Bour , N. E. Odling , P. Davy , I. Main , P. Cowie , and B. Berkowitz , Rev. Geophys. 39 , 347 ( 2001 ). REGEEP 8755-1209 10.1029/1999RG000074 [7] 7 B. Berkowitz and H. Scher , Phys. Rev. Lett. 79 , 4038 ( 1997 ). PRLTAO 0031-9007 10.1103/PhysRevLett.79.4038 [8] 8 M. W. Becker and A. M. Shapiro , Water Resour. Res. 36 , 1677 ( 2000 ). WRERAQ 0043-1397 10.1029/2000WR900080 [9] 9 L. C. Meigs and R. L. Beauheim , Water Resour. Res. 37 , 1113 ( 2001 ). WRERAQ 0043-1397 10.1029/2000WR900335 [10] 10 R. Haggerty , S. W. Fleming , L. C. Meigs , and S. A. McKenna , Water Resour. Res. 37 , 1129 ( 2001 ). WRERAQ 0043-1397 10.1029/2000WR900334 [11] 11 S. Painter , V. Cvetkovic , and J.-O. Selroos , Geophys. Res. Lett. 29 , 20-1 ( 2002 ). GPRLAJ 0094-8276 10.1029/2002GL014960 [12] 12 S. Geiger , A. Cortis , and J. T. Birkholzer , Water Resour. Res. 46 , W12530 ( 2010 ). WRERAQ 0043-1397 [13] 13 P. K. Kang , T. Le Borgne , M. Dentz , O. Bour , and R. Juanes , Water Resour. Res. 51 , 940 ( 2015 ). WRERAQ 0043-1397 10.1002/2014WR015799 [14] 14 J. D. Hyman , M. Dentz , A. Hagberg , and P. Kang , J. Geophys. Res. 124 , 1185 ( 2019 ). JGREA2 0148-0227 10.1029/2018JB016553 [15] 15 B. Berkowitz , A. Cortis , M. Dentz , and H. Scher , Rev. Geophys. 44 , 1 ( 2006 ). REGEEP 8755-1209 10.1029/2005RG000178 [16] 16 R. Benke and S. Painter , Water Resour. Res. 39 , 1324 ( 2003 ). WRERAQ 0043-1397 [17] 17 S. Painter and V. Cvetkovic , Water Resour. Res. 41 , W02002 ( 2005 ). WRERAQ 0043-1397 10.1029/2004WR003682 [18] 18 B. Noetinger , D. Roubinet , A. Russian , T. Le Borgne , F. Delay , M. Dentz , J.-R. De Dreuzy , and P. Gouze , Transp. Porous Media 115 , 1 ( 2016 ). TPMEEI 0169-3913 [19] 19 P. K. Kang , M. Dentz , T. Le Borgne , and R. Juanes , Phys. Rev. Lett. 107 , 180602 ( 2011 ). PRLTAO 0031-9007 10.1103/PhysRevLett.107.180602 [20] 20 R. Haggerty and S. M. Gorelick , Water Resour. Res. 31 , 2383 ( 1995 ). WRERAQ 0043-1397 10.1029/95WR10583 [21] 21 J. Carrera , X. Sánchez-Vila , I. Benet , A. Medina , G. A. Galarza , and J. Guimerá , Hydrogeol. J. 6 , 178 ( 1998 ). 10.1007/s100400050143 [22] 22 O. Huseby , J.-F. Thovert , and P. Adler , Phys. Fluids 13 , 594 ( 2001 ). PHFLE6 1070-6631 10.1063/1.1345718 [23] 23 V. V. Mourzenko , J.-F. Thovert , and P. M. Adler , Phys. Rev. E 84 , 036307 ( 2011 ). PRESCM 1539-3755 10.1103/PhysRevE.84.036307 [24] 24 Y. Edery , S. Geiger , and B. Berkowitz , Water Resour. Res. 52 , 5634 ( 2016 ). WRERAQ 0043-1397 10.1002/2016WR018942 [25] 25 J. D. Hyman , S. L. Painter , H. Viswanathan , N. Makedonska , and S. Karra , Water Resour. Res. 51 , 7289 ( 2015 ). WRERAQ 0043-1397 10.1002/2015WR017151 [26] 26 P. K. Kang , M. Dentz , T. Le Borgne , S. Lee , and R. Juanes , Adv. Water Resour. 106 , 80 ( 2017 ). AWREDI 0309-1708 10.1016/j.advwatres.2017.03.024 [27] 27 J. D. Hyman , S. Karra , N. Makedonska , C. W. Gable , S. L. Painter , and H. S. Viswanathan , Comput. Geosci. 84 , 10 ( 2015 ). CGEODT 0098-3004 10.1016/j.cageo.2015.08.001 [28] 28 J. D. Hyman , C. W. Gable , S. L. Painter , and N. Makedonska , SIAM J. Sci. Comput. 36 , A1871 ( 2014 ). SJOCE3 1064-8275 10.1137/130942541 [29] 29 T. Walmann , A. Malthe-Sørenssen , J. Feder , T. Jøssang , P. Meakin , and H. H. Hardy , Phys. Rev. Lett. 77 , 5393 ( 1996 ). PRLTAO 0031-9007 10.1103/PhysRevLett.77.5393 [30] 30 Svensk Kärnbränslehantering AB, Data report for the safety assessment SR-Site (TR-10-52), Tech. Rep. Svensk Kärnbränslehantering AB, 2010. [31] 31 J. R. De Dreuzy , P. Davy , and O. Bour , Phys. Rev. E 62 , 5948 ( 2000 ). PLEEE8 1063-651X 10.1103/PhysRevE.62.5948 [32] 32 J.-F. Thovert , V. V. Mourzenko , and P. M. Adler , Phys. Rev. E 95 , 042112 ( 2017 ). PRESCM 2470-0045 10.1103/PhysRevE.95.042112 [33] 33 J. Bear , Dynamics of Fluids in Porous Media ( Dover Publications , New York, 1988 ). [34] 34 R. W. Zimmerman and G. S. Bodvarsson , Transport Porous Med. 23 , 1 ( 1996 ). [35] 35 See Supplemental Material at http://link.aps.org/supplemental/10.1103/PhysRevLett.123.248501 for details on DFN simulations, and the derivations for the scaling form of the stable distribution. [36] 36 LaGriT, Los Alamos Grid Toolbox, (LaGriT) Los Alamos National Laboratory, http://lagrit.lanl.gov ( 2013 ). [37] 37 P. Lichtner , G. Hammond , C. Lu , S. Karra , G. Bisht , B. Andre , R. Mills , and J. Kumar , Los Alamos National Laboratory Technical Report No. LA-UR-15-20403 , 2015 . [38] 38 N. Makedonska , S. L. Painter , Q. M. Bui , C. W. Gable , and S. Karra , Comput. Geosci. 19 , 1 ( 2015 ). CGEODT 0098-3004 [39] 39 S. L. Painter , C. W. Gable , and S. Kelkar , Computat. Geosci. 16 , 1125 ( 2012 ). 10.1007/s10596-012-9307-1 [40] 40 J. D. Hyman , G. Aldrich , H. Viswanathan , N. Makedonska , and S. Karra , Water Resour. Res. 52 , 6472 ( 2016 ). WRERAQ 0043-1397 10.1002/2016WR018806 [41] 41 R. Nordqvist , E. Gustafsson , P. Andersson , and P. Thur , Swedish Nuclear Fuel and Waste Management Co. Technical Report No. R-08-103 , 2008 . [42] 42 G. Samorodnitsky and M. S. Taqqu , Stable Non-Gaussian Random Processes: Stochastic Models with Infinite Variance ( Chapman & Hall Ltd, London , New York, 1994 ). [43] 43 J. P. Nolan , Stable Distributions–Models for Heavy Tailed Data ( Birkhauser , Boston, 2018 ), Chap. 1, http://fs2.american.edu/jpnolan/www/stable/stable.html . [44] 44 M. Dentz , P. K. Kang , A. Comolli , T. Le Borgne , and D. R. Lester , Phys. Rev. Fluids 1 , 074004 ( 2016 ). PLFHBR 2469-990X 10.1103/PhysRevFluids.1.074004 [45] 45 G. Srinivasan , J. D. Hyman , D. A. Osthus , B. A. Moore , D. O’Malley , S. Karra , E. Rougier , A. A. Hagberg , A. Hunter , and H. S. Viswanathan , Sci. Rep. 8 , 11665 ( 2018 ). SRCEC3 2045-2322 10.1038/s41598-018-30117-1 [46] 46 A. Comolli and M. Dentz , Eur. Phys. J. B 90 , 166 ( 2017 ). EPJBFY 1434-6028 10.1140/epjb/e2017-80370-6 [47] 47 A. Koponen , M. Kataja , and J. Timonen , Phys. Rev. E 54 , 406 ( 1996 ). PLEEE8 1063-651X 10.1103/PhysRevE.54.406 [48] 48 V. V. Uchaikin and V. M. Zolotarev , Chance and Stability: Stable Distributions and Their Applications ( Walter de Gruyter , London, 1999 ). [49] 49 S. Kullback and R. A. Leibler , Ann. Math. Stat. 22 , 79 ( 1951 ). AASTAD 0003-4851 10.1214/aoms/1177729694
Publisher Copyright:
© 2019 American Physical Society.