Accurate simulation of multiphase flow is essential for predicting pressure buildup and CO₂ plume migration in geological carbon storage. The two-point flux approximation (TPFA) method is widely used in reservoir simulation due to its simplicity, but it is only consistent on K-orthogonal grids. Realistic subsurface models, however, often involve non-K-orthogonal grids caused by geological features such as faults and fractures, where TPFA may introduce significant errors. This study systematically quantifies the numerical error associated with using TPFA on non-K-orthogonal grids in the context of CO₂ storage in saline aquifers. Based on the formalism of the TPFA method, a numerical index is proposed in this work for the first time to quantify the degree of grid non-K-orthogonality. A series of numerical experiments are presented to compare the TPFA and multi-point flux approximation (MPFA) methods, with high-resolution reference solutions used as benchmarks. Results show that grid non-K-orthogonality can impact the accuracy of both TPFA and MPFA solutions with larger degree of grid non-K-orthogonality leading to larger solution errors in general. The magnitude of solution errors of the MPFA method is significantly smaller than that of TPFA in terms of both pressure solution and CO₂ saturation as it is much more robust against grid non-K-orthogonality effect. In particular, the TPFA method can produce substantial deviations in CO₂ saturation distributions when grid non-K-orthogonality is present, indicating the necessity of more advanced discretization methods such as MPFA for modeling the CO₂ plume migration behavior more accurately. These findings highlight the importance of selecting appropriate discretization methods for geologically complex reservoirs and the proposed grid non-K-orthogonality index can help evaluate the solution accuracy in a general simulation study. Our results offer practical guidance on the tradeoffs between computational efficiency and physical accuracy in carbon storage modeling.

Modeling flow and transport in fractured porous media has been a topic of intensive research for a number of energy- and environment-related industries. The presence of multiscale fractures makes it an extremely challenging task to resolve accurately and efficiently the flow dynamics at both the local and global scales. To tackle this challenge, we developed a computational workflow that adopts a two-level hierarchical strategy based on fracture length partitioning. This was achieved by specifying a partition length to split the discrete fracture network (DFN) into small-scale fractures and large-scale fractures. Flow-based numerical upscaling was then employed to homogenize the small-scale fractures and the porous matrix into an equivalent/effective single medium, whereas the large-scale fractures were modeled explicitly. As the effective medium properties can be fully tensorial, the developed hierarchical framework constructed the discrete systems for the explicit fracture–matrix sub-domains using the nonlinear two-point flux approximation (NTPFA) scheme. This led to a significant reduction of grid orientation effects, thus developing a robust, applicable, and field-relevant framework. To assess the efficacy of the proposed hierarchical workflow, several numerical simulations were carried out to systematically analyze the effects of the homogenized explicit cutoff length scale, as well as the fracture length and orientation distributions. The effect of different boundary conditions, namely, the constant pressure drop boundary condition and the linear pressure boundary condition, for the numerical upscaling on the accuracy of the workflow was investigated. The results show that when the partition length is much larger than the characteristic length of the grid block, and when the DFN has a predominant orientation that is often the case in practical simulations, the workflow employing linear pressure boundary conditions for numerical upscaling give closer results to the full-model reference solutions. Our findings shed new light on the development of meaningful computational frameworks for highly fractured, heterogeneous geological media where fractures are present at multiple scales.