Enhancing Paleoreef Reservoir Characterization through Machine Learning and Multi-Attribute Seismic Analysis: Silurian Reef Examples from the Michigan Basin

: Historically, Silurian pinnacle reef complexes in the Michigan Basin have been largely identiﬁed using 2D seismic with very little research on the reservoir characterization of these reefs using 3D seismic data. By incorporating a high-resolution 3D dataset constrained by a well-studied and data-rich paleoreef reservoir, the Puttygut reef, seismic attributes were correlated to petrophysical properties through machine learning and self-organizing maps (SOMs). A suite of structural and frequency-based attributes was calculated from pre-stack time migrated (PSTM) seismic data, with only a subset of them selected as SOM inputs. Structural attributes enhanced details in the reef but frequency attributes were overall more useful for correlating with reservoir quality. A strong relationship between certain combination percentages of attributes and certain sections of the reef with porosity and permeability was found after the SOM results were compared to wireline log and core analysis data. Areas with high permeability and porosity correlated with the average frequency and spectral decomposition at 29 and 81 Hz. Areas with high porosity and varying permeability correlated with the average frequency and spectral decomposition at 29, 57, and 81 Hz. Areas with intermediate porosity correlated with the average frequency and spectral decomposition at 29 and 57 Hz. The efﬁcacy of the procedure was then demonstrated on two nearby reefs with very similar results.


Introduction
The focus of this study was to characterize features in ancient reef reservoirs that are at and below conventional vertical seismic resolution by creating a machine learning and multi-attribute analysis workflow. A seismic attribute, first defined by Taner and Sheriff 1977 [1], is a measurement derived from seismic data, most commonly based on time, amplitude, frequency, and/or attenuation. The primary use for attributes is that they can help an interpreter visualize features, relationships or patterns that might go unnoticed [2,3]. An ideal area to develop and test such a workflow was in the Michigan Basin on Silurian Niagara-Lower Salina pinnacle reefs. The reefs form a sublinear trend on the northern and southern rims of the basin at depths between 900 and 2200 m [4]. The basin encompasses the entire lower peninsula and parts of the upper peninsula of Michigan, as well as adjacent areas of eastern Wisconsin, northeastern Illinois, northern Indiana, northwestern Ohio, and western Ontario [5]. As of 2019, these reefs have produced over 500 million barrels of oil and over 3 trillion cubic feet of natural gas [6]. In addition to what has already been produced, there is still great potential for secondary oil recovery through water, CO 2 or gas injection. Michigan reef reservoirs have 26% average primary recovery and 12.5% average secondary recovery with only around 5% of the discovered drilling future gas storage wells. Overall, if you have seismic data, a way to calculate attributes, and a method of machine learning, you will be able to explore that seismic data with increased detail and clarity which was not possible before with just the data and attributes. From this study, a relationship was found between the SOM results and the porosity and permeability within the Puttygut reef. Furthermore, applying this workflow on the Ira and Lenox reefs confirmed the repeatability of this method by finding very similar results.   (Brown) Formation. An A-A' line is drawn through the four wells used in this study. There are also all wells in and near the reef marked by the white and black dots. The black dots represent the wells that were cored. Image (D) is a generalized stratigraphic scheme of the reefs in the southern trend of the Michigan Basin from [8]. In this region, the Clinton formation is refered to as the Cordell formation.
It is important to increase the level of detail seen in seismic data to advance our understanding of the subsurface. With machine learning and multi-attribute analysis, it is possible to re-explore pre-existing seismic data through a new lens and see features on a much finer scale than previously possible. Machine learning can open new doors to exploring our geologic past and to increase our understanding of smaller scale events and processes. It can aid us in identifying and characterizing features of interest that were previously hidden in the data, and it is possible to use this method on any geologic regime.
More practically, the machine learning method can be used to increase the level of detail seen in reservoirs which makes the task of characterization much easier. With the increased level of detail within the reservoir from the machine learning workflow, utilizing the reef for gas storage or carbon sequestration will be more precise, accurate, and easier overall. The workflow from this research was developed in an attempt to increase confidence when interpreting seismic data and the lower risk associated with planning and drilling future gas storage wells. Overall, if you have seismic data, a way to calculate attributes, and a method of machine learning, you will be able to explore that seismic data with increased detail and clarity which was not possible before with just the data and attributes. From this study, a relationship was found between the SOM results and the porosity and permeability within the Puttygut reef. Furthermore, applying this workflow on the Ira and Lenox reefs confirmed the repeatability of this method by finding very similar results.

