Exploration of the Brookian sandstone reservoirs in the Nanushuk and Torok Formations on the North Slope of Alaska is a hot topic and presents opportunities to the oil and gas community because of their shallow depth, vast extent, and scope of development. The consecutive hydrocarbon discoveries announced by Repsol-Armstrong, Caelus Energy, and ConocoPhillips in 2015, 2016, and 2017 have indicated the presence of the vast recoverable resources on the North Slope in the Nanushuk and Torok reservoirs. We have investigated the detailed geophysical and petrophysical characteristics of these reservoirs. Our goal is to detect dominant geologic features in these formations using a combination of seismic attributes at the regional scale and analyze critical petrophysical and rock physics properties to evaluate formation heterogeneities and identify the reservoir targets by integrating well log and core data at the well scale. The Nanushuk Formation is expressed as topset reflections, whereas the Torok and gamma-ray zone formations are expressed as foresets and bottomsets on the seismic reflection data. Using seismic attributes, we mapped the extent of different geomorphological features, including shelf edges, channels, slides, and basin-floor fans, all with significant amplitude anomalies. The shelf edges continue for tens to hundreds of miles along the north/northwest and east–west directions, depending on the areas. The internal characters of these formations delineated by conventional well logs and advanced petrophysical analysis reveal their vertical heterogeneities and complexities, in terms of reservoir properties. We conclude that the reservoirs are vertically and laterally heterogeneous. These are thin-bedded low-resistivity reservoirs. Only a few zones in the parasequences are oil-saturated. We find that a combination of low VP/VS ratio and acoustic impedance can be a useful proxy to detect the hydrocarbon-bearing sand intervals in these formations.

The North Slope of Alaska has been the primary target for oil production in the state since 1977. Although the deep sedimentary formations (such as the Kekiktuk, Lisburne, Ivishak, Sag River, and Kuparuk Formations) have served as major hydrocarbon-producing reservoirs, the volume of the remaining hydrocarbon in them in the existing oil fields is significantly lower now (Bhattacharya and Verma, 2019). This has either forced the companies to cut expenditures on long-term exploration in Alaska and sell their assets or explore previously overlooked shallow hydrocarbon reservoirs in the state. However, the consecutive discoveries (Pikka-Horseshoe, Smith Bay, and Willow) of the hydrocarbon reservoirs in the Cretaceous Nanushuk-Torok Formations in the National Petroleum Reserve in Alaska (NPRA) announced by Repsol-Armstrong, Caelus Energy, and ConocoPhillips in 2015, 2016, and 2017 have demonstrated their massive exploration potential, which has attracted several oil and gas companies and investors around the world (Decker, 2018). A recent study by Houseknecht et al. (2017) has estimated that there are at least 8.7 billion barrels of oil and 25 trillion cubic feet of natural gas reserves in the Nanushuk and Torok Formations on the North Slope, which may be able to reinvigorate the declining oil output from the resource-rich state for some time. Based on the statistics, the deep Ellesmerian reservoirs (i.e., Kekiktuk, Lisburne, Ivishak, and Sag River) contributed to 74% of total hydrocarbon production in the early 1980s–1990s, and the Brookian reservoirs had close to 0% production, which has changed to 48% for the Ellesmerian reservoirs and 12% for the Brookian reservoirs (including the topsets and turbidites) in 2017 (Gregersen, 2019). Therefore, a detailed geophysical and petrophysical characterization of these emerging Nanushuk-Torok reservoirs is needed prior to oilfield development.

The Nanushuk Formation is a thick fluvial-deltaic-shelf succession, whereas the Torok Formation is its basinward equivalent. These two formations are genetically related and composed of clastic rocks. Hydrocarbons are held in these clastic rocks predominantly by stratigraphic traps. Over the years, several studies have focused on these formations using 2D seismic lines, limited wells with conventional well logs (mostly gamma-ray), core data, and outcrop studies (Houseknecht et al., 2009; Frierson, 2013; LePain et al., 2017; Ramon-Duenas et al., 2018; Houseknecht, 2019a). Although most of these studies provided geologic insights to sedimentology, stratigraphy, and depositional processes, they lack detailed geophysical and petrophysical characterization of these formations from the perspectives of the reservoir distribution, properties, and heterogeneities. Recently, Houseknecht (2019a) shows basic seismic stratigraphic interpretations on a proprietary 3D seismic survey. Helmold (2016) publishes a paper on the reservoir quality using petrographic thin sections and well logs. However, detailed quantitative interpretations and 3D architecture of different geologic features in these formations are missing from the published literature. Few studies focused on the hydrocarbons aspect from the petrophysical and rock physics standpoints.

