When the permeability field of a given porous medium domain is heterogeneous by the existence of randomly distributed fractures such that numerical investigation becomes cumbersome, another level of upscaling may be required. That is such complex permeability field could be relaxed (i.e., smoothed) by constructing an effective permeability field. The effective permeability field is an approximation to the real permeability field that preserves certain quantities and provides an overall acceptable description of the flow field. In this work, the effective permeability for a fractured rock system is obtained for different coarsening scenarios starting from very coarse mesh all the way towards the fine mesh simulation. In all these scenarios, the effective permeability as well as the pressure at each cell is obtained. The total flux at the exit boundary is calculated in all these cases, and very good agreement is obtained.

1. Introduction

The first level of upscaling encountered when studying phenomena occurring in porous media is related to adopting the continuum hypothesis to such complex formations. That is moving from pore scale simulations in which variables are only defined within the fluid region to the continuum description in which macroscopic variables are defined everywhere in the simulation domain. In doing so, certain conditions and length scale constraints need to be fulfilled in order to generate macroscopic quantities that are scale independent [1]. While these length scale constraints may be easily adopted in statistically homogeneous porous media, it may be difficult to satisfy particularly closer to the interface boundaries between media with very different textures [2]. Unfortunately, these sharp interface boundaries are ubiquitous in rock systems which apparently contain several fractures and fissures. These fractures have significant effects on both flow and transport in porous media. They provide fluids and solutes with significantly less resistance paths compared with the surrounding matrix and therefore can lead solutes to reach further distances than what simulations ignoring them would normally predict. There are two approaches to dealing with these fractures depending on their density and distributions in rock systems. For rock systems with few penetrating fractures, it is customary to consider the discrete fracture model (e.g., [3, 4]). In this model, fractures are assumed as two-dimensional surfaces impeded within the given rock matrix domain. These fracture surfaces have their own properties and equations which include interaction terms between the rock matrix and the fractures. When the numbers of the fractures are large and their geometrical distributions are complex such that it becomes exceedingly difficult to consider them in simulation, another alternative is to construct an equivalent porous medium domain by further increasing the averaging representative volume. In this approach, one abandons the fractures and considers another continuum, coarser in nature, but may be acceptable for our engineering application and design.

In terms of simulation, although the discrete fracture model provides an elegant way to deal with fractures when their density is manageable, the drawback of this approach, however, is the fact that the interaction term between the fracture and the surrounding rock matrix is, in most cases, not known and ends up as a fitting parameter. The other approach is to refine the mesh all the way to the fractures level, which apparently poses enormous difficulties in large systems particularly when modeling complex transport problems like multiphase compositional flows. It has been, therefore, proposed that an even higher level of upscaling may be required to generate field properties to reduce the required computing resources and provide certain requirements satisfaction. These requirements are related to the fact that the solution of the upscaled system should produce approximately close results to those obtained from simulating the real system. That is if the permeability field may be relaxed such that the sharp interface boundaries may be smoothed, in a sense, the requirements for denser mesh closer to the sharp interface boundaries may be abandoned, and it may be possible to significantly decrease the number of degrees of freedom of the system. Different approaches have been considered in the literatures based on either deterministic approaches (e.g., [5–13], etc.) or stochastic approaches (e.g., [14–18], etc.). In the next section we introduce a quick review of the previous work emphasizing on two methods which are related to this work closely.

2. Investigation of the Upscaling Techniques

