Research Article  Open Access
Adaptive DoubleDiffusion Model and Comparison to a Highly Heterogeneous MicroModel
Abstract
Doublediffusion model is used to simulate slightly compressible fluid flow in periodic porous media as a macromodel in place of the original highly heterogeneous micromodel. In this paper, we formulate an adaptive twogrid numerical finite element discretization of the doublediffusion system and perform a comparison between the micro and macromodel. Our numerical results show that the micromodel solutions appear to converge to the macromodel linearly with the parameter ε of periodic geometry. For the twogrid discretization, the a priori and a posteriori error estimates are proved, and we show how to adapt the grid for each component independently.
1. Introduction
When modeling phenomena in highly heterogeneous media, one frequently finds that the coefficients of differential equations describing these phenomena vary by several orders of magnitude between closeby locations. Consequently, the solutions to the exact, that is, micromodels, vary at multiple separate spatial and time scales, and it is convenient to work with their spatial and temporal averages, that is, with the solutions to macromodels. Various multiscale modeling techniques have been introduced with the aim to derive, analyze, and approximate micromodels by the macromodels. In particular, the mathematical framework for construction and analysis of the average, that is, homogenized, models for multiscale media is well understood, see, for example, [1–3]. However, there are few results devoted to the comparison between the micro and macromodels for evolution equations and to their adaptive numerical discretizations.
Consider the following transient diffusion model which describes the density of singlephase slightly compressible flow in a porous domain : where the coefficients , denote, respectively, the porosity and conductivity.
In this paper, we propose a twogrid adaptive finite element approximation of the following multiscale averaged model for (1.1) in highly heterogeneous media such as porous media with fracturesproposed in [4, 5]. The system (1.2a) and (1.2b) describes the average behavior of fluids in media with two or more subregions of distinct features, where the subregions may or may not be connected globally at microscale. At macroscale the interaction of fluids associated with these subregions is modeled with the exchange term , and the coefficients , are computed from microscale geometry and coefficients. The system (1.2a) and (1.2b) is useful in applications in porous media as well as in description of other multiscale evolution phenomena in heterogeneous media, for example, heat conduction. In addition, its special, degenerate, case of (1.2a) and (1.2b) with known as the WarrenRoot model [6]is a popular way of modeling of subscale diffusion accompanying advection and adsorption, see [7–10].
The macromodel (1.2a) and (1.2b) of interest to this paper was introduced in [4] based on heuristic arguments. The model, while intuitively clear and useful for practical purposes, for a long time was lacking a rigorous multiscale analysis and interpretation via an associated appropriate micromodel. Only recently in [11] it was shown rigorously by twoscale convergence that (1.2a) and (1.2b) is indeed a homogenized limit of (1.1) in certain highly heterogeneous media. Most interesting is that the original microstructure which leads to (1.2a) and (1.2b) is composed of three rather than two media as in [12]. Also, the derivation in [11] demonstrates clearly that in twodimensional geometries and more generally for quasistatic models in disconnected microscale subregions we have .
However, the analysis in [11] does not include a direct comparison of (1.1) and (1.2a) and (1.2b), and this motivated our research. We compare the computational solutions to (1.2a) and (1.2b) and (1.1) quantitatively in function of the microstructure parameters. In addition, since typically , it is natural to suggest that in numerical simulation and should be approximated on separate grids, since they evolve differently and since this can lead to more efficient numerical algorithms. In fact, the spatial grids should be chosen adaptively guided by some a posteriori estimators and should be allowed to vary in time. The idea of using two grids extends that in [13] where we considered a stationary analogue of (1.2a) and (1.2b) and proved upper and lower bounds for residual type a posteriori estimators. The estimators proposed in [13] were shown to be robust, that is, insensitive to and to other parameters. In this paper we extend those to the evolution system (1.2a) and (1.2b) and present an a posteriori error estimator of residual type based partly on the work in [14–16] for scalar equations.
The aim of this paper is thus twofold, and the paper is organized as follows. In Section 2 we introduce the notation and the generic flow model as well as its numerical approximation. In Section 3 we develop the details of micro and macro models and compare their discrete solutions in function of a characteristic parameter of the microstructure . In Section 4, we define the numerical approximations to these models and formulate the apriori and aposteriori estimates for the approximation error. We then illustrate the use of aposteriori estimators with an adaptive example. The appendix provides the technical proofs of some of the results derived in this paper.
2. Notation and Preliminaries on the Flow Model
Here we introduce the notation and preliminaries. We follow, for example, [17, 18]. Let the domain of flow , , be a bounded polygonal region with piecewise smooth boundary , where , are disjoint regions of the boundary where Dirichlet and Neumann boundary conditions are prescribed, respectively.
For any subset , we use the standard notation for Lebesgue and Sobolev spaces , . These are equipped with the usual seminorms , norms and the usual scalar product (s) . We denote . If , the subscript will be omitted. We also let under the standard norm in , .
For any and any vector space , the space consists of all squareintegrable functions with values in such that is finite. The space is defined similarly.
2.1. Variational Form of MicroModel
Consider (1.1) as a model for singlephase slightly compressible flow in porous media in which the gravity terms have been ignored; see [19] for derivation. Here and denote the nonnegative coefficients of porosity and conductivities, respectively, and accounts for external sources. We assume for simplicity that is a scalar quantity.
It is well known that (1.1) is well posed if appropriate boundary and initial data are imposed. We will use For simplicity we assume that and that the flow is driven by . However, in simulation examples we use , .
The variational formulation of problem (1.1) and (2.1) reads
Find so that For any , , there is a unique solution to (2.7); see [20, Theorem 11.1.1, page 366]).
2.2. Variational Form of MacroModel
Now consider (1.2a) and (1.2b) and assume that , , , , and are bounded uniformly from above and below by positive constants and are independent of time. We consider separately the particular case of (1.3a) and (1.3b) when .
Our variational formulation for (1.2a) and (1.2b) uses the theory developed in [17, 18] for evolution problems with maccretive operators by the HilleYosida Theorem. For analysis of classical solutions see [21, 22].
For the macromodel (1.2a) and (1.2b) and if , we assume the following initial and boundary conditions defined for : As before, we assume that for now.
Let , , . The variational form of (1.2a), (1.2b), and (2.3) is as follows:
Find , , such that The well posedness of this system (2.4) follows immediately by showing that the operator is maccretive. The accretiveness is straightforward by showing that in the product space. The maximal property follows from that of the operators , .
Now, if , we do not impose boundary conditions on . The variational formulation of (3) and (2.3), reads
Find , such that The well posedness of (2.6) follows similarly from the HilleYosida theorem. Also, has the minimum of regularity of its initial data and of .
2.3. Notation on Numerical Discretization
For computations of solutions we use conforming (Galerkin) finite elements for spatial discretization and implicit Euler time stepping. We partition the time interval into subintervals , , such that . We denote and for any function .
For spatial discretization, we adopt standard finite element nomenclature that can be found in textbooks such as [23–25]. At each time step , we denote by , , a family of admissible and shaperegular partitions of into a finite number of elements. For any element , we let be the diameter of and denote . Also, is the set of all edges of . For any edge , we let denote the diameter of the edge . Finally, we denote by the space of polynomials of degree in .
Numerical Approximation of MicroModel
At each time step , we define
and seek an approximation such that for all Here denotes an interpolation or projection operator into .
It is standard [20] that there is a unique solution for (2.8a) at each , . A priori estimates for the error between the solution to (2.2) and that of (2.8a) and (2.8b) can be found in the literature, see, for example, [20, 25, 26]. For a posteriori error estimates for parabolic problems see [15, 27–29]; the methods in the latter two papers are most relevant to our paper.
Numerical Approximation of MacroModel
Consider now a fully discrete approximation to the solutions to (2.4) or (2.6). The formulation with a fixed grid immediately extends (2.8a) and (2.8b), and appropriate error estimates can be proven. More generally, each of can be associated with its own grid ; see details in Section 4.
3. Comparison of Micro and MacroModels
The difficulties associated with numerical approximation (2.8a) and (2.8b) to (1.1) arise when , are highly varying coefficients. In particular, consider a porous domain composed of three disjoint regions , , see Figure 1, with positive constants so that Here denotes the characteristic function of .
An accurate numerical approximation to (1.1) with (3.1) requires that the grid is very fine so that it lines up with the interfaces . The computational complexity of the associated numerical implementation is very large whenever the geometry of interfaces is complex, but this can be overcome by a modeling approximation. In [11] it is shown that (1.2a) and (1.2b) provides solutions close to those of (1.1), with the coefficients which do not vary locally. Thus a grid for (1.2a) and (1.2b) can be chosen depending on the global macrodynamics rather than local microdynamics and, consequently, it can be much coarser than .
3.1. Derivation of MacroModel
One successful avenue to reduce the complexity of a highly heterogeneous model is to consider a homogenized or upscaled problem, whose coefficients are derived from the original coefficients [1, 2, 30], and which can be simulated on a coarse grid quite accurately. However, for timedependent problems, this leads to inaccurate representation of the global dynamics, especially if and/or .
The doublediffusion models proposed in [4], the relaxation model from [6], and doubleporosity models derived in [12] are families of successful modeling strategies for parabolic PDEs with periodic highly heterogeneous media. The recent general derivation in [11] is most general and includes those from [4, 6, 12] as special cases.
The doublediffusion model developed in [4] was shown in [11] to be a limit of a certain micromodel with three regions arranged periodically. The rigorous derivation in [11] allows for each of the regions to be connected. Furthermore, it assumes that separates and , that is, that . The model in [11] develops a general nonlocal exchange term acting across . Its special case equivalent to the model in [4] assumes thereby making a quasistatic assumption and allows to compute the coefficient explicitly.
Doubleporosity models proposed in [12] assume that only one of the regions is connected and so effectively one can set . However, . In the macromodel the second equation (1.2b) is replaced by a local micromodel equation, the exchange follows through macromicro boundary conditions and is nonlocal in nature. We refer to [31–35] for various analyses, models, approximations, and extensions of doubleporosity models.
Of course, all three regions , , and can be connected only in . In only one of these regions, say , can be connected; see Figure 1(a). As we point out below, this implies , and the macromodel (1.2a) and (1.2b) in becomes the WarrenRoot relaxation model (3) as in [6].
In this paper we are interested in the doublediffusion models with a constant from [4, 6] and its corresponding quasistatic micromodel from [11]. The exchange region is assumed to be small and acts like a thick skin separating from . In consequence, the fluid flows effectively only in two regions and ; this is represented by the two subequations in (1.2a) and (1.2b). Details, comparison, and simulations are shown in what follows.
3.1.1. MacroModel Parameters
To compute parameters , , , , and of the macromodel (1.2a) and (1.2b), one has to calculate local averages and compute auxiliary solutions of differential equations. These depend on , and on the geometry of . We follow closely [11].
First we formalize the notion of periodicity and heterogeneity in . Without loss of generality we assume that is globally connected and ; all other coefficients in (3.1) are positive constants. We follow the usual structure [1, 11, 12]. Let the unit cube be divided into three distinct regions , , and , such that and as illustrated by Figure 1(b). Denote by the characteristic function of , , extended periodically to all . We assume that the domains , , have a smooth boundary. Also, we define the periodic characteristic functions , , . Thus, for a fixed , is subdivided into three distinct regions , . Another subdivision of is into periodic cells.
Macroporosities
From [11, equation 16, page 206] we have .
Macroconductivities
From [11, equations 16, 14, page 205206] we have
where , for ;, is the solution to an auxiliary PDEOne can prove [1, 3, 11] that is a symmetric matrix.
Now consider . Here the region cannot be connected. Thus, the boundary condition (3.3c) does not apply and (3.3a), (3.3b), and (3.3c) admit the trivial solution for . Consequently, is the null matrix, and the model (1.2a) and (1.2b) becomes the WarrenRoot model (3) proposed in [6].
Exchange Term Parameter
In general, the flow between and across has transient character, and an appropriate term describing it must be nonlocal in nature.
However, for very small , the flow is of quasistatic nature as discussed in [4]. The exchange term then has the form . To compute , we consider the solution ofThe parameter see [11, equation 19, page 209] is given by
We note that the general nonquasistatic case is when is included in (3.4d) thereby changing the character of exchange term to nonlocal. Conversely, (3.4a), (3.4b), (3.4c), and (3.4d) can be considered as its special case when .
3.1.2. Numerical Computation of MacroModel Parameters
The solutions of the auxiliary PDEs (3.3a), (3.3b), and (3.3c) and (3.4a), (3.4b), (3.4c), and (3.4d) are approximated using finite elements with and . See [36] for treatment of periodic boundary conditions.
In what follows we ignore the distinction between the exact and numerical values of macromodel parameters , and provide examples of their calculations for various values of and choices of geometry of .
In all examples hence . Also, in all examples and the geometry of is axisymmetric, and thus is a scalar constant.
Example 3.1. Here , , , and . We obtain
Example 3.2. Let , . , and . We obtain
Example 3.3. Here , , , and . We obtain
Example 3.4. Let , , , and . We get
Example 3.5. , , , and . We get
3.2. Comparison of the Models
Now we can compare the solutions of the micromodel (1.1) to those of the macromodel (1.2a) and (1.2b). In [11] it was shown that the latter twoscale converge to the former as . In addition, the analysis in [11] suggests that it is appropriate to compare to . However, this notion of convergence does not give information about the rate at which may converge, and it involves special periodic test functions.
In this paper we estimate this rate by comparing their numerical approximations directly without the use of any test functions but rather in a certain metric of interest. To the best of our knowledge such comparison or convergence rate was not discussed elsewhere.
Strictly speaking, the twoscale convergence proof considered in [11], and a similar proof for the doubleporosity model in [12], includes scaling of with . This scaling is a formal device necessary to preserve certain parts of the boundary value problem under investigation in the limit as . However, the homogenization limit is intended to serve only as an approximation of the true model, which has a given fixed set of parameters. We do not include this scaling in our computations. Rather, we treat each example, for a given , as a data set in its own merit, rather than as an element of a sequence intended to twoscale converge to the limit.
3.2.1. Setup of Computational Experiments
We set up simulations for comparison using compatible data for micro and macromodels; we set and choose initial and boundary data driving the flow with interesting dynamics. For porosities and conductivities, we use the values computed in Section 3.1.2.
Let and . On we define . On the lateral sides of the Neumann noflow condition is imposed. Also, let . Thus the flow in the micromodel goes from left to right, and the solution evolves towards the stationary solution . Due to the incompatibility between initial and boundary data, we have high gradients of the solution close to for small .
We fix the final time and use uniform time stepping with . The solution of the micromodel depends on the number of cells in , with . For example, if , the domain is composed of cells. We denote the numerical solution of (2.8a) and (2.8b) at by where denotes the grid parameter for the micromodel.
For the macromodel, we use . Also, we use for when . If , then the boundary condition for is not prescribed.
3.2.2. Qualitative Comparison
In our first comparison we use parameters from Example 3.1. We solve the macromodel (4.2) for and the micromodel (2.8a) and (2.8b) for . The plots of and are in Figure 2, and a different view is shown in Figure 1(d).
(a)
(b)
(c)
(d)
The heterogeneous structure of is well visible from the behavior of the micromodel solution . Large gradients of solution are visible on cell boundaries due to the large difference between , , and .
However, the behavior in the fast region is well approximated by which envelopes , while envelopes well. This is very well seen in a side view in Figure 2(d).
Quite interesting is the behavior of . Since in our examples , no boundary condition at is imposed on . Thus, for small , evolves from its initial constant value of according to (1.3b) with the input from . The latter satisfies, however, the homogeneous boundary condition at . In particular, with constant , , we have Thus, for small and small , is away from , but, as time increases, its magnitude decreases proportionally to .
Next, we use parameters from Example 3.2 and plot them in Figure 3. In this example the conductivity is much larger than that of previous case, the local gradients in the micromodel are smaller, and the solutions to the macromodel achieve a smoother profile faster.
(a)
(b)
(c)
(d)
3.2.3. Quantitative Comparison—Fine Grid in MacroModel
Now we are ready to discuss a quantitative comparison between solutions to the micro and macromodels. For this, we set up a family of cases with , . We use uniform mesh with for now. See Section 3.2.4 for and Section 4 for adaptive gridding.
From [11, Theorem 18, page 208] and from our plots it follows that we can actually compare the micro and macrosolutions, as long as the characteristic functions are involved. While the results in [11] use twoscale convergence with , that is, rely on special periodic test functions, we compare and directly. We see that it is easy to do so only in the connected region .
We define a quantity to be used in comparison. We also define, for a fixed and some function , which will be applied to . Clearly is not a norm but is a useful quantity of interest.
In Tables 1, 2, 3, and 4 we present for different and . For each case we consider convergence of with .




