Numerical modelling of f l ow and transport in rough fractures

Scott Briggs,Bryan W.Karney,Brent E.Sleep

University of Toronto,Toronto M5S 1A4,Canada

Numerical modelling of f l ow and transport in rough fractures

Scott Briggs*,Bryan W.Karney,Brent E.Sleep

University of Toronto,Toronto M5S 1A4,Canada

A R T I C L E I N F O

Article history:

Received 30 June 2014

Received in revised form

5 October 2014

Accepted 8 October 2014

Available online 7 November 2014

Hydrogeology

Fracture f l ow

Solute transport

Computational f l uid dynamics

Lattice Boltzmann method(LBM)

Random walk(RW)

Simulation of f l ow and transport through rough walled rock fractures is investigated using the lattice Boltzmann method(LBM)and random walk(RW),respectively.The numerical implementation is developed and validated on general purpose graphic processing units(GPGPUs).Both the LBM and RW method are well suited to parallel implementation on GPGPUs because they require only next-neighbour communication and thus can reduce expenses.The LBM model is an order of magnitude faster on GPGPUs than published results for LBM simulations run on modern CPUs.The f l uid model is veri f i ed for parallel plate f l ow,backward facing step and single fracture f l ow;and the RW model is veri f i ed for pointsource diffusion,Taylor-Aris dispersion and breakthrough behaviour in a single fracture.Both algorithms place limitations on the discrete displacement of f l uid or particle transport per time step to minimise the numerical error that must be considered during implementation.

?2014 Institute of Rock and Soil Mechanics,Chinese Academy of Sciences.Production and hosting by Elsevier B.V.All rights reserved.

1.Introduction

As computer modelling has evolved,so too have two diametrically opposed requirements of those models in modelling. Certainly,it is desirable to improve a model’s computational speed and,yet,it is also desirable to improve the model’s accuracy and capability,and attributes that usually require an increase in computational time.Improvements can take the form of enhanced or additional physics or an increase in model resolution to improve accuracy.Often it is a combination of both factors,for example,a new physics model may require higher resolution,or f i ner meshes to simulate new phenomena.Therefore,one must balance model accuracy and computational time.

To address these computational limitations over the past few decades,more computational power was used,often simply taking advantage of new processor technology.However,single processors,or CPUs,are reaching a performance limit due to manufacturing constraints.Therefore,to continue improving performance,the CPU industry has moved toward using multiple CPUs in parallel.The challenge with this approach then becomes the implementation of conventional numerical algorithms and methods on parallel architectures,including clusters of CPUs and graphics processing units(GPUs).GPUs have evolved over time with more complex computing capabilities,similar to a conventional CPU,and are referred to as general purpose GPUs(GPGPUs).

The lattice Boltzmann method(LBM)is increasingly used for the simulation of fl uid fl ows in complex geometries(Stockman et al., 1998;Eker and Akin,2006;Yan and Koplik,2008;Dou et al., 2013).However,its engineering applications have been limited by the required computing power.The local nature of LBM,where only next-neighbour node communication is required,is suitable for parallelization.Previous work has shown that an increase of an order of magnitude in performance can be expected when implementing LBM on a GPGPU(Bailey et al.,2009;T?lke,2010).However,such work did not show the applicability and validation for fl ows in rock fractures that are of interest to hydrogeology.

Other computational f l uid dynamic(CFD)methods begin with the continuum Navier-Stokes equations governing the macroscopic movement of f l uids,and then discretize these equations with a suitable numerical method(Eker and Akin,2006).In the LBM model,the microscopic interaction of particles on a grid and the averaging of those interactions emerge into the macroscopic continuum of a f l uid.These interactions include two main steps: streaming and collision.The streaming step is a translation of particles from one node on the grid to the next.The collision step conserves momentum by redirecting particles that‘collide’or occupy the same node.

This study demonstrates a veri f i ed GPGPU code for simulating two-dimensional(2D)laminar f l ow through rock fractures using a D2Q9 LBM.

Effective understanding of solute transport in fractures is underpinned by the need for accuracy in the simulation of fl uid fl ow.To account for the effects of tortuosity(Tsang,1984)and Reynolds number above unity,a CFD approach is used.A CFDapproach,like the LBM,provides local velocities throughout the model domain which are used to move particles through the process of advection.A diffusive process is also included using a random walk(RW)algorithm which is shown to accurately capture the complex geometries associated with single fractures.

In this study,high performance GPGPU numerical methods are developed,validated and shown to be capable of modelling f l ow and transport in synthetic and real fractures.It is increasingly important for projects of all scales to conduct modelling studies which may require meshes with millions of nodes or large parametric searching,for example,variations in Reynolds number, boundary conditions and fracture geometries.The large meshes,or grids,become computationally and f i nancially expensive on conventional CPUs or clusters.However,a single GPGPU can bridge the gap to high performance computing if the required algorithms are well implemented on the GPGPU with suf f i cient performance advantages.

2.Model implementation

2.1.Lattice Boltzmann method(LBM)

The LBMs have been used in a variety of engineering applications,including the f i eld of porous and fractured media f l ow(Sukop et al.,2013).Additional development of LBM theory can be found in the literature(Succi,2001;Sukop and Thorne,2006;Latt,2007).For the purpose of modelling f l ow in fractures,a 2D LBM code was developed using nine velocity directions ei,also known as D2Q9. The LBM can be summarised in the following form:

where the left hand side of the equation represents the streaming step and the right hand side represents the collision step,andτis the relaxation parameter that governs the rate at which the f l uid tends towards equilibrium.For the LBM model presented,τtakes the following form:

whereυLis the numerical viscosity de f i ned by the discretization of the system into lattice units.

