A benchmark study of the multiscale and homogenization methods for fully implicit multiphase flow simulations

Accurate simulation of multiphase ﬂow in subsurface formations is challenging, as the formations span large length scales (km) with high-resolution heterogeneous properties. To deal with this challenge, diﬀerent multiscale methods have been developed. Such methods construct coarse-scale systems, based on a given high-resolution ﬁne-scale system. Furthermore, they are amenable to parallel computing and allow for a-posteriori error control. The multiscale methods diﬀer from each other in the way the transition between the diﬀerent scales is made. Multiscale (ﬁnite element and ﬁnite volume) methods compute local basis functions to map the solutions (e.g. pressure) between coarse and ﬁne scales. Instead, homogenization methods solve local periodic problems to determine eﬀective models and parameters (e.g. permeability) at a coarser scale. It is yet unknown how these two methods compare with each other, especially when applied to complex geological formations, with no clear scale separation in the property ﬁelds. This paper develops the ﬁrst comparison benchmark study of these two methods and extends their applicability to fully implicit simulations using the algebraic dynamic multilevel (ADM) method. At each time step, on the given ﬁne-scale mesh and based on an error analysis, the fully implicit system is solved on a dynamic multilevel grid. The entries of this system are obtained by using multiscale local basis functions (ADM- MS), and, respectively, by homogenization over local domains (ADM-HO). Both sets of local basis functions (ADM-MS) and local eﬀective parameters (ADM-HO) are computed at the beginning of the simulation, with no further updates during the multiphase ﬂow simulation. The two methods are extended and implemented in the same open-source DARSim2 simulator (https://gitlab.com/darsim2simulator), to provide fair quality comparisons. The results reveal insightful understanding of the two approaches, and qualitatively benchmark their performance. It is re-emphasized that the test cases considered here include permeability ﬁelds with no clear scale separation. The development of this paper sheds new lights on advanced multiscale methods for simulation of coupled processes in porous media.


Introduction
Geological formations span large (km) length scales, having heterogeneous properties characterized at high resolutions (cm and below). As for the uncertainty within the integrated field data, typically, several equiprobable realizations of the property fields are generated to study and simulate the fluid flow and transport. Classical simulation approaches are too expensive for such studies. Therefore, advanced simulation methods are required to allow for an accurate representation of the heterogeneous properties. At the same time, they should pro-These approaches are different in the sense that the multiscale method deals with crossing the solution, e.g., the pressure, across the scales Aarnes and Hou (2002) ; Jenny et al. (2003) ; Hajibeygi et al. (2008) ; Chung et al. (2015) , whereas in the latter effective, lower-resolution parameters and functions like the permeability or the transmissibility, are derived Weinan and Yue (2004) ; Abdulle et al. (2012) ; Weinan et al. (2007) ; Li et al. (2020) ; Singh and Wheeler (2018) ; Vasilyeva et al. (2020) . Moreover, while the multiscale basis functions have been expressed in a purely algebraic formulation Wang et al. (2014) , the same does not hold for the homogenization approach. Specially the integration of homogenized parameters within the fully implicit framework in an algebraic manner has not yet been developed so far. The present work is a first step in this direction.
At the same time, the two methods have many similarities. Both find their mapping strategy via local solutions of the original governing equations with local boundary conditions. Multiscale basis functions often employ reduced-dimensional boundary conditions Tene et al. (2015) ; Møyner and Lie (2016) , while homogenization schemes impose periodic boundary conditions on local problems, and consider local representative micro-structures even in the case of non-periodic properties Allaire (1992) ; Abdulle and Weinan (2003) ; Arbogast and Xiao (2013) ; Bastidas et al. (2019) ; Brown et al. (2013) . Both methods are effective for global equations within the fully coupled system of local-global unknowns, e.g., the global pressure and the local saturation. Both have been extended to nonlinear and geologically complex models Amanbek et al. (2019a) ; HosseiniMehr et al. (2018) ; Singh et al. (2019a) . Recent developments of these two classes of approaches have introduced a fully-implicit dynamic multilevel simulation framework (ADM) in which heterogeneous detailed geo-models are mapped into adaptive dynamic coarser mesh Cusini et al. (2018) ; Faigle et al. (2014) ;  ; Carciopolo et al. (2020) .
The ADM method develops a fully-implicit discrete system for coupled flow and transport system of equations, in which each equation can be represented at a different resolution than the defined fine-scale one. More importantly, the procedure can be done fully algebraic, with the dynamic mesh resolution defined based on a front-tracking strategy. In contrast to the rich existing literature of Adaptive Mesh Refinement (AMR) methods Pau et al. (2009Pau et al. ( , 2012 ; Berger and Oliger (1984) ; Schmidt and Jacobs (1988) ; Edwards (1996) ;Sammon (2003) ;  , ADM can be defined as an adaptive mesh coarsening strategy which is conveniently applicable for heterogeneous and nonlinear coupled systems Cusini et al. (2016) .
Irrespective of the choice of the dynamic mesh strategy, it is always a challenge to construct adaptive multiscale entries of the implicit systems. The ADM method so far has included multiscale basis functions Cusini et al. (2016) . In addition, homogenization methods have also been developed for multiphase simulations on dynamic grids Amanbek et al. (2019a) ; Cusini et al. (2019) . In this context, two aspects can be of interest: the study of the homogenization-based coarser system entries and the development of a benchmark study of the quality of the two approaches of ADM-multiscale (ADM-MS) and ADM-homogenized (ADM-HO) for coupled implicit multiphase flow scenarios.
This paper develops such a unified framework in which the ADM method is extended to account for both multiscale and homogenization schemes for multiphase flow simulations. This development makes it possible to allow for different coarse-scale entries for dynamic simulations, and importantly to benchmark the two classes of multiscale and homogenization strategies. Noteworthy is that, once the effective parameters are computed, all other homogenization procedures are implemented algebraically. This is done by introducing constant unity local basis, with the support of primal (non-overlapping) coarse-scale partitions. The multiscale ADM is implemented fully algebraic since local basis functions are also solved algebraically over the overlapping (dual) coarse grid domains Zhou and Tchelepi (2012) . The outcome of this development is made available to the public via an open-source DARSim2 simulator, https://gitlab.com/darsim2simulator. Numerical test cases are considered for the challenging, highly heterogeneous SPE10 Christie and Blunt (2001) . The number of active grid cells, pressure and saturation errors, and the solution maps are all reported in detail. The development of this paper sheds new lights in the application of multiscale and homogenization approaches in advanced next-generation environments for field-relevant simulation scenarios.
The paper is structured as follows. Next, in Section 2 , the mathematical model is stated briefly. Section 3 presents the computational framework for both multiscale and homogenization ADM methods. Section 4 presents the test cases, and conclusions are drawn in Section 5 . The Appendix gives more details on the multiscale and the homogenization approaches.

Governing equations
We consider flow of two immiscible and incompressible phases of and through a heterogeneous porous medium. At the Darcy scale, mass balance for the phase i ∈ { , } reads (1) Here, is the porosity of the medium, i [kg/m 3 ] and S i are the density and saturation of the phase i , respectively. The phase mobility tensor i is equal to ∕ , where K [m 2 ] is the rock absolute permeability tensor, and = ( ) is the saturation-dependent relative permeability of phase i . Moreover, [Pa.s] is the phase viscosity. For the ease of presentation, the two phase pressures are assumed equal, = = [Pa] (see e.g. Aziz and Settari (2002) ). However, the extension to models involving a saturation dependent capillary pressure is also possible. In addition, g [m/s 2 ] is the gravitational acceleration which acts on ∇ z direction, and q [1/s] is the phase source term.
Here it is assumed that the two fluids are occupying completely the pore space, and no other fluid phase is present. This gives the constraint + = 1 , which reduces the number of unknowns in the above equations to two: S (in short from here on, S ) and p . Finally, the model is completed by initial conditions for the saturation, and with boundary conditions. We do not specify them explicitly since none of them play a role in the multiscale strategy.
The fully-implicit coupled simulation approach Aziz and Settari (2002) estimates all the parameters at next time step ( + 1) . As such, the semi-discrete nonlinear residual for the phase i ∈ { , } reads (2) For finding the solution pair ( +1 , +1 ) one needs to employ a linearization scheme. Here we restrict the discussion to the Newton scheme, which is 2nd-order convergent but requires a starting point that is close enough to the solution. In other words, the time step may be subject to restrictions also depending on the mesh size. Alternatively, one may consider approaches like the modified Picard Celia et al. (1990) or the L-Scheme Radu et al. (2017) , which are less demanding from the computational point of view, or more robust w.r.t. the starting point and mesh resolution, but converge slower than the Newton scheme Bastidas et al. (2019) . Applied to (2) , the Newton linearization reads which can be expressed algebraically as +1 = − , i.e., In each time step, the linear Eq. (4) is solved iteratively (inner loop) several times until nonlinear convergence (outer loop) is reached. The overall computational complexity of the simulation depends highly on the complexity of the solution of this linear system. Advanced multiscale and homogenization methods aim at solving this linear system on a dynamic multilevel mesh. Note that, as shown before Cusini et al. (2018) , the overall efficiency of any advanced method should include not only the speedup of solving the linear Eq. (4) but also the count of the Newton (outer) loops. Next, the ADM method based on multiscale and homogenization formulations is presented.

ADM Framework formulation
The fully-implicit linear system (4) is too expensive to be solved for real field scenarios. A multilevel dynamic mesh, as shown in Fig. 1 , is generated within the ADM framework, based on an error estimate strategy. The error estimate is developed based on a front tracking criterion, which applies fine-scale grids only at sub-regions with sharp gradients. The fine-scale system is then algebraically reduced into this multilevel grid, through sequences of restriction and prolongation operators. To obtain the ADM grid, first, sets of = × hierarchically nested coarse grids are imposed on the fine mesh. Here, l indicates the coars-ening level. Moreover, l is the coarsening ratio which is defined as for two-dimensional (2D) domains. The ADM grid is constructed by assembling a combination of cells at different resolutions within the computational domain. By using the sequence of restriction ( R ) and prolongation ( P ) operators, one can express the ADM system as Here, ̂ −1 is the restriction operator which maps the parts of the solution vector that are at level ( − 1) to level l . Similarly, the prolongation operator ̂ −1 maps the parts of the solution vector that are at level l to level − 1 . Once the ADM system (6) is solved, the approximated fine-scale solution ′ 0 can be acquired by prolonging the ADM solution ̂ ADM , i.e.
The ADM Restriction ̂ −1 and prolongation ̂ −1 operators are assembled using the static multilevel multiscale restriction −1 and prolongation −1 operators, respectively. They are constructed only at the beginning of the simulation and are kept unchanged throughout the entire simulation.
The static prolongation operator −1 is constructed as an assembly of the locally computed basis functions at each coarsening level l and reads Here, ( ) −1 and ( ) −1 are the two main diagonal blocks corresponding to main unknowns (i.e., pressure p and saturation S ). In the case of using the homogenization scheme, i.e. ADM-HO, as will be described in Section 3.3 ), constant basis functions for pressure are used. However, for the multiscale-based ADM, i.e. ADM-MS, as will be described in Section 3.2 , locally-computed basis functions are used. Note that the saturation prolongation operator for both approaches is constant unity function at all coarsening levels, which represents the conservative finite-volume integration.
The static restriction operator −1 reads In this work, a finite-volume restriction operator is used to guarantee local mass conservation, i.e.

ADM Using multiscale (ADM-MS)
In the ADM-MS method, the prolongation operator for pressure is found based on multiscale basis functions. These local basis functions are computed algebraically Wang et al. (2014) , based on the steady-state pressure equation. In this study, the incompressible flow equation (elliptic pressure equation) is used to construct the multiscale basis functions Tene et al. (2015) . An example of a basis function is shown in Fig. 2 .
The coarse grid construction and computation of multiscale basis functions are explained in more details in Appendix A .

ADM Using homogenization (ADM-HO)
Homogenization is another method that can be applied to problems involving multiple scales. In this method, one uses the mathematical models at micro (fine) scale (1) to derive effective upscaled models  and parameters in which the rapidly oscillating characterstics are averaged out. In doing so, the upscaled model may have a different structure than the ones at the fine scale. We refer to Amaziane et al. (2017) ; Bourgeat et al. (1996) ;Hornung (1997) ; van Duijn et al. (2007) for theoretical details.
The goal of this work is to build a unified ADM platform where the multiscale and homogenization methods can be compared. Therefore, here, the homogenization method is used only to construct effective properties at the dynamic multilevel mesh. In this setup, the homogenized properties of ADM-HO at multilevel mesh are found as in ADM-MS by solving local flow (pressure) equations based on an incompressible (elliptic) equation.
More precisely, one assumes that a scale separation holds and doubles the spatial variable into a fast and a slow one. The method relies on the homogenization ansatz , meaning that all quantities in (1) Henning et al. (2015Henning et al. ( , 2013 to develop effective numerical simulation schemes even in case of non-periodic media. More details about the homogenization procedure can be found in Appendix B . In the present context, for a given fine-scale effective permeability K and for each coarsening level l , an effective permeability tensor K l is computed locally in a pre-processing step. First the domain Ω ⊂ ℝ 2 is divided into coarse cells Ω l that correspond to a partition of the domain Ω as shown in Fig. 3 .   Fig. 4. Example of the local solutions 1 (top right, for x -direction) and 2 (bottom right, for y -direction) for a coarse cell inside a 2D domain. The heterogeneous permeability field is also shown for the entire domain (left).
For each coarse cell Ω l at level l , the components of the effective permeability tensor are calculate as for , = 1 , 2 . Here j are the periodic solutions of the pressure equation on local domains (known as micro-cell equation in HO literature), i.e.
We remark that { } 2 =1 is the canonical basis of dimension 2 and K is the above mentioned permeability tensor. To guarantee the uniqueness of the solution j one assumes that its average value over the local coarse cell Ω l is 0.
To determine the value of the effective permeability tensor at each coarse cell Ω l , two local (micro-cell) problems (12) are solved for each spatial direction in 2D. Fig. 4 provides an illustration of these local solutions for a coarse element.
Note that the local problems (12) capture the rapidly oscillating characteristics within a coarse element, completely decoupled from other coarse elements. The homogenized parameters, like multiscale bases, are computed at the beginning of the simulation. Fig. 5 illustrates the calculation of the effective permeability at different levels.
The homogenized parameters are used to construct the coarse system entries. More precisely, the homogenized value in a coarse cell is distributed equally to the fine cells constructing it. Then the fine-scale Jacobian and residual are computed with the fine-scale saturation field. This system is then mapped to the ADM resolution by setting prolongation operators in (6) to unity. This is a convenient procedure, developed in this work, to integrate the numerical homogenization method with an existing advanced simulator.
Notice that based on the features of the permeability tensor K , the resulting effective parameter K may depend on the macro-scale location and the size of the coarse-scale partition. Nevertheless, one can show that in practice, the adaptive refinement of the mesh is an important aspect that improve the calculation of the effective parameters (see Bastidas et al. (2019) and Fig. 5 ).
The ADM procedure is sketched in Fig. 6 . More details about the role of the homogenization in the offline stage and the complete algorithm of ADM-MS and ADM-HO can be found in Algorithm 1 .

Simulation results
To benchmark the homogenization and multiscale based solutions for the dynamic mesh on heterogeneous media, two heterogeneous nonperiodic permeability fields from the top and bottom layers of the SPE 10th Comparative Solution Project Christie and Blunt (2001)    injection well, while the reservoir fluid is produced from the production well. The locations of the injection and production wells are specified in each test case. Table 1 shows the input parameters of the fluid and rock properties used in all test cases. Note that the density and viscosity ratios are assumed to be 1. Since compressibility and gravitational forces are both neglected, the density values have no influence on the results.
The numerical results provided by the ADM-MS and ADM-HO methods are compared to those obtained from simulation at fine scale (reference). Both ADM methods employ the coarsening ratio of 3 × 3 with two coarsening levels. This is set according to the size of the domain.

Test case 1: SPE10 top layer
In this test case, one injection well and one production well are placed in the bottom left corner and top right corner of the domain, respectively. The simulation time is = 1000 [days] and the results are reported on 100 equidistant time intervals. The permeability distribution of the SPE10 top layer is shown in Fig. 7 . Fig. 8 shows the homogenized version of the permeability at two different levels. We highlight that the homogenized permeability at both coarse levels preserves the structure of the original fine-scale permeability. The high and low permeable zones remain clearly detectable.    The saturation and pressure fields at the final time step are shown in Fig. 9 and Fig. 10 , respectively.
From these results, it is understood that ADM-HO on a coarse cell containing high and low permeable fine cells can lead to a higher flow leakage, as compared to fine-scale and ADM-MS approaches. This effect can be seen in Fig. 9 . Fig. 11 , illustrating the adaptive mesh at 2000 days after injection. Notice that the refinement of the permeability is most dominant at the saturation front, due to the chosen mesh refinement criterion. For this figure, the coarsening threshold value is Δ = 0 . 3 , i.e., a cell is successively coarsened if ΔS is lower than 0.3.
The error history maps for both ADM-MS and ADM-HO are shown in Fig. 12 . The relative errors, presented in Fig. 12 and Fig. 14 , are expressed in terms of the L 2 norm over the entire medium, calculated with respect to the fine-scale solution as  The results indicate that the homogenization-based simulations have higher errors compared with the multiscale-based simulations. They both have similar average usage of active grid cells, with ADM-MS having slightly fewer grid cells. This is shown in Fig. 13 . Note that the grid cells around wells are kept at the fine-scale resolution permanently. Furthermore, for tighter error tolerance values, the quality of both approaches become comparable.
Fig. 14 provides the average pressure and saturation errors together with the average percentage of active grid cells during the whole simulation time as functions of the coarsening criterion threshold.

Test case 2: SPE10 bottom layer
In the second test case the permeability distribution of the SPE10 bottom layer, presented in Fig. 15 , is considered. The location of the injection and production wells are the top left and the bottom right corners, respectively. The simulation time is 20 [days]. All other simulation parameters remain unchanged. Fig. 16 shows the homogenized permeability values at two different levels. Due to the many high contrast channels, more active cells are employed compared with the SPE top layer, as shown in Fig. 17 .
The saturation and pressure maps at the final time step are shown in Fig. 18 and Fig. 19 , respectively.
Similar to the previous test cases, Fig. 20 compares the error between the two ADM approaches. Moreover, in Fig. 21 , the percentage of active grid cells per each time-step is shown.   The results indicate a noticeable difference in the errors of ADM-MS and ADM-HO. The pressure error in ADM-HO is significantly higher since ADM-HO uses homogenized effective parameters. This aspect can be improved by employing first order corrections. However, such an      Fig. 19. Pressure profiles at 20 days. The threshold value for the front tracking criterion is Δ = 0 . 3 . approach would deviate from the ADM framework, and requires more computational effort, therefore it is not adopted here. ADM-MS instead employs multiscale basis functions. Due to the more accurate pressure calculations, the ADM-MS saturation error is also lower than that of ADM-HO. The difference in the percentage of active grid cells used in the two approaches is less noticeable than the difference in the errors. However, the ADM-HO uses more active grid cells, especially in this SPE10 bottom layer test case.

Conclusion
Homogenization and multiscale methods have been developed and evolved during the past decade as promising advanced simulation approaches for large-scale heterogeneous systems. In this work, the two methods were investigated, extended into a unified fully-implicit framework, and benchmarked for simulation of multiphase flow in porous media. It was shown that the two methods allow the construction of coarser level systems, and both rely on local solutions to find their corresponding maps. While homogenization methods deliver effective models and parameters, multiscale methods find an interpolation of the solution (pressure) across scales. This is the main difference between the two approaches.
For highly heterogeneous test cases, it was shown that the two approaches provide accurate solutions. With the developed multiscale numerical strategies, the ADM-MS solutions are more accurate when compared to ADM-HO. The use of a constant effective parameter instead of local multiscale basis functions results in relatively higher errors. In addition, using constant unity prolongation operator along with the effective coarse-scale parameters allows for straightforward implementation of the ADM-HO method for domains with non-periodic permeability fields. The study of this paper sheds new light on the application of multiscale and homogenization methods for real-field simulation of multiphase flow in porous media. Note that the computational costs of the two approaches were comparable, as they applied almost the same active cells during the simulation. Ongoing study includes benchmark studies of ADM-HO and ADM-MS for 3D fractured porous media, on compilable simulation platform, which allows scientific CPU comparison study.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix A. ADM based on multiscale method
In the ADM-MS approach, multiscale finite volume method (MSFV) Jenny et al. (2003) ; Cortinovis and Jenny (2014) is used to compute local basis functions at multiple coarsening levels. The computation of basis functions Φ is done by solving the incompressible fluid flow equation (elliptic part of the mass balance) Cortinovis and Jenny (2017) which reads The incompressible basis functions are found to be the most efficient ones, compared with the compressible and more complex formulations Tene et al. (2015) . The first step is to impose coarse grids on top of the fine mesh, for the coarse level 1. Here, to simplify the visualization, a 2D 15 × 15 discrete domain is considered (see Fig. 23 ). By connecting the centers of the coarse cells, the dual-coarse grid is obtained. The dual grid makes an overlapping partitioning of the fine-scale domain, with 3 categories of interior (white), edge (green), and vertex (blue) cells. The coarsening ratio in the illustrated example of Fig. 23 is 5 × 5.
The Eq. (15) is solved at each dual coarse grid h and for each coarse node (vertex) k , i.e., −∇ ⋅ ( ⋅ ∇Φ ℎ ) = 0 . In order to solve this local system, Dirichlet boundary conditions of 1.0 (for the corresponding coarse node) and 0.0 (for the other three coarse nodes) are imposed. These Dirichlet values allow to solve the basis functions on the edges, if a reduced dimensional (1D) elliptic problem is considered. The solution at the edge and vertex cells are then imposed as Dirichlet boundary condition for the full 2D problem. The solution of this well-posed system is the basis function of the corresponding coarse node at the corresponding dual coarse grid. Fig. 24 shows a schematic of the mentioned dual coarse grid h and an example of a basis function belonging to the bottom left coarse node ( Φ ℎ 1 ). Fig. 25 shows all the four basis functions for the mentioned dual coarse grid h .
The combination of the basis functions at all the dual coarse grid cells surrounding the corresponding coarse node forms the basis function belonging to that coarse node. Fig. 26 illustrates an example of a    To obtain the basis functions at higher coarsening levels, the hierarchically nested coarse grid is constructed on the same domain. The same procedure is followed to compute the basis functions at higher coarsening levels. Fig. 27 shows the coarse grid construction at 2 consequent coarsening levels, for a 2D domain with 75 × 75 fine cells.
Note that, according to the vast multiscale literature, construction of basis functions can be done purely algebraic, once the wire-basket decomposition of the fine cells into Vertex, Edge, Face, and Interior is known Tene et al. (2016) . A partitioning method should be applied for complex mesh Møyner and Lie (2016)

Appendix B. ADM based on homogenization theory
The main idea of the ADM-HO is to use a homogenized version of the permeability K instead of volume-averaged permeabilities at differ- ent coarse levels. In doing so, two assumptions are commonly made: the permeability is periodic at each of the coarsening levels, and the scales are well separated. We refer to Allaire (1992) for the rigorous mathematical support of this approach. Although the test cases considered here do not satisfy the two assumptions stated before, the homogenization idea can still be considered for developing multiscale simulation tools, and in this sense we refer to Bastidas et al. (2019)  At each coarsening level l , we call micro-scale cell (i.e., local coarse cell) the region Ω l wherein the parameters change rapidly. For each Ω l , the characteristic length is , where L is the characteristic length for the macro-scale domain Ω. The factor ∶= reflects the scale separation. To identify the fast changes in the parameters we double the variables and define the fast variable y ∶= x . In the non-dimensional setting, Ω can be written as the finite union of the local cells Ω l . We let Ω = ∪  Ω for some set of indices .
For calculating the homogenized permeability we consider an auxiliary elliptic problem: Auxiliary problem. Given a fine-scale permeability K , find a function u that satisfies ∇ ⋅ ( ⋅ ∇ ) = 0 in Ω, = 0 on Ω.
Here the boundary conditions are specified for completeness. We employ the homogenization ansatz , meaning that the unknown u can be written as where each ̂ ( x , y ) is periodic. Due to the doubling of variables, the gradient and divergence operators become ∇ = ∇ + 1 ∇ and div = div + 1 div .
Collecting the terms with factor −2 we obtain for each domain Ω l ∇ ⋅ ( ∇ ̂ 0 ) = 0 for all y ∈ Ω , Clearly, any 0 = 0 ( x ) which does not depend on the fast variable y is a solution of the problem (18) . Thus, one can prove that all solutions depend only on x . The function u 0 is in fact the macro-scale approximation of the original unknown u . On the other hand, for the −1 terms one has ∇ ⋅ ( ( ∇ ̂ 0 + ∇ ̂ 1 )) = 0 for all y ∈ .
One can determine ̂ 1 as a function of ̂ 0 and eliminate it from the system. To this end, ̂ 1 is written as a linear combination of functions j , and with ̂ 0 as coefficients The functions j are solutions of the micro-cell (local domain) problems defined on each Ω l : = 0 for all y ∈ , isperiodicin .
Here, { ⃗ } d =1 is the canonical basis of dimension 2. To guarantee the uniqueness of the solution one requires that ∫ Ω = 0 y , for all Ω .
The homogenized permeability in the auxiliary problem is obtained by considering the terms of order 0 , in which ̂ 2 appears. Averaging these terms and using the periodicity of the functions on the right hand side of (17) , one obtains that the auxiliary unknown ̂ 0 ( x ) solves the homogenized problem ∇ ⋅ ( ⋅ ∇ ̂ 0 ) = 0 in Ω, Here the matrix valued function ∶ Ω → ℝ 2×2 has the elements Note that these steps are carried out only to determine the effective permeability tensor K l . The solution ̂ 0 of the effective problem (19) is not of interest here. In other words, we use the auxiliary problem for the sole purpose of defining an effective parameter that could be used at each level of the dynamic multilevel algorithm.