Integrated fracture characterization of Asmari reservoir in Haftkel field

The Asmari reservoir in Haftkel field is one of the most prolific naturally fractured reservoirs in the Zagros folded zone in the southwest of Iran. The primary production was commenced in 1928 and continued until 1976 with a plateau rate of 200,000 bbl/day for several years. There was an initial gas cap on the oil column. Gas injection was commenced in June 1976 and so far, 28% of the initial oil in place have been recovered. As far as we concerned, fracture network is a key factor in sustaining oil production; therefore, it needs to be characterized and results be deployed in designing new wells to sustain future production. Multidisciplinary fracture evaluation from well to reservoir scale is a great privilege to improve model’s accuracy as well as enhancing reliability of future development plan in an efficient manner. Fracture identification and modeling usually establish at well scale and translate to reservoir using analytical or numerical algorithms with the limited tie-points between wells. Evaluating fracture network from production data can significantly improve conventional workflow where limited inter-well information is available. By incorporating those evidences, the fracture modeling workflow can be optimized further where lateral and vertical connectivity is a concern. This paper begins with the fracture characterization whereby all available data are evaluated to determine fracture patterns and extension of fracture network across the field. As results, a consistent correlation is obtained between the temperature gradient and productivity of wells, also convection phenomenon is confirmed. The findings of this section help us in better understanding fracture network, hydrodynamic communication and variation of temperature. Fracture modeling is the next step where characteristics of fractures are determined according to the structural geology and stress directions. Also, the fault’s related fractures and density of fractures are determined. Meanwhile, the results of data evaluation are deployed into the fracture model to control distribution and characteristics of fracture network, thereby a better representation is obtained that can be used for evaluating production data and optimizing development plan.


Introduction
The Asmari reservoir in Haftkel field is one the best examples of the NFRs in the world where the desired reservoir engineering elements are found during the 93-year production history of this field. It is located in the Northern Dezful Embayment in the southwest of Iran (Fig. 1) in the center of the chain of anticlinal structure of Asmari which extends Mamatain in the southeast to the Naft-Safid in the northwest. It is the main producing reservoir (Tertiary) of Early Miocene age which is overlain by the Gachsaran deposits of Miocene age which act as a seal. The Asmari Formation consists of the well-fractured limestone/dolostone sequences. SI Metric Conversion Factors: bbl/6.2898 = m 3 , Ft 3 / 35.3150 = m 3 , Ft / 3.2808 = m, (°F-32)/1.8 =°C, Psi × 6894.76 = N/m 2 (Pa), mD/1.0132E + 15 = m 2 1 3 The primary production was commenced in 1928 and continued until 1976 with a plateau rate of 200,000 bbl/ day for several years. There was an initial gas cap on the oil column. Massive oil production was caused a major pressure decline also changes in the contact levels. Thus, gas injection was commenced in 1976 to maintain reservoir pressure and stabilize contact levels.
The Haftkel anticline is an asymmetrical anticline with the southwest vergence which is controlled by thrust faults. The folding of the southwest flank is somewhat steeper than the northeast flank due to the northeast direction of thrust faults which caused the folding (Ghanadian et al. 2017). Based on the inter-limb angle, the structure can be classified as an open fold (between 70° to 120°). According to the results of 2D seismic interpretation, fold style of the Haftkel structure is a fault propagation fold as shown in Fig. 2 (Sherkati et al. 2006;Ghanadian et al. 2017;Solgi et al. 2017). Therefore, this type of folding along with the limited faults are considered into the present fracture study. Unfortunately, the quality of current seismic data is not adequate to reflect the whole structural features properly, particularly in the center of structure and plunges.
Asmari in Haftkel contains low porous carbonates with a moderate fluid capacity, which has been enhanced by several sets of fracture network. Fractures have been distinguished in the geological samples and image logs (Ahmadhadi et al. 2008;Ghanadian et al. 2017). The most important factors controlling the fracture distribution are the geometry and the folding mechanism (Sherkati et al. 2006;Ghanadian et al. 2017;Solgi et al. 2017). Recognition of fracture patterns on the reservoir structure derived attentions because of the enhanced production potential whereby the maximum density of the fractures found over the plunging noses and central parts of the structure (McQuillan 1973). On the other hand, history matching together with variations of contact levels was a challenge. This problem was further crucial in the southeast of field. Thus, one of the objectives of this study was constructing a robust reservoir model including fracture network able to predict reservoir behaviors with minimal deviations. It can also be used as a reliable tool for well design to reduce failure risks.
In this paper, we present fracture characterization at first, where all available data including petrophysical logs, oil production, productivity index, saturation pressure, as well as temperature and pressure surveys are reviewed to obtain a holistic view of fracture network. The findings of this study help us in better understanding of fracture network, hydrodynamic communication and variation of temperature across the field. Then, fracture modeling is presented, where results of fracture characterization are deployed into the fracture model to get the fracture map and the corresponding characteristics (i.e., permeability and porosity). The last section concludes this paper.