The model runs on a GPGPU using CUDA,a proprietary programming model developed by NVIDIA.One of the drawbacks of GPGPU implementations is the discrepancy between 32-bit and 64-bit f l oating point precision as current hardwares have limited support for 64-bit,or double precision calculations.Typical GPGPU implementations offer double precision performance that is approximately one third or one quarter that of single precision performance,depending on the model and manufacturer.Without double precision calculations,the numerical error,or the code complexity required to compensate for error,increases.

Another type of numerical error in CFD models,conventionally referred to as numerical dissipation,describes the arti f i cial dissipation of momentum in f l uid.Since the LBM is essentially a f i nite difference approximation to the Boltzmann equation,it is subjected to the same numerical truncations as other f i nite difference methods.The numerical error can cause dissipation of the advection term which by de f i nition should be free of dissipation(Zhu et al.,2006).

To minimise the potential for numerical instabilities in the LBM and maintain the second order accuracy of the LBM,the model parameters are de f i ned using the method presented by Latt and Krause (2008)as part of the OpenLB User Guide.The process involves selecting physical units then converting to lattice units to f i nally obtain the relaxation parameterτ.The relaxation parameter plays an important role in the collision term of the LBM.It controls the tendency of the system to move towards local equilibrium.In the literature,the relaxation parameter has been found to cause numerical instabilities when it approaches 0.5(τmust be greater than 0.5 for physical viscosities).Stable values ofτclose to unity are preferred for numerical accuracy of the LBM and can be found using the method outlined below(Sukop and Thorne,2006;Sukop et al.,2013).

In this study,water is the physical f l uid being simulated with a kinematic viscosity,υ,in a fracture of aperture,2a,and with physical velocity,u.This leads to an expression for the Reynolds number:

The dimensionless expression for Reynolds number is then used to convert from the physical units of the system to lattice units.The fracture width is discretized into lattice nodes of lengthδxwith discrete timeδt.In order to minimise the slightly compressible nature of the LBM and maintain the second order accuracy,the following constraints are used respectively when determining system discretization:

Practically,to ensure stability and numerical accuracy,these constraints are addressed by limiting the numerical velocity to a maximum of 0.1 lattice units per time step,which minimises the partial compressibility of the LBM(Sukop and Thorne,2006).The lattice viscosity(υL)is calculated based on the discretization of the system and the dimensionless Reynolds number.Finally,the relaxation parameter is calculated according to Eq.(2)and is keptas close to unity as possible by adjusting the mesh size and maximum lattice velocity.

2.2.Boundary conditions

One of the distinct advantages of the LBM comes from its discrete nature.It is ef fi cient for modelling complex geometries (Chen et al.,1994;Eker and Akin,2006;Lammers et al.,2006; Brewster,2007)that arise in the analysis of rock fractures.Within the modelling domain,each node may represent a rock mass or fl uid node.At the solid boundaries,a no-slip condition is used to create a zero velocity boundary along the surface.A different set of collision equations are used at the solid boundary and are referred to as mid-plane bounce-back boundary conditions(Succi,2001). The name arises from the applied boundary rules where particles entering a boundary at time t are sent back out with equal velocity magnitude and opposite direction at time t+Δt.This effectively puts the boundary at a distance midway between a fl uid and solid node.

Constant fl ux,pressure and gravity-driven boundary conditions can be used to drive the fl uid through the fracture.Solid and no-slip boundaries are used along the fracture surfaces while periodic boundary conditions are used atentry and exitof the fracture where fl uid and solutes leaving the fracture are re-injected with equivalent velocity and direction.Periodic boundary conditions simulate an in fi nite domain with periodically repeated geometry.A periodic,or wrapped boundary,in combination with applied gravity boundary conditions,removes entry or exit effects which would otherwise arise under conventional constant fl ux boundaries.

Gravity-driven f l ow is used to minimise computational complexity and entry effects as mentioned above.The force of gravity is converted to a velocity and added to the velocity component parallel to the primary fracture axis.Gravity-driven conditions are implemented according to the method described by Sukop and Thorne(2006).Gravity-driven f l ow acts on each cell of the LBM independently,therefore it is unnecessary to use conventional pressure or velocity boundary conditions as this would only add an arti f i cial constraint into the system,possibly creating entry or exit effects.The acceleration due to gravity is converted to a velocity term from the following relationship:

where F is the external force added into the LBM calculations in the form of a local velocity.The mass(m)is proportional to the density (ρ)and the relaxation parameter(τ)can be substituted for time arriving at:

whereΔu represents a discrete velocity increment that is added to the velocity component parallel to the fracture plane used to calculate the equilibrium distribution function.

2.3.Solute transport

Solute transport is simulated by modelling the discrete movement of particles through the fracture.Advective processes are known from the LBM simulation and the local velocity information is used to displace particles over each time step.Diffusion follows a RW process to model discrete particle movement.The RW group of methods has been developed and used extensively for the purpose of solute transport in porous and fractured media(Ahlstrom et al., 1977;Tompson and Gelhar,1990;Wels et al.,1997;James and Chrysikopoulos,1999;Delay et al.,2005;Nowamooz et al.,2013).

For the purposes of studying solute transport and the effects of roughness in a single fracture,particles are assumed to be neutrally buoyant and exhibit no decay or matrix diffusion.Particle-particle interactions are not modelled,nor do particles affect the f l ow solution.The only forces acting on the particles are due to advection and diffusion.Thus,transport is described by the following Fokker-Planck equation(James and Chrysikopoulos,2011): where Dmis the molecular diffusion coef f i cient,uiis the local velocity at the location x of the particle at time t,Z(0,1)is a normally distributed random number for each dimension i with mean zero and a standard deviation of unity.For a more detailed development of the Fokker-Planck approach to RW,the reader is referred to Delay et al.(2005).

