Three-Dimensional Modeling and Fluid Flow Simulation for the Quantitative Description of Permeability Anisotropy in Tidal Flat Carbonate

: Three-dimensional (3D) facies and petrophysical models were generated from previously published data of carbonate strata in the Dam Formation (eastern Saudi Arabia) to quantitatively investigate, describe, understand, model, and predict the permeability anisotropy in tidal ﬂat carbonate on the basis of a sequence stratigraphic framework. The resulting 3D models were used to conduct ﬂuid ﬂow simulations to demonstrate how permeability anisotropy inﬂuences the production of hydrocarbons and ultimately a ﬀ ects decisions concerning future drilling in the exploration and development of carbonate reservoirs with tidal ﬂat strata. The constructed 3D facies model consists of four lithofacies associations, two of which are grain-dominated associations and two of which are mud-dominated associations. These lithofacies associations vary spatially in four reservoir zones (zones 1 to 4), which represent two fourth-order sequences in the uppermost part of the Dam Formation. Zones 1 and 3 consist of transgressive parasequences, and zones 2 and 4 consist of the regressive parasequences of these sequences. The 3D porosity and permeability models have a coherent match with the distribution of the lithofacies and the stratigraphic framework of the Dam Formation. The results suggest that the permeability anisotropy in zones 1 and 3 is controlled by the occurrence of the grain-dominated lithofacies associated with tidal ﬂat channels. This lithofacies association overlies the sequence boundaries of sequences 1 and 3, forms reservoir bodies with relatively high permeability values, and is elongated perpendicular to the shoreline of the depositional environment. In contrast, permeability anisotropy in zones 2 and 4 is thought to be controlled by the occurrence of the grain-dominated lithofacies associated with the oolitic shoal. This lithofacies association overlies the maximum ﬂooding surface of sequences 2 and 4, forms reservoir bodies with relatively high permeability values, and is elongated parallel to the shoreline of the depositional environments. Fluid ﬂow simulation results suggest that the trend in hydrocarbon production from the constructed 3D models depends on permeability anisotropy in each reservoir zone. Thus, recognizing trends in permeability anisotropy, which can be predicted using sequence stratigraphy, could help to identify potential areas for future drilling.


Introduction
In terms of the properties of porous media, anisotropy is defined as the variation in those properties based on direction [1]. Thus, such properties depend on the location and orientation in a porous medium [2,3]. Isotropy, the opposite term, implies that the physical properties of the porous medium are not directionally dependent [2,3]. The permeability (the ease with which a fluid flows through rocks) of sedimentary strata is usually considered to be an anisotropic physical property because of the considerable vertical and horizontal variation in its values [3]. Recognizing and analyzing the controls on permeability anisotropy is very important to optimizing the fluid flow performance in hydrocarbon reservoirs and water aquifers and to locating sites for the storage of carbon dioxide.
In tidal flat carbonate, permeability anisotropy results from the interfingering of the heterogeneous facies mosaic [4][5][6][7][8][9][10]. Therefore, distinguishing the spatial and temporal patterns of sedimentary strata in the rock record provides insights into the controls on the permeability anisotropy in subsurface reservoirs [11][12][13][14][15][16]. Previous works in modern tidal flat environments and in their equivalent strata in the rock record have recognized various geological parameters that can be used to develop predictive facies models that facilitate the understanding of permeability anisotropy in these environments [12,[17][18][19]. Most of the previous works were centered on qualitative descriptions, such as sedimentological, stratigraphic, and diagenetic characteristics, but quantitative descriptions have received less attention.
In this context, the objective of this study is to investigate illustrative outcrops of tidal flat carbonate (Dam Formation) in the Lidam area of eastern Saudi Arabia to enhance the quantitative descriptions and to understand the controls on permeability anisotropy in carbonate reservoirs of tidal flat settings. Through generating 3D facies and petrophysical models and running fluid flow simulations on these models, this effort recognized trends in the spatial variability in sedimentary strata, permeability, porosity, and hydrocarbon production within a sequence stratigraphic framework. Therefore, the results provide perceptions that facilitate the understanding of the controls on permeability anisotropy that influence fluid flow behaviors in subsurface reservoirs with tidal flat carbonate.

Study Area and Geological Background
The Early Miocene (Burdigalian) Dam Formation cropped out in the Al-Lidam area (26 • 21 42" N, 49 • 27 42" E) approximately 80 km west of Dammam in the Eastern Province of Saudi Arabia (Figure 1). The formation was described first by Steineke and Koch [20] as strata that consist of carbonate rocks, and based on the lithological characteristics, fossil contents, and stratigraphic pattern, these strata were interpreted to have been deposited in tidal flat settings [21][22][23][24][25][26]. The strata of the Dam Formation are exposed in the study area as a series of connected and disconnected escarpments with a general trend in the NNW-SSE direction ( Figure 1C). Stratigraphically, the Dam Formation (which is~100 m thick) overlies the Hadruk Formation and is overlain by Hofuf Formation the ( Figure 1B) [21,27]. In the study area, only the uppermost part of the Dam Formation is exposed ( Figure 1C, [21,22]). Recent studies [23,24] have divided the exposed strata of the Dam Formation in the Al-Lidam area into several high-frequency sequences (see Section 3 for more details).  [21]); (C) geological map of the Lidam area showing outcrops of the uppermost part of the Dam Formation. Note that the outcrops with numbers in blue boxes trend~NW-SE parallel to the presumable paleoshoreline (studied by Ali [24]). The outcrops with numbers in black circles trend approximately E-W, perpendicular to the presumable paleoshoreline (studied by Bashri et al. [25]). The red box indicates the area of interest (AOI) where the 3D modeling was performed.

