Three-dimensional analysis of spreading and mixing of miscible compound in heterogeneous variable-aperture fracture

Zhi Dou*,Zhi-fang Zhou,Jin-guo Wang

aDepartment of Civil Engineering,University of Toronto,Toronto M5S 1A4,Canada

bSchool of Earth Sciences and Engineering,Hohai University,Nanjing 210098,China

Three-dimensional analysis of spreading and mixing of miscible compound in heterogeneous variable-aperture fracture

Zhi Doua,b,*,Zhi-fang Zhoub,Jin-guo Wangb

aDepartment of Civil Engineering,University of Toronto,Toronto M5S 1A4,Canada

bSchool of Earth Sciences and Engineering,Hohai University,Nanjing 210098,China

As mass transport mechanisms,the spreading and mixing(dilution)processes of miscible contaminated compounds are fundamental to understanding reactive transport behaviors and transverse dispersion.In this study,the spreading and dilution processes of a miscible contaminated compound in a three-dimensional self-affine rough fracture were simulated with the coupled lattice Boltzmann method(LBM). Moment analysis and the Shannon entropy(dilution index)were employed to analyze the spreading and mixing processes,respectively.The corresponding results showed that the spreading process was anisotropic due to the heterogeneous aperture distribution.A compound was transported faster in a large aperture region than in a small aperture region due to the occurrence of preferential fl ow.Both the spreading and mixing processes were highly dependent on the fluid fl ow velocity and molecular diffusion.The calculated results of the dilution index showed that increasing the fl uid flow velocity and molecular diffusion coeffi cient led to a higher increasing rate of the dilution index.

Mixing;Spreading;Solute transport;Three-dimensional fracture;Self-affinity;Hurst exponent

1.Introduction

The investigation of spreading and mixing processes in the subsurface is important in many scientific disciplines,technical applications,and engineering practices,including contaminant hydrology,nuclear waste storage,groundwater remediation,and reservoir engineering(Qian et al.,2011; Yang et al.,2012;Mondal and Sleep,2013;Cai et al., 2010).In groundwater environments,the mass transport of amiscible contaminated compound is complex due to the geological formations in the subsurface.Characterization of the spreading and mixing processes of miscible compounds is signifi cant to quantifi cation and description of the chemical and biological reactions in both porous media and fractures. According to Rolle et al.(2009)and Kitanidis(1994), spreading caused by the variability of advection velocity refers to the change of compound cloud shape,while dilution refers to the change of water volume occupied by the solute,which is the only process allowing mass exchange between different streamlines and resulting in the decrease of the peak concentration in breakthrough curves.Insufficient or incomplete mixing of a compound often has a certain infl uence on the overalland localrates of chemicalreaction.In particular,since only the mixing process transverse to the fl ow directions is dispersive,the mixing process plays an important role in transverse dispersion of the compound cloud.Therefore,it is important to understand and quantify the spreading and mixing processes of compounds in the subsurface.

A large number of numericaland experimentalstudies have shown that the heterogeneous nature of geological formations has a signifi cant influence on the spreading and mixing processes of compounds.Rolle and Kitanidis(2014)used a porescale numerical model to study the effects of compoundspecific dilution on transient transport and solute breakthrough.Because characterizing a dilution process from integrated measurements of solute breakthrough is challenging, they modifi ed the dilution index,introduced by Kitanidis (1994),and proposed a transient fl ux-related dilution index to measure the dilution breakthrough curve.Hochstetler et al. (2013)studied the effects of compound-specific transverse mixing on steady-state reactive plumes by comparing porescale simulations with Darcy-scale experiments.Ballarini et al.(2014)evaluated the effect of heterogeneities on transverse mixing in bench-scale tank experiments by analyzing the spatialmomentand dilution index of the compound tracer.The results showed that the magnitude of permeability had a signifi cant influence on the spreading and dilution processes;as in heterogeneous porous media,a non-uniform fl ow velocity field occurred in the natural fracture due to the heterogeneous aperture distribution.Cirpka et al.(2015)described the quantitative assessment of transverse mixing and its enhancement in three-dimensional(3D)heterogeneous anisotropic porous media.They used a quantitative method to analyze transverse mixing based on the spreading and dilution processes.Although many studies have shown that the heterogeneity of porous media has a significant infl uence on the spreading and dilution processes,and,consequently,transverse dispersion is highly dependent on the heterogeneity of the geological formation,studies of the spreading and dilution processes in rough fractures are still limited.Recently,Dou and Zhou(2014)described the miscible compound transport in a two-dimensional(2D)single rough fracture.The results showed that the slow mass exchange through solute molecular diffusion between mobile and immobile regions,which was caused by the roughness of the fracture,resulted in a long breakthrough tailing.However,they did not analyze the spreading and mixing processes.