Numerical error is minimised by ensuring a particle moves a maximum of one halfΔx per time step(Maier et al.,2000).The left hand side of Eq.(9)represents the maximum advection,due to umax,and the maximum likely diffusion distance during any given time step:

With a suf f i cient number of particles,a RW method can accurately model the process of Brownian motion used to model diffusion.In general,more than 100,000 particles are required to suf f i ciently model diffusion using RW(Hassan and Mohamed, 2003).Such a large number of particles,running over many iterations,requires a high quality pseudo random number generator (PRNG).It is important that the period of the PRNG,among other statistical properties,is of suf f i cient capability to ensure that the random numbers generated do not affect the model results.While robust PRNGs have been available for some time(Matsumoto and Nishimura,1998),parallel implementations are an area of ongoing research.Only recently have developments by Saito and Matsumoto(2013)signi f i cantly improved the quality of random number generation on GPGPU hardware.

The particles are injected at the inlet to the fracture as an instantaneous injection.Given the synthetic fracture is only 100 mm long and the length constraint obtained from Eq.(9),the particles must wrap around the system until the minimum length is reached to ensure that dispersion is allowed to fully develop. The resident times,t,for all particles are tracked and plotted as a histogram representing concentration against time.To generalise the presented data,resident times are non-dimensionalized using the relative fracture properties and are expressed by the pore volume(PV):

where Q is the f l ow rate in 2D f l ow through the fracture calculated from the LBM velocity data,and L is the total length of the fracture in which the particle travels.

2.4.Fracture pro f i les

Flow through a single rock fracture is modelled for two cases and a third fracture set is used to explore solute transport.The f i rst consists of a one-sided fracture aperture collected by Boutt et al. (2006)and the second consists of a 2D slice through a dolomite fracture generated in the laboratory(Mondal and Sleep,2012). Finally,synthetic fractures are generated and used for solute transport.

The dolomite fracture was created from a block approximately 350 mm long,250 mm wide and 70 mm thick.The rock sample contained stylolites,which are the planes of weakness,parallel to the length of the block.A fracture was introduced in the rock block using the method described in Reitsma and Kueper(1994)resulting in f i nal dimensions of 280 mm?210 mm?70 mm.A threedimensional(3D)stereo-topometric measurement system,the advanced topometric sensor(ATOS)II manufactured by GOM mbH, was used to determine the surface pro f i le of the fracture walls and its aperture distribution.For more details on the preparation of the sample and the ATOS II system,see Mondal and Sleep(2012)and Tatone and Grasselli(2009),respectively.A 16 mm 2D slice through the 3D surface created by the ATOS II was used in the LBM.Using a 2D approximation of the fracture to represent the 3D surface saves signi f i cant computational resources.However,a 2D system cannot capture or quantify the effects of contact points in a fracture and the impact of reducing effective apertures and increasing tortuosity (Zimmerman and Bodvarsson,1996).Tortuosity in fractures refers to the circuitous path that a f l uid particle will travel due to the small and large scale roughness of a rock fracture.The 2D modelling is able to capture the effects of 2D tortuosity and can be an effective means of providing insight into the hydraulic behaviour of rough fractures.

For the purpose of transport in fractures,synthetic fractures are used.Synthetic fracture generation creates systems with controlled surface properties.To capture the variation of the crucial roughness of a fracture,a series of similar fractures with increasing roughnesswas generated.Fracture surfaces are generated with the software package SynFrac developed by Ogilvie et al.(2006)who expanded based on earlier work(Brown et al.,1995;Glover et al.,1998a,b)to capture the complex nature of natural fractures with synthetic approximations.Glover and Hayashi(1997)demonstrated that modelling a synthetic fracture at the centimetre scale applied directly to f i eld f l ow measurements at the 100 m scale.An important consideration when generating synthetic fractures is capturing the fracture properties at all wavelengths.The top and bottom of a single fracture willhave correlated geometry and surface properties at long wavelengths but are mostly independent at short wavelengths.The threshold separating long and short wavelength is called the‘mismatch length’.SynFrac has multiple methods for determining the mismatch length.In this study,the SynFrac implementation of the Brown et al.(1995)mismatch length was set to 15 mm.

Using SynFrac for 2D fracture generation,two studies were conducted.First,a series of fractures were generated with increasing roughness by increasing the value of the fractal dimension input parameter.Second,to analyse random variations that may occur in the fracture surfaces generated by SynFrac,multiple fractures with identical characteristic parameters were created by varying the seeds of the Park and Miller pseudo-random number generator (SRNG)in SynFrac.A 100 mm 2D pro f i le is manually extracted from the data and selected such that it has no contact point.Since each 3D fracture generated with SynFrac is adjusted so that the relative separation of the top and bottom surfaces creates a single contact point,an equivalent adjustment was needed in 2D for consistency between fracture studies.Therefore,the mean aperture of each fracture was kept constant for each study by manually adjusting the pro f i le separation.Other SynFrac settings are kept constant including the resolution(1024?1024),standard deviation(1 mm) and anisotropy factor(1.0).The fracture length of 100 mm and a 1024 element resolution is expanded to a 2048 element resolution using interpolation,resulting in a 48.8μm element size.

Fractures were generated by specifying a fractal dimension for the 3D surface in SynFrac ranging from 2.00 to 2.35.It is recognised that using a fractal dimension for de f i ning roughness is incomplete as fractal dimensions are not unique to an object;two similar but unique objects may have the same fractal dimension.It has also been shown that the direction of f l ow in a fracture yields varying results(Boutt et al.,2006)whereas the fractal dimension of a surface is independent of the direction of measure.These fractures with varying roughness were used and compared with an equivalent parallel plate system and a real dolomite fracture(Fig.1). Ogilvie et al.(2006),developers of SynFrac,used complementary software,ParaFrac,to analyse real rock fractures and found that sandstone and granodiorite samples had fractal dimensions approaching 2.35.