The main goals of this study are to identify and interpret the dominant and prospective geologic features present in the Nanushuk and Torok Formations and analyze their reservoir properties using multiple high-resolution 3D seismic surveys, conventional and advanced petrophysical logs, and core data. We selected the study areas strategically based on the access to publicly available multiscale subsurface data and their proximity to the hydrocarbon discovery fields (such as the Pikka-Horseshoe and Smith Bay discoveries). In this study, we used a top-down approach to identify the dominant geologic features in these formations using 3D seismic data and seismic attributes at the regional scale followed by core- and well-log-based petrophysical analysis and rock physics studies to understand their complexities and possibly detect hydrocarbons at the well scale.

Our study is novel and complementary to the existing studies from multiple aspects. First, the study reveals the presence of different types of stratigraphic features and their 3D distribution at the regional scale identified using several 3D seismic data sets, which would be important for a better understanding of the depositional processes and oilfield development. This is one of the few examples of the application of advanced seismic attributes, such as coherent energy, Sobel-filter similarity, and spectral decomposition on the North Slope to delineate important stratigraphic features (such as basin-floor fans, shelf edges, and channels), elucidate their complexities, and infer the potential distribution of the sand and shale sequences. We think the petrophysical modeling and rock physics relations, obtained by this study, will be crucial to identify the complex hydrocarbon reservoirs when enough core samples are not available. Joint petrophysical inversion modeling, followed by acoustic velocity-based rock physics analysis, has never been published on these formations. Although we do not show a detailed regional variation of rock properties (such as rock physics) in the Nanushuk and Torok Formations, we think the rock physics relationships, obtained by this study, can be used for predicting the prospective sandstone reservoirs in these formations using seismic inversion techniques. We conclude that this study will provide an integrated quantitative workflow that can be used to identify and map these complex formations and delineate potential reservoir targets in the future. Our workflow can be useful to other areas outside Alaska as well, such as the Gulf of Mexico, Permian Basin, and Northwest Shelf of Australia.

Our study area is from the NPRA, northern Alaska (Figure 1). This is a large area and covers a significant portion of the Colville foreland basin. The basin is bound by the Brooks Range in the south and Barrow Arch (a rift shoulder) in the north. In terms of tectonostratigraphy, the main tectonic episodes on the North Slope can be summarized, to some extent, by four major sequences including the Frankilinian, Ellesmerian, Beaufortian, and Brookian sequences (Alvey et al., 2008; Grantz et al., 2011; Helwig et al., 2011) (Figure 2). The Nanushuk and Torok Formations are part of the Brookian sequence. These rocks are Early to Late Cretaceous and Cenozoic progradational deposits, mostly composed of sand and shale sequences. Houseknecht et al. (2009) and Houseknecht (2019a, 2019b) suggest that the voluminous sediment in the Nanushuk-Torok Formations was derived from the Chukotka highlands in the southwest and Brooks Range orogens in the south. According to Houseknecht (2019a), sediment derived from the Chukotka highlands was routed to the east and northeast, whereas sediment from the Brooks Range was routed to the north, filling a significant portion of the Colville Basin. The Nanushuk Formation transitions seaward to the slope and basinal facies in the coeval Torok Formation (Bird and Andrews, 1979; 0; Houseknecht, 2019a). Sherwood et al. (2002) and Houseknecht (2019a) document a change in the clinothem character in the basin. In the western part of the basin, predominantly shallow marine to nonmarine topsets, with minor foresets, have been found to thin westward and overlap against the Chukotka highlands. In the east-central part, the clinothem displays the progradational nature of the shelf margins, with topsets, foresets, and bottomsets. The Nanushuk Formation consists of sandstone, siltstone, mudstone, and subordinate coal representing facies deposited in the alluvial plain, fluvial, deltaic, and other environments (LePain et al., 2009). The shelf-edge deltas of the Nanushuk Formation are considered as major hydrocarbon-bearing facies. The Torok Formation includes shale, silty mudstone, and sandstone that comprise facies deposited in the outer shelf, marine slope, and basin-floor fan environments (Houseknecht et al., 2009).

Although drilling in the Brookian sequence, including turbidites and topsets, is a hot topic now, it started back in 1946 with the U.S. Navy. The U.S. Navy drilled 13 wells in the Brookian sequence in the Umiat field in 1944–1952. After a hiatus in drilling in the Brookian for a few years, the early 1970s saw an uptick in the drilling activities in this sequence. The Tarn and Meltwater turbidite fields were discovered in the early 1990s. The Nanuq field, discovered in 2000, is a basin-floor setting (Torok Formation). The Qannik field was discovered and put in production in 2006. It represents the topset component of the progradational sequence. A few wells in the Schrader Bluff Formation were drilled and developed in 2004–2010. And then, with the new 3D seismic and drilling technology, new discoveries (Pikka-Horseshoe, Smith Bay, and Willow) have been made in the Brookian sequence since 2013. Since 2012, more than 35 exploration wells have been drilled in the Brookian topsets and turbidites. A few wells targeting the Brookian sequence, such as the Grizzly 1, Heavenly 1, Kokoda 1, Kokoda 5, and Winx 1, were found to be either dry or not containing movable hydrocarbons (Houseknecht, 2019a; Edison Investment Research, 2019). With all of the new exploration, discoveries, ongoing appraisal, and very limited development, the oil and natural gas liquids production from the Brookian reservoirs is close to 60 MBBL/day (Gregersen, 2019).

