Naturally fractured hydrocarbon reservoir simulation by elastic fracture modeling

Petroleum Science, May 2017

Accurate fluid flow simulation in geologically complex reservoirs is of particular importance in construction of reservoir simulators. General approaches in naturally fractured reservoir simulation involve use of unstructured grids or a structured grid coupled with locally unstructured grids and discrete fracture models. These methods suffer from drawbacks such as lack of flexibility and of ease of updating. In this study, I combined fracture modeling by elastic gridding which improves flexibility, especially in complex reservoirs. The proposed model revises conventional modeling fractures by hard rigid planes that do not change through production. This is a dubious assumption, especially in reservoirs with a high production rate in the beginning. The proposed elastic fracture modeling considers changes in fracture properties, shape and aperture through the simulation. This strategy is only reliable for naturally fractured reservoirs with high fracture permeability and less permeable matrix and parallel fractures with less cross-connections. Comparison of elastic fracture modeling results with conventional modeling showed that these assumptions will cause production pressure to enlarge fracture apertures and change fracture shapes, which consequently results in lower production compared with what was previously assumed. It is concluded that an elastic gridded model could better simulate reservoir performance.

A PDF file should load here. If you do not see its contents the file may be temporarily unavailable at the journal website or you do not have a PDF plug-in installed and enabled in your browser.

Alternatively, you can download the file locally and open with any standalone PDF reader:

Naturally fractured hydrocarbon reservoir simulation by elastic fracture modeling

