Application of a conservative, generalized multiscale finite element method to flow models

We consider a heterogeneous oil reservoir which is confined in a domain $\Omega$. The reservoir is equipped with an injection well, from which water is discharged to displace the trapped oil towards the production wells, situated on the perimeter of the domain. The dynamics of the movement of the fluids in the reservoir are categorized as an immiscible two-phase system with water and oil (denoted by $w$ and $o$, respectively) that is incompressible. Capillary pressure is not included in the model. Further simplifying assumptions that we use are a gravity-free environment and that the two fluids fill the pore space. Then, the Darcy’s law combined with a statement of conservation of mass allow us to write the governing equations of the flow as

and

where $\boldsymbol{v}$ is the Darcy velocity, $S$ is the water saturation, and $k$ is the permeability coefficient. The total mobility $\lambda(S)$ and the flux function $f(S)$ are respectively given by:

where $k_{rj}$, $j=w,o$, is the relative permeability of the phase $j$.

Notice that the elliptic part (2.1) and the transport part (2.2) of the system are coupled through the total mobility. In order to solve this problem numerically, we use an operator splitting technique [3], where saturation at the previous time step is used when solving the elliptic part of the system to obtain the velocity $\boldsymbol{v}$. This velocity is then used in an explicit time stepping scheme for the transport equation. This velocity is held constant for a predetermined number of time steps, which yields a new saturation. This new saturation is then used to solve the elliptic problem again and the process is continued until the final time is reached. A schematic of the operator splitting is shown in Fig. 1.

To discretize the transport equation, we first integrate (2.2) with respect to time and over some $C_{z}\subset\Omega$. Here we apply the left end point quadrature rule to its second term in time and use integration by parts to arrive at the approximation

where we have neglected the error terms and

The transport part of the system is solved with an explicit time stepping as seen above and an upwinding scheme is used on the term $\int_{\partial C_{z}}{\boldsymbol{v}}\cdot\boldsymbol{n}f(S)\,\text{d}l$. A review of upwinding on a rectangular mesh can be found for example in [25]. Typically in this situation, it is imperative that numerical approximation of ${\boldsymbol{v}}$ satisfies certain local conservation property. In particular, given ${\boldsymbol{v}}_{h}\approx{\boldsymbol{v}}$, it is desirable to have

A natural way of obtaining this local conservation property is to seek the solution of $(p,{\boldsymbol{v}})$ from the first order system (2.1). At the approximation stage, one ends up with the mixed finite element formulation [23]. A more common approach is to transform (2.1) into a second order equation that governs $p$. The approximate solution for $p$ is sought after which the approximate solution of ${\boldsymbol{v}}$ is calculated using the relation $\boldsymbol{v}=-\lambda(S)k(\boldsymbol{x})\nabla p$, and hence a postprocessing procedure is used. Unfortunately, a standard technique such as the Galerkin finite element method does not allow for a straightforward postprocessing to obtain a ${\boldsymbol{v}}_{h}$ that is locally conservative.

Furthermore, aside from the issues pertaining to local conservation, it is almost always impossible to conduct numerical simulations at the finer scales if one is to include the ever increasing details of geological information. A viable alternative is the use of multiscale methods, which is the subject of the next section.

The multiscale finite element method (MsFEM) was introduced in [17] and further analyzed in [18]. Similar to the approximation method described earlier, MsFEM is a continuous Galerkin finite element method that is based on solving (3.2) in a finite dimensional subspace of $H^{1}_{\text{D}}$. Its distinct feature is in the careful choice of a multiscale finite dimensional subspace that allows calculation of the approximate solution on $\mathcal{T}_{h}$ without directly resolving the fine scale heterogeneity globally. This means $h$ is much larger than the characteristic fine scale of $k({\boldsymbol{x}})$. Since much of the fine scale heterogeneity has its source from $k({\boldsymbol{x}})$, this information should be ingrained in this subspace. This is done through the construction of the multiscale basis functions that characterize the subspace.

