This study is applicable to sandstones reservoirs only. The data used to create each geological model is not originally from a specific oilfield; instead, a dataset of virtual rock samples is created using e-Core technology (Numerical Rocks, Trondheim, Norway) from the geological model. These models are based on outcrop models such as Fontainebleau, Bentheim and Berea. The geological modeling procedure uses a random porosity within a range given by equation 3.1. Equations 3.1 and 3.2 show porosity and permeability ranges from samples used in the work of Smits and Jing [2] and Cense [1].
The present work data set contains 34 digital rock models, subdivided in 19 Fontainebleau models, 11 Bentheim models, 2 Berea and 2 North Sea reservoir models.
Fontainebleau sandstone consists of monocrystalline quartz grains that have been eroded and transported for long periods before being deposited in dunes along the sea shore during the Oligocene, roughly 30 million years ago. The rock is well sorted containing grains of around 200 µm in diameter. During its geological evolution, the sand was cemented by silica crystallizing around the grains. Fontainebleau sandstone exhibits intergranular porosity ranging from 0.03 to roughly 0.3 [13] the permeability varies widely, which is mainly caused by the variation of volume of quartz cement overgrowth. Depending on the weathering process of the sandstone, kaolinite, illite and smectite can be found in Fontainebleau sandstone (Thiery et al. [48]).
Bentheim, is a water wet sandstone, well sorted in fine to medium coarse grains with subrounded to rounded angularity and mainly composed of quartz (70-80%), feldspar (20-25%) and authigenic clays (2-3%) [13]. The weathered feldspar and the quartz overgrowth are related; feldspar weathers to kaolinite and often to quartz, illite and water. The formed quartz precipitates over the quartz grains and yields the overgrowth structure. The weathering feldspar can produce secondary porosity; heavily weathered feldspars yield an intra-granular porosity, whereas totally weathered out feldspar results in secondary porosity.
Berea sandstone is a sedimentary rock that generally presents angular grains instead of rounded and relatively fine grain size. This stone has been used for different purposes in construction, mainly as a building rock.
The types of sandstone digital rocks mentioned above were reconstructed previously by Øren et at. [13, 14] using petrographical information obtained from image analysis of 2D thin sections of sandstone. The end results showed a fair agreement with experimental data from the predicted primary drainage and waterflood relative permeability.
The digital rocks type used in this work were built by means of the e-Core technology (Numerical Rocks, Trondheim, Norway), which offers a set of different model types, e.g. Fontainebleau, Bentheim, Berea, laminated model (refer to appendix A). In table 3.1 can be seen the digital rock models created and their respective properties.
Table 3.1 Digital rock models created and their properties.
Below, figure 4.1 show the cummulative frequency distribution of the properties obtained from the geological models created that were used to generate the pore network models.
b)
c) d)
Figure 3.1. The cummulative frequency distribution obtained from the 34 geological models separate in a) Porosity (%), b) Clay (%), c) Quartz (%) and Feldspar (%).
An existing grain size distribution can be either predefined or configured by the users. The examples of predefined grain size distribution correspond to rocks that were used in this work, such as Berea, Bentheim and Fontainebleau. All of them are based on predefined distributions extracted from thin section data. The other alternative is to customize the grain size distribution of a model in which the user could provide minimum, mean and maximum grain size and their respective standard deviation as an input, the reservoir rock properties is shown in table 2.1.
Figure 3.2 Grain size distributions for Fontainebleau, Bentheim, Berea and North Sea reservoir models.
Figure 4.2 shows relative normal frequency distribution with respect to the grain size for a Fontainebleau, Bentheim, Berea and North Sea reservoir geological models. Bentheim distribution accounts for the biggest grain size of the models whereas the North sea reservoir accounts for the smallest grain size of the models. Notice that dashed lines observed in the left side of the Bentheim and North Sea Reservoir curves have not been taking into account for the creation of their geological models. The reason is not fully understood but the fact is that e-Core software truncated the curves.
3.2. Basic Assumptions
In this section the geological model, pore network model, simulations as well as quality control and verification of the model are discussed. Furthermore assumptions and considerations will be explained.
3.2. 1. Digital pore scale rock modeling
The geological model should not only represent the physical domain of interest but also; avoid the variation of macroscopic properties such as permeability and porosity with the scale.
The representative elementary volume (REV), should be large enough that fluctuations on a molecular level are not taken into account but small enough that the properties of the porous medium can be described and not change significantly with the scale [49].
In order to select the proper REV, an autocorrelation function a could be deployed. The autocorrelation function is known as a classic statistical tool for rocks "microgeometry" and is defined as a characteristic length scale of the porous rock which provides semi-quantitative information about the overall sample textures. It was demonstrated that for absolute permeability, the REV should be larger than 10 autocorrelation lengths and for relative permeability the REV should be larger than 20 autocorrelation lengths [50].
As well the error made selecting the proper REV for relative permeability and absolute permeability have to be taking into account and will not be significant when the grid spacing d, defined as the side length of the cube model, is smaller than a over ten [50]. The parameter a is the correlation length that can be calculated fitting the two-point correlation function. The definition of the two-point correlation is viewed as a characteristic function dichotomizing a rock image into pore and grain, with an exponential function M(u).
3.2.2. Pore network modeling
Highly irregular geometry of the pore bodies and pore throats are replaced by an equivalent idealized geometry, such as cubic, cylindrical and irregular triangles shapes. Those idealized geometries lead to mathematical tractable problems but hopefully retains the essential features of the pore space which are relevant to the flow [14].
Pore throats are modeled as long narrow channels whose cross-section shape is determined by the dimensionless shape factor G. In the same way, the cross-section of the pore bodies is determined using the shape factor.
Microporosity of clays is usually not properly quantified even with more advance scanning techniques. If clay is present in the model, microporosity needs to be included. In this work a microporosity of 0.5 is assumed. Empirical studies [13] have shown that 0.5 is an average value for microporosity in the most common clay minerals. Therefore, in the model, every voxel of clay is associated with a porosity of 0.5 and is furthermore assumed to be filled with water and not available for flow. The total porosity of the pore network model is given by:
where is the total porosity, is the intergranular porosity, is the total fraction of clay and is the clay microporosity.
3.2.3. Network flow modeling
Two phase flow in the pore network is modeled taking into account the following assumptions and considerations:
The flow is laminar and all fluids are Newtonian, incompressible and immiscible.
Fluid contents in a interface node change with time.
Fluid pressures are only defined in pore bodies.
For transport properties it is assumed that the flow is dominated by capillary forces, therefore viscous forces are insignificant. Thus, the capillary Number Nc should be smaller than 10-5
where µ is the dynamic viscosity in Pa.s, u is the Darcy flow velocity in m/s and σ is the interfacial tension in N/m.
Fluid density of wetting value , non-wetting phase value as well as interfacial tension are kept constant.
Predictive modeling of reservoir rocks requires an accurate characterization of wettability. Based on the triangular pore model for oil and water (refer to theory section) wettability fluid configuration is assumed. In order to characterize wettability effects, the contact angle, oil-wet pore fraction and wettability index are used as follows
A) The wettability contact angles for drainage and imbibition are divided in three different conditions, given by equations 3.6, 3.7 and 3.8. Note that the receding contact angle was kept between zero to thirty for all of them.
0< θr<30, wetting fluid: 30<θa<60, non wetting fluid: 100< θa<150 (3.6)
0< θr<30, wetting fluid: 60<θa<75, non wetting fluid: 105< θa<120 (3.7)
0< θr<30, wetting fluid: 30<θa<75, non wetting fluid: 105< θa<180 (3.8)
where θr is the receding contact angle and θa is the advancing contact angle. Both vary randomly within the limits.
B) The oil wet pore fraction is varied from zero to one while the spatial distribution of these pores was kept constant at 2. This assumption was put forward by Øren et al. [14] when oil recoveries obtained from their work were shown to be fairly insensitive to the way the oil wet pore fraction was distributed.
where is the oil wet pore fraction.
C) Wettability index Amot-Harvey [21] and USBM [22] are used and calculated during the network flow modeling. The Amott index is based on the amount of spontaneous imbibition of a certain phase. For water, the Amott index Iw is defined as
where Sspw is the water spontaneous imbibition, Scw is the connate water saturation and Sor is the residual oil saturation. Similarly, but a reverse operation with oil is carried out; the sample is brought to residual oil saturation state, where two volumes are produced. The first volume occur when drainage spontaneously displaces the water whereas the force drainage (flooding or centrifuging) produces the rest of the volume, the Amott index for oil Io is defined as
where Sspw is the oil spontaneous imbibition. For an extremely water-wet system Io will be zero, while for an oil-wet system Iw equals zero. Clearly, the shape of the capillary pressure curve is not taken into account in the Amott index.
The value of IAH by definition is between -1 and 1. A value of 1 corresponds to the case of perfect water wettability and -1 for oil.
The normalized Amott-Harvey index, goes from 0 to 1, is obtained by
The USBM (U.S. Bureau of Mines) method [22] is based on the observation that in order to increase the saturation in non-wetting fluid of a porous space, energy must be put into the system (conversely, an increase in wetting fluid releases energy). This observation is clear in case of strong wettability. This observation can be generalized by measuring the energy required (area between the curve and the x-axis) to increase the saturation of each fluid successively and estimate that the ratio of these areas is an indication of wettability.
Figure 3.3 Imbibition and drainage capillary pressure curve, required to determine the Amott and USBM wettability indices. Reference [2].
The area Aw between the secondary drainage curve and the x-axis (at +Pc) represents the energy required for the oil to penetrate into the system. In the same way, the area Ao between the primary imbibition curve and the x-axis (at -Pc) represents the energy require for the water to penetrate the system.
For an extremely water-wet system USBM is very large and positive, for an intermediate-wet USBM lies around zero and for an extremely oil-wet system USBM will be very large and negative.
The USBM index can be normalized in order to have the same range as the Amott-Havey Index from -1 to 1.
The wettability index USBM after previous normalization, can be again normalized using equation 3.14 where the index goes from 0 to 1. The conversion is given by
D) Wettability classification based on the average difference between IUSBM and IAH relationship is adopted in this work. This classification was introduced by Dixit et al. [51] using analytical relationships that later Skauge et al. [52] reported some experimental confirmation of wettability classes, where a set of cores from sandstone reservoirs of North Sea were grouped in three subclasses, fractionally-wet (FW), mixed-wet large (MWL), and mixed-wet small (MWS). Considering the first one (FW) in the wet-state where the oil-and-water-wet sites are random with respect to size and the other two (MWL and MWS) in the mixed-wet state water-and-oil-wet pores are sorted by size.
The subclasses are distinguished from the following analytical relationships derived between IUSBM and IAH:
where, Rmax and Rmin are radius of pore size distribution and v is the volume exponent. The theoretical line of FW is as well defined as :
From this point onwards, wettability will be referred to as wettability condition 1, 2 and 3.
Another condition established over the network is the boundary conditions; periodic boundary conditions are imposed across the planes perpendicular to the direction of the applied pressure gradient.
Hysteresis value is assumed to be constant and equal to 20% or 0.2.
3.2.4. Statistical quality control & Verification of the model
Statistical quality control
One important step of the work is to verify that the digital pore scale model is properly built before extracting the pore network which will predict constitutive relations. This procedure is done after the geological model is completed.
The aim is to use geometrical quantities to characterize and distinguish statistically the microstructures. The goal is to generate a 3D configuration that matches stochastic functions of a reference medium.
The statistical tools introduced within theory section 3, such as: two-point correlation function, Lineal path function, local porosity probability and local percolation probability are used to characterize the medium (microstructure) and verify if the virtual rock created complies with the minimum requirements for the network extraction. This fact was first demonstrated for two dimensional images [25].
Verification of the model
Once absolute permeability from the pore network is known, it is considered that for the pore network model to be representative for the geological model, their ratio should be close to 1 or no more than 2. The relation is given by,
where is the absolute permeability from the geological model and is the absolute permeability from the pore network model.
The difference between total porosity from the digital pore scale model (geological model) and the porosity from the pore network should not be larger than 3 %. This values has been observed in order during several works.
3.3. Reconstructed Fontainebleau model example
3.3.1. Geological model (Voxelized 3-D model)
The reconstruction procedure consists of three main modeling steps: sedimentation, compaction and diagenesis [12].
Sedimentation
The first step consists of defining the model size and resolution. The size of the reconstructed model is 8003 voxels with a resolution of 7 µm. The process itself commences with the generation of a grain-size distribution curve, is predefined in e-Core for Fontainebleau, Bentheim and Berea (meaning there is no need to specify parameters such as mean, maximum and minimum grain size with their respective standard deviation). Then, each grain is dropped into a bounding box forming the sandbed, using a sequential deposition algorithm [12].
This process continues in a low energy environment or in a high energy environment. In this work the latter was selected, where the sandbed is usually influenced by lateral forces like streams and waves as well as gravity, placing the sand grains at the lowest available position (global minimum).
Figure 3.4. 2-D section after sedimentation process is completed. Grains forming the sandbed are shown.
Compaction
This process is not used for the Fontainebleau model but is covered in details within the theory section of this report.
Diagenesis
Diagenetic processes in real reservoirs are commonly rather complex, difficult to model. They can completely change the heterogeneities of the primary sandstone structure. Therefore, few diagenetic processes are modeled, in this case, quartz cement growth and authigenic clay growth.
To model the quartz cement growth process, two options are available in the software, the algorithm described by Schwatz and Kimminau [24] called voronoi method and the Pilloti's method [23], which is more frequently used due to good in-house experience.
Table 3.2 shows a set of parameters require by Pilloti's technique and clay content used in the present model.
Table 3.2. Parameters used for Pilotti's technique (all these parameters are user defined): target φ is the desired porosity after the cementation process is completed, cutting planes (used from 15-30 because Fontainebleau is well rounded), sphericity factor, deviation of a particle from a sphere and clay content, clay pore filling and clay pore bridging.
The target porosity φ is used as termination criteria for either the Pilotti processing or voronoi cementation. This value will be reached through a binary search for the inner radius.
Figure 3.5. 2 D section after diagenesis process is completed. Porosity showed in color blue, quartz in color yellow and clays in brown.
After the completion of the three previous modeling steps, the results are shown within the table 3.3
Table 3.3. Results from the geological model. The porosity, clay content and quartz content are the results of the diagenetic process. Total φ is total porosity is calculated from equation 3.4.
After the geological model is completed, absolute permeability and formation factor are calculated, refer to table 3.4. Absolute (single phase) permeability is calculated using lattice Bolztmann simulation on the voxelized binary image, derived from fluid velocity field. The average absolute permeability is defined as the average of the directional dependent absolute permeabilities whereas the formation factor is determined from the electrical potential solving the Laplace equation.
Table 3.4. Absolute permeabilities kx, ky, kz and formation factors, Fx, Fy,,Fz calculated in the x-, y- and z- direction respectively.
3.3.2. Statistical Verification
From the two point correlation function, the autocorrelation length, a , is 100 µm. This value is obtained when the two-point correlation function reaches zero. In figure 3.5 (a) some fluctuations are observed in the tail of the function which can be related to the presence of anisotropy in the microstructure. However, nothing more can be said because this function does not contain connectedness information [50]. It is demonstrated that for both absolute permeability and relative permeability, the REV is representative because the value is greater than 10a and 20a.
Lineal path curves, as can be seen in figure 3.6 (a), show the probability of finding a line segment of 260 µm entirely in the pore space.
Figure 3.5 (b) shows the 3D local porosity distribution µ(φ,L) [53]. The µ(φ,L) measures the empirical probability of finding the local porosity φ in a cubic cell of linear dimensions L. A measuring box with side length L = 220 µm is set as an input to find the local porosity distribution curve. The width of the distribution indicates the fluctuation in local porosities (variance in porosity) and the most probable peak, situated around 0.23-0.24. The most probable peak is in good agreement with the value of total porosity obtained in table 3.3.. During the generation of the geological model, unfortunately the characterization of the sample in terms of µ(φ,L) does not allow to quantify the connectivity of the model. Transport properties depend critically on the connectivity of the pore space.
The local percolation curve shown in figure 3.6 (b) represents the connectivity of the pore microstructure in a 2D section. The model is globally connected in 2 directions, however, does not imply that they have a similar connectivity; in fact, x-direction has a better connectivity than y-direction which corresponds to the results obtained from the two-point correlation curve.
In general the microstructure shows a well connected pore structure with slight differences between the x and y direction (anisotropy). Wide range of local porosities distribution is present in which the predominant porosity circa around 0.23-0.24.
b)
Figure 3.6. Two-point correlation function (a) and a 3D Local Porosity distribution µ(φ,L=220 µm) (b)
b)
Figure 3.7. Lineal path probability function(a) ,local percolation probability function in 2D (b)
3.3.3. Pore network modeling
Pore network extraction comes from the inter-granular pore space that is defined by sedimentary grains and quartz cement. The procedure is called ultimate dilation of the grains and is explained by Øren and Bakke [12]. The skeleton is used as a basis for mapping pore throats (links) and pore bodies (nodes), please refer to theory section for further details.
An important parameter which influences the behavior of network models is the coordination number and is defined as a number of links that are connected to the nodes [12]. It is shown in the table 4.5 as connectivity. Every pore body and throat is described in terms of dimensionless shape factor defined as where A is the average area and P perimeter length. This factor approximates the irregular pore bodies or throats by an equivalent cubic, irregular triangular or cylindrical shape [14].
Table 3.5. Data obtained from the pore network extraction from the Fontainebleau virtual rock. Total porosity (clay fraction and total network), average connectivity of the network, nodes and links radius divided in maximum, minimum and mean with their respective total number.
Figure 3.8. Measured pore body (left side) and throats (right side) size distributions for the reconstructed Fontainebleau. Mean values can be observed.
3.3.4. Transport Properties from pore network
Three flood cycles are carried out within the pore network, i.e. primary drainage, primary imbibition and secondary drainage. Considering a water-wet model sample initially filled with brine, the following processes are carried out:
During primary drainage, oil is injected to the model sample until connate water saturation is reached.
Waterflooding process starts when water displaces oil until it reached irreducible oil saturation of the model. During this process two cycles could be present, i) spontaneous water imbibition occurred when the capillary pressure reached zero, therefore allowing water to spontaneously imbibe into water-wet pores in the network. ii) Force water drive occurred when the capillary pressure further decreased to negative to enable waterflooding oil-wet pores.
Finally, during the secondary drainage oil displaces brine until the connate water saturation is obtained. During the process two cycles could be present, i) spontaneous oil imbibition occurred when the capillary pressure increase from negative values up to zero. ii) Force oil drive occurred when capillary pressure is increased from zero until it reached the highest value obtained during primary flooding in which the water is displace from water-wet pores.
It is assumed that oil-water displacements are mainly governed by capillary pressure. Variation in the advancing and receding contact angles are taking into account in the simulation as well as various pore level mechanism such as trapping of wetting and non-wetting phase by snap-off and bypassing, film flow via wetting phase.
Taking into account assumptions and consideration for fluid flow, simulations results can be visualized in the table 3.6. Properties depending on saturation can be calculated by solving Navier Stokes equation [54].
Table 3.6. Fluid Flow results from the pore network model.
3.3.5. Verification of the model
To verify the correct reconstruction of the model, the results have to be in tolerance or within the ranges specified previously. Note that in table 3.7, the porosity difference, established from equation 3.22 shows good agreement, as well as absolute permeability ratio, showing a value of less than two. Also, the formation factor difference is relatively small close to 4. The formation factor verification will be discussed further.
Table 3.7. Results obtain after verification of the model, where porosity difference, permeability ratio are based on equation 3.21 and 3.22
The pore network model was verified to be representative of the geological model which was characterized by statistical tools.
Pore Network data
A total of 34 different pore network models were extracted from every geological model described on table 3.1. The following table contains the average, maximum and minimum values of each category. Note that one of the pore networks is not shown due to confidentiality protection.
Table 3.8. Pore network data per each type of geological rock model.
Analysis & discussion
During the reconstruction of the geological model of the Fontainebleau sandstone, some constrains were found using a model resolution of 6 µm. It was found, please refer to figure 3.9, that when target porosity was set below 15%, a porosity overestimation was obtained from the geological model with respect to the target porosity (desire porosity from the geological model without taking into account clays). In consequence, after the pore network is extracted, the porosity of the network will present an overestimated porosity in comparison from the desire porosity set for the geological model.
Figure 3.9 Porosity overestimation from diagenesis process (end results of the geological model) with respect to the target porosity (desire geological model porosity).
Another constrain or limitation, reader refer to figure 3.10, was found by the time the ratio between the absolute permeability from the geological model and from the pore network model was desired to stay in between 0.5 and 2 (Eq.3.21). Using a constant resolution of 6 µm for the model was not possible to keep the ratio within limits. Therefore, when porosity was desired to be more than 25%, the ratio started to be greater than two, on the other hand, when decreasing porosity below 15% the ratio started to be less than one. Figure 3.10. shows how varying the resolution of the model at different desired porosities, the value of the absolute permeability ratio is obtained.
Figure 3.10. Changes in resolution and their effect on absolute permeability ratios at different desire porosities.
The data shown in figures 3.9 and 3.10 are related, for porosities lower than 15%, some details in the microstructure were missing and therefore the resolution was needed to decrease. On the other hand , for porosities greater than 25% too much details from the microstructure was taken into account which effect the absolute permeability ratio negatively.
Sensitivity Analysis (from voxalized model)
Figure 3.11 shows a qualitative comparison, generated from maximum and minimum cases of different parameters such as: clays (pore filling and pore bidging), grain sphericity, compaction factor (even thought was not used for Fontainebleau model), Cluster algorithm for pore filling clays, spherical type of grains and depositional energy environment and their effect on the difference in model porosity and absolute permeability ratio. The sensitivity analysis is not an exact quantitative result and is not intended to be taken as such. Instead, this result could be used as a qualitative reference to build a Fontainebleau model. Figure 4.10 clearly indicates the effect of pore filling clays in the model either for absolute permeability ratio or porosity differences. Other parameters that influence permeability and porosity are displayed in the figure 3.11.
Figure 3.11 Sensitivity analysis, left side with respect to absolute permeability ratio and the right side with respect to differences in porosity
Formation factor relationship
No literature was found related to how the formation factor could be verified in order to be representative from the pore network model. Therefore, from all the data obtained from different models a relationship was found. It is shown in figure 3.12 how the formation factor presented as a ratio, defined as the formation factor calculated from the geological model and the formation factor calculated from the pore network, stayed between 1 and 2. Note that this relationship was found after the previous verification of the model for either porosity or permeability was done.
A linear regression with a negative slope was found when plotted with porosity. For higher porosity, the ratio tends to be closer to one and at lower porosities, the differences between the geological model formation factor and the pore network model formation factor increase, thus the ratio increase.
.
Figure 3.12 Formation Factor ratio vs Porosity obtained from the geological model