A review on the upscaling of hydraulic conductivities in heterogeneous media may be found in the work of Farmer [11] and in the work of Wen and Gomez-Hernandez [8]. In both reviews the authors pointed out that the need for upscaling stems from the fact that the scales at which measurements are taken and the scale at which reservoirs are discretized are usually different. Therefore, it is suggested that a transformation of hydraulic conductivities from the scale of the measurements into a coarser grid of block conductivity tensors amenable for input to a numerical flow simulator may be needed. Holden and Nielsen [19] pointed out that the scale at which permeability in a reservoir may change is much finer than what is possible to use in a reservoir simulator. They proposed a methodology in which permeability is formulated as a minimization problem. In the method, the upscaled permeability is defined as that which minimizes the difference between the pressure and the velocity fields generated by the fine and coarse scale pressure equations. In this work, however, we apply an averaging over a coarser grid of the porous media domain such that certain quantities are conserved. The problem may be casted as follows: consider a three-dimensional rectangular domain, Ξ© with boundary πœ•Ξ©, if 𝐯 represents Darcy’s velocity and 𝑝 the pressure field of the fluid saturating the porous medium domain, the governing equations describing this system may be written as βˆ‡β‹…π―=π‘ž,(2.1)𝐯=βˆ’πŠβˆ‡π‘.(2.2) From which one obtains an equation only in the pressure given asβˆ’βˆ‡β‹…πŠβˆ‡π‘=π‘ž.(2.3) The upscaling problem, as described in [11], is to find a corresponding coarse hydraulic conductivity field, ξ‚πŠ such that the new velocity and pressure fields, ̃𝐯,̃𝑝 obtained from solving the equivalent system of equations, are close, in some sense, to the fine scale velocity and pressure fields, 𝐯,𝑝. The system of equations that need to be solved in this case takes the formΜƒΜƒξ‚βˆ‡β‹…π―=Μƒπ‘ž,𝐯=βˆ’πŠβˆ‡Μƒπ‘.(2.4) The general methodology is to seek an equivalent (effective) hydraulic conductivity field ξ‚πŠ, calculated from a set of boundary conditions, and it may be used afterwards to yield accurate results with different boundary conditions. To do so, the rectangular domain is divided into a number of segments in each direction, that is, 𝑁π‘₯, 𝑁𝑦, and 𝑁𝑧 to form boxes of sizes, 𝐻π‘₯, 𝐻𝑦, and 𝐻𝑧. These cells are called the coarse cells. Furthermore, each of the boxes is divided into 𝑛π‘₯, 𝑛𝑦, and 𝑛𝑧 finer boxes with sizes β„Žπ‘₯, β„Žπ‘¦, and β„Žπ‘§. The hydraulic conductivity tensor field 𝐊, generally, takes a different value in each of the fine cells. There have been two approaches highlighted by Farmer [11], which may be used to construct the coarser hydraulic conductivity field, namely, the global (overall) strategy and the local methods. The first technique is performed in two stages. In the first stage one needs to perform one or more fine grid experiment which may cover the whole domain or at least most of it. Three test cases for each scenario are performed, one in each coordinate direction on the fine grid within each of the coarse cells. In the second stage, the fine grid simulations are used to extract the coarse scale permeability, ξ‚πŠ. In the second method, which is sometimes called the local-local or the subdomain method, each coarse cell is considered separately, as shown in Figure 1, and the three independent flow experiments in the three directions are performed. In each of these experiments, no flow conditions are assigned to the four surfaces parallel to the flow direction and constant pressure conditions on the two surfaces normal to the flow direction. If the total flux may be computed through one of the faces that were assigned pressure boundary condition, an equivalent permeability for this cell may be computed [20, 21].

There are, however, certain requirements for the correct upscaling as defined in [1, 2]. For the case of single phase incompressible fluids, it suffices to require that the flux across any area in the domain should be the same for all levels of upscaled description. This requirement ensures a divergence-free flow field in case there is no source/sink term. This condition can be defined asξ€œπ΄Μƒξ€œπ―(𝐫)⋅𝐧𝑑𝐴=𝐴𝐯(𝐫)⋅𝐧𝑑𝐴,(2.5) where ̃𝐯 and 𝐯 are the velocity vectors at different level of upscaling, and 𝐫 is the position vector which scans the area 𝐴. Applying this requirement on the local-local method mentioned earlier in which each coarse cell is considered separately, the above equation reduces toΜƒξ€œπ―(𝐱)⋅𝐧𝐴=𝐴𝐯(𝐫)⋅𝐧𝑑𝐴,(2.6) where ̃𝐯(𝐱) represents the velocity vector at the center of the coarse cell face, and 𝐱 is the position vector which defines the center of the cell face. Likewise for volumetric quantities (e.g., the source/sink term, π‘ž), Salama and Van Geel [1] pointed out thatξ€œπ‘£ξ€œΜƒπ‘ž(𝐫)𝑑𝑣=π‘£π‘ž(𝐫)𝑑𝑣=Μƒπ‘ž(𝐱)𝑣.(2.7)

3. The Subdomain Technique