As in the standard Galerkin finite element method, the multiscale basis functions are associated with the nodes/vertices in the discretization of $\Omega$. In the interest of offering a straightforward presentation, in what follows we assume that $\mathcal{T}_{h}$ is a collection of rectangles $\tau$ such as depicted in Fig. 2. Extension of the method to triangles follows the same line of arguments.

For a vertex $z\in Z_{\text{in}}\cup Z_{\text{n}}$, the corresponding multiscale basis functions $\chi_{z}$ are defined in such a way that $\chi_{z,\tau}=\chi_{z}\big{|}_{\tau}$ is governed by

where $\phi_{z,\tau}({\boldsymbol{x}})$ is the standard bilinear function in $\tau$ and $z$ is a vertex of $\tau$. The multiscale finite dimensional space is defined as

The MsFEM solution is $p_{\text{ms},h}$ with $(p_{\text{ms},h}-p_{\text{D},h})\in\mathcal{V}_{\text{ms},h}$ satisfying

Similar to the standard Galerkin finite element method, the linear system resulting from (3.7) is

where $\mathbf{A}_{\text{ms}}$ is a square matrix with entries $a(\chi_{\zeta},\chi_{z})$, $\mathbf{f}_{\text{ms}}$ is a vector with entries $(f,\chi_{\zeta})$ for $\zeta\in Z_{\text{in}}$ and $(f,\chi_{\zeta})+\langle g_{\text{N}},\chi_{\zeta}\rangle$ for $\zeta\in Z_{\text{n}}$, and $\mathbf{p}$ is as in (3.4).

To this end, we emphasize on the significant advantage of MsFEM. If one is to solve the pressure equation (3.1) to the extent that the fine scale heterogeneity of $k({\boldsymbol{x}})$ is directly resolved in the finite element formulation (3.3), then ideally the level of mesh resolution on which the flow and transport are simulated has to be in a comparable scale with the level of subsurface resolution exhibited by $k({\boldsymbol{x}})$. In turn, the resulting linear system such as expressed in (3.4) can have a very large dimension whose inversion poses a very challenging task. Employing MsFEM on the other hand, avoids this drawback as there is no need to pose (3.7) on a mesh having comparable size to the characteristic scale of $k({\boldsymbol{x}})$. This is because the multiscale finite element space in (3.6) already contains this information.

Still, potentially there can be a significant error associated with MsFEM approximation that stems from the imposition of linear boundary condition on $\chi_{z}$ on $\partial\tau$, which causes a mismatch as compared to the true solution. This mismatch is even more pronounced when $k({\boldsymbol{x}})$ exhibits channelized features and/or scattered inclusions in which the the values of $k({\boldsymbol{x}})$ are orders of magnitude higher (or lower) compared to the neighboring regions. Within the context of modeling and simulation of multiphase flow and transport, there have been several research directions aiming toward alleviating this situation. For example, oversampling techniques [5], adoption of global information [11], and the use of a local-global iterative approach [8] have been used for error reduction. More recent investigations involve the alternative of systematically enriching the original coarse space $\mathcal{V}_{\text{ms},h}$. This is the subject of next subsection.

The Generalized Multiscale Finite Element Method (GMsFEM) is based on a systematic enrichment of $\mathcal{V}_{\text{ms},h}$ [10, 14, 15]. This enrichment is made available by taking advantage of the knowledge of the spectral properties of the original differential operator governing the multiscale basis functions $\chi_{z}$; namely, the left hand side of (3.5). These types of enriched spaces yield pressure solutions whose errors decrease with respect to the localized eigenvalue behavior. In related work, it is shown that the errors typically depend on the first eigenvalue that is not included in the space construction. We refer the interested reader to [10] for rigorous error estimates.

The GMsFEM solution is $p_{\text{ems},h}$ with $(p_{\text{ems},h}-p_{\text{D},h})\in\mathcal{V}_{\text{ems},h}$ satisfying

