Practical Aspects of Upscaling Geocellular Geological Models for Reservoir Fluid Flow Simulations: A Case Study in Integrating Geology, Geophysics, and Petroleum Engineering Multiscale Data from the Hunton Group

Optimal upscaling of a high-resolution static geologic model that reflects the flow performance of the reservoir is important for reasons such as time and calculation efficiency. In this study, we demonstrate that honoring reservoir heterogeneity is critical in predicting accurate production and reducing the time and cost of running reservoir flow simulations for the Hunton Group carbonate. We integrated three-dimensional (3D) seismic data, well logs, thin sections, outcrops, multiscale fracture studies, discrete fracture networks, and geostatistical methods to create a 100 × 150 × 1 ft gridded representative geologic model. We calibrated our gridded porosity and permeability parameters, including the evaluation of fractures, by history matching the simulated production rate and cumulative production volumes from a baseline fine-scale model generated from petrophysical and production data obtained from five wells. We subsequently reperformed the simulations using a suite of coarser grids to validate our property upscaling workflow. Compared to our baseline history matching, increasing the horizontal grid cell sizes (i.e., horizontal upscaling) by factors of 2, 4, 8, and 16 results in cumulative production errors ranging from +0.5% for two time (2×) coarser to +74.22% for 16× coarser. The errors associated with vertical upscaling were significantly less, i.e., ranging from +0.5% for 2× coarser to +10.8% for 16× coarser. We observed higher production history matching errors associated with natural fracture size. Results indicate that greater connectivity provided by natural fracture length has a larger effect on production compared to the higher permeability provided by larger apertures. We also estimated the trade-off between accuracy and run times using two methods: (a) using progressively larger grid cell sizes; (b) applying 1, 5, 10, and 20 parallel processes. Computation time reduction in both scenarios may be described by simple power law equations. Observations made from our case study and upscaling workflow may be applicable to other carbonate reservoirs.


Introduction
Time is money. Therefore, despite the availability of unlimited computer resources, a reservoir analyst may have a larger impact on field performance by quickly history matching dozens of wells Energies 2020, 13, 1604 2 of 27 using a more economical upscaled model (coarser model) than by slowly history matching a single well using a higher accuracy fine-scale model. Upscaling is the process of reducing a large number of the fine-scale grid blocks of a geological model to a smaller number of coarse-scale grid blocks while retaining the underlying geological properties [1][2][3][4]. There are copious upscaling methods that include arithmetic, harmonic, geometric, power law, pressure solver, and weighted averages algorithms, where the choice depends on the rock property [4].
The necessity of the upscaling process arises from the limitation on computing capacity and amount of available time required to process large geologic models. Although high cell resolution is preferred, such bulky models containing millions of cells are prohibitive for iterative production history matching, sensitivity analyses, production forecasting, reservoir optimization, uncertainty assessment, proposing new or infill drilling locations, and risk analysis using a probabilistic approach [5][6][7][8]. When a large number of discrete fractures is modeled, the model becomes even more cumbersome [9][10][11]. Therefore, upscaling is necessary in most cases. However, upscaling a fine-scale reservoir model must be performed with care. Excessive upscaling leads to loss of the detailed heterogeneities in the geologic model that directly affects the result accuracy [4].
The flow from thin bed conduits and thin bed baffles is particularly sensitive to inaccurate upscaling. Santacruz et al. [12] simulated the flow of a leveed-channel system, assigning different permeabilities to the thin slump deposits. They concluded that slump deposits acted as barriers to waterflood penetration into the oil zones, whereby the lower permeability of the slumps significantly reduced the production rate. In many reservoir simulation studies, upscaling is performed only in the vertical direction, as the areal grid cell size of the reservoir model typically approximates the heterogeneity provided by three-dimensional (3D) seismic data [13]. Furthermore, Ma et al. [13] found a 50% error in the simulated production volumes after upscaling. Figure 1 illustrates the change in production as a function of upscaling permeability. strations of permeability upscaling processes that lead to loss of the detailed heterogeneities in the geologic model that directly affects the rofile. The fine scale permeability model, which consists of more than 13 million cells, has cell sizes of 20 ft in X-direction, 30 ft in Y-direction irection. This model was upscaled by enlarging the grid size (reducing the number of cells to 2,530,134) and averaging the permeabilities. Sim for the excessively upscaled model resulting in only 87, 000 cells. Note the permeability heterogeneities along A-A' section are not the same, e al shape. The right figure shows the change in production rate over time for the two upscaled models. The fine scale model could not be handling 13 million cells is beyond their capacity. Notice the 30% increase in production rate after upscaling (Upscaled 1 to Upscaled 2).