Petrophysical evaluation
Advanced and conventional well logs can be used for qualitative assessment of the fracture networks. Generally, conventional logs are not accurate enough to be used for quantitative characterization of fracture network due to log resolution issue. However, some quantitative assessment can be done on the basis of pad-based logs where log resolution is sufficiently acceptable. At larger scale, non-pad tools such as gamma ray, sonic and resistivity can be affected by the presence of fracture network depending on the pore filling materials/fluids or the geometry of fracture plans (Saboorian-Jooybari et al. 2016). Many works have been devoted to establish proper links between the presence of fracture network and their impacts on the conventional well logs; however, they have not been conclusive yet. Thus, currently conventional logs are considered as the qualitative assessment tools. In the current study, we focus on the image logs alone because quality of log data is acceptable in this field.
Image log is one of the advanced log technologies which can provide important information about the fracture appearance and its characteristics (Aghli et al. 2014;Xian et al. 2018;Shafiabadi et al. 2021;Nabiei et al. 2021). There are eight sets of image logs across the Haftkel field, where the locations of the wells are highlighted with red boxes in Fig. 3.
Here, log data and the results of interpretation of two wells at the northwest (A) and the southeast (B) are presented in Figs. 4 and 5, respectively. Majority of wells are located in the southeast flank where focus of the recent development campaign was there. It should be highlighted that image data mostly covers the area that majority of production comes from. Based on the results, fracture networks have been systematically developed in reservoir sections without any meaningful relationship to specific litho-facies or reservoir sublayers. Figure 4 represents log data of Well A whereby distribution and intensity of fracture networks are shown in last tracks. Two typical sections from the highly fractured intervals are presented in the zoomed plots distinguished by the green and red boxes. Stereonet plot in "DIP Track" shows fracture setting within 50 m of reservoir section whereby dominant direction is NW-SE with an average dip angle of 70˚. Figure 5 presents the results of image log interpretation in conjunction with core data in Well B. Likewise, fracture networks are distinguished across the whole reservoir section almost with the same intensity trend. Apart from the image logs, core permeability data also demonstrate some signs of fracturing particularly in the low porosity and tight facies (Fig. 5). It seems that fracturing has impacted remarkably on the weak reservoir units rather than good ones.
Although intensity has remained constant at well scale, significant changes in dip direction are observed across the field. Dominant dip direction is significantly influenced by the active compressional and tension regime moving from NW to SE flanks as shown in Fig. 6. Fracture intensity is also reduced in the same direction due to the presence of gentle structure trends. Based on the findings, dip direction is generally controlled by the folding mechanism rather than fault system. However, some scatter patterns are also expected where the major faults are present.