3.Validation

3.1.Flow between parallel plates

The cubic law conventionally used to describe f l ow through rock fractures assumes the fractures can be modelled as parallel plates and is used for comparison with the LBM.A plate separation,2a,is de f i ned by cubic law calculations for f l ow and conventionally a fracture aperture is approximated using the arithmetic mean of the aperture.The f l ow through parallel plates as described by the cubic law is

where W is the width of the fracture and in all simulations is taken as unity for the 2D models,L is the length of the fracture,andΔh represents the change in head over the length of the fracture.In the case of gravity-driven f l ow where gravity acts along the length of the fracture,Δh=L,and Eq.(11)becomes

Fig.1.Fracture pro f i les with varying roughness through dolomite fracture generated using a synthetic fracture generator called SynFrac.Total fracture length is 100 mm and each fracture has a mean aperture of 1.7 mm,only the fractal dimension(FD)variable is changed in SynFrac.Fracture pro f i le(a)represents a parallel plate system with an equivalent 1.7 mm aperture.Fracture pro f i le(j)represents a 16 mm long strip from a dolomite fracture with mean aperture 0.1 mm.

The f l ow rates calculated by the cubic law are compared to those calculated by the LBM model.

LBM model results of f l ow between parallel plates using incompressible f l uids have been compared with available analytical solutions.For laminar f l ow conditions,the Hagen-Poiseuille equation can be used to describe the horizontal velocity through a cross-section:

where x is the distance from the centreline and G is the driving force.This analytical solution yields a parabolic velocity pro f i le.For the case of gravity-driven f l ow,G=ρg,which de f i nes a local pressure gradient throughout the system.The maximumvelocity occurs at the centreline where x=0 and the average velocity is 2/3 of the maximum velocity.Substituting for these changes gives:

where uavgis the average velocity.

Using the dimensionless Reynolds expression,physical parameters can be converted to equivalent lattice parameters.Lattice spacing is determined by the geometry and discretization of the physical system.Lattice velocity is limited to a maximum 0.1 lattice unit per time step which arises from the approximations used in the LBM(Sukop and Thorne,2006).Viscosity is determined by constraining the relaxation parameter to unity(τ=1)to ensure numerical stability.

The force of gravity in Eq.(14)can be altered and is used to drive fl ow in the model.For the case of parallel plates,when the numerical model reaches steady state,it compares well to the analytical solution as shown in Fig.2.Fig.2 shows horizontal velocity pro fi les versus the ratio of velocity(u)to the maximum velocity(Umax)at all nodes across the model domain.The parallel plate boundaries are con fi gured for a 256-node spacing.The only locations where the Poiseuille pro fi le and model pro fi le differ are at the two closest nodes to the boundary and this can be attributed to the implementation of the bounce-back rules and is typical of the LBM approach(Sukop and Thorne,2006).

3.2.Backward facing step

Modelling parallel plates is a crucial step in the validation of an algorithm as it allows for comparison with analytical solutions for the same problem.However,fracture geometries introduce complex geometries with changing apertures,including constrictions and openings.To simulate the effects of an abrupt aperture opening,a simple backward facing step is modelled and compared with experimental results.

The geometry used to simulate a backward facing step consists of parallel plates with an upstream step that is half the plate separation and a length at least 20 times the height to minimise the interference from out f l ow boundary conditions.Flow is gravitydriven and the boundary conditions are periodic.For the upstream of the step,the f l ow f i eld is given suf f i cient distance to develop a parabolic velocity pro f i le and similarly,for downstream of the step,the f l ow becomes parabolic.Finally,the f l ow is gradually ramped back to the step height for the periodic boundary.All other boundaries use standard bounce-back rules.

The geometry backward facing step is well studied but conventionally at higher Reynolds number than that necessary for the study of f l ow in fractures.Results presented by Armaly et al. (1983)for Re?200 are used with Re=100 being the lowest of the available data.The relationship of the dimensionless Reynolds number(Re=2au/υ)de f i nes the f l ow and the upstream height or step height is used.

Generally,f l ow over a backward facing step at Re?200 consists of three segments:f l ow separation directly after the step consisting of an area of recirculation,followed by a bulk f l ow reattachment point,and f i nally development of a parabolic velocity pro f i le downstream of the step.The reattachment point refers to the end of the recirculation zone.Fig.3 illustrates three cases,Re=100, Re=150 and Re=200 where f l ow is from left to right and relative colouration from blue to red represents slow to fast f l uid velocities. Reynolds number is varied by adjusting lattice parameters,specifically the force of gravity for the system.Using graphical means of measurement,the reattachment points of the model results are in good agreement with those reported by Armaly et al.(1983)of 3,4 and 5 times the step height for the Reynolds number of 100,150 and 200,respectively.In each case the velocity pro f i le becomes parabolic again far downstream.

Fig.2.Horizontal velocity pro f i les comparing the analytical results of a Poiseuille pro f i le and the model result.

3.3.Flow in a single fracture

It is not trivial to validate f l ow in real or synthetic fractures because it would require comparison against other CFD implementations.Many CFD codes are not freely available or those codes that are free,or open source,thus require their own implementation,compilation and validation.Alternatively,one can compare the model results with those from the literature, including experimental results.Once such comparison arises from a method developed by Zimmerman et al.(1991)and Renshaw(1995),which uses a statistical measure of fracture roughness.

Fig.3.Flow,from left to right,over a backward facing step(leftmost edge).Shown as red vertical lines,the reattachment lengths are approximately 3,4,5 times of step heights for Re=100,150 and 200,respectively.The step height is half the downstream width.Velocity is plotted with red representing the fastest velocities and blue the slowest.The velocity pro f i le is parabolic immediately upstream and far downstream of the step while the zones outside of this region are omitted for clarity.

