A multiscale Robin-coupled implicit method for two-phase flows in high-contrast formations

In the presence of strong heterogeneities, it is well known that the use of explicit schemes for the transport of species in a porous medium suffers from severe restrictions on the time step. This has led to the development of implicit schemes that are increasingly favoured by practitioners for their computational efficiency. The transport equation requires knowledge of the velocity field, which results from an elliptic problem (Darcy problem) that is the most expensive part of the computation. When considering large reservoirs, a cost-effective way of approximating the Darcy problems is using multiscale domain decomposition (MDD) methods. They allow for the pressure and velocity fields to be computed on coarse meshes (large scale), while detailed basis functions are defined locally, usually in parallel, in a much finer grid (small scale). In this work we adopt the Multiscale Robin Coupled Method (MRCM, [Guiraldello, et al., J. Comput. Phys., 355 (2018) pp. 1-21], [Rocha, et al., J. Comput. Phys., (2020) 109316]), which is a generalization of previous MDD methods that allows for great flexibility in the choice of interface spaces. In this article we investigate the combination of the MRCM with implicit transport schemes. A sequentially implicit strategy is proposed, with different trust-region algorithms ensuring the convergence of the transport solver. The method is assessed on several very stringent 2D two-phase problems, demonstrating its stability even for large time steps. It is also shown that the best accuracy is achieved by considering recently introduced non-polynomial interface spaces, since polynomial spaces are not optimal for high-contrast channelized permeability fields.


Introduction
Multiscale domain decomposition methods are a suitable choice to deal with the huge elliptic problems arising from the discretization of the equations governing multiphase flows in oil reservoirs [1]. They allow for the pressure and velocity fields to be computed on a coarse mesh (large scale), while detailed basis functions (locally defined for each subdomain) incorporate rock heterogeneity on a much finer grid (small scale) [2]. The local problems can be solved simultaneously in state-of-the-art parallel machines, making the simulation of huge problems feasible [3].
Several multiscale methods have been presented in the context of the finite volume method [4,5], the finite element method [6,7,8], and mixed finite elements [9,10,11,12,13]. We consider here the Multiscale Robin Coupled Method (MRCM [14,15]) for the solution of two-phase flow problems. It is a generalization of the Multiscale Mixed Method (MuMM [12]) that allows for the independent choice of pressure and flux interface spaces through local Robin boundary conditions. For improved accuracy when the permeability field has high contrast and is channelized, we incorporate adapted interface spaces which behave better than classical polynomials [16,17].
Our focus in this work is the solution of nonlinear two-phase flow models. In the literature the coupling of multiscale flow and transport problems has been treated by explicit operator splitting techniques [18,19] and implicit formulations [20,21]. We have considered operator splitting techniques with explicit approximations for the transport problem in previous works [15,17]. Here, we propose to combine the MRCM with a Sequential Implicit (SI) scheme [22] that allows for the use of large time steps for the coupled flow and transport problem, in contrast to explicit time integration approaches where a CFL-type condition restricts the size of the time step for the transport calculation. The SI algorithm solves at each time step one update of the pressure equation followed by the implicit solution of the transport by the Newton method, considering the velocity field fixed in time.
To ensure convergence of the nonlinear loop we consider trust-region methods to guide the Newton iterations. Specifically, our new solver has three options of trust-region algorithms. The first is the strategy of Jenny et al. [23], where two successive iterations cannot cross any trust-region boundary delineated by the inflection-point of the analytic flux function. The other two algorithms are based on least-square methods, that automatically select the trust-regions to define the iterative step of the implicit solver which are the trust-region reflective and the trust-region dogleg algorithms [24,25,26]. We remark that many improvements have been proposed to the inflection-point strategy [23], for example, the extension to flows with buoyancy and capillary forces [27], compositional simulations [28], and the development of a numerical trust-region solver based on the discretized flux function [29]. Other developments based on trust-region methods are available in the literature. For instance, in [30] the authors developed a flux-search solver for three-phase flow problems.
In this paper we investigate the behavior of domain-decomposition-based multiscale methods (specifically, the family of methods parameterized by the MRCM) when coupled with implicit transport solvers for the simulation of two-phase flows. The approximation is assessed in terms of accuracy and computational efficiency. Additionally, we compare three nonlinear iterative schemes considering different trust-region algorithms so as to determine the one that provides the best performance for the proposed multiscale SI method.
Though satisfactory results in the studied flows were obtained, the semi-implicit treatment of the velocity can generate material balance errors in more complex models (e.g., compositional flows) [31].
The Sequential Fully-Implicit (SFI) scheme, developed in the Multiscale Finite Volume Method (MSFV) framework [20], is an option to deal with this issue. See [32,33,34,35] for extensions to compositional flow simulations, and [36] for developments on nonlinear acceleration techniques. We also study in this work the MRCM combined with the SFI scheme to approximate a challenging two-phase flow problem, with strong fingering instabilities [37].
Summarizing, the main contributions of this work are: • We show that the family of multiscale methods considered here can be efficiently coupled, through both SI and SFI schemes, with implicit transport solvers for the simulation of two-phase flows.
• We identify the trust-region algorithm, with the inflection-point strategy, as the best-performing nonlinear iterative method for the aforementioned coupling.
The rest of the paper is organized as follows: In Section 2, the two-phase flow model is presented. We recall the MRCM in Section 3 and present the sequential implicit formulation in Section 4. Numerical simulation results are presented in Section 5. Finally, our conclusions are presented in Section 6.