The Navier-Stokes equation has typically been used to solve the non-uniform fl ow velocity fi eld.Recently,lots of numerical methods have shown that,as an alternative to directly solving the Navier-Stokes equation,the lattice Boltzmann method(LBM)can simulate the non-uniform fl ow velocity by solving the discrete Boltzmann equation(Chen et al.,2013).Much effort has been devoted to designing different boundary conditions for the LBM,such as the periodic boundary condition,pressure boundary condition,and bounce-back boundary condition.For the geometric model used in this study,the LBM has natural advantages over other conventional computational fluid dynamics(CFD)methods, especially in dealing with the boundary of a self-affi ne rough fracture wall.In addition,Dou and Zhou(2014)have proved that a coupled LBM model is capable of simulating the miscible compound transport.

The main objective of this study was to simulate the spreading and mixing processes of a compound in a 3D rough fracture and investigate the influence of non-uniform flow velocity and molecular diffusion on the spreading and mixing processes.The non-uniform fl ow velocity and compound concentration fi elds were obtained by solving the Navier-Stokes equation and the advection-diffusion equation through a coupled LBM model,respectively.A 3D rough fracture was generated using the technique of successive random additions. The spreading and mixing processes of the compound were analyzed by calculating the spatial moment and original dilution index for different molecular diffusion coeffi cients in different average fl uid fl ow velocity fields.

2.Method

2.1.Numerical model

In this study,a coupled LBM model was developed to simulate the miscible compound transport in a 3D rough fracture.The coupled LBM model used two particle distribution functions to represent the pure water and the miscible compound,respectively:

where fi(X，t)and gi(X，t)are the distribution functions of fl uid particles and compound particles,respectively,at spatial position X and time t with velocity vector ei,respectively;the subscript i is the number of particles;τandτDare the nondimensional relaxation times related to the kinematic viscosity,and the molecular diffusion coefficient,respectively;csis the lattice sound speed;and feqi(X，t)and geqi(X，t)are the equilibrium distribution functions for fluid particles and compound particles, respectively,defined as

where ueqis the local equilibrium velocity,ρfis the fl uid density,and C is the compound concentration.In this study, the D3Q19 model and corresponding weightωiwere employed to simulate the particle movement for each node. The fl uid density,compound concentration,fluid pressure p, and fluid velocity u in the absence of any additionalforces are given by the following:

Using the Chapman-Enskog expansion,the Navier-Stokes equation can be recovered with second-order accuracy from the coupled LBM model for the fl uid flow velocity field(Eq.(1)). Similarly,the compound particle distribution function(Eq.(2)) can be recovered with the macroscopic advection-diffusion equation for the compound concentration field.Thus,governing equations in the model that couple the flow and concentration fi elds are as follows:

whereμis the dynamic viscosity of water.The coupled LBM model was developed as a C++code.The validation of this coupled LBM model has been introduced in Dou and Zhou (2014).Any lattice node in the computational domain represented either a solid node or a fl uid node.For the solid node, the collision step in the half-way bounce-back algorithm was implemented,instead of the usual collision step,to provide a non-slip wall boundary condition.

2.2.Geometric model

In this study,a 3D self-affine fracture surface was generated with the technique of successive random additions(Voss, 1988),an efficient and fast algorithm.The height of the selfaffi ne rough fracture wall surface,h(x，y),satisfi es the conditionλHh(x，y)=h(λx，λy),whereλis the scaling factor,and H is the Hurst exponent.Unlike self-similar structures,the scaling transformation for the self-affi ne surface is isotropic in the horizontal plane(xy plane),but anisotropic in the vertical direction(z direction),with the Hurstexponent H varying from 0 to 1.Fig.1 shows the 3D self-affine surfaces with H=0.85 and H=0.40,respectively,where lu represents one lattice spacing.

Fig.1 shows that the fracture surface with a lower H value is relatively rough,and the variation of the surface height is remarkable.A higher H value results in a higher spatial correlation and a relatively smooth fracture surface.However, there is ample experimental evidence that the natural fracture surface has the characteristics of self-affi nity,with the Hurst exponent of H≈0.8±0.05(M?l?y et al.,1992;Schmittbuhl et al.,2008).