We used four 3D seismic surveys (Nanuq South 3D, Harrison Bay 3D, NE-NPRA 3D, and Puviaq South 3D) and several wells (Grizzly 1, Horseshoe 1, Itkillik River Unit 1, Cronus 1, Flat Top 1, South Harrison Bay 1, Cape Halkett 1, and Amaguq 2) with well logs in this study (Figure 1). The 3D seismic surveys cover a large area of approximately 1201  mi2 on the North Slope. The 3D seismic surveys are widely distributed in the NPRA, so we expect that they would capture many heterogeneities in the Nanushuk and Torok Formations at a high resolution. Although 2D seismic surveys (acquired back in the 1970s) are available throughout the NPRA, we did not use those data sets because of their limited resolution to offer ideas about the 3D distribution pattern of the geologic features and perform geophysically meaningful quantitative analysis. It is worth noting that ice lake static corrections were applied to the 3D seismic data sets because tundra lakes and permafrost generate high-magnitude and high-frequency static anomalies. Table 1 summarizes the specific features of the 3D seismic surveys used in this study (Alaska Geologic Materials Center, 2019).

In this study, we used poststack data for all the 3D seismic surveys for seismic horizon mapping and attribute computations. In addition, we used near-, mid-, and far-offset stack data for the Nanuq South 3D survey to show amplitude differences. We used conventional well logs (gamma-ray, resistivity, bulk density, neutron porosity, photoelectric, and P-wave sonic) as available from the wells, all of which are either inside or nearby the 3D seismic surveys in this study. In addition, we used advanced petrophysical logs such as the pulsed neutron spectroscopy (PNS), nuclear magnetic resonance (NMR) spectroscopy, and sonic scanner (providing compressional and shear sonic logs) available from some wells to better understand the reservoir complexities. The Horseshoe 1 in the Pikka-Horseshoe area (in the vicinity of the Nanuq South 3D survey) and CT-1 and CT-2 wells in the Smith Bay area (on the western side of the NE-NPRA 3D survey) are considered hydrocarbon discovery wells (GEOExPro, 2017). Although well data from the Horseshoe 1 well in the Pikka-Horseshoe area are publicly available and we used them in this study, the data from the CT-1 and CT-2 wells in the Smith Bay discovery area may be released in the future.

We observed the different geologic features corresponding to the Nanushuk and Torok Formations in different 3D seismic surveys. To study all the different geologic features that could be present in these formations, we used several 3D seismic surveys for horizon mapping and seismic attribute analysis. The use of several 3D seismic surveys provides a higher resolution regional picture compared to publicly available 2D seismic lines. To start our study, first, we identified the important geologic reflectors in all of the 3D seismic surveys and prepared the corresponding time-structure maps (horizons). We did this by performing well-seismic ties to identify the key formations (Figure 3). We extracted statistical wavelets from the 3D seismic data and used them in generating the synthetic seismograms. We applied reasonable time shifts and minimal stretch-squeeze to match the real seismic trace with the synthetic trace, constructed using sonic and density logs. The identification of the geologic horizons was easy because it is well-known that the Nanushuk Formation is expressed as the topset, whereas its time-equivalent Torok Formation is expressed either as the foreset or bottomset on the seismic reflection data. Figures 4, 5, 6, 7, and 8 show the clinoforms mapped and interpreted on the 3D seismic data sets.