Naturally fractured hydrocarbon reservoir simulation by elastic fracture modeling Mehrdad Soleimani 0 0 Faculty of Mining, Petroleum and Geophysics, Shahrood University of Technology , Shahrood , Iran Accurate fluid flow simulation in geologically complex reservoirs is of particular importance in construction of reservoir simulators. General approaches in naturally fractured reservoir simulation involve use of unstructured grids or a structured grid coupled with locally unstructured grids and discrete fracture models. These methods suffer from drawbacks such as lack of flexibility and of ease of updating. In this study, I combined fracture modeling by elastic gridding which improves flexibility, especially in complex reservoirs. The proposed model revises conventional modeling fractures by hard rigid planes that do not change through production. This is a dubious assumption, especially in reservoirs with a high production rate in the beginning. The proposed elastic fracture modeling considers changes in fracture properties, shape and aperture through the simulation. This strategy is only reliable for naturally fractured reservoirs with high fracture permeability and less permeable matrix and parallel fractures with less cross-connections. Comparison of elastic fracture modeling results with conventional modeling showed that these assumptions will cause production pressure to enlarge fracture apertures and change fracture shapes, which consequently results in lower production compared with what was previously assumed. It is concluded that an elastic gridded model could better simulate reservoir performance. Reservoir performance model; Naturally fractured reservoir; Elastic gridding - Mehrdad Soleimani 1 Introduction To handle the complexity of reservoir heterogeneity which comes from natural fractures, the nature of the reservoir should be determined in advance (Agar and Hampson 2014). The complexity in reservoir modeling either comes from reservoir characterization (e.g., heterogeneity and anisotropy in permeability) or from the process of oil recovery (e.g., capillarity, gravity and phase behavior), or from both (Kresse et al. 2013). To resolve this complexity, a dual-porosity model (Nie et al. 2012) and recently a triple-porosity model (Huang et al. 2015; Sang et al. 2016) have been introduced to simulate fractured reservoirs. Despite many useful features of the dualporosity model, it cannot provide reliable results in reservoirs in which fractures do not intersect (Karimi-Fard and Firoozabadi 2003). Applicability of the classic double-porosity model with a constant shape factor to low-permeability reservoir simulation is questionable (Cai et al. 2015). It also does not describe discrete fractures, which are the main source of challenge in naturally fractured reservoirs (Chen et al. 2008; Presho et al. 2011; Soleimani 2016a). History matching of naturally fractured reservoirs becomes more challenging particularly when these models are represented using discrete fracture network (DFN) models in carbonate rocks (Bahrainian et al. 2015). DFN models were also represented using generally unstructured grids to achieve high degrees of geological realism. Karimi-Fard et al. (2004) introduced an efficient discrete fracture model applicable for general-purpose reservoir simulators. Accurate representation of each individual fracture requires use of unstructured grids which can consider effects of the fracture aperture (Mi et al. 2014). Jiang and Younis (2015) implemented a lower-dimensional DFN model based on unstructured gridding for handling complex fracture geometry in simulated formation. Although DFN models have several advantages, they cannot be used directly with standard history matching techniques for strongly heterogeneous reservoirs with complex fracture geometry (Yan et al. 2016). Figure 1a shows a real structure and shape of different porosities in a limestone cube. Figure 1b, c shows conventional dual-porosity and DFN approaches that model the realistic fractures in a limestone fractured reservoir. In such models, if in the course of history matching, more fractures are added, moved or deleted, the model must be regridded from the beginning (Salimi and Bruining 2010). In this study, we try to increase realism of the fracture model by considering shape of the fracture and the variation of its properties through production and pressure regime change in fluid flow procedure. Fracture characteristics which change through time of production are included in elastic properties of fractures. The elastic gridding scheme was introduced in this study into the fracture modeling scheme to handle heterogeneity of a complex reservoir. This approach was then applied on a fractured limestone reservoir in southwest Iran. Results of the application of the proposed strategy for fracture modeling in the study field show production of more water in comparison with conventional fracture modeling result. 2 Fracture simulation methods for carbonate rocks The intersection relationships of fractures are probably very complex in a realistic DFN model (Zhang 2015). Bisdom et al. (2016) investigated impact of in situ stress and outcrop fracture geometry on hydraulic apertures in reservoirs. They stated that each fracture network containing fractures is created of at least one fracture set, but is not necessarily limited to it. By proposing a fracture propagation model using multiple planar fractures with a mixed model, Jang et al. (2016) stated that there is a large discrepancy in reservoir volume stimulation, because of a number of intersections of fracture connectivity. Heffer and King (2006) introduced a spatial correlation function of fractures as displacement strain vectors using renormalization techniques in representation of stochastic tensor fields for strain modeling. Masihi and King (2007) applied this method to generate fracture networks based on the assumption that the elastic energy in the fractured media follows a Boltzmann distribution. Koike et al. (2012) used geostatistical fracture distribution and fracture orientation (strike and dip) in simulation of the fracture system to estimate the hydraulic conductivity. Bisdom et al. (2016) also proved that the fracture orientation and the associated hydraulic aperture distribution have stronger impact on equivalent permeability than length or spacing. Thus, spatial correlation of fractures is the most important parameter in any fracture networks model or gridding scheme. To make this correlation between fractures, various methods are introduced for assigning precise values from fracture characteristics to the model of fractures. Among them, various approaches of using outcrop fracture characteristics such as fracture spatial distribution, length, height, orientation, spacing and aperture are widely used for model regularization (Wilson et al. 2011; Hooker et al. 2012). Lapponi et al. (2011) used outcrop data to construct a 3D model in a dolomitized carbonate reservoir rock from the Zagros Mountains, southwest Iran. Lee et al. (2011) studied the spatial fracture intensity effect on hydraulic flow in fractured rock. They used outcrop for simulation of spatial fracture intensity Conventional dual porosity model 3D discrete fracture model Fig. 1 Dual-porosity conventional model a actual porosity model, b sugar cube representation and c discrete fracture model for a fractured reservoir (Biryukov and Kuchuk 2012) distribution. They have built three spatial fracture intensity distribution models and showed that flow vectors are strongly affected by spatial fracture intensity. They also proved that the higher the fracture intensity, the higher flow velocity. Boro et al. (2014) presented a workflow to construct an upscaled fracture model based on outcrop studies in a carbonate platform. It is important to note that fractures in the outcrop might have been affected by surface processes like weathering and stress release. The overall analysis, however, can help to constrain possible scenarios on fracture populations that may be relevant to the subsurface reservoir. Malinouskaya et al. (2014) illustrated a method to rapidly estimate permeability of a fracture network, using fracture data from outcrops of a Jurassic carbonate ramp. The method proposed by Maffucci et al. (2015) used outcrop data combined with a discrete fracture network (DFN) model to increase the reliability of fracture system characterization in the case of limited data for carbonate rocks. Connectivity of fracture networks in carbonate rock is dependent on orientation, size distribution and densities of the different fracture sets. These parameters also define size of the blocks enveloped by fractures (i.e., the matrix block size), which is generally used to model transfer of fluids between the matrix and fractures (Wennberg et al. 2016). To incorporate interaction between the matrix and fractured media, Huang et al. (2014) divided the fractured porous medium into two non-overlapping subdomains. One domain has a continuum model in the rock matrix, and the other, in deep fractured and fissure zones, is described by a DFN model. Then, they coupled these domains to simulate groundwater flow in their case study. However, Boro et al. (2014) stated that in general, fracture intensities, apertures and their intrinsic permeability would have more significant impact on the permeability of the field. Fracture shape and orientations are more important in affecting connectivity. Bisdom et al. (2016) also showed the strong importance of fracture orientation and associated hydraulic aperture distribution on equivalent permeability. In general, it is widely believed that fluid flow is affected by heterogeneities at all scales, from millimeter scale (porosity) to kilometer scale (Shekhar et al. 2014). In the present work, field evidence of different solutions and fractures in limestone is used to certify the nature of elasticity of fractures through time of production. The proposed strategy also exposes incorrect assumptions in conventional fracture modeling, which has a great impact on history matching studies. Afterward, the concept of allocating each fracture to a fracture set and subsequently to a fracture network was considered by comparing the formation outcrop and considering the elasticity nature of the fractures, which comes from formation fluid pressure and/or regional stress in modeling. 3 Elastic fracture modeling Dennis et al. (2010) stated that the physical structure properties and complexity of fracture characterization both have a significant effect on fluid flow in fractured rock. Thus, they proposed that the fracture zone should be characterized fully before simulation. Wang et al. (2016) studied the flow stress damage and reservoir responses to injection rate under different DFN-connected configuration states. Their results proved significant influence of the hydraulic pressure flow on the properties of hydraulic fractures. Generally, Gan and Elsworth (2016) stated that in any fracture modeling for complex reservoirs, the upcoming assumptions should not be neglected: Fractures initiate from flaws and the process is controlled by the elastic stress around them, the material surrounding the flaws can be viewed as continuum media, and individual flaws are spaced widely enough so that stress anomalies associated with each do not overlap. These assumptions are necessary for analytical formulation and therefore will be retained in the analysis, although a slight degree of plasticity is allowed (Wang and Shahvali 2016). However, in any elastic gridding, it should be kept in mind that a fracture is initiated when the maximum stress concentration occurring on the critical flaw boundary reaches the strength of the material which surrounds the flaw. Fracture extension also occurs from the tensile and not the compressive stress concentrations under both tensile and compressive loading. Fan et al. (2012) stated that microcracks induced by the excess oil/gas pressure may propagate and form an interconnected fracture network. This indicates also that during production, different pressure regimes could change characteristics of fractures. Not only might they be closed in the case of pressure drop or reservoir depletion, but they also might be widened due to high production rates and excessive pressure of fluid flow to the walls of fractures or cracks. It also might create fractures which make connections between vugs, while cracks could be widened and/or become fractures. Fan et al. (2012) investigated mechanism of fracture propagation by a linear elastic model. They have shown that critical crack propagation takes place if the intensity of the induced stress reaches the fracture toughness of the reservoir rock. On the other hand, subcritical crack propagation occurs in the rock when the stress intensity has not reached the fracture toughness of the reservoir rock, but exceeds a threshold value, which is usually a fraction (e.g., 20%– 50%) of fracture toughness of the reservoir rock. This conclusion states how important it is to consider fracture– matrix interaction and/or boundary condition of fractures in accurate flow simulation. Subsequently, Hassanzadeh and Pooladi-Darvish (2006) considered the time variability of the fracture boundary condition by the Laplace domain analytical solutions of the diffusivity equation for different geometries of fractures in constant fracture pressure through a large number of pressure steps. Guerriero et al. (2013) proposed an analytical model that could pave the way to a full numerical model allowing one to calculate the pore pressure within fractures, at several scales of observation, in a reasonable time. They also suggested that the model also allows one to obtain a better understanding of the hydraulic behavior of fractured porous rock. However, not only the pressure regime change, but other factors affecting fracture shape and apertures could be accounted for by considering elastic behavior of fractures in the model construction during gridding and throughout production history matching investigation. Bisdom et al. (2016) stated relationships between the fracture geometrical parameters and some other parameters such as the stress applied in the medium. However, for this specific case, accurate relationship between degree of elasticity and fracture’s aperture would be defined only by core analyses in a wide range of applied stresses from pore fluid to the walls of fractures. However, previous studies have shown that this relationship is linear in a narrow range of applied stress (Bisdom et al. 2016). Elastic gridding is a variant of the grid optimizationbased technique on a length functional with a non-Euclidean metric tensor. In creation of the grid, it is considered to be a system of springs connecting neighboring grid vertices along the grid lines. On the other hand, the problem of grid optimization is thereby reduced to a problem of elasticity and the problem of translating grid optimization criteria into criteria for assigning spring constants to grid lines. After having assigned values to the spring constants of the elastic grid, the grid vertex positions can be found by solving the equilibrium equation for the elastic system. In case of reservoirs with parallel flow channels, strong pressure regime change will cause fracture apertures to become wider to some extent. However, not all the fractures in carbonate rocks, regardless of their sizes, are responsible for fluid flow. Observations of fractures in core and outcrop indicate that flow in open fractures in carbonate rock tends to be channeled rather than through fissures. Most of the flow takes place along a few dominating channels in the fracture plane, whereas most of the fracture plane is not effective for fluid flow. Wennberg et al. (2016) stated that the effect of channeled flow should be taken into account during evaluation of fractured carbonate reservoirs and building dynamic flow models. However, the consequence of the aperture variation is that the fluid flow in fractured carbonate reservoirs will tend to be channeled instead of that of fissure-type flow, which is the assumption in most flow simulators (Wennberg et al. 2016). Figure 2 shows different solutions and fractures of the Sarcheshmeh limestone formation outcrops in the Shahrood area, Iran. Major fractures shown by black lines (which are considered as the pathway of fluid) are oriented normal to the maximum stress direction, while red lines show minor fractures crosscut the major ones. The nature of major fractures is that these are the main fluid flow channels, while minor fractures are not necessarily important for that purpose. This is a typical fracture pattern in limestone reservoirs in Iran. This pattern would make only those parallel fractures, which are the only pathway of fluid flow, undergo effects of pressure regime change through high rate production. This effect will cause aperture widening, which changes the shape of the fracture, the parameter that is going to be considered in fracture modeling in this study. Another high production rate effect in such reservoirs is connecting vugs by propagation and/or widening of fractures. Lapponi et al. (2011) showed that vuggy porosity seems to increase porosity only locally and to a limited extent, developing a non-connected pore network. However, fracture propagation under pressure or high production rate would make connections between vugs. This is not an important effect in production, since these connecting fractures show small apertures, not appropriate for fluid flow. Figure 3 shows an outcrop of the same formation in the same location and an example of vugs connected by fracture propagation under pressure. However, in this study reservoir, the thermal fracturing is not planned in the master development plan of the field and the production history of the reservoir also shows some degree of overpressure fluid in the first periods of initial production, when accurate data were not available. Thus, neither the thermal fluid injection nor the fluid overpressure was considered here. The majority of fractures that we had in the study reservoir were of the type of fractures shown in Fig. 2. As was previously mentioned, some fractures connect vugs, which results in fracture propagation and/or fracture shape. All the main fractures responsible for fluid flow are modeled as flat planes in fracture modeling and are fixed through the production simulation procedure. However, in the proposed strategy, these fractures will change in shape through history matching. This is what we called the elastic gridding modeling. Figure 4 shows an example of an elastic grid and the elastic nature of fractures in modeling. As it was previously mentioned, fracture properties might change through high pressure regime change. The more the pressure regime changes in a short time, the more this will cause more change in the shape of fractures. Figure 4a shows a conventional model of fractures. Figure 4b, c illustrates shape change in the same fracture after high pressure regime change, modeled by elastic fracture modeling. Figure 4d also shows change in fracture aperture modeled by elastic fracture modeling. Fig. 2 a Outcrop of parallel fractures in the limestone Sarcheshmeh Formation, Iran. Only fractures shown by black lines are responsible for fluid flow. Fractures defined by red lines do not make flow paths. b The same formation in another part of the study area Fig. 3 Outcrop of the Sarcheshmeh Formation, same location as in Fig. 2, Iran. a Various porosity types in the outcrop, b connected and nonconnected vugs and worm channels (red lines) However, the most important parameter that changes the permeability of the reservoir rock is the dilation of fractures. Unlike other fracture parameters, the dilation degree of a fracture under stress cannot be described from core sample tests. Clearly, long fractures cause the core to fall apart during the experiment. Core recovery is also very poor in intensely fractured intervals, and the stress release when the core is taken to the surface will affect the observed apertures (Wennberg et al. 2016). Use of electrical image logs also suffers from uncertainty if the absolute aperture is large. This problem would be boosted if we want to accurately measure the dilation of the fracture. Therefore, it necessitates that the relationship between aperture dilation and applied stress is defined by numerical analysis and other permeability tests in various conditions. Taron et al. (2014) tested fracture dilation in a geothermal system. However, Min et al. (2004) and Farahmand et al. (2015) derived relationships between aperture dilation and applied stress. Min et al. (2004) studied completely all cases and we used the results that they derived in their complete study. Min et al. (2004) have stated that the exact extent of shear dilations of fractures can only be identified through numerical experiments. Their experiments showed that on the one hand, equivalent permeability decreases with increase in stresses, when the differential stress is not large enough to cause shear dilation of fractures. On the other hand, the equivalent permeability increases with the increase in differential stresses, when the stress ratio was large enough to cause continued shear dilation of fractures. Fig. 4 a Conventional model of matrix and fractures (Karimi-Fard 2013). b an example of change in shape of a fracture modeled by elastic fracture modeling, c another change in shape and d change in fracture aperture modeled by elastic fracture modeling In this case, shear dilation is the dominating mechanism in characterizing the stress-dependent permeability. They have proved that the maximum contribution of dilation is more than one order of magnitude in permeability. Thus, we have used this role in our study for fracture dilation. Since the elastic system is attached to the boundary of the model, it would not collapse. For both convex and concave domains, the elastic method is also very flexible in fracture geometry and robust to system collapse (Fig. 5). Figure 5 shows an example of a conventional grid, and Fig. 5b illustrates only a schematic of the same system but with elastic grids. Since we do not exactly know how the grids would differ during simulation, the applied pressure and elastic properties of rock would define that thus Fig. 5b is an example of any shape of grids with the unique geometry. Figure 6a depicts a simple case of conventional fracture modeling in a medium. Figure 6b shows example of a low allowed degree of elasticity and/or low pressure regime change, and Fig. 6c exhibits a sample of high degree of elasticity allowed in the model due to possible high production rates with high fluctuation in the pressure regime. Every parameter that is going to model the elastic behavior of a fracture should be defined by numerous experimental tests on core samples. The pore fluid effects needed to calculate elastic properties of fluid-saturated rock could be obtained also in each study from the simulation case. However, in this study, we did not have enough data to derive an implicit model for elastic behavior simulation of a fracture. However, it is not only the matter of data, but it is about the matter of accurate and implicit relationship between the applied stress and strain of the medium and elastic properties and behavior of fractures. Thus, in this study, we used an explicit relationship between the stress applied to a fracture and the allowed degree of elasticity of the fracture. Consequently, we should define a maximum degree of change in shape that we consider for a fracture and the value of changes in curvature of the fracture. The former is defined based on the Bulk and/or Young’s modulus of the rock. For hard rocks, the lowest grade of change in shape is allowed and it is vice versa for soft rocks. The latter needs more explanation. At first, it should define whether curvature has any effect on permeability change or it is only a fracture shape change, without effect on permeability. According to Fig. 4, it changes the permeability only if both sides of the fracture experience convex curvature, which increases the permeability. In case of same convexity or concavity of fracture walls, no changes in permeability would happen. S p rn ig S p rn ig Fig. 5 a Example of a conventional grid and b one out of thousands of examples of an elastic grid In this study, we assumed the first case for our fractures, as shown in Fig. 4d. The degree of curvature is defined based on the toughness of rock. Although the exact degree of convexity and curvature of the fracture wall should be defined by core sample test, it could be defined as a linear function of applied stress for medium to rocks. 4 The study reservoir The study reservoir is an extensively faulted anticline 74 km long and 6–8 km wide, located in the Dezful Embayment, southwest Iran (Fig. 7a). The study anticline is an asymmetrical anticline with a NW–SE trend and a sinuous axis. The target formation is a prominent carbonate unit of Oligocene and early Miocene age called the Asmari Formation. This formation in the study field contains limestone, dolomite and minor marl and shale (Abdollahie Fard et al. 2006). The study field was divided into five sectors based on geological and engineering data. Most of the boundaries of these sectors are in agreement with faults. These sectors are depicted in the underground contour (UGC) map of the target formation shown in Fig. 7b. It should be mentioned that the UGC map in Fig. 7b was obtained from time–depth conversation of 3D migrated seismic data with well top adjustment supervision. However, an advanced migration algorithm is needed for imaging in such complex media (Soleimani 2016b). Limited lateral communication, the presence of faults acting as barriers, hydrodynamic system and imbalanced offtake have caused the Asmari Formation to be divided into several compartments (Sherkati and Letouzey 2004). The oil production from the study reservoir has been very imbalanced mainly because early development took place in the central and northwestern part of the field, where oil column was thicker and the wells had higher productivity. In spite of extensive fracturing of the Asmari Formation, a pressure differential has developed in the gas zone across the field due to imbalanced oil and gas production. 4.1 Petrophysical data Limited cores were cut from the Asmari Formation for routine and special analyses. Total length of cores cut was 550 m out of which 421 m was recovered. The mean porosity of plugs cut from 421 m of cores in six wells was 9.8%, while 21% of samples have porosities less than 4% and only 3% of samples have porosity more than 20%, implying that the Asmari Formation is a low porosity reservoir (Hoseinzadeh et al. 2015). The median permeability of cores was calculated to be 0.43 mD with 60% of samples having permeability of less than one milli-darcy, indicating a very low permeable carbonate rock matrix. Production logging tools (PLT) logs were recorded in 15 oil wells and in 5 gas wells. Reviewing of the PLT log results indicates that distance between flowing intervals varies from 1 m to the maximum of 44 m. This indicates an active mechanism of fracture production in this field. 5 Zonation and fracture study Porosity distribution in the Asmari Formation is very diverse. Therefore, division of this formation into several zones and subzones is unavoidable. Based on petrophysical and petrographic characteristics, the target formation was divided into seven different oil-bearing zones with the water column as zone 8. Consequently, zones 1, 2, 6 and 7 were divided into two subzones (Table 3 of Appendix). From the study of cores, some parameters were extracted Fig. 6 a Conventional modeling in which fractures were modeled by flat planes with no change in shape through production (Karimi-Fard 2013), b Low degree of elasticity allowed, and c high degree of elasticity allowed in elastic fracture modeling in the proposed strategy such as type of fracture; spacing (number of fractures per meter); dip of fractures (relative to the core axis in the core description tables and relative to horizontal in the interpretation); width of fractures measured in laboratory conditions; and filling mineral type. Among these parameters, type of fracture and filling mineral type are classified and coded as shown in Tables 1 and 2, respectively. Occurrence of anomalous losses of mud during drilling is often the first indirect indication of the existence of a naturally fractured interval. Comparison of the mud loss rate in different wells may give the intensity of fracturing in different parts of the reservoir (Xia et al. 2015). In spite of the fact that there are many parameters that can affect the mud loss rate and inherent errors in the results, such as mud weight, weight on bit, stroke per minute of pumps and rock type, the mud loss rate is one of the most important parameters in indirect observations of fractures. Other information that could be useful in fracture studies includes productivity index (PI), PLT data and maximum daily flow rate. Production rate is usually dependent on the extent of fracturing of the formation in this low-matrix permeability reservoir (Xia et al. 2015a). Therefore, higher production rates correspond to a more extensive fracturing of the strata in the well location (Fig. 8). To determine the fractured zones within the reservoir, a radius of curvature method (Yoshida et al. 2016) was used. The 3D curvature model was not shown here, but it was obtained from 3D migrated seismic data and application of maximum and minimum curvature attributes. After applying this method to the target formation, we have concluded that the most common types of fracture are type 1 (cross-axial tensional and conjugate shear fractures) and type 2 (axial tensional and conjugate shear fractures). Interpretation of mud losses, daily flow rate, PI, PLT and other data suggest that sector 1N and partially sector 4 are highly fractured. Fracturing in sector 3 is relatively intensive and not extensive in sector 1S. Figure 9a shows the fracture quality map, and Fig. 9b illustrates strain distribution map on top of the target formation. The strain regime here is of shortening type. The mud loss, PI and PLT integration was performed by finding the relationship between fracture quality index and these three parameters in GIS media. 6 Elastic fracture model generation All the required data for elastic fracture and geological model construction were collected. Petrophysical data contain two series of net-to-gross (NTG) ratio, porosity and water saturation related to each horizon. Monte Carlo statistical simulation was used to obtain a 3D distribution of fracture spacing density and, consequently, matrix block size. Input probabilistic frequencies were scaled in such a way as to have the mode of matrix block height distributions changing from approximately 3 m at the top to approximately 6 m at the bottom of the Asmari Formation, respectively. Matrix–fracture communication was defined by considering elasticity behavior of fractures. At that stage, matrix and fracture properties were considered as first approximations to be modified during history matching. After construction of the model, a geological model containing 633 9 451 9 12 (3,425,796 cells) mesh cells was prepared including all faults. Configuration of the elastic gridding model was based on fault traces and a contour line of -2600 m from the structure map on top of the Asmari Formation. The elastic grid has 133 blocks in Fig. 7 a Location of the study field is shown by the rectangle (Soleimani and Jodeiri-Shokri 2015). b Zonation of the study field and underground contour map (UGC) of top of the Asmari Formation Table 1 Classifying different types of fracture in the cores Type of fracture Closed or hairline Table 2 Classifying different types of filling minerals in the cores one direction and 17 blocks in the other direction with 30 layers. The elastic grid structural elevations were obtained from underground contour maps and verified against well picks. Gross thickness of zones (subzones) was based on isopach and isochore maps. Porosity versus normalized depth profiles were used as guides to split zones (subzones) into 30 layers. Figure 10 displays a map view of the elastic grids on the top of the target formation. Zone (subzone) matrix porosity and NTG ratio were obtained from geologic maps and then downscaled to layers based on the porosity–normalized depth profiles. 1 7 0 0 9 0 0 Fig. 8 Production index based on production of the wells in each zone Fig. 9 a Fracture quality map based on mud loss, PI and PLT and b strain distribution map of the study field. Small symbols show well location on both maps Fig. 10 Map view of the grids on the top of the target formation Matrix permeability was based on matrix porosity and permeability–porosity correlations. A fracture porosity distribution was generated based on the fracture porosity map. This map was horizontally scaled according to the qualitative map of fracture intensity. It was vertically scaled also in accordance with empirical correlations reflecting reduction in fracture intensity from top to bottom of the Asmari Formation (Fig. 11). Well trajectories and perforation, production, and static pressure histories were incorporated in the model. Sporadic information about salt and/or water production was only available for a limited number of oil wells. Individual well gas/oil ratios were found to have resulted from prorating rather than from measurements; on this basis, they were excluded from history matching. 6.1 Model initialization Providing the model with initial distributions of water saturation, reservoir pressure and hydrocarbon components were the subject of this part of the study. To obtain distribution models of these characteristics, we took seven available drainage capillary pressure curves (not shown here) and scaled them in accordance with the following equation: where Rw is the formation water resistivity, Rt is the formation resistivity, and w = f(Rw, R, ø). The right side of Eq. (1) was derived from well logs and averaged for different regions and zones (subzones). To have a reasonable value for the right side of equation, a primary blocking step of the target formation was performed for averaging in each block. Subsequently, a weighted averaging based on the thickness of each zone/subzone was performed. Due to the large number of obtained maps with different blocking approaches and based on different interpretation of the averaged map result with the help of geological and hydrological data, the maps are not shown here. Then, Swi was estimated from Eq. (1) for each simulation model elastic grid block based on its matrix porosity. The resulting Swi values were introduced in the model as connate water saturation. To obtain different oil–water contacts (OWC), two aquifers and six equilibration regions were introduced into the model (Fig. 12). Two analytical aquifers were introduced instead of one aquifer, while the northwest part of the reservoir was completely disconnected from any analytical aquifer. Aquifer 2 is stronger 2 7 3 8m 2 1 4 1m No analytical aquifer connections Fig. 11 a Fracture porosity and b fracture intensity map of top of the target formation Fig. 12 a Two different aquifers in the field and b six different OWC zones, all introduced into the model than the other aquifers. The aquifers are separated by faults, while six OWC zones are controlled by the pressure regime of the field. The water saturation in the matrix also could affect the OWC in different zones. Thus for model initiation, the pressure regime and the matrix water saturation (besides the fracture water saturation) should be introduced as the initial condition for model running. Figure 13 shows the pressure regime and matrix water saturation in the model. 6.2 Preliminary simulation The study field’s production history is complicated, and the aquifer strength varies from the southeast to the northwest of the field. Thus, introducing the elastic behavior of fractures into the model was proposed here to realize the history matching of the reservoir production. However, running time was a significant issue right from the beginning of the matching process. Only using advanced hardware (a total of 12 CPUs of 3 GHz each) and implementing extensive debugging allowed the reduction of the running time for a history match to nearly 75 h. In the first step, the conventional history matching objectives in the traditional gridded model showed the entire field’s solution of gas production and water production only in the oil wells located southwest of the field. To introduce the elastic behavior of fractures into the model, initial global modifications involved elevating the aquifer permeability, elevating fracture permeability in x direction by a factor of 10 and reducing fracture permeability in z direction by a factor of 10 to compensate for the elasticity of microfracturing observed in the cores. Examination of elastic model runs showed that oil flow from matrix to fractures should be increased to match water and gas production from the oil wells. Elevating of this Matrix pressure, psia 7000 5ia000 sp , R 3P000 F 1000 0 Fig. 13 a Matrix water saturation and b the pressure regime both introduced into the model for initiation Fig. 14 Fracture water saturation in a conventional fractures and b elastic fracture modeling Fracture water saturation Fracture water saturation Fracture oil saturation Fracture oil saturation Fig. 15 Fracture oil saturation in a conventional fracture and b elastic fracture modeling oil transfer was done by modifying the imbibition capillary pressure. Individual oil well modifications included changing capillary pressure curves assigned to the matrix grid blocks surrounding the well, varying the height of matrix blocks, modification of the pore volume of initially oil- and gas-saturated grid blocks and modification of the fracture permeability in x, y and z directions. Figure 14 shows the conventional and elastic fracture water saturation model. As it can be seen on the elastic fracture model, wells located in the southeast (SE) of the field may produce more water. This is due to the increase in cracks and/or fracture apertures, which means previous small cracks were filled with water. Figure 15 shows elastic and conventional models of fracture oil saturation. Again, the southeast (SE) part of the field does not produce more oil after a while according to the elastic modeling results. This phenomenon is less observed for gas saturation. Figure 16 shows result of elastic and conventional fracture modeling. As a result, change in shape of fractures, (increasing aperture) would replace oil by gas in the vicinity of oil-gas contact. This is obvious in Fig. 16. Fracture gas saturation Fracture gas saturation Fig. 16 Fracture gas saturation in a conventional fracture and b elastic fracture modeling 7 Conclusions There are many problems in fractured reservoirs that may require new approaches in fracture modeling based on advanced concepts. Some problems cannot be simplified beyond a certain limit. The complexity, scale and uncertainty of natural systems are the primary reasons for most of the modeling difficulties which have to be handled. The objective of this work was to develop a prototype workflow for history matching in naturally fractured reservoirs. This concern is also with the problem of generating more realistic computational grids for reservoir simulations. Here we face the particular problems connected to the complexity of the reservoir. The proposed strategy combines the accuracy of fracture modeling with the efficiency of elastic gridding. Elastic gridding is very simple, and in principle, the method should be well suited for weighing between mutually conflicting optimization criteria with highly nonlinear and discontinuous cost functions. The strategy was applied on a complex reservoir from southwest Iran. The strategy was observed to be reasonably effective in achieving field-level agreement in total oil and water production rates. Aperture increase in Table 3 Zonation of the target formation in the study reservoir Zone It is 28 m thick. Mainly consists of 1–2 dolomite, limestone, sandy dolomite, thin layer of shale and sandstone cracks will mean they will be filled by water, while it replaces gas by oil in the vicinity of the gas oil contact. Thus, elastic fracture modeling shows that more water will be produced compared to what was assumed by conventional fracture models, which do not take into account fracture properties change through production. Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creative, which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Table 3 shows zonation of the target formation. As was mentioned in the text, the target formation is divided into seven zones with the aquifer as zone 8. This zonation is based on the geological and petrophysical information. In each zone and subzone, lithology porosity, permeability and NTG of the zone are fully described. This zonation was used as the base of the reservoir modeling and simulation. Porosity and permeability Most of pore volume is filled with cements. The average porosity is 10%. Water saturation increases from 10% in the southern flank to 30% in northern flank as well as from west to east Average porosity is 8% in sector 4 and varies from 10% to 14% in both flanks. In eastern portion, it ranges from 11% in the southern limb to 7% in the northern flank. Water saturation varies from 20% to 25% in the southern flank to 55%–60% in the northern flank The NTG is about 90%–95% in the western half and is variable in the eastern half; it varies from 80% to 100% in the southern limb, from 60% to 42% in the northern flank The calculated NTG ratio of the subzone is about 90% in the central portion of the field, 60% in the eastern portion of the northern flank and 45% in the sector 4 Table 3 continued Zone This subzone predominantly contains 2–1 dolomite and limestone. On average, it is 41 m thick Zone Its thickness is 27 m. Composed of 2–2 limestone, dolomitic limestone, calcareous dolomite and partly thin shaly layers. Thickness of 40 m. Contains alternation of limestone, dolomitic limestone and dolomite 18 m thick and is composed of dense limestone and partly of dolomite. A thin shaly bed is observed. Thickness is 60 m. Composed of limestone, partly of dolomite and of thin shaly layers Zone 18.5 m thick and composed of 6–1 dolomite and intercalation of argillaceous layers Zone Thickness of 33 m and composed of 6–2 limestone and partly of dolomite. Dies out at the western part Zone 37 meters thick and consists of 7–1 limestone and partly of dolomite and of a few thin shaly layers Zone Thickness 38 m. Composed of 7–2 limestone, marly limestone and shale Porosity and permeability NTG In most of the wells throughout the field, Average NTG ranges from 85 to 95% in porosity varies from 8% to 10%. Water the southern flank with exception of a saturation in this subzone is distributed few wells having NTG about 65%. In similarly to the previous subzones with the northern flank of the central and the minimum saturation of 20% at the eastern portions, NTG averages 65% southern flank and the maximum and 35% accordingly saturation of 60% at the northern flank. Porosity decreases from 12% in the southern flank to 7%–10% in the northern limb. In sector 4, porosity decreases to 5%. Water saturation increases from 15% in the south to 55% in the western part and to 80% in the eastern part of the northern flank Porosity varies form 15% in the southern limb to 9%–11% in the northern flank Water saturation increases from 10–20% in the southern flank to 50% in the western part of the northern flank The average porosity of the zone ranges from 8% to 14%. In the central and eastern portions of the field water, saturation increases uniformly from 20% in the southern limb to about 70%–80% in the northern flank Average water saturation ranges between 25% and 35%. A general trend shows a south to north increase in the water saturation The average porosity in the eastern part is 7%. It is 16% in the southern and 8% in the northern flank The water saturation is 40% in the south flank and 80% in the northern limb The average calculated porosity of the subzone reduces from 16% in the southern limb to 9%–10% in the northern flank. Water saturation changes from 25% in the southern flank to 50% in the northern flank The average porosity is 12%. In the southern flank, the porosity is 13% decreasing to 10% in the northern limb. As a result, the water saturation value of 75% was reached In sector 1, porosity is 12% in the southern flank and 7% in the northern flank. Average porosity values of 8% to 10% were assumed In sector 4, water saturation varies from 45% in the southern flank to 90% in the northern limb NTG decreases from 90% to 95% in the southern flank to 25%–50% in the northern flank. In sector 4, NTG value averages 80% Since rock quality is relatively good, NTG ratio is almost 100% in the southern flank and decreases to 85% in the northern flank Except sector 4, NTG ratio decreases from 75% to 95% in the southern flank to 25% and to 35% in the western and eastern portions of the northern flank, respectively The NTG ratio varies between 80% and 95% throughout the field NTG ratio ranges between 13 and 54%. It is 95% in the south flank and 60% in the north flank of the central part of the field NTG ratio ranges from 90 to 95% in sector 1, between 50 and 60% in sectors 2 and 3, and from 66% to 100% in the southern flank of sector 4 NTG ratio varies from 80% to 100% decreasing to 34%. In the eastern part of the field, the average NTG ratio ranges between 75 and 80% Table 3 continued Thickness 72 m and composed of shale, marl and marly limestone Porosity and permeability Porosity of 9%–10% was estimated. Water saturation in the southern and northern flanks was evaluated to be 40% and 55%, respectively The NTG ratio in sector 4 (37%) and in sector 1 (53%), decreasing from about 60% in the southern flank to about 35% in the northern flank Abdollahie Fard I , Braathen A , Mokhtari M , Alavi A. Interaction of the Zagros fold-thrust belt and the Arabian-type, deep-seated folds in the Abadan plain and the Dezful embayment SW Iran . Pet Geosci . 2006 . doi:10.1144/ 1354 - 079305 - 706 . Agar SM , Hampson GJ . Fundamental controls on flow in carbonates: an introduction . Pet Geosci . 2014 . doi:10.1144/petgeo2013- 090 . Bahrainian SS , Daneh Dezfuli A , Noghrehabadi A. Unstructured grid generation in porous domains for flow simulations with discretefracture network model . Transport Porous Media . 2015 . doi:10. 1007/s11242- 015 - 0544 -3. Biryukov D , Kuchuk FJ . Transient pressure behavior of reservoirs with discrete conductive faults and fractures . Transport Porous Media . 2012 . doi:10.1007/s11242- 012 - 0041 -x. Bisdom K , Bertott G , Nick HM . The impact of in situ stress and outcrop-based fracture geometry on hydraulic aperture and upscaled permeability in fractured reservoirs . Techtonophysics . 2016 . doi:10.1016/j.tecto. 2016 .04.006. Boro H , Rosero E , Bertotti G. Fracture-network analysis of the Latemar Platform (northern Italy): integrating outcrop studies to constrain the hydraulic properties of fractures in reservoir models . Pet Geosci . 2014 . doi:10.1144/petgeo2013- 007 . Cai L , Ding DY , Wang C , Wu YS . Accurate and efficient simulation of fracture-matrix interaction in shale gas reservoirs . Transport Porous Media . 2015 . doi:10.1007/s11242- 014 - 0437 -x. Chen Y , Cai D , Fan Z , Li K , Ni J. 3D geological modeling of dual porosity carbonate reservoirs: a case from the Kenkiyak pre-salt oilfield Kazakhstan . Pet Explor Dev . 2008 . doi:10.1016/S1876- 3804(08) 60097 - X . Dennis I , Pretorius J , Steyl G. Effect of fracture zone on DNAPL transport and dispersion: a numerical approach . Environ Earth Sci . 2010 . doi:10.1007/s12665- 010 - 0468 -8. Fan ZQ , Jin ZH , Johnson SE . Gas-driven subcritical crack propagation during the conversion of oil to gas . Pet Geosci . 2012 . doi:10. 1144/ 1354 - 079311 - 030 . Farahmand K , Baghbanan A , Shahriar K , Diederichs MS . Effect of fracture dilation angle on stress-dependent permeability tensor of fractured rock . In: 49th U.S. Rock mechanics/geomechanics symposium, San Francisco, 2015 , ARMA- 2015 - 542 . Gan Q , Elsworth D. A continuum model for coupled stress and fluid flow in discrete fracture networks . Geomech Geophy Geo Energy Geo Resour . 2016 . doi:10.1007/s40948- 015 - 0020 -0. Guerriero V , Mazzoli S , Iannace A , Vitale S , Carravetta A , Strauss C. A permeability model for naturally fractured carbonate reservoirs . Mar Pet Geol . 2013 . doi:10.1016/j.marpetgeo. 2012 .11. 002. Hassanzadeh H , Pooladi-Darvish M. Effects of fracture boundary conditions on matrix-fracture transfer shape factor . Transport Porous Media . 2006 . doi:10.1007/s11242- 005 - 1398 -x. Heffer KJ , King PR . Spatial scaling of effective modulus and correlation of deformation near the critical point of fracturing . Pure Appl Geophys . 2006 . doi:10.1007/978- 3 - 7643 - 8124 - 010 . Hoseinzadeh M , Daneshian J , Moallemi SA , Solgi A. Facies analysis and depositional environment of the Oligocene-Miocene Asmari Formation, Bandar Abbas hinterland Iran . Open J Geol . 2015 . doi:10.4236/ojg.2015.54016. Hooker JN , Gomez LA , Laubach SE , Gale JFW , Marrett R. Effects of diagenesis (cement precipitation) during fracture opening on fracture aperture-size scaling in carbonate rocks . Geol Soc Lond Spec Publ . 2012 . doi:10.1144/SP370.9. Huang T , Guo X , Chen F. Modeling transient flow behavior of a multiscale triple porosity model for shale gas reservoirs . J Nat Gas Sci Eng . 2015 . doi:10.1016/j.jngse. 01 .022. Huang Y , Zhou Z , Wang J , Dou Z. Simulation of groundwater flow in fractured rocks using a coupled model based on the method of domain decomposition . Environ Earth Sci . 2014 . doi:10.1007/ s12665- 014 - 3184 -y. Jang A , Kim J , Ertekin T , Sung W. Fracture propagation model using multiple planar fracture with mixed mode in naturally fractured reservoir . J Pet Sci Eng . 2016 . doi:10.1016/j.petrol. 02 .015. Jiang J , Younis RM . Numerical study of complex fracture geometries for unconventional gas reservoirs using a discrete fracture-matrix model . J Nat Gas Sci Eng . 2015 . doi:10.1016/j.jngse. 2015 .08. 013. Karimi-Fard M , Firoozabadi A. Numerical simulation of water injection in fractured media using discrete-fracture model and the Galerkin method . SPE Reserv Eval Eng . 2003 . doi:10.2118/ 83633 -PA. Karimi-Fard M , Durlofsky LJ , Aziz K. An efficient discrete fracture model applicable for general-purpose reservoir simulators . SPE J . 2004 . doi:10.2118/ 88812 -PA. Karimi-Fard M. Modeling tools for fractured systems: gridding , discretization, and upscaling. Stanford University. 2013 . http:// Koike K , Liu C , Sanga T. Incorporation of fracture directions into 3D geostatistical methods for a rock fracture system . Environ Earth Sci . 2012 . doi:10.1007/s12665- 011 - 1350 -z. Kresse O , Weng XW , Gu HR , Wu RT . Numerical modeling of hydraulic fractures interaction in complex naturally fractured formations . Rock Mech Rock Eng . 2013 . doi:10.1007/s00603- 012 - 0359 -2. Lapponi F , Casini G , Sharp I , Blendinger W , Ferna´ndez N , Romaire I , Hunt D. From outcrop to 3D modelling: a case study of a dolomitized carbonate reservoir , Zagros Mountains Iran. Pet Geosci . 2011 . doi:10.1144/ 1354 - 079310 - 040 . Lee CC , Lee CH , Yeh HF , Lin HI . Modeling spatial fracture intensity as a control on flow in fractured rock . Environ Earth Sci . 2011 . doi:10.1007/s12665- 010 - 0794 -x. Maffucci R , Bigi S , Corrado S , Chiodi A , Di Paolo L , Giordano G , Invernizzi C. Quality assessment of reservoirs by means of outcrop data and discrete fracture network models: the case history of Rosario de La Frontera (NW Argentina) geothermal system . Tectonophysics . 2015 . doi:10.1016/j.tecto. 2015 .02.016. Malinouskaya I , Thovert JF , Mourzenko VV , Adler PM , Shekhar R , Agar S , Rosero E , Tsenn M. Fracture analysis in the Amellago outcrop and permeability predictions . Pet Geosci . 2014 . doi:10. 1144/petgeo2012- 094 . Masihi M , King PR . A correlated fracture network: modeling and percolation properties . Water Resour Res . 2007 . doi:10.1029/ 2006WR005331. Mi L , Jiang H , Li J , Li T , Tian Y. The investigation of fracture aperture effect on shale gas transport using discrete fracture model . J Nat Gas Sci Eng . 2014 . doi:10.1016/j.jngse. 09 .029. Min KB , Rutqvist J , Tsang CF , Jing L. Stress-dependent permeability of fractured rock masses: a numerical study . Int J Rock Mech Min Sci . 2004 . doi:10.1016/j.ijrmms. 2004 .05.005. Nie RS , Meng YF , Jia YL , Zhang FX , Yang XT , Niu XN . Dual porosity and dual permeability modeling of horizontal well in naturally fractured reservoir . Transport Porous Media . 2012 . doi:10.1007/s11242- 011 - 9898 -3. Presho M , Woc S , Ginting V. Calibrated dual porosity, dual permeability modeling of fractured reservoirs . J Pet Sci Eng . 2011 . doi:10.1016/j.petrol. 2011 .04.007. Salimi H , Bruining H. Upscaling in vertically fractured oil reservoirs using homogenization . Transport Porous Media . 2010 . doi:10. 1007/s11242- 009 - 9483 -1. Sang G , Elsworth D , Miao X , Mao X , Wang J. Numerical study of a stress dependent triple porosity model for shale gas reservoirs accommodating gas diffusion in kerogen . J Nat Gas Sci Eng . 2016 . doi:10.1016/j.jngse. 2016 .04.044. Shekhar R , Sahni I , Benson G , et al. Modelling and simulation of a Jurassic carbonate ramp outcrop , Amellago, High Atlas Mountains. Morocco. Pet Geosci . 2014 . doi:10.1144/petgeo2013- 010 . Sherkati S , Letouzey J. Variation of structural style and basin evolution in the central Zagros (Izeh zone and Dezful Embayment) Iran . Mar Pet Geol . 2004 . doi:10.1016/j.marpetgeo. 01 . 007. Soleimani M , Jodeiri-Shokri B. 3D static reservoir modeling by geostatistical techniques used for reservoir characterization and data integration . Environ Earth Sci . 2015 . doi:10.1007/s12665- 015 - 4130 -3. Soleimani M. Seismic imaging by 3D partial CDS method in complex media . J Pet Sci Eng. 2016a . doi:10.1016/j.petrol. 2016 .02.019. Soleimani M. Seismic image enhancement of mud volcano bearing complex structure by the CDS method, a case study in SE of the Caspian Sea shoreline . Russ Geol Geophs. 2016b . doi:10.1016/j. rgg. 2016 .01.020. Taron J , Hickman S , Ingebritsen SE , Williams C. Using a fully coupled, open-source THM simulator to examine the role of thermal stresses in shear stimulation of enhanced geothermal systems . In: 48th US Rock mechanics/geomechanics symposium held in Minneapolis , 2014 , ARMA 14 - 7525 . Wang Y , Shahvali M. Discrete fracture modeling using Centroidal Voronoi grid for simulation of shale gas plays with coupled nonlinear physics . Fuel . 2016 . doi:10.1016/j.fuel. 2015 .09.038. Wang Y , Li X , Tang CA . Effect of injection rate on hydraulic fracturing in naturally fractured shale formations: a numerical study . Environ Earth Sci . 2016 . doi:10.1007/s12665- 016 - 5308 -z. Wennberg OP , Casini G , Jonoud S , et al. The characteristics of open fractures in carbonate reservoirs and their impact on fluid flow: a discussion . Pet Geosci . 2016 . doi:10.1144/petgeo2015- 003 . Wilson CE , Aydin A , Karimi-Fard M , et al. From outcrop to flow simulation: constructing discrete fracture models from a LIDAR survey . AAPG Bull . 2011 . doi:10.1306/03241108148. Xia Y , Jin Y , Chen M. Comprehensive methodology for detecting fracture aperture in naturally fractured formations using mud loss data . J Pet Sci Eng. 2015a . doi:10.1016/j.petrol. 10 .017. Xia Y , Jin Y , Chen M , et al. Hydrodynamic modeling of mud loss controlled by the coupling of discrete fracture and matrix . J Pet Sci Eng. 2015b . doi:10.1016/j.petrol. 2014 .07.026. Yan X , Huang Z , Yao J , Li Y , Fan D. An efficient embedded discrete fracture model based on mimetic finite difference method . J Pet Sci Eng . 2016 . doi:10.1016/j.petrol. 2016 .03.013. Yoshida N , Levine JS , Stauffer PH . Investigation of uncertainty in CO2 reservoir models: a sensitivity analysis of relative permeability parameter values . Int J Greenh Gas Control . 2016 . doi:10. 1016/j.ijggc. 2016 .03.008. Zhang QH . Finite element generation of arbitrary 3-D fracture networks for flow analysis in complicated discrete fracture networks . J Hydrol . 2015 . doi:10.1016/j.jhydrol. 2015 .08.065.

This is a preview of a remote PDF:

Mehrdad Soleimani. Naturally fractured hydrocarbon reservoir simulation by elastic fracture modeling, Petroleum Science, 2017, 286-301, DOI: 10.1007/s12182-017-0162-5