Regional Geology
The Michigan Basin is a circular, intracratonic basin that covers an area of~316,000 km 2 shown in Figure 1 [5]. The reefs in this basin are Wenlock (Homerian) in age [8] and form a sublinear trend on the northern and southern rims at depth between 900 and 2200 m [4]. The Niagaran Guelph Formation reefs are unconformably overlain by, and laterally encased in, the cyclic limestones and evaporites of the Silurian Salina Group, and underlain by the Lockport Dolomite Formation [12]. See Figure 1D for a generalized stratigraphic scheme from Rine 2019 [8].
A summary of the time-step basin wide depositional evolution of the Michigan Basin developed by Rine 2019 [8] can be generalized as follows: Through a time of sea level oscillations the first formation to be deposited in the basin is the Lockport formation. Eventually, a sea level fall exposes the Lockport and ends its deposition. Following is a sea level rise which coincides with the initiation of the bioherm mounds and marks the beginning of the lower Guelph Formation. Next, the sea level falls again which exposes the bioherm mounds and causes significant wasting. That is what creates the flat tops seen today. A rapid sea level rise follows initiating the pinnacle reef growth in the upper Guelph Formation. Tidal flat facies begin to be deposited as the basin experiences restriction. The relative sea level fall in the basin results in the A-1 Evaporite deposition. The sea level rises, falls, and rises again resulting in the deposition of the lower A-1 Carbonate, the rabbit ear anhydrite (REA), and the upper A-1 Carbonate, respectively. Moreover, at this point in time, the reefs are completely covered in sediment. The A-2 Evaporite section is then deposited with another fall in sea level. Finally, through another rise and fall in sea level, the A-2 Carbonate and the B-Salt are deposited, respectively, which marks the end of the relevant sequence for this study. While the primary depositional model for these reefs is relatively well known, the postdepositional diagenetic history of these reefs is still relatively unknown. Since the reefs in SE Michigan have been 100% dolomitized and are partially salt plugged, the primary depositional facies are not the only factors that contribute to the petrophysical properties, hence the need for utilizing seismic attribute analysis for better understanding reservoir quality distribution spatially.
The reservoir interval for these reef reservoirs extends from the top of the A-1 Carbonate (Ruff Formation) and continues through the Brown Niagaran (Guelph Formation) down to the base of the reservoir or Gray Niagaran (Lockport Dolomite Formation). While these formation's thicknesses vary widely throughout the reef, at well P-201 the thickness of the A-1 Carbonate is~44 feet, the tidal flat cap is~43 feet thick, the reef core is 115 feet thick (with a reservoir interval which is 105 feet thick), and the bioherm is 151 feet thick (with a reservoir interval which is 38 feet thick). See Figure 1 for locations of well P-201 and the formations. Since these formations are all comprised of dolomite, they do not exhibit good contrasts in acoustic impedance and are near or below seismic resolution. This general homogeneous lithology affected acoustic isotropy and contributes to the difficulty in determining where to pick these horizons by inspection of the seismic data alone. Figure 4 Geosciences 2021, 11, 142 5 of 23 shows an inline (east-west direction) through the Puttygut reef. The seismic horizons were selected by an expert geologist (M. Rine) in the area with substantial apriori geologic and structural knowledge of the Michigan Basin in order to achieve the most accurate horizons picks by correlating the seismic with the available well data and cross sections. The A-2 Carbonate and Clinton Formation are easily visible owing to their acoustic impedance contrast to their surroundings. However, the A-2 Anhydrite, A-1 Carbonate, Brown, and Gray Formations are not so easy to resolve. These horizons are detectable at points along the line but are not confidently chosen without the prerequisite geologic and structural knowledge in the area. Unfortunately, the horizons which are more difficult to pick are the horizons that are of most interest in reservoir quality. The reef core and bioherm between the A-1 Carbonate and the Gray Formation are commonly the areas where the best reservoirs are found, as these areas of high permeability and porosity tend to exist in the skeletal wackestone and reef boundstone of the reef core, as well as the foreslope reef talus conglomerates that flank the windward (eastern) side of the pinnacles. The thrombolitic bindstone facies within the A-1 Carbonate that caps the reef complex in the reef crest position is also an area of high permeability and porosity [8].

Materials and Methods
The initial workflow consisted of running a large suite of attributes on the pre-stack 3D volume around the Puttygut reef. Geologically, we are primarily interested in improving the mapping of continuous areas of high porosity and permeability found in the subseismic formations that comprise the reef core and bioherm. As these areas are seismically difficult to detect, areas of interest with high permeability and porosity are initially chosen using core and well log data. There is one limitation with this approach, however, as the core and well log data are in depth and the seismic data are in time. There are no reliable sonic or density logs within the reef proper, only sonic logs in nearby off-reef wells. Therefore, producing a reliable well tie within the reef is particularly challenging without making assumptions that drastically affect the reliability of the results. The workflow for picking sub-seismic horizons was as follows: (1) Pick formation tops from geophysical logs and cores; (2) tie the wells to the seismic using a VSP from the P-201 well (see Figure 1C for well location) and synthetic seismograms generated from wells with DT and RHOB logs outside of the survey; (3) shift the wells to time using key horizons (i.e., A-2 Carbonate, A-2 Salt, A-1 Carbonate, Clinton. See Figures 2 and 3 for the VSP results. See Supplementary Figures S1-S3 for uncropped versions); (4) pick sub-seismic horizons that cross-cut reflectors based on well control (i.e., Brown Niagaran, Bioherm, Gray Niagaran). From there, the smaller, sub-seismic horizons were mapped through apriori structural geologic knowledge and using the pseudo-tied wells to aid in picking locations. While this workflow is not ideal, there is a very high level of confidence with the accuracy of these provided horizons.
Attributes were generated and evaluated by the interpreter to determine how effective they were at displaying minute details within the target reef. Attributes that exceled at revealing the internal, less detectable aspects, of the reef were then selected for multi-attribute analysis in the form of self-organizing maps (SOMs) to achieve sub-seismic resolution.
tors". These output detectors are then fit to the input characteristics such that they both attempt to classify the data and are influenced by the other detectors. What is made in the end is a map which is organized in a meaningful manner that can be used for potential pattern recognition. These inputs can be anything from variations of textures in an image, to the variations in the quality of life between nations (an example created by Kohonen 1995 [19]), to the different attributes that can be created using seismic data. As mentioned previously, these SOMs analyze the input data at the sample interval level.    In machine learning, for a desired anthropologically reasoned outcome, rather than choose characteristics randomly, it is optimal to take advantage of the interpreter's knowledge when choosing SOM attributes and parameters. In this case, two categories were decided upon when initially analyzing single attributes: Structural based attributes and frequency-based attributes. The structural based attributes were subjectively graded on how well they illuminated the internal structure of the reefs. The frequency-based attributes were similarly graded on how well they distinguished unique sections of the reef. Highlighting the sub-seismic internal structure of the reef and identifying different sections of the reef via frequency content potentially create the key to isolating areas of higher porosity and permeability. Narrowing down the attributes to a smaller set is also useful in the SOM stages later in the workflow. It is not necessarily beneficial to have more attributes when conducting multi-attribute analysis as many attributes are calculated in similar ways and can show similar features that can be unnecessary [3]. If all the attributes that were calculated are used in one SOM then there would be a great amount of unnecessary redundancy. In addition, using more data in the calculations will increase the The theory of using SOMs and multi-attribute analysis to image and interpret subseismic resolution has been successfully tested and proven in the past [11]. Most commonly SOMs were calculated over unconventional plays, such as in the Denver-Julesburg Basin and in the Eagle Ford shale [13][14][15]. They have also been tested in other areas, such as for interpreting direct hydrocarbon indicator (DHI) characteristics, classifying carbonate facies, and recognizing geologic patterns with varying degrees of success [16][17][18]. A SOM was explored in detail by Teuvo Kohonen [19] in his book where a very mathematical description of SOMs can be found. Kohonen describes SOMs as "nonparametric regressions" where a number of ordered discrete reference vectors are fit to a distribution of vectorial input samples. In order to approximate continuous functions though, the reference vectors are used to define nodes or neurons, in a hypothetical elastic network. Local interactions between these neurons along the neural network create that essence of elasticity. By setting up the map in this way, the neurons develop into specific detectors of signals in the input space. These detectors are formed in the map in a meaningful order as if falling into place on a feature coordinate system. These resultant feature maps are used for the preprocessing of patterns for recognition [19]. Stated simply, SOMs require some number "m" of input characteristics and a desired number "n" of output "detectors". These output detectors are then fit to the input characteristics such that they both attempt to classify the data and are influenced by the other detectors. What is made in the end is a map which is organized in a meaningful manner that can be used for potential pattern recognition. These inputs can be anything from variations of textures in an image, to the variations in the quality of life between nations (an example created by Kohonen 1995 [19]), to the different attributes that can be created using seismic data. As mentioned previously, these SOMs analyze the input data at the sample interval level.
In machine learning, for a desired anthropologically reasoned outcome, rather than choose characteristics randomly, it is optimal to take advantage of the interpreter's knowledge when choosing SOM attributes and parameters. In this case, two categories were decided upon when initially analyzing single attributes: Structural based attributes and frequency-based attributes. The structural based attributes were subjectively graded on how well they illuminated the internal structure of the reefs. The frequency-based attributes were similarly graded on how well they distinguished unique sections of the reef. Highlighting the sub-seismic internal structure of the reef and identifying different sections of the reef via frequency content potentially create the key to isolating areas of higher porosity and permeability. Narrowing down the attributes to a smaller set is also useful in the SOM stages later in the workflow. It is not necessarily beneficial to have more attributes when conducting multi-attribute analysis as many attributes are calculated in similar ways and can show similar features that can be unnecessary [3]. If all the attributes that were calculated are used in one SOM then there would be a great amount of unnecessary redundancy. In addition, using more data in the calculations will increase the amount of time required to complete a SOM. With these criteria in mind, the chosen set of attributes are defined below: Structural attributes: The following structural defining attributes were chosen since they all excelled at defining some form of structure (internal or boundary) and since they are all calculated in reasonably different ways, which is important later during the multi-attribute analysis step.