Two-phase flows
The governing system of equations for the two-phase problem consists of a second-order elliptic equation for pressure coupled to a hyperbolic conservation law for the saturation of one of the phases [38].
The phases considered are water and oil (denoted by w and o, respectively), and the sum of their saturations is equal to one since we assume a fully saturated porous medium. We consider an immiscible and incompressible two-phase flow in a reservoir containing injection and production wells. For simplicity, capillary pressure effects are not taken into account. However, the proposed method could incorporate these effects by means of an operator splitting algorithm [39]. Adopting a mixed formulation for the pressure equation and a finite volume scheme for the hyperbolic conservation law, the unknowns of the two-phase flow problem are the Darcy velocity upx, tq, the fluid pressure ppx, tq and the water saturation spx, tq. The pressure and velocity are related by Darcy's law so that the elliptic problem can be written where Ω Ă R d , d " 2 or d " 3 is the domain, BΩ " BΩ p Y BΩ u , BΩ p X BΩ u " H; Kpxq is the symmetric, uniformly positive definite absolute permeability tensor; q " qpx, tq is the source term; p b " p b px, tq is the pressure boundary condition at boundary BΩ p ; u b " u b px, tq is the normal velocity boundary condition (n is the outward unit normal) at the boundary BΩ u ; λpsq is the total mobility, given by the sum of the mobilities of the phases: where k rj psq and µ j , j P tw, ou, are, respectively, the relative permeability and viscosity of phase j. The water saturation problem is governed by the transport equation where s 0 ands are, respectively, the initial and injection conditions for the water saturation. Here, BΩ´" tx P BΩ, u¨n ă 0u represents the inlet boundaries. The function f psq is the nonlinear fractional flow of water, given by For simplicity, we assume a constant porosity (scaled out by changing the time variable).
We combine the MRCM to solve equation (1) with an implicit approximation of the hyperbolic conservation law for the water saturation (3) in a sequential fashion, aiming at a compromise between accuracy and computational efficiency of numerical simulations.

The Multiscale Robin Coupled Method
The Multiscale Robin Coupled Method has been introduced to solve elliptic equations accurately, presenting advantages when compared to other existing multiscale mixed methods for flows in highly heterogeneous porous media [14,15]. The MRCM solves equation (1) by a decomposition of the domain Ω into non-overlapping subdomains Ω i , i " 1, 2,¨¨¨, N . Weak continuity of normal fluxes and pressures are imposed to the multiscale solution pu h , p h q at the skeleton Γ of the decomposition (the union of the interfaces Γ ij " Ω i X Ω j ), whose characteristic size H is significantly larger than the fine-scale of the discretization (H " h). The weak continuities are enforced by the following compatibility conditions ż Γ puh´uh q¨ň ψ dΓ " 0 and for all pφ, ψq P U HˆPH , that are the interface spaces defined over the edges E h of the skeleton Γ as subspaces of In Eq. (5), the`and´superscripts represent the solution on each side of the interface Γ, whileň is a fixed normal vector to the skeleton Γ (pointing outwards from the subdomain with the smallest index).
The formulation of the MRCM consists of finding local solutions pu i h , p i h q within each subdomain Ω i , and interface unknowns pU H , P H q satisfying along with the compatibility conditions for all pφ, ψq P U HˆPH , where κ " λpspxqqKpxq, andň i is the normal vector to Γ pointing outwards of Ω i . The parameter for the local Robin boundary conditions can be written as where αpxq is a dimensionless algorithmic function. We remark that the variation of this dimensionless parameter, along with a suitable choice of interface spaces, result in a family of different methods, and by setting this function to extreme values (α Ñ 0 and α Ñ`8) well-known multiscale mixed methods can be recovered [14].
The implementation of the MRCM consists of solving the equations (7) independently to obtain a local set of multiscale basis functions for each subdomain. The global solution is then given by a linear combination of the basis functions and coefficients obtained from the global interface system (8).

