Estimation of hydraulic anisotropy of unconsolidated granular packs using finite element methods

The effect of particle shape and heterogeneity on hydraulic anisotropy of unconsolidated granular packs is hereby investigated. Direct simulation was carried out on synthetically generated spherical, aspherical, ellipsoidal (aspect ratio of 2 and 3) and lenticular samples. Single phase Stokes equation was solved on models discretised on finite element geometries and hydraulic permeability computed in the horizontal and vertical directions to estimate the degree of anisotropy. The spherical and aspherical packs with varying degrees of particle shapes and heterogeneities are virtually isotropic. Ellipses with aspect ratios 2 and 3 have higher anisotropy ratios compared to the spherical and aspherical geometries while the lenticular geometry is the most anisotropic. This is attributable to the preferential alignment of the grains in the horizontal flow direction during random dynamic settling under gravity.

Estimation of hydraulic anisotropy of unconsolidated granular packs using finite element methods

INTRODUCTION
Anisotropy (K h /K v ) is a measure of the directional change in permeability of a rock sample.Anisotropy controls single and multi-phase fluid flow effective mobility through porous systems.Small scale anisotropy in geologic formations may occur due to: stratification, directional rock fractures or orientation of non-spherical grains.
The hydraulic anisotropy of soils and sedimentary rocks has a great effect on fluid flow and contaminant transport.Hence, an understanding of the anisotropy ratio is important for many subsurface fluid-associated problems such as the design of oil or water well fields, underflow beneath dams and dykes, internal erosion in soil masses or settlement rates of consolidating clays.
In the field, anisotropy is mainly due to rock fracturing and sediments forming strata.Induced anisotropy by continuous or discontinuous isotropic layers has been studied extensively [1][2][3][4][5].The geo-statistics approach to the permeability of heterogeneous media was developed by [6,7] and used in numerical simulations [8].Small scale anisotropy may occur from stratification due to the process of formation of the rock [9], directional rock fractures [10] or orientation of non-spherical grain particles during deposition [11].
Anisotropy also affects production from hydrocarbon reservoirs; the larger the degree of anisotropy the higher the productivity index.When vertical permeability is low, wells completed horizontally may become economically unattractive.Wells drilled and completed normal to the larger horizontal permeability will potentially be a much better producer than one drilled and completed arbitrarily or normal to the smaller horizontal permeability.It is therefore pertinent to measure permeability before the horizontal section is drilled [30].
A lot of work has been done on the anisotropy ratio of cohesive soils and sedimentary rocks using probe (mini) permeameter (e.g.[31,32]).However, little has been done on how grain shape and packing heterogeneity influence the hydraulic anisotropy of unconsolidated granular materials.
In order to investigate the effect of particle shape on hydraulic anisotropy, several numerical experiments were conducted on synthetic granular packs consisting of spherical, aspherical, ellipsoids; with aspect ratios 2 and 3 and lenticular geometries.During the dynamic settling of the granular packs, the aspect ratio was varied and the degree of particle alignment was determined using the nematic order parameter [33].

CONSTRUCTION OF SYNTHETIC GRANULAR PACKS AND SETTLING PROCEDURE
A detailed description of how the synthetic granular packs used in this study were constructed has been presented in [33].Basically, the shapes of particles are constructed using the information derived from the surface mesh obtained via scanning technologies.A model particle is then constructed by sequentially adding spheres to the volume enclosed by the mesh.As more spheres are added, the enclosing surface of the cluster approximates the surface of the real particle represented in the original mesh.
The grain packs are then constructed by allowing the grains to fall under gravity into a confined box while allowing interaction between the grains and the bottom wall.When the kinetic energy of the system reaches a negligible value and every grain contacts either the wall or another grain, new batch of grains settles in the box above those previously deposited.This process repeats itself until the initial box is filled; the grains are not frozen at any stage.Hence, each new subset of grains entering the box interacts with previously deposited grains, allowing for future relocation of already settled grains.This extra dynamics allow the construction of denser packs of real sands.The porosities of the resulting packs depend on the particle shape and the friction coefficient used for the grain-grain interaction.Figures 1(a-d) show packs of spherical and aspherical grains with varying degree of heterogeneity.Figures 2(a-b) show ellipsoidal grain packs of aspect ratio of 2 and 3 and Figure 2(c) shows a lenticular pack geometry.