From homogenization theory it is not clear what order of convergence should be expected. We observe that decreases linearly with for both like quantities in all cases.
Next we discuss the dependence of the results on conductivity and other data. Consider Example 3.2 as a reference example with its corresponding Table 1. Compare it now with that from Table 2 where is larger than that in the reference case to see how this influences the errors. As expected, and as shown in Example 3.3, changes by the same factor, and this means that the coupling between and is stronger which corresponds to faster flow across in the micromodel. However, the influence on the convergence of the micromacro difference, from Table 2 is rather weak.
Next consider Example 3.4 which uses a different geometry than that in the reference case, with thicker region . This produces a slightly smaller since the gradient of is smaller. Also, and are predictably different. The error quantity is influenced insignificantly.
The opposite effect is seen in Table 4 when is selected as thinner than in the reference Example 3.2.
3.2.4. Coarse Macrogrid and Parameter Grid
In the examples presented in Section 3.2.3 we used . This choice of compatible grids is convenient for visualization purposes. However, for the idea of the macromodel to be useful, the computational complexity of the macromodel needs to be lower than that of the micromodel. Since the macromodel is a system of two equations, its lower complexity can be only achieved if the grid in the macromodel is significantly coarser than that of micromodel. Below we show that one can use in the macromodel.
Example 3.6. We use data from Example 3.2 with . Consider a fixed microgrid with and a few cases of macrogrid with . Now we compare the micro and macrosolutions and , see Table 5. We also compare the run time of the macromodel to that of the reference micromodel which is 120 seconds.