•
Energy ratio similarity (ERS): An attribute which is designed for edge-detection. ERS is similar to other coherence attributes except that it is calculated along the structure and that it uses analytic traces rather than just the real traces. ERS is the ratio of the energy of the Karhunen-Loève filtered data over the total unfiltered energy of the input data within the analysis window [20]. For the purposes of this research, ERS discriminated best the reef margin as well as some internal features that other attributes were unable to detect. The idea for reef edge and internal structure detection by the use of coherence or similarity attributes has been tested extensively in the past. Chopra and Marfurt 2007a [20] and Skirius et al. 1999 [21] successfully used ERS to aid in carbonate reef interpretations. • Aberrancy: Defined as the third derivative of structure, it is the measure of the asymmetry of a curve about its normal [22]. Aberrancy measures the lateral changes in the curvature along a surface and is also able to detect faults that are below seismic resolution [23]. Qi and Marfurt 2018 [23] calculated aberrancy on a data volume over the Barnett shale gas reservoir in Fort Worth Basin, Texas and were able to identify small karst features unseen with coherency. Aberrancy was able to define smaller and more internal reef structures than other structural attributes for the purpose of this study.   [28] later used this attribute to help extract more information on channel interiors while using a coherence attribute to distinguish the edges of the channels. GLCM was useful for identifying additional features within the reef core that were unnoticed by the other structural attributes.
Frequency attributes: The following frequency attributes were chosen based on how well they differentiated the internal reef structure.

•
Cosine of Instantaneous phase: Instantaneous phase attribute provides the interpreter with the phase of the wavelet at a certain location in the data. Yet, instantaneous phase is cyclical in nature, wrapping around from −180 to +180 • , and is also mathematically discontinuous. The cosine of the phase is implemented in order to remove those discontinuities [29]. Taking the cosine of the instantaneous phase also removes the cyclical nature of the attribute, creating a scale from −1 to +1, which is advantageous later on when being used in a SOM. Sukmono et al., 2006 [30] and Huang et al., 2011 [31] incorporated the generation of cosine of instantaneous phase (or just instantaneous phase) in order to aid in carbonate reef interpretations. • Average frequency: Instantaneous frequency is a measure of the frequency of the wavelet of the seismic trace at a given location in the data [29]. Toelle and Ganshin 2018 [32] found a relationship between the instantaneous frequency of seismic data and the amount of porosity present in reefs located in the northern Michigan Basin. They found that the instantaneous frequency decreased from 73 to 45 Hz as the porosity increased from 5% to 20%. They interpreted this decrease in frequency to be associated with frequency attenuation from the increase in porosity. Huang et al., 2011 [31] and Sarhan 2017 [33] also used instantaneous frequency to aid in carbonate interpretations. Unfortunately, instantaneous frequency tends to have spikes and noise that can mask important information. One way to counter this interference is to take the average of that instantaneous frequency. Consequently, the reduction in noise and spikes in this attribute from averaging outweigh the potential loss of resolution [34]. The average frequency from the Puttygut volume produced pockets of low frequency near wells with higher porosities and sections of high frequency in other areas. • Spectral decomposition (continuous wavelet transform (CWT)): This is an attribute that separates the amplitude seismic data cube into different frequency cubes. Although there are more than one way to calculate spectral decomposition, the CWT is a linear form which is based on the short time Fourier transform (STFT). Nejad et al., 2009 [35] and Saadatinejad et al., 2012 [36] used spectral decomposition successfully to aid in reef interpretation and exploration.
One of the ultimate goals of this research is to attempt to map/identify areas of high permeability and porosity within the reef reservoirs using machine learning and multi-attribute analysis. Yet, these two characteristics of reservoirs are challenging to unequivocally map with absolute accuracy [37][38][39][40][41]. These past studies have one thing in common when attempting to map or predict permeability or porosity: The attributes chosen are mostly derived from or related to the frequency and/or phase of the seismic data. In order to determine if the frequency-based attribute approach at identifying permeability and porosity is effective for the study area of this research, all attributes will be looked at initially on the provided high confidence horizons that represent the reef lithofacies. These detailed horizons make it possible to accurately estimate the location of high permeability and porosity zones surrounding wells that penetrate the reef. From there, it should be possible to work out from each well into the reef to determine how much variability in permeability and porosity there is. Four value ranges were set to differentiate "great" from "poor" permeability and porosity. The ranges for porosity are: 0-3% is poor, 3-6% is intermediate, 6-10% is good, and greater than 10% is great. The ranges for permeability are: <1 mD is poor, 1-10 mD is intermediate, 10-50 mD is good, and greater than 50 mD is great.

