Shales will play an important role in the successful transition of energy from fossil-based resources to renewables in the coming decades. Aside from being a significant source of low-carbon intensity fuels, like natural gas, they also serve as geologic seals of subsurface formations that may be used to isolate nuclear waste, sequester CO₂, or store intermittent energy (e.g., solar hydrogen). Despite their importance, shales pose significant engineering and environmental challenges due to their nanoporous structure and extreme heterogeneity that spans at least ~10 orders of magnitude in spatial scale. Two challenges inhibit a system-level understanding: (1) the physics of fluid flow and phase behavior in shales are poorly understood due to the dominant molecular interactions between minerals and fluids under confinement, and (2) the apparent lack of scale separation that prevents a reliable (closed) description of the physics at any single scale of observation. In this review, we focus on the latter issue and discuss scale translation, which in its broadest sense is transforming data or simulations from one spatiotemporal scale to another. While effective scale translation is not exclusive to shales, but all geologic porous media, the need for it is especially acute in shales given their high degree of heterogeneity. Classical theories like homogenization, while indispensable, fail when scales are not separated. Other methods, like numerical upscaling, scale-translate in only one direction: small to large, but not the reverse, called downscaling. However, the confluence of advances in three areas are bringing challenging problems such as shales within reach: increased computational power and scalable algorithms; high-resolution imaging and multi-modal data acquisition; and machine learning to process massive amounts of data. While these advances equip geoscientists with a wide array of experimental and computational tools, no individual tool can probe the entire gamut of heterogeneity in shales. Their effective use, therefore, requires an ability to bridge between various data types obtained at different scales. The aim of this review is to present a coherent account of computational and experimental methods that may be used to achieve just that, i.e., to perform scale translation. We provide a broader definition of scale translation, one that transcends classical homogenization and upscaling methods, but is consistent with them and accommodates notions like downscaling and data translation. After a brief introduction to homogenization, we review hybrid methods, numerical upscaling and its recent extensions, multiscale computing, high-resolution imaging, and machine learning. We place particular emphasis on multiscale computing and propose an algorithmic framework to bridge between the pore (micro) and Darcy (macro) scales. Throughout the paper, we draw comparisons between the various methods and highlight their (often hidden) similarities, differences, benefits, and pitfalls. We finally conclude with two case studies on shales that exemplify some of the methods presented.

Enhanced-oil-recovery (EOR) processes involve complex flow, transport, and thermodynamic interactions; as a result, compositional simulation is necessary for accurate representation of the physics. Flow simulation of compositional systems with high-resolution reservoir models is computationally intensive because of the large number of unknowns and the strong nonlinear interactions. Thus, there is a great need for upscaling methods of compositional processes. The complex multiscale interactions between the phase behavior and the heterogeneities lie at the core of the difficulty in constructing consistent upscaling procedures. We use a mass-conservative formulation and introduce upscaled phase-molar-mobility functions for coarse-scale modeling of multiphase flow. These upscaled flow functions account for the subgrid effects caused by the absolute permeability and relative permeability variations, as well as the effects of compressibility. Upscaling of the phase behavior is performed as follows. We assume that instantaneous thermodynamic equilibrium is valid on the fine scale, and we derive coarse-scale equations in which the phase behavior may not necessarily be at equilibrium. The upscaled thermodynamic functions, which represent differences in the component fugacities, are used to account for the nonequilibrium effects on the coarse scale. We demonstrate that the upscaled phase-behavior functions transform the equilibrium phase space on the fine scale to a region of similar shape, but with tilted tie-lines on the coarse space. The numerical framework uses K-values that depend on the orientation of the tie-lines in the new nonequilibrium phase space and the sign of upscaled thermodynamic functions. The proposed methodology is applied to challenging gas-injection problems with large numbers of components and highly heterogeneous permeability fields. The K-value-based coarse-scale operator produces results that are in good agreement with the fine-scale solutions for the quantities of interest, including the component overall compositions and saturation distributions.