We apply the local-local method to obtain the effective permeability field for a given complex permeability field in fractured rock system. The simulation domain represents a rectangular two-dimensional area with a set of randomly generated fractures. We use the cell-centered finite difference scheme (CCFD) to approximate the governing equations because it satisfies mass conservation law locally. Furthermore, we consider the fractures oriented horizontally and vertically, as shown schematically in Figure 2. The local-local method deals with each macroscopic cell separately. For this purpose, we define two different meshes; the first mesh is the coarse mesh and it is defined globally over the whole simulation domain, and the second mesh is fine according the complexity of the permeability field within each coarse cell. That is we construct the coarse rectangular mesh over the whole simulation domain, {𝑋0,𝑋1,𝑋2,…,𝑋𝑖,…,π‘‹π‘š}Γ—{π‘Œ0,π‘Œ1,π‘Œ2,…,π‘Œπ‘—,…,π‘Œπ‘›}, and for each coarse mesh we construct a finer mesh given as {π‘₯0,π‘₯1,π‘₯2,…,π‘₯𝑖,…,π‘₯𝑝}Γ—{𝑦0,𝑦1,𝑦2,…,𝑦𝑗,…,π‘¦π‘ž}.

4. Discretization of the Governing Equations

For the generic cell shown in Figure 3, the governing equations for this problem may be approximated using the cell-centered finite difference scheme with the various velocity components at the mid edges of the cell obtained as follows: 𝑒π‘₯ξ‚€1𝑖,𝑗+2=βˆ’πΎπ‘₯π‘₯ξ‚€1𝑖,𝑗+2𝑃(𝑖+1/2,𝑗+1/2)βˆ’π‘ƒ(π‘–βˆ’1/2,𝑗+1/2),𝑒π‘₯(𝑖+1/2)βˆ’π‘₯(π‘–βˆ’1/2)(4.1)π‘₯ξ‚€1𝑖+1,𝑗+2=βˆ’πΎπ‘₯π‘₯ξ‚€1𝑖+1,𝑗+2𝑃(𝑖+3/2,𝑗+1/2)βˆ’π‘ƒ(𝑖+1/2,𝑗+1/2)𝑒π‘₯(𝑖+3/2)βˆ’π‘₯(𝑖+1/2),(4.2)𝑦1𝑖+2,𝑗=βˆ’πΎπ‘¦π‘¦ξ‚€1𝑖+2,𝑗𝑃(𝑖+1/2,𝑗+1/2)βˆ’π‘ƒ(𝑖+1/2,π‘—βˆ’1/2)𝑒𝑦(𝑗+1/2)βˆ’π‘¦(π‘—βˆ’1/2),(4.3)𝑦1𝑖+2,𝑗+1=βˆ’πΎπ‘¦π‘¦ξ‚€1𝑖+2,𝑗+1𝑃(𝑖+1/2,𝑗+3/2)βˆ’π‘ƒ(𝑖+1/2,𝑗+1/2).𝑦(𝑗+3/2)βˆ’π‘¦(𝑗+1/2)(4.4) Equation (2.2) is discretized, likewise, in the following form:𝑒π‘₯(𝑖+1,𝑗+1/2)βˆ’π‘’π‘₯(𝑖,𝑗+1/2)+𝑒π‘₯(𝑖+1)βˆ’π‘₯(𝑖)𝑦(𝑖+1/2,𝑗+1)βˆ’π‘’π‘¦(𝑖+1/2,𝑗)ξ‚€1𝑦(𝑗+1)βˆ’π‘¦(𝑗)=π‘žπ‘–+21,𝑗+2.(4.5) Now substituting (2.7)–(4.3) into (4.4), one obtains an equation in the pressure only, as explained earlier, which may be solved to obtain the pressure field from which the velocity field may be determined. The permeability at cell edges is obtained based on the harmonic mean relationship which is given asπ‘˜π‘₯π‘₯ξ‚€1𝑖,𝑗+2=π‘₯(𝑖+1)βˆ’π‘₯(π‘–βˆ’1)ξ€·π‘˜(π‘₯(𝑖)βˆ’π‘₯(π‘–βˆ’1))/π‘₯π‘₯ξ€Έξ€·π‘˜(π‘–βˆ’1/2,𝑗+1/2)+(π‘₯(𝑖+1)βˆ’π‘₯(𝑖))/π‘₯π‘₯ξ€Έ,π‘˜(𝑖+1/2,𝑗+1/2)𝑦𝑦1𝑖+2=,𝑗𝑦(𝑖+1)βˆ’π‘¦(𝑖)ξ€·π‘˜(𝑦(𝑖)βˆ’π‘¦(π‘–βˆ’1))/π‘¦π‘¦ξ€Έξ€·π‘˜(𝑖+1/2,π‘—βˆ’1/2)+(𝑦(𝑖+1)βˆ’π‘¦(𝑖))/𝑦𝑦.(𝑖+1/2,𝑗+1/2)(4.6)