We see that, as predicted above, the macromodel run time at is not competitive with that of the micromodel; this is exacerbated by the overhead of various adaptive procedures to be discussed below. However, when , the ability of the macromodel solutions to approximate those of micromodel appears reasonable and is associated with a much lower cost. This is true even when is by a factor of 10 coarser than resulting in a computational time decreased by two orders of magnitude.
Clearly, the macromodel can run much faster than the micromodel especially when adaptive nonuniform grids with are implemented, see Section 4.
Influence of Parameter Grid
Last we address the effect of the computed parameters , on the closeness of the macromodel solution to that of micromodel. These coefficients are precomputed as discussed in Section 3.1.1.
Example 3.7. Let , and let the micromodel grid be . Use . Compute numerical parameters , by approximating solutions to (3.4a), (3.4b), (3.4c), and (3.4d) and (3.3a), (3.3b), and (3.3c), respectively.
Example 3.7 and Table 6 show that the grid used for computing these parameters only mildly affects the quality of the solution, thus one can use coarse grid with confidence.

4. Error Estimates and TwoGrid Adaptivity
We present now the details of discrete formulation of (2.4) using independent grids for each component and adaptive gridding. The adaptive twogrid algorithm further reduces the computational cost of the macromodel. The formalism of the twogrid solution connects loosely to [37, 38] where their benefit was considered for a nonlinear scalar equation.
In what follows we assume that , that is, that the first component varies in space faster than the second. As a special case, this includes the case . We consider two triangulations , , which will be used for , respectively. In principle, these can be chosen independently. In fact, we allow the triangulations to vary in time and thus consider and the associated spaces as in (2.7). To avoid the loss of accuracy due to excessive interpolation and intergrid projections and because , we assume that is a refinement of the partition .
We need intergrid operators to handle two components that live on separate grids. Let be the interpolation operator and the projection defined by
Now we define the discrete solutions. At each time step , we seek , , which satisfy the discrete problem, that is, (2.4) restricted to the finite dimensional subspaces so that, for all , for all If , the system (4.2) is modified appropriately so that instead of , and we do not impose boundary conditions on this component. In fact, instead of (2.7) we define (It holds that but it is not necessary). It is straightforward that (4.2) is uniquely solvable.
4.1. A Priori Error Estimate
We have the following convergence result proved for the error in energy norm. In what follows and . We denote and define the energy norm for the product space
Assuming that is sufficiently smooth we have the following apriori estimate proven in the appendix.
Theorem 4.1. Let be the smooth solutions of (2.4) and (4.2), respectively. Then one has
where , , and are independent of .
Furthermore, if , then
Now we notice that (4.5) suggests the following choice of grid parameters , . If , then to balance the components of the error one can use . Similarly, if , then (4.6) implies that can be chosen to be of the order of . The use of coarse grid for the second component reduces the size of the linear system to be solved and decreases the computational time of the macromodel.
Also, the apriori estimates (4.5), (4.6) are global. One can further reduce the computational cost while maintaining accuracy by using local grid adaptivity guided by a posteriori estimators. This is pursued below.
4.2. APosteriori Error Estimates
Estimators for timedependent problems can be defined in many ways including the now classical spacetime element and adjoint approaches [27, 28]. In this paper we follow the residual estimator framework extending [15, 29] to the doublediffusion system using the ideas in [13, 39, 40] originally formulated for an elliptic system.
The estimator for (4.2) is composed of the temporal part which adapts the time stepping and the spatial part which guides the spatial grid adaptivity. We define and , where , and where, as usual (see [23]), one defines The element and edge residuals in , are given by Here denotes the jump of the flux across the edge of an element. Usually , . Also, if , then .
The scaling constants take into account the contribution of the exchange term and are defined so that the estimators remain robust when the coefficients of the problem change substantially: This a posteriori estimator works well with twogrid discretizations as shown in the examples to follow. It is well formulated also for when . We have the following result on reliability of the estimator proven in the appendix.
Theorem 4.2. Let , be the solution of (1.2a) and (1.2b), (2.3) and (4.2), respectively. Then
4.3. Adaptive TwoGrid Discretization and Implementation
Consider a fixed time step for which some triangulation is chosen, the solution is found, and the spatial error indicators are computed. To apply the twogrid spatial adaptivity we use the following refinement algorithm.
Adaptive TwoGrid Algorithm
(1)Select the triangulation(s) to be refined: If , refine only . If , refine only . Otherwise, refine both. (2)In , selected in Step 1, refine any for which . (3)Enforce the requirement that is a refinement of by adding extra elements as needed.
The steps in should be repeated a certain number of times until an ideal grid is found. See [41] for analysis of whether the iterative process of refinement and coarsening is, in general, convergent.
Example of TwoGrid Adaptivity
Now we show an example on how this adaptive algorithm works for the doublediffusion system. We use data and setup from Example 3.2, with and uniform time stepping with . Other data is as in Section 3.2.1.
Consider a fixed with a uniform triangulation with . Then apply times the steps in , see the solution and the refined meshes in Figure 4. The process works as expected: since at the time the solution has a high gradient near , the refinement occurs there. The refinement affects the grid for the first component only, because the gradients of the second component are not included due to . This can be compared to the apriori estimate in Theorem 4.1 from which we know the convergence is .
In Table 7 we provide the details on and compare the effectiveness of the local twogrid refinement by algorithm with uniform grid refinement. With adaptive refinement we get with unknowns, while the uniform refinement needs unknowns to achieve comparable .