Location
The input data for the geostatistical models in this study came from previous works [23][24][25][26]. These works investigated 26 outcrops in the study area ( Figure 1C), and 14 outcrops were selected to be measured in detail. These 14 outcrops represent two transects ( Figure 1C). The first transect ( Figure 1C) includes six measured sections and represents a trend that runs approximately E-W, presumably perpendicular to the paleoshoreline of the depositional environment [25]. The second transect ( Figure 1C) includes eight measured sections and represents a trend that runs approximately NW-SE, presumably parallel to the paleoshoreline of the depositional environment [24]. For the 3D geostatistical modeling, we selected an area encompassing eight of these studied outcrops ( Figure 1C). The location that was selected was found to have one side approximately 2.40 km long in the N-S direction and another side approximately 3.64 km long in the E-W direction ( Figure 1C). This location was enclosed by a rectangular polygon that had an area of 8.73 km 2 . In the subsequent discussion, we will refer to this polygon as the area of interest (AOI).

Stratigraphy
Data from the measured sections in the study area include details about stratigraphic units and the sequence stratigraphic framework, which were presented as 2D cross-sections ( Figure 2) by Ali [24] and Bashri et al. [25]. Stratigraphic surfaces, such as sequence boundaries and maximum flooding surfaces ( Figure 2), were used in this study to construct surfaces and constrain the zones in the 3D geostatistical models. Parasequences, bed sets, and bed thickness were used to define the layers in the 3D geostatistical models.

Lithology
The data from the measured sections in the study area include detailed lithological descriptions (Table 1). Although there were slight differences in the lithological descriptions provided by Ali [24] and Bashri et al. [25], particularly in the naming and grouping of the facies, the differences were not substantial, and the lithological information on the study area from these two studies can be unified. Based on grain characteristics (type, texture and size), sedimentary structure, fossil content, and stratigraphic pattern, Ali [24] identified 17 lithofacies, whereas Bashri et al. [25] identified 15 lithofacies (Table 1). These lithofacies were grouped in lithofacies associations and used as input for modeling the facies (Table 1; see Section 4.2).

Petrophysics
Data from the measured sections in the study area include the results of 105 paired porosity and permeability measurements conducted on core plugs with diameters of 1.5 inches [26]. The porosity and permeability data were integrated with the lithofacies data and used to generate 3D porosity and permeability models, which were used for mapping permeability anisotropy and running the flow simulation conducted in this study.

Conceptual Model
The data from the measured sections include conceptual models for the Dam Formation in the study area [23][24][25][26]. The conceptual model demonstrated that the Dam Formation in the Lidam area represents facies belts of tidal flat lithofacies associations crossed with tidal channels. The paleoshoreline of this conceptual model was interpreted to trend approximately NW-SE. These conceptual models were used to verify the final geostatistical model by comparing the output of the 3D facies model and the conceptual model in terms of the distribution of lithofacies belts.

Workflow
The workflow of this study comprised five steps, all of which were conducted in Petrel™ 2019.
Step (1): generating a 3D structural grid for the AOI to serve as a structural framework for the 3D facies and petrophysical models.

