Three-dimensional morphometric analysis and statistical distribution of the Early Kimmeridgian Hanifa Formation stromatoporoid/coral buildups, central Saudi Arabia

The high porosity and permeability stromatoporoid/coral facies partly characterize some of the most productive supergiant Late Jurassic carbonate reservoirs in the Middle East. Reservoir models often overlook detailed morphology and distributions of this facies due to the limited resolution of subsurface data. 3D morphological architecture and distribution of the Late Jurassic stromatoporoid/coral buildups are poorly understood as not many studies specifically investigate these aspects. This study performs a full three-dimensional outcrop investigation of the stromatoporoid/coral complex part of the Late Jurassic Hanifa reservoir analog in Wadi Birk, Saudi Arabia. We investigated spatial correlation, clustering tendencies, scaling relationships, and comprehensive statistical description over a 1.2 × 1 km area analog to a small-scale sectoral model of a typical Middle Eastern oilfield. This study utilizes integrated datasets including measured sections, drone-based digital outcrop model (DOM), core, ground penetrating radar (GPR), and shallow seismic to perform a full three-dimensional outcrop investigation. We identify and map three levels of the scaling hierarchy of the buildups: S1, S2, and S3. S1-scale buildups (a few meters in length) combine to form S2 clusters (~35 m in length), which in turn may grow into larger S3 complexes (~128 m in length). S1-scale buildups are near-circular and isotropic. However, as the S1 buildups cluster and grow into S2 and S3, their circular shape evolves into pseudo-ellipsoids with a strong preferential orientation towards northwest-southeast. The S2 and S3-scale amalgamated buildups are associated with asymmetric flank accumulation towards the northwest. We propose that the orientation of the S2 and S3-scale buildups and associated flanks are attributed to hydrodynamic currents driven by paleo-trade winds. This study also performs a brief static connectivity assessment to demonstrate the subsurface implications of these facies to field development planning. Results show that a vertical 5-spot pattern with a 1 km spacing can, in the best case, only access 14% of the Gross Rock Volume of the build-up facies due to the extreme lateral heterogeneity. Some 84% of the build-ups would be connected via flank and other facies with significantly different reservoir quality. This result underlines the extreme heterogeneity of the late Jurassic shallow-water platform carbonates and the potential for substantial quantities of bypassed hydrocarbons remaining in subsurface analog reservoirs.