Pressure surveys
There are numerous static pressure surveys that have been collected through production and observation wells. They are used to find out the reservoir pressure at datum depth and evaluating reservoir performance. Figure 7 compares the profile of pressure surveys with the field oil production and gas injection profiles. The profile of field oil production rate is presented in the middle panel (Fig. 7). The first oil in HaftKel was produced in 1928 with an initial production rate of about 5,700 stb/day. Oil production was increased by infill drilling and reached to the plateau rate of about 200,000 stb/d in 1944 and maintained for almost 7 years. Production was dropped to about 14,000 stb/d along with rising up water-oil contact (WOC) in 1976. In order to maintain reservoir pressure and WOC, gas injection campaign in gas cap was commenced in June 1976 at an average rate of 400 MMscf/Day for several years then regulated at 100 MMscf/Day to replace reservoir voidage due to production (Saidi 1996;Alavian and Whitson 2005) as shown in the bottom panel. Two major field shutdowns took place during 1951-1954 and 1980-1986 periods. Consistent behaviors are seen on pressure, oil production and gas injection profiles (Fig. 7), thereby reservoir pressure has accordingly increased gently after the first shutdown but sharply after gas injection.
The profile of static bottom-hole pressure surveys (top plot) indeed superimposes the profiles of pressure surveys of 50 wells across the field. It demonstrates almost similar behavior with minimal deviation, which deduces high interwell connectivity and communication due to extensive and conductive fracture network across the field.

Oil production and productivity index
Figures 8 and 9 present the bubble maps of cumulative oil production and productivity index (PI) of wells, respectively. PI is an indicator of efficiency of a well to produce fluid and is obtained by dividing production rate (bbl/day) by drawdown pressure (psi). Although, production wells have not started to produce at the same time; however, it is assumed that they have almost depleted their own drainage area after such a long time. In addition, since all wells are vertical so high well's productivity (> 10) is an indication of extension of fracture network around the wells. Comparing Figs. 8 and 9 reveals almost a consistent correlation between the cumulative oil production and productivity of wells whereby the most productive wells have produced vast majority of field production.
In addition, the most productive wells are within the fractured zones due to high reservoir connectivity and deliverability.

Production GOR data
The field GOR follows the Asmari reservoir pressure (Saidi 1996) as shown in Fig. 10. It decreases when reservoir pressure (top panel) declines at the early production times then it increases during shutdown times and re-pressurizing periods. Saidi (1996) attributed this variation to diffusionconvection process. Figure 10 presents the field GOR profile (bottom panel), whereby variation is relatively small before 1979 but suddenly jumps during 1979-1980 due to gas cap expansion and gas breakthrough into the production wells. High gas production problem has been mitigated by changing the perforation intervals afterward.
Chan's diagnostic plots (Chan 1995) are extensively used for diagnosing water production and explain mechanisms behind the excessive water production in oil wells (Firouz et al. 2019;AlOtaibi et al. 2019;Garcia et al. 2019). Although there is not much water production in this field; however, this method was used to distinguish free gas production mechanisms. Figures 11, 12, 13, 14, 15, 16 present GOR and GOR' data of typical wells on log-log scale where the blue curves shows GOR data and the red filled circles show derivative of GOR data (or GOR'). Gas coning at the early production period is distinguished by the diverging profiles (increasing GOR and decreasing GOR' trends). In contrast, gas channeling is distinguished in the middle and/ or in the late periods by the positive gradients or increasing both GOR and GOR' trends.

Pressure transient analysis
Fracture network can also be characterized by pressure transient analysis (PTA). It is often distinguished by a valley shape on the pressure derivative profile on the log-log scale (Gringarten 1984). The PTA profile of well HK-056 as a typical case is shown in Fig. 17. Pressure data are shown by the green plus signs and pressure derivative data by the red circles. Valley on the pressure derivative plot (red circles) is a signature of fracture network or dual porosity effect in the vicinity of this well.
Dual porosity model (Bourdet et al. 1989) is often used for characterizing fracture network whereby two parameters (ω and λ) are used as regressing parameters. The parameter ω is the ratio of fracture storage to total storage (sometimes called storativity contrast), and the parameter λ, which often called the inter-porosity flow represents the communication between the matrix and the fracture network. The predicted pressure profile by model is shown by red line and the derivative profile by the black line (Fig. 17).
Besides, investigations on the pressure transient data of five other wells across the field showed that pressure in these wells has built up to 90 percent of the final pressure in less than 10 min, which are additional evidences of extension of fracture network in the vicinity of these wells.