Generating the 3D Structure Grid
First, we constructed a 3D structural grid for the AOI (approximately 2.4 km by 3.64 km, i.e., an area of approximately 8.73 km 2 ) ( Figure 3A). The 3D grid was generated using a simple grid method since no faults were found in the AOI. The location of the 3D grid was constrained by four geographic points, i.e., two latitude points (i.e., 25 • 13'55.65 and 25 • 15'15.20) and two longitude points (i.e., 49 • 29' 52.15 and 49 • 32'4.66) ( Figure 3A). In the 2D horizontal dimension, the 3D grid had 86,400 cells, i.e., 240 cells on the side with a north-south orientation and 360 cells on the side with an east-west orientation ( Figure 3B). Each cell had dimensions of 10 × 10 m ( Figure 3B). These dimensions were designed to capture the lateral lithological heterogeneity of the tidal flat setting, and they were intended to generate a 3D facies model that is significantly finer than the subsurface models, which usually have cells that have dimensions of 250 × 250 m. We used a stratigraphic framework from previous studies to generate surfaces within the 3D grid that was constructed ( Figure 3C,D). Five stratigraphic intervals were defined from previous studies, and they were used to generate five surfaces in the 3D grid ( Figure 3D). These stratigraphic intervals are as follows (from bottom to top): (1) the first sequence boundary (SB 1); (2) the first maximum flooding surface (MFS 1); (3) the second sequence boundary (SB 2); (4) the second maximum flooding surface (MFS 2); and (5) the top of the outcrops (which we considered as the third sequence boundary [SB 3]). These five surfaces defined the four reservoir zones (zones 1-4) in the constructed model ( Figure 3C,D).
The layers in the reservoir zones that were constructed ( Figure 3E) determine the thickness of the cells in the 3D grid. These layers are considered the main architectural elements in each zone, and they must adequately replicate the architecture of the stratigraphic units in the outcrops. Thus, layers within the constructed reservoir zones were assigned based on lithofacies thickness. The thickness of the layer in the 3D model was 0.2 m, and there were 65 layers in the model. After constructing these 65 layers, the total number of cells in the 3D grid was 360 × 240 × 65 cells, i.e., 5,616,000 cells.

Generating the 3D Facies Model
In 3D facies modeling, we populated the lithofacies that were defined by previous studies as discrete data into the 5,616,000 cells of the 3D grid. This step consisted of five tasks, i.e., (1)    The first task is to group the lithofacies of the studied outcrop into genetically related lithofacies associations (Table 1, Figure 5A,B). The suggested lithofacies associations were as follows (from offshore to onshore, Figure 4): (1) the grain-dominated skeletal oolitic shoal; (2) the mud-dominated skeletal oolitic shoal fringe; (3) the mud-dominated tidal flat, and (4) the grain-dominated tidal channel.
The second task is to digitize the lithofacies of the measured sections and insert them as facies logs in pseudo-wells that represent the locations of the outcrops ( Figure 5C,D). The facies data of the measured sections in the AOI were digitized manually. Each lithofacies association was assigned a facies code (e.g., 1, 2, 3). Each of these codes was given a certain stratigraphic interval based on the outcrop data. Thus, the facies log can be used as input data for 3D modeling.
The third task was to scale up the digitized facies log to the cells of the 3D grid ( Figure 5C,D). In this sense, each pseudo-well includes 65 layers of different facies. The upscaled facies data ( Figure 6D) represent the hard data of the 3D model at the well location and were used for data analysis ( Figure 5E-G). The upscaled data were used to perform quality control of the modeling results in the sense that the resulting 3D facies model should honor the upscaled data. The fourth task was to perform the data analysis ( Figure 5E-G). The functionality of the analysis of the statistical data is determined by the discrete (indicator) variogram analysis ( Figure 5E), the variogram map ( Figure 5F), the proportion of the vertical facies ( Figure 5G), and the thickness of the facies. This functionality allowed us to check the quality of the input of data to understand statistical variations in the lithofacies. The results of the data analysis and the conceptual sedimentological model were used to construct a geologically sensible 3D facies model.
In the fifth task, we generated 2D trend maps (horizontal maps) using the proposed conceptual model for the Dam Formation in the study area ( Figure 5H,I). In these trend maps, the depositional environments of the lithofacies associations (shoal, shoal fringe, tidal flat, and tidal channel) were shaped to present the spatial arrangement of their geographic trends ( Figure 5H,I). These trend maps determined the spatial probability of each lithofacies association.
In the final task of this step (task 6, Figure 5J), the digitized lithofacies in the pseudo-wells were distributed in the cells of the models using sequential indicator simulation (SIS). The SIS function was constrained by the parameters of the variograms that were constructed in the data analysis step and guided by the trend maps.

Generating the 3D Property Model
Two important properties, i.e., porosity and permeability, were modeled in this study ( Figure 6). The values of porosity and permeability were taken from Abdelkarim et al. [26] ( Figure 6C). Porosity and permeability were populated into the 5,616,000 cells of the 3D grid using a facies-based function in which the distribution of these properties was performed distinctly for each facies in the 3D facies model using the Gaussian random function simulation (GRFS) algorithm in Petrel™2019 ( Figure 6D,E). The data required to perform GRFS in Petrel™2019 were the range, minimum, maximum, mean, and standard deviation values of the simulated porosity and permeability. These values were used to constrain the random distribution of porosity and permeability in each distinct lithofacies association throughout the 3D grid ( Figure 6F,G). Average porosity and permeability maps for each zone were generated to detect the overall trend of these properties in the zones. These types of maps are useful because they filter out the noise that results from the unrealistic distribution of facies (algorithm errors). Variogram maps of the average maps were generated to understand the permeability anisotropy of each reservoir zone in the 3D model.

