Robust solvers for fully coupled transport and flow in saturated/unsaturated porous media

Abstract

Numerical simulations and laboratory studies are our main tools to comprehend better processes happening in the subsurface. The phenomena are modeled thanks to systems of partial differential equations (PDEs), which are extraordinarily complex to solve numerically due to their often highly nonlinear and tightly coupled character. After decades of research on new and improved solving algorithms, there is still a need for accurate and robust schemes. In this work, we investigate linearization schemes and splitting techniques for fully coupled flow and transport in porous media.

A particular case of multiphase flow in porous media, the study of water flow in variably saturated porous media, modeled by the Richards equation, is studied here. An external component, e.g., a surfactant, is transported by the water phases. The resulting system of equations is fully coupled and nonlinear.

In this work, we investigate three different linearization schemes, the classical Newton method, commonly used throughout the industry, the modified Picard method, and the L-scheme. Even though only linearly convergent, the latter appears to be the most robust scheme for some particular cases. The convergence of the L-scheme has also been studied theoretically. Both the Newton method and modified Picard are faster, in terms of numbers of iterations, but fail to converge in cases of unsaturated-saturated porous media or complex phenomena such as hysteresis effects.

The rate of convergence of the L-scheme can be improved by combining it with the Newton method or by using the Anderson acceleration. The scheme resulting by combining the L-scheme and the Newton method appears practically to be both quadratically and globally convergent. It requires fewer iterations than the L-scheme, it is more robust than the Newton method, and it also converges for larger time steps. Using larger time steps can considerably decrease the total number of iterations over the full simulation. Alternatively, by applying the Anderson acceleration, one avoids the computation of any derivatives, and thus the implementation of the algorithm is less invasive. The rate of convergence of the L-scheme depends on user-defined parameters. Finding the optimal L values can be tedious. Optimizing the Anderson acceleration is more straightforward, and the improvements obtained are remarkable.

The equations investigated in this work are not only characterized by nonlinear terms but are also fully coupled. The external component, dissolved into the water phase, directly influences the flow. We investigate two splitting approaches, the canonical nonlinear splitting, and an alternate linearized splitting. We compare them in terms of the numbers of iterations required to achieve the convergence, and the condition numbers of the systems to be solved within each iteration. After all, the latter appears to be a better alternative; it requires fewer iterations and achieves equally accurate results. A monolithic or fully implicit formulation is also investigated. The solution algorithm is computationally slower than the splitting ones but appears to be more robust.

Finally, a global random walk approach for solving the coupled nonlinear problem is also studied. The scheme results in being free of numerical diffusion. Using vast numbers of particles, almost as many as the molecules involved in the reaction, it produces an intuitive representation of the process. The global random walk algorithms are explicit and thus often more straightforward than the typical finite volume/element schemes.

The models investigated in this work represent a particular case of two-phase flow in porous media. Still, we are confident that similar results, concerning the linearization schemes and solving algorithms, can also be achieved in the case of more canonical models.