Puttygut Data
The Puttygut reef reservoir has a recent vintage 3D seismic acquisition and processing and seven cored wells with petrophysical properties and facies descriptions but has a relatively older vintage wireline log dataset. These log suites are lacking sonic and density logs which are needed for tying the wells to seismic. From the Puttygut seismic survey, two processed volumes were generated in 2019 (one is pre-stack time migrated and the other is post-stack time migrated). The inline and crossline spacing is 440 feet, with shot and receiver spacing of 110 feet, a final bin size of 55 by 55 feet, and a sample rate of 1 ms. The data is zero-phase with American polarity, where the positive impedance contrast is recorded as a peak. There are 30 wells in the area drilled within and around the reef ( Figure 1C). All 30 wells had at least a gamma ray log, 14 wells had neutron porosity logs, two had resistivity logs, one well had an SNP log, and one well had a DT log. There were also 11 wells that had core data that were compiled and transferred into the LAS format. Eleven wells had permeability and porosity data, nine wells had oil and water saturation data, and four wells had horizontal permeability measured at 90 • to the maximum permeability direction. In addition to all this well and 3D seismic data, there is a vertical seismic profile (VSP) at well P-201.
Generating well ties in the Puttygut seismic area was an iterative process. Owing to the lack of geophysical logs (sonic and density) within the Puttygut reef complex proper, fortunately a vertical seismic profile (VSP) was conducted in conjunction with the Puttygut 3D seismic acquisition in the P-201 wellbore. Since the reservoir interval (A-1 Carbonate to Gray Niagaran) was perforated throughout, a packer was set just above the A-1 Carbonate resulting in the A-2 Anhydrite being the deepest formation recorded during the VSP survey. That VSP survey (see Figures 2 and 3 for the VSP results. See Supplementary Figures S1-S3 for uncropped versions) was used to generate time-depth tables for the P-201 and those time-depths tables were applied to the surrounding reef crest wells. In the off-reef position, the Buszek well (PN 23843) had a sonic log that was used to generate a synthetic seismogram. The synthetic seismogram tied well to the 3D seismic (Figure 4), and thus was used to generate a time-depth table for that well and was applied to the surrounding off-reef wells. Since the VSP was unable to be conducted within the reef proper, the top (A-1 Carbonate/Brown Niagaran) and base (Gray Niagaran) horizons of the reef reservoir had to be interpreted by comparing the seismic to structural cross sections. Furthermore, the lack of acoustic impedance on the reef crest between the A-2 Carbonate (dolomite), A-2 Anhydrite (anhydrite), A-1 Carbonate (dolomite), Brown Niagaran (dolomite), and Gray Niagaran (dolomite), results in a lack of strong reflectors to correlate horizons from off-reef to on-reef. Therefore, while those horizons were picked by closely correlating with structural cross-sections from the available well control, there is some uncertainty associated with them.

Results
A wide range of single attributes were calculated and compared by the interpreter to decide which attributes best showed features of interest in the reef. Therefore, after grading according to the criteria previously described, the following attributes were deemed most worthy. Four structural attributes (energy ratio similarity, negative curvature, aberrancy, and GLCM homogeneity) and four frequency attributes (cosine of instantaneous phase, average frequency, continuous wavelet transform (CWT), spectral decomposition at frequencies of 29, 57, and 81 Hz) were extracted from the seismic volume. The approximate thickness ranges of these frequencies that correspond to velocities at the reef interval are 107, 55, and 38 feet for 29, 57, and 81 Hz, respectively.