We present a reservoir simulation framework for coupled thermal-compositional-mechanics processes. We use finite-volume methods to discretize the mass and energy conservation equations and finite-element methods for the mechanics problem. We use the first-order backward Euler for time. We solve the resulting set of nonlinear algebraic equations using fully implicit (FI) and sequential-implicit (SI) solution schemes. The FI approach is attractive for general-purpose simulation due to its unconditional stability. However, the FI method requires the development of a complex thermo-compositional-mechanics framework for the nonlinear problems of interest, and that includes the construction of the full Jacobian matrix for the coupled multi-physics discrete system of equations. On the other hand, SI-based solution schemes allow for relatively fast development because different simulation modules can be coupled more easily. The challenge with SI schemes is that the nonlinear convergence rate depends strongly on the coupling strength across the physical mechanisms and on the details of the sequential updating strategy across the different physics modules. The flexible automatic differentiation-based framework described here allows for detailed assessment of the robustness and computational efficiency of different coupling schemes for a wide range of multi-physics subsurface problems.

We propose a two-stage preconditioner for accelerating the iterative solution by a Krylov subspace method of Biot’s poroelasticity equations based on a displacement-pressure formulation. The spatial discretization combines a finite element method for mechanics and a finite volume approach for flow. The fully implicit backward Euler scheme is used for time integration. The result is a 2 × 2 block linear system for each timestep. The preconditioning operator is obtained by applying a two-stage scheme. The first stage is a global preconditioner that employs multiscale basis functions to construct coarse-scale coupled systems using a Galerkin projection. This global stage is effective at damping low-frequency error modes associated with long-range coupling of the unknowns. The second stage is a local block-triangular smoothing preconditioner, which is aimed at high-frequency error modes associated with short-range coupling of the variables. Various numerical experiments are used to demonstrate the robustness of the proposed solver.

Miscible gas injection is one of the most effective enhanced oil recovery techniques. There are several challenges in accurately modeling this process, which occurs in the near-miscible region. The adjustment of relative permeability for near-miscible processes is the main focus of this work. The dependence of relative permeability on phase identification can lead to significant complications while simulating near-miscible displacements. We present an analysis of how existing methods incorporate compositional dependence in relative permeability functions. The sensitivity of the different methods to the choice of reference points is presented with guidelines to limit the modification of the relative permeabilities to physically reasonable values. We distinguish between the two objectives of reflecting near-miscible behavior and ensuring smooth transitions across phase changes. We highlight an important link that combines the two objectives in a more general framework. We make use of Gibbs free energy as a compositional indicator in the generalized framework. The new approach was implemented in an automatic differentiation general purpose research simulator and tested on a set of near-miscible gas-injection problems. We show that including compositional dependencies in the relative permeability near the critical point impacts the simulation results with significant improvements in nonlinear convergence.

A novel multiscale method for discrete fracture modeling on unstructured grids (MS-DFM) is developed. To this end, the DFM fine-scale discrete system is constructed using unstructured conforming cells for the matrix with lower-dimensional fracture elements placed at their interfaces. On this unstructured fine grid, MS-DFM imposes independent unstructured coarse grids for the fracture and matrix domains. While the conservative coarse-scale system is solved over these coarse-grid cells, overlapping dual-coarse blocks are also formed in order to provide local supports for the multiscale basis functions. To increase the accuracy, but maintaining the computational efficiency, fracture-matrix coupling is considered only for the basis functions inside the matrix domain. This results in additional (enriching) fracture basis functions in the matrix. By construction, basis functions form the partition of unity for both fracture and matrix sub-domains. Furthermore, to enable error reduction to any desired level, a convergent iterative strategy is developed, where MS-DFM is employed along with a fine-scale smoother in order to resolve low- and high-frequency modes in the error. The performance of MS-DFM is assessed for several 2D and 3D test cases. The proposed method achieves accurate results for several test cases even without iterations, and for challenging ones with only a few iterations. MS-DFM is the first of its kind, and thus extends the application of multiscale methods to unstructured discrete fracture models. As such, it provides a promising framework for real-field application of unstructured DFM.

