Determine the Reservoir Characterization, Using Seismic Data, Well Logs, Attribute Interpretation, and Static Modelling in the Duhok Province, Kurdistan Region, Iraq

Abstract


Introduction
Today, the global demand for oil and gas rising continuously, and meeting these needs will require the discovery of new hydrocarbon reserves.The northeastern margin of the Arabian plate (Zagros-Taurus Folds and Thrust belt) extends from northern Iraq to Oman and forms one of the most fertile regions for recoverable hydrocarbon reserves in the world in terms of global hydrocarbon reserves.(Beydoun, 1998;Bretis et al., 2011).The area under study in the "WN" oil field forms part of this belt, which is located in Duhok Province to the northwest of Duhok City within the High Folded Zone in the Kurdistan Region of Iraq (Fig. 1).The belt was formed and affected by Alpine Orogeny during the closure of the New-Tethys Ocean between the Arabian and the Eurasian plates and the subsequent continental collision between these plates which was started at the Late Cretaceous and reached its culmination at the Pleistocene (Strevanovic and Markovic 2003;Sissakian and Al-Jiburi, 2014).In Iraq, more than 99% of the discovered oil reserves exist in the Cretaceous and Tertiary hydrocarbon reservoirs (Aqrawi et al., 2010).Worldwide, carbonate rocks are regarded as a significant host rock for oil and gas reserves (Durrani et al., 2020).The "WN" oil field consists of different reservoir horizons and among them, horizon No.5 is considered for the current study, which forms an important Tertiary carbonate reservoir in the other oil fields of the region.Compared with other reservoirs, due to heterogeneity and pore structure complexity, modelling and characterization of the carbonate reservoirs are more complicated and challenging for geoscientists (Bueno et al., 2014;Tawfik et al., 2019).Evaluating those heterogeneities may be accomplished using interpretation of seismic data, seismic attributes and welllog data, which are crucial for building accurate geologic models (Bueno et al., 2014;Khalid et al., 2015).Integrating seismic data interpretation and well log data analysis have been widely employed in significant studies over the past several decades and have been successful in the evaluation of hydrocarbon reservoirs around the world (e.g., Eshimokhai and Akhirevbulu, 2012;Abe and Olowokere, 2013;Niri, 2015;Fajana, 2020;Waziri et al., 2022), and in Iraq, (Abdulateef, 2018).In this study, the computer-processed interpretation of well logs was used for the two boreholes in the WN oil field.(Khawaja et al., 2023) used seismic attributes as criteria for direct hydrocarbon indicators (DHI) to identify the presence of hydrocarbons.Envelope (reflection strength) and root-mean-square (RMS) amplitude of the surface seismic attributes computed from the seismic data, high amplitudes anomalous suggest the presence of hydrocarbon and improve reservoir characterization.There is no previous petrophysical analysis, 3D modelling, or any published academic studies using 3D seismic survey data and well logs regarding this reservoir.Thus, this work is aimed towards the determination of the reservoir characterization by combining various data types including a 3D seismic cube, well logs data, check shots of four wells and other techniques such as seismic interpretation, attribute computation, 3D static modelling and 3D property modelling.

Regional Tectonic Setting and Geological Formations
Tectonically, the "WN" oil field is located in the High Folded Zone (Fig. 1) within the unstable shelf as a part of the Arabian platform.The basic structural framework of the High Folded Zone is mainly the result of the closure of the Neo-Tethys Ocean and the subsequent continental collision between the Arabian and Eurasian (Turkish and Iranian) Plates (Jassim and Goff, 2006).In general, this zone is characterized by high-amplitude anticlines and narrow synclines, the anticlines cored by Mesozoic carbonates with outcropped Cenozoic clastic on their limbs.They are trending almost in the NW-SE direction parallel to the Zagros structures and swinging trend here to the E-W direction parallel to the Taurus Belt in the southern part of the Turkey Plate (Jassim and Goff, 2006).The outcrop geological formations from older to younger age include Gercus (M Eocene) consisting of mudstones, sandstone, sandy and gritty marl; Avanah Limestone (M-U Eocene) is composed of limestone, usually dolomitized and recrystallized; Pila Spi (M-U Eocene) is formed of bituminous dolomitic and chalky limestone; Al-Fatha (M Miocene) is composed of anhydrite, gypsum, marl, and limestone; Injana (U Miocene) consists of thin-bedded sandstone and claystone; Al-Muqdadiya (L Pliocene) is composed of sandstone, mudstone, and siltstone; Bai Hassan (U Pliocene) consists conglomerate, siltstone, claystone, and sandstone (Bellen et al., 1959).A few parts of the area are covered by surficial Quaternary alluvial, terraces and slope deposits (synclinal trough filling of polygenetic origin).