Fig.1.Self-affine rough fracture surfaces with different Hurst exponents.

3.Results and discussion

In this study,a computational domain,with lattice units of 67,67,and 22 in the x,y,and z directions,respectively,was used to represent a 3D rough fracture.H=0.86 and H=0.82 were set for the top and bottom self-affi ne fracture surfaces, respectively.The corresponding statistical histogram of aperture distribution in the 3D rough fracture is shown in Fig.2.It can be seen that the aperture distribution in the 3D rough fracture is a Gaussian-type distribution.

In the LBM,the lattice spacing and time step are usually dimensionless and typically set to unity.In this study,one lattice spacing was set to be 3.6μm and one lattice time step was set to be 2.2×10-6s,while the kinematic viscosity was set to be 1/6.

Fig.2.Distribution of apertures in 3D rough fracture.

For the initial simulation conditions,a miscible contaminated compound,with an initialconcentration of C0=1.0,was injected instantaneously into the region of 7.2μm＜x＜18μm, and an initial concentration C0was set to be 0.001 in other computational domain to avoid numerical instability.Two different average fl ow velocities were used,hereafter referred to as the high-flow velocity and low-flow velocity.In order to simulate the high-and low-flow velocities,two constant body forces,G=1×10-5and G=1×10-4,were introduced to substitute for the pressure gradients of fluid and accelerate the fluid flow velocity field,respectively.For the fluid fl ow velocity and concentration fi elds,the periodic boundary condition was imposed at the inlet and outletto generate a conservative mass transfer throughout the computationaldomain.If the variations of the spatial moment and dilution index over 4000ts(ts represents one lattice time step)were less than 0.1%,it was assumed that the spreading and dilution processes had arrived at a quasi-equilibrium state.The time step for each numerical case was about 0.42 s.This computation took about 8 h to run on an Intel CPU(2.1 GHz)computer for a computational domain of 241.2μm×241.2μm×79.2μm.

3.1.Flow field in 3D rough fracture

The fl ow fi eld in the 3D rough fracture was simulated by adding a constant body force along the assumed flow direction,as shown in Fig.3.For clarity of illustration,the velocity slice planes parallel,respectively,to the xy,yz,and zx planes, were extracted.The transparent white contour represented the self-affine surface.

Fig.3 shows the non-uniform flow fi eld in the 3D rough fracture.There is a preferential flow path in the flow direction,where the flow velocity is higher than those in other regions of the 3D rough fracture.Although the non-uniform flow velocity fi eld was solved with the Navier-Stokes equation in this study,the cubic law was used to analyze the relationship between the spatial aperture distribution and local fluid flow velocity.The cubic law for the single rough fracture was defined as where Q is the bulk flow rate,J is the hydraulic gradient,g is the gravitational acceleration,and b is the mean fracture aperture.According to the cubic law,the local fluid fl ow velocity is highly dependent on the localaperture.For a constant body force,a smaller aperture resulted in a lower fluid flow velocity.The Reynolds number for the performed 3D fracture was defi ned as follows(Qian etal.,2015;Brush and Thomson, 2003):

Fig.3.Flow fi eld in 3D rough fracture with body force of 1×10-5.

where lvis the characteristic length of viscous force,and UIis the characteristic velocity due to the inertial force.In this study,lvwas defined as the mean fracture aperture b,and UIwas defi ned as UI=Q/bW,where W is the fracture width.In the simulation,the maximum Reynolds number during a periodic fl ow fluctuation was Re=53.Assuming that there was Darcy flow in all of the cases studied,for the compound transport,the dimensionless Peclet number(Pe)is

In this study,four numerical cases with different fluid flow velocities and molecular diffusion coeffi cients were investigated.For the low-fl uid flow velocity,the values of Pe were 24 and 70 for the high and low molecular diffusion coefficients, respectively.For the high-fl uid flow velocity,the values of Pe were 240 and 700 for the high and low molecular diffusion coeffi cients,respectively.

3.2.Spatial moments

The general expression of the spatial moments for the concentration distributionρg(x，y，z，t)is

where Mijk(t)is the spatial moment with a sum of i,j,and k equal to 0,1,or 2,respectively;and n is the porosity of the porous media(Freyberg,1986).For a single rough fracture,n was assumed to be 1.0.Physically,the zeroth spatial moment forthe concentration distribution is a measure of the totalmass of the compound.The first spatial moment for the concentration distributionρg(x，y，z，t),normalized by the zeroth spatial moment,defines the center of mass of the compound (xρg，yρg，zρg):