5. Application of the Local Upscaling Technique to Fractured Rock System

In this work we consider the local-local upscaling technique to generate an effective permeability field which would require less computational load and in the same time would result in reasonable match with the fine scale modeling. As explained earlier the scheme follows the following steps.(i)Discretize the whole domain with the appropriate fine mesh based on the complexity of the permeability field and fractures structure.(ii)Discretize the whole domain with the required coarse mesh that conforms with the fine mesh (i.e., the mesh coincides with the fine mesh).(iii)For each coarse cell, solve the boundary value problem given by (2.3) over the finer cells subject to predefined boundary conditions. That is to get the effective permeability in the π‘₯ and 𝑦-directions, 𝐾π‘₯π‘₯ and 𝐾𝑦𝑦, we consider two cases as shown in Figure 4. For example, in Figure 4(a), we consider a no-flow boundary condition at the top and bottom boundaries and prescribed pressure boundary condition (say 𝑝1=1 and 𝑝2=0) at the other two sides. Similar boundary conditions are suggested in Figure 4(b). The velocity at the finer scale could then be calculated.(iv)Using (2.6), one may calculate the π‘₯ and 𝑦-velocity at the center of the faces of the coarse cell which encompass the finer mesh.(v)These velocities can then be used to calculate the effective permeability components in π‘₯ and 𝑦 directions with the Darcy’s equation (2.2).

6. Numerical Examples

Several scenarios with different fractures numbers and distributions have been considered. The simulation domain in all these scenarios is [0, 0.6] Γ— [0, 0.6]. For the sake of simplicity, the width and the length of the fractures in all the simulation scenarios are considered identical as 0.00015 and 0.12, respectively (i.e., 0.2 times the length of the domain). In the 2-fracture scenario, the length of the fractures is chosen 0.3, 0.5 times of the size of the domain. The locations of the fractures are random. The flow rate in the fracture is usually described with the cubic law [22] in which the flux through the fracture is given as a cubic function of the fracture aperture, therefore𝑄=βˆ’πœŒπ‘”π‘312πœ‡Ξ”β„ŽΞ”πΏ,(6.1) where Q is the flux across the fracture section, 𝜌 is the density, πœ‡ is the viscosity, 𝑔 is the gravity, and 𝑏 is the fracture aperture. Generally, the hydraulic conductivity in the fracture region will be much higher than that in the matrix and plays a dominant role in the flow in the fracture-matrix system. In the paper, the hydraulic conductivity of the fracture region is specified as 3Γ—105 times that of the matrix. The mesh size of the fine mesh is 0.0006 for the matrix and 0.00015 for the fractures. Four scenarios with different number of fractures are considered. In the first scenario, only one fracture in the horizontal direction and one fracture in the vertical direction, respectively, are assumed. The other three scenarios consider, respectively, 5, 20 and 400 fractures in both the horizontal and vertical directions. At first, fine meshes are generated with the number of the mesh cells in horizontal and vertical directions are respectively as follows: 1002Γ—1002 for the 2-fracture scenario, 1008Γ—1008 for the 10-fracture scenario, 1032Γ—1035 for the 40-fracture scenario, and 1480Γ—1500 for the 800-fracture scenario. The flow equations are solved with the fine mesh, and the pressure distribution is obtained for each scenario. As will be seen later, the flow in the whole domain is very much affected by the fractures. To consider different levels of upscaling, the domain is discretized with a coarse mesh of different resolutions. In this work we consider very coarse resolutions and moderately fine resolutions. In all these levels of resolutions, effective permeability is obtained based on (2.2).