Interface spaces
It is well known that the classical polynomials are not optimal for high-contrast channelized permeability fields [16]. Here we recall novel interface spaces based on physics to deal with permeability fields containing highly-permeable channels and barriers, that has been introduced recently [17,40]. These interface spaces are particularly relevant when the channels and barriers are relatively large as happens in karst reservoirs [41,42].
Let Γ i,j Ă Γ be an interface with support in the line segment ra, bs through which N high highpermeability channels pass. Let ra k , b k s Ă ra, bs, k " 1,¨¨¨, N high denote the respective support of each channel. The pressure space contains the following basis functions that mimic the behavior of the pressure solution across the channels: for k " 1,¨¨¨, N high . We assume b 0 " a and a N high`1 " b in Eq. (12). Therefore, the total number of pressure basis functions at Γ i,j is 2`N high .
Next, let Γ i,j Ă Γ be an interface with support in ra, bs through which N low low-permeability structures pass. Let ra k , b k s Ă ra, bs, k " 1,¨¨¨, N low the respective support of each low-permeability structure. The flux space mimics the behavior of the flux across the barriers and contains the following basis functions: for each low-permeability structure k " 1,¨¨¨, N low , where we assume a N low`1 " b. Therefore the total of flux basis functions at Γ i,j is 1`2N low .
We use the MRCM with the physics-based interface spaces to capture the geometry of the highpermeability channels and low-permeability structures at each interface. We also consider the adaptivity in the αpxq function according to permeability variations [15]. Thus, one can control the relative importance of each interface space at each location. As proposed in [15], we take a small value (pressure is favored) at the high-permeability channels and a large value (flux is favored) for the remaining areas.
The interface spaces at the interfaces without high-permeability channels or low-permeability structures are chosen as linear polynomials.