Temperature surveys
Temperature surveys often demonstrate reservoir cooling or warming. One of the potential reasons for the reservoir cooling in the conventional reservoirs is the Joule-Thomson effect (JTE). It can take place around the choke and sampling valves and creates noticeable temperature gradient as well (Steffensen and Robert 1973;App 2009). Cooling effect often takes place around the low permeability wellbores with high drawdown pressures, i.e., ~ 6000 psi, so oil mobility is reduced accordingly due to increasing oil viscosity (App 2008). Thus, temperature gradient data are important for reservoir performance study. One of the important implication of temperature gradient is evaluating potential convection process in NFRs (Peaceman 1976; Rees and Pop 2000;Ghorayeb and Firoozabadi 2000). Figure 18 shows variations of reservoir temperature during 1933-1950 and 1993-2011 periods with a significant gap during 1950-1992 period. Cooling effect is seen during both 1933-1950 and 1993-2011 periods. High drawdown pressure and releasing gas around the production wells could be a potential reason. However, high throughput of gas injection and the associated friction effect could be a potential reason for the reservoir warming effect during 1999-2001. In order to reduce uncertainties pertaining to cooling and warming effects while gas injection, only before gas injection data are used for fracture characterization purpose. Figure 19 shows the variation of temperature across the reservoir thickness between gas-oil contact (GOC) and WOC levels. The red dashed line represents the average geothermal profile and the black dashed line is the fitted curve through available data. Figure 19 shows that reservoir temperature around the GOC is somehow warmer than the nominal geothermal profile. Whereas, the average regional geothermal gradient is 1.66˚F/100ft (Heidarifard et al. 2008), the temperature gradient across the GOC is 0.46˚F/100ft. This implies rather uniform reservoir temperature across the GOC. This phenomenon can be explained by the convective motions through fracture system particularly vertical fractures. Therefore, a more uniform vertical temperature profile is established, so reservoir temperature is rather cooler at the bottom and warmer at top when it is contrasted with the nominal geothermal profile (Saidi 1996). Similar experience has been reported in Kirkuk Oilfield by Freeman and Natanson (1959). Figure 20 illustrates the areal distribution of temperature gradient records on the field UGC map. Three color codes are used to show variation of temperature gradients, so yellow means 0.3-0.7˚F/100ft, light blue means 0.7-1.0˚F/100ft and red means 1.0-2.0˚F/100ft. The low temperature gradient records are seen in the central and crestal wells where the maximum fracturing and connectivity are anticipated.
In turn, high temperature gradient records are found at the edges where the lowest fracture density is anticipated.
This important conclusion leads us to incorporate temperature gradient data into the fracture modeling workflow. As such, higher fracture connectivity and transmissibility are applied at the crestal area to address low temperature gradient also enable convection process there.
According to the above finding, a correlation between the temperature gradient and productivity index, similar to the oil production and productivity index (Figs. 8 and 9) is 1 3 expected. Figure 21 compares the temperature gradient data with the productivity index records in the same wells. A negative exponential relationship is seen, thereby the lowest temperature gradient records are seen in the most productive wells. It means, the most productive wells are in the proximity of the highly interconnected fractures where temperature gradient is minimal due to strong convective motions. Gas channeling and fracture connectivity in some of the key wells (HK-003, HK-034, HK-43, HK-044, HK-048 and HK-050) have been demonstrated recently on the GOR and GOR' diagnostic plots (Figs. 11,12,13,14,15,16). Therefore, this plot can provide a new insight about the proximity of wells to the fracture network, or how fracture network has been developed across the field.

Saturation pressure surveys
Vertical projection of the earliest saturation pressure data represents a well-mixed oil in the Asmari reservoir. Figure 22 shows the profile of 11 saturation pressure data Although production was commenced in 1928; however, production uplift was not that much to change the saturation profile substantially. Figure 22 shows that the saturation pressure records are very close to 1400 psi when reservoir temperature was at 112˚F. Negligible difference in fluid characteristic between top and bottom of reservoir is attributed to the continuous convection process through the vertical fractures. Areal projection of the above data on the reservoir top view represents a parabolic distribution as shown in Fig. 23. Higher saturation pressures are seen in the flanks, perhaps due to less depletion and higher reservoir pressure. In contrast, lower saturation pressures are seen in the middle due to higher depletion. However, saturation pressure records in the central area are fairly uniform (Fig. 23).
This implies a well-mixed reservoir fluid due to establishing high convective motions through fracture system in the central area.
In summary, reservoir fluid in the middle of reservoir is rather uniform and homogeneous; hence, the fracture network in the model should be well developed and interconnected to enable natural convection process and mixing reservoir fluid.