Since both equations(Eqs.(15)and(16))depend only on the variance of the lognormal aperture distribution,they can be combined and are shown in Fig.4.The ratio of hydraulic aperture to mechanical aperture tends towards unity as either the fracture aperture increases or the standard deviation decreases as the walls become smoother.Experimental data by Zimmerman and Main (2003),not shown on the graph,f i t well with the theoretical data. Other numerical work by Patir and Cheng(1979)and Brown(1987) also compared similarly with the theoretical data.The model predictions plotted in Fig.4 were calculated using f l ow results from a single fracture while increasing the mechanical aperture(dm)and shows that,as the separation between the rock fracture walls increases,the lognormal variance also changes.

Comparison to cubic law f l ow rates is another form of validation, speci f i cally when compared to results in the literature.Flow through a single rock fracture is modelled for two cases to demonstrate the utility of the LBM.The f i rst is a one-sided fracture aperture collected by Bouttet al.(2006)and the second consists of a 2D slice through a fracture generated in the laboratory.

For the f i rst fracture data set,the cubic law deviates 8.4%from the f l ows rates determined by the LBM model at Re=6.Flow rates are chosen based on the discretization of the fracture and the need to maintain a relaxation parameter close to unity.The fracture has an equivalent aperture of 359μm which is translated into a fracture velocity of 1.67?102m/s for an equivalent parallel plate system.The same 8.4%deviation from the cubic law is found at Reynolds number of 0.06,0.6,6 and 60.These results are in line with observations reported by other researchers such as Brush and Thomson(2003)who found the cubic law to be within 10%of their Stokes Law simulations for Reynolds number less than unity.

The second data set analysed in our study has an equivalent aperture of 100μm and at Re=6,the velocity through an equivalent parallel plate system is 5.97?102m/s.In this case the deviation from the cubic law at Re=6 is 50%and the same approximate deviation holds for Re from 0.06 through 6.Again, variations in literature can be found.For example,Brown(1987) who used the Reynolds equation to describe f l ow between slightly rough non-planar surfaces,found the cubic law to hold within a factor of 2,while Tsang(1984)showed an order of magnitude variation from the cubic law if tortuosity was ignored.

The two fracture data sets show different deviations from the cubic law that could be due to their different physical attributes. The f i rst data set consists of the bottom pro f i le of a fracture and the top pro f i le is a smooth plate.This geometry would more closely approximate that of parallel plates then the second data set,where both sides are represented by a fracture pro f i le.The equivalent aperture is also much larger in the f i rst data set,359μm versus 100μm in the second data set,potentially affecting hydraulic properties and deviations from the cubic law.

The two data sets highlight the dif fi culty in using the cubic law for fracture fl ow.That is, fl ow in fractures cannot always be simulated by the cubic law.The simulations show the advantages of using the LBM model which can be discretized to account for roughness intrinsically.The main advantage of LBM models compared to other approaches is that they can be used to resolve small scale effects such as an abrupt change in aperture where there is a signi fi cant change in velocity streamlines and potential secondary fl ows.Fig.5 shows the fl ow streamlines calculated at two different locations along the same fracture at three different Reynolds numbers.The fi rst location on the left is an area of a large change in aperture while the second location is of a small depression in the fracture.Even at low Reynolds number,Re=0.6,the fl ow has zones of recirculation,creating areas of the fl uid that do not actively contribute to bulk fl ow.As the Reynolds number is increased(6 then 60),the recirculation zones become larger and appear in more places.

Fig.4.The ratio of hydraulic aperture to mechanical aperture is plotted against statistical roughness of the fracture as described by Renshaw(1995).The model f i ts well with theoretical data.The model predictions are plotted from a single fracture by increasing the mechanical aperture dm.

These results demonstrate how the roughness of a fracture can affect fl uid fl ow within a fracture even at low Reynolds numbers and provide some hints to explain the discrepancy between fl ow rates expected from the cubic law and results from the LBM.Fig.6 compares the results from fl ow in the fi rst data set compared to fl ow between parallel plates with an equivalent mechanicalaperture.The left hand side of Fig.6 consists of a rock fracture along the base of the model with a no-slip smooth top boundary,constant gradient outlet and parabolic inlet boundaries.The right hand side of Fig.6 models f l ow through parallel plates spaced at an equivalent aperture calculated using the arithmetic mean.It can be seen that the actual rock fracture compresses the velocity pro f i le much more than that of the equivalent fracture.It is the peaks of the rock fracture that signi f i cantly change the velocity distribution,leading to an apparently smaller equivalent aperture.The f l ow distribution is clearly different fromthat predicted by simple parallel plates,and although it cannot be seen in Fig.6,there are areas of recirculation downstream of each fracture constriction(see Fig.5).Since this is a complex phenomenon,it would be dif f i cult to create a single variable that could be adjusted for such effects.Rather,it is important that a given system be simulated with a model such as the presented LBM model.

Fig.5.The streamlines are plotted as the Reynolds number increases from 0.6 to 60.Secondary f l ows develop in the form of eddies and grow to f i ll a larger cross-section of the aperture.Each node represents approximately 2μm.

3.4.Solute transport

For a point-source in 2D space,the analytical solution for diffusion as developed by Crank(1975)in the form shown by Sukop and Thorne(2006)is

Fig.6.Left hand side:Flow through a fracture.Right hand side:Flow through parallel plates with the mechanical aperture equivalent to the fracture aperture on the left.Relative velocity is plotted with yellow representing the fastest velocities and dark blue the slowest.