Figure 1.
Illustrations of permeability upscaling processes that lead to loss of the detailed heterogeneities in the geologic model that directly affects the accuracy of production profile. The fine scale permeability model, which consists of more than 13 million cells, has cell sizes of 20 ft in X-direction, 30 ft in Y-direction, and 1 ft in vertical (Z) direction. This model was upscaled by enlarging the grid size (reducing the number of cells to 2,530,134) and averaging the permeabilities. Similar process was repeated for the excessively upscaled model resulting in only 87,000 cells. Note the permeability heterogeneities along A-A' section are not the same, especially in the black oval shape. The right figure shows the change in production rate over time for the two upscaled models. The fine scale model could not be run on our computers as handling 13 million cells is beyond their capacity. Notice the 30% increase in production rate after upscaling (Upscaled 1 to Upscaled 2).
Meddaugh [14] examined the impact of vertical and areal upscaling on fluid flow, using multiple geostatically generated porosity model realizations. He concluded that the changes in the cumulative water injection are more sensitive to the number of realizations than to the level of vertical upscaling, and that the recovery factor decreases when the areal upscaling increases. A challenge in upscaling geologic reservoir model parameters is to determine a sufficiently coarse reservoir grid cell size [2,3,8] without sacrificing the detail of the original structures, rock properties, layered nature of the reservoir, and reservoir heterogeneity. Regardless of the upscaling process used, a proper production history match using fluid flow simulation validates the upscaling process. The purpose of this study is, therefore, to find the optimal level of upscaling for the Hunton carbonate oil wells, beyond the production history match grid cell size, at which the most significant geological and petrophysical properties and the simulation accuracy is retained ( Figure 1). Moreover, this study demonstrates the relative error magnitudes arising from the numerical uncertainty versus reservoir heterogeneity that it is essential for modelers to understand the sources of error.
We begin our study with a review of the geology of the Hunton carbonate of Oklahoma as it appears in the three study areas. We then describe our methodology, based on the well-established workflow of defining and then history matching a finely gridded, high resolution model ( Figure 1). We construct a fine-grid model as our reference and use of it for simulating production from available Hunton wells. We then use a series of homogenous models to simulate production; thus, removing the effect of geologic heterogeneity while divulging the numerical errors associated with each subsequent increase in grid cell sizes. Next, we show the effect of upscaling both the horizontal and vertical grid cell properties on the computation run times and accuracy, constructing power law relations that can be used to predict the run times for future models. We conclude with a summary of our results and recommendations that the reservoir analysis requires a balance between the need for accuracy of a given well simulation to the value of having less accurate results for multiple wells in the oil field.
The significance of this work lies in: 1. demonstrating seamless integration of lithologic, natural fracture, petrophysical, 3D seismic, reservoir, and production data to obtain a comprehensive reservoir description; 2. dividing reservoir property upscaling errors (during reservoir simulation and history matching) into that associated with numerical instability and that associated with geologic heterogeneity in a quantifiable manner (i.e., through empirical equations) for different grid cell sizes; 3. quantifying the reservoir simulation time required for different number of parallel processes in a predictable manner through empirical equations; 4. quantifying the simulation time vs. number of grid cells (or size) through empirical equations; 5. describing the effect of upscaling in different directions on reservoir simulation accuracy; 6. describing the effect of natural fracture size and aperture on production profile.

Geology of the Hunton Reservoir
The areas of study are in Pottawatomie, Pontotoc, and Murray Counties, Oklahoma ( Figure 2). We conducted field measurements of the Hunton Anticline exposure as well as the Ada and Fittstown carbonate outcrops of the Chimneyhill subgroup to build a realistic fracture model for our upscaling and simulation studies. The subsurface data (logs, core, and seismic) comes from a well located in Pottawatomie County.
In the studied area (Figure 2), the Hunton Group Limestone is highly heterogeneous as inferred from the Dykstra Parsons coefficient [15]. This quantity is used to quantify the level of heterogeneity in a reservoir using permeability values from core data. Twenty-five Klinkenberg corrected core permeability values ( Figure 3) available in our study reveal a Dykstra Parson's coefficient of 0.902 (calculated using Equation (1)) which represents a very high reservoir heterogeneity. Lake and Jensen (1991) describe various other methods for the estimation of reservoir heterogeneity. surface with less than one degree dip [18,19]. In the area of study (Figure 2), the stratigraphic thickness ranges from 50 to 120 ft. Three primary lithologies of the Chimneyhill in the study area include wackestone, mudstone, and mud-dominated wackestone [20]. The Hunton Group facies comprise a series of progradational and aggradational parasequences (shallowing-upward cycles) [21]. Fractures in the Hunton Group are considered as one of the key components to enhance the hydrocarbon production [22]. Ghosh et al. [23] and Ghosh et al. [24] provide further details on the outcrop and subsurface fractures, current tectonic stress regime, and paleo tectonic stress regimes in the study area and surrounding areas in Oklahoma.    The Hunton Group is a prolific oil and gas-producing reservoir in the Midcontinent [16]. The cumulative oil production is 290 Million Barrels of Oil (MMBO) and the cumulative gas production is 5 Trillion Cubic Feet of Gas (TCFG) [17], with undiscovered 9 MMBO and 38 Billion Cubic Feet (BCF) of oil and gas, respectively [16]. The Hunton Group carbonates were deposited during Late Ordovician-Silurian-Early Devonian time in Oklahoma and is bounded by the underlying conformable Sylvan Shale and the overlying unconformable Woodford Shale ( Figure 4). The Hunton Group, including the Chimneyhill subgroup, is associated with a shallow-warm water carbonate ramp of a gently sloping surface with less than one degree dip [18,19]. In the area of study (Figure 2), the stratigraphic thickness ranges from 50 to 120 ft. Three primary lithologies of the Chimneyhill in the study area include wackestone, mudstone, and mud-dominated wackestone [20]. The Hunton Group facies comprise a series of progradational and aggradational parasequences (shallowing-upward cycles) [21]. Fractures in the Hunton Group are considered as one of the key components to enhance the hydrocarbon production [22]. Ghosh et al. [23] and Ghosh et al. [24] provide further details on the outcrop and subsurface fractures, current tectonic stress regime, and paleo tectonic stress regimes in the study area and surrounding areas in Oklahoma. Dykstra