Materials and Methodology
The dataset which is used in this study includes an SEG-Y format 3D seismic cube that consists of 536 inlines (3500-4035) and 491 crosslines (2150-2640), with line spacing of 12 m covering a total area of about 37 km 2 and well logs data with their four wells of check shots survey.The seismic lines and oil well locations have been plotted by which the base map generated as shown in Fig. 3. Schlumberger's Petrel software (version 2017.4) and Interactive Petrophysics (IP) software (version 2018) have been applied for the 3D seismic interpretation and geological model of the studied horizon and well logs analysis All data used in this study was obtained from International Oil Companies (IOC) working in the Kurdistan Region under licenses granted by the Ministry of Natural Resources of the Kurdistan Regional Government and with the Ministry's permission The seismic data, well tops, check shots and well logs of the four wells were loaded into Petrel software, and then, the seismic data was realized for fast handling.The seismic data condition was achieved by applying some post-stack filters and volume attributes such as frequency, automatic gain control (AGC) and structural soothing filters to improve signals and remove noise from seismic data even to aid structural interpretation.A synthetic seismogram for the well WN.1 was generated based on a calibrated sonic log, density log, check shots and an appropriate wavelet extraction from the stacked seismic traces beside the well.Identifying and locating the top and bottom of the target reservoir (horizon-5 and horizon-6) on the seismic sections was performed by a well-to-seismic tying process (Fig. 4 and Fig. 5).Picking horizons for the whole area has been done manually by conducting a threedimensional horizon interpretation at a regular interval of all fifth lines in both inlines and crosslines, while one-by-one line intervals in the fractured and weak reflector zones.Manually, nineteen minor faults of different sizes have been identified and determined.Picking of the top and bottom of the target reservoir throughout the area has been conducted, and surfaces of the picked two-way travel time (TWT) values have been made by using them to build a 3D static grid model.The bottom of the reservoir is not necessary for interpretation anymore which is just used in the construction 3D model, thus, only the isochron map for the top of the reservoir has been created (Fig. 6).Many methods exist for the depth conversion process, one of the popular ways is creating the simple velocity model.A layer cake velocity model has been selected using the equation V= V0 + K (Z-Z0), which has been implemented by many researchers (e.g., Qadri et al., 2017;Osinowo et al., 2018).Where V0 is the velocity (reference velocity) at the top of the layer at the reference depth Z0, K is the velocity gradient (compaction factor) while Z is the calculated depth.
This velocity model has been applied to generate a depth map of the top of the reservoir (Fig. 7) using the isochron (TWT) surface map, and its corresponding well top.In addition, an isopach map between the top and bottom of the reservoir has been constructed (Fig. 8).The 3D static grid model (Fig. 9) has been built to indicate the geometry of the reservoir and determine the distribution and visualization of the petrophysical properties across the 3D geological model based on fault and horizon modelling.The model involves incorporating all of the minor faults, the top and bottom of the reservoir depth maps, their corresponding well tops and the isochore thickness map (constructed just for zonation) between its top and bottom.In the next step, its unique zone (reservoir thickness) has been divided into a hundred layers.The petrophysical analysis has been performed for the available composite logs in the wells WN.2 and WN.3 applying IP software to derive and estimate the essential parameters such as effective and secondary porosity, permeability, clay volume (VCL), and water saturation (SW).Computer-processed analysis for well logs from the same two wells (Fig. 10 and Fig. 11), Neutron-Density cross plot, M-N cross plot and Buckles plot have been done to check well logs and determine the lithology and porosity of the reservoir, which are shown in Figs. 12, 13, 14, 15, 16 and 17, respectively.The extracted properties from the petrophysical (well-logs) analysis have been upscaled to the scale of the constructed 3D static grid model.Then, these upscaled properties were inserted into the grid skeleton to populate throughout the reservoir, and in such a way 3D property model for each of the essential petrophysical parameters was constructed, as presented in Figs. 18, 9, 20, 21, and 22, respectively.Seismic attributes are any quantifiable aspect extracted from seismic data (Roden et al., 2015).Seismic attribute interpretation has been performed on the isochron map (top of the reservoir surface) using both the root mean square (RMS) amplitude and the envelope surface attributes to produce the attributes maps (Figs.23 and 24) and to detect hydrocarbon presence.

