Dynamic effects of fault modeling on stair-step and corner-point grids

Fault modeling has become an integral element of reservoir simulation for structurally complex reservoirs. Modeling of faults in general has major implications for the simulation grid. In turn, the grid quality control is very important in order to attain accurate simulation results. We investigate the dynamic effects of using stair-step grid (SSG) and corner-point grid (CPG) approaches for fault modeling from the perspective of dynamic reservoir performance forecasting. We have performed a number of grid convergence and grid-type sensitivity studies for a variety of simple, yet intuitive faulted flow simulation problems with gradually increasing complexity. We have also explored the added value of the multipoint flux approximation (MPFA) method over the conventional two-point flux approximation (TPFA) to increase the accuracy of reservoir simulation results obtained on CPGs. Effects of fault seal modeling on grid-resolution convergence and grid-type sensitivity have also been briefly examined. For simple geometries, both SSG and CPG can be used for fault modeling with similar accuracy in conjunction with the pillar-grid approach. This is evidenced by the fact that simulation results from SSG and CPG converge to identical solutions. SSG and CPG yield different results for more complex geometries. Simulation results approach to a converged solution for relatively fine SSGs. However, a SSG only provides an approximation to the fault geometry and reservoir volumes when the grid is coarse. On the other hand, non-orthogonality errors are increasingly evident in relatively more complex faulted models on CPGs and such errors cannot be addressed by grid refinement. It has been observed that MPFA partially addresses the discretization errors on non-orthogonal grids but only from the total flux accuracy perspective. However, transport related errors are still evident. Grid convergence behaviors and grid effects are quite similar with or without fault seal modeling (i.e., dedicated fault-zone modeling by use of scaled-up seal factors) for simple geometries. However, in more complex test cases, we have observed that it is more difficult to achieve converged results in conjunction with fault seal modeling due to increased heterogeneity of the underlying problem.