Generate Upscaled Models
In the upscaling process, we generated a coarser 3D grid model from the fine grid model, which had 5,616,000 cells. The upscaling process is essential to allow our current computer lab capacity to run the flow simulation, which involves applying complicated mathematical algorithms to the generated models. We created coarser grid cells with a cell dimension of 100 × 100 m 2 (a grid with 23,125 cells). This coarser grid reduced the total number of cells from approximately 5,616,000 to 23,125 cells. The 3D property models of the fine grid cells (i.e., the 5,616,000-cell grid) were transferred to these coarser grid cells (23,125 cells). Although the number of cells was reduced substantially in the upscaled grid, this reduction did not considerably affect the petrophysical distribution in the 3D models.

Fluid Flow Simulation
The fluid flow simulation was run using the black-oil simulator (Eclipse 100) in Eclipse-Petrel™2019. The simulation was run on the upscaled model, which was assumed to be a heavy oil reservoir with impermeable sides and a cap as boundary conditions. The pressure support was assumed to come from a water aquifer below the oil-water contact. The production of heavy oil from this hypothetical outcrop reservoir comes from vertical wells deployed to represent three transects.
The first transect (Wells A1-A3) and the second transect (Wells A1, A4, and A5) were designed to test the variability in oil production with permeability anisotropy in zone 1, along and across the tidal flat channel (first and second transects, respectively). Oil-water contacts for both transects were placed at the base of the tidal channel lithofacies. The perforation interval of the wells in both transects included the entire interval between the base of the channel and the upper boundary of zone 1 (MFS1). Thus, the production of oil was ensured to come from the tidal flat channels.
The third transect included four wells (B1-B4) and was designed to test the variability in oil production with permeability anisotropy in zone 2. The oil-water contact in Wells B1-B4 was placed at the base of zone 2 (MFS 1). Perforation of these four wells included the entire interval in this zone. Thus, the production of oil was ensured to come only from the lithofacies in zone 2.
Nine fluid flow simulation runs were performed with the production of oil coming from only one single well along the three transects in each simulation run (one production well at a time). Each well in these three transects had a single control model with a production rate of 10,000 STB/day. The initial water saturation of the hypothetical reservoir of the outcrop was 0.5, whereas the preliminary heavy oil formation-volume factor was 1.1. Water saturation and formation factors are dynamic parameters that continuously fluctuate with time, production, and changes in pressure. The software considers these dynamic changes and adjusts these parameters through time. It should be noted that the fluid flow simulation was run only for zones 1 and 2 with the notion that the results of zone 1 would be similar to the result in zone 3 and that the results of zone 2 would be similar to the result of zone 4 because of the similarity in lithofacies distribution and petrophysical data of these zones.

Analysis of the Data
The results of the data analysis ( Figure 5E-G) provide constraints on the vertical facies proportion, facies thickness, and facies probability of the 3D facies model. These constraints were computed from the upscaled lithofacies logs as vertical proportion curves and as percentages of the volume to the total model ( Figure 5E-G). The results of the data analysis also include experimental and modeled variogram parameters (nugget, sill, and major and minor ranges, Table 2) and variogram maps ( Figure 5F). The results of the data analysis identified two different directions of the major range of the variograms. The major ranges for the grain-dominated skeletal oolitic shoal, the mud-dominated skeletal oolitic shoal fringe, and the mud-dominated tidal flat were found to be oriented approximately NW-SE (Table 2). In contrast, the major range for the grain-dominated tidal channel was found to be oriented approximately NE-SW ( Table 2). The minor range of the constructed variograms for all of these lithofacies associations was markedly lower than their corresponding major ranges (Table 2).

Distribution of the Lithofacies
The results showed both lateral and vertical variations in the lithofacies of the 3D facies model (Figure 7; Table 3). Each reservoir zone in the 3D facies model had a distinct set of lithofacies with a proportion comparable to the proportion of lithofacies in the outcrops in the study area ( Figure 7A). Zones 1 and 3 had extensive tidal channels (Figure 8) that trended in the same direction as that suggested by the outcrop data (NE-SW). The average facies map of each zone in the 3D model (Figure 9) smoothed the results of facies modeling and provided the general trends of the lithofacies in each zone. Whereas zones 1 and 3 ( Figure 9A,C) showed a general NE-SW trend, zones 2 and 4 ( Figure 9B,C) showed a NW-SE trend. Descriptions of the lithofacies in each reservoir zone are provided in the following sections.    The grain-dominated tidal channel lithofacies association is clustered as a single channel volume that comprises~21.9% of the total volume of the lithofacies in the zone (Figure 8). The long axis of this channel volume trends~NE-SW ( Figure 9A), and it cuts across the other two lithofacies associations in the zone, which are oriented NW-SE ( Figure 9A).