Next, we computed a variety of seismic attributes from the original seismic surveys (poststack time-migrated data) and extracted them on the prominent horizons. Seismic attributes are quantities derived from original seismic volume using different mathematical equations for structural and stratigraphic analysis and rock property estimation. In this study, we used geometric and amplitude-based attributes to identify geologic features of interest, such as shelf edges, channels, and basin-floor fans. After testing an ensemble of seismic attributes and extracting them on a few mapped geologic horizons in the Nanushuk-Torok interval, we found that only a subset of attributes is useful for this study. We used the root-mean-square (rms) amplitude, coherent energy, Sobel-filter similarity, and spectral decomposition attributes. The rms amplitude is an amplitude-based attribute that computes the square root of the sum of the squared amplitudes divided by the number of samples in a window. The rms amplitude attribute has been historically used as one of the direct hydrocarbon indicators (DHI) in several areas; however, this attribute also gets affected by several other factors, including laminated sand-shale sequences, stratigraphic continuity, and chaotic slumping. The other factors may result in significant amplitude anomalies, which can be mistaken as hydrocarbon accumulation. However, this attribute can be used as a first step to identify the presence or absence of any anomalies. The results from this attribute depend on the noise in the data. We used 9 ms as the window length for the rms attribute in this study. Next, we computed the coherent energy volume. Coherent energy estimation involves singular-value decomposition of the seismic data in different principal components, with different windows aligned along the structural dip. In our coherent energy estimation, we take the energy estimate of the first principal component (Verma et al., 2016). Although it is similar to rms amplitude, it is less susceptible to noise because the first few principal components of the seismic data contain the signal, whereas the latter components contain noise (Chopra and Marfurt, 2007). Coherent energy permits estimates of effects due to the subtle differences among lithology, fluid, porosity, and thickness (Johnston, 2010; Sarkar et al., 2016). Next, we computed Sobel-filter similarity, which is a geometric attribute for image enhancement. It is an edge-detection attribute used to identify channels, faults, and other types of breaks in reflectors. It computes the square root of the sum of the squares of derivatives of the amplitude in the inline and crossline directions along the true dip. The classical Sobel filter, S, on a flat image can be expressed as
where x and y represent different axes and a represents the amplitude. In the case of 3D seismic data, we normalize equation 1 by the rms amplitude within the same seismic window. The normalization process enhances lower amplitude edges when applied to a coherence input volume (Chopra and Marfurt, 2007). If two neighboring traces are similar or contiguous to each other in the zone of interest, then the Sobel-filter similarity value is one; if they are completely different, then the value of the similarity is zero. Please see Luo et al. (1996) and Marfurt (2006) for detailed mathematics of Sobel-filter similarity.

We performed spectral decomposition to delineate channels of different thicknesses. It is a useful technique to decompose the seismic data into different frequency elements that correspond to thickness (Peyton et al., 1998; Partyka et al., 1999). As channels become very thin (below ¼ wavelength), their waveform shape becomes constant (the amplitude still changes), such that similarity or coherence attributes based on waveform similarity cannot detect the channels anymore (Liu, 2006; Marfurt, 2006). The spectral decomposition images provide complementary information to the similarity attributes because they are more sensitive to channel thickness rather than the lateral changes in seismic waveform or amplitude (Liu, 2006; Lyu et al., 2019). Spectral decomposition-based maps are typically interpreted qualitatively, using seismic geomorphologic pattern recognition or semiquantitatively, to infer the relative thickness variation of beds with acoustic impedance contrast. Spectral decomposition can be used in reservoir characterization, estimating temporal bed thickness, identifying stratigraphic traps, and facies classification (Ye et al., 2019). Because subtle stratigraphic traps mostly hold hydrocarbons in these reservoirs, we think the application of the abovementioned seismic attributes is appropriate.

After interpretation of the geologic features of interest using 3D seismic data and seismic attributes at a high level, we used processed conventional well logs for the identification of lithology and calculation of basic petrophysical and reservoir properties. We used the Simandoux equation to estimate water saturation for the Nanushuk and Torok Formations, rather than the well-known Archie equation. The use of the Simandoux equation is justified compared to simple Archie’s equation in these laminated sand-shale reservoirs. The Simandoux equation corrects the effect of shale in the zone of study for fluid saturation calculation, whereas Archie’s equation does not take into account the effect of shale. The Simandoux equation can be expressed as follows:
where Sw is the water saturation, a is the tortuosity, Rw is the formation water resistivity, ϕ is porosity, m is the cementation factor, Vsh is the shale volume, Rsh is the shale resistivity, and Rt is the formation resistivity. We also used advanced petrophysical logs, such as PNS, NMR, and sonic scanner tools in conjunction with conventional well logs to identify hydrocarbon-rich sandstone reservoir zones. In places where the advanced petrophysical logs such as PNS were not available, we performed joint inversion of conventional well logs to generate a statistical multimineral solution to identify hydrocarbon-rich sandstone layers. This type of petrophysical modeling was useful to better understand the heterogeneities of these reservoirs. In essence, the statistical multimineral models are generated using a linear inversion technique (Moses and Harrison, 1985; Mitchell and Nelson, 1988; Kulyapin and Sokolova, 2014; Bhattacharya and Carr, 2019). The input to the inversion process is the conventional well logs, such as the gamma-ray, resistivity, bulk density, neutron porosity, and photoelectric logs, and the output mineral solution is composed of varying proportions of quartz, clay, carbonate, pyrite, oil, and water. Such a solution can be validated against available PNS logs and/or core data. After petrophysical modeling, we performed rock physics analysis to understand the relation between the acoustic impedance, VP/VS ratio, hydrocarbon saturation, and lithology. It is known that the combination of the VP/VS ratio and acoustic impedance can be used to detect hydrocarbons in rocks. The presence of hydrocarbons decreases the VP/VS ratio. The analysis was based on petrophysical properties derived from conventional and advanced well logs (such as the sonic scanner and PNS tools). Such an analysis will be significantly helpful to predict the distribution of the hydrocarbon-charged (HC) sand bodies in the reservoirs using 3D seismic data.

3D seismic reflection data and attribute-assisted geologic feature detection