Results
A wide range of single attributes were calculated and compared by the interpreter to decide which attributes best showed features of interest in the reef. Therefore, after grading according to the criteria previously described, the following attributes were deemed most worthy. Four structural attributes (energy ratio similarity, negative curvature, aberrancy, and GLCM homogeneity) and four frequency attributes (cosine of instantaneous phase, average frequency, continuous wavelet transform (CWT), spectral decomposition at frequencies of 29, 57, and 81 Hz) were extracted from the seismic volume. The approximate thickness ranges of these frequencies that correspond to velocities at the reef interval are 107, 55, and 38 feet for 29, 57, and 81 Hz, respectively.
The following section details observations of the single attributes results extracted onto horizons picked by an expert geologist that works in the Michigan Basin. The single attributes that were calculated are extracted onto the Brown horizon slice and initial observations were recorded. The Brown horizon is the top of the Guelph reef complex (the top of the Reef Core, Reef Apron, and Reef Talus geobodies), which made it the most important horizon for displaying the attribute maps which is the most representative of the morphology of the Guelph reef complex. Then, six total SOMs were calculated. All six SOM calculations were output with a 5 × 5 neuron grid. Three of the six SOM runs used the aforementioned structural attributes in addition to the PSTM amplitude, while the other three used the aforementioned frequency attributes without the PSTM amplitude. The amplitude was removed from the frequency SOM runs since it tended to dominate the results. Following removal, there was a greater variability in which the attribute contributed the most for each cluster of data. The first two (one structure and one frequency) SOM runs were restricted between the A-2 Carbonate and Clinton layers and utilized all the inlines and crosslines. The next two (one structure and one frequency) were restricted between the same two horizons but the inlines and crosslines were restricted to isolate the reef from the surrounding geology. The final two SOM runs (one structure and one frequency) used the same inline and crossline restrictions as the second set of SOMs but were further confined vertically to be between the A-1 Carbonate and the Gray layers. The goal of this arbitrary section is to identify areas of interest along well bores where there are permeability and porosity data. Four type wells that had the most complete suite of porosity and permeability data were chosen for this task: P-106, P-201, P-102, and P-103.  Figure 5A is a two-way time (TWT) map of the Brown horizon. The reef complex is a clear structural high compared to the off-reef areas, with a slightly higher structure observed in the northern section of the reef complex compared to the southern. The slope is gentler on the western edge of the reef and steeper on the eastern edge, which corresponds to the asymmetric reef model proposed by Rine et al., 2017 [9]. Figure 5B shows pre-stack time migrated (PSTM) amplitude data. There is a consistent negative amplitude that composes most of the reef top, followed by a ring of positive amplitude that contours the reef outline. Figure 5C shows the energy ratio similarity (ERS). The yellow arrows filled with black are showing areas of lower similarity along the eastern side of the reef, as well as one spot in the southern portion. The yellow lines outline some sinuous features, potentially talus deposits, that surround the main reef body. Figure 5D shows the negative curvature attribute. Outlined in yellow is a connected feature that runs throughout the reef body. There is also a broad trend of negative curvature that wraps around the reef outline itself everywhere except for the southern tip. The area that surrounds the reef complex is comprised of alternating carbonate (A-1 Carbonate) and halite (A-2 Salt) layers, which have a much greater acoustic impedance than the reef complex interior which is comprised of only carbonate. Figure 5E shows aberrancy. An interconnected feature is present here as well, but it takes up more of the reef. Within that interconnected feature there are small pockets of little to no aberrancy. At the bioherm level, these features might be associated with coral growth patterns within the reef. Figure 5F shows GLCM homogeneity. The yellow circle shows an area of high homogeneity which is concentrated in the northern section of the reef. There is also a low homogeneous zone surrounded by a tight ring of high homogeneity on the reef outline. The tight ring of homogeneity may be areas in which reef debris or talus collected as the reef was exposed and eroded during the sea-level fall. The yellow arrows filled with black are pointing to areas surrounding the reef which are high homogeneity concentrations. These concentrations might indicate where smaller bioherm mounds initially formed but were abandoned/drowned as the sea level rose during the time of reef growth. Figure 5G shows the average frequency. The black shapes outline areas of particularly low average frequency in the main reef body, approximately 20 Hz on average. The rest of the reef reveals a higher average frequency, approximately 70 Hz on average. The black arrows are pointing to areas of anomalously high or low average frequency which lie outside of the reef body. The varying frequencies found within the reef core potentially indicate a transition from a higher porosity and permeability to a lower porosity and permeability lithofacies. Figure 5H shows that the cosine of the instantaneous phase accentuates what the PSTM amplitude shows. There is a negative cosine value in the reef body, surrounded by a positive ring that hugs the reef body. The attribute also shows that there is another ring with a negative cosine value that tightly curves the inner features to the east and more broadly curves to the west. Outside of that ring lies mostly positive values speckled with negatives. Figure 5I-K displays the CWT spectral decomposition at frequencies of 29, 57, and 81 Hz, respectively. Here, 29 Hz corresponds to approximately 107-foot thicknesses, 57 Hz to 55-foot thicknesses, and 81 Hz to 38-foot thicknesses. Figure 5I shows that the CWT spectral decomposition at 29 Hz shows the main reef body has a lower presence of 29 Hz relative to its surroundings everywhere except where the yellow circle is in the southern portion of the reef. The presence of 29 Hz in the southern portion of the reef might indicate that it is slightly thicker than the rest of the reef. Figure 5J is the CWT spectral decomposition at 57 Hz. Around the reef there is a ring of higher presence of 57 Hz. The southern half of the ring exhibits a higher presence than the northern half. Most of the reef body is near the void of 57 Hz, except where the yellow lines are outlining. There appears to be small shapes or features that have some presence of 57 Hz within the reef itself. Going from 29 to 57 Hz, we begin to pick up on thinner features. In this case, we go from frequencies which are associated with 107-to 55-foot thicknesses. Seeing more detail within the reef itself as the frequency increases could imply that the internal growth patterns of the reef happened on a smaller scale. The ring surrounding the reef seen with the 57 Hz spectral decomposition may be associated with an area of reef debris that can be found surrounding the reef core. Figure 5K is the CWT spectral decomposition at 81 Hz. There is a seemingly random distribution at first glance within the reef outline but upon closer inspection, areas where the attribute displayed similar values are usually connected. The yellow lines outline the more drastic regions where the higher presence occurred within the reef. Looking out from the reef the same pattern can be seen. The highly interconnected areas seen with the 81 Hz spectral decomposition could be associated with the more complex internal coral growth patterns found within the reef core. Figure 5L is a traditional RGB co-blending of the three spectral decomposition volumes displayed at 430 ms. Through this co-blending the same internal features, with perhaps a little more detail, can be better visualized in order to interpret each spectral decomposition image individually. Yet, using the spectral decomposition volumes in a SOM, there is potentially more detail that can be extracted. Each of these attributes brings something different to the table that will become more important when combining them into a SOM. In addition, all SOM results output the data from the wavelet scale to the voxel scale.  After comparing the set of SOM runs, the third SOM run with the frequency attributes was next used. Figure 6 shows a comparison of the second and third SOM runs. The second SOM run was created with more input data than the third. Some notable differences are noted. The white circles show an area in the second SOM characterized by 3 neurons but in the third SOM that same area is characterized by at least 5 neurons. The white rectangle shows an area characterized by a few neurons in the second SOM and a significant amount more in the third SOM. Finally, the arrows point to areas where the general trends and patterns are the same. The gold and purple neurons in the second SOM are similar to the brown and green neurons in the third SOM. However, there are still some slight changes in these features, as well. These differences are caused by the varying amount of input data into each SOM. The third SOM was used due to the higher level of detail seen in areas of interest within the reef.

Discussion
The third SOM run with the frequency attributes was chosen to explore the reef further. The input data and parameters for the third SOM were best suited to achieve the goal of identifying areas of poor to great permeability and porosity. Similar graded ranges as used previously were used for both porosity and permeability. Horizons with high levels of confidence were used to bridge the gap between the well logs and the time volumes. A qualitative interpretation on the efficacy of the frequency SOM runs to identify subseismic areas of varying permeability and porosity was conducted using the available data. Figure 7 shows the A-A' cross section from Figure 1. The top image in Figure 7 reveals the third SOM run from the frequency attributes with the wiggle traces of the seismic amplitude overlain. Clearly the SOM results are able to define and characterize the reef on a much finer scale. The size of the voxels for this particular dataset were 55 by 55 feet laterally and 1 ms vertically.
Areas of varying permeability and porosity were located on the well logs and then the neurons associated with those locations were noted. At well P-106, the best-case scenario was found with great permeability and porosity on average. Neurons associated with this case are numbers 19 and 23. Two areas of greater than 10% (good to great) porosity were found with widely varying permeability at wells P-201 and P-102 that correspond to neuron 19. Finally, well P-103 had locations where the porosity is roughly the same, but with vastly different permeability were found. The two neurons associated with this scenario are 10 and 25. The bottom image of Figure 7 shows these aforementioned neurons which are associated with the desirable porosity and permeability in blue, while the rest of the neurons are shown in grey. The well logs are also displayed over the SOM.