The components of the second spatial moment in the x,y, and z directions atthe location of(x,y,z)are quantified by theprincipal diagonal elements(σxx，σyy，andσzz)in the spatial covariance tensor,respectively,given by the following:

Eqs.(19)through(21)show the spreading of the compound around the center of the compound mass in the x,y,and z directions,respectively.

Figs.4(a)through(c)show the first spatial moment of the compound in the main flow direction,transverse direction,and vertical direction.The normalized first spatial moment representsthe centerofmassofthe compound.The resultsin the main fl ow direction show that the compound is transported faster in the high-fl uid fl ow velocity fi eld than in the low-fluid flow velocity fi eld.Forthe cases in the high-fluid flow velocity fi eld,the infl uence of molecular diffusion on the fi rst spatial moment in the main flow direction is negligible.However,for the cases in the low-fluid flow velocity field,it can be seen from Figs.4(a) through(c)thatthe spreading ofthe compound in the main fl ow, transverse,and verticaldirections is significantly influenced by moleculardiffusion.Ahighermoleculardiffusion leads to faster transport in the main flow,transverse,and vertical directions. Moreover,itshould be noted thatthe spreading ofthe compound in this study arrived atthe quasi-equilibrium state while the fi rst and second spatialmoments appeared constant.This was due to the limited computational domain and periodic boundary conditions.Under such conditions,the spreading process of the compound lasts a longer time in the transverse direction than in the main fl ow and vertical directions.

The second spatial moment indicates the compound spreading around its center of mass.In Figs.4(d)through(f), the corresponding results show thatboth the fl uid fl ow velocity and molecular diffusion have signifi cant influences on the second spatial moment.At the beginning of the spreading of the compound,the second spatial moment decreases in the transverse direction,and increases in the main flow and verticaldirections.This is because the compound spreading in the transverse direction is larger in the initial state than in the quasi-equilibrium state.In general,the process through which the second spatial moment arrives at the quasi-equilibrium state is faster in the high-fl uid flow velocity fi eld than in the low-fl uid fl ow velocity field.For the cases in the low-fl uid fl ow velocity field,the effects of molecular diffusion in the main flow,transverse,and vertical directions are remarkable. However,the effects weaken in the high-fl uid fl ow velocity field.It can be seen from Fig.4 that there are slight fl uctuations of the fi rst and second spatial moments in the high-fluid flow velocity fi eld.This is due to the re-equilibrium process of the compound concentration distribution among the apertures in the 3D rough fracture.The re-equilibrium process is sensitive to the heterogeneous aperture distribution.The frequency of the fl uctuation is lesser in the low-fl uid fl ow velocity field than in the high-fl uid fl ow velocity field.

3.3.Dilution index

Fig.4.First and second spatial moments of compound in different directions.

Although the spreading of the compound can be quantifi ed with the fi rstand second spatialmoments in both homogenousand heterogeneous fl ow velocity fields,the mixing of the compound in a heterogeneous flow velocity field,such as in a self-affine rough fracture,cannot be measured with the second spatialmoment.On the one hand,the mechanism of spreading is differentfrom thatofdilution;the spreading in heterogeneous systems is inadequate to represent a measure of dilution or mixing.On the other hand,when a compound is transported inside a preferential channel,its second spatial moment along the transverse direction decreases while the actualdilution does not decrease because dilution is an irreversible process(Rolle et al.,2009;Ballarini et al.,2014).In this study,the original dilution index(Shannon entropy),introduced by Kitanidis (1994),was used to quantify the global compound dilution in the 3D self-affi ne rough fracture,given by the following:

where pg(x，y，z，t)is the concentration probability function of the compound,defi ned as follows:

The maximum value of the originaldilution index,E(t)max, is technically equal to the volume of the entire fracture.Then, the dimensionless dilution index,representing a non-uniform concentration distribution in the fracture,can be obtained through division of the dilution index by E(t)max.

Fig.5.Dilution index for compound with different molecular diffusion coefficients in different fl uid flow velocity fi elds.

Fig.5 shows the variation of the dimensionless dilution index in the compound mixing process.It can be seen from Fig.5 thatboth the fl uid flow velocity and molecular diffusion have signifi cant infl uences on the global dilution index.In the high-fluid flow velocity field,the increasing rate of the dilution index increases with the molecular diffusion coeffi cient. The same results can be seen in the low-fluid fl ow velocity fi eld.The influence of molecular diffusion on the mixing process is greater in the low-fl uid fl ow velocity field than in the high-fluid fl ow velocity field.The increasing fl uid flow velocity results in the increased deformation of the compound,which enhances the mass exchange between different apertures.