Upon comparison, it is obvious that $\mathcal{V}_{\text{ms},h}\subset\mathcal{V}_{\text{ems},h}$ and thus $\dim(\mathcal{V}_{\text{ems},h})>\dim(\mathcal{V}_{\text{ms},h})$. Consequently, the resulting linear system $\mathbf{A}_{\text{ems}}\mathbf{p}=\mathbf{f}_{\text{ems}}$ has a larger dimension than its counterpart in (3.8), where here we view $\mathbf{p}$ as a vector containing the coordinates (rather than the pressure nodal values) in the linear combination of $\Phi_{z,l}$ that span the approximate pressure. However, this system is still significantly smaller compared to the fine scale analogue in (3.4) that is solved on mesh that has size of the order of the characteristic scale of $k({\boldsymbol{x}})$. Thus, the enriched coarse system offers a suitable reduced-order alternative for obtaining approximate pressure solutions while maintaining an acceptable level of accuracy.

We reiterate that the construction of the enriched basis functions in Eq. (3.12) is performed on a respective $\omega_{z}$ (refer back to Fig. 2). However, the postprocessing technique offered in the next section is localized to $\tau$. As such, it is important to note that the enriched basis functions need to be restricted into respective element contributions to implement the proposed procedure. So while we refer to the enriched basis functions synonymously (whether they are posed on a $\omega_{z}$ or $\tau$), we make this distinction for additional clarity.

Since we desire independent element-by-element calculations, we proceed with a formulation of an auxiliary boundary value problem in $\tau$ (see the right plots of Fig. 3) having vertices $v(\tau)=\{w,x,y,z\}$. The boundary value problem governs $\widetilde{p}_{\tau}$ with

Here we designate $\partial\tau=\cup_{\zeta\in v(\tau)}E^{\tau}_{\zeta}$, where $E^{\tau}_{\zeta}=\partial\tau\cap\partial t_{\zeta}$ (i.e. half of each element edge containing the vertex $\zeta$). Furthermore, we set $\widetilde{g}_{\tau}$ as piecewise function on $\partial\tau$ such that

where

The existence and uniqueness of the above problem is stated below

We note that (4.8) implies

which shows that the solution of (4.5) recovers the flux of $p$ averaged over $\partial\tau$, i.e., a local conservation property. We use (4.5) as a governing principle to derive the postprocessing technique for calculating a locally conservative flux from $p_{\text{ems},h}$.

This elemental calculation is based on discretization of $\tau$ into quadrilaterals $t_{\zeta}$, i.e., $\tau=\cup_{\zeta\in v(\tau)}t_{\zeta}$, each of which yields $t_{\zeta}=C_{\zeta}\cap\tau$, see the right plots of Fig. 3. We set the local solution space as $\mathcal{V}(\tau)=\text{span}\{\Phi_{\zeta,1}\}_{\zeta\in v(\tau)}$. The numerical solution associated with (4.5) is to find $\widetilde{p}_{\tau,h}\in\mathcal{V}(\tau)$ satisfying

The following four equations result from (4.10):

where

and $E^{\tau}_{\zeta\eta}=\partial t_{\zeta}\cap\partial t_{\eta}$, for $\zeta,\eta=w,x,y,z$, and $\zeta\neq\eta$. We note that similar equations will be created if $\tau$ is triangle (see [4]).

Because $\widetilde{p}_{\tau,h}=\sum_{\zeta\in v(\tau)}\alpha_{\zeta}\Phi_{\zeta,1}$ with unknown $\alpha_{\zeta}$, the above equation yields a linear system of the form $\widetilde{\mathbf{A}}\boldsymbol{\widetilde{\alpha}}=\widetilde{\mathbf{f}}$ where

Here the system is of dimension 4. We note that in the case of $\tau$ adjacent to $\Gamma_{\text{N}}$, a similar linear system can be derived that takes into account the effect of $g_{\text{N}}$.

Since this linear system is obtained from numerical approximation of (4.5), which is a Neumann problem, the matrix $\widetilde{\mathbf{A}}$ is singular. On the other hand, because the system has a small dimension, we may specify the value at one of the vertices to remove the null space of $\widetilde{\mathbf{A}}$, for instance by adding a constant to one of the entries of $\widetilde{\mathbf{A}}$, and then invert the modified matrix. The fact that $\boldsymbol{\widetilde{\alpha}}$ is not unique is irrelevant because the desired final result from the system is the flux as governed by $q^{\tau}_{\zeta\eta}$ which is unique.