Results and Discussion
The seismic wavelet is the most required element for synthetic seismogram generation.Wavelet extraction is the most crucial link between seismic and well data in the context of forward seismic modelling (Bisaso, 2011).The extraction was performed by applying the statistical technique wavelet extraction algorithm which then properly fit between the traces and seismic data was obtained (Figs. 4  and 5).Seismic data-based statistical wavelets provide best-fit models representing the typical behaviour of actual parameter characteristics (de Morais, 2018).So, any stretching or squeezing of the synthetic seismogram was avoided through the seismic-well tie procedure and showed a perfect matching.Depth conversion is one of the most crucial steps in seismic interpretation and reliable geological modelling.The layer cake velocity model for time-depth conversion has been used which extracts the parameters from the well logs data for each zone.The moving average algorithm for interpolation has been run to interpolate and distribute the velocity across the model.This algorithm calculates the average of the input data and adjusts the weighting based on the distance away from wells for every unsampled area (Schlumberger, 2010).The best fit was made between the depths of the reservoir in the produced depth map and its depths at the wells (Fig. 9b).The isochron (TWT) map in Fig. 6 shows contour line values with an interval of 50 milliseconds (ms) ranging from 700 to 1350 ms.The depth map in Fig. 7 contour lines exhibits relatively varied values with a contour interval of 50 m ranging from 1050 to 1750 m.At the uppermost of the horizon, the crest of the anticline occurs at 1011m and spreads towards the lowermost part of 1786 m on the northern flank which suggests a hydrocarbon structural prospect.The interesting contour lines in both isochron and depth map of the reservoir show a gradual increase in twoway time and depth towards the north and south of the area revealing an asymmetrical subsurface structural closure that is depicted as a double-plunging rollover anticlinal structure.Numan (2000) suggested that listric faults were created as normal listric faults to accommodate the rifting and extensional tectonic regime during the detachment of the Iranian and Turkish plates from the Arabian Plate due to the opening of the Neo-Tethys Ocean in the Triassic.
During the Cretaceous, the Neo-Tethys began to close, the Arabian Plate subducted beneath the Iranian and Turkish plates, and subsequent continental plate collision gave rise to an inversion of the regional tectonic regime that changed from extensional to compressional tectonism (Numan, 2000).The hanging-wall block of listric normal faults could buckle as it adapts to the curving fault surface, which might lead to anticlines, referred to known as rollover anticlines that form important traps in the hydrocarbon fields (Suppe, 1985;Tearpock et al., 2020).It is commonly accepted that the folds in the Taurus-Zagros Belt within northern Iraq are decollement buckle folds (Ameen, 1991).The length and width of this anticline are about 6.4km and 3.6km, respectively.However, it has an east-west trend axis which suggests north-south compression and coincides with the general structural series of the Taurus fold-thrust orientation in the region.But its southern limb appears gentler than the northern limb which is inconsistent with the common Taurus structural series.This conflict could be probably interpreted as the result of the presence and effect of a hidden deepseated fault (Ameen, 1991).Basement faulting is responsible for the asymmetrical anticlines, their steeper limbs are likely to occur directly over the basement faults.The northern limb of this anticlinal closure is bisected by nineteen of the minor reverse faults with small throws.Fault throws range between 4.5 and 43 m and fault lengths between 170 m and 2.75 km.Most of these faults are trending in the W-E directions parallel to the anticlinal axis, and the resets are trending toward the WSW and ENE directions which are perpendicular to the anticlinal  and 7).The pre-existing normal faults, known as the foreland listric fault were reactive and change to reverse (thrust) faults due to the tectonic inversion from extension phase to compression (Numan, 2000).The current structural pattern in the study area forms part of the Zagros-Taurus foldthrust belt which developed and resulted from the opening and closure of the Neo-Tethys Ocean and the progressive subsequent continental collision between the Arabian and Eurasian (Iranian and Turkish) plates from the Permo-Triassic period to the present (Sharland et al., 2001;Jassim and Goff, 2006).This type of anticlinal structure and the forming of these minor reserve faults on its northern flank has been suggesting a response to this tectonic system of folding, faulting, and thrusting on the originally inherited normal listric faults along the older obduction zone during the prevailing compressive (convergent) phase.
The isopach map (Fig. 8) displays a spatial variation of thickness values of the reservoir with a contour interval of 20 m, which has a relatively variable range from 160 to 330 m with an increase of values in some places as indicated by the blue colour.Its maximum thickness within the area is concentrated at the location of the well WN.4.This variation might suggest changes in sedimentation or removal after sedimentation.A structural framework (Fig. 9) forms the basis of the 3D static (geological) modelling, which provides borders to reservoir geometry, allows distribution properties throughout it, and offers a practical method for reservoir heterogeneity.
In other words, it precisely represents the subsurface conditions of the reservoir that the wells have encountered (Ojo et al., 2018).To represent the tiny vertical heterogeneity in a suitable level of detail, by layering process, the reservoir was divided into a hundred layers.The thickness of each layer is around 10m on average along the vertical small-scale of the wells to show a more precise spatially populating the reservoir properties.
The computer-processed interpretation (CPI) of wells WN.2 (Fig. 10) and WN.3 (Fig. 11) display the relationship of the varied petrophysical properties such as lithology, calculated effective and secondary porosity, permeabilities, water saturation (SW) and clay volume (VCL) with gamma-ray (GR), calliper and the depths of the wells.The neutron-density cross plot for the well WN.2 (Fig. 12) and well WN.3 (Fig. 13) for the reservoir show that the majority of plotted data points tend to be clustered close and aligned to the dolomite trend line and a few of the points tend to scatter around the limestone trend line.The M-N cross plot for well WN.2 (Fig. 14) and well WN.3 (Fig. 15) depicts the composition of the reservoir, the dolomite trend line and part of the calcite trend line were blurred and obscured by the majority of the plotted cluster points confirming that the reservoir consists of dolomite and dolomitic limestones.Buckles plot is a hyperbolic relationship between porosity and irreducible water saturation made up of a series of hyperbolic curves of fixed Bulk Volume of Water (BVW).The plot can be used as a direct technique for determining the bulk volume of the irreducible water saturation zone (Jumaah, 2021).The majority of the points in the well WN.2 cluster between the BVW hyperbolic curves (blue lines) 0.02 and 0.04, and in the well WN.3 cluster between the BVW hyperbolic curves 0.02 and 0.06.Tracking the majority of points around the BVW 0.04 curve indicates that the reservoir zone is at some irreducible water saturation perspectives and oil production with little or without water (Hussain, et al., 2022;Rashid et al., 2023).
The effective porosity model of the reservoir in Fig. 18 shows the spatial distribution of the effective porosity which has a minimum value of 0.97%, and a maximum value of 23%. with an average value of 14%.suggests good reservoir potential for this horizon.The 3D secondary porosity model in (Fig. 19) illustrates the variable distribution of secondary porosity values which ranged between 0.43% and 11.7% and has an average of 6.5%.The knowledge of permeability distribution is the key parameter to reservoir evaluation (Soto et al., 2001).The permeability in the 3D model (Fig. 20) shows the spatial distribution of the permeability that ranged between 0.00 millidarcy (mD) and 238.8 mD and has an average permeability of 18.9 mD.The reservoir has heterogeneous permeability, with the blue, green, and yellow colours dominating the highest values and pinkish the lowest values.Permeability from the Buckles model, in the well WN.2, the values of permeability (Fig. 16) are observed to be more than (1000) millidarcy (mD), while permeability (Fig. 17) in the well WN.3 reaches up to (500) mD.Carbonate reservoirs are complex, extremely heterogeneous, and have limited porosity and permeability (Jassam and Al-Fatlawi, 2023), and organizing carbonate reserves is a challenge due to their compact and heterogeneous nature (Faidhllah and Hamd-Allah, 2023).In the same reservoir with the same porosity, both high and low permeability zones can be present (Salman et al., 2022).Computing the porosity and permeability for two samples of the cored interval in one well and has been found that they have the same porosity value of 15.8% with different permeability values of 75.8 mD and 858 mD (Dou et al. 2011).Permeability computed for the carbonate packstones (Holocene) from Florida and the Bahamas, showed ranges from 31.5 and 9300 mD (Enos and Sawatsky, 1981).Makhankova et al. (2020) studied a carbonate rock (M.Miocene) in Luconia, South China Sea, and found extraordinarily high permeability of up to 3,500 mD with an average of 800 mD.The permeabilities of reservoirs can be broadly classified (Tiab and Donaldson, 2015) as follows: poor if k < 1mD; fair if 1mD < k < 10 mD; moderate if 10 mD < k < 50 mD; good if 50 mD < k < 250 mD; very good if k>250, thus, accordingly this reservoir with average permeability of 18.9 mD has corresponded with the moderate type reservoirs.The clay volume (VCL) 3D model of the reservoir (Fig. 21) reflects the clay volume value ranges from 0.00% to 9.8% and the dominant lithology is dolomite-rich carbonate rocks characterized by a very low average volume of clay that hardly reached 1.6%.
The 3D water saturation (SW) model (Fig. 22) reveals water saturation distribution within the reservoir varies from 12.44% to 99.8% with an average of 33.67%.The average hydrocarbon saturation of the reservoir can be computed using the equation Sh = 100 -SW% (Sanuade et al., 2017) where Sh=hydrocarbon saturation; SW=water saturation Thus, the average hydrocarbon saturation will be equal to 66.33%.This saturation (Sh) average and the results of the Buckles plots in Fig. 16 and Fig. 17 suggest that the reservoir might be more hydrocarbon-productive than water.
Root-mean-square (RMS) amplitude and envelope amplitude attributes can be run to check the existence of hydrocarbons and aid in reservoir characterization (Thota et al., 2022).The envelope (reflection strength) attribute, also called the instantaneous amplitude attribute is one of the complex trace attributes (Taner et al., 1979).Root-mean-square (RMS) amplitude attribute computes the square root of the sum of squared amplitudes divided by the number of samples within a defined window used (Alotaby, 2015).Kalu et al. (2020) used envelope and RMS amplitude seismic attributes in their research as hydrocarbon direct indicators (HDI) tools.In this study, both the envelop and RMS attributes have been extracted from the reservoir horizon surface that exhibited similar amplitude characteristics.The brighter amplitude reflection can be seen on the both envelope and RMS seismic attribute maps in Fig. 23 and Fig. 24 respectively.

