An accurate prediction of transport in sub-surface formations is challenging for Eulerian schemes, especially if advection dominates. Lagrangian particle-tracking schemes, not being hindered by numerical dispersion, are compelling alternatives and have increased relevance in the modelling of non-linear transport, e.g., saturation transport in a multi-phase setting [1].

In this work, we have focused on the development of a particle-tracking scheme for multiphase flows in fractured media with a permeable matrix. To this end, we adopted an Embedded Discrete Fracture Model where fractures are treated as lower dimensional manifolds [2]. The fracture-matrix interfaces are not resolved at the sub-grid level, thereby necessitating an alternative modeling strategy for inter-continuum interactions. In [3], we presented a stochastic particle-tracking scheme for advective solute transport in single-phase flow. A particle’s continuum state is modeled as a two-state Markov chain with transition probabilities for inter-continuum particle transfer. The probabilities are pathline-specific and scale with the particle’s modeled travel time through the grid cell.

In a first step, we aim to improve the efficiency of the above-mentioned scheme, especially for the scenarios involving high contrast in matrix-fracture pore-volumes. Then, the scheme is extended to model saturation evolution in two-phase immiscible flows. Following the work of [1], saturation is modeled as a statistical quantity and is estimated with two distinct particle ensembles, i.e., one for each phase.

In the absence of dispersive effects, e.g., in flows with high Péclet numbers and without capillary pressure differences, saturation discontinuities are a typical feature. Therefore, near the fronts, inaccurate estimates and instabilities may result due to finite-sized particle ensembles in cells of the computational domain. To render such inaccuracies and instabilities insignificant, an adaptive diffusivity is added to the system, which selectively acts near the front and attains negligible values away from it. To quantify the diffusivity, a Smagorinsky-type [4, 5] model is proposed, which scales with the magnitude of the saturation gradient.

In the pursuit of efficient yet accurate large-scale modeling, the proposed Lagrangian scheme can be extended to connect dynamic and locally unresolved sub-grid processes, e.g., phase dissolution, with the macroscopically observed effects.

References:

1. Tyagi, M., Jenny, P., Lunati, I. and Tchelepi, H.A. (2008). JCP, 227(13): 6696-6714.
2. Deb R. and Jenny P. (2017), IJNAMG 41, 18: 1922-1942.
3. Monga, R., Deb, R., Meyer, D.W. and Jenny, P. (2020), In ECMOR XVII ,EAGE.
4. Smagorinsky, J. (1963). MWR, 91(3): 99-164.
5. Lilly, D.K. (1966). NCAR manuscript: 123.

Governing equations for multicomponent porous media flow typically exhibit strong and unbalanced nonlinearities and orders of magnitude variations in time constants. Together with the complexity of field-scale geomodels, with orders-of-magnitude variations in rock properties and cells with high aspect ratios and contorted geometries, this requires implicit discretizations with robust, global nonlinear solvers. We discuss a nonlinear domain decomposition preconditioning method for fully implicit simulations based on an additive Schwarz preconditioned exact Newton method (ASPEN) that is applicable to multicomponent porous media flow simulation. The method efficiently accelerates nonlinear convergence by resolving unbalanced nonlinearities in a local stage and long-range interactions in a global stage. In previous research, we have shown that use of ASPEN can improve robustness and (significantly) reduce the number of global iterations compared with standard Newton. However, each global iteration is more expensive because of the extra work introduced in the local stage.

Herein, we focus on implementation aspects for the local and global stage, including linear and nonlinear solver settings, and scalability. We show how the global-stage Jacobian can be transformed to have the same structure as the fully implicit system, enabling use of highly efficient linear preconditioners, and discuss how the global-stage update can be recast as a correction to the local-stage solutions.

We compare the computational performance of ASPEN to standard Newton on a series of test cases, ranging from conceptual cases with simplified geometry or flow physics to cases representative of real assets. Our overall conclusion is that ASPEN is outperformed by Newton when this method works well and converges in a few iterations. On the other hand, ASPEN avoids time-step cuts and has significantly lower runtimes in time steps where Newton struggles. A good approach to computational speedup is therefore to adaptively switch between Newton and ASPEN throughout a simulation. A few examples of such strategies are outlined.

At geo-storage conditions, carbon-dioxide is a supercritical low-viscosity fluid that is buoyant with a mobility ratio >> 1, highly prone to unstable displacement.
The best candidate storage formations consist of highly permeable and porous fluvial – intertidal and deltaic siliciclastics with prominent bedforms and laminations. The impact of these heterogeneities on multiphase flow and trapping already gained considerable attention (e.g., Pickup & Sorbie 1996; Trevisan et al. 2017; Ringrose and Bentley 2021). The USGS (Ruben & Carter 2005) created a software tool that permits realistic geometric modelling of cross bedding. Output quadrilateral surfaces are converted into a boundary representation (BREP) with a water-tight subdivision into distinct rock types contributing to representative elementary volumes. The resulting heterogeneous sandstone models are periodic, and are meshed with tetrahedra, providing the flow grids for the present analysis aimed at determining dynamic relative permeability using a full physics approach, accounting for viscous, gravitational, and capillary forces.
To accurately model complex coarse-fine interface flow dynamics, extra degrees of freedom / discontinuities are introduced into the flow models, employing the new nonlinear interface transfer algorithm of Tran et al. (2020). This resolves the formation of dynamic capillary pressure and saturation discontinuities.
Our results reveal that 1) at the >10-cm scale, gravity influences dynamic barrier formation, accelerating capillary sealing. 2) The latter can choke the flow at small fluid supply rates. 3) Increasing rate will lead to episodic drainage-imbibition cycling or continuous flow. 4) Ensemble relative permeability is flow direction dependent and a function of flow rate. 5) saturation distributions are strongly influenced by bedforms. A new method is presented to extract ensemble relative permeability curves for cross-bedded sandstones. It does not rely on uniform inlet fractional flows raising the physical realism of this multiphase flow upscaling.


Cited References

Pickup, G.E. and Sorbie, K.S., 1996. The scaleup of two-phase flow in porous media using phase permeability tensors. SPE Journal, 1(04), 369-382.

Ringrose, P. and Bentley, M., 2021. Models for Storage. In Reservoir Model Design, Springer, Cham.

Rubin, D.M., Carter, C. L. 2005, Bedforms and Cross-Bedding in Animation, USGS.

Tran, L.K., Kim, J.C. and Matthäi, S.K., 2020. Simulation of two-phase flow in porous media with sharp material discontinuities. AWR 142, p. 103636.

Trevisan, L., Krishnamurthy, P.G. and Meckel, T.A., 2017. Impact of 3D capillary heterogeneity and bedform architecture at the sub-meter scale on CO2 saturation for buoyant flow in clastic aquifers. International Journal of Greenhouse Gas Control, 56, 237-249.