where C is the concentration,Cois the initial concentration,Dmis the molecular diffusion coef f i cient,Mois the initial mass,t is thetime,and r is the spatial coordinate in any given linear direction. Fig.7 shows the results from the RW algorithm at three different time increments and the corresponding analytical solutions from Eq.(17).The f i t between the RW and the analytical solution is excellent with some variation from the analytical solution due to the random nature of the RW method.

Taylor-Aris dispersion between parallel plates is de f i ned as follows(Stockman,1997):

where Deffis the effective dispersion coef f i cient.The above equation holds for a range of Peclet numbers(Pe)of

where the Peclet number is de f i ned as

Fig.8 indicates a good f i t between the discrete RW and the analytical dispersion equation for parallel plate applications.Values chosen are similar to those found in Sukop and Thorne(2006)as an additional measure of comparison and validation.

Finally,the model for solute transport is validated by measuring resident time in a single fracture to capture the breakthrough curve (BC).BC data from the model results are compared with the analytical solution for a parallel plate system with an equivalent mechanical(arithmetic mean)aperture using Eq.(21)and labelled as the analytical solution in Fig.8.The analytical solution for the concentration of a solute subjected to uniform f l ow u at a location x at time t for a one-dimensional(1D)instantaneous injection is (Hunt,1978):

where M is the initial concentration of particles in the system,σis the porosity and taken at unity.Fig.9 shows the BCs for Dm=7? 1010m2/s at Re=1 for the synthetic fractures generated in Syn-Frac with fractal dimensions from 2.00 to 2.35.Also plotted is a parallel plate system with the same mechanical aperture of the rougher fractures modelled using the same LBM and RW method (labelled as a slit).From Fig.9,the modelled slit results are very close to the expected analytical solutions,while non-Fickian behaviour was apparent if the inequality of Eq.(19)was not maintained as expected(Cardenas et al.,2009;Qian et al.,2011). Results from fractures of varying roughness show some deviation from the analytical solution due to the inherent geometry of the system that is captured by the discrete approach of the LBM and RW method.

Fig.7.Point source diffusion in 2D and the relative concentrations at a given radius from the source.Analytical results for time step ts=1000,2000 and 10,000 are shown respectively.

Fig.8.Effective dispersion for the values:uavg=0.0038 lu/ts and Dm=0.0013（lu）2/ts after Sukop and Thorne(2006).The input values are given in terms of lattice unit(lu)and time step(ts),typical for LBM applications.

3.5.Transport sensitivity analysis

The susceptibility of the transport code to numerical error is estimated by Eq.(14).The constraint de f i nes the maximum diffusion and advection distance that a particle may move during a single time step.To maintain the inequality,the discrete time step can be reduced while the discretization of space is f i xed for a given model.For a given Reynolds number,the maximum time step will change.The example of Re=50 is shown in Fig.10.A total of four models are run,two above and two below the empirical limit expressed by Eq.(19).

Fig.9.Breakthrough curves for Dm=7?1010m2/s at Re=1 for synthetic fractures generated from a 2D slice of a 3D surface with fractal dimensions from 2.00 to 2.35. The‘Slit’represents a parallel plate system modelled in the same way;f i nally the analytical solution is shown for comparison.

Fig.10.Data shown are for Re=50 for synthetic fractures generated from a 2D slice of a 3D surface with fractal dimensions(FD)from 2.00 to 2.35.Cases 1 and 2 do not meet the constraint for minimising numerical error while Cases 3 and 4 do satisfy the constraint.

Case 1 represents the largest discrete time while Case 4 represents the smallest time discretization,absolute values range over two orders of magnitude.

The sensitivity analysis illustrated in Fig.10 shows that considerations of the minimum time step must be taken into account to reduce numerical error that is most prevalent in Case 1.Cases 3 and 4 begin to reach convergence and for the purpose of accuracy balanced with computation limits the system discretization of Case 3 would be suf f i cient for most engineering applications.In terms of particle count,215particles are suf f i cient for most models (Fig.11).While more particles allow a better averaging of the RW method,the resident times converge to a solution for particle counts around 215.

4.Model performance

Performance of the presented LBM on the GPU is approximately an order of magnitude faster than a comparable LBM model run on a CPU,which is consistent with the f i ndings of T?lke(2010).Typically,performance in LBM codes is measured in million lattice updates per second(MLUPS).Single CPU codes typically perform around 6.2 MLUPS(Bailey et al.,2009)and more recently 88 MLUPS (Habich et al.,2013).The GPU model in this study,the NVIDIA GeForce GTX Titan,achieves over 1000 MLUPS for a grid size of 2048 by 128 nodes using double precision calculations.The comparison should be taken as a rough estimate as this is not intended to compare directly between models which would require the equivalence of grid size,LBM implementation,optimisations and other factors affecting the performance of computer code.

5.Conclusions

The LBM model is well suited for simulating laminar fl uids through single rock fractures where complex fl ow patterns are produced by the irregular bulk rock surfaces.Even under laminar fl ow conditions,tortuous fl ow paths and surface roughness result in unique fl ow conditions that the LBM model can effectively capture.The LBM model presented in this study agrees with results of modelling of fl ow in rock fractures(Tsang,1984;Brown,1987)and is also consistent with the statistical roughness model described by Zimmerman et al.(1991)and Renshaw(1995).

The RW model is a simple implementation of advection and diffusion,however,it is able to successfully capture particle transport in systems with varying geometry.The RW model coupledwith the LBM matches well with the analytical results of f l ow and transport between parallel plates and single fracture geometries with varying synthetic properties,for example using fractal dimensions to distinguish varying surface roughnesses.