Zone 2
Zone 2 (between MFS 1 at the base and SB 2 at the top) consists of two lithofacies associations ( Figure 7C; Table 3), i.e., (1) the mud-dominated skeletal oolitic shoal fringe (comprising the majority of the zone, 51.8%) and (2) the grain-dominated skeletal oolitic shoal (comprising 48.2% of the total volume of the zone). Unlike the zones above and below, this zone does not contain the grain-dominated tidal channel lithofacies (Figure 9). The grain-dominated skeletal oolitic shoal represents the downdip lithofacies for the mud-dominated skeletal oolitic shoal fringe ( Figure 9B), and in a few areas, the grain-dominated skeletal oolitic shoal is surrounded by the mud-dominated skeletal oolitic shoal fringe ( Figure 9B).

Zone 3
Zone 3 (between SB 2 at the base and MFS 2 at the top) consists of two lithofacies associations ( Figure 7A; Table 3): (1) grain-dominated tidal channel (comprising 18.6% of the total volume of the zone) and (2) mud-dominated skeletal oolitic shoal fringe (comprising the majority of the zone, 81.4%).
In this zone, the grain-dominated tidal channel trends NE-SW (Figures 8 and 9) and comprises~15.2% of the total volume of lithofacies in this zone. Although visually appearing as isolated bodies in the average map of the zone (Figure 9C), the tidal channels show high connectivity (Figure 8). Of the 15.2% of the total lithofacies volume composed of the tidal channel lithofacies, 12% is part of a connected volume (Figure 8).

Zone 4
Zone 4 (between MFS 2 at the base and the top of the outcrops) consists predominantly of grain-dominated oolitic shoal (which comprises 81.5% of the total volume of the zone) (Figure 7; Table 3). There are patches of mud-dominated skeletal oolitic shoal fringe (comprising a minority of the zone, 18.5%), which occur as elongated bodies and tend to be parallel to the shoreline ( Figure 9D). The grain-dominated tidal channel lithofacies is also absent from this zone (Figure 9).

Architecture of the Lithofacies
Previous studies on the Dam Formation have resulted in the construction of shoreline-controlled sequence stratigraphy with the suggestion that the paleoshoreline of the Dam Formation had a strike direction trending approximately NW-SE. The 3D model of the facies in this study showed lithofacies architectures that were comparable to those suggested by the sequence stratigraphic framework of the Dam Formation in previous studies ( Figure 7). Therefore, the internal architecture of the 3D facies model can be classified into transgressive and regressive parasequences.
Zones 1 and 2 comprise sequence 2 of the Dam Formation ( Figure 2). The basal part of zone 1 (Figure 7) is a mud-dominated tidal flat, and the grain-dominated tidal channel lithofacies associations, which occur as transgression-related facies, were deposited retrogradationally on top of SB 1 (Figure 7). The occurrence of these lithofacies associations in zone 1 is consistent with the early transgressive parasequences of Bashri et al. [25]. The grain-dominated oolitic shoal occupies the upper part of zone 1 ( Figure 7) and occurs extensively proximal to the proposed shoreline as late transgressive parasequences that are consistent with the Bashri et al. [25] sequence stratigraphic model. Zone 2, which overlies MFS 1 (Figure 7), is dominated by the extensive grain-dominated skeletal oolitic shoal, which reflects the suggested high-stand system tract of sequence 1. Remarkably, in zone 2, the grain-dominated tidal channel lithofacies is absent, suggesting shoreline transgression and retrogradation of the shoal lithofacies in the landward direction.
Zones 3 and 4 comprise sequence 3 of the Dam Formation ( Figure 2). The mud-dominated tidal channel lithofacies association with tidal flat channel bodies in zone 3 is interpreted to represent the transgressive parasequences of sequence 3 ( Figure 7). However, the extensive grain-dominated skeletal oolitic shoal of zone 4 is interpreted to represent the regressive part of sequence 3. This interpretation is consistent with the sequence stratigraphic framework of the Dam Formation proposed by Bashri [25]. The 3D facies model revealed some similarities in trends and distribution of the facies between zones 1 and 3 ( Figure 7B) and between zones 2 and 4 ( Figure 7C). This similarity can be linked to the above-defined sequence stratigraphic framework of the 3D model. Both zones 1 and 3 ( Figure 7B) represent the transgressive facies (TF) of sequences 2 and 3, respectively, and are referred to in the following discussion as zones with TF. Zones 2 and 4 ( Figure 7C), on the other hand, represent the regressive facies (RF) of these sequences, and these zones are referred to in the following discussion as zones with RF. The zones with TF have wedges that thicken to the SW and thin to the NE ( Figure 7B). In contrast, the zones with RF have wedge shapes with the opposite trend. In other words, they thicken to the NE and thin to the SW ( Figure 7C).