(a)
(b)
(c)
(d)
Last, we describe the process with which , , may vary between time steps. For a fixed denote by the final triangulation obtained by the adaptive algorithm . For the new time step , we use as the initial triangulation. If the algorithm suggests that is modified, we need to project to the new grid so it can be used as initial data for the new step, but differs from their interpolation only in the elements where the triangulation is coarsened. Another way to deal with different triangulations , is to compute the a posteriori error estimator in a common refinement of the two triangulations , [29, 42], but this will not be pursued further.
Implementation
The decision to use two grids implies that we have to interpolate between the finite element spaces and , see (4.1). Suppose that , , is a basis for . Then the operator in (4.1) is represented by the matrix whose components are given by .
Consider now and compare the computational cost associated with onegrid and twogrid approaches. To solve (4.2) with onegrid we need to assemble the stiffness matrix , the mass matrix , and the block matrix
At each time step we solve the linear system with ; this is referred to as .
To solve (4.2) with two grids and two spaces , we need to assemble the stiffness matrix , the mass matrices , , the interpolation matrix , and the block matrix
To finish we solve the linear system which we refer to as .
Suppose now we solve (4.2) and wish to maintain for some tolerance over 1,000 time steps. With data as in Section 3.1.2 we run the simulations using one grid and find that we need . Using the twogrid approach we only need .
The computational time needed to assemble the matrices and solve the systems in our MATLAB implementation is displayed in Table 8. Clearly the assembly process takes more time for twogrid case than for onegrid discretization. However, the economy in the cost of solving the linear system makes up for that cost.

5. Conclusions
We compared the solutions to the micromodel to those of the macromodel and have shown that the latter is a good approximation to the former. This remains true also when a very coarse computational grid is used. Moreover, we established a linear convergence rate in function of the periodicity parameter in a certain quantity of interest which further shows that the doublediffusion model is an excellent approximation to the micromodel, at least for the considered scenarios.
To make the doublediffusion model computationally efficient, we proposed an algorithm for local grid adaptivity which allows each component to live on its own grid. The grid refinement is guided by an a posteriori error estimator for which we proved theoretical results.
Appendix
A. Proof of Theorem 4.1
Proof. The proof follows the techniques presented in [25, 26] for a scalar PDE. We present an outline of the proof for the doublediffusion system. By (4.1), if is a refinement of , we can rewrite (4.2) as
Next, for , as in [26] we define the elliptic projection of into via
Applying (A.2) for into the weak formulation (2.4), with arbitrary we arrive at
Next we define , , subtract (A.3) from (A.1), add the resulting equations, and use as test functions and to get