We present a new fully-implicit, mixed-hybrid, finite-element (MHFE) discretization scheme for general-purpose compositional reservoir simulation. The locally conservative scheme solves the coupled momentum and mass balance equations simultaneously, and the fluid system is modeled using a cubic equation-of-state. We introduce a new conservative flux approach for the mass balance equations for this fully-implicit approach. We discuss the nonlinear solution procedure for the proposed approach, and we present extensive numerical tests to demonstrate the convergence and accuracy of the MHFE method using tetrahedral elements. We also compare the method to other advanced discretization schemes for unstructured meshes and tensor permeability. Finally, we illustrate the applicability and robustness of the method for highly heterogeneous reservoirs with unstructured grids.

Enhanced Oil Recovery (EOR) processes involve complex flow, transport, and thermodynamic interactions; as a result, compositional simulation is necessary for accurate representation of the physics. Flow simulation of compositional systems with high-resolution reservoir models is computationally intensive because of the large number of unknowns and the strong nonlinear interactions. Thus, there is a great need for upscaling methods of compositional processes. The complex multiscale interactions between the phase behavior and the heterogeneities lie at the core of the difficulty in constructing consistent upscaling procedures. We employ a mass-conservative formulation and introduce upscaled phase molar mobility functions for coarse- scale modeling of multiphase flow. These upscaled flow functions account for the sub-grid effects due to the absolute and relative permeabilities variations, as well as the effects of compressibility. Upscaling of the phase behavior is performed as follows. We assume that instantaneous thermodynamic equilibrium is valid on the fine scale, and we derive coarse-scale equations, in which the phase behavior may not necessarily be at equilibrium. The upscaled thermodynamic functions, which represent differences in the component fugacities, are employed to account for the non-equilibrium effects on the coarse scale. We demonstrate that the upscaled phase-behavior functions transform the equilibrium phase space on the fine scale to a region of similar shape, but with tilted tie- lines on the coarse space. The numerical framework uses K-values that depend on the orientation of the tie-lines in the new non-equilibrium phase space and the sign of upscaled thermodynamic functions. The proposed methodology is applied to challenging gas-injection problems with large numbers of components and highly heterogeneous permeability fields. The K-value based coarse-scale operator produces results that are in good agreement with the fine-scale solutions for the quantities of interest, including the component overall compositions and saturation distributions.

A multiscale method for Discrete Fracture Modeling (DFM) using unstructured grids is developed. The fine-scale discrete system is obtained by imposing tetrahedron (triangular for 2D domains) shaped grid cells, while lower-dimensional fractures are imposed at the grid interfaces. The DFM approach is then used to describe the transmissibility coefficients for all the interfaces, including those with the lower-dimensional fractures. On this fine-scale discrete system, a new algebraic multiscale formulation is developed, which first imposes two sets of coarseand dual-coarse grids. The former grid is essential for conservative multiscale formulations, while the second one is used for the calculation of local multiscale basis functions. The coarse-scale partitioning of the fracture and matrix domains is flexible and totally independent. Moreover, the multiscale basis functions are constructed for both the matrix and fracture domains. By construction, the basis functions are a partition of unity. For 2D and 3D test cases, the performance of the multiscale method is systematically assessed. It is shown that the method (with no multiscale iterations) provides accurate results, even for complex fractured systems. The presented multiscale method is a promising framework for real-field application of DFM models.

We present a fully implicit mixed hybrid finite-element (FE) formulation for general-purpose compositional reservoir simulation. The formulation is locally conservative, and the momentum and mass balance equations are solved simultaneously; including Lagrange multipliers on element interfaces. The method utilizes automatic differentiation for the Jacobian construction. This hybrid FE approach accommodates unstructured grids, and we present black-oil and compositional test cases with permeability tensors. We also discuss the accuracy and computational efficiency for the new formulation. For all tests, we compare the performance and accuracy of the proposed approach with the Multi-Point Flux Approximation (MPFA-O) method.