Conclusions
The current study has been carried out using the integration of 3D reflection seismic data, well logs, and check shots, with attribute interpretation, petrophysical analysis, cross plots and 3D property models to characterize the reservoir at the WN oil field.Isochron and depth contour maps of the reservoir are successful in identifying and confirming the presence of asymmetrical double-plunging subsurface rollover anticlinal structure with length and width of about 6.4km and 3.6km that forms a favourable structure for a hydrocarbon trap.It has an axis trending to the west-east direction parallel and corresponds with the Taurus Mountain series, while its northern flank is steeper than its southern flank, which is discordant with the common structural pattern in this region.This incompatibility might reveal the existence of a concealed deep-seated fault.A total of nineteen minor reverse faults with throws ranging between 4.5 and 43m and lengths between 170m and 2.75km striking the northern limb of the anticline structure in E-W and ENE-WSW directions have been determined and then used  for building and gridding 3D modelling.The presence of these reverse faults suggests that this rollover anticline is subjected to compressional stress causing inversion in the faulting system in this area.Different cross-plot analyses allowed for determining and computing lithology, porosity, permeability and fluid content for the reservoir, which enhanced the interpretation.The average thickness of the reservoir is 249m, the effective porosity ranges between 0.97-23%, the secondary porosity between 0.43-11.7%, the permeability between 0.00-238.8mD, the water saturation between 12.44-99.8%,and the clay volume between 0.00-9.8%,with the average of 14%, 6.5%, 18.96mD, 33.67%, and 1.6%, respectively, and the hydrocarbon saturation is 66.33% suggesting a potential hydrocarbon reservoir.Envelope and root-mean-square (RMS) amplitude surface seismic attributes results successfully revealed high amplitude areas that identify the presence of hydrocarbon zones and improved reservoir characterization.3D property models were beneficial to confirm that the reservoir is of moderate quality reservoir.