Trends in Petrophysical Properties
The 3D models of porosity and permeability have a coherent match with the facies distribution and architecture, showing consistent vertical and lateral correspondence ( Figure 10). Unsurprisingly, the grain-dominated lithofacies associations have better porosity and permeability values than the mud-dominated lithofacies association ( Figure 10). The porosity and permeability values of the zones with TF (zones 1 and 3) appear to be influenced by the presence of the grain-dominated tidal channel lithofacies association (Figure 10). The average maps of the porosity and permeability in zones with TF (zones 1 and 3) (Figure 11) showed relatively high values of porosity (ranging from approximately 15% to 20%) and permeability (ranging from approximately 100 to 200 mD) compared to the values of the background lithofacies (which range from 5% to 10% and from 2 to 20 mD for porosity and permeability, respectively). These relatively high values organize in trends that are consistent with the channel body trends (~NE-SW) ( Figure 11). Generally, the porosity and permeability are highest at the center of these channel bodies (~20 and~200 mD, respectively) and decrease to the flank of the channel body (~15 and 100 mD, respectively) ( Figure 11). The variogram map of permeability in zone 1 indicated an elliptical shape with NE-SW direction in the variance value and indicated that geometric anisotropy occurred perpendicular to the proposed shoreline of the depositional environment ( Figure 12A). In zone 3, this trend is not well defined ( Figure 12C).   In contrast to the zones with RF, the average porosity and permeability maps for zones with RF have areas trending NW-SE and arranged in a parallel manner ( Figure 13). These maps ( Figure 13) show a contrast between areas of a relatively high porosity and permeability (ranging from 10% to 15% and from 150 to 300 mD, respectively) and areas with a relatively low permeability (ranging from 2% to 5% and from 2 to 25 mD, respectively) ( Figure 13). Generally, the porosity and permeability values have intermediate values in the transition areas between these contrasting areas ( Figure 13). The variogram map of permeability in zones 2 and 4 (zones with RF) shows an elliptical shape in the NW-SE direction for the variance value, indicating that geometric anisotropy occurs parallel to the proposed shoreline of the depositional environment ( Figure 12B,D).

Trends in Fluid Flow Behavior
In the simulation runs, we used three indicators for oil production that changed through the specified time for the simulation runs (Figures 14-16), i.e., (1) oil in place (OIP); (2) cumulative oil production (COP), and (3) oil recovery efficiency (ORE). In each of the three well transects, these indicators varied based on the location of the well. At each location of a well (i.e., before the production of oil in the static model), the permeability in the reservoir zone varied, whereas the OIP for the entire model stayed the same. Thus, the variation in production indicators can be attributed to the lateral variation in permeability (i.e., permeability anisotropy). The simulation time of all of the wells in the three transects was set to 80 years ( Figure 13B). For example, at time zero in zone 1, the OIP in the model was the same for all of the wells, i.e., 32.85 × 10 6 STB, and started to decrease through time with the production of oil ( Figure 14).   . Note that these production indicators suggest that oil production increased to the SW. Note that these production indicators suggested that oil production increased in the direction of the basin (NE).

First Transect
The first transect (Wells A 1-3, Figure 14A) has an updip-downdip trend perpendicular to the interpreted shoreline of the depositional environment of the Dam Formation. This transect was designed to test the variation in the production of oil along the axis of the tidal flat channel. In the simulation, the following results were obtained after 20 years: • Well A 1, which was located in the updip direction ( Figure 14A), showed the highest ORE (0.2) ( Figure 14C) with a COP of 6.9 × 10 6 (STB) ( Figure 14B) and a decrease in OIP to 25.9 × 10 6 (STB) ( Figure 14A). • Well A 2, which fell between Wells A 1 and A 2 ( Figure 14A), had an ORE of 0.1 ( Figure 14C), a COP of 2.2 × 10 6 (STB) (Figure 14B), and a decrease in OIP to 25.9 × 10 6 (STB) ( Figure 14A). • Well A 1, which was located in the downdip direction ( Figure 14A), had the lowest ORE (0.0005) ( Figure 14C), a COP of 0.004 × 10 6 (STB) (Figure 14B), and an insignificant decrease in OIP tõ 32.79 × 10 6 (STB) ( Figure 13A).