Heterogeneous Model
The general workflow includes the following steps: 1. creation of the geologic model including the discrete fracture network (DFN), 2. upscaling the geologic and DFN model for subsequent flow simulation, 3. production history matching of the fine grid baseline model (Upscaled 1 in Figure 1), and repeated forward modeling of the production rate and volume at increasingly coarser grid sizes and fracture dimensions. Figures 5 and 6 describe the workflow used to populate the 3D volume, where we: (1) construct the relationship between the acoustic impedance logs, lithology (using thin sections), and porosity logs at well locations, (2) derive core porosity and vertical/horizontal permeability relationships (derived using laboratory experiments, (3) calculate fluid saturation from well logs, (4) invert the seismic data to estimate a compressional wave (P-wave) impedance volume, (5) correlate the P-wave impedance to the well log electrofacies and porosity measurements to construct a horizontal variogram to geostatistically construct 3D estimates of lithology and porosity, (6) build DFN models based on the outcrop field measurements, and finally, (7) upscale the fracture and matrix properties together for production simulations. The original geologic model using a 20 × 30 × 1 ft geologic grid resulted in 13,292,228 cells which exceeded the capabilities of our computer resources. We therefore upscaled the initial grid to 100 × 150 × 1 ft to form our baseline model, which was modified to match the production history. This baseline model was then further upscaled to coarser 200 × 300 × 1 ft, 400 × 600 × 1 ft, 800 × 1200 × 1 ft, and 1600 × 2400 × 1 ft grid cell sizes with a constant bed thickness of 1 ft for the horizontal upscaling to compare difference in production with the 100 × 150 × 1 model. Doubling the grid cell size in the X and Y directions (constant 1 ft bed thickness in Z direction) reduces the lateral heterogeneity and modifies the production rate and cumulative production. To understand the effect of vertical upscaling on production we coarsened the cell size to 100 × 150 × 2 ft, 100 × 150 × 4 ft, 100 × 150 × 8 ft, and 100 × 150 × 16 ft. Finally, to understand the effect of natural fractures on production, we increased the fracture lengths and widths by factors of two and four for the baseline 100 × 150 × 1 ft grid used in history matching.   and 1600 x 2400 x 1 ft grid cell sizes with a constant bed thickness of 1 ft for the horizontal upscaling to compare difference in production with the 100 x 150 x 1 model. Doubling the grid cell size in the X and Y directions (constant 1 ft bed thickness in Z direction) reduces the lateral heterogeneity and modifies the production rate and cumulative production. To understand the effect of vertical upscaling on production we coarsened the cell size to 100 x 150 x 2 ft, 100 x 150 x 4 ft, 100 x 150 x 8 ft, and 100 x 150 x 16 ft. Finally, to understand the effect of natural fractures on production, we increased the fracture lengths and widths by factors of two and four for the baseline 100 x 150 x 1 ft grid used in history matching.

Figure 6.
Workflow depicting history match. The reference fine grid has cells that measure Δx = 20 ft, Δy = 30 ft, and Δz = 1 ft. Flow was simulated for the coarser horizontally and vertically upscaled models, as well as for the models with upscaled fracture lengths and apertures.

Homogenous Model
We constructed homogenous models in order to quantify the production errors, resulting from the upscaling processes, which are associated with reservoir rock heterogeneity. The homogenous models were derived from the already-existing heterogeneous models using average values for the matrix and fracture parameters that were obtained from the fine grid 20 x 30 x 1 model ( Table 1).
The partial differential equations used for the reservoir simulators include the material balance equation, flow equations for three-phase flow of oil, water, and gas, and formation and fluid compressibility equations. The finite difference is a numerical method that is used to solve partial differential equations. However, the numerical method only gives an approximate solution to the discretized grid system of partial differential equations resulting in a truncation error [3,25,26]. Therefore, the production errors can be associated with both reservoir heterogeneity and numerical errors.
We quantified the effect of truncation error (numerical error) in the different simulations, first by making a series of homogenous vertical and homogenous horizontal upscaled models with different grid cell sizes (horizontal upscaling: 100 x 150 x 1 ft, 200 x 300 x 1 ft, 400 x 600 x 1 ft, 800 x 1200 x 1 ft, and 1600 x 2400 x 1 ft; vertical upscaling: 100 x 150 x 2 ft, 100 x 150 x 4 ft, 100 x 150 x 8 ft, and 100 x 150 x 16 ft grid cell sizes). We then performed the calculations as shown below. We used the 150 x 100 x 4 grid cell size as an example: Figure 6. Workflow depicting history match. The reference fine grid has cells that measure ∆x = 20 ft, ∆y = 30 ft, and ∆z = 1 ft. Flow was simulated for the coarser horizontally and vertically upscaled models, as well as for the models with upscaled fracture lengths and apertures.

Homogenous Model
We constructed homogenous models in order to quantify the production errors, resulting from the upscaling processes, which are associated with reservoir rock heterogeneity. The homogenous models were derived from the already-existing heterogeneous models using average values for the matrix and fracture parameters that were obtained from the fine grid 20 × 30 × 1 model (Table 1). Table 1. Reservoir simulation parameters for the homogenous models used in all vertical and horizontal upscaling and production simulation cases. Where: Φ = porosity; Kv = matrix permeability in vertical direction; Sw: water saturation; K I , K J , K K : fracture permeability in x, y, and z directions. Sigma factor: link between the matrix and fracture properties that describes fluid flow between the matrix and the dual-porosity/permeability model in a porous medium. If Sigma factor = 0, no communication between the matrix and natural fractures occurs.

Parameters Average Values Units Parameters Average Values Units
The partial differential equations used for the reservoir simulators include the material balance equation, flow equations for three-phase flow of oil, water, and gas, and formation and fluid compressibility equations. The finite difference is a numerical method that is used to solve partial differential equations. However, the numerical method only gives an approximate solution to the discretized grid system of partial differential equations resulting in a truncation error [3,25,26]. Therefore, the production errors can be associated with both reservoir heterogeneity and numerical errors.
We quantified the effect of truncation error (numerical error) in the different simulations, first by making a series of homogenous vertical and homogenous horizontal upscaled models with different grid cell sizes (horizontal upscaling: 100 × 150 × 1 ft, 200 × 300 × 1 ft, 400 × 600 × 1 ft, 800 × 1200 × 1 ft, and 1600 × 2400 × 1 ft; vertical upscaling: 100 × 150 × 2 ft, 100 × 150 × 4 ft, 100 × 150 × 8 ft, and Energies 2020, 13, 1604 9 of 27 100 × 150 × 16 ft grid cell sizes). We then performed the calculations as shown below. We used the 150 × 100 × 4 grid cell size as an example: Horizontal upscaling production error due to geologic heterogeneity = Horizontal upscaling production error due to numerical uncertainty = Horiz_upsc_prod_ err_ homog 150×100×4 , i.e., same as (3) Percentage errors for all grid cell sizes are calculated with respect to the history match grid cell size (i.e., 150 × 100 × 1 ft). Similar calculations may be performed for the error in vertical upscaling. However, this limitation is overcome using percentage errors instead of actual errors as shown in the equations above.

Estimation of Reservoir Parameters
Electrofacies form one of the key components of the models described in Figure 5. We begin with Principal Component Analysis (PCA) of the electric logs. These values serve as input to subsequent Self-Organizing Mapping (SOM) to construct electrofacies logs. We corroborate the geologic interpretation of the electrofacies predictions using thin sections and borehole images. We constructed the relationship between the porosity and vertical and horizontal permeabilities from cores. Porosity logs were correlated with core porosity measurements which in turn are linked to the facies distribution ( Figure 7e). We estimated fluid saturation by applying Archie's law (Equation (6)) on the well logs. We measured porosities as well as vertical and horizontal permeabilities from core plugs facilitating the construction of porosity-permeability relationships. Following previous work by Milad et al. [11], Milad and Slatt [20], Milad et al. [27], we constructed vertical variograms from the well logs to populate the petrophysical data spatially.

Seismic Data
We used seismic data to delineate formation boundaries, to define major faults, and to provide a volumetric estimation of compressional (P-wave)-and shear (S-wave) impedances. These impedances were crucial in differentiating fossiliferous wackestone from dolomitized mudstone rocks. The low frequency background model used in impedance inversion used five wells (Figures 2c and 7a) with P-wave velocity, S-wave velocity, and density logs and ranged from the top of the Woodford surface to the top of the Viola surface ( Figure 4). Six common angle gathers ranging from 0 o to 30 o at 5 o increments were used for prestack inversion using © Hampson-Russell Software. Detailed descriptions about our inversion methodology and its results are available in Milad et al. [28] and Milad [29]. The main purpose of the inversion is to provide horizontal variograms for porosity distribution. Another purpose is to find the major and minor directions of lithologic variation (Figure 7b), impedance (Figure 7c), and the facies distribution (Figure 7d). Figure 7f shows a DFN model for one of the four fracture sets. We considered four main parameters in the creation of the fracture network: intensity, geometry, orientation, and aperture. Parameters defining the fracture intensity and dimension can be found in Milad et al. [27]. The mean dip and mean dip azimuths were assigned for each fracture set based on the borehole measurements. We assigned the statistical distributions of the borehole-measured fracture apertures (calibrated to the maximum aperture value observed in the core) to the model. We used the stochastic technique described by Schlumberger [30] in © Petrel Software manual to create a DFN volume for the entire seismic cube with a distribution guided by the fracture intensity map (grid cell intensity values) for each set, where the Fracture porosity = fracture aperture/fracture spacing (7) Natural fracture permeabilities were calculated using a cubic aperture (y) law

Upscaling
We follow Oda [31] to upscale the fracture porosity, the fracture sigma length in the x, y, and z directions (spanned by J, K, and I grid indices), and the fracture permeability. Assuming Darcy's Law for laminar flow through an isotropic porous medium, Oda [31] shows how to upscale the geometry and properties of the fractures within a given cell to construct the permeabilities of an equivalent grid cell. Projecting the fracture isotropic permeability onto the fracture plane and then scaling it to the ratio between the fracture volume and grid cell volume provides a net permeability tensor for each grid cell. Eigenvectors of the net permeability tensor provide the principal permeabilities.
The simulation software requires two simulation cells to model a dual porosity/dual permeability system where adjacent grid cells communicate through both the fractures and the matrix. This construct allows one to define porosity and permeability where the transmissibilities are connected or coupled to each other using a quantity called the sigma factor [30].

Geocellular Gridding and History Matching
The grid cells are aligned with an orthogonal X, Y, Z Cartesian coordinate system. We used small vertical cell sizes to accommodate the heterogeneities observed in the well logs. The relative grid size is a function of the heterogeneity, where coarser spacing can be used if the reservoir properties are more homogeneous in one direction (Figure 7b). We used proportional slicing between the top and base correlation surfaces to divide the reservoir into appropriate layers. Based on P-impedance from seismic inversion, we rotated the computational grid to align with the direction of the depositional trend (30 ft in the major direction of North 15 East (N15E) [x] and 20 ft in the minor [y] direction) (Figure 7b).
The history matching process was undertaken to adjust the uncertain model parameters of fluid reservoir properties, including well production rates, reservoir pressure and reservoir temperatures ( Table 2) without changing the properties of static geological models. Once the fluid model parameters and well data constraints were refined and individual well production histories were matched in the 100 × 150 × 1 ft model, we used the same reservoir simulation parameters for the subsequent upscaled  Figure 7a shows the trajectories of five wells used in this study. Geologic models in Figure 7b-f were obtained using the method described in Figure 5, i.e., by utilization of core, well log, outcrop, and seismic data. Figure 8a shows a 3-dimensional porosity model that has been upscaled to 150 × 100 × 1 ft grid from a 20 × 30 × 1 ft (Figure 7e) grid. Figure 8b shows the water saturation calculated from the well logs, upscaled onto the 20 × 30 × 1 ft grid, and then subsequently upscaled onto the 150 × 100 × 1 ft grid. Figure 8c,d show the upscaled combined fracture and matrix permeabilities.  Table 2 provides detailed information and properties of the reservoir.  Table 3 shows total errors resulting from the truncation computation and geologic heterogeneity. Table 4 shows numerical errors due to computation truncation from homogenous models. Table 5 shows the effect of geologic heterogeneity on the effect of production using Equations (1-3) as mentioned earlier. The effect of horizontal upscaling on production is highly variable. Table 3. Cumulative production comparisons in horizontal upscaling cases for the heterogeneous models. Well by well actual cumulative production, history-matched (100 × 150 × 1 grid) cumulative production, and comparisons of larger grid cell sizes with history-matched cumulative production are shown. Errors in the history match column are with respect to the actual cumulative production and those in remaining columns are with respect to the history-matched production.   Table 5. Cumulative production comparisons in horizontal upscaling cases after subtracting the numerical errors in Table 4 from total errors in Table 3  For the total (numerical plus heterogeneity) error cases, Well 3 ( Figure 9c) shows little variation (-4%) in cumulative production (Table 3) between the 100 × 150 × 1 (history match) green line and the next larger grid (200 × 300 × 1 ft). In other wells (Figure 9a,b,d,e), the next larger grid (200 × 300 × 1 ft) shows a significant change (total error) in production rate (+22 to + 131%) ( Table 3). The remaining larger grid cell sizes (400 × 800 × 1 ft, 800 × 1600 × 1 ft, and 1600 × 2400 × 1 ft) induced significant differences in cumulative production and production rates when compared to the history-matched grid cell size (i.e., +22% to +537%) [ Figure 9a-f, Table 3]. Note that progressively larger grid cell sizes do not necessarily result in a progressively larger departure from the history match (green curve) for individual wells. However, such a correlation can be observed for the cumulative production in all wells combined (Table 3, Table 4, and Table 5). For example, the orange (800 × 1200 × 1 ft) and black (1600 × 2400 × 1 ft) curves swap position from one well to the other (Figure 9a-e).

Actual
Energies 2020, 13, x FOR PEER REVIEW 7 of 29 grid cell sizes do not necessarily result in a progressively larger departure from the history match (green curve) for individual wells. However, such a correlation can be observed for the cumulative production in all wells combined (Tables 3, 4, and 5). For example, the orange (800 x 1200 x 1 ft) and black (1600 x 2400 x 1 ft) curves swap position from one well to the other (Figures 9a-9e).   Table 3. Cumulative production comparisons in horizontal upscaling cases for the heterogeneous models. Well by well actual cumulative production, history-matched (100 x 150 x 1 grid) cumulative production, and comparisons of larger grid cell sizes with history-matched cumulative production are shown. Errors in the history match column are with respect to the actual cumulative production and those in remaining columns are with respect to the history-matched production.   Table 6 shows total errors, Table 7 shows numerical errors, and Table 8 shows errors due to geologic heterogeneity for the vertical upscaling cases. For the total (numerical plus heterogeneity) error cases, Figure 10a-c,e-f show little difference between the history match (green curves) and larger grid cell curves. In addition, essentially all vertically upscaled matches may be considered as reasonable history matches to the observed data (Figure 10a-c,e-f). Larger grid cell size always leads to a jump in the cumulative production (Table 6) and production rate compared to the base history-match grid cell size of 100 × 150 × 1 ft. Vertical upscaling makes a noticeable difference in only well 4 (Figure 10d) which exhibits a noticeable 17% deviation in production rate and cumulative production ( Table 6)  (f) Sum of production rates for all wells (field production). The difference in the field production dates is due to different start time of production for these five wells. Dots: actual production rates; production rates for different grid cell sizes-green curves: 100 x 150 x 1 (approximate history match); blue curves: 100 x 150 x 2; red curve: 100 x 150 x 4; orange curves: 100 x 150 x 8; black curves: 100 x 150 x 16 of cumulative production for individual wells. Table 6. Cumulative production comparisons in vertical upscaling cases for the heterogeneous models showing total (numerical uncertainty plus geologic heterogeneity). Well by well actual cumulative production, history match (100 x 150 x 1) cumulative production, and comparisons of larger grid cell size with history match cumulative production are shown. Errors in the history match column are with respect to the actual cumulative production and those in remaining columns are with respect to the history match (i.e., not actual production) case (100 x 150 x 1).  (f) Sum of production rates for all wells (field production). The difference in the field production dates is due to different start time of production for these five wells. Dots: actual production rates; production rates for different grid cell sizes-green curves: 100 × 150 × 1 (approximate history match); blue curves: 100 × 150 × 2; red curve: 100 × 150 × 4; orange curves: 100 × 150 × 8; black curves: 100 × 150 × 16 of cumulative production for individual wells. Table 6. Cumulative production comparisons in vertical upscaling cases for the heterogeneous models showing total (numerical uncertainty plus geologic heterogeneity). Well by well actual cumulative production, history match (100 × 150 × 1) cumulative production, and comparisons of larger grid cell size with history match cumulative production are shown. Errors in the history match column are with respect to the actual cumulative production and those in remaining columns are with respect to the history match (i.e., not actual production) case (100 × 150 × 1).  Table 7. Cumulative production comparisons in vertical upscaling cases for the homogenous models showing errors due to numerical uncertainty. Well by well actual cumulative production and comparisons of larger grid cell sizes with (100 × 150 × 1 grid) cumulative production are shown.  Table 8. Cumulative production comparisons in vertical upscaling cases after subtracting the numerical errors in Table 7 from total errors in Table 6

Comparison
The difference in production rates from vertical upscaling is not nearly as dramatic as that observed from horizontal upscaling (Figure 11a,b), though both follow a logarithmic relationship.  Figure 11. Comparison between vertical and horizontal upscaling production error. Total errors that are associated with numerical errors and rock heterogeneity due to upscaling (orange), numerical (uncertainty) error resulting from the truncation computational errors (black), and upscaling errors (rock heterogeneity) resulting from changes in grid cell sizes (green). Some orange points are overlapping black points in the upper plot (a) and black points are overlapping each other on the lower plot (b). a) Semilog plot depicting logarithmic best-fit equations (inversely proportional) for the three horizontal upscaling curves (orange, black, and green). Note numerical production errors due to the computational truncation errors (black) are higher than the production errors due to the upscaling of geologic heterogeneity (green) for grid cell sizes. The exact error values for the horizontal upscaling are found in the last columns ("all wells combined") in Table 3, Table 4, and Table 5. b) Semilog plot depicting logarithmic best-fit equations (inversely proportional) for the total errors and numerical errors vertical upscaling. Note numerical errors with the vertical upscaling is very low (black). Note much lower numerical errors with vertical upscaling (black) and higher errors with horizontal upscaling (black). The exact error values for the vertical upscaling are found in the last columns ("all wells combined") in Table 6, Table 7, and Table 8. Table 9 and Figure 12a-f show that increase in aperture has a lesser effect on production than the length. For example, Well 1 shows a difference of 11% in cumulative production compared to the history match using 2x aperture, and 301% using 2x length (Table 9, also see Figure 12a solid green, solid blue, and solid red curves). Well 1 also shows 190% increase compared to the history match using 4x aperture and 673% increase using 4x length (Table 9, also see Figure 12a solid green, broken blue, broken red curves). For Well 2, 2x fracture aperture results in a 63% increase, while 2x fracture length results in 861% increase in production. Similarly, for 4x increase in aperture and length, the increase is 169% and 1595% in production respectively. With 4x aperture and 4x length, the jump in production rate with respect to the respective 2x cases is much higher than the difference between the respective 2x and history match cases. The remainder of the wells show a similar increase for 2x cases. The 4x values indicated by the arrows in Figure 12c-e (constant values) should be ignored as the model boundaries were most likely reached. In addition, the cumulative production of the whole field for 4x aperture and length (arrows in Figure 12f) should be ignored as that includes all five wells. Table 9. Cumulative production comparisons for history match, 2x, and 4x fracture aperture and fracture lengths. Errors in the history match column are with respect to the actual cumulative production and those in the remaining columns are with respect to the history match (i.e., not actual) production case (100 × 150 × 1). upscaling of geologic heterogeneity (green) for grid cell sizes. The exact error values for the horizontal upscaling are found in the last columns ("all wells combined") in Tables 3, 4, and 5. b) Semilog plot depicting logarithmic best-fit equations (inversely proportional) for the total errors and numerical errors vertical upscaling. Note numerical errors with the vertical upscaling is very low (black). Note much lower numerical errors with vertical upscaling (black) and higher errors with horizontal upscaling (black). The exact error values for the vertical upscaling are found in the last columns ("all wells combined") in Tables 6, 7, and 8. Table 9 and Figure 12a-f show that increase in aperture has a lesser effect on production than the length. For example, Well 1 shows a difference of 11% in cumulative production compared to the history match using 2x aperture, and 301% using 2x length (Table 9, also see Figure 12a solid green, solid blue, and solid red curves). Well 1 also shows 190% increase compared to the history match using 4x aperture and 673% increase using 4x length (Table 9, also see Figure 12a solid green, broken blue, broken red curves). For Well 2, 2x fracture aperture results in a 63% increase, while 2x fracture length results in 861% increase in production. Similarly, for 4x increase in aperture and length, the increase is 169% and 1595% in production respectively. With 4x aperture and 4x length, the jump in production rate with respect to the respective 2x cases is much higher than the difference between the respective 2x and history match cases. The remainder of the wells show a similar increase for 2x cases. The 4x values indicated by the arrows in Figure 12c-12e (constant values) should be ignored as the model boundaries were most likely reached. In addition, the cumulative production of the whole field for 4x aperture and length (arrows in Figure 12f) should be ignored as that includes all five wells.