We present a new framework for solving coupled multi-physics problems. The objective is to develop a platform where different coupling strategies for the simulation of complex physical processes can be employed with great flexibility in order to find an optimal - in terms of robustness and computational efficiency - strategy for a given problem. The new simulator is modular; each module represents a particular physics process, such as compositional, thermal, poromechanics, reactions, wells, and surface facilities. The platform provides seamless coupling between the physics modules without resorting to conditional branches and intermediate interface and treats terms that are coupled across multiple physics modules efficiently. The different modules can be coupled with each other in a sequential or fully coupled manner, and different solution strategies can be applied to different modules. This allows investigation of complex coupling strategies that have not been studied before. Examples of target problems include modeling compositional-thermal EOR processes in both conventional and unconventional resources with tight coupling to nonlinear poromechanics. The paper addresses the design of the framework and provides details of its implementation.

A multiscale method for Discrete Fracture Modeling (DFM) using unstructured grids is developed. The fine-scale discrete system is obtained by imposing tetrahedron (triangular for 2D domains) shaped grid cells, while lower-dimensional fractures are imposed at the grid interfaces. The DFM approach is then used to describe the transmissibility coefficients for all the interfaces, including those with the lower-dimensional fractures. On this fine-scale discrete system, a new algebraic multiscale formulation is developed, which first imposes two sets of coarseand dual-coarse grids. The former grid is essential for conservative multiscale formulations, while the second one is used for the calculation of local multiscale basis functions. The coarse-scale partitioning of the fracture and matrix domains is flexible and totally independent. Moreover, the multiscale basis functions are constructed for both the matrix and fracture domains. By construction, the basis functions are a partition of unity. For 2D and 3D test cases, the performance of the multiscale method is systematically assessed. It is shown that the method (with no multiscale iterations) provides accurate results, even for complex fractured systems. The presented multiscale method is a promising framework for real-field application of DFM models.

In this paper, we propose a strategy to bypass the phase identification of fluid mixtures that can form three, or more, phases. The strategy is used for reservoir simulation of multicomponent, three-phase, thermal compositional displacement processes. Since the solution path in compositional space is determined by a limited number of “key” tie-simplexes, the proposed “bypass” method uses information from the parameterized tie-simplexes and their extensions. The tie-simplex parameterization is performed in the discrete phase-fraction space. Once the phase-fraction space is discretized, a conventional three-phase flash is used adaptively to compute the phase states at the discretization nodes. If all discretization vertices of a given discrete cell, in phase-fraction space, have the same phase state, then this state is assigned to the entire cell and expensive flash calculations are bypassed. We demonstrate the robustness and efficiency of our phase identification bypassing strategy for several cases with three-phase flow, including a six-component ES-SAGD (enhanced solvent SAGD) model.

This paper presents a detailed numerical analysis of two methods for phase-behavior computations associated with thermal compositional reservoir simulation. Specifically, we analyze and compare the standard K-values approach with an Equation of State (EoS) model. We study steam and steam–solvent injection into a one-dimensional reservoir saturated with heavy oil in order to identify the impact of using different methods to describe the phase behavior. We then present results for complex displacement processes in multidimensional reservoir models. The investigation shows the differences in the predictions between K-values- and EoS-based simulations can vary from small to quite significant depending on the compositional interactions between the injected and resident fluids in the modeled displacement process.

Miscible gas injection is one of the most effective enhanced oil recovery (EOR) techniques. There are several challenges in accurately modeling this process that mostly occur in the near-miscible region. The adjustment of relative permeability for near-miscible processes is the main focus of this work. The dependence of relative permeability on phase identification can lead to significant complications while simulating near-miscible displacements. We present an analysis of how existing methods incorporate compositional dependence in relative permeability functions. The sensitivity of the different methods to the choice of reference points is presented with possible guidelines to limit the modification of the relative permeabilities to physically reasonable values. We distinguish between the objectives of reflecting near miscible behavior and ensuring smooth transitions across phase changes in the existing methods. We highlight an important link that combines the two objectives in a more general framework. We make use of Gibbs free energy as a compositional indicator to honor the generalized framework. The new approach was implemented in the Automatic Differentiation General Purpose Research Simulator (ADGPRS) and tested on a set of near-miscible gas injection problems. We show that including compositional dependencies in the relative permeability near the critical point impacts the simulation results with significant improvements in nonlinear convergence