Both the LBM and the RW are well suited to implementation on the parallel computing architecture of GPGPU hardware,as the local nature of their respective algorithms reduces communication overhead between processing nodes.Additionally,with the development of new parallel PRNG,RW can be effectively implemented on GPGPU hardware with,statistically,high quality random numbers.The LBM model developed allows for the effective simulation of fracture f l ow and is capable of simulating 2D systems at the micron scale over a 100 mm global domain and can be used to simulate systems an order of magnitude faster than some CPU-based codes,allowing for signi f i cantly faster analysis of single fracture f l ow and transport.However,both LBM and RW are susceptible to numerical error unless proper constraints are maintained.

Fig.11.For a set of bin size when calculating the histogram,a larger number of particles give a more accurate description of the dispersion of particles through the fracture without changing the overall behaviour.Data shown are for Re=50 for synthetic fractures generated from a 2D slice of a 3D surface with fractal dimensions(FD)from 2.00 to 2.35.

Con f l ict of interest

The authors wish to con f i rm that there are no known con f l icts of interest associated with this publication and there has been no signi f i cant f i nancial support for this work that could have in f l uenced its outcome.

Ahlstrom SW,Foote HP,Arnett RC,Cole CR,Serne RJ.Multicomponent mass transport model:theory and numerical implementation(discrete-parcelrandom-walk version).Tech.Rep.BNWL-2127.Richland,USA:Battelle Paci fi c Northwest Laboratories;1977.

Armaly BF,Durst F,Pereira JCF,Sch?nung B.Experimental and theoretical investigation of backward-facing step fl ow.Journal of Fluid Mechanics 1983;127: 473-96.

Bailey P,Myre J,Walsh SDC,Lilja DJ,Saar MO.Accelerating lattice Boltzmann fl uid fl ow simulations using graphics processors.In:Proceedings of the 2009 International Conference on Parallel Processing(ICPP 2009).Piscataway,USA:IEEE; 2009.p.550-7.

Boutt DF,Grasselli G,Fredrich JT,Cook BK,Williams JR.Trapping zones:the effect of fracture roughness on the directional anisotropy of fl uid fl ow and colloid transport in a single fracture.Geophysical Research Letters 2006;33(21).http:// dx.doi.org/10.1029/2006GL027275.

Brewster JD.Lattice-Boltzmann simulations of three-dimensional fl uid fl ow on a desktop computer.Analytical Chemistry 2007;79(7):2965-71.

Brown SR.Fluid fl ow through rock joints:the effect of surface roughness.Journal of Geophysical Research 1987;92(B2):1337-47.

Brown SR,Stockman HW,Reeves SJ.Applicability of the Reynolds equation for modeling fl uid fl ow between rough surfaces.Geophysical Research Letters 1995;22(18):2537-40.

Brush DJ,Thomson NR.Fluid f l ow in synthetic rough-walled fractures:Navier-Stokes,stokes,and local cubic law simulations.Water Resources Research 2003;39(4).http://dx.doi.org/10.1029/2002WR001346.

Cardenas MB,Slottke DT,Ketcham RA,Sharp JM.Effects of inertia and directionality on f l ow and transport in a rough asymmetric fracture.Journal of Geophysical Research:Solid Earth 2009;114(B6).http://dx.doi.org/10.1029/2009JB006336.

Chen S,Doolen GD,Eggert KG.Lattice-Boltzmann f l uid dynamics:a versatile tool for multiphase and other complicated f l ows.Los Alamos Science 1994;22:98-109.

Crank J.The mathematics of diffusion.2nd ed.Oxford:Clarendon Press;1975.

Delay F,Ackerer P,Danquigny C.Simulating solute transport in porous or fractured formations using random walk particle tracking:a review.Vadose Zone Journal 2005;4(2):360-79.

Dou Z,Zhou Z,Sleep BE.In f l uence of wettability on interfacial area during immiscible liquid invasion into a 3D self-af f i ne rough fracture:lattice Boltzmann simulations.Advances in Water Resources 2013;61:1-11.

Eker E,Akin S.Lattice Boltzmann simulation of f l uid f l ow in synthetic fractures. Transport in Porous Media 2006;65(3):363-84.

Glover PWJ,Hayashi K.Modelling f l uid f l ow in rough fractures:application to the Hachimantai geothermal HDR test site.Physics and Chemistry of the Earth 1997;22(1-2):5-11.

Glover PWJ,Matsuki K,Hikima R,Hayashi K.Fluid f l ow in synthetic rough fractures and application to the Hachimantai geothermal hot dry rock test site.Journal of Geophysical Research:Solid Earth 1998a;103(B5):9621-35.

Glover PWJ,Matsuki K,Hikima R,Hayashi K.Synthetic rough fractures in rocks. Journal of Geophysical Research:Solid Earth 1998b;103(B5):9609-20.

Habich J,Feichtinger C,K?stler H,Hager G,Wellein G.Performance engineering for the lattice Boltzmann method on GPGPUs:architectural requirements and performance results.Computers and Fluids 2013;80:276-82.

Hassan AE,Mohamed MM.On using particle tracking methods to simulate transport in single-continuum and dual continua porous media.Journal of Hydrology 2003;275(3-4):242-60.

Hunt B.Dispersive sources in uniform ground water f l ow.Journal of the Hydraulics Division,ASCE 1978;104(1):74-85.

James SC,Chrysikopoulos CV.Transport of polydisperse colloid suspensions in a single fracture.Water Resources Research 1999;35(3):707-18.

James SC,Chrysikopoulos CV.Monodisperse and polydisperse colloid transport in water-saturated fractures with various orientations:gravity effects.Advances in Water Resources 2011;34(10):1249-55.

Lammers P,Beronov KN,Volkert R,Brenner G,Durst F.Lattice BGK direct numerical simulation of fully developed turbulence in incompressible plane channel f l ow. Computers and Fluids 2006;35(10):1137-53.