Fig.6.Mixing processes of compound under different conditions.

Fig.6 shows the non-uniform concentration distribution of the compound in the 3D rough fracture due to theheterogeneous aperture distribution.The geometry of each slice for the concentration distribution in the 3D rough fracture is dependent on the rough fracture surface.In Fig.6,the concentration distribution is highly sensitive to apertures.For different Pe values,the compound transport is faster in larger aperture regions.Since the initial concentration and total compound mass are the same for all numerical cases,a higher peak concentration is associated with a smaller magnitude of dilution.Thus,for a local aperture region,we can simply compare the magnitude of dilution depending on the peak concentration.In Fig.6,it can be seen that the peak concentration decreases with time.Furthermore,comparing Figs.6(a)and(c)shows that the mixing process of the compound is dominated by the fluid fl ow velocity.The compound is more diluted in the large aperture region than in the small aperture region.Preferentialflow occurs in the large aperture region and enhances the mixing process.Comparing Figs.6(b)and(d)shows that the concentration is evenly distributed over the whole fracture for the larger Pe value. However,it can be seen from Fig.5 that the dilution index is below 1.0,indicating the occurrence of incomplete mixing.In fact,although increasing both the fluid fl ow velocity and molecular diffusion coefficient leads to the increase of the compound concentration in the small aperture region,the concentration distribution is non-uniform in large and small aperture regions.The non-uniform concentration distribution caused by the heterogeneous aperture distribution results in incomplete mixing throughout the rough fracture.

4.Conclusions

In this study,the spreading and mixing processes of a compound in a 3D rough fracture were simulated to investigate the influence of non-uniform flow velocity and molecular diffusion in these processes.The non-uniform flow velocity and compound concentration fields were obtained by solution of the Navier-Stokes equation and the advection-diffusion equation,respectively,through the coupled LBMmodel.Some conclusions are as follows:

The presented coupled LBMmodelis capable of simulating and analyzing the spreading and mixing processes of a miscible contaminated compound in the 3D self-affi ne rough fracture.Although both the fluid flow velocity and molecular diffusion have signifi cant influences on the spreading and mixing processes,the heterogeneous aperture distribution is the dominantfactor in the processes.It leads to a non-uniform fl uid flow velocity fi eld in the rough fracture,anisotropy in the spreading of the compound,and incomplete mixing.Due to the occurrence of preferential flow,the compound is transported faster and is more deformed and stretched in the large aperture region than in the small aperture region,resulting in differences in the spreading and mixing among apertures.In addition,increasing the fluid flow velocity and molecular diffusion coefficient leads to a higher increasing rate of the dilution index.

Ballarini,E.,Bauer,S.,Eberhardt,C.,Beyer,C.,2014.Evaluation of the role ofheterogeneities on transverse mixing in bench-scale tank experiments by numerical modeling.Groundwater 52(3),368-377.http://dx.doi.org/ 10.1111/gwat.12066.

Brush,D.J.,Thomson,N.R.,2003.Fluid fl ow in synthetic rough-walled fractures:Navier-Stokes,Stokes,and local cubic law simulations.Water Resour.Res.39(4),1085.http://dx.doi.org/10.1029/2002WR001346.

Cai,J.C.,Yu,B.M.,Zou,M.Q.,Mei,M.F.,2010.Fractal analysis of surface roughness of particles in porous media.Chin.Phys.Lett.27(2),024705. http://dx.doi.org/10.1088/0256-307X/27/2/024705.

Chen,Q.,Zhang,X.,Zhang,J.,2013.Improved treatments for general boundary conditions in the lattice Boltzmann method for convectiondiffusion and heat transfer processes.Phys.Rev.E 88(3),033304.http:// dx.doi.org/10.1103/PhysRevE.88.033304.

Cirpka,O.A.,Chiogna,G.,Rolle,M.,Bellin,A.,2015.Transverse mixing in three-dimensional nonstationary anisotropic heterogeneous porous media. Water Resour. Res. 51(1), 241-260. http://dx.doi.org/10.1002/ 2014WR015331.