Second Transect
The second transect (Wells A 1, A 4, and A 5, Figure 15A) had a trend approximately parallel to the interpreted shoreline of the depositional environment of the Dam Formation. This transect was designed to test the variation in oil production across the axis of the tidal flat channel in zone 1. In the simulation, the following results were obtained after 20 years: • Well A 1, which represented the central axis of the channel ( Figure 15A), had the highest ORE (0.2) ( Figure 15C) with a COP of 6.9 × 10 6 (STB) ( Figure 15B) and a decrease in OIP to 25.9 × 10 6 (STB) ( Figure 15A). • Wells A 4 and A 5, which represent production from the flanks of the channel, showed markedly lower production indicators, i.e., ORE values of 0.04 and 0.006, respectively; COP values of 1.5 × 10 6 (STB) and 0.21 × 10 6 (STB), respectively; and OIP values of 31.2 × 10 6 and 32.6 × 10 6 , respectively ( Figure 15).

Discussion
The results of the 3D geostatistical modeling made it possible to recognize trends and patterns in the spatial variability of the sedimentary strata and the petrophysical properties of the Dam Formation within a sequence stratigraphic framework. These results showed sedimentological, stratigraphic, and petrophysical patterns similar to the Dam Formation in the outcrops and similar to the modern analogs of tidal flat settings in different locations worldwide. The results of the fluid flow simulation suggested trends in the production of oil that were consistent with the trend in the anisotropy of permeability. Most importantly, these recognized trends and patterns in the results can be described quantitatively. Therefore, the resulting quantitative data can be utilized extensively as an analog to provide a robust assessment of permeability anisotropy and reservoir quality of strata deposited as tidal flat carbonate elsewhere.

Similarity of the 3D Model to the Modern Analog and Stratigraphic Record
Both the conceptual (Figure 4) and the 3D facies model (Figure 7) in this study showed a grain-dominated, skeletal, oolitic shoal surrounded by a mud-dominated skeletal oolitic shoal fringe that was cross-cut by tidal channels and flanked on the onshore side by the mud-dominated tidal flat lithofacies association. Interestingly, a similar arrangement of this pattern of lithofacies occurs in modern analogs [28][29][30]. For example, the documented lateral facies patterns observed in the Little Bahama Bank in the northern Bahamas show shoal facies surrounded by a mud-dominated background similar to what the 3D facies model has produced (Figures 4 and 7) [31][32][33]. This is especially significant given that the Bahamas have different carbonate physiography than the Dam Formation, i.e., an isolated platform [31][32][33] versus a carbonate ramp [34]. The 3D facies models also seem to be consistent with the modern facies pattern of the offshore UAE (Figure 7) in that the grain-dominated skeletal oolitic shoal is surrounded by a fringe of mud-dominated skeletal oolitic shoals [28].
The main difference between the stratigraphic record and the modern analogs of tidal flat settings is that the modern analogs represent one snapshot of the distribution of facies over time, whereas the stratigraphic record provides spatial and temporal perspectives for the distribution of facies [28]. Comparing the results of the 3D models ( Figure 7) with a stratigraphic record of a tidal flat in the Dam Formation suggests that the 3D models successfully reproduced the pattern of the tidal flat lithofacies association. This similarities between the 3D model and the modern analog and stratigraphic record of tidal flats suggests that we constructed a 3D model that can represent the setting of a tidal flat that ultimately can be used to investigate the controls on permeability anisotropy.

Understanding Permeability Anisotropy Using Variogram Parameters
The results of the 3D modeling suggested that the permeability anisotropy of the tidal flat settings in the study area can be understood from the variogram parameters. The results showed that the variograms of the permeability in the four zones have elliptical shapes with a major range substantially longer than the minor range, which suggests geometric anisotropy ( Figure 12).
Previous modeling results of the Dam Formation in the Lidam area [26] showed similar results; namely, the variograms of the petrophysical properties have a major range that is remarkably longer than the minor range. However, it is important to note that the major and minor ranges of the variograms determined in this study (Table 1), i.e., approximately 2 km and approximately 1 km, respectively, are considerably longer than the ranges of the variograms reported by Abdelkarim et al. [26], i.e., approximately 0.1 km and 0.03 km, respectively. This difference could be attributed to the limited area that was modeled by Abdelkarim et al. [26], who modeled an area of approximately 0.23 km by 0.06 km, i.e., approximately 0.014 km 2 , using measured sections that were closely spaced, which means that a full representation of the spatial extent of the lithofacies was not available. Another possibility could be that, in this study, we modeled lithofacies associations by grouping many individual lithofacies, whereas Abdelkarim et al. [26] separately modeled the variograms of these lithofacies. In either of these two cases, the pattern of the variograms suggests geometric anisotropy for the modeled lithofacies, which resulted in permeability anisotropy, which can be explained in a sequence stratigraphic framework and by sea-level change.