Latt J.Hydrodynamic limit of lattice Boltzmann equations.Geneva,Switzerland: University of Geneva;2007.

Latt J,Krause J.OpenLB user guide.Institute of Mechanical Engineering,Ecole Polytechnique Federale de Lausanne(EPFL);2008.

Maier RS,Kroll DM,Bernard RS,Howington SE,Peters JF,Davis HT.Pore-scale simulation of dispersion.Physics of Fluids 2000;12(8):2065-79.

Matsumoto M,Nishimura T.Mersenne twister:a 623-dimensionally equidistributed uniform pseudo-random number generator.ACM Transactions on Modeling and Computer Simulation 1998;8(1):3-30.

Mondal PK,Sleep BE.Colloid transport in dolomite rock fractures:effects of fracture characteristics,speci f i c discharge,and ionic strength.Environmental Science and Technology 2012;46(18):9987-94.

Nowamooz A,Radilla G,Fourar M,Berkowitz B.Non-Fickian transport in transparent replicas of rough-walled rock fractures.Transport in Porous Media 2013;98(3):651-82.

Ogilvie SR,Isakov E,Glover PWJ.Fluid f l ow through rough fractures in rocks.II:a new matching model for rough rock fractures.Earth and Planetary Science Letters 2006;241(3-4):454-65.

Patir N,Cheng HS.Application of average f l ow model to lubrication between rough sliding surfaces.Journal of Tribology 1979;101(2):220-9.

Qian J,Zhan H,Chen Z,Ye H.Experimental study of solute transport under non-Darcian f l ow in a single fracture.Journal of Hydrology 2011;399(3-4):246-54.

Reitsma S,Kueper BH.Laboratory measurement of capillary pressure-saturation relationships in a rock fracture.Water Resources Research 1994;30(4):865-78.

Renshaw CE.On the relationship between mechanical and hydraulic apertures in rough-walled fractures.Journal of Geophysical Research 1995;100(B12): 24629-36.

Saito M,Matsumoto M.Variants of Mersenne twister suitable for graphic processors.ACM Transactions on Mathematical Software(TOMS)2013;39(2). http://dx.doi.org/10.1145/2427023.2427029.

Stockman HW.A lattice gas study of retardation and dispersion in fractures: assessment of errors from desorption kinetics and buoyancy.Water Resources Research 1997;33(8):1823-31.

Stockman HW,Glass RJ,Cooper C,Rajaram H.Accuracy and computational ef f iciency in 3D dispersion via lattice-Boltzmann:models for dispersion in rough fractures and double-diffusive f i ngering.International Journal of Modern Physics C 1998;9(8):1545-57.

Succi S.The lattice Boltzmann equation for f l uid dynamics and beyond.Oxford: Clarendon Press;2001.

Sukop MC,Huang H,Alvarez PF,Variano EA,Cunningham KJ.Evaluation of permeability and non-Darcy f l ow in vuggy macroporous limestone aquifer samples with lattice Boltzmann methods.Water Resources Research 2013;49(1):216-30.

Sukop MC,Thorne DT.Lattice Boltzmann modeling:an introduction for geoscientists and engineers.Springer;2006.

Tatone BSA,Grasselli G.A method to evaluate the three-dimensional roughness of fracture surfaces in brittle geomaterials.Review of Scienti f i c Instruments 2009;80(12).http://dx.doi.org/10.1063/1.3266964.

T?lke J.Implementation of a lattice Boltzmann kernel using the compute uni f i ed device architecture developed by NVIDIA.Computing and Visualization in Science 2010;13(1):29-39.

Tompson AFB,Gelhar LW.Numerical simulation of solute transport in threedimensional,randomly heterogeneous porous media.Water Resources Research 1990;26(10):2541-62.

Tsang YW.The effect of tortuosity on f l uid f l ow through a single fracture.Water Resources Research 1984;20(9):1209-15.

Wels C,Smith L,Beckie R.The in f l uence of surface sorption on dispersion in parallel plate fractures.Journal of Contaminant Hydrology 1997;28(1-2):95-114.

Yan Y,Koplik J.Flow of power-law f l uids in self-af f i ne fracture channels.Physical Review E 2008;77(3):036315.

Zhu H,Liu X,Liu Y,Wu E.Simulation of miscible binary mixtures based on lattice Boltzmann method.Computer Animation and Virtual Worlds 2006;17(3-4): 403-10.

Zimmerman RW,Kumar S,Bodvarsson GS.Lubrication theory analysis of the permeability of rough-walled fractures.International Journal of Rock Mechanics and Mining Sciences and Geomechanics Abstracts 1991;28(4):325-31.

Zimmerman RW,Bodvarsson GS.Hydraulic conductivity of rock fractures.Transport in Porous Media 1996;23(1):1-30.

Zimmerman RW,Main I.Hydromechanical behavior of fractured rocks.In:Yves G, Maurice B,editors.International geophysics.Academic Press;2003.p.363-421.

Scott A.Briggs is a post-doctoral fellow at the University of Toronto.His current work focuses on the numerical modelling of sulphide diffusion in bentonite with applications for the long-term storage of nuclear waste.He obtained his doctorate in Environmental Engineering from the University of Toronto with a focus on contaminant hydrogeology.His Ph D research included modelling the impact of single fracture roughness on the f l ow,transport and development of bio f i lms.

*Corresponding author.Tel.:+1 416 978 5969.

E-mail address:scott.briggs@mail.utoronto.ca(S.Briggs).

Peer review under responsibility of Institute of Rock and Soil Mechanics,Chinese Academy of Sciences.

1674-7755?2014 Institute of Rock and Soil Mechanics,Chinese Academy of Sciences.Production and hosting by Elsevier B.V.All rights reserved.

http://dx.doi.org/10.1016/j.jrmge.2014.10.004