PORE-SPACE CAD MODEL CONSTRUCTION
The digitised binary representation of the porous media are processed and differentiated into GRAINS and PORES using Computer-Aided (CAD) geometrical package accomplished with Non-Uniform Rational B-splines (NURBS) curves and surfaces (e.g.Figures 1 and 2).
The NURBS curves and surfaces allow to capture the geometries with a tolerance-based (typical absolute tolerance is 1e −7 ) level of detail, independent of scale allowing a purpose dependent adaptation of the mesh to smooth geometry.Volumetric objects are then defined by grouping curve-delimited surfaces together by a technique called Boundary Representation (BREP).
This refers to a hierarchical, internally consistent tree structure of points (nodes), holes and surfaces (loops) and surface enclosed pore volumes (body), recording their relations to each other.Fourth order splines are used to represent grains and outer boundaries are resolved into surfaces with six side boundaries; four no-flow and two Dirichlet conditioned boundaries, for further import to the geometry editor of the mesh generation code.For all the CAD models constructed here, the knots in the NURBS curve are non-uniform and of degree 3.

BOUNDARY CONDITIONS FOR UNCONSOLIDATED GRANULAR PACKS
A sub-volume far from the boundaries of the porous media packs is extracted to avoid boundary effects such as distortions in the solid skeleton or zones of abnormal porosity.The dimensions of the granular packs are approximately 6 − 7 times the mean grain diameter in all directions.The part of the grains that are not contained in the sub-volume are discarded.The packs are then differentiated into two domains: Grains and Pores.In order to simulate flow, two Dirichlet and four no-flow boundary conditions (Figure 3), were defined in the 'CAD' geometrical package accomplished with NURBS curves and surfaces.

FINITE ELEMENT DISCRETISATION OF THE PORE-SPACE
In order to discretise the pore-space of all the geometric samples investigated, an unstructured hybrid mesh consisting of tetrahedral, hexahedral, prism and pyramid elements (Figure 4) was used.The unstructured grids can fit free-form geometrical entities, such as NURBS with spatially variable refinement and they can also be generated automatically.For realistic hybrid meshes of free-form geometry, the quality of the resulting mesh can be evaluated by using the element to node ratio.A value close to 2 can be obtained for hybrid meshes when compared with 5 − 6 for pure tetrahedral meshes [34].
In the geometrically unconstrained regions, hexahedral dominated elements were used.More shape-adaptive tetrahedral elements are used to capture geometric complexities and intentional refinement variations.These perform well even if they have a large aspect ratio [35].This dramatically reduces the number of elements required to represent pore geometries.

FLUID FLOW EQUATIONS
The partial differential equations governing the flow of an incompressible Newtonian fluid are the Navier -Stoke's equation (1) and the law of mass conservation for incompressible fluids can be written as (2) The variables ∇P, ρ, µ and reduced pressure gradient, fluid density, fluid viscosity and velocity vector respectively.
is the partial derivative with respect to time.The boundary conditions include: no slip at any boundary between the fluid and the solid; this implies that at the grain walls, the normal and tangential components of the velocity are zero.The only body force accommodated is gravity.This is defined through the reduced pressure gradient [ 36,37] Where p is the pressure and is a unit vector in the vertical direction.
Since the relationship between permeability and pore geometry is most readily studied using steady-state flow, transient effects are ignored and fluid density is assumed to be constant.This implies that at any fixed point in space, the velocity does not vary with time and the equations reduce to (3) The presence of the advective acceleration (second term on the left hand side) causes the equations to be nonlinear and consequently very difficult to solve.
In the case where the fluid is incompressible and this term is far smaller than the viscous term (first term on the left), which is a typical situation in flows where the fluid velocities are very slow (creeping flow), the viscous forces are very large, or the length-scales of the flow are very small so that the Reynolds number (Re = uL/µ, where L is a characteristic length) is small (typical Re values for flow in geological porous media are 10 −5 or lower), this equation reduces to the Stoke's equation [e.g.38, 39]: (4)