Introduction
Faults are associated with reservoirs of oil and gas in many basins of the world. A fault is a break or planar surface in brittle rock across which there is observable displacement. Faults can strongly alter fluid flow and can act both barriers and conduits. According to the stress state that produces faults, there are different types of faults: normal faults, reverse faults, strike-slip faults, dip-slip faults, and obliqueslip faults. Structural settings differ globally, and thus different types of faults are incorporated into geological models to realistically define the structural framework. A normal fault is shown in Fig. 1a and a reservoir model with fault groups is illustrated in Fig. 1b, respectively. The fault density varies greatly in reservoirs from none to a very high degree. When flow simulation models are constructed, it is relatively simpler to generate a grid when there are only a few faults in the model. It becomes increasingly more challenging to model and grid structurally complex reservoirs with intersecting and bifurcating faults such as the one in Fig. 1b.
Representation of faults in static modeling packages can have a significant impact on the structure of the resulting simulation grid. In common static modeling packages, there are different gridding options available to customize the 3D grid for geomodeling and flow simulation purposes, e.g., the corner-point and the stair-step gridding options (Fig. 2). Construction of the structural model involves a process called the pillar gridding (e.g., Hoffman et al. 2003Hoffman et al. , 2006. This is a concept in static modeling packages where the faults in the fault model are used as a basis for generating the 3D grid used for static reservoir modeling. Accurate fault modeling either using fault sticks or polygons provide the first step to pillar gridding. The process of representing the fault planes as pillars of the grid could result in a 3D grid containing a large number of (severely) non-orthogonal gridcells. A gridcell is non-orthogonal when it is skewed. In other words, grid non-orthogonality arises in hexahedral grids (commonly used by a large number of reservoir simulators) when the gridcells depart from being a rectangular prism.
Severely non-orthogonal grids lead to large discretization errors in the solution of flow equations by reservoir simulators that utilize the TPFA (e.g., Wu and Parashkevov 2009). TPFA is commonly used as the default discretization method in commercial and proprietary reservoir simulators owing to its computational efficiency. Here, it is important to highlight that there are discretization techniques that accurately solve flow equations on non-orthogonal grids by taking into account the diagonal flow, e.g., MPFA (Aavatsmark et al. 1996), mimetic method (Alpak 2010), and others. However, these discretization schemes are usually computationally more intensive than TPFA. On the other hand, they effectively eliminate the errors from non-orthogonality. A variant  of the MPFA scheme, namely the O-method (Aavatsmark et al. 1996(Aavatsmark et al. , 1998, is available in addition to TPFA in the proprietary reservoir simulator used in the studies described in this paper. Orthogonal grids can be generated by using a method where a fault is represented by a number of stair-steps in the vertical direction. This SSG approach automatically fulfills the orthogonality assumptions of TPFA. However, this approach leads to faults in the model that do not accurately follow the surfaces as defined by seismic interpretation (Fig. 2a).
On the other hand, the CPG approach leads to modeled faults that more accurately follow the interpreted surfaces but results in a non-orthogonal grid. The severity of the non-orthogonality depends on the number of faults, their orientations, and the skill of the geomodeler. Another disadvantage of CPG is that it is difficult to model geometries encompassing fault groups because it has a fixed number of six-sided gridcells. In this circumstance, a combination of SSG and CPG is generally used in gridding.
A summary of comparison between SSG and CPG is documented in Table 1. For SSG, the fault is modeled as a zigzag surface. As a result, the ensuing grids are uniform Cartesian grids, which are orthogonal. The flux calculation through the fault surface is accurate (for TPFA), however, the geometry and the volumes are only approximate. As a remedy, grid refinement can be used to render SSGs follow fault surfaces more closely. Grid refinement in the static modeling package in general provides a more effective solution than performing a global or local grid refinement in the simulator because the static modeling package also allows the refinement of the flow properties (permeability, porosity, net-to-gross ratio, etc.) populating these gridcells for a more accurate representation. Results from SSG are also affected by grid orientation. The gridlines used for constructing a SSG are recommended to be parallel to the line connecting the injector and the producer, not diagonal, which means the resulting grids provide more accurate solutions when aligned with the main flow directions. Another disadvantage of SSG is its inflexibility in the representation of well locations.
In the literature, there are many comparisons between simulation results of SSG and CPG. For example, for cases with and without faults, He and Durlofsky (2006) demonstrated that the results of CPGs with grid quality control are more accurate than Cartesian grids when total grid numbers are almost the same. Wu and Parashkevov (2009) showed that SSGs result in a smaller error than the deviated or verticalized grids. The usage of SSG and CPG is case dependent. Therefore, it is beneficial to understand and predict the circumstances that one gridding approach is better than the other. The contents of our paper are as follows: The paper proceeds with the description of the fault modeling workflow using the pillar gridding approach that give rise to SSGs and CPGs. The rest of the paper is dedicated to simple and intuitive numerical studies that involve SSGs and CPGs. Our investigations feature grid converge studies, comparisons of simulation predictions using SSGs and CPGs, investigations of the effect of MPFA on CPGs, and fault seal modeling. The paper is closed with the main conclusions reached as a result of our study.

Fault modeling protocol
The pillar-grid fault models and static reservoir models were constructed using a commercial static modeling package (Schlumberger Information Solutions 2004). The resulting static model was exported into a proprietary static-todynamic model processing package, which can be used for both upscaling and downscaling operations (e.g., Alpak and Vink 2018). Using the geological model as input, reservoir simulation was conducted in a proprietary simulator (Por et al. 1989;Regtien et al. 1995). Both TPFA (e.g., Aarnes et al. 2007) and MPFA O-method (Aavatsmark et al. 1996(Aavatsmark et al. , 1998) discretization techniques were investigated. Fault modeling and gridding were performed with the following steps: 1. Faults traces and horizons (reservoir top and bottom) are imported into the static modeling package. 2. A structural model is defined using the geometric modeling capabilities. 3. Fault modeling: Fault traces are converted to key fault pillars in the fault modeling stage. 4. Pillar gridding: Pillar gridding is the process of defining the 3D grid, which is a 2D mesh extended into the third dimension. A 2D grid is defined by rows and columns oriented in the horizontal ( x andy ) directions. On the other hand, a 3D grid is defined by rows, columns, and pillars oriented in thex , y , and z directions. In other words, a 3D grid is a series of 2D grids stacked on top of one another. Pillars are the lines connecting corresponding grid nodes from each of the 2D grids. The pillar gridding process starts with a set of rows and columns evenly spaced at a specified grid increment. The pillars are vertical lines passing through each row-column intersection at this stage. The previously defined key pillars delineating the fault geometries guide the reorientation of these initial pillars to create the final 3D grid. Through a series of iterations guided by an automated optimization process, the pillars, rows, and columns are made sub-parallel to the key pillars in the final 3D grid. SSG or CPG grid types can be specified in the pillar gridding process. Both grid types are generated using the fault trend as guideline. 5. Horizons generation: Depth information is added to the 3D grid by associating it with inputs such as depth maps and/or well tops. 6. Layering: The geological model is generated with one layer (top-to-bottom). It is later divided into 10 layers through downscaling in the static-to-dynamic model processing package. 7. Property modeling is performed in the simulator because we use simple properties in the studies reported in this paper. 8. The model is quality checked by carefully inspecting the cross sections and editing problem areas in the static model and repeating steps #3 to #7. 9. Optionally, seal factors (transmissibility multipliers) are generated for every gridcell on fault surfaces by use of a fault seal modeling algorithm, which is an effective property-based modeling approach for incor-

Numerical experiments
We describe the results of a simulation-based study encompassing three cases: (1) a single fault with flat horizons, (2) a single fault with non-flat horizons, and (3) two interior faults with flat horizons. Complexity was built into geological models incrementally. For each case, geomodels were set up in the static modeling package, and both SSGs and CPGs were generated. Subsequently, grid convergence tests were conducted through flow simulations on all cases because there was a need to identify a converged solution and determine an appropriate grid density that strikes the best balance in terms of delivering both accuracy and efficiency. The effect of gridding was investigated by comparing the results from simulations on SSGs and CPGs. TPFA was used in all the cases except for those in a dedicated section, in which MPFA was introduced to improve the results of CPG cases with relatively severe grid non-orthogonality. Note that fault seal modeling was not performed in the initial test cases. It was only used at a later stage, where the effect of gridding with fault seal function was investigated for different flow scenarios.

Case #1: model with a single fault and flat horizons
The geological model for this test case has a north-south normal fault and flat horizons (Fig. 3a). The fault has a small throw compared to the reservoir thickness as shown in the cross section in Fig. 3b. For this faulted reservoir, SSGs and CPGs were generated (Fig. 4). The 3D reservoir has homogeneous properties: isotropic permeability = 100 mD, porosity = 0.25, net-to-gross ratio (NTG) = 0.5. The size of the reservoir is about 7800 ft × 5250 ft × 100 ft. Our base case is CPG with 26 × 15 × 10 gridcells and gridcell sizes of Δx = Δy = 50 ft and Δz = 10 ft. The static model was converted to a simulation model through one-to-one upscaling. All simulation models have ten layers. Two wells, specifically one vertical producer and one vertical injector, were placed in the model as shown in Fig. 4. The wells penetrate through the entire model and were perforated in all of the layers. Absolute coordinates were used to specify the well locations so that they do not change for different gridding systems, such as during grid refinement. Both wells operate under bottom-hole pressure (BHP) constraints. The maximum BHP constraint for the injector is 6000 psi and the minimum BHP constraint for the producer is 4020 psi. The initial reservoir pressure is 5000 psi. Given that the bubble point is 4014.7 psi, the reservoir oil remains undersaturated during production. The maximum rate of oil production is 5000 bbl/day. The maximum water injection rate is 5000 bbl/day. Since fault seal modeling was not included, wherever there is sand-to-sand contact, the seal factor is 1.0 (open to flow); otherwise, it is 0.0 (closed to flow).
Oil production rate versus. time and average reservoir pressure versus. time results are plotted in Fig. 5 for the base case. The average reservoir pressure drops during the first year of production, because the water injection is not sufficient to replace production. As a result, the oil production decreases quickly in the first year after which it continues to drop slowly.
The saturation map after a 32-year production shows the water from the injector moves across the fault to the producer (Fig. 6a). Pressure declines around the producer and increases around the injector (Fig. 6b). As mentioned previously, the fault seal factor is 1.0 when there is sandto-sand contact. For this case, the fault seal factor is 1.0 everywhere and the fault throw is much less than the reservoir thickness. Therefore, the impact of the fault is negligible, which is reflected by smooth color maps of saturation and pressure fields.
An analysis was performed on grid convergence. The grid used in the base case was only areally refined, not vertically. The coarse base-case grid of size 26 × 15 × 10 was used as the starting point, and it was refined in each horizontal direction by a factor of approximately 2 each time. In Fig. 7, oil production rate versus. time profiles are documented for four levels of refinement: level 1 (2 0 times refinement, grid size: 26 × 15 × 10), level 2 (~2 1 times refinement, grid size: 51 × 31 × 10), level 3 (~2 2 times refinement, grid size: 103 × 62 × 10), and level 4 (~2 3 times refinement, grid size: 205 × 125 × 10). The oil production rate converges to an almost unified solution when the grid is refined. Average reservoir pressure versus. time profiles are plotted for four The results indicate that pressure is not sensitive to grid size, because both the injector and the producer are at BHP control. Using the results from the finest grid as reference, relative errors were calculated at the end of the year ( Table 2). The magnitude of the error in oil production rate decreases from 8% to 4%, and then to 1% with each refinement. Errors in average reservoir pressure are within 1% for all the cases. Therefore, the simulation results for the use of CPGs exhibit convergence for Case #1.  Fig. 7 a Oil production rate and b average reservoir pressure as a function of time at four levels of refinement for Case #1 using CPGs. The orange line is the base case; green line is ~2 times refinement; red line is ~4 times refinement; and blue line is ~8 times refinement. c Zoomed-in view of the average reservoir pressure profiles Similarly, a convergence test was performed on SSGs. Oil production rate converges very fast (Fig. 8). Average reservoir pressure is insensitive to the grid size. The errors decrease as the grids are refined (Table 3). Both the ~2 times refined-grid (in each horizontal direction) and the ~4 times refined-grid yield reasonably accurate results. The ~2 times refined-grid results were selected for further analysis because it sufficiently strikes the best balance in terms of simulation accuracy and turnaround time.
After the grid convergence test, solutions for different gridding approaches were compared. The solutions are based on the second-coarsest grids (i.e., 51 × 31 × 10 for CPG and 52 × 30 × 10 for SSG) in each case. In Fig. 9, oil production rate as a function of time profiles are plotted for both SSG and CPG. It was observed that the results are almost the same for CPG and SSG. It indicates that these cases converge to the "true" solution, which is an expected result since this test case is very simple with one north-south oriented fault and a very small throw. It can be concluded that the gridding-approach effect for Case #1 is negligible.  . 8 a Oil production rate and b average reservoir pressure as a function of time at four levels of refinement for Case #1 using SSGs

Case #2: model with a single fault and non-flat horizons
Case #2 is more complicated than Case #1. The geological model has a moderately dipping normal fault with a varying throw (Fig. 10). The horizons are shown in Fig. 10a, and they are non-flat. Two cross sections, shown in Fig. 10b and Fig. 10c, illustrate an example completely offset part of the fault zone and an example part of the fault zone with sand-to-sand contact in the model, respectively. It can be seen that the lower part of the fault zone is closed to cross-fault flow, and the upper part permits cross-fault flow. Two types of grids were generated in the static modeling package, namely CPGs and SSGs (Fig. 11). The base case uses the CPG with 22 × 36 × 10 gridcells and gridcell sizes of Δx = Δy = 100 ft and Δz = 20 ft.
The 3D reservoir has homogeneous properties: isotropic permeability = 100 mD, porosity = 0.25, net-togross ratio = 0.5. The size of the reservoir is about 5000 ft × 6700 ft × 200 ft. The well constraints are the same as in Case #1.
A single injector -producer pair was placed in the reservoir as shown in Fig. 11. Oil production rate versus. time and average reservoir pressure versus. time profiles are reported in Fig. 12. Similarly, as in the Case #1, during the first year of production, both oil production rate and average reservoir pressure decrease as a function of time.
The saturation map after 32 years of production shows the water from the injector finds a path where the fault is connected and then moves toward the producer (Fig. 13a). The reservoir pressure decreases around the producer and increases around the injector (Fig. 13b). Similarly, as in Case #1, a grid sensitivity analysis was conducted. In Fig. 14, oil production rate versus. time and average reservoir pressure versus. time profiles are shown for four levels of refinement. The error in oil production rate decreases from − 19% to − 4%, and then to − 1% with each refinement for the CPG approach (Table 4). Errors in average reservoir pressure are within 1% for all investigated cases. Therefore, the results converge after the grid sizes become increasingly smaller for Case #2 on CPGs. A very large error is observed for a grid size of 22 × 36 × 10, which corresponds to gridcell sizes of 100 ft × 100 ft × 20 ft. The error decreases dramatically after refining the grid by a factor of ~2, from − 19% to − 4% in oil production rate. A grid convergence study is a crucial process in reservoir simulation as shown in this investigation.
A convergence test has been performed for using the SSG approach as well (Fig. 15). The oil production rate converges very fast. Average reservoir pressure is insensitive to the grid size. Errors decrease as the grid is refined (Table 5). As in  . 15 a Oil production rate and b average reservoir pressure as a function of time at four levels of refinement for Case #2 using SSGs the Case #1, the ~2 times refined-grid results were selected for further analysis. A comparison was performed between SSG and CPG simulations using the second coarsest grid in each case in terms of the oil rate versus. time profile (Fig. 16). Results from these two types of grids are very close to each other. Thus, even with a dipping fault, the gridding effect is small for Case #2. Both SSG and CPG yield accurate results.  Fig. 19 Coarse-scale and fine-scale oil production rate profiles for Case #3 using SSG (ortho) and CPG (non-ortho) The geological model of Case #3 is similar to Case #1 as shown in Fig. 17. Both of these cases have flat horizons (Fig. 17a). The difference is that there is an additional fault in Case #3. One of the faults was made almost parallel to the right-hand side boundary, the other fault has an angle of 45°. The cross section plot shows that the fault throw is small compared with the reservoir thickness (Fig. 17b). Again, two types of grids were generated in the static modeling package, namely SSGs and CPGs (Fig. 18). In CPGs, gridcells are non-orthogonal in a large area around the fault.
The 3D reservoir has the same homogeneous properties, including porosity, permeability, and net-to-gross ratio, as in Case #1 and Case # 2. The size of the reservoir is about 7800 ft × 5250 ft × 100 ft. Well constraints are the same as in Case #1.
A grid convergence study was performed for Case #3 in a similar manner to the ones conducted for Case #1 and Case #2. In Fig. 19, it could be observed that the simulation results for SSG and CPG approaches both converge. However, they converge to different solutions. As observed earlier, SSG results converge to the "true" solution" after refinement. On the other hand, CPG results do not converge to the "true" solution due to relatively more severe non-orthogonality errors. An additional investigation was performed to improve the results from CPG via MPFA. This effort is explained in the section below.

MPFA versus TPFA investigation on Case #3
A comparison of using TPFA and MPFA (O-method) on SSG and CPG for Case #3 is shown in Fig. 20 for the field oil production rate versus. time prediction. MPFA does not make any difference on SSG because both TPFA and MPFA methods operate accurately, as would be expected. While the corrective effect of MPFA on CPG is evident, there is still some difference between SSG and CPG results after MPFA is used. This is because while MPFA largely acts on the total flux (flow) part of the problem, there is also a transport component to multiphase flow and MPFA alone is not sufficient to resolve the mobility discretization errors on non-orthogonal grids. In order to prove this hypothesis, we have designed a single-phase flow variant of the Case #3 problem on the same grids, which eliminates the multiphase transport component of the simulation problem. We have injected oil at the injection well instead of water to reduce the two-phase flow problem into a single-phase flow Fig. 21 Comparison of oil production rate profiles for TPFA and MPFA schemes with faults for modified (single-phase flow) Case #3

Fig. 22
Oil production rate profiles at four levels of refinement for Case #1 including fault seal calculations: a CPGs and b SSGs problem for testing purposes. We have operated the reservoir above bubble point pressure during the entire simulation period to ensure single-phase flow under subsurface conditions. The wells were driven by the same BHP constraints used in Case #1. Oil production rate versus. time profiles for the single-phase flow variant of Case #3 are illustrated in Fig. 21. Once multiphase transport effects are eliminated from the test problem, MPFA prediction on CPG coincides with TPFA and MPFA predictions on SSG while the TPFA prediction on CPG differs significantly from them.

Effect of fault seal modeling
In this section, we study the effect of gridding algorithms in conjunction with a fault seal modeling capability, in which fault permeability and thickness are transferred into transmissibility multipliers (or so called seal factors) at the fault Fig. 23 Oil production rate profiles using SSG and CPG for Case #1 including fault seal calculations

Case #1: model with a single fault and flat horizons
For Case #1 with the fault seal calculations, field oil production rate profiles for CPGs and SSGs individually exhibit converge with grid refinement, as shown in the panels of Fig. 22, and converged profiles agree with each other, as illustrated in Fig. 23. Relative errors are documented in Tables 6 and 7. The second coarsest CPG and SSG yield almost identical results (Fig. 23). The results are also similar to the results obtained for Case #1 without fault seal modeling. Thus, for simple cases, such as Case #1, the grid convergence behavior and gridding effects are quite similar with or without fault seal modeling.

Case #2: model with a single fault and non-flat horizons
For Case #2 including fault seal modeling, the field oil rate versus. time profiles for the CPG approach do not exhibit converge after ~16 times refinement in each of the horizontal directions (Fig. 24a). Relative errors decrease after each refinement. However, the error between the finest-grid result and second-finest-grid result is still − 13% (Table 8), which is relatively large. Convergence behavior of SSG simulations is illustrated in Fig. 24b. The SSG approach exhibits convergence, yet at a slower rate compared to the Case #2 without fault seal modeling. Relative errors for the SSG approach are reported in Table 9. The grid convergence is more difficult because fault seal modeling brings another source of heterogeneity. We tried to further reduce the grid size, however, found out that the static modeling package does not work properly when asked to generate a grid with ~32 times refinement in each of the horizontal directions (with a corresponding gridcell size is about 3 ft in x and y directions). This is the result of the fact that the level of refinement is restricted by static modeling package's maximum number of gridcells or minimum grid sizes. The solutions are slightly different for the ~16-times refined mesh with CPG and SSG approaches (Fig. 25). The difference at the end of the first year is 14%.

Discussion
We use simple conceptual layer-cake models to better isolate and quantify the dynamic effects of gridding for fault modeling purposes. Having stated that, key learnings of our study carry to more complex and heterogeneous reservoir models as well. This is because the investigated simpler models encompass crucial fundamental aspects of fault modeling. We recommend, as much as the model geometry permits, to utilize the SSG approach due its straightforward grid convergence properties. If using the CPG approach is inevitable, special attention must be paid to the grid orthogonality and minimizing the number of highly non-orthogonal grids. An iterative grid construction and quality-checking protocol may be considered to strike an optimal balance between accuracy in fault representation and quality of the resulting grid. We strongly recommend the use of MPFA in conjunction with CPG as it cures part of the grid non-orthogonality related errors, namely the flow part of the coupled flow and transport problem that underlies the simulation of multiphase flow. One can consider a grid-adaptive implementation of MPFA, where the MPFA application is restricted to the relatively more severely non-orthogonal gridcells of a given model (typically corresponding to the vicinity of the fault zone), while the remaining parts of the model with more orthogonal gridcells operate with TPFA. It is important to note that, in this paper, grid convergence results are interpreted from the viewpoint of reservoir performance forecasting purposes. On the other hand, other uses of reservoir models such as well-testing and well-productivity studies focus more on the initial production. As such, more stringent criteria and focusing on the earlier part of the production time profiles are required. For such studies local grid refinement can be considered, which is kept outside the scope of our investigations.

Conclusions
For simple geometries, both the SSG and CPG approaches can be used for fault modeling with similar accuracy in conjunction with the pillar-grid approach. This is evidenced by the fact that simulation results from SSG and CPG approaches converge to identical solutions. SSG and CPG approaches yield different results for more complex geometries. Simulation results for relatively fine SSGs approximate the "true" solution accurately. However, the SSG approach only provides a rough approximation to the geometry and reservoir volumes when the grid is coarse. Non-orthogonality errors occur for the CPG approach and such errors cannot be addressed by grid refinement. MPFA partially addresses the problem of non-orthogonal grids that arise from the CPG approach.
The grid convergence behaviors and gridding effects are quite similar with or without fault seal modeling for simple geometries. However, in a more complex test case examined in this paper, we have observed that it is more difficult to achieve converged results in conjunction with fault seal modeling due to increased heterogeneity of the underlying problem.