We found low-angle clinothems to be distributed over a large area (Figures 4, 5, and 6). The clinoforms mapped in the Nanuq South 3D survey area illustrate the distribution of the prograding shelf edges mostly along the north/northeast–south/southwest, which change its orientation toward the northwest near the NE-NPRA 3D survey. In the NE-NPRA 3D and Puviaq South 3D survey areas, the clinothems are larger and mostly oriented along the northwest–southeast. This pattern is commensurate with the shelf-trajectory maps by Houseknecht (2019a, 2019b) based on 2D seismic lines. Mapping of the clinoforms in all the 3D seismic surveys reveals that the shelf edges continue over tens of miles. Traditional seismic stratigraphic interpretations are out of the scope of this study. Please see Houseknecht (2019a) and Frierson (2013) for details on seismic stratigraphy in the NPRA. In this study, we rather focus on the quantitative aspects using basic and advanced seismic attributes. The seismic attributes reveal the presence of these features: channels, shelf edges, beach ridges, canyons, slump, and basin-floor fans.

The rms attribute extracted on the seismic section shows amplitude anomalies, which could be related to the presence of sand bodies along the shelf edge, slope apron, basin floor, cemented horizon, and mass transport deposit. We found amplitude anomalies at various places in the clinothems, especially near the shelf edges and basin floor (Figures 4 and 5). The extraction of the coherent energy attribute shows the evidence of high-energy (dominantly sand bodies based on drilling results) deposits in two dominant forms: (1) laterally continuous features along the shelf edges and (2) sporadically present on the basin floor (Figures 7 and 8). This will have an implication on the exploration program in this area in the future. It is worth noting that the Horseshoe 1 well was drilled near one of the shelf edges in the Nanuq South 3D survey area (Figure 3), whereas the CT-1 and CT-2 wells were drilled on the basin-floor fans in the Torok Formation in the west of the NE-NPRA 3D survey area (GEOExPro, 2017). The coherent energy attribute also shows the evidence of possibly beach ridges (Figure 6). Beach ridges are wave-swept or wave-deposited ridges that run parallel to a shoreline. These are commonly composed of sand as well as sediment reworked from underlying beach material. Although they have less potential for preservation, beach ridges have been identified using advanced 3D seismic attributes in other areas, such as the Australian Northwest Shelf (Bourget et al., 2014; Paumard et al., 2017, 2019) and North Sea (Jackson et al., 2010). To the best of our knowledge about the published studies on these Brookian Formations on the North Slope, beach ridges have not been described in the subsurface, which could be due to their subtle seismic expression or poor seismic imaging at places. Another possible explanation of the presence of such features at this particular location can be attributed to the shoreline trajectory and progradation. Shoreface and delta-front sand bodies can be observed on the coherent energy attribute (Figure 8). A core-based interpretation by LePain et al. (2018) confirms their presence in the NPRA.

We corendered two seismic attributes (coherent energy and Sobel-filter similarity) to further the identification of the geologic features present in these formations. In the NE-NPRA 3D survey area, it shows the evidence of distributary channels feeding sand (indicated by the coherent energy values) to the shelf edges and canyons on the slope (Figure 9). We interpreted the presence of large basin-floor fan-type deposits based on seismic attributes. In a seismic reflection section, these deposits are expressed as positive-relief and undeformed mounds on the basin floor. High coherent energy anomalies in the fan deposits could be indicative of shelf-winnowed sands or mass transport deposits or a hard shale. The dimension of the basin-floor fans varies in the study area. Apart from the basin-floor fans, we observed parallel-to-near-parallel elevated ridges on the basin floor in the Torok Formation (Figure 9). This pattern is indicative of slump deposits, formed due to the continued migration of the shelf trajectories and oversteepening, which might have resulted in the collapse of the unstable upper slope material. At places, the foresets (in the vertical seismic section) show contorted geometries and concave-up features near the slump deposits, which we interpreted as slump scars. Their presence influences the depositional pattern of the sediment and can restrict the flow of sand to the deeper portion of the basin. These submarine slides affect the lower portion of the Torok Formation. Most of these slides overlie and downlap onto either the gamma-ray zone (GRZ) or Lower Cretaceous unconformity. An offset 3D seismic survey in the Harrison Bay area (southeast of the NE-NPRA 3D survey) confirms the presence of giant slides in this area, some of which are related to the well-known Fish Creek Slide (Homza, 2004). The seismic sections from the NE-NPRA 3D and Harrison Bay 3D surveys show several slide features (Figure 10 and 10). The slide blocks are expressed on the seismic sections by the mound-like steep ramp/escarpment characteristics, associated with axial lows that meet the detachment surface. It is interesting to note that the presence of the GRZ can generate a velocity pushdown or traveltime delay effect in places (e.g., in the Harrison Bay 3D survey area) because of its low velocity, which was also noted by Weimer (1987) and Homza (2004). The slide block near the northern portion of the NE-NPRA 3D survey could be related to the Dalton Slide, whereas the ones present in the Harrison Bay 3D survey are mostly related to the Fish Creek Slide (Homza, 2004; T. Homza, personal communication, 2019). Some of the slide blocks are narrow and steeply dipping with significant amplitude contrasts. The geometry of the slides in the subsurface of the North Slope of Alaska is complex, which can be distributed in an organized and disorganized fashion, and sometimes these are just erosional remnants. In the Harrison Bay 3D survey, there is a large erosional feature, which appears to have affected the Fish Creek Slide blocks (possibly younger based on seismic reflector terminations) (Figure 10). The seismic expression of the erosional feature is similar to a broad valley, bound by steep escarpment surfaces against the Fish Creek Slide blocks. One can alternatively interpret the parallel- to near-parallel, narrow, laterally extensive, concave-up features (which could be indicative of erosion by bottom currents) on the stratal slice in the NE-NPRA 3D survey, with high coherent energy anomalies as contourite sand deposits. Contourites are common in the deepwater environment. According to Brackenridge (2014), there could be at least five types of sand-rich contourite deposits, with variable reservoir quality. Weimer (1991) suggests the lower likelihood of contourite deposits in this area because it would require the lower slope material to be eroded significantly prior to the deposition of the contourites. We interpreted several channels on the slope in the NE-NPRA 3D survey, most of which have low to moderate sinuosity (Figure 9).