Geological storage of carbon dioxide (CO₂) is a promising technology for reducing atmospheric emissions. The large discrepancy in the time- and length-scales between up-dip migration of buoyant supercritical CO₂ and the sinking fingers of dissolved CO₂ poses a challenge for numerical simulations aimed at describing the fate of the plume. Hence, several investigators have suggested methods to simplify the problem, but to date there has been no reference solution with which these simplified models can be compared. We investigate the full problem of Darcy-based two-phase flow with gravity-current propagation and miscible convective mixing, using high-resolution numerical simulations. We build on recent developments of the Automatic Differentiation - General Purpose Research Simulator (AD-GPRS) at Stanford. The results show a CO₂ plume that travels for 5000 years reaching a final distance of 14 km up-dip from the injection site. It takes another 2000 years before the CO₂ is completely trapped as residual (40%) and dissolved (60%) CO₂. Dissolution causes a significant reduction of the plume speed. While fingers of dissolved CO₂ appear under the propagating gravity current, the resident brine does not become fully saturated with CO₂ anywhere under the plume. The overall mass transfer of CO₂ into the brine under the plume remains practically constant for several thousands of years. These results can be used as a benchmark for verification, or improvements, of simplified (reduced-dimensionality, upscaled) models. Our results indicate that simplified models need to account for: (i) reduced dissolution due to interaction with the plume, and (ii) gradual reduction of the local dissolution rate after the fingers begin to interact with the bottom of the aquifer.

Steam Assistant Gravity Drainage (SAGD) is widely used to recover heavy oil and bitumen reservoirs. Typical SAGD operations involve a pair of horizontal wells separated vertically. Steam, or a steam-solvent mixture (e.g., Expanding- Solvent SAGD), is injected into the upper well to form a steam chamber and mobilize oil, which drains to the lower production well. Significant mechanical stresses associated with this process can increase the risk of fracturing the reservoir, or cap-rock. We perform a fully-coupled thermal-compositional-mechanical numerical simulation of SAGD and ES-SAGD processes for a typical bitumen reservoir in the Fort McMurray region of Alberta, Canada. A mixed finite-volume approximation for the flow and a Galerkin finite-element approximation for the mechanics are used, and the resulting set of nonlinear equations are solved using a fully implicit formulation. The two discretizations share the same unstructured grid. We demonstrate that thermo-mechanical effects can be quite significant in SAGD operations. The sharper gradients associated with the standard SAGD process increase the risk of damaging of the cap-rock. On the other hand, ES-SAGD operations lead to more dispersed temperature and pressure distributions, which decreases the possibility of damaging the cap-rock.

In this paper, we present a novel strategy for phase-state identification that can be used to bypass the need for full Equation-of-State Computations in multicomponent, multiphase thermal-compositional displacement processes. Analysis based on the Method of Characteristics (MOC) indicates that the displacement path in compositional space is determined by a limited number of tie-simplexes. Our "bypass" method uses information from the parameterized extensions of these "key" tie-simplexes. The parameterization is performed in the discrete phase-fraction space. In this space, the phase fractions can be negative, or greater than unity; this allows for parameterization of the hyperplane that corresponds to a tiesimplex extension. Once the tie-simplexes extension is discretized, a conventional three-phase flash is used adaptively to compute the phase states at the discrete points in compositional space. If all discretization nodes for a given discrete cell that lie on the hyperplane have the same phase state, this state is assigned to the entire computational cell, and flash calculations are bypassed for the compositions of that cell. Here, we make use of the continuity of the tie-simplex parameterization, which was proven in earlier work. We demonstrate the efficiency and robustness of our "bypass'" strategy for the simulation of flow and transport in thermal, three-phase compositional models of heterogeneous reservoirs.