Introduction
Stromatoporoids, a hypercalcified benthic demosponge, and corals are prominent metazoans reef builders in the Middle Palaeozoic (Ordovician -Late Devonian) and the Jurassic (Scholle and Ulmer-Scholle, 2003;Leinfelder et al., 2005;Kershaw et al., 2013Kershaw et al., , 2018. In the Arabian Late Jurassic strata, stromatoporoids are usually co-occurring with scleractinian corals and form two styles of stratigraphic reefs, biostromes and bioherms (Leinfelder et al., 2002(Leinfelder et al., , 2005Hughes et al., 2008;Fallatah and Kerans, 2018;Al-Mojel et al., 2020). Interpretation of growth characteristics and morphologies as a proxy for assessing environmental factors are primarily similar between stromatoporoids and their associated scleractinian corals (Nose and Leinfelder, 1997;Leinfelder et al., 2002Leinfelder et al., , 2005. Moreover, the Late Jurassic Arabian stromatoporoids are interpreted to have been adapted to overheated and slightly hypersaline waters (Leinfelder et al., 2005). Therefore, investigating morphologies of the Late Jurassic Arabian stromatoporoids and their associated corals might give an insight into the adaptability of the past reef community and their interaction with in-situ environmental conditions. Laterally heterogenous high porosity and permeability stromatoporoid/coral facies partly characterize giants and supergiants hydrocarbon carbonate reservoirs in the Middle East (Al-Khalifah and Makkawi, 2002;Afifi, 2004;Lindsay et al., 2006). This stromatoporoid/coral facies exhibits sub-seismic interwell meter-scale depositional heterogeneities often overlooked or overly homogenized by most subsurface models. The subsurface data neither have the resolution nor data density to resolve this level of facies heterogeneity. A field development plan based on an overly "homogenous" model will lead to sub-optimum ultimate recovery due to underestimating bypassed oil, poorly predicted water breakthrough dates, and incorrect estimation of recovery rates from simulations (Vahrenkamp et al., 2019). Including the inter-well depositional heterogeneity of the stromatoporoid/coral facies in the reservoir models are even more critical for CO 2 -based Enhanced Oil Recovery (EOR).
Despite their importance as a reservoir component, almost no study explicitly investigates the 3D distribution and morphology of the Late Jurassic Arabian stromatoporoid/coral buildups. Fallatah and Kerans (2018) and Al-Mojel et al. (2020) performed a regional (~300 km) sequence stratigraphy study of the Oxfordian-Early Kimmeridgian Hanifa Formation in central Saudi Arabia. Both studies reported the presence of stromatoporoid/coral biostromes/bioherms within the Ulayyah member of the Hanifa Formation. El-Sorogy et al. (2018) performed depositional architecture and sequence stratigraphy study of the Ulayyah Member of the Hanifa Formation in Al-Abbakayn, Sadous, and Maashabah mountains, Saudi Arabia. However, their studies are limited to vertical stacking patterns and two-dimensional facies profiles from interpolated measured outcrop sections extrapolated into 3D. Elzain et al. (2020) performed lithofacies modeling and constructed an outcrop-based geocellular model of the upper Ulayyah reservoir unit of the Hanifa Formation in Jabal Abbakayn, central Saudi Arabia. However, their facies distribution, including the stromatoporoid/coral facies, is propagated stochastically using sequential indicator simulation based on several 1D measured sections. Observation-based 3D statistical distribution of stromatoporoid/coral facies was not reported or explicitly modeled in either study. Ground penetrating radar (GPR), along with the digital outcrop model (DOM), has been utilized as a methodology to perform three-dimensional outcrop investigation (Adams et al., 2005;Asprion et al., 2009;Nielsen et al., 2009;Jorry and Bievre, 2011;Janson et al., 2015). Outcrop exposures of the Hanifa Formation in Wadi Birk (coordinate: 23.28124 • ; 46.64420 • ), central Saudi Arabia, were included in Al-Mojel et al. (2020) regional studies. They reported that the expression of stromatoporoid/coral facies in Wadi Birk is among the most prominent. The Hanifa Formation outcrop in this wadi is also part of the multi-attribute studies by Ramdani et al. (2021). Therefore, this wadi provides an analog to document and investigate km-scale spatial distribution and 3D morphology of the Late Jurassic Arabian stromatoporoid/coral facies.
The objective of this study is to quantitatively investigate and document the 3D morphology and spatial distribution of the Late Jurassic stromatoporoid/coral facies part of the Ulayyah member of the Hanifa Formation in Wadi Birk, Saudi Arabia (Fig. 1A). We investigate spatial correlation, clustering tendencies, and scaling relationships of this facies over a 1.2 × 1 km field area analogous to a small-scale sectoral model of typical Middle Eastern carbonate reservoirs. This study reports statistical findings and discusses relationships between this facies' 3D organizational and morphological heterogeneity with paleohydrodynamic energy. This study also performs a brief static connectivity assessment to demonstrate the subsurface implications and the

AI
Acoustic impedance DOM Digital outcrop model EOR Enhanced oil recovery GNSS-IMU Global navigation satellite system -Inertial measurement unit GPR Ground penetrating radar GRV Gross rock volume HCS Hanifa composite sequence  (Ziegler, 2001;Scotese, 2003;Markello et al., 2008;Fallatah and Kerans, 2018). The red circle shows the location of Wadi Birk (study area) during this time. (1C) Regional north-south sequence stratigraphic correlation (~370 km) of the Hanifa Formation from Al-Mojel et al. (2020) where the Hanifa formation is divided into two 3rd order sequences: Hanifa Composite Sequence-1 (HCS-1) and Hanifa Composite Sequence-2 (HCS-2). The dashed rectangle highlights the location of the study area (Wadi Birk) within this sequence stratigraphic framework. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.) importance of correctly assessing the architecture of stromatoporoid/ coral facies in the field development.

Geological and stratigraphic context
The Late Jurassic Hanifa Formation (Fig. 1A) is exposed within the 800 km-long Tuwaiq Mountain Escarpment in central Saudi Arabia with a regional structural dip of ~1 • to the east (Powers et al., 1968). A sharp lithological contrast between cliff-forming Tuwaiq Mountain limestone and the base Hanifa shale forms the basal contact of the Hanifa Formation in the outcrop (Al-Mojel et al., 2020). An extensive erosional and iron-stained surface marks the top contact between the Hanifa Formation and the overlying Jubaila Formation (Powers et al., 1968;Al-Mojel et al., 2020). This contact is also associated with abrupt facies changes from coral/stromatoporoid-dominated facies of the upper Hanifa Formation to dark-colored quartz-rich grainstones found at the base of the Jubaila Formation (Al-Mojel et al., 2020). The Hanifa formation in the outcrop is lithostratigraphically divided into two members (Vaslet et al., 1983): the lower Hawtah and the upper Ulayyah members (Fig. 1A). Calcareous shale interbedded with multiple carbonate grainstone beds and a thin coral/stromatoporoid bearing layer (Fig. 1C) characterizes the Hawtah member (Al-Mojel et al., 2020). Carbonate grainstones with coral/stromatoporoid biostromes/bioherms layer at the top ( Fig. 1C) characterize the Ulayyah member (Al-Mojel et al., 2020). The upper part of the Ulayyah member in Wadi Birk is notably characterized by coral/stromatoporoid biostromes/bioherms layer associated with oncolithic grainstone-packstone (Fig. 1C).
During Late Oxfordian to Early Kimmeridgian, the Arabian Plate was situated in the northeastern part of Gondwana (Scotese, 2003), facing the western part of the Neo-Tethys Ocean on a relatively stable passive margin (Fig. 1B). The Arabian Plate during this time was located around the equatorial area with arid to semi-arid paleoclimate (Scotese, 2003;Al-mojel, 2017). A vast shallow-water epeiric sea flooded most of the Arabian Plate in this period, leading to the formation of an extensive (~10 3 km) shallow-water marine carbonate platform (Ziegler, 2001;Lindsay et al., 2006). Concurrently, differential intraplate subsidence partly enhanced by peripheral aggradation of shallow-marine carbonates led to the development of the Gotnia, Arabian, and Rub'Al-Khali intrashelf basins (Ziegler, 2001;Lindsay et al., 2006). The study area (Wadi Birk) was located approximately 10 • -15 • south of the paleo-equator and ~1000 km inwards (Fig. 1B) from the Neo-Tethys continental margin (Ziegler, 2001;Al-Mojel et al., 2020). The prevailing trade wind vector was approximately northwest-southeast in this location (Thierry et al., 2000;Markello et al., 2008;Fallatah and Kerans, 2018).
Depositional settings and sequence stratigraphy of the Hanifa Formation are subjects of ongoing discussions (Hughes et al., 2008;El-Sorogy et al., 2018;Fallatah and Kerans, 2018;Al-Mojel et al., 2020). This study follows depositional settings and sequence stratigraphic architecture of the Hanifa Formation proposed by Al-Mojel et al. (2020). In this model, the Hanifa Formation was deposited on a flat-topped shelf in an inner-platform setting and vertically divided into two third-order sequences (Fig. 1C): The Hanifa Composite Sequence-1 (HCS-1) and the Hanifa Composite Sequence-2 (HCS-2). The entire Hawtah Member is a part of HCS-1, while the Ulayyah Member includes the upper part of HCS-1 and the whole HCS-2. The interval of interest in this study is HCS-2 which represents the upper part of the Ulayyah member and is analogous to the subsurface Hanifa reservoir (Hughes et al., 2008;Al-Mojel et al., 2020). An interval of laterally heterogeneous stromatoporoid/coral facies dominates the upper ~25 m part of HCS2 (Fig. 1C).

Field data
This study utilizes outcrop sedimentological sections, aerial drone surveys, GPR, seismic, and core data from the previously described study area (Fig. 2). We drilled a 50 m-long core (~4.5 inches in diameter) in the study area, with the wellhead located ~2 m below the Jubaila-Hanifa contact stratigraphically. A total of ~200 m of sections were measured in the outcrop and the core. The sections started from slightly below the base of HCS-2 up to the Hanifa-Jubaila contact. This study also utilizes two sections located ~1.5 km west and ~1.7 km southeast to provide a general overview of the facies heterogeneity away from the main focus area. The sections include a description of depositional texture (Dunham and Ham, 1962;Embry and Klovan, 1971), sedimentary structures, bioturbation, bedding contact, and observable macrofossils. We calibrated the sections with 242 thin sections made from hand samples and core plugs. The thin sections are impregnated with blue-dyed resin to visualize porosity and stained with Alizarin red to differentiate calcite and dolomite (Dickson, 1966). Measured sections also incorporate spectral gamma-ray measurements using Radiation Solution™ RS-230 gamma-ray spectrometer. We measured the percentage of K (Potassium), ppm of U (Uranium), and ppm of Th (Thorium) at 50 cm spacing from the outcrop and 20 cm spacing from the core. We calculated the total gamma-ray count in the API unit using the empirical formula of Luthi (2001). Core thin sections and gamma-ray data are available in Ramdani et al. (2022a).
We collaborated with FalconViz™ (Thuwal, Saudi Arabia) to conduct aerial photography surveys over the entire study area using Wingtra™ fixed-wing, DJI™ M600 hexacopter, and DJI™ Phantom 4 Pro quadcopter drones. The surveys acquired more than two terabytes of outcrop images with 90% overlaps. The images are georeferenced by ground control points and onboard GNSS-IMU instruments connected with a Real-time Kinematic GPS. We constructed the digital outcrop model (DOM) using Pix4D® software, following the method reported in Khanna et al. (2020) and Menegoni et al. (2019). For data interpretation and further analysis, we imported the DOM into VRGS software®. The processed drone datasets are available in Khanna et al. (2022).
This study utilizes GPR and seismic data as the geophysical methods to image the stromatoporoid/coral facies "behind-the-outcrop". We acquired eight km-long 2D 50 MHz/100 MHz GPR lines and three grids (60 m × 50 m; 50 m × 20 m; 50 m × 36 m) of 3D 50 MHz GPR. We utilized PulseEKKO™ Smart Cart Acquisition System, where both transceiver and receiver antennas are mounted on a draggable cart allowing for rapid data acquisition. The spatial sampling of GPR acquisition is 0.1 m for 100 MHz and 0.25 m for 50 MHz data to avoid spatial aliasing. We processed both 50 MHz and 100 MHz GPR data through a similar processing flow (Table 1) with different parameters representing their antenna frequency in Reflexw™. The raw GPR data is available in Ramdani et al. (2022b).
We recorded three 2D seismic lines that vary in length, receiver, and shot intervals in the study area. The longest profile is a 300 m-long 120 shots-120 receivers line with 2.5 m shot and receiver interval. The shortest profile is a 144 m 96 shots-96 receivers line with a 1.5 m shot and receiver interval. The seismic source is a 500 kg accelerated weight drop mounted on a 4 × 4 extreme utility vehicle, while the receivers are 30 Hz vertical P-wave geophones. The recording length for all lines is 500 ms with a 0.5 ms sampling interval and a dominant frequency of ~120 Hz. The raw seismic data is available in Ramdani et al. (2022c). We processed the seismic reflection data (Table 1) in Reflexw™. In addition to the regular reflection processing, we also performed colored inversion (Lancaster and Whitcombe, 2000) to convert reflection amplitude into acoustic impedance. Colored inversion offers some advantages over regular seismic reflection amplitude, such as removing the wavelet smearing effect and enhancing thin beds and discontinuities (Lancaster and Whitcombe, 2000) applicable to our objective.

Laboratory measurements
This study utilizes laboratory-scale acoustic velocity, density, and relative dielectric permittivity measurements from the core and hand samples of the same study area reported in Ramdani et al. (2022d). The acoustic velocity was measured using an ultrasonic measurement apparatus integrated into Epslog's Wombat™ core scratch tester. A rectangular steel cutter coated with the polycrystalline diamond compact is "scratching" the entire 50 m core length, creating continuous 1 cm-wide and 5 mm-deep, flat grooves on the surface of the core. The acoustic velocities were measured along this groove using piezoelectric crystal probes capable of transmitting and receiving 300 KHz acoustic waves at 4 cm resolution. Bulk densities were measured from 192 cylindrical core plugs using the Metarock PDP-300™ helium gas porosimeter. Using the open-ended coaxial probe technique, Ramdani et al. (2022d) measured the relative dielectric permittivity of 24 selected hand samples and core plugs (Zheng and Smith, 1991) at a 50 MHz measurement frequency. They calibrated the apparatus using short circuit measurement, de-ionized water at 24 • C (room temperature) using the Cole-Cole equation (Cole and Cole, 1941), and permittivity measurement of air. We utilized acoustic velocity and density of the core to perform the well-seismic-tie and colored inversion. The relative dielectric permittivity measurements are mainly utilized to determine GPR velocity during data acquisition and processing.

Sedimentary facies of the Hanifa Composite Sequence-2 (HCS-2) in the study area
This study adopts and expands the established sedimentary facies classification within HCS-2 by Al-Mojel et al. (2020) and Al-Mojel and Razin (2022). Six distinct facies, including three sub-facies (eight in total), are described based on their depositional texture, sedimentary structures, bioturbation, components, and geophysical properties (Table 2, Figs. 3 and 4).
• FC1: Sharp-based intraclastal peloidal skeletal grainstone FC1 is characterized by ~10-30 cm tabular cross beds grainstone crudely stratified with poorly-sorted sub-rounded muddy intraclastal with occasional imbrication structure at the base (Figs. 3 and 4). This facies is occasionally bioturbated at the top, associated with parallel laminated beds of ~10-20 cm. The main components of this facies are poorly-sorted fine to medium-sized peloids, granule to pebble-sized mud intraclast, reworked clasts of brachiopods, gastropod shells (Nerinea sp.), and echinoderms.
• FC3: Cross-bedded composite-grain and peloidal grainstone This facies is characterized by ~10-~30 cm thick trough crossbedded grainstone with densely-packed medium-sized compositegrains as the primary components (Figs. 3 and 4). The composite grains are mostly aggregate of peloids, skeletal mollusks, and benthic foraminifera (Nautiloculina oolithica, Kurnubia palastiniensis, and Redmondoides lugeoni) that agglutinated into larger grains. The main components are accompanied by fragments of brachiopods, gastropods, and echinoderm plates as minor components. Occasional hummocky crossstratification is also observed in the slabbed core and the outcrop.
FC4A is in-situ stromatoporoids and corals bioherms developing buildup morphologies (Figs. 3 and 4). This facies is also associated with brachiopods, gastropods, lagoonal benthic foraminifera (Nautiloculina oolithica, Redmondoides lugeoni, Quinqueloculina sp), and echinoderm. FC4A represents the "core" of the stromatoporoid/coral buildup. FC4B is characterized by the abundance of pebble to cobble-sized stromatoporoid and coral detritus floating within wackestone to mud-lean packstone matrix (Figs. 3 and 4). The skeletal matrix components are lagoonal benthic foraminifera (Nautiloculina oolithica, Kurnubia palastiniensis, Redmondoides lugeoni), sponge spicules (clusters and Table 2 Depositional Facies identified in HCS2. The color code utilized in this table serves as facies color legend for all figures in this paper. Please refer to the text for a detailed explanation of each facies.. monoaxon), echinoderms, and minor clasts of dasyclad algae. This facies onlaps the flank of the buildup core (FC4A) and downlaps onto the previously deposited strata. FC4B is occasionally bioturbated and grades into FC1 away from the buildup core.

Stromatoporoid/coral facies mapping
This study interprets and maps the outline of stromatoporoid/coral facies (FC4A, FC4B, and FC4C) in the spatially georeferenced digital outcrop model (DOM) calibrated with outcrop measured sections ( Fig. 5A-C). As defined by Janson et al. (2015), the buildup core is expressed as massive, bulbous, and structureless with a lack of bedding features compared to the surrounding relatively smooth-texture and bedded strata ( Fig. 6A and B). The stromatoporoid/coral floatstones (FC4B and FC4C) are characterized by strata that exhibit small-scale internal granular bulbous features corresponding to the floating stromatoporoid/coral detritus. The floatstones onlap the buildup core as flanks and downlap onto the previously deposited strata (Figs. 5C and 6B).
This study utilizes the methodology outlined in Ramdani et al. (2022d) to assess, process, interpret, and integrate outcrop data (i.e., drone imageries and sedimentological sections) with GPR and seismic surveys in the same study area. The radargrams and acoustic impedance (AI) profiles calibrated with core data are utilized to map the stromatoporoid/coral facies outline "behind-the-outcrop" (Fig. 7A-G). The radargram facies of the stromatoporoid/coral buildups is convex to domal outline, chaotic, with or without internal structure/bedding reflection patterns (Fig. 7A -"D" symbol). Occasionally, convex reflection within an otherwise concordant background also characterizes the radargram expression of this facies (Fig. 7A -"C" symbol). This radargram facies expression is analogous to the studies done by Neal (2004) (2011), and Janson et al. (2015). In the 3D GPR depth slices (Fig. 7C, D, and 7G), the buildup is circular to pseudo-ellipsoid anomalies that outline internal chaotic reflection patterns. The stromatoporoid/coral facies in the acoustic impedance (AI) profiles ( Fig. 7B and F) are seen as patches of tabular to convex geometry with anomalously high acoustic impedance This study utilizes lab-measured acoustic velocity and relative dielectric permittivity of facies associated with the buildup complex (Table 2) performed by Ramdani et al. (2022d). These properties are utilized to estimate the maximum depth of penetration and resolution of the geophysical methods. The maximum penetration depth of the GPR is ~25 m for 50 MHz and ~13 m for 100 MHz data. The vertical resolution limit of the stromatoporoid/coral buildups is ~45 to 30 cm in 50 MHz radargram and ~15 to 12 cm in 100 MHz radargram. The vertical resolution limit of the stromatoporoid/coral complex in seismic data is at ~6.5 m. This maximum depth penetration and resolution are considered when interpreting the geophysics data.

Ripley's K analysis
Ripley's K analysis (Ripley, 1981;Dixon, 2002) is a spatial statistic method commonly applied to investigate the clustering tendencies of spatially distributed point features. The method examines how the distribution of points changes with distance by comparing points located inside the circles of variable radius (centered at each data point) with points from a completely random distribution (Khanna et al., 2020). This method was utilized in outcrop studies such as by Khanna et al. (2020) and Jacquemyn et al. (2018) to investigate the clustering tendencies of microbial buildups and rudstone-filled scoured channels. This study employs a similar approach to calculate Ripley's K-function of the mapped buildup centroids. The results are interpreted as the distance at which the buildups are likely clustered or dispersed. The distance is also employed to investigate the association between clustering and scaling relationship.

Static connectivity
Static connectivity is a simple method to infer reservoir unit association with some entities, typically wells or sets of perforations (Snedden et al., 2007). This study performs non-directional static connectivity, calculating the number and volume of connected cell clusters corresponding to one or more facies codes with vertical wells. The connected clusters indicate the likelihood of accessible volume and are a proxy for heterogeneity.

Brief vertical facies evolution of the Hanifa Composite Sequence-2 in the study area
In the Al-Mojel et al. (2020) sequence stratigraphic model, HCS-2 represents a third-order composite sequence with a long overall transgression followed by a short regression sequence. The HCS-2 in the study area (Fig. 4) is initiated by the deposition of intraclastal peloidal skeletal grainstone facies (FC1). The abundant presence of peloids and mud intraclastal, low diversity of well-preserved skeletal components (Fig. 3) and cross-bedding structure of FC1 indicate a transgressive lag deposit (Al-Mojel et al., 2020). This deposit overlay an iron-stained surface that   (Fig. 4) is also high (~30-40 API), indicating closer proximity with the terrestrial argillaceous influx (Ehrenberg and Svånå, 2001). The early transgressive lag deposit is overlain by mud-dominated FC2 intercalated with grainstones (FC3 and FC6). The association between these facies with a rich diversity of lagoonal benthic foraminifera (Hughes et al., 2008;Al-Mojel et al., 2020) suggests a general lagoonal depositional setting. West-east sedimentary log sections across the study area (Fig. 4) show limited lateral facies variations in the lower ~20 m of HCS-2. The facies variation at this level is not as extreme as the top 25 m of HCS-2, which comprises the stromatoporoid/coral complex.
Laterally heterogenous stromatoporoid/coral complex (FC4) capped by oncoidal and skeletal-foraminiferal grainstone (FC5 and FC6) characterizes the upper ~25 m of HCS-2 (Fig. 4). The stromatoporoid/coral complex of the HCS-2 in the study area is associated with the backbarrier patch reef environment and part of the general backstepping architecture of the stromatoporoid/coral bodies on the platform (Al-Mojel et al., 2020). The first appearance of stromatoporoid/coral facies is preceded by a negative shift of gamma-ray value (Fig. 4), which indicates a cleaner shallow marine carbonate environment devoid of argillaceous influx (Ehrenberg and Svånå, 2001). In general, the gamma-ray value of the stromatoporoid/coral complex, including the contemporaneous FC1 and the overlain oncoidal and skeletal-foraminiferal grainstone, is the lowest of all HCS-2 (Fig. 4). This observation suggests that the stromatoporoid/coral complex is located farthest from the shoreline (Ehrenberg and Svånå, 2001) and likely the most distal facies encountered in the study area. This proposition is also supported by Al-Mojel et al. (2020), that placed the maximum flooding surface of HCS-2 and the entire Hanifa Formation within this interval. The highstand system tracts within HCS-2 are indicated by an increase (both in number and diversity) of lagoonal benthic foraminifera within the stromatoporoid/coral complex, signaling the progradation of lagoon into the back reef environment. The abundant appearance of oncoid grains (Fig. 4) also points to the limitation of accommodation space in this area (Védrine et al., 2007;Al-Mojel et al., 2020). The oncoidal packstone-grainstone facies (FC5) is intercalated with cross-bedded peloidal skeletal-foraminiferal grainstones (FC3), forming a larger-scale foreset.
An extensive iron-stained hardground surface marks the boundary between the top of Hanifa and the base of the Jubaila Formation in this area (Al-Mojel et al., 2020). The gamma-ray values of HCS-2 around this boundary remain low (Fig. 4), indicating that the shoreline was very far from this area or the regressive part of HCS-2 has been eroded (Al-Mojel et al., 2020). Due to its extreme lateral heterogeneity, the high-frequency cycles within the stromatoporoid/coral complex are challenging to define solely based on 1D sedimentary log sections. Hence, we perform a physical correlation from the DOM. Three high-frequency sequences (cycle-1, 2, and 3/Top of Hanifa) are identified within the stromatoporoid/coral complex (Fig. 5A) and utilized for the subsequent investigation.

Stromatoporoid/coral buildup scaling relationships and parameterizations
We observe a scaling relationship between the different sizes of the buildups in the outcrop and DOM ( Fig. 5B and C). A buildup cluster with 10's of m in length is constructed from an amalgamation of smaller meter-scale buildups. A buildup cluster with 10's of meter-scale in length may also become a part of an even larger 100's of meter complex buildup assembly. Some buildup clusters grow from multiple smaller initiation points that, with time, merge with the adjacent clusters forming a more extensive complex (Fig. 6A and B). These scaling relationships are also observed in the GPR images (Fig. 7A, E, 7C, and 7G). Radargram facies of the stromatoporoid/coral buildups show smaller buildups congregating into a larger cluster. An extensive circular to pseudo-ellipsoid envelope outlines multiple buildup clusters forming a larger buildup complex.
We identify three levels of scaling relationship for the stromatoporoid/coral buildups (Fig. 8A, B, and 8D). S1, dubbed small buildup, is the smallest mappable buildup unit (meter to 10's of meter in length). The boundary between one S1 buildup with another is sometimes difficult to differentiate as these buildups tend to grow very closely or build competitively on top of each other. S1 buildup might also comprise even smaller sub-S1 scale buildup units that are too difficult to be mapped individually in this study. We mapped 628 S1-scale buildups across three high-frequency cycles within the stromatoporoid/coral complex. S2, dubbed buildup cluster, is an amalgamation of multiple S1 units constructing 10's of meter-sized buildup structures. S2 is the easiest mappable buildup in the GPR data since there is a clear radargram facies distinction between the interbuildup strata with the S2 buildup cluster. We mapped 140 S2-scale buildups (Fig. 8C) across three high-frequency cycles within the stromatoporoid/coral complex. Almost all S1 buildups are part of S2 clusters. Individual solitary S1 units are rarely found in the datasets. Buildup clusters mapped right below the cycle boundary exhibit flatter, thinner, and more elongated morphology.
S3, dubbed buildup complex, is the largest mappable buildup unit within the study area and the only buildup scale identifiable in the acoustic impedance profiles (Fig. 7B and F). S3 buildup complex evolves from multiple S2 clusters that grow vertically and merge from different initiation points (Fig. 8A and B). We mapped 37 S3 buildup complexes within the HCS-2. These 37 mapped S3 buildup complexes are assembled from 118 S2 buildup clusters, leaving 22 solitary S2 clusters that "failed" to merge into an S3 complex. Most solitary S2 units are contemporary with the oncoidal grainstone facies (FC5) in the last HCS-2 cycle. It is important to mention that the S3 complex is likely not the largest buildup scale in HCS-2. Instead, it is the largest buildup scale that can be mapped within the areal extent of our 1 km × 1.2 km dataset (Fig. 3). We speculate that S4 (and even larger scales) likely exists outside our mapping range, forming 100s of m to km-length buildup complexes. However, outcrop datasets larger than ~1 km 2 are needed to identify and map this buildup scale.
The buildups exhibit semi-spherical to pseudo-ellipsoid morphology ( Fig. 8C and D). Therefore, we employ four parameters: RL, RS, R3, and Azimuth, to simplify and describe the morphology of mapped buildups (Fig. 8E). RL and RS are the ellipsoids' long and short axis lengths. Azimuth is the angle of RL measured clockwise relative to the north. For example, a circular buildup will have a similar length of RL and RS and uniform Azimuth orientation. R3 is the maximum thickness of the buildup measured from the lowest basal contact elevation up to the maximum peak of the buildup (Fig. 8F). Most buildups, especially the S2 clusters, have a relatively flat base in a 1:1 vertical to lateral scale radargram display and allow easy measurements of the R3. On the other hand, many S3 complexes have rugged basal contact as most of them are initiated from multiple S2 clusters that merge into S3. We measure the R3 from the lowest elevation point of S2 clusters that form the S3 complex for such buildups. Notice the consistency in morphology between buildups interpreted in the GPR and DOM data. Three scales of buildups are identified: S1, S2, and S3. (8C) Select S2 buildups interpreted in the study area across three high-frequency cycles. (8D) Illustration of a unit of S3 complex that develops from clusters of S2 and S1 buildups. (8E) Pseudo-ellipsoid buildup parameterization. RL: long axis length, RS: short axis length, Azimuth: Azimuth of long axis (8F) Illustration of R3 parameter (thickness).

The direction and orientation of the flank strata
The radargram facies clearly shows the stromatoporoid/coral buildup flank strata (Figures as they are primarily identifiable from their prominent onlap towards the buildup core and downlap stratal termination away from the core (Fig. 9C and D). These flanks consist of stromatoporoid/coral floatstone facies (FC4B and FC4C). However, these stratal termination features are mostly below the resolution of the acoustic impedance (AI) profiles ( Fig. 9A and B). Instead, in the AI image, these flanks are inferred from high AI envelopes extending from the buildup of core anomalies (Ramdani et al., 2022d). A closer examination of AI sections parallel and perpendicular to the regional depositional dip shows an asymmetric AI envelope ( Fig. 9A and B). The AI envelopes associated with flank strata have a significantly longer "tail" to the northwest of the buildup core. These envelopes extend as long as ~30 m to the northwest, while the southern, western, and eastern envelopes are roughly similar in length (less than ~10 m-long). In the 3D GPR data (Fig. 9C and D), the downlapping flank strata are also asymmetrical and significantly elongated towards the northwest.

Statistical description and clustering of the buildup complex
We present a compiled statistical summary of RL, RS, R3, and Azimuth for S1, S2, and S3 scale buildups across all three high-frequency cycles within the stromatoporoid/coral complex of HCS-2 as boxplots ( Fig. 10A-E). We observe an exponential increase in RL, RS, and R3 when S1 buildups evolve into S2 clusters and S3-scale buildup complexes. We observe a log-linear relationship in the areal expansion ( Fig. 10E) from S1 to S2 and S3. We present the aspect ratio (RL/RS) and the Azimuth of S1, S2, and S3 scale buildups as boxplot and roseplot for data visualization (Fig. 10C, F, 10G, and 10H). The S1 buildups are the most circular of all three buildup scales, as their range of aspect ratio is between 1.02 and 1.7, with a mean of 1.3 and a standard deviation of 0.2. The S2 buildup clusters are more elongated than S1 in shape. The range in aspect ratios of S2 clusters is between 1.2 and 2, with a mean of 1.5 and a standard deviation of 0.4. The ellipsoids of the S3 buildup complexes are the most elongated of all three buildup scales. They have a range of aspect ratios between 1.3 and 2.5, a mean of 1.9, and a standard deviation of 0.3. The roseplot of S1 Azimuth (Fig. 10F) shows a uniform distribution with no preferential orientation. Roseplots of S2 and S3 scale buildups (Fig. 10G and H) reveal a strong NW-SE orientation (~N33W), which is exceptionally well developed for buildups with an aspect ratio >~1.8. ~80% of S2 and S3 buildups are elongated and oriented between N45W to N20W (northwest-southeast orientation).
We utilize Ripley's K-function to investigate the clustering tendencies of the buildups. Ripley's K-function of S1 buildups shows a statistically significant clustering below ~60 m (Fig. 10I). The function line of observed-K between ~60 m and ~120 m is below the expected K, indicating that S1 starts being dispersed at this range. However, the degree of dispersion is relatively weak compared to the significant dispersion at a distance larger than 180 m. Therefore, a weak clustering can still be inferred for S1 buildups up to ~140 m. Ripley's K-function of S2 buildup clusters (Fig. 10J) shows a statistically significant clustering up to ~150 m. However, unlike S1 buildups that exhibit weak clustering at a significant distance beyond the intersection, S2 buildups immediately dispersed beyond ~150 m. Ripley's K-function of S3 buildup complexes (Fig. 10K) shows a statistically significant clustering below ~250 m. S3 buildups immediately dispersed with no significant clustering at the intersection point. Buildup mapping and statistical observations have established a hierarchical scaling relationship where S1 buildups congregate into S2 clusters, and S2 clusters grow into S3 complexes. Ripley's K function of each buildup scale confirms the long axis of the following buildup scale in the hierarchy. S1 buildups (mean long axis of ~8 m) have a clustering distance of ~60 m, corresponding to the maximum long axis of S2 buildups at ~68 m. S2 buildups have a clustering distance of ~150 m, corresponding to the maximum long axis of S3 buildups (~158 m).

Morphometric description of the S2 buildups
Based on their radargram facies characteristics, the size of our dataset, and the definitive clustering discrepancy presented by Ripley's K analysis, S2 buildup clusters are the easiest buildup scale to map and differentiate in 3D. Therefore, we perform a more detailed quantitative analysis of this buildup scale (Fig. 11A-C). In the study area, we mapped 72 S2 buildups in cycle-1, 49 in cycle-2, and 19 in cycle-3. Histograms of RL, RS, and R3 show a lognormal distribution across the three cycles. However, their magnitude decreases from cycle-1 to cycle-3. For example, the mean RL for cycle-1 is ~40 m, while the mean RL for cycle-2 and cycle-3 is ~35 m and ~19 m, respectively. We observe a positive linear correlation with a strong correlation coefficient (~0.87) between RL and RS (Fig. 11D). On the other hand, RL and RS have a moderate positive correlation with R3 ( Fig. 11E and F).

Statistical description of the S2 flank strata
Histograms of flank distributions for all three high-frequency cycles follow lognormal distribution (Fig. 12C-E). Roseplots of flank Azimuth color-coded by their length (Fig. 12F-H) reveal that flank distribution in cycle-1 shows a significant asymmetric elongation towards the northwest (~N22W) with some minor elongation towards the west (~N112W). This asymmetric NW-SE elongation of the flanks is also observed in cycle-2 (~N35W). However, a considerable percentage of flank facies deviate from this orientation with elongations towards the When sorted by their hierarchical scale relationships (S1, S2, and S3), the buildups show an exponential or logarithmic parameter increase when one buildup scale forms the next buildup scale in the hierarchy. (10F) The distribution of Azimuth color-coded by aspect ratio shows that S1 buildups are circular, and their Azimuth is uniformly-distributed. (10G and 10H) Distribution of Azimuth color-coded by aspect ratio shows that S2 and S3 buildups are more elongated (more ellipsoid) and more strongly oriented towards northwest-southeast. Ripley's K function of S1 (10I), S2 (10J), and S3 (10K) show the distance at which the buildups at each scale tend to be dispersed or clustered. The resulting K-functions also confirm the long axis of the next buildup scale in the hierarchy. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.) northeast, southeast, and southwest. In cycle-3, a consistent flank orientation is not observed anymore (Fig. 12H).

Environmental condition during stromatoporoid/coral facies initiation
A negative shift in gamma-ray value distinctly characterizes the first appearance of the stromatoporoid/coral facies in the study area. This clean gamma-ray signature indicates a more distal depositional environment in this extensive lagoon, free of argillaceous input and well-lit with a well-functioning carbonate factory (Ehrenberg and Svånå, 2001;Al-Mojel et al., 2020). The low gamma-ray signature is maintained throughout the rest of HCS-2 deposition, indicating that the stromatoporoid/coral buildup is built during the late transgressive to early highstand system tracts of HCS-2 when argillaceous influx is relatively absent. The association of stromatoporoid/coral facies in our study area with diverse lagoonal benthic foraminifera, echinoderm, and lack of dasycladacean algae domination also suggests seawater was of near-normal salinity (Leinfelder et al., 2005).

Stromatoporoid/coral buildup organization and the role of hydrodynamic energy
The stromatoporoid/coral buildups exhibit scaling relationships. The areal expansion of S1 to S2 and S3 follows a linear logarithmic trend while vertical expansion follows an exponential increase (Fig. 10). Detailed histograms of S2-scale buildups sorted by their high-frequency cycles show a consistent lognormal distribution for all parameters (Fig. 11). This exponential growth from small buildups to larger complexes and their logarithmic frequency distribution is consistent with studies by Purkis et al. (2007) and Purkis (2018), where reef size is best described by power law, logarithmic, or negative exponential relationship. Ripley's K analysis shows that clustering distance from each buildup scale confirms the long axis of the next buildup scale in the hierarchy. S1 scale buildups remain relatively clustered up to the clustering distance of ~140 m, almost similar to the mean long axis of S3 complexes (~128 m). These clustering consistencies across multiple scales strongly indicate the self-organized nature of the stromatoporoid/coral buildups in our study area. This spatially self-organized scaling pattern is found in both modern and ancient reef builders, as reported in Purkis et al. (2015), Rameil et al. (2010), and Khanna et al. (2020). S1 buildups are more circular and have no preferential orientation. However, once they cluster into a larger scale (S2), the buildups evolve into more anisotropic and preferentially oriented towards northwestsoutheast. This preferential orientation is even more prominent for S3scale complexes as they are the most anisotropic (average aspect ratio of ~2), and most of them are oriented between N45W to N20W. During the deposition of HCS-2, the study area was located 10 • -15 • south of the equator with a northwest-southeast paleo-trade wind direction (Fig. 1B). The present-day Arabian Plate has rotated approximately 15 • clockwise (Ziegler, 2001;Scotese, 2003) from where it was during HCS-2 deposition. Rotating the Azimuth of S2 and S3 buildups 15 • degree counter-clockwise results in the consistency between the direction of the paleo-trade wind and the buildups orientation (Fig. 15B). Therefore, we interpret that currents oriented parallel to the trade wind direction exerted a dominant control over the orientation and eccentricity of the larger-scale buildup clusters and complexes. Once buildups were initiated, their growing topography and clustering started to restrict current-driven water flow into the remaining inter-buildup space leading to the amplification of current speed. This amplification caused the buildups to develop elongated shapes with an orientation parallel to the currents.
The development of flank facies with a preferred orientation approximately parallel to the proposed NW-SE current direction further supports this interpretation. The asymmetric elongation of flanks to the NW indicates bi-directional flow with a stronger component towards the NW. The minor orientation components of flanks in the southern, eastern, and western directions away from the buildups are attributed to erosion by tidal currents. Additionally, we associate the division between the mud-supported (FC4B) and grain-supported (FC4C) flank strata with the influence of hydrodynamic energy. The stromatoporoid/ coral floatstone with wackestone-packstone matrix (FC4B) is associated with sheltered flanks, while stromatoporoid/coral floatstone with grainstone matrix (FC4C) indicates agitated and current exposed flanks. The role of local currents and tidal influence in reef detritus distribution is also reported by Khanna et al. (2020) in the Cambrian microbial buildup in Texas. Hence, currents generated by trade winds and tidal forces may have contributed to shaping sediment distribution and maintaining near-normal marine conditions in this Late Jurassic giant epicontinental lagoon.
Based on the three-dimensional morphology, clustering, and scaling relationships of the stromatoporoid/coral buildups, we propose the following model (Fig. 13) in order to explain the lateral and vertical growth of the stromatoporoid/coral buildup complex in the three stromatoporoids dominated subcycles of HCS-2: • Initiation stage (Fig. 13D, H, and 13I) The stromatoporoid/coral buildups initiate at the nucleation points. The buildup cores are primarily small and circular, with no preferential orientation at this initial stage. Due to their small size and lack of resulting current restriction, flank strata do not develop at this stage.
• Congregation stage ( Fig. 13C and G) The buildup cores that grow close enough within the clustering distance (e.g., ~60 m for S1-scale buildups) assemble into larger buildup clusters. The dominant hydrodynamic energy starts to influence and be influenced by the growth and assembly of these clusters, creating a selfimposing restriction. The buildup clusters become elongated and oriented parallel to the paleo-tradewind direction. The flank strata at this • Maximum growth stage (Fig. 13B, F, 13J) With increasing accommodation, the congregated buildup clusters form even larger complexes from multiple initiation points. The buildup complexes reach maximum volume expansion at this stage. Similarly, the asymmetric development of muddy floatstone facies as flank strata achieve their maximum extent at this stage.
• Termination stage (Fig. 13A, E, and 13K) The accommodation is limited towards the end of a cycle. The vertical growth of the buildup complex is constrained. Hence, the buildup complex expands more laterally. Consequently, the buildup complex at this stage becomes flat-topped. Due to the relatively shallowing environment, the inter-buildup strata switch from predominantly mudsupported facies into skeletal peloidal grainstones. Hence, the flank strata at this stage are primarily stromatoporoid/coral floatstone facies with a grain-supported matrix.

Modeling method for the stromatoporoid/coral complex in HCS-2 based on quantitative observation
This study presents a comprehensive quantitative observation of the buildups and their associated flank strata. This information is utilized to perform simple modeling that attempts to build a realistic geocellular outcrop model of the upper part of HCS-2, which honors the complex lateral heterogeneity of the stromatoporoid/coral facies observed in this study (Fig. 14). The modeling methodology described here is analogous to the surface-based method by Janson and Madriz (2012), where the buildups are constructed by first recreating the overall three-dimensional morphology as surfaces. The model is constructed based on the S2-scale buildups using ArcGIS™ and JewelSuite™ software. The model's dimension is 1.2 km by 1 km by 25 m, representing the areal extent of the study area and the thickness of three stromatoporoid/coral-dominated cycles in HCS-2. The model's grid size is 4 m × 4 m x 25 cm.
Step-by-step modeling methods are explained below.
Step 1. Four stratigraphic surfaces (or grids) representing the base of the stromatoporoid/coral complex, top of cycle-1, top of cycle-2, and top of Hanifa (top of cycle-3) are utilized as inputs for stratigraphic zonation. The base of the stromatoporoid/coral complex is inferred from the measured sections at the elevation of the first significant negative shift in gamma-ray value. The top of cycle-1, cycle-2, and cycle-3 are imported from drone and geophysics interpretation. After a structural correction of ~1 • to the east, these surfaces are flat on a 1.2 km by 1 km scale (Fig. 14A).
Step 2. Random buildup nucleation points are generated on the base of the stromatoporoid/coral complex surface (Fig. 14B). A conditional statement is enforced to limit the spread of nucleation points. The maximum number of points is dictated by the number of mapped buildups in the area for each cycle. For example, the maximum points in cycle-1 are 288 buildups per 1.2 km 2 . The minimum distance between points is taken from Ripley's K analysis (i.e., 60 m for S1 to S2 clusters, Fig. 10I).
Step 3. A value of RL and Azimuth is assigned for each nucleation point. These values are sampled from normalized RL distribution (i.e., Fig. 11A) and Azimuth roseplot (i.e., Fig. 10G). Since RL and RS have a strong positive correlation coefficient (0.87), we assign RS from RL by an empirical equation (i.e., Fig. 11D). A 2D ellipse with a major axis, minor axis, and orientation based on RL, RS, and Azimuth is generated from each nucleation point. Overlapping ellipses are merged into a more extensive complex (Fig. 14C).
Step 4. A value of R3 is assigned for each nucleation point. Since R3 moderately correlates with RL (0.54) and RS (0.63), this parameter is determined by an empirical equation derived from multiple linear regression using RL and RS as variables. The equation is R3 = 0.98 + 0.07*RL + 0.06*RS -0.0004*RL*RS. Using RL, RS, R3, and Azimuth, a 3D ellipsoid representing the buildup for each nucleation point is constructed ( Fig. 14D) at this stratigraphic level.
Step 5. Each buildup ellipsoid is expanded using the value of FN, FS, FE-FW, and their Azimuth sampled from normalized flank distribution (i.e., Fig. 12). The expansion value is proportional to the size of the buildup ellipsoid. The resulting surface (Fig. 14E) represents flank strata associated with the buildup ellipsoid. The overlapping flanks are merged into a larger complex.
Step 6. A surface from the elevation of R3 is constructed and serves as the base for the subsequent buildup nucleation points (Fig. 14F). A density map is generated based on the previous nucleation points. The unit of this density map is normalized into a probability of S2 clusters growing into an S3 complex. For example, 118 out of 140 S2 clusters mapped in the study area grow into S3 complexes. Therefore, the maximum probability of nucleation points in level n+1 to grow on top of the nucleation point n is 0.8 (Fig. 14G). A similar conditional statement (buildup density and minimum cluster distance) is implemented to limit the spread of nucleation points in the next level (n+1).
Step 7. Steps 3 to 6 are repeated, filling the available accommodation space within one cycle with modeled buildups. Buildups that grow beyond the top of the cycle boundary are truncated ( Fig. 14H and I). Afterward, steps 1 to 7 are repeated for a new cycle using a different set of nucleation points, different distributions, and parameters representing their cycles. The steps are repeated until the entire HCS-2 (up to the top of cycle-3) are populated with buildup ellipsoid and their flanks strata.
Step 8. Facies FC4A is assigned for all buildup ellipsoids in all cycles. The inter-buildup strata (background) is filled with facies proportion identified from measured sections. The inter-buildup strata in cycle-1 are FC2. The inter-buildup strata in cycle-2 are mostly FC2, with the occurrence of FC6 on the upper part. The background strata in cycle-3 are intercalation between FC5 and FC6. Flank strata with FC2 background are assigned as FC4B, while flank strata with FC6 background are FC4C (Fig. 15A).

Subsurface implications related to the heterogeneity of the stromatoporoid/coral buildup facies
A detailed 3D statistical investigation of the stromatoporoid/coral buildup and its' associated flank in the study area allows for a quantitative description and characterization of the 3D morphological architecture of these facies. The collective influence of these facies on development decisions in an analog reservoir requires an integrated reservoir modeling approach that includes dynamic reservoir simulation, which would be an individual study by itself. Instead, this study briefly analyzes the impact of the heterogeneity of the stromatoporoid/ coral facies in an analog subsurface reservoir through computing the static reservoir connectivity. We construct 3D outcrop-based porosity and permeability models from the depositional facies model (Fig. 15A) and populate it with porosity and permeability from subsurface reservoirs (Table 3)   1993; Lucia et al., 2001;Al-Khalifah and Makkawi, 2002), assuming depositional facies controls the distribution of petrophysical properties (Pringle et al., 2006). These models ( Fig. 15C and D) serve as a visual example of property heterogeneities within the stromatoporoid/coral complex of the HCS-2. The patchiness of high porosity and permeability stromatoporoid/coral buildups and their flank strata within an otherwise lower porosity-permeability background is visible (Fig. 15C and D). Thus, connectivity between the stromatoporoid/coral facies might play an essential part in establishing the optimum field development concept in this type of reservoir.
To briefly assess the connectivity between the buildup complexes in the context of the reservoir without running a dynamic simulation, we perform a static connectivity exercise. This process calculates the number and volume of connected cell clusters corresponding to one or more facies codes. We found over 250 buildup cell clusters within the 1.2 km × 1 km x 28 m model size. The Gross Rock Volume (GRV) of individual cell clusters' varies from ~1% to ~8% of the total GRV of the stromatoporoid/coral complex. This high number of clusters indicates the lack of coalescence and quantifies the extreme 3D spatial heterogeneity of the buildups. The top 20 largest buildup clusters present in the model contain ~30% of the total GRV (Fig. 16A). The exercise continued by placing five vertical wells into the model with ~1 km spacing as an buildups from cycle-1 show a preferential northwestsoutheast orientation consistent with the direction of paleo trade wind (after 15 • rotation). This consistency indicates a dominant influence of trade wind associated current to the lateral expansion of the buildups. 3D porosity (15C) and permeability (15D) models show lateral property heterogeneities of the stromatoporoid/coral facies.

Table 3
Subsurface equivalent porosity and permeability of the modeled facies (Kompanik et al., 1993;Lucia et al., 2001;Al-Khalifah and Makkawi, 2002 analog to the typical state of Middle Eastern field developments (Vahrenkamp et al., 2019). We investigate the extent of connected stromatoporoid/coral facies volume through these five vertical wells. The "best" scenario is when all five vertical wells at ~1 km spacing drill through clusters of connected stromatoporoid/coral facies. However, this "best" scenario would only access ~14% of the total GRV occupied by the FC4A facies (Fig. 16B). The "worst" scenario is when all five wells do not encounter stromatoporoid/coral facies. Another "bad" scenario is when 3 out of 5 wells encounter these facies, yet due to the extreme lateral heterogeneity, these wells would only access a mere ~2% of the total GRV available (Fig. 16C).
One of the key uncertainties in the connectivity of buildup clusters is the influence of the flank strata onlapping against the buildup cores. The flank facies are either stromatoporoid/coral floatstones within a mudsupported matrix (FC4B) or with a grain-supported matrix (FC4C). The equivalent subsurface porosity/permeability of the FC4C approaches that of the buildup cores (FC4A), making it property-wise almost an extension of the buildups. Thus, we performed static connectivity analysis, assuming that both depositional facies form bodies of an extended reef complex with similar properties. In this case, we still detect more than 250 clusters. However, the grainy flank strata connect several major buildup clusters and form an extended "finger-like" buildup mega cluster elongated northwest-southeast (Fig. 16D). This mega cluster contains ~37% of the total GRV present in the extended reef complex (Fig. 16D) and is the most significant volumetric "pay zone" in the area. The "best" case scenario using five vertical wells The "best" case scenario of five vertical wells at ~1 km spacing can only access ~12% of the total stromatoporoid/coral facies' GRV. (16C) The second "worst" five vertical wells scenario at ~1 km spacing can only access 2% of the total GRV of the stromatoporoid/coral facies. (16D) The top 20 most significant extended buildup clusters assuming stromatoporoid/coral buildup (FC4A) and grainy flank strata (FC4C) are connected. Notice the northwest-southeast elongated mega cluster outlined by the dashed line. (16E) The "best" case scenario of five vertical wells in the connected buildup-grainy flank setting can access 43% of the total GRV. (16F) The second "worst" five vertical wells scenario in the connected buildupgrainy flank setting can only access 4% of the total GRV. (Fig. 16E) can access 43% of the available GRV present in the extended reef complex. 86% of the connected reefal GRV comes from the one mega cluster (light green in Fig. 16D). In this setting, a "worse" case scenario would only access 4% of the available GRV present in the extended reef complex (Fig. 16F). This order of magnitude difference emphasizes the importance of well planning and well design in a field development project to access the clustered and non-clustered highquality reservoir rock in an analog reservoir.

Conclusions
This study presents a quantitative investigation of meter-scale 3D morphology and spatial distribution of the stromatoporoid/coral buildups and their associated flank strata in a level of detail that was not previously known in the Hanifa Formation. This study provides a novel and valuable insight into the growth morphology of the stromatoporoid/ coral buildups and their relationship with environmental factors. The statistical findings compiled in this study can be directly utilized to model similar facies in the subsurface reservoirs stochastically. We also have shown that the static connectivity of analog subsurface reservoirs is directly correlated with the lateral heterogeneity of the stromatoporoid/ coral buildups and their flank strata.
The key findings of this study are: 1. Integrated drone-based digital outcrop model, measured sections, GPR, seismic, and core data can accurately map the 3D morphology and distribution of laterally heterogeneous stromatoporoid/coral buildups in the upper Hanifa Formation. 2. Stromatoporoid/coral buildups mapped in this study exhibit threelevel of hierarchical scaling relationships: S1, S2, and S3.
• S2-scale clusters are developed from S1 buildups. They exhibit pseudo-ellipsoid morphology (~35 m long axis and ~25 m short axis length) and are preferentially elongated towards northwestsoutheast. • S3-scale complexes grow from multiple S2 clusters. They are the most elongated (~128 m long axis and ~92 m short axis length) and oriented towards northwest-southeast. 3. S2 and S3-scale buildups are accompanied by asymmetric flank strata elongated to the northwest of the buildups. 4. Decreasing accommodation caused by buildup growth likely caused the focusing and strengthening of currents to the inter-build-up space, which self-restricts the orientation of the S2 and S3 buildup complex into a direction parallel to the dominant hydrodynamic regime. In the study area, northwest-southeast orientation is observed parallel to the paleo trade wind. 5. The development of asymmetric flanks is taken as evidence that trade winds and tidal currents impacted build-up orientation and flank facies deposition. The combination of wind and tidal currents likely caused the maintenance of near-normal marine conditions indicated by the biota in this vast epicontinental lagoon.

Declaration of competing interest
The authors declare the following financial interests/personal relationships which may be considered as potential competing interests:

Data availability
Data will be made available on request.