Dou,Z.,Zhou,Z.F.,2014.Lattice Boltzmann simulation of solute transportin a single rough fracture.Water Sci.Eng.7(3),277-287.http://dx.doi.org/ 10.3882/j.issn.1674-2370.2014.03.004.

Freyberg,D.L.,1986.A natural gradient experiment on solute transport in a sand aquifer:2.Spatial moments and the advection and dispersion of nonreactive tracers.Water Resour.Res.22(13),2031-2046.http:// dx.doi.org/10.1029/WR022i013p02031.

Hochstetler,D.L.,Rolle,M.,Chiogna,G.,Haberer,C.M.,Grathwohl,P., Kitanidis,P.K.,2013.Effects of compound-specific transverse mixing on steady-state reactive plumes:Insights from pore-scale simulations and Darcy-scale experiments.Adv.Water Resour.54,1-10.http://dx.doi.org/ 10.1016/j.advwatres.2012.12.007.

Kitanidis,P.K.,1994.The concept of the dilution index.Water Resour.Res. 30(7),2011-2026.http://dx.doi.org/10.1029/94WR00762.

M?l?y,K.J.,Hansen,A.,Hinrichsen,E.L.,Roux,S.,1992.Experimental measurements of the roughness of brittle cracks.Phys.Rev.Lett.68(2), 213-215.http://dx.doi.org/10.1103/PhysRevLett.68.213.

Mondal,P.K.,Sleep,B.E.,2013.Virus and virus-sized microsphere transport in a dolomite rock fracture.Water Resour.Res.49(2),808-824.http:// dx.doi.org/10.1002/wrcr.20086.

Qian,J.,Ma,L.,Zhan,H.,Luo,Q.,Wang,X.,Wang,M.,2015.The effect of expansion ratio on the critical Reynolds number in single fracture flow with sudden expansion.Hydrol.Process.30(11),1718-1726.http:// dx.doi.org/10.1002/hyp.10745.

Qian,J.Z.,Zhan,H.B.,Chen,Z.,Ye,H.,2011.Experimental study of solute transportunder non-Darcian flow in a single fracture.J.Hydrol.399(3-4), 246-254.http://dx.doi.org/10.1016/j.jhydrol.2011.01.003.

Rolle,M.,Eberhardt,C.,Chiogna,G.,Cirpka,O.A.,Grathwohl,P.,2009. Enhancement of dilution and transverse reactive mixing in porous media: Experiments and model-based interpretation.J.Contam.Hydrol.110(3), 130-142.http://dx.doi.org/10.1016/j.jconhyd.2009.10.003.

Rolle,M.,Kitanidis,P.K.,2014.Effects of compound-specific dilution on transient transport and solute breakthrough:A pore-scale analysis.Adv.Water Resour.71,186-199.http://dx.doi.org/10.1016/j.advwatres.2014.06.012.

Schmittbuhl,J.,Steyer,A.,Jouniaux,L.,Toussaint,R.,2008.Fracture morphology and viscous transport.Int.J.Rock Mech.Min.Sci.45(3), 422-430.http://dx.doi.org/10.1016/j.ijrmms.2007.07.007.

Voss,R.F.,1988.Fractals in nature:From characterization to simulation.In: Barnsley,M.F.,Devaney,R.L.,Mandelbrot,B.B.,Peitgen,H.O.,Saupe,D., Voss,R.F.,eds.,The Science of Fractal Images.Springer New York,New York,pp.21-70.

Yang,Z.,Niemi,A.,Fagerlund,F.,Illangasekare,T.,2012.A generalized approach for estimation of in-plane curvature in invasion percolation models for drainage in fractures.Water Resour.Res.48(9),W09507. http://dx.doi.org/10.1029/2012WR011829.

Received 5 July 2016;accepted 15 September 2016

Available online 21 January 2017

This work was supported by the National Natural Science Foundation of China(Grant No.41602239),the Natural Science Foundation of Jiangsu Province(Grant No.BK20160861),the Fundamental Research Funds for the Central Universities(Grant No.2016B05514),the International Postdoctoral Exchange Fellowship Program from the Office of China Postdoctoral Council (Grant No.20150048),and the“333 Project”of Jiangsu Province(Grant No. BRA2015305).

*Corresponding author.

E-mail address:Dz.uriah@gmail.com(Zhi Dou).

Peer review under responsibility of Hohai University.

?2016 Hohai University.Production and hosting by Elsevier B.V.This is an open access article under the CC BY-NC-ND license(http:// creativecommons.org/licenses/by-nc-nd/4.0/).