DISCRETE FRACTURE NETWORK (DFN) MODELLING OF FRACTURED RESERVOIR BASED ON GEOMECHANICAL RESTORATION, A CASE STUDY IN THE SOUTH OF IRAN

Fractured reservoirs have always been of interest to many researchers because of their complexities and importance in the oil industry. The purpose of the current study is to model the fractured reservoir based on geomechanical restoration. Our target is the Arab Formation reservoir which is composed of seven limestone and dolomite layers, separated by thin anhydrite evaporate rock. First of all, in addition to the intensity, the dip, and the azimuth of the fractures, the magnitude and the direction of the stresses are determined using wireline data e.g. photoelectric absorption factor (PEF), sonic density, neutron porosity, a dipole shear sonic imager (DSI), a formation micro imager (FMI), and a modular formation dynamics tester (MDT). Then, the seismic data are interpreted and the appropriate seismic attributes are selected. One of our extracted attributes was the ant tracking attribute which is used for identifying large-scale fractures. Using this data, fractures and faults can be identified in the areas away from wells in different scales. Subsequently, the initial model of the reservoir is reconstructed. After that, the stress field and the distribution of fractures are obtained using the relationship between the stresses, the strains, and the elastic properties of the existing rocks. The model is finely approved using the azimuth and the intensity of fractures in the test well. Our findings showed that the discrete fracture network (DFN) model using geomechanical restoration was positively correlated with real reservoir conditions. Also, the spatial distribution of fractures was improved in comparison to the deterministic-stochastic DFN.