Run Time and Accuracy
There is a fourth order negative relationship between solution time vs. grid cell size (Figure 13a), and solution time vs. number of processes (Figure 13b). In other words, using five parallel processes does not reduce runtime to one-fifth of the time taken by one processor. For example, for the history match case, five processors consume nearly 13 hours which is one-third, rather than one-fifth, the time taken by a single processor. For the 200 × 300 × 1 grid cell size, this ratio is different as for the remaining grid cell sizes. The time-gap closes with an increase in the number of processes (Figure 13b).  x 300 x 1 (blue), 400 x 600 x 1 (red), and 800 x 1200 x 1 (orange). All curves were best fit by negative power law equations.  Figure 13. Time required for flow-simulation completion. (a) Log-log plot depicting simulation time required using a different number of grid cells (or grid cell size) while applying different numbers of processes: a single (blue), five (orange), ten (red), twenty (black) processes. All curves were best fit by a power law. (b) Plot (both x and y-axes are linear) depicting simulation time using different number of processes (curves depicting single grid cell sizes) while using various grid cell sizes: 100 × 150 × 1 (green), 200 × 300 × 1 (blue), 400 × 600 × 1 (red), and 800 × 1200 × 1 (orange). All curves were best fit by negative power law equations.

Discussion
The aforementioned data integration at different scales from core to seismic volume were used to develop a high-resolution stratigraphic framework and complex structures to build detailed fine grid cell size static geological models consisting of lithofacies, natural fractures, porosities, permeabilities, and saturation. The resulting detailed (20 × 30 × 1 ft grid cell size) reservoir model was upscaled to 100 × 150 × 1 ft grid cell size for flow (reservoir) simulation. We used a dual porosity/dual permeability model for higher accuracy even though the calculations require higher time for upscaling and flow simulation.
A non-linear relationship between the grid cell size (and consequently the number of grids) and run duration was observed. The runtime gap closes with larger grid cell size as the number of grid cells decreases. Therefore, if the grid cell size is large enough, parallel processes may not be required. The increase in grid cell size, however, results in unacceptably large errors in the simulation results. Therefore, the optimal grid cell size and the run time taken based on the number of available processes requires optimization. Comparison between vertical and horizontal upscaling shows greater error with horizontal upscaling. In this study, the vertical grid cell can be reasonably upscaled to 100 × 150 × 2 ft or 100 × 150 × 4 ft without incurring more than 6% error compared to the history match (100 × 150 × 1 ft). Larger vertical grid cell sizes, i.e., up to 10 ft may be used without incurring more than 11% error. In other words, a reasonable history match and forward modeling can be obtained with 1-10 ft vertical grid cell size showing that thin beds do not affect production in the current study. Increase in lateral grid cell size, however, results in a substantial total error (0.5% to 166%) that includes upscaling processes and numerical errors. The magnitude of numerical truncation (i.e., numerical uncertainty or instability) errors ranges from 4.77% to 91.78% and is inversely proportional to the number of grid cells, or directly proportional to the grid cell size, in the model. Thus, the effect of numerical errors can be eliminated by decreasing the grid cell sizes in the model. These results agree with Ezekwe [3] and Lantz [25] where the numerical errors increase with increases in the areal grid cell sizes. Increase in lateral grid cell size due to upscaling results in errors in the range of 0.5%-74.2% compared to the vertical grid cell sizes errors which are in the range of 0.5%-10.8%. This observation is counterintuitive, given more frequent changes in rock properties in the vertical direction than in the lateral directions. However, it is also noteworthy that the vertical cell size is much smaller compared to horizontal cell size. Moreover, each step in vertical upscaling increases the grid cell volume by a factor of two, while each step in horizontal upscaling increases grid cell volume by a factor of four (i.e., twice in the x-direction and twice in the y-direction). Based on the overall production, a lateral grid cell size of more than the minimum required for successful simulation is not recommended.
The equation in Figure 11, i.e., the one defining a monotonous (logarithmic) increase in error with grid cell size (or decrease in the number of grid cells), is presented only for reference and likely to vary in different areas based on geology. The error can also decrease with an increase in grid cell size as observed with individual wells. However, the equation in Figure 13 (grid cell size vs. time), i.e., power law relationship between grid cell size and solution time, is likely to hold true at other places in other geological settings.
Regarding the effect of natural fracture geometry on production, natural fracture length has a higher impact on the production compared to the aperture. Therefore, assuming a constant natural fracture intensity or surface area, greater production may be expected with longer natural fractures. This observation likely results from a higher reservoir connectivity through longer natural fractures. Since fracture length can affect reservoir production significantly [32], its accurate estimation is paramount, even though it is very difficult to achieve [10]. An increase in natural fracture aperture, on the other hand, does not result in higher natural fracture connectivity but still increases production. This observation indicates that natural fracture conductivity is not infinite. Moreover, constant high production rates with 4x apertures imply a shift from finite conductivity (1x and 2x aperture cases) to infinite conductivity leading to model boundary dominated flow.