Discussion
The third SOM run with the frequency attributes was chosen to explore the reef further. The input data and parameters for the third SOM were best suited to achieve the goal of identifying areas of poor to great permeability and porosity. Similar graded ranges as used previously were used for both porosity and permeability. Horizons with high levels of confidence were used to bridge the gap between the well logs and the time volumes. A qualitative interpretation on the efficacy of the frequency SOM runs to identify subseismic areas of varying permeability and porosity was conducted using the available data. Figure 7 shows the A-A' cross section from Figure 1. The top image in Figure 7 reveals the third SOM run from the frequency attributes with the wiggle traces of the seismic amplitude overlain. Clearly the SOM results are able to define and characterize the reef on a much finer scale. The size of the voxels for this particular dataset were 55 by 55 feet laterally and 1 ms vertically. Once each SOM was completed, the resultant neurons were filtered based on their attribute contribution percentages. If there were any blends similar to the successful blends developed from the Puttygut trials then they were kept, if not they were discarded. There were no exact matches with the attribute percentages in neurons calculated at Ira, but there were three results that were very close. Lenox was primarily composed of three neurons within the reef body, so the filtering and selection process was very short. Again, there were no exact matches at this reef either, but two of the three neurons that were found within the reef were very close. The three blends of attributes that were developed at Puttygut and looked for in the Ira and Lenox reefs were: (1) Areas of >50 mD permeability and >10% porosity had a blend of CWT spectral decomposition at 81 Hz, average frequency, and CWT spectral decomposition at 29 Hz in order from the highest contribution to the lowest. (2) Areas with 6% to >10% porosity had a blend of CWT spectral de- Areas of varying permeability and porosity were located on the well logs and then the neurons associated with those locations were noted. At well P-106, the best-case scenario was found with great permeability and porosity on average. Neurons associated with this case are numbers 19 and 23. Two areas of greater than 10% (good to great) porosity were found with widely varying permeability at wells P-201 and P-102 that correspond to neuron 19. Finally, well P-103 had locations where the porosity is roughly the same, but with vastly different permeability were found. The two neurons associated with this scenario are 10 and 25. The bottom image of Figure 7 shows these aforementioned neurons which are associated with the desirable porosity and permeability in blue, while the rest of the neurons are shown in grey. The well logs are also displayed over the SOM. The permeability log (left) has a logarithmic scale ranging from 1 to 1000 mD and the porosity log (right) has a linear scale ranging from 0% to 20%. Blue areas of each respective log are associated with greater than 50 mD permeability or greater than 10% porosity. The green areas are associated with 10-50 mD permeability or 6-10% porosity. Yellow areas are associated with 1-10 mD permeability or 3-6% porosity. Finally, the red areas are associated with <1 mD permeability or 0-3% porosity. Neuron 23 was only present at the one location where there is the highest values of permeability and porosity. Neuron 23 is also the only neuron with a high percentage of CWT spectral decomposition at 81 Hz at 40.6% contribution. The other top contributors for this neuron are the average frequency at 35.3% and CWT spectral decomposition at 29 Hz with a 22.4% contribution. Neuron 19 can be found where there was a wide range of permeability but intermediate to great porosity. Neuron 19 is also found next to neuron 23 at well P-106 where the best properties were found. The top three attribute contribution percentages for neuron 19 are CWT spectral decomposition at 29 Hz with 39.3%, average frequency with 29.9%, and CWT spectral decomposition at 57 Hz with 24.3%. The two neurons found at well P-103 may seem confusing at first, but upon closer inspection the attribute contribution percentages for neurons 10 and 25 are almost identical. The porosity in this well has slight variations from good to great through the whole section, while the permeability is more highly varied and does not seem to follow any specific trend. That information implies that these two neurons are mostly driven by the porosity and not as affected by the permeability. The blends of neurons 10 and 25 have attribute contribution percentages of~38% for CWT spectral decomposition at 29 Hz,~31-35% for CWT spectral decomposition at 57 Hz,~14-15% for average frequency, and~11% for CWT spectral decomposition at 81 Hz. The banding seems to come from the 3.2% difference in the cosine of instantaneous phase attribute contribution. There were also areas of undesirable rock properties that were found to be associated with two other neurons, 18 and 20. These two neurons only appear at well locations associated with poor porosity and permeability conditions. These two neurons are dominated by two attributes. Neuron 18 has a 51.2% contribution from CWT spectral decomposition at 29 Hz and a 30.5% contribution from average frequency. Neuron 20 has a 44.5% contribution from average frequency and a 39.3% contribution from CWT spectral decomposition at 29 Hz. Essentially, areas of high permeability and porosity are found where CWT spectral decomposition at 81 Hz, average frequency, and CWT spectral decomposition at 29 Hz are the highest contributing attributes. Areas that have high porosity, but widely varying permeability, were found to be characterized by two sets of attributes: (1) CWT spectral decomposition at 29 Hz, average frequency, and CWT spectral decomposition at 57 Hz. (2) CWT spectral decomposition at 29 Hz, CWT spectral decomposition at 57 Hz, average frequency, CWT spectral decomposition at 81 Hz. Finally, areas of poor porosity and permeability were characterized by attribute combinations dominated by average frequency CWT spectral decomposition at 29 Hz.
The next phase of this research applies the workflow developed over the Puttygut reef and is applied to other nearby reef complexes to test the repeatability and robustness of the method. These reefs have less porosity and permeability control, so this is a qualitative study in the end. Figure 1 shows the locations of the reef complexes in relation to each other. The following discussion will consist of a description of each reef, the parameters used to create the attributes and SOMs, and the process by which the neurons were filtered and selected. Then, an interpretation of the results and supplementary evidence of well information to confirm the claims follows.
The input data for the Ira and Lenox SOMs were the average frequency, cosine of instantaneous phase, CWT spectral decomposition at 29, 57, and 81 Hz. All these frequency attributes were created from the 3D PSTM amplitude volumes. Horizons for each reef were also provided from the A-2 Carbonate down to the Clinton. Only the A-1 Carbonate, the Brown, and the Gray horizons were used for this exercise. These attributes and horizons were imported, and the SOMs were created. The boundaries for each SOM were slightly different but followed the same principles as with the third Puttygut SOM run. They were vertically confined between the A-1 Carbonate and Gray horizons and confined laterally by a set of inlines and crosslines.
Once each SOM was completed, the resultant neurons were filtered based on their attribute contribution percentages. If there were any blends similar to the successful blends developed from the Puttygut trials then they were kept, if not they were discarded. There were no exact matches with the attribute percentages in neurons calculated at Ira, but there were three results that were very close. Lenox was primarily composed of three neurons within the reef body, so the filtering and selection process was very short. Again, there were no exact matches at this reef either, but two of the three neurons that were found within the reef were very close. The three blends of attributes that were developed at Puttygut and looked for in the Ira and Lenox reefs were: (1) Areas of >50 mD permeability and >10% porosity had a blend of CWT spectral decomposition at 81 Hz, average frequency, and CWT spectral decomposition at 29 Hz in order from the highest contribution to the lowest.
(2) Areas with 6% to >10% porosity had a blend of CWT spectral decomposition at 29 and 57 Hz in higher percentages and average frequency and CWT spectral decomposition at 81 Hz in lower, but still significant, percentages. (3) Areas with porosity ranging from 3% to >10% had a blend of CWT spectral decomposition 29, average frequency, and CWT spectral decomposition at 57 Hz in order from the highest contribution to the lowest. A close proxy to each of these cases was found at the Ira reef, while proxies to cases one and two were found at the Lenox reef. The attribute blends were surprisingly similar to what was found at the Puttygut reef.
Ira 3D is located to the south of Puttygut. The reef has an area of 100 acres and is approximately 190 feet in height with a capacity of 6.3 BCFG. The seismic data that were used as the base for all the attribute and SOM analysis was the PSTM amplitude data shot in 2018. The inline and crossline spacing was 440 feet, the shot and receiver spacing was 110 feet, and bin size was 55 by 55 feet, and the sample rate was 1 ms. All the frequency-based attributes were calculated from the PSTM amplitude seismic volume. These attributes are the cosine of the instantaneous phase, the average frequency, and the CWT spectral decomposition at 29, 57, and 81 Hz. The input parameters for the SOM mimicked the Puttygut SOM as close as possible, with the lateral extent confined to the reef as much as possible, and the vertical extent between the same two horizons (the A-1 Carbonate and the Gray horizon). After filtering through the results of the SOM, three neurons were found that had attribute contribution percentages very similar to those found at Puttygut when looking for areas of varying permeability and porosity. The right two images of Figure 8 show crossline 58 at the Ira reef SOM which goes through five wells in the area. The top image is with all the neurons displayed, while the bottom image is with the neurons that are similar to the ones found at Puttygut displayed. Neuron 6 had the highest contributions from CWT spectral decomposition at 57 Hz of 48.6%, average frequency of 33.9%, and CWT spectral decomposition at 29 Hz of 17%. Neuron 16 has the highest contributions from CWT spectral decomposition at 57 and 81 Hz, average frequency, and CWT spectral decomposition at 29 Hz of 36.2%, 26.2%, 23.4%, and 13.4%, respectively. Finally, neuron 23 has the highest contributions from average frequency of 35.3% and CWT spectral decomposition at 29 and 81 Hz of 34.9% and 28.4%, respectively. When comparing the neuron attribute blends of Ira to Puttygut, it is observed that neuron 6 from Ira is most similar to the third case at Puttygut where there is a wider range of porosities from intermediate to great but changes in permeability are not captured. Neuron 16 is most similar to case two from Puttygut where there are areas of 6% to >10% porosity, but again changes in permeability are not well captured. Finally, neuron 23 is most similar to case one at Puttygut where there are areas of >50 mD permeability and >10% porosity present. From the images on the right side of Figure 8 we see that neuron 6, which had a similar attribute contribution percentage as case three from Puttygut, was not very abundant overall. There is a higher concentration of neurons 16 and 23 within the Ira reef. Neuron 16 is localized near the center of the reef, while neuron 23 is seen more in larger clumps on the outskirts of the reef, as well as in a single location in the center of the reef, in between a section of neuron 16. The vertical white lines on the crossline image represent the locations of wells that penetrate the reef. Three of the five wells land within areas of higher porosity and permeability (ranges mentioned earlier) indicated by the active neurons, with the other two wells to the north only coming close. Those three wells, I-201, I-111, and I-108 also happened to be the three highest ranked wells in the field. The other two wells that did not penetrate the neurons but only came close, I-104 and I-202, were ranked ninth and twelfth out of twelve in the field. Lenox 3D is located to the west and slightly south of Puttygut. The reef has an area of 35 acres and is approximately 260 feet tall with a capacity of 3.2 BCFG. The seismic data that were used as the base for all the attribute and SOM analysis was the PSTM amplitude data shot in 2018. The inline and crossline spacing was 660 feet, the shot and receiver spacing was 110 feet, and bin size was 55 by 55 feet, and the sample rate was 1 ms. All the frequency-based attributes were calculated from the PSTM amplitude seismic volume. These attributes are the cosine of the instantaneous phase, the average frequency, and the CWT spectral decomposition at 29, 57, and 81 Hz. The input parameters for the SOM mimicked the Puttygut SOM as close as possible, with the lateral extent confined to the reef as much as possible, and the vertical extent between the same two horizons (the A-1 Carbonate and the Gray horizon). The filtering process for the Lenox reef was much easier, as there were only three neurons that were within the reef body, one of which did not match any of the attribute trends or relationships from Puttygut. That left two that were rela- Lenox 3D is located to the west and slightly south of Puttygut. The reef has an area of 35 acres and is approximately 260 feet tall with a capacity of 3.2 BCFG. The seismic data that were used as the base for all the attribute and SOM analysis was the PSTM amplitude data shot in 2018. The inline and crossline spacing was 660 feet, the shot and receiver spacing was 110 feet, and bin size was 55 by 55 feet, and the sample rate was 1 ms. All the frequency-based attributes were calculated from the PSTM amplitude seismic volume. These attributes are the cosine of the instantaneous phase, the average frequency, and the CWT spectral decomposition at 29, 57, and 81 Hz. The input parameters for the SOM mimicked the Puttygut SOM as close as possible, with the lateral extent confined to the reef as much as possible, and the vertical extent between the same two horizons (the A-1 Carbonate and the Gray horizon). The filtering process for the Lenox reef was much easier, as there were only three neurons that were within the reef body, one of which did not match any of the attribute trends or relationships from Puttygut. That left two that were relatively similar to the cases discovered at Puttygut. Neuron 23 had the highest contributions from CWT spectral decomposition at 81 and 29 Hz, and average frequency of 34.1%, 30%, and 25.3%, respectively. In addition, neuron 24 which had the highest contributions from CWT spectral decomposition at 57 and 29 Hz, average frequency, and CWT spectral decomposition at 81 Hz of 29.1%, 23.5%, 23.4%, and 15.7%, respectively. It is very clear from the left image in Figure 8 that the dominant neurons within this reef are numbers 23 and 24. When comparing the neuron attribute blends of Lenox to Puttygut it can be seen that neuron 23 from Lenox is most similar to case one at Puttygut where there are areas of great permeability and porosity present, and neuron 24 at Lenox is most similar to case two from Puttygut where there are areas of good to great porosity, but changes in permeability are not well captured. The left image in Figure 8 displays the frequency SOM through the center of the reef at inline 56. It is very clear that neuron 24 is most of the reef interior, while neuron 23 makes up more of the reef boundary. There is some slight intrusion of another neuron, number 9, in the reef core, but its attribute contribution percentages did not closely resemble anything from the three Puttygut cases, so it was excluded from further analysis.
These results hold promise for future applications of machine learning to carbonates. The relationship found between the neurons and the porosity and permeability information demonstrates that this method can successfully locate areas of desirable rock properties. There are, however, some limitations with the current workflow. The main limitations being that the workflow was created from one reef, the Puttygut reef. Subsequently, that workflow was only tested on two other nearby reefs, the Ira and Lenox reefs. Reefs reservoirs are created by many different processes throughout the world, so it would be ideal to test this developed workflow in a new location and then to tailor it slightly to suit the needs of the geology in the area. Those changes to the workflow could include, but are not limited to, increasing or decreasing the amount of input data into the SOM, altering or changing the input attributes or potentially looking for different attribute contribution percentages in order to represent areas of higher porosity and permeability.