Mud loss data
A common problem while drilling the NFRs is severe mud loss, which often leads to the nonproductive time (Dokhani et al. 2019) and formation damage (Lietard et al. 1999). Various parameters can control mud loss rate such as rate of penetration, mud weight, viscosity, weight on  (Sun et al. 2015). Mud loss rate indeed is a dependable variable and therefore, it cannot be used for determining fracturing index solely, unless its relationship with the fracture index is verified appropriately by a supportive evidence. Unfortunately, due to shortcoming and inaccuracy of data, available mud loss data cannot be used for determining fracturing index here.

Production logging tool (PLT) data analysis
The Production Logging Tool (PLT) mostly is used to infer contribution of reservoir layers in production. PLT data were reviewed to obtain fracturing index through the Asmari reservoir. There are three sets of PLT data in three wells; however, data are not complete and/or consistent enough with other data to be used for determining contribution of Asmari's zones and characterizing fracture network. For instance, the available PLT and mud loss data do not belong to the same wells to search for any correlation.

Fracture modeling
Over the last decades, considerable efforts have been devoted to understand the physical processes occurring in the fractured reservoirs. The naturally fractured reservoirs may contain many fractures that permeate different regions of the reservoir, occurring with greater intensity near faults or folds. Fractures serve as the highly conductive flow paths for the reservoir fluid, increase the reservoir's effective permeability significantly. The fracture volume is usually small, so the porous rock matrix serves as the primary source of hydrocarbons. Understanding fault and fracture network can enhance hydrocarbon recovery process significantly through optimizing well location and directional drilling (Panien et al., 2010).
In order to incorporate fracture network in the geological model, an integrated fracture modeling study was initiated by re-interpreting 2D seismic data and incorporating image logs as well as dynamic data including oil production and productivity index data. Additionally, the results of pressure and temperature surveys, gas-oil ratio (GOR) and saturation pressure data were used as guides. The fracture  Fig. 24. The 3D fracture model was built by using FracaFlow® software then results were deployed into the static model by using Petrel® software.
Fracture sets were defined through interpreting image logs in eight wells, as mentioned earlier. Considering fracture studies in the nearby fields (Ahmadhadi et al. 2008) and geological data (Sherkati et al. 2006), the mechanism of natural fracturing in Haftkel field is different, and classified as joints or extensional fractures (tension gashes, veins). Fracture sets that are striking approximately parallel to the fold axis with coordinates of N25W-S25E, and those are striking approximately perpendicular to the fold axis with coordinates of N64E-S64W. Shear fractures exhibiting angles of approximately 30° to the local shortening axis and anticline axis (SET N85W-S85E and N52W-S52E) as shown in Fig. 25.
The small-scale and large-scale (seismic) faults are two types of faults that are considered in this analysis. Small scale faults have been reported in the middle and southeast parts of the Haftkel structure. Four small-scale normal faults with NE-SW direction have been reported previously in several wells (NISOC 2017). The presence of small-scale faults would explain the broken interval in the intra-Asmari. Although the side and geometry of small-scale faults are not still clear but well production and well test results confirm the presence of fault and fracture swarms.
The large-scale faults; however, have been detected in seismic study. The fault zone in the NW is composed of Fig. 21 Comparing PI and temperature gradient data Fig. 22 The profile of earliest saturation pressure data against the sampling depth seven parallel faults which are dying out upward the Gachsaran and are shown by black lines in Fig. 26. Moreover, in case of scatter 2D lines, the chance of coming across with whole faults is too low. However, in NW part of reservoir a corridor of normal faults was detected in the new seismic interpretation. Furthermore, possible faults based on the curvature map of top of Asmari are shown by the light blue dotted lines in Fig. 26.
Furthermore, the curvature, production index, temperature gradient, maximum production maps were used to quantify the contribution of each fracture set to flow, which then can be mapped on the reservoir basis with the more widely distributed log and seismic data. As results, fracture density driver sets in N85W-S85E, N64E-S64W, N25W-S25E and N52W-S52E directions were determined, which are shown in Figs. 27,28,29,30,respectively. The resultant set of maps then were deployed into the model for constraining the fracture model.
Next step, the parameters of discrete fracture model (DFM) such as permeability, porosity and block sizes were   Incorporating dual porosity model and calibrating the associated parameters with the test data can be considered as a reliable method to take matrix-fracture interactions into account; however, use of analog data in the untested wells imposes uncertainties in estimating fracture parameters. Therefore, the fracture modeling approach presented here is recommended for the mature NFRs, in which sufficient production history along with the test data are available to describe fracture effects on the fluid flow and reservoir characteristic. On the other hand, it is not recommended for green NFRs where production and test data are not sufficiently available to characterize fracture model thoroughly. Even though all available data have been reviewed in the current case; however, there are inadequate data for determining aperture and width of fractures. Thus, they have been used from the nearby fields. In addition, seismic attributes were not available in the entire areas (i.e., southeast of field) to control fracture distribution in the area with inadequate data.
In summary, four fracture set, which were determined by the image logs and core samples, were incorporated into the fracture model whereby their characteristics were allocated in the light of test and production data. Finally, the model was used to generate the corresponding fluid flow properties like porosity, permeability and block size of fractures.