Conclusions
This study combines several geologic, geophysical, and engineering data and parameters methodically to establish the feasibility of such effort in understanding simulation accuracy and practicability. We demonstrated that it is essential for modelers to understand the sources of error (i.e., numerical uncertainty vs. reservoir heterogeneity) and acquire ability to estimate the relative error magnitudes arising from these sources. Additionally, we found that vertical upscaling is a better option than horizontal upscaling for preservation of computational accuracy. For horizontal upscaling, we showed that the lateral grid cell size should be kept at the original bin size (~110 × 110 ft) of the seismic data used in its construction in order to balance runtime and accuracy. Moreover, reservoir connectivity attributable to natural fracture length was found to have a significantly greater impact on the production response compared to natural fracture aperture. A less intuitive phenomenon (i.e., overhead of inter process communication) causing a non-linear relationship between the number of parallel processes and the total simulation time was demonstrated. Our workflow of using a fine grid and comparing the history match to those of a coarser grid exhibits a means of estimating the accuracy of a given simulation. Funding: This research was funded by STACK-MERGE-SCOOP-consortium members at the Institute of Reservoir Characterization at the University of Oklahoma. Check carefully that the details given are accurate and use the standard spelling of funding agency names at https://search.crossref.org/funding, any errors may affect your future funding.