NUMERICAL SOLUTION SEQUENCE
Eqn 4 is solved numerically in three-dimensional pore-spaces by first solving for a function ψ (x, y, z) and then computing the pressure gradient ∇P.Considering a Newtonian fluid; which exhibits a linear stress-strain relationship with constant viscosity, the velocity can be written as: (5) The function ψ (x, y, z) can be obtained simply by solving the Poisson's equation Estimation of hydraulic anisotropy of unconsolidated granular packs using finite element methods with ψ(x, y, z) = 0 at the grain boundaries (homogeneous boundary conditions).
Given that the parabolic profile within the pore space is adequately captured [40]; with no-slip boundary condition at the grain surface, this two-step approach is less restrictive and basically approximates the Stoke's equation.
To compute the pressure field, the equation ( 7) is solved once ψ is known.This ensures that for an incompressible fluid considered here, the divergence of is zero thus: Since the computed velocity is piece-wise constant from finite-element to finite-element and discontinuous across the element boundaries; a complementary node-centered finitevolume (NCFV) is used to measure .This approach is detailed in [41].

FINITE ELEMENT DISCRETISATION OF THE FLOW EQUATION
In order to compute the velocity field on each finite element, a linear finite element method with linear interpolation function N was adopted where the pressure field (Ph) is approximated by expansions of the form: where P k is scalar of nodal value of pressure for nodes k = 1...n.The integral form of Eqn 8 for an incompressible fluid over a typical element Ωe bounded by Γe is now written as (10) In this Bubnov-Galerkin method [42,43], the weighting functions are the same as the interpolation functions.Introducing the expansion for P (Eqn 9), we have (11) The left hand side of the equation is expanded by partial integration after which Green's theorem is applied to obtain (12) The above equation can be written in a n × n system of linear equations of the form The solution of the above equation requires that the value of the pressure is fixed at the inlet and outlet of the model (Dirichlet boundary conditions) while no slip boundary conditions are assigned on all the other faces.The resulting matrix of linear algebraic equations is solved with the algebraic multigrid solver [44].The results are stored at the finite element nodes.Given that shape function derivatives are constant (since linear shape functions are used) within each element, and the parabolic function Ψ as well as the fluid viscosity µ are constant then, the product of these and the nodal values of pressure gives a constant value of velocity at the barycenter of the element.For a range of other finite element shapes, the interpolation functions can be found elsewhere [45,46].

HYDRAULIC PERMEABILITY AND ANISOTROPY
Having computed the velocity fields as described in section V, the permeability, κ, is computed thus: (14) where ∆P is the pressure difference per unit length along the direction of main flow, A is the cross-sectional area open to flow and q is the cross-sectional flux.
The flux q, (in m 3 s −1 ) through a finite volume cell is defined as (15) and is accummulated by summing the fluxes across the facets of each sector of the FV: (16) for n sectors and m facets.Thus, the local piece-wise constant velocity is projected onto facet normals, n.These dot products are multiplied by facet area and summed over all finite volume facets to complete the surface integral.Green's theorem establishes the equivalence of volume (V ) and surface (S) integrals of the form (17) The effective horizontal hydraulic permeability (K h ) and the effective vertical hydraulic permeability (K v ) are computed as shown in Figure 3 and the anisotropy ratio (A r ) determined from the ratio K h /K v and results summarised in Table I.

POROSITY
The porosity is established as a by-product of the discretisation of the pore-space and is computed using the following relations: Pore Volume (Vp) / Bulk Volume (Vb) where, Vb = (Grain Volume + Pore Volume), i.e. q An u = ⋅ κ µ = q A P ∆

DISCUSSION
The spherical (sample C) and aspherical (samples A, B, E) with varying degrees of particle shapes and heterogeneities and aspect ratio of ≈ 1, are virtually isotropic (Table I).As the aspect ratio increases, the preferential alignment increases resulting in higher hydraulic anisotropy.Porosity also plays a great role in the manner in which grains are aligned and it is therefore considered to have a secondary influence on hydraulic anisotropy.Elongated geometries (F, G) with aspect ratios 2 and 3 have higher A r compared to the spherical and aspherical grain packs (samples A, B, C, and E) with aspect ratios 1 (Figure 8).[14] and sand [20] are also plotted.