We also observed the presence of large, moderately sinuous channels in the Puviaq South 3D survey in the western part of the NE-NPRA 3D survey (Figure 11). Most of these channels are oriented along northeast–southwest. It is interesting to note that the amplitude anomalies along the channels are not uniform, which could be indicative of the varying sand/clay ratio along the channel (Figure 11). Most of these channels are along the slope (toward the northeast in the Puviaq South 3D survey), which can explain their varying amplitudes due to marine influence. We performed spectral decomposition on the channels to reveal the relative thickness variation (Figure 11). We found that not all of these channels have a uniform thickness. High-frequency components (i.e., the blue color, 50 Hz) correspond to thin channels, and low-frequency components (i.e., the red color, 20 Hz) indicate thick channels. The dominant frequency is approximately 38 Hz, which is close to the central frequency used in the spectral decomposition process. In the Puviaq South 3D survey area, we observed thick channels near the western side, whereas thin channels are near the northeastern direction.

It is well-known that several factors, such as laminated depositional pattern, stratigraphic continuity, chaotic slumping, tuning thickness, and fluid, can impact the amplitude anomalies. Sometimes, rms amplitude can serve as a DHI; however, the use of amplitude anomalies as a DHI on the North Slope can be hard because some dry holes look significantly good on poststack seismic data. In general, we can observe the effects of fluid, especially if oil-bearing sand, on far-stack seismic or amplitude variation with offset (AVO) gradient volumes. Figure 12 shows the comparison of amplitude anomalies among the near-, mid-, and far-offset stack seismic data around the Horseshoe 1 discovery well. We found a subtle difference between the near- and far-offset seismic signatures in terms of amplitude anomalies. We observed a blob of a small amplitude anomaly just below the N2 reflector on the far-offset data, corresponding to low gamma-ray readings. This is the zone of hydrocarbon discovery in this well. In this case, the far-offset stack data have some issues related to frequency absorption. Future work will involve a detailed analysis of angle gathers, AVO, and seismic inversion to better decompose the impact of lithology and fluid from the seismic data. Because of these uncertainties, we performed a detailed petrophysical characterization to analyze the complexities of these formations to detect hydrocarbons at the well scale.

Core and well-log-based petrophysical characterization

Well-log correlation shows multiple coarsening (or cleaning) upward sequences punctuated by a sharp increase in gamma-ray response in the Nanushuk Formation, which is indicative of progradation and flooding surfaces at the top of the parasequence set, respectively (Figure 13). Interpretation of the conventional well logs (e.g., density porosity) calibrated to core data in the Horseshoe 1 well shows the presence of sandstone, with porosity in the range of approximately 9%–21%. These values fall in the general range of core-measured porosity values reported in these sandstones from other areas on the North Slope; see Houseknecht (2019a) for further details. It is worth noting that the resistivity responses in these laminated sand-shale reservoirs are mostly low in the study area. In general, a high-resistivity response is indicative of hydrocarbon-bearing clean sand; however, the resistivity reading (approximately 810  Ωm) is suppressed in these formations in several areas on the North Slope. This is due to the interbedded nature of sand-shale sequences in the formation, which can significantly affect the petrophysical log responses and reservoir estimates. These types of low-resistivity reservoirs have been found in other parts of the world, including the Gulf of Mexico, North Sea, West Africa, and Indonesia (Darling and Sneider, 1992; Worthington, 2000; Audinno et al., 2016). The presence of conductive clay minerals affects the resistivity response (Worthington, 2000; Hamada and Al-Awad, 2000). In addition, the presence of a small amount of pyrite and siderite can reduce the resistivity of the rock (Pratama et al., 2017). Also, the dispersed clay can affect the oil saturation and fluid mobility. However, there are exceptions to this observation, such as the Umiat area near the Brooks Range, where the Nanushuk reservoir exhibits high resistivity in places, which may be due to the increase in quartzose components and the presence of permafrost. The water saturation based on the Simandoux equation matches well with core-based water saturation data from the Horseshoe 1 well (Figure 14). The main parameters we used in the Simandoux equation for this well are as follows: water resistivity, 0.24  Ωm; tortuosity, 1; cementation factor, 1.4; saturation exponent, 2; and shale resistivity, 3 Ωm. These values are approximate and vary across the flow units. We also boosted cementation factor and water saturation when the shale volume and porosity values cross certain thresholds.