A sequential implicit solver for two-phase subsurface flows
In the presence of strong heterogeneity, explicit schemes for the transport of saturation suffer from severe time-step restrictions. To solve the coupled equations (1) and (3) we consider the sequential implicit method [22], which allows for the use of large time steps, improving the computational efficiency of the simulation.
In the SI algorithm, each time step consists of a sequential update for the flow and transport problems, where a (nonlinear) Newton loop is used to solve the transport equation implicitly. We denote by ∆t the time step used to update the coupled problems of flow and transport at times t n " n∆t, for n " 0, 1, . . . .
Let p n pxq, u n pxq and s n pxq denote the pressure, velocity and saturation approximations for ppx, t n q, upx, t n q and spx, t n q respectively, at time t n . To find the updated variables, one first computes s n`1 pxq (as detailed below) and then solves (1) for the pressure p n`1 pxq and velocity u n`1 pxq keeping the saturation frozen at s n`1 .
The saturation s n`1 pxq is computed through Eq. (3) by using a simple implicit Euler time integration considering u constant in time as follows s n`1´sn ∆t`∇¨`f ps n`1 qu n˘" 0.
Upon a finite volume discretization, the problem for s n`1 reads where I refers to a computational cell of an orthogonal, uniformly spaced (by directions) grid identified as an index (I " pi, jq in 2D and I " pi, j, kq in 3D), V I represents the volume of cell I, and F n`1 I is a function of f ps n`1,ν`1 q and u n , that represents the balance of fluxes at the faces of cell I.
We solve Eq. (17) by variants of Newton's iterative method. Let ν refer to the iteration level of the Newton loop for saturation and set s n`1,0 " s n . In a pure Newton scheme, the next iterate s n`1,ν`1 is defined by the linear system H 1 ps n`1,ν q d ν "´Hps n`1,ν q, where Hps n`1,ν q " H 1 ps n`1,ν q is the Jacobian matrix of H and s n`1,ν`1 " s n`1,ν`dν . The solution at the new time level is achieved when the change in the saturation between two successive iterations is less than a specified tolerance denoted by η. In other words, s n`1 " s n`1,ν`1 if s n`1,ν`1´sn`1,ν ď η.
We consider here a first-order upwind scheme to define F n`1,ν I . In the 2D case we have with discrete fluxes F n`1,ν i´1{2,j and G n`1,ν i,j´1{2 on respective interfaces x i´1{2 and y j´1{2 given by where u x " u x px, yq and u y " u y px, yq denote the x and y components of the velocity field u, and " spx i , y j , t n`1 q represents the saturation (assumed to be a piecewise constant over each computational cell) at time t " t n`1 and at Newton iteration ν.
To ensure the convergence of the nonlinear loop we consider trust-region algorithms instead of pure Newton iterations. In the next subsection we present more details about the algorithms that we have implemented and compared.

Trust-region algorithms
Nonlinearity is a challenging issue for reservoir simulations. Convergence failures for the transport problem are related to the nonlinearity of the flux function. The Newton method is not guaranteed to converge for large time steps, and it can be sensitive to the initial guess [23]. To ensure convergence of the nonlinear loop we consider trust-region algorithms to guide the Newton iterations. Specifically, our solver has three options of trust-region algorithms: the inflection-point strategy of Jenny et al. [23], the trust-region reflective algorithm [24], and the trust-region dogleg algorithm [25].
The inflection-point strategy was introduced to deal with Newton's initial guesses that are on the opposite side of the saturation inflection point with respect to the saturation solution [23]. To ensure the convergence of the Newton iterative process for any time step size, two successive saturation updates are made on the same side of the saturation inflection-point. Hence, if an update would cross the inflection point, selective under-relaxation is applied, i.e., if f 2 ps n`1,ν`1 qf 2 ps n`1,ν q ă 0, then s n`1,ν`1 " ps n`1,ν`1`sn`1,ν q{2. Additionally, it is necessary to enforce the constraint 0 ď s n`1,ν`1 ď 1 after every iteration, that is justified by the physics of the problem [23].
The inflection-point strategy can be seen as a trust-region method that defines different saturation regions delineated by the inflection-point. The updates are performed such that two successive iterations cannot cross any trust-region boundary. A more general trust-region Newton method, that includes saturation trust-regions delineated by the unit-flux and endpoints, is presented in [27]. Another extension of the inflection-point strategy was the development of a numerical trust-region solver, that is based on the discretized flux function [29]. However, for the problem at hand (two-phase flows without gravity and capillary effects) the inflection-point strategy is a particular case of the methods of [27] and [29].
The other two trust-region algorithms considered here are quite popular for nonlinear least-squares problems [43]. They define the iterative update d ν by minimizing a model function in a selected region [44]. To explain this approach, consider the unconstrained minimization problem where x represents the vector with unknowns s n`1 i,j and ϕ : R n Ñ R is the objective function to be minimized, i.e., The trust-region algorithm at iteration ν defines d ν by solving the following sub-problem where m ν is a model function that represents ϕ near the current point x ν and ∆ ν is the trust-region radius, that is adjusted at each iteration to produce a sufficiently decreasing approximation (ϕpx ν`dν q ă ϕpx ν q).
The reflective algorithm uses as model function the quadratic function where B ν is the Hessian matrix ∇ 2 ϕpx ν q or an approximation to it [43]. The minimization problem is restricted to d belonging to the two-dimensional subspace spanned by the gradient direction ∇ϕpx ν q and the Newton direction B ν d "´∇ϕpx ν q [24]. For this algorithm, the Newton system is solved by applying the preconditioned conjugate gradient method [45].
The trust-region dogleg algorithm, on the other hand, adopts [26] m ν pdq " Hps ν q`H 1 ps ν qd The update d is computed as a linear combination of the Cauchy and Newton steps, as presented in [25].
The Cauchy step r d minimizes the model function m ν in Eq. (27) along the steepest descent direction.
The Newton step is the unrestricted minimum of m ν given by H 1 p d "´H (i.e., the update defined in Eq. (18)). The dogleg algorithm chooses d " r d`χp p d´r dq, where χ is the largest value in r0, 1s such that Our solver for the transport problem considers the three trust-region algorithms mentioned above.
We perform comparisons of the approximations provided by them in the section with numerical results.
Additionally, we consider for comparison the Newton method with a global under-relaxation factor of 0.5, in line with [23]. The Newton method with global under-relaxation is stable but requires significantly more iterations to converge when compared to the Newton method with the inflection-point strategy, where the under-relaxation is only applied locally.

Numerical experiments
In this section we present numerical experiments to study the performance of the sequential implicit solver using the MRCM for the approximation of two-phase flows. We compare the saturation approximations provided by the Newton method combined with the trust-region algorithms mentioned in subsection 4.1.
The tolerance for the Newton step size is set to η " 10´6 in the L 2 norm, and the time is expressed in PVI (Pore Volume Injected) [1]. The relative permeabilities are given by k ro " p1´sq 2 and k rw " s 2 , such that the fractional flow of water can be written as This is the configuration considered in the numerical studies unless stated otherwise.
We first compare the convergence of the saturation solution provided by the upwind method in the implicit and explicit versions. Then, we investigate the MRCM combined with the sequential implicit solver for different choices of permeability fields. Finally, we close our numerical experiments with an example that considers gravity effects.

Implicit versus explicit
The Newton method combined with the trust-region algorithms considered is unconditionally convergent, allowing for arbitrary sizes of time steps. Thus, the choice for the size of the time step is based only on accuracy requirements. One can take much larger time steps by using the implicit method instead splitting scheme (see [39] for additional discussion about the operator splitting framework). Here, the explicit approach fixes the same time step size for both elliptic and hyperbolic equations.
We consider a high-contrast permeability field in the domain Ω " r0, 1sˆr0, 1s (with 30ˆ30 fine grid cells). We show in Fig. 1 the log-scaled permeability filed (left), the saturation solution at the final time T PVI " 0.0625 (center), and the relative L 1 pΩq error for saturation as a function of the time step size (right). The time steps considered have size varying from ∆t " 9.75ˆ10´6 to ∆t " 1.25ˆ10´3 (in PVI), while the reference solutions consider ∆t " 10´6 (in PVI, and satisfying the CFL condition). One can note that the errors are essentially the same for the time step sizes smaller than the CFL restriction (∆t ď ∆t CF L « 7.8ˆ10´5). The explicit scheme cannot handle ∆t ą ∆t CF L , while the implicit method maintains the same behavior (linear slope) for the larger sizes of time steps. The average number of required Newton iterations per time step for the cases simulated vary from 2.36 for the smallest time step size to 6.5 for the largest one. Therefore, the implicit solver allows for approximating accurate solutions with large time step sizes, being a good choice to improve the efficiency of two-phase flow simulations.

A Gaussian permeability field
In this example we consider a permeability field given by Kpxq " e 4.5ξpxq , where ξpxq is a self-similar Gaussian distribution having zero mean and covariance function given by Cpx, yq " |x´y|´1 {2 [46]. For this field, the permeability contrast is K max {K min « 10 6 and the computational grid has 64ˆ64 cells in Ω " r0, 1sˆr0, 1s. Concerning the MRCM, we consider linear interface spaces and set αpxq " 1. The domain decomposition considered has 4ˆ4 subdomains, each one containing 16ˆ16 fine grid cells. In order to recover continuous fluxes at the interfaces of the skeleton we use the stitch downscaling procedure, we refer the reader to [47] for details about this scheme. The convergence results for the Newton method for a choice of ∆t " 2ˆ10´2 « 100∆t CF L (that corresponds to a total of 10 time steps) are shown in Fig. 3. In this study, we compare the size of the Newton step provided by the trust-region algorithms. We also consider in the comparison the Newton method with global under-relaxation. The fine grid and MRCM procedures are considered for the elliptic updates. We show the size of the iterative Newton step computed at times T PVI " 0.04, 0.08, 0.12, 0.16, 0.2, i.e., T PVI " 2∆t, 4∆t, 6∆t, 8∆t, 10∆t. We note that the number of iterations needed by the trust-region algorithms are significantly smaller than the required by the underrelaxation technique in all cases. The advantage of the the inflection-point strategy with respect to under-relaxation has been investigated in [23]. Here, we show that the trust-region dogleg and reflective algorithms can also be competitive to approximate the transport problem. Moreover, we find that by using the MRCM to update the velocity field we obtain similar number of Newton iterations.
The behavior of the methods over time is shown in Fig. 4, where we present the number of Newton iterations required in a simulation with three different time step choices. We present results for ∆t chosen as ∆t " ∆t CF L , ∆t " 10∆t CF L and ∆t " 100∆t CF L , that generate, respectively, a total of 1000,  Both in Fig. 4 and Fig. 5 we note that the number of iterations required by the procedures that use the fine grid and MRCM to compute the velocity field is essentially the same. Therefore, the SI solver does not suffer from an increase in the number of iterations when combined with the MRCM.
We perform a convergence study by setting ∆t " 0.1∆t CF L as the reference time step. We estimate the errors of each method for the hyperbolic equation considering each reference solution computed with ∆t " 0.1∆t CF L . In space, we consider the fine grid (with the fine grid velocity field approximation to be the reference) in Fig. 6 (left), the MRCM (with the fine grid velocity field approximation as reference) in Fig. 6 (center), and the MRCM (with the MRCM approximation as reference) in Fig. 6 (right). In the convergence study reported at the center of Fig. 6, the multiscale inaccuracies in the velocity field seem to be relevant when ∆t ă 2ˆ10´2 " 100∆t CF L , while the error of the transport process is dominant for the largest time step choices. We observe linear slope in all methods for the hyperbolic solver when

A channelized permeability field
The next example considers the layer number 36 of the SPE10 project [48], that has a highly-permeable channel and permeability contrast of K max {K min « 10 6 , see Fig. 8

(left). The domain for this example
is Ω " r0, 11{3sˆr0, 1s with 220ˆ60 fine grid cells. For this high-contrast channelized formation, we apply the physics-based interface space for pressure to better represent the solution in the highpermeability channel. We investigate the accuracy of the MRCM combined with the SI approach. The domain decomposition considered contains 11ˆ3 subdomains with 20ˆ20 cells in each one of them.
The flux interface space is linear, as well as the pressure space at the interfaces that do not cross the high-permeability channel. We use the adaptive version of the MRCM [15] by setting αpxq " 10´2 at the interfaces that cross the high-permeability channel and αpxq " 10 2 at the remaining interfaces. Figure 8 (right) shows a map of the absolute permeability variations at the boundaries of the subdomains. The red color identifies the high-permeability channel, where αpxq " 10´2 is set and the physics-based interface spaces for pressure are defined. Figure 9 shows the total accumulated of Newton iterations until time T PVI " 0.11 for different sizes of time step taken as multiples of ∆t CF L « 1.3ˆ10´5. We start with ∆t " 64∆t CF L (that corresponds to a total of 128 time steps) and multiply by four until ∆t " 4096∆t CF L (that corresponds to a total of 2 time steps). We show results for the SI scheme combined with the fine grid velocity field and MRCM.
We    and MRCM is displayed in Fig. 10. The transport problem considers the Newton method using the trust-region algorithm with the inflection-point strategy. We remark that after the convergence of the Newton method all the algorithms considered (the under-relaxation technique, inflection-point, trustregion dogleg, and trust-region reflective) provide similar accuracy. We present saturation maps given by different sizes of time steps. For each choice of ∆t, we note that the approximations that use the fine grid velocity field and MRCM are closely related.
The study of the convergence in time is reported in Fig. 11, where we consider as reference the approximation computed with ∆t " 32∆t CF L . Here, we choose the Newton method using the trust-region algorithm with the inflection-point strategy to report the convergence in time since all the hyperbolic solvers presented the same behavior. Linear slope is attained when the error of each procedure (fine grid velocity field or MRCM) considers its corresponding spatial reference solution. The MRCM errors with respect to the fine grid solution show that the multiscale inaccuracies of the velocity field are relevant when ∆t ă 2.5ˆ10´2 (that corresponds to ∆t ă 512∆t CF L ), while the error of the transport procedure is dominant for larger time step choices. If we compute only the multiscale error, i.e., the error of the MRCM by considering as reference the respective fine grid velocity field with the same time discretization, we obtain an error of the order 3ˆ10´2 for all ∆t, which is consistent with the dominant multiscale error observed on the convergence curve.
Finally, we report in Fig. 12

A homogeneous medium with fingering instability
In this subsection we test the sequential implicit solver for a challenging problem with a fingering instability in a homogeneous medium. In contrast to the previous examples, here we do not consider heterogeneity, however the high nonlinearity of the coupling of flow and transport generates an unstable oil-water interface.
The numerical set-up considers a slab geometry with flow established from left to right by imposing pressure p " 0 on the left and p "´10 4 on the right along with no-flow at top and bottom. The domain Ω " r0, 3sˆr0, 1{2s, with 300ˆ50 fine grid cells, has an initial water front at left, while the rest of the reservoir is filled with oil. The water front contains a small perturbation at the center, as shown in Fig.   13. No source terms are considered. This geometry generates a finger that evolves in time, characterizing a 2D Riemann problem [37]. The nonlinearity, and hence, the physical instabilities of this problem are connected to the viscosity ratio value. Here, we choose M " 4 in line with [37] and [49],   Although the Newton method in the SI solver is unconditionally convergent, the accuracy for this problem with fingering instability drops quickly when ∆t increases. We intend to choose large sizes of time step, aiming at computational efficiency. Therefore, we combine the MRCM with the Sequential Fully-Implicit (SFI) scheme [20], which is adequate to simulate more complex models in which the semi-implicit treatment of the velocity can generate low accuracy [31].

Sequential fully-implicit approximation
The Sequential Fully-Implicit [20] method consists of an outer loop to solve the coupled problems of flow and transport at each time step and the inner (nonlinear) Newton loop to solve the implicit transport problem.
Concerning the SI solver explained in section 4, an additional external loop is added to advance from  previously considered η " 10´6. We remark that if a single iteration of the external loop is performed, the SI algorithm is recovered. Figure 15 shows a comparison of the saturation profiles at time T PVI " 0.6 approximated by the SI and SFI schemes for different sizes of time steps. In this figure, we show results for the procedure that uses the fine grid velocity field. The reference solutions for each case are obtained by taking ∆t " ∆t CF L .
Note that the reference solutions provided by the SI and SFI solvers are similar. For each choice of ∆t ą ∆t CF L , the SFI approximations are more accurate than the SI ones. The SFI scheme presents inaccuracies for choices of ∆t Á 10∆t CF L as well as the SI scheme, however the oil-water interface is better captured by the SFI solver. The corresponding comparison between the SI and SFI schemes combined with the MRCM is shown in Fig. 16. We note that the MRCM works properly when combined with the SFI scheme. The same observations about the relation between the SI and SFI approximations can be drawn if we combine them with the fine grid solution or MRCM.
We show a convergence study for the SI and SFI schemes by setting their respective solutions computed with ∆t " ∆t CF L as references. The fine grid and MRCM are used to compute the velocity field, the latter considering its corresponding reference (in space) and the fine grid reference solution. Essentially the same behavior (linear slope) is observed in all the curves. The SFI solver produces errors slightly lower than the SI scheme. However, the differences between the SI and SFI solutions are better observed in the saturation maps, where we note more detailed information when compared to the global error norms.
From the saturation maps we can conclude that the SFI solver is more accurate, and captures well the fingering instabilities when compared to the SI scheme. In terms of computational efficiency, the SFI is clearly more expensive than the SI (see [50,51,40] that would significantly decrease the cost of the SFI procedure). However, larger time steps can be chosen when using the SFI solver. We show in Fig. 18 the number of iterations required by the Newton method to approximate the 2D Riemann problem for the intermediate size of time step ∆t " 16∆t CF L . The number of Newton iterations for the SI and SFI schemes are presented, the latter being composed of five external iterations. Note that the first external iteration of the SFI requires a similar number of Newton iterations to the SI scheme. The subsequent SFI iterations require a smaller number of Newton iterations: the second requires around 5; the third requires around 3; the fourth requires around 2; the fifth is only necessary at the beginning of the simulation.
These observations are essentially the same for both procedures in the case of the fine grid velocity field Fig. 18 (left) and MRCM Fig. 18 (right). Thus, to attain better accuracy than that provided by the SI scheme, some external iterations in the SFI solver are necessary. For this 2D Riemann problem, a maximum of four external iterations was required when ∆t " ∆t CF L , while a maximum of seven external   iterations was required in the case of ∆t " 16∆t CF L .
To close this discussion, Fig. 19 shows the total accumulated of Newton iterations for approximating the 2D Riemann problem until the final time T PVI " 0.6. We report the total of iterations required by the SI and SFI schemes combined with the fine grid velocity field and MRCM for all the different sizes of time steps previously considered. We note the high cost (in terms of the number of Newton iterations) of the SFI scheme when compared to the SI.
Despite the fact that the SFI scheme is more expensive than the SI, the former presents more accurate solutions. We also test the SFI scheme to approximate the previous examples (in subsections 5.2 and 5.3) and we only noticed very small differences between the SFI and SI approximations (in terms of saturation maps) for choices of ∆t Á 200∆t CF L . The same accuracy for both SFI and SI schemes is observed for the mentioned examples in terms of global errors. This is an expected result since the SI algorithm produces satisfactory approximations in the context of two-phase flows considered [31]. However, for the problem with physics instabilities, we notice the SFI scheme performing better than the SI. We have shown that the MRCM works properly when combined with both SI and SFI schemes. Therefore, the MRCM combined with the SFI algorithm can be further applied to simulate more complex models.

An example with gravity
Our last example considers an application of the sequential implicit solver to approximate a two-phase flow problem with gravity. The elliptic equation (1) when incorporating the gravity effects is given by where g is the gravitational acceleration, h is the height, and λ g is the gravitational mobility, that is given by where ρ w and ρ o are, respectively, the densities of the water and oil. The transport equation (3) in the case with gravity is given by

31)
We solve the transport problem by the Newton method as defined in Eq. (18). Here, the balance of fluxes F n`1,ν i,j at a cell pi, jq is a function of f ps n`1,ν q, λ o ps n`1,ν q, u n , Kpxq, ρ w , ρ o , g, and ∇h. We use the Implicit Hybrid Upwinding (IHU) method [52,53] to define where V n`1,ν i,j represents the viscous part (i.e., the term f psqu in Eq. (31)), while G n`1,ν i,j represents the gravity part (i.e., the term f psqKpxqλ o psqpρ w´ρo qg∇h in Eq. (31)). The viscous and gravity parts are treated separately in the IHU framework.
The upwinding of the viscous term is based on the velocity field, and hence, the balance of the fluxes is the same defined in Eq. (20) with discrete fluxes given by Eqs. (21) and (22). On the other hand, the definition of the gravity term is based on density differences, being the balance of the gravitational fluxes given by where the discrete fluxesG n`1,ν i,j˘1{2 on respective interfaces y j˘1{2 are defined based on the fact that the heavier fluid goes down and the lighter fluid goes up as follows: qλ o ps n`1,ν i,j˘1 q λ w ps n`1,ν i,j q`λ o ps n`1,ν i,j˘1 q pρ w´ρo qg if py j˘1´yj qg ą 0 ∆x K i,j˘1{2 λ w ps n`1,ν i,j˘1 qλ o ps n`1,ν i,j q λ w ps n`1,ν i,j˘1 q`λ o ps n`1,ν i,j q pρ w´ρo qg otherwise (34) where K i,j˘1{2 is the harmonic average of Kpx i , y j q and Kpx i , y j˘1 q. Note that the gravity effect acts along the y direction.
To approximate the problem with gravity we use the extension of the trust-region algorithm with the inflection-point strategy that includes trust-regions delineated by the unit-flux and endpoints [27].
This extension also treats kinks in the Newton paths generated by the gravity term. We use an underrelaxation factor of 0.5 to ensure that the solution update does not extend the new state beyond the trust regions delineated by inflection-points and kinks.
The use of large time steps is essential for computational efficiency, especially for examples with gravitational effects that introduce restrictive CFL condition for stability, given by: By considering the relative permeabilities k ro " p1´sq 2 and k rw " s 2 , and hence, f psq " M s 2 M s 2`p 1´sq 2 , with M " µ o {µ w , we have the following estimate max |F 1 psq| ď max |f 1 psqu|`maxˇˇf 1 psqKpxqλ o psqpρ w´ρo qg∇h`f psqKpxqλ 1 o psqpρ w´ρo qg∇hˇď Such severe restriction highlights the importance of transport implicit methods, that allow for the use of large time steps when compared to explicit time integration approaches.
The geometry considered for this example is a classical quarter of a 5-spot problem [1], where a point source for water injection is placed at the bottom left corner and a production well is placed at the top right corner of a square domain. The boundary conditions are no-flow at all boundaries. We consider the dimensionless form of the two-phase flow problem with gravity as presented in Appendix A. Table 1 shows the specification of the model, and Fig. 20 (left) presents the high-contrast permeability field considered, that is part of the top layer of the SPE10 project [48]. The MRCM approximation for this problem uses the adaptivity of the Robin parameter by setting αpxq " 10´2 at the interfaces that cross highly-permeable regions and αpxq " 10 2 at the remaining interfaces. We denote by aMRCM-PBS  the MRCM version that uses the adaptive Robin parameter and the physics-based interface spaces, and include for comparison the adaptive MRCM with linear interface spaces denoting by aMRCM-POL. In all cases a domain decomposition of 3ˆ3 subdomains with 20ˆ20 cells into each one is considered. Figure   20 (right) shows a map of the absolute permeability variations at the boundaries of the subdomains.
The red color identifies the high-permeability regions, where αpxq " 10´2 is set and the physics-based interface spaces for pressure are defined. The cyan color identifies the low-permeability regions, where the physics-based interface spaces for flux are defined.   The MRCM produces accurate and robust results for the simulation of two-phase flows with gravity when combined with a sequential implicit scheme. Implicit schemes are fundamental to simulate practical reservoirs. The trust-region Newton method considered for the problem with gravity is a significant improvement to the pure Newton method [27]. However, to overcome the severe restrictions on the time step size that still appear, more specialized methods should be used, for example, the numerical trust-region solver proposed in [29].   We note that the simulation converges consistently in the case of the SI scheme. On the other hand, the SFI scheme is quite inefficient for this problem, where the strong coupling between the flow and transport equations causes slow convergence of the outer-loop. One way to deal with such difficulty is presented in [36], where the authors use nonlinear acceleration techniques to improve the convergence of the outer-loop.

Conclusions
The MRCM has been combined with the SI solver for approximating two-phase flows, allowing for the use of large time steps, in contrast to conditionally-stable explicit time integration approaches. The numerical experiments provide strong evidence that the results produced by the MRCM combined with the SI solver are accurate and efficient. Therefore, we can replace fine grid procedures by the MRCM keeping the same parameters of the sequential implicit hyperbolic solvers. This is a promising result in terms of computational efficiency since the MRCM can take advantage of state-of-the-art parallel machines and produce two-phase flow simulations at a reduced computational cost.
To ensure convergence of the nonlinear loop, we have tested the SI solver with different trust-region algorithms. We find that the trust-region reflective and dogleg schemes are appropriate when the size of the time step is chosen of the order of 10∆t CF L , while the inflection-point strategy is adequate to handle sizes of time step of the order of 100∆t CF L to 1000∆t CF L . We have combined the MRCM with an extension of the trust-region algorithm with the inflection-point strategy for solving problems with gravity. In the case with gravity, severe restrictions to the time step appear, making the use of transport implicit methods essential for computational efficiency. The MRCM combined with the SI scheme has shown accurate results for approximating the problem with gravitational effects. Moreover, we have also shown that the best accuracy is achieved by considering the recently introduced physics-based interface spaces for high-contrast channelized permeability fields.
The MRCM has also been combined with the SFI scheme for approximating an example with physical instabilities (fingering). We show that the MRCM works properly when combined with both SI and SFI schemes. The SFI algorithm, as well as nonlinear acceleration techniques, are currently being considered by the authors in order to further apply the MRCM to more complex flow models. Future works also include the implementation of the sequential implicit solver in multi-core machines. Note that each quantity with superscript˚denotes a dimensionless quantity. Now, using these quanti- and the dimensionless form for the transport equation

A.4)
where t˚" u ref L t represents the dimensionless time. However, we report the time variable in PVI (also dimensionless) in our experiments.