EFFECT OF ORIENTATIONAL ORDER ON HYDRAULIC ANISOTROPY
The possibility of the grains to rotate and rearrange during settling under gravity favours preferential orientation of the elongated particles.It has been shown [33] that though particle shape and packing heterogeneity have small but noticeable effects on permeability, porosity is the most important parameter affecting it.The results of simulation (Figure 8) clearly show that, for a fixed porosity, the higher the aspect ratio, the greater the preferential alignment; indicating the influence of grain shape and preferential alignment on hydraulic anisotropy.Packings with the same degree of heterogeneity have higher anisotropy at higher aspect ratio.In the bedding plane of sedimentary rocks, the A r is usually lower than 1.5 which means that these rocks are nearly isotropic in that plane while it is usually lower than 4 for clay, with a few exceptions up to 42 [20].Few results are available for granular materials, their A r is not always higher than 1.Experimental values for sands and gravels are in the range of 0.75 to 4.1.
It is known that Clays, either flocculated or dispersed, have different hydraulic anisotropies at the same void ratio [14,47].Similarly, the same granular material, either statically or dynamically compacted, has different hydraulic anisotropies at the same void ratio.The Ar obtained from our simulation for the spherical (e.g.pack C with regular grains and aspect ratio 1) and aspherical (e.g.packs A, B, E with irregular grains and aspect ratio 1) geometries are similar to values reported for sands and gravels at high porosity in [20].However, the higher values observed for the ellipsoids and lenticular geometries (e.g.packs F, G, H) suggest that particle shapes and preferential alignment in the horizontal flow direction have significant influence on the Ar (Figure 9).This further explains why clays in a 'loose' condition like the sensitive clays of Eastern Canada [21] or Sweden [48] are nearly isotropic, whereas, the over-consolidated homogeneous London clay has a value close to 2 Estimation of hydraulic anisotropy of unconsolidated granular packs using finite element methods

CONCLUSIONS
Direct simulation was carried out on synthetically generated spherical, aspherical, ellipsoidal (aspect ratio of 2 and 3) and lenticular samples to estimate the degree of anisotropy.As the porosity of the samples ϕ decreases; the degree of anisotropy increases significantly for the ellipsoidal (aspect ratio of 2 and 3) and lenticular samples.The spherical and aspherical geometries of varying degree of heterogeneities are virtually isotropic.Homogeneous spherical packs appear to be hydraulically isotropic at their highest porosity and their anisotropy ratio slightly increases when they become more heterogeneous.
The results of the simulation experiments indicate that particle shape and preferential alignment influence hydraulic anisotropy.Packings with the same degree of heterogeneity have higher anisotropy at higher aspect ratio.It can therefore be concluded that, as with geological history, the placement method influences hydraulic anisotropy.

Figure 1
Figure 1 CAD models of (a) spherical homogeneous grain pack (b) spherical but heterogeneous grain pack (c) aspherical homogeneous grain pack (d) aspherical and heterogeneous grain pack.

156Figure 3 AFigure 4
Figure 3 A geometric representation of the model set-up showing four no-flow and two Dirichlet (inflow and outflow) boundaries. .

Figure 5 Figure 6
Figure 5 Typical finite element meshes of the CAD models shown in Fig. 1.Hybrid mesh consisting of tetrahedral, hexahedral, prism and pyramid elements are used and number of elements ranges between (1.4 − 1.8) × 10 6 .

Figure 7
Figure 7 (a) Pressure and (b) velocity fields within the pore space of the geometry (c) Zoom into the velocity field of a region in (b).

Figure 9
Figure 9  Plots of anisotropy ratio A r versus grain particles aspect ratio R for some of the measured samples.

TABLE I
Typical anisotropy ratio (A r ) of synthetic granular packs.K h and K v are the horizontal and vertical permeabilities in Darcy (1 Darcy ≈ 1 × 10 −12 m 2 ), respectively.Several other A r of packs A, B, C, E, F, G and H are plotted in Figure8.