Understanding Permeability Anisotropy with a Sequence Stratigraphic Framework
The results suggested that changes in sea level could be the first-order control of permeability anisotropy in tidal flat settings. Zones with TF (zones 1 and 3) (Figure 7) represent deposition at the time of decreases in sea level. In such cases, permeability anisotropy is controlled by tidal channels (Figures 9-12). These channels usually flow perpendicular to the shoreline of marine depositional environments (Figures 9-12). In this case, high permeability areas concentrate around the channel belts, and low permeability areas concentrate away from the channel belts ( Figure 11). In the zones with RF (Figures 9 and 13), which represent deposition at the time of the sea level high stand, permeability anisotropy was controlled by location with respect to the oolitic shoal lithofacies associations, which were arranged parallel to the shoreline of the marine environments ( Figure 13). In this case, areas with high permeability are concentrated in areas in the downdip direction and around areas occupied by shoal deposits (Figure 13).

Predicting Permeability Anisotropy Using Production Data and Sequence Stratigraphy
The hydrocarbon reserve of the presumable oil reservoir in the constructed 3D model is a function of the total porosity in each zone (Figures 14-16). However, the oil production of this presumable reservoir is a function of the permeability in the 3D model (Figures 14-16). Since the production in each simulation run comes from a single well (Figures 14-16), the oil reserves would not be affected by the location of the well in the model. In contrast, oil production indicators (i.e., OIP, COP, and RF) through time (as functions of permeability) are affected significantly by the location of the well in the model because of the permeability anisotropy (Figures 14-16).
The results from fluid flow simulation showed that the oil production indicators (OIP, COP, and RF) of wells that penetrate zone 1 a higher closer to the tidal channel axis and in the updip direction and lower at the flank of the channel belt and in the downdip direction (Figures 14 and 15). The tidal flat channel lithofacies association on top of sequence boundaries usually occurred extensively in the updip direction and sporadically in the downdip direction from the shoreline of the depositional environment ( Figure 7). This lithofacies association also has a wedge shape with a thick interval in the updip direction and a thin interval in the downdip direction (Figure 7). The flanks of these channels are usually associated with mud-dominated lithofacies, resulting in poorer permeability in these areas than in the axis of the channel ( Figure 11). Thus, better reservoir quality rocks are most likely to be located at the axis of the tidal channel lithofacies in the updip direction in strata representing transgressive parasequences in a tidal flat setting (Figures 14 and 15).
In contrast to the oil production trend in zone 1, the oil production indicators (OIP, COP, and RF) of wells that penetrate zone 2 increase basinward, i.e., in the downdip direction ( Figure 16). The association of the oolitic shoal lithofacies on top of a maximum flooding surface (Figure 7) usually occurs extensively in the downdip direction, and it pinches out in the updip direction toward the shoreline of the depositional environment ( Figure 13). This lithofacies association also has a wedge shape with a thick interval in the downdip direction and a thin interval in the updip direction ( Figure 7). Thus, better quality reservoir rocks are most likely to be associated with oolitic shoal lithofacies in the downdip direction in strata representing regressive parasequences of the tidal flat setting (Figures 13 and 16).

Implications and Limitations
This study provides interesting implications when integrating sequence stratigraphy with facies modeling, petrophysical modeling, and fluid flow simulation to understand permeability anisotropy in a sequence stratigraphic framework. The results provide insights into several aspects of the heterogeneity of lithofacies, petrophysical variations and reservoir connectivity in tidal flat lithofacies. Ultimately, the results offer quantitative data on how permeability anisotropy varies in tidal flat settings and how the variations affect the production of hydrocarbons. This information is useful in predicting areas for future drilling in carbonate reservoirs with tidal flat lithofacies.
It should be noted, however, that this study did not take into account the impact of the diagenesis of the petrophysical properties of reservoirs. Rather, the study considered mainly the variation in the depositional texture of lithofacies. This could be one limitation when using the results of modeled outcrops to understand the equivalent subsurface reservoirs. It would be interesting if data on diagenesis in the Dam Formation from the subsurface were integrated with the results of this outcrop study.

Conclusions
The 3D outcrop modeling (facies and property modeling) and the fluid flow simulation of tidal flat strata in the Dam Formation (in eastern Saudi Arabia) resulted in the following key findings: 1.
In the 3D models, porosity and permeability vary vertically and laterally with the variation in lithofacies associations, suggesting depositional controls on permeability anisotropy of the studied tidal flat strata.

2.
Fluid flow simulation results suggest that hydrocarbon production is influenced significantly by permeability anisotropy. 3.
The 3D models and the results of the fluid flow simulation reveal patterns and trends in permeability anisotropy and hydrocarbon production that could be linked to the sequence stratigraphic framework.
Quantitatively, this study shows that the permeability anisotropy in carbonate reservoirs deposited in tidal flat settings is controlled significantly by trends in the depositional environment (shoal trend versus tidal channel trend). In the stratigraphic record, such trends are predictable in the sequence stratigraphic framework. Therefore, permeability anisotropy could also be predicted based on sequence stratigraphy.