6.1. Example 1: One Fracture in Both the Horizontal and Vertical Directions

In this scenario, two fractures are generated randomly as shown in Figure 5(a). The pressure field for this system obtained by solving the governing equations using the fine mesh is shown in Figure 5(b), from which we can see the apparent influence of the fractures in distorting the pressure field.

Several coarsening scenarios are considered, namely, 2Γ—2, 5Γ—5, 20Γ—20, 40Γ—40, 100Γ—100, and 500Γ—500. In all these scenarios, the requirement that the flux across any surface calculated using the fine mesh and any other coarse mesh is the same. The generated effective hydraulic conductivity of each cell, Μƒπ‘˜π‘₯π‘₯ and Μƒπ‘˜π‘¦π‘¦, for different resolution scenarios is shown in Figure 6. Apparently, the permeability of the cells containing fractures is larger. Furthermore, as the resolution of the mesh increases, the permeability contrast sharpens around the fractures. Figure 6 also shows the pressure at the center of each cell for the different resolution scenarios. It is apparent that the pressure field gets closer to the fine scale simulation as the resolution increases. On comparing the flux at the exit boundary (Table 1), it is clear that the flux for different coarsening scenario is close to one another.

6.2. Example 2: Five Fracture in Both the Horizontal and Vertical Directions

In this scenario, five fractures are generated randomly in both the horizontal and vertical directions (i.e., 10 fractures in total), as shown in Figure 7(a). The pressure field for this system as obtained by solving the governing equations using the very fine mesh is shown in Figure 7(b). Again, the influence of the fractures in distorting the pressure field is significant.

Effective permeability in each cell is shown in Figure 8, which shows that the effective permeability approaches that of the fine mesh simulation as the number of grid cells increase. Again the flux at the exit boundary is very close to the different coarsening scenarios as depicted in Table 1.

6.3. Example 3: Twenty Fractures in Both the Horizontal and Vertical Directions

In this scenario, 20 fractures in π‘₯ and 20 in the 𝑦 directions are considered, as shown in Figure 9(a). The pressure field obtained by solving the governing equations based on the actual permeability field is shown in Figure 9(b). The effect of the fractures on the pressure field is apparent.

On the other hand, the effective permeability field obtained from solving the equivalent equations as explained earlier is shown in Figure 10. It is clear that as the number of grid cells increases, the system approach that based on the actual permeability field. Table 1 confirms again that the flux at the exit boundary is almost identical for the different coarsening scenarios.

6.4. Example 4: 400 Fractures in Both the Horizontal and Vertical Directions

As shown in Figures 11 and 12, similar behavior is obtained for this scenario which assumes 400 fractures in both the horizontal and vertical directions for both the effective permeability and the pressure fields. In this scenario, we consider several degrees of coarsening including, 2Γ—2, 10Γ—10, 20Γ—20, 40Γ—40, 300Γ—300, and 500Γ—500.

Apparently in this scenario, the large number of fractures will increase the overall effective hydraulic conductivity field which will result in a significantly larger flux compared with the other scenarios. Indeed this is the case as shown in Table 1.

7. Conclusions

In this work several numerical examples have been considered to calculate effective hydraulic properties for a given fractured porous medium domain. Several scenarios of fractured systems have been considered starting with two fractures up to 800 fractures. Cell-centered finite different method has been used, and therefore the fractures orientations are considered only in the horizontal and vertical directions. Different coarsening scenarios have been considered for each fracture scenario starting with very coarse mesh (2Γ—2) all the way to the fine mesh. The fine mesh used in this work is that associated with the actual permeability field of both the matrix and the rock. Effective permeability for each cell is obtained for each coarsening scenario based on the requirement that the flux at the face of each cell is the same as that considered with the actual permeability field. It is found that the spatial distribution of effective permeability of the medium changes with the change in the degree of coarsening. That is for coarser mesh the effective permeability spans a larger domain while it reduces to the actual permeability field as the mesh approaches the fine mesh. The usefulness of this approach stems from the significantly less computational effort and memory requirement that would be required when using the upscaled permeability field and in the same time the satisfaction of conservation laws at the new scale.


The work is funded by KAUST’s GRP-CF (Global Research Partnership Collaborative Fellows) Program, and KAUST-UTAustin AEA project (7000000058).