Introduction
Reservoir characteristics, such as initial porosity, secondary porosity (e.g. fractures), permeability, heterogeneity, anisotropy, fluid properties, fluid saturation, geothermal gradient, and oil recovery mechanisms (Dandekar, 2013;Safari et al., 2017Safari et al., , 2020 are factors that affect the ultimate oil recovery. Fractures as secondary porosity are vital for geological aspects (Pavičić et al., 2017), mining engineering (Elkarmoty et al., 2017;Hajibagherpour et al., 2020), and hydrogeological potential (Pavičić et al., 2017). Research studies have shown that adding geomechanical considerations improves the estimation capability and provides more realistic fault and fracture models (Maerten et al., 1999;Bourne et al., 2000;Mohd Nor et al., 2016). The procedure is started with calculating the stress field using the rock characteristics and a structural model. Then, after modifying the stress components based on time and structural geometry, the final DFN models containing intensity, azimuth, and dip angle of fractures are derived (Maerten, 2010). The final model is useful for drilling aspects (Rad et al., 2021), and production in heterogeneous porous media (Shiri and Hassani, 2021).
Naturally fractured reservoirs, due to their complexity, have complex production and recovery mechanisms. To handle this complexity, dual porosity models (Nie et al., 2012), and triple porosity models (Sang et al., 2016) are developed. Despite their favourable features, they suffer from a lack of realism for modelling naturally fractured reservoirs. The solution could be making deterministic and stochastic discrete fracture network (DFN) models (Fang et al., 2017), using geometrical (Soleimani, 2017) and geomechanical (Maerten, 2010) elastic laws. In numerous studies, the use of geomechanical data for structural modelling shows that these models are quite valuable for estimating natural fractures and faults (Bourne et al., 2000;Bourne and Willemse, 2001;Davatzes et al., 2005). However, numerical models based on continuum mechanics are important tools for describing structural geometry in oil exploration and production. Using a combination of structural and ge-Rudarsko-geološko-naftni zbornik i autori (The Mining-Geology-Petroleum Engineering Bulletin and the authors) ©, 2021, pp. 151-162, DOI: 10.17794/rgn.2021.4.12 omechanical restoration leads to the improvement of the model, which is mainly due to calculating the deformation and the stress distribution over time (Maerten, 2010). This method is based on physical laws and the linear elastic theory that can be replaced with geometrical and dynamic methods (Fenz, 1962). The following studies illustrate the development of this modelling: Martha introduced two-dimensional physical and linear elasticity laws for rock deformation and cross-sectional reconstruction (Martha, 2003). De et al. performed 3D restoration of a normal fault model using the finite element method. In this method, tensile properties of rocks such as Young's modulus and Poisson's ratio were used to move the fault surfaces by coupling the dynamic relaxation method and the finite element method (De et al., 2003). Mace et al. obtained a total strain field by combining the well data and the stress history that was extracted by the discrete element method (Mace et al., 2005). Shackleton performed geomechanical restoration for estimating fractures of the Santa Cornell anticline in the Pyrenees in Spain (Shackleton, 2009). Laurent et al. used geomechanical restoration to determine paleo-geographical parameters and paleo-coordinates of the layers in their depositional time (Laurent et al., 2010). Shukla and Jayakumar performed geomechanical and geometrical restoration for cap rock modelling of a basaltic rift reservoir in India. Following this research, the stress and strain components, the reservoir parameters, the spacing, the intensity, the dip, and the azimuth of fractures were extracted, and the DFN model was constructed based on the statistical distribution of fractures (Shukla and Jayakumar, 2011). Durand-Riard et al. performed 3D geomechanical restoration using geological boundary conditions on strike-slip and dip-slip faults in the Niger Delta (Durand-Riard et al., 2013). Phillips et al. applied geomechanical restoration to find the optimal location of wells in fractured reservoirs. The geomechanical restoration was performed using all types of data, and hence, stress variations and all the relative impacts on fractures were modelled accurately (Phillips et al., 2014).
The deterministic-stochastic DFN model has both large-scale fractures and small-scale fractures. Seismic data is used to identify and determine faults and build a large-scale deterministic DFN model based on the data. Core data and wireline logging data are used to find small-scale fractures. However, they are in the vicinity of wells. Therefore, a stochastic model is used to perturb fracture properties far away from well locations and build a small-scale stochastic DFN model. The final DFN model is a combination of the deterministic and stochastic models (Fang et al., 2017).
The forward geomechanical methods use mechanical concepts and physical laws to explain the formation of faults and fractures. Moreover, the basis of the method involves calculating the stress distribution at the time of fracture generation using all structural data including: faults, fractures and folds, rock characteristics, and tectonic events. Also, secondary activities, such as mineralization and rock inclusions are essential, especially for mining exploration (Mahdavi et al., 2015). The main components of the forward geomechanical restoration model are: (1) the initial condition of paleo-geometry, (2) the constitutive laws of sedimentary units, and (3) boundary constraints. The restoration technique is a deformation mechanism of layers from the present structural and geomechanical characteristics back to their original state in paleo time. In this current study, we applied this methodology to a synthetic model of an Iranian oil field. In the oil field, we increase the reality of fractures of the Arab Formation reservoir by updating the fracture characteristics and their changes through geological time. This formation is a fractured limestone reservoir in the south of Iran which is folded by a vertically upward force of a salt dome during its depositional time. The present study aims to build an improved DFN model of folded and faulted reservoir formation using the finite element method by the forward geomechanical restoration method. The reversibility of the modelling allows us to see past and future conditions in the restoration model. The results showed improvement of fracture modelling by the geomechanical restoration DFN model with respect to the deterministic-stochastic DFN model.

Data and Methodology
In this study, data was collected from an Iranian offshore oil field which is located in the Persian Gulf (X oil field in Figure 1). This oil field was discovered in 1967. It is located around 80 km to the southwest of the Lavan oil processing centre and 40 km southeast of South Pars Field. The most important oil reservoir of the current oil field, which is the target of our study, is the Arab Formation. This formation has seven limestone and dolomite layers that are separated by thin anhydrite evaporite rocks. It is around 150 m in thickness and is saturated by 40.1 API oil. Interpreted FMI logs (secondary porosity indicators, Figure 2) indicated a significant amount of secondary porosity in the Arab Formation. Figure 3 shows the stereo net plot of extracted fault patches in the Arab Formation. The predominant dip azimuth in well data is to the northwest and to the southeast. It indicates that most of the fractures in the field are oriented northwest. Upper formations of Khatiyah, Shuaiba, and Khuff have a mean thickness of 100, 70, and 400 m of argillaceous limestone, chalky limestone, and dolomite-limestone rocks, respectively. 2D seismic data, 97 km 2 high resolution 3D seismic data, and data from 15 wells are used to make the geo-model including the top horizons and fault surface of the Khatiyah, Shuaiba, Arab, and Khuff formations. The well data are leak-off test (LOT), wireline logging as follows: calliper log (CALI), sonic log (DT), gamma ray log (GR), micro spherically fo- cused log (MSFL), shallow resistivity log (LLS), deep resistivity log (LLD), neutron porosity log (NPHI), density log (RHOB), formation micro imager (FMI), and Modular Formation Dynamics Tester (MDT). Aside from the conventional well logging data, which are used for geo-modelling, LOT is used to determine stress field magnitude and orientation, MDT is used to determine pore pressure, and FMI log is used for fracture intensity and orientation.
The studied oil field is located in the Persian Gulf (see Figure 1) and folded by different forces. The major force is exerted vertically upward by a salt dome. It made an anticline during the depositional time of upper layers. The topological map of the investigated formation, the position of normal faults, the oil field geometry, and the salt dome in the seismic survey are shown in Figure 4. Aside from the upward stress of the salt dome, horizontal stresses are also applied by tectonic activities. An analysis of thickness variations in the field indicates moderate and consistent increases in interval thickness in all directions away from the crest of the dome. This shows that deformation was syn-depositional, and the structure is being actively uplifted through depositional time.
To reduce the risk of exploration, drilling, and production of oil reservoirs, it is necessary to build a threedimensional model as close as possible to the reality of petrophysical and production characteristics of oil reservoirs (e.g. porosity, permeability, oil saturation, and pore pressure). The first step in this process is the use of a three-dimensional constructional approach for paleogeometry, and then, the use of the forward geomechanical simulation by specifying the boundaries. The threedimensional constructional model is the definition of structural framework using the fault modelling process. Fault modelling is the process of constructing a representative fault model based on the complex fault interpretations from the seismic data. The next step is the incorporation of stratigraphic surfaces. The objective of horizon and zonation modelling is to generate gridded surfaces based on seismic horizon interpretation and the formation tops.
The equivalent continuous model and the DFN model are two different strategies to study the effects of frac- tures for further fluid flow simulation (see Figure 5). In the equivalent continuous model (see Figure 5(b)), the highest storage capacity belongs to the matrix, and the highest fluid flow is through fractures. In this model, the relationship of fluid flow between matrix and fractures is defined by a transfer function (Dai et al., 2007). The deterministic-stochastic DFN model (see Figure 5(c)) is another way to describe the fracture distribution. In this model, fractures longer than 200 m are generated deterministically by the ant tracking attribute, and fractures shorter than 200 m are generated stochastically based on the spatial distribution of fracture intensity. The final DFN model is dependent on the distribution of fracture intensity, fracture orientation or azimuth, fracture length, and fracture dip angle. The flowchart of this method is shown in Figure 6.
The three-dimensional reconstruction, which is based on the finite element method (FEM), allows all physical laws to analyse the geological strain. In this process, the gradual deformation of various geological structures including faults, joints, folds, beddings, and also the elastic/plastic rock parameters, the internal stress and plastic strain distribution, the volumetric strain, and the plastic shear strain are calculated. Then the FEM is run with moving boundary conditions for transforming the layers to their original horizontal condition. Plastic shear strain from this methodology is responsible for the present fault and fracture distributions.
The present study is conducted on a fractured carbonate reservoir with several faults. The data used in this study are: (1) image logs for dip, azimuth, and intensity of the fractures; (2) sonic and other logs to determine the elastic properties of rocks; (3) 3D post-stack seismic data for reservoir modelling, extraction and integration of seismic attributes for identifying faults and fractures. The geomechanical restoration algorithm is shown in Figure 7. This method uncovers the creation and development of faults as well as showing the hidden faults. The method gives petroleum engineers information on oil migration pathways and on the geometry of the oilfields to reduce risk and to improve decision making (Maerten, 2010). Based on this method, once the rock mass is subjected to stress, if the stress yields the tensile strength of the rock, the strain and displacement of the rock will occur irreversibly. The rock mass is subdivided into small elements and each element has its own properties, and its own reaction toward other elements. The reaction depends on the applied forces, elasticity, displacement, and contact area. The purpose of the restoration is a better interpretation of the geological structures and their development through time.
The forward geomechanical simulation was used to evaluate and construct the structural model from the time of the deposition to the current state of the layers and also, to obtain the structure's geomechanical and lithological properties. According to the algorithm (see Figure 7), the data from wells, beddings, faults, and petrophysical properties along with pore pressure and other data were used to run these analyses.
In geological mechanics-based restoration, rocks are mechanically deformed backward to their paleo geometry.
The fundamental equation of the geomechanical restoration method by combining mass and momentum conservation laws based on FEM formulation is as follows (Belytschko et al., 2014): (1) Where: is the divergence operator in the unrestored space, P is the nominal stress tensor, is the unrestored rock density, b is body force, is acceleration vector. This equation is in the strong form state. For solving this equation by the finite element method, it is necessary to have it in the weak form state. The weak form state of this equation is as follows: (2) After assembling all element matrices to global matrices, it is time to solve the following system of equations by the dynamic relaxation method of finite element method: (3) Where: M is the mass matrix, F int is the internal force vector, F ext is the external force vector, c is a damping factor, is the acceleration vector, is the velocity vector. The dynamic relaxation method is based on explicit time integration. There are two terms of explicit and im- plicit methods for solving the equation for each node of n. In explicit, the solution in time t is only dependent on the previous time steps. However, in implicit, the solution of the equation is only dependent on time t. Explicit time integration is based on the central difference scheme as follows: (4) By inserting these relationships into the above system of the equation we have: (5) model in size and physical properties (see Figure 8(b)). The third synthetic model is a layer similar to the second model with an inclusion with Poisson's ratio of 0.35 and Young's modulus of 100 GPa (see Figure 8(c)).

Image log and stress field model
Core and well log data for mechanical rock properties Conventional well logs 3D Seismic Data: Extraction of seismic attributes, folds, faults and structural interpretation Petrophysical and structural 3D model including folds and faults interpretation, quality control, extraction of elastic parameters from mechanical earth model, linking of major tectonic events to paleo-geometry DFN modeling by geomechanical forward modeling and restoration method, including applied gravity loading, stress-strain relationship constraining boundary conditions, and distribution of plastic strain through geological time scale.
Azimuth, dip, and intensity of fractures In Figure 8, three synthetic symmetrical folds are examined, similar in structure and thickness, but different in detail. The first model consists of three layers, in which the interfaces have neither friction nor spacing between them, but they are always in contact. They are mechanically interconnected and have movement, deformation, and rotational capabilities. The layers are homogeneous, with Poisson's ratio of 0.25 and Young's modulus of 10 GPa (see Figure 8(a)). The second synthetic model is a homogeneous layer similar to the first

Results and Discussion
To analyse the methodology, we applied the finite element method to restore three simple synthetic 3D geo-logical structures. These structures illustrate the impact of contact interfaces and heterogeneity on restoration process. The restoration results of synthetic models are shown in Figure 9. The first synthetic model shows the layers are horizontally stretched to both the left and the right sides, and the deformation is spread between the layers. It causes compression at the top of the layers and tension in the bottom of the layers that leads the lower layers to become narrower than the upper layers (see Figure 9(a)).
In the second model, a homogeneous layer, there are compressive conditions at the top of the layer and tensile conditions at the bottom. In this case, the bottom part of the layer is narrower, and the top part of the layer is wider (see Figure 9(b)).
The purpose of the third model is to investigate the effects of heterogeneity in restoration (see Figure 9(c)). As shown in the figure, inclusions will obey the stress during restoration. It applies high compression at the top of the inclusions and high tension at the bottom of them. The geometry of the restoration model is similar to the second model.
These synthetic models are introduced and examined to show the capability of geomechanical restoration for building stress-strain constitutive laws and for flattening synthetic anticlines. In the following content, we evalu-  ate the strength of this methodology for real data obtained from an Iranian offshore oil field. The 3D seismic reflection data with 600 xlines and 520 inlines from the surface to the time of 2016 ms is used for faults and fracture extraction by the ant tracking algorithm. This data is used for making DFN models. The stereo net of extracted fractures and faults is shown in Figure 3. Interpretation of image logs is an accurate and advanced methodology for identifying and evaluat-ing fracture intensity. These logs were used for determining fracture intensity, fracture orientation, and fracture dip angle in well locations. The validation of fracture intensity by 3D geomechanical restoration model is shown in Figure 10 and the visual outputs of fracture intensity and orientation are shown in Figure 11(b) and Figure 12. Figure 10 shows the cross-plot of fracture intensity between well data and the geomechanical restoration DFN model in the Arab Formation, with R = 0.87. Compared to the results of the deterministic-stochastic DFN model with R=0.72 in the validation well, it showed better correlation with the fracture intensity. Therefore, the 3D geomechanical restoration model was able to reliably predict fracture intensity around the well.
The results of DFN modelling using the geomechanical restoration algorithm on the fractured dolomitic limestone reservoir (the Arab Formation) are shown in Figure 11. The fractures are better estimated (based on correlation coefficients) using the DFN model by the geomechanical restoration (see Figure 11(b)), than the Deterministic-stochastic DFN model (see Figure 11(a)) which is obtained based on the flowchart (see Figure 7). Simulation of fluid flow using this model, due to its high accuracy, will improve the porosity and permeability variations in all directions, and lead to better decisionmaking.
The current methodology is important due to the conservation of the physical laws, such as momentum, mass, and energy that contributed to reservoir modelling. Therefore, the main advantage of this current study is that the laws of physics and linear elasticity theory are replaced with the geometrical and dynamic properties in the methodology for structural reconstruction.
Lack of direct information about paleo structure and paleo data, and using simple constitutive laws are general limitations of geomechanical restoration. Also, Limited FMI and BHI logs data for fracture intensity and orientation are limitations of the current study. No core fracture data is available. No mud loss data is available in the Arab Formation that is often associated with drilling in fractured formations. The rate of mud loss during drilling can therefore be used to indicate the presence of fractures and fracture intensity in a well. The geomechanical restoration model has different advantages. For example, the constitutive models (stressstrain relationships) are known and therefore a complete strain history is available from sedimentation toward final deformation. Also, it has an acceptable match to the current observation and helps to reveal the spatial distributions of porosity, fractures, stress, geometry, and so on. However, as a disadvantage, this methodology suffers from hard criteria for the selection of the deformation constraints or from the flexibility in movement.

Conclusions
The estimation of fractures farther away from wells, based only on the deterministic-stochastic DFN approach, has ambiguities. This limitation is due to the resolution of data for fracture identification. For example, for seismic and well logging data, the minimum resolution for fracture identification should be 100, and 0.1 meters, respectively. The geomechanical restoration technique is a solution that improves the accuracy of fracture spatial distribution. In the present study, this method is applied to the reservoir of an Iranian oilfield. Then, the modelling is compared to the deterministicstochastic DFN method, and the following results are obtained: • Determination of azimuth and intensity of the fractures in the reservoir using the geomechanical model improves the DFN model. Supervised data of faults and fractures during the geomechanical restoration increases fracture generation and identification in all scales that correspond to the reality of the reservoirs. • Fractures obtained by the geomechanical restoration have higher accuracy and resolution than using only seismic data at distances far away from the wells. • Compared to the deterministic-stochastic DFN method, the 3D geomechanical restoration model is able to predict fracture intensity reliably around the well.