Fig. 1 .
Fig. 1.Tectonic map of the Kurdistan Region of Iraq (from Abdula et al., 2015) showing the location of the study area

Fig. 2 .
Fig. 2. The geological map of the study area and its surroundings (from Bamerni et al., 2021)

Fig. 3 .
Fig. 3. Base map showing the seismic lines geometry and well locations

Fig. 4 .Fig. 5 .
Fig. 4. Synthetic seismogram for the well WN.1, shows identifying the top and bottom of the reservoir that correlated with the original seismic section of the inline No. 3876

Fig. 6 .
Fig. 6.Isochron (TWT) map of the top of the reservoir (H-5), minor faults positions, and well locations

Fig. 7 .
Fig. 7. Depth map of the top of the reservoir (H-5), minor faults positions, and well location

Fig. 8 .
Fig.8.Isopach map of the gross reservoir thickness throughout the area with well locations

Fig. 9 .Fig. 11 .Fig. 12 .
Fig. 9. 3D skeleton of the reservoir, (a) The side view shows its well top at the wells and fault patches.(b) The front view displays a depth map of the top of the Horizon-5 (Reservoir top), its well top at the wells, and fault patches

Fig. 13 .
Fig. 13.Neutron-density cross plot of the reservoir zone at the well WN.3 Figs. 16 and 17 represent Buckles plots for the reservoir in the wells WN.2 and WN.3 respectively.

Fig. 14 .
Fig. 14.M-N cross plot for the reservoir zone in the well WN.2

Fig. 15 .
Fig. 15.M-N cross plot for the reservoir zone in the well WN.3

Fig. 16 .
Fig. 16.Buckles plot of the reservoir zone for the well WN.2.The grey lines represent lines of constant permeability and the blue lines represent the constant values of the product of p

Fig. 17 .
Fig. 17.Buckles plot of the reservoir zone for the well WN.3.The grey lines represent lines of constant permeability and the blue lines represent the constant values of the product of porosity

Fig. 20 .
Fig. 20.3D permeability model for the reservoir

Fig. 21 .
Fig. 21.3D clay volume CLV model for the reservoir

Fig. 23 .
Fig. 23.Surface envelope amplitude attribute map of the top of the reservoir

Fig. 24 .
Fig. 24.Surface RMS amplitude attribute map of the top of the reservoir