Conclusions and Future Work
From this work, there does seem to be a trend or relationship between porosity and permeability in the seismic/well data and the SOM calculations (see Figure 7). When comparing the well logs at the Puttygut reef to the SOM results, areas with >50 mD permeability and >10% porosity present coincide with attribute blends with high percentages of CWT spectral decomposition at 81 Hz, average frequency, and CWT spectral decomposition at 29 Hz. The core profile for well P-106 associates this depth with a heavily karsted reservoir zone. Other attribute combinations were also found that were mainly affected by porosity and not so much by permeability. Areas of 6% to >10% porosity were associated with attribute blends of CWT spectral decomposition at 29 and 57 Hz in higher percentages and average frequency and CWT spectral decomposition at 81 Hz in lower, but still significant, percentages. The relationship between those four attributes and the porosity was discovered at well P-103 and stretched from the A-1 Carbonate down almost to the base of the bioherm layer. These two neurons exhibit a similar vertical extent in multiple locations throughout the reef, creating pillars of desirable porosity. Finally, areas with a wider range of porosity, from 3% to >10%, had attribute blends in decreasing percentages of CWT spectral decomposition 29, average frequency, and CWT spectral decomposition at 57 Hz. These locations were near the upper section of the bioherm. Areas with attribute blends dominated with only CWT spectral decomposition at 29 Hz and average frequency should be avoided as they were associated with areas of lower porosity and permeability. A Puttygut reef model created by Garrett 2016 [42] was compared to the four best porosity and permeability case neurons from the Puttygut SOM from this study. It is possible to make a qualitative comparison between the chosen neurons and the geology of the reef. The attributes associated with neuron 23 seem to correspond to areas where the reef talus is present. Neuron 19 corresponds to areas where the reef core is located. Neurons 10 and 25 correspond to areas where the reef core and bioherm are located. Garrett 2016 [42] also created the same Puttygut model populated with calculated porosity and permeability values. It can be seen that areas of higher porosity and permeability also correlate well with the chosen neurons from the SOM created from this workflow.
After repeating this study at the Ira and Lenox reefs, it became apparent that there was a relationship between certain attribute contribution percentages and rock properties. While there is no case between the attribute contribution percentages that are exactly the same, there is most definitely a consistent trend/relationship when looking solely at the highest contributing attribute combinations. For the case of >10% porosity and >50 mD permeability, CWT spectral decomposition at 81 Hz, average frequency, and CWT spectral decomposition at 29 Hz were always the highest contributors by a significant margin. In the second case where there is a range of 6% to >10% porosity but no noticeable effect on permeability, every attribute significantly contributes except for the cosine of the instantaneous phase. Finally, in the case where a wider range of 3% to >10% porosity are found, CWT spectral decomposition at 29 Hz, average frequency, and CWT spectral decomposition at 57 Hz are found to be the most prevalent contributors. From the results shown, it seems that the workflow created at Puttygut was able to translate over the Ira and Lenox wells. The similar attribute blends being found at similar geological locations within each reef increase the confidence that this workflow was able to accurately distinguish areas of higher porosity and potentially permeability.
In conclusion, this described SOM set up and search for neurons where the attribute contributions are heavily skewed towards CWT spectral decomposition at 81 Hz, average frequency, and CWT spectral decomposition at 29 Hz, provide the best quantification of permeability and porosity. While this workflow was only tested in a localized area, it is an initial demonstration on how the multi-dimensional, multi-attribute analysis can be applied to reefs. Reefs are formed by different processes depending on where they are geographically located. As such, different attributes or workflow modifications may be needed as this method is applied to other reefs around the world.