Summary and conclusions
The Asmari reservoir in the Haftkel field is one of the most prolific naturally fractured reservoirs in the southwest of Iran where fracture network is a key factor in sustaining oil production. For the sake of characterizing fracture network; therefore, all the related data were evaluated and the following conclusions were derived: • Dominant dip direction has significantly been influenced by the active compressional and tension regimes moving from NW to SE flanks. Fracture intensity has also reduced in the same direction due to the presence of gentle structure trend. Dip direction has generally been controlled by the folding mechanism rather than fault system. • A consistent correlation between oil production and productivity of wells is seen whereby the most productive wells have produced vast majority of total field oil production; therefore, the most productive wells are within the fractured zones. • Comparing the static bottom-hole pressure surveys of 50 wells across the field demonstrates almost similar behaviors, which deduces high inter-well connectivity and communication due to extensive and very conductive fracture network across the field. • Gas coning followed by the gas channeling signatures are seen on the GOR-GOR' plots of most of the productive wells, which represent gas production through the fracture network in the central area. • The profile of temperature data across the reservoir thickness is rather uniform, particularly across the GOC due to strong convective motions through vertical fractures. In addition, low temperature gradient records are seen in the central and crestal wells where maximum fracturing and connectivity are anticipated. In turn, high temperature gradient records are found at the edges where the lowest fracture density is anticipated. Additionally, a correlation is seen between the temperature gradient and productivity index of wells, which provide insights about the proximity of wells to the fracture net- work as well as extension of fracture network across the field. • Negligible difference in the initial fluid characteristic across the reservoir thickness is attributed to continuous convection process through the vertical fractures. The uniformity of saturation pressure records in the central area implies a well-mixed reservoir fluid due to establishing convective motions through the fracture system in this area. • Fracture modeling was initiated by identifying the largeand small-scale faults using 2D seismic lines and well data, respectively. Furthermore, image logs, production, productivity index, pressure and temperature surveys were used to determine fracture sets as well as their pattern. After determining fracture density driver sets, the parameters of discrete fracture model (DFM) such as permeability, porosity and block sizes were calibrated in the light of well test results. • Fracture modeling can be improved by running further logs and tests to quantify aperture, width and conductivity of fractures particularly in the undeveloped areas. Furthermore, a seismic acquisition to supply seismic attributes in the southeast of field would be a great opportunity to address uncertainties of large-scale faults.