Because conventional well logs have certain limitations (e.g., low vertical resolution and a combination of mineral and fluid effects, etc.), we used advanced petrophysical logs (e.g., PNS and NMR) and modeling techniques (e.g., joint petrophysical inversion) to analyze the formation heterogeneities and identify the oil-saturated sandstone reservoirs (Figure 14). Joint inversion of conventional logs produces a multimineral solution, which is composed of varying proportions of quartz, clay (i.e., illite and chlorite), carbonate (i.e., calcite and dolomite), oil, and water. The interpretive petrophysical inversion model shows the vertical heterogeneity of these formations and the relative increase of quartz content near the top of the parasequences (Figure 14). The quartz content increases from 40% to 80% from the base to the top of some of the productive parasequences in the Horseshoe 1 well. Figure 14 shows an increasing oil saturation near the top of the parasequence. At least 15%–20% of the pore volume appears to be saturated with oil near the top of the parasequence, which goes down to 5% near its base. The water content (including clay bound-water) increases toward the base of the parasequence. We also had access to the downhole fluid sampling data (modular formation dynamics tester or MDT) from the Horseshoe 1 well. Two MDT test results at 4760.69 and 4770.95 ft show the presence of oil in the interval, which matches the petrophysical-inversion-based multimineral solution. This is a final validation test of the petrophysical inversion model. In addition to the petrophysical model, we used NMR-log-based permeability estimates in reservoir heterogeneity evaluation. NMR log measurements show the gradational increase in the permeability from the bottom to the top of the parasequences. Based on these observations, we can divide most of the parasequences into at least two to three petrophysical zones (zone 1, 2, and 3) with varying mineralogy, porosity, permeability, and hydrocarbon saturation. In the NMR log, we observed a bimodal distribution of the T2 time (the proton precession relaxation time), which is indicative of oil in these clastic rocks (Figure 15). This can be validated with the core photographs showing the oil-stained Nanushuk Formation. Figures 14 and 15 show the reservoir heterogeneity. Based on our analysis and Houseknecht (2019a), permeability in the Nanushuk sandstones ranges between approximately 0.001 and 1000 mD at places; however, there are a few areas in the NPRA where these formations have reported higher permeability than that. The Torok Formation is thicker than the Nanushuk Formation. The Torok Formation has a higher proportion of shale than the Nanushuk Formation; however, there are thick sandstone packages present in this formation. The Cronus 1 well, near the Nanuq South 3D survey area, has a few sandstone packages (probably turbidites) near the bottom of the Torok Formation (Figure 16). The petrophysical inversion model shows that the general mineralogy of the Torok Formation is similar to the Nanushuk Formation; however, the former contains illite and chlorite type of clay minerals. The presence of chlorite can affect reservoirs in different ways. In deeply buried clastic reservoirs, chlorite can form coatings on the quartz grains, which impedes quartz-overgrowth and it can lead to good reservoir quality (Wilson et al., 2014; Worden et al., 2018). On the other hand, if the amount of chlorite is high, it can block pore throats and reduce permeability significantly (Worden et al., 2018). The availability of point-count data from the Cronus 1 well and X-ray diffraction results from other wells indicate the presence of different types of clay minerals in varying proportions. The Torok Formation also contains high-density minerals, including dolomite fragments at places. These sandstones have porosity varying between approximately 6% and 14%, with permeability ranging between approximately 0.001 and 4 mD in this well. The porosity and permeability values in the Torok Formation in the Cronus 1 well fall in the general range of core-measured values reported from other wells on the North Slope (Houseknecht, 2019a). Similar to the Horseshoe 1 well, the availability of an NMR log in the Cronus 1 well revealed a typical bimodal T2 amplitude distribution (across the sandstones shown in Figure 16), which is indicative of hydrocarbons. These sandstone packages could be good reservoir targets, with hydraulic stimulation. It is worth noting that the mineralogy and reservoir properties (porosity and permeability) in these formations change spatially to an extent due to several reasons (such as provenance, clay type, thermal maturity, burial depth, and pore fluid chemistry); however, the multiwell petrophysical analysis is out of the scope of this study.