The next lemma verifies that the aggregation of elemental calculations satisfies (2.4) for $z\in Z_{\text{in}}$. We note that the same is true for $z\in Z_{\text{n}}$ whose proof we omit for simplicity.

As we see from the above description, what the postprocessing has been able to accomplish is a reconstruction of upscaled (or averaged) conservative flux in $C_{z}$ whose diameter is on the order of $h$. However, MsFEM (and for that matter GMsFEM) always utilizes $h$ that is far greater than the characteristic scale of $k({\boldsymbol{x}})$. In this context, the upscaled flux gathered in Lemma 2 is comparable to entries of $\mathbf{p}$ coming from solving (3.7) (respectively from solving (3.14) for GMsFEM). Since $\mathcal{V}_{\text{ms},h}$ and $\mathcal{V}_{\text{ems},h}$ are built from the multiscale basis functions containing fine scale information of $k({\boldsymbol{x}})$, having this $\mathbf{p}$ allows for calculation of the approximate pressure down to the level of characteristic scale of $k({\boldsymbol{x}})$. However, the current stage of development does not yet allow the upscaled flux to have this capability.

At many practical levels, obtaining upscaled locally conservative fluxes already allows for simulation of multiphase flow and the transport problem in Sect. 2. However, the transport equation (2.2) must then be discretized at the same mesh level as where the approximate pressure is solved, i.e., with $h$ that is far greater than the characteristic scale of $k({\boldsymbol{x}})$. Indeed, this practice is sufficient and suitable for some target applications. Nonetheless, for a large class of problems, such as modeling flows in channelized subsurface with potential localized features, achieving an acceptable accuracy does require the simulation to be performed with a discretization parameter that is comparable with the scale of $k({\boldsymbol{x}})$. For this to occur, we need the flux to be conservative at that finer scale. In the next subsection, we offer a downscaling procedure that allows this capability by taking advantage of the upscaled conservative flux above.

The main driving fact for the downscaling procedure is the realization that after the postprocessing in Sect. 4, we have

which can be thought of a statement of as compatibility condition in $C_{z}$. This means we can proceed with formulating a boundary problem

Here $\widetilde{\boldsymbol{v}}_{h}=-\lambda k\nabla\widetilde{p}_{\tau,h}$ that is evaluated pointwise on segments of $\partial C_{z}$ that are in $\tau$. In fact, this is the same calculation that we did in order to derive (4.11). In this way the compatibility condition is readily satisfied, and thus the solution of (4.14) exists. Similar to the postprocessing in Sect. 4, nonuniqueness of the solution is of no concern since our interest is to gather $-\lambda k\nabla\widetilde{p}_{{}_{C_{z}}}$ from (4.14).

In our implementation, we numerically solve (4.14) for every $C_{z}$ in $\Omega$ with the discretization of $C_{z}$ using the same mesh configuration as that of $\tau$ when numerically solving (3.5). Obviously, the associated mesh parameter should be smaller than $h$ and comparable to the characteristic scale of $k({\boldsymbol{x}})$. We use the standard Galerkin finite element method (3.3) for both of these problems. In addition, the numerical solution of (4.14) is further postprocessed (see [4] for the description) to get a locally conservative flux on the finer scale. Once this is done for all $C_{z}$, we obtain a the downscaled locally conservative flux in $\Omega$, which is in turn used in the simulation of transport equation (2.2).

We should like to emphasize that the main advantage of the proposed series of upscaled and downscaled postprocessings is due to their independence of each other. The procedure is immediately parallelizable, fits well in the framework of CPU-GPU clusters, while the cost for communication is minimal. In this respect, they are indeed in the same spirit as the multiscale basis functions calculations.

In this subsection, we compare the velocities obtained from postprocessing Eq. (2.1) to be coupled to the transport problem in Eq. (2.2). To solve the pressure equation in (2.1) we use Dirichlet conditions $p_{L}=1$ and $p_{R}=0$ on the left and right boundaries of $\Omega$, along with no flow (zero Neumann) conditions on the top and bottom boundaries. We also assume that there is no external forcing, i.e., that $q=0$. Table 1 shows the relative error of the dowscaled velocities obtained by the postprocessing procedure in Sect. 4. We observe that for the deterministic permeability, an appropriate choice of $L_{z}$ (the number of basis functions chosen for each respective node) significantly reduces the error associated with the downscaled velocities. For the random field, we see that the error reduction is less pronounced, yet that a clear error decline is still evident. This difference is due to the nature of the spectrum associated with the eigenvalue problem in Eq. (3.10) and, in particular, the eigenvalue behavior corresponding to the deterministic field lends itself to a more pronounced error decline (see, e.g., [10]). For this initial set of results all MsFEM/GMsFEM solutions were computed on $10\times 10$ coarse mesh.

For a visual representation of the improvement, the computed velocity is plotted on global and localized regions for the deterministic and random permeability fields in Figs. 5 and 6, respectively. In either figure, the left plots show the reference fine scale velocity on the global domain, and on a subregion where the permeability has a large change in value. The center plots show the resulting downscaled velocity on the same regions resulting from the MsFEM (i.e., the case when $L_{z}=1$). We note that significant differences are noticeable between the reference values and this initial case. Finally, the right plots show the computed velocity corresponding GMsFEM with $L_{z}=4$ on the same regions. Figs. 5 and 6 clearly illustrate how the enrichment produces a more accurate representation of the velocity on the fine scale. In addition, this increase in accuracy is particularly evident in regions where the permeability contrast is most extreme. We also note that there are some negative values for the horizontal velocity in the random permeability due to the channelized nature of the permeability in some regions.

In solving Eq. (2.2) we use quadratic relative permeability curves $k_{rw}=S^{2}$ and $k_{ro}=(1-S)^{2}$, along with $\mu_{w}=1$ and $\mu_{o}=5$ for the fluid viscosities. For the initial condtion, the value at the left edge is set as $S=1$ and we assume $S(x,0)=0$ elsewhere. Results for application of the method to two-phase flow are shown in Figs. 7 and 8 using the same permeabilities from Fig. 4. Each of the plots shows the reference saturation in the first row at three time levels, $L_{z}=1$ on the second row at the same three time levels, $L_{z}=2$ in the third row, and $L_{z}=4$ in the last row. The improvement in the deterministic case is significant as can be seen in the Fig. 7 and verified with the relative error shown in Fig. 9. And as expected, for the random field we see a noticeable (but less pronounced) error decline from the profiles in Figs. 8, and the relative error plots in Fig. 9. We mention that in Figs. 7 and 8, the solutions corresponding to the case when $L_{z}=4$ are essentially indistinguishable from the fine scale reference solutions. We reiterate that the size of the resulting linear systems from the pressure equation in Eq. (2.1) will depend on $L_{z}$ and the coarse mesh size. The system for the reference cases are of size $10201^{2}$ for the deterministic permeability and $14641^{2}$ for the random permeability. Since the coarse mesh is $10\times 10$ for both cases, using $L_{z}=1,L_{z}=2,$ and $L_{z}=4$ results in linear systems of size $121^{2}$, $202^{2},$ and $364^{2}$, respectively.

The $L^{2}$ error of the saturation profiles offered in Figs. 7 and 8 is presented in Fig. 9. As was evident from the velocity results, the error is improved by increasing the number of enrichment functions up to a certain threshhold, and then the reduction is minimal as more functions are added. This suggests the number of enrichment functions should be chosen to minimize the overhead in the calculation. As seen above, the size of the linear system to be solved grows quickly as the number of functions is increased and, at some point, there is little to be gain by adding to the enrichment.

For a comparison on the mesh size, the deterministic permeability was used to simulate two-phase flow with $L_{z}=4$ using a $10\times 10$ and $25\times 25$ coarse mesh. Analogous results are presented for the random field, except that a $30\times 30$ coarse mesh was used for the refinement. The results are shown in Fig. 10. As expected, the error is reduced with a finer coarse mesh at the same level of enrichment. This is due to the fact that the coarse mesh is finer and the behavior of the permeability is captured to an acceptable degree by fewer basis functions.