Rock physics analysis

In the last step, we performed rock physics analysis on well-log data obtained from sonic scanner tools (e.g., dipole sonic) in a quantitative way to detect the hydrocarbon-bearing sand zones. These well-log-based relations can later be used in the seismic inversion process to predict the sandstone bodies and potential reservoirs. We plotted the VP/VS ratio against the acoustic impedance of P-wave (Figure 17). This plot indicates that if the VP/VS ratio and the acoustic impedance of the sandstone units are significantly low (less than 1.7, or at least less than 2.0 due to the significant overlaps of the data points), there is a chance of finding hydrocarbons in certain areas. One should keep in mind that this proxy based on a low VP/VS ratio and low acoustic impedance is a combination of several factors, including lithology, reservoir quality, porosity, fluid saturation, and burial depth. The VP/VS ratio in these formations changes across the NPRA. In the Pikka-Horseshoe discovery area, the VP/VS ratio can be high; however, the acoustic impedance still remains low (less than approximately 30,000 g.ft/[cm3.s]) in the hydrocarbon-bearing sandstones. In the Umiat area, near the Brooks Range, the VP/VS ratio in these sandstones is low. Based on the separations among the data points color-coded by the gamma-ray curve, it also appears that this plot is a good lithology indicator. Figure 17 shows a similar plot in the wells tested to be dry (such as the Grizzly 1 well) in the studied interval. The HC-sand cluster is completely missing from the dry hole. We can use such rock physics relations to highlight and detect hydrocarbon-charged sand bodies in the reservoir using 3D seismic and petrophysical logs. We recommend taking significant precautions before using this proxy or just the deterministic cutoff values to detect and predict hydrocarbons in this area.

This work focuses on 3D seismic geomorphology and petrophysical characterization to identify lithology and potentially fluid in the Nanushuk and Torok Formations. As already stated, more advanced work is needed to better characterize these formations, specifically, detection and prediction of hydrocarbon reservoir units, which cannot be done on poststack seismic data without uncertainties. Our future work will include data reprocessing, multiwell petrophysical studies, AVO analysis of angle (offset) gather data, detailed rock physics, and prestack seismic inversion to derive acoustic impedance, VP/VS, and lambda-rho and mu-rho parameters from seismic data to predict the sweet spots.

The analysis of the publicly available 3D seismic surveys in conjunction with petrophysical and rock physics analyses provided a detailed understanding of the complex Nanushuk and Torok Formations. We summarize the key conclusions from this study here.

  1. 1)

    Seismic-attribute-assisted mapping of the Nanushuk and Torok Formations reveals the presence of channels, prograding shelf edges, canyons, slumps, slide blocks, and basin-floor fans, many of them with significant amplitude anomalies. The shelf edges continue for tens of miles. We found advanced seismic attributes, such as coherent energy, Sobel-filter similarity, and spectral decomposition, useful in identifying these geologic features and quantifying their extent.

  2. 2)

    The internal character of the Nanushuk and Torok Formations delineated by well logs and petrophysical inversion modeling shows their vertical heterogeneities, which affect the reservoir properties, such as porosity and permeability. These are laminated, low-resistivity reservoirs. The distribution pattern of sand and shale layers (interbedded or massive) in these formations can affect the petrophysical log responses, fluid mobility, and reservoir estimates. This also depends on the clay mineral types. It is important to keep in mind that the heterolithic nature of these formations is not always accurately represented on the conventional log responses with their limited resolution; hence, the character of these formations can be oversimplified. Therefore, detailed core analysis and advanced high-resolution petrophysical logging are critical.

  3. 3)

    The results from the core data and calibrated petrophysical inversion model, in conjunction with the high-resolution NMR logs, show that only a few zones in the same parasequence are oil-saturated. We think that this will have implications on reservoir estimates, targeting, deliverability, and field development in the future. Because of the low effective porosity and permeability in these formations, hydraulic stimulation could be a good strategy to produce hydrocarbons.

  4. 4)

    The combination of low VP/VS ratio and low acoustic impedance is a good indicator of the hydrocarbon-bearing sand in these formations.

We thank the Alaska Department of Natural Resources, Division of Oil and Gas, for making the tax-credit 3D seismic data available to the University of Alaska Anchorage and the University of Texas Permian Basin (State of Alaska, 2017). We also thank Schlumberger and CGG for donating the Petrel and PowerLog software to the University of Alaska Anchorage and the University of Texas Permian Basin. We computed seismic attributes using AASPI Consortium software. We appreciate the three anonymous reviewers whose comments helped to strengthen the manuscript.

Interested persons can request the 3D seismic data from the Alaska Geologic Materials Center (http://dggs.alaska.gov/gmc/seismic-well-data.php) and well log data from the Alaska Oil and Gas Conservation Commission.

Biographies and photographs of the authors are not available.

Freely available online through the SEG open-access option.