A semidiscrete Galerkin scheme for a coupled two-scale elliptic–parabolic system: well-posedness and convergence approximation rates

In this paper, we study the numerical approximation of a coupled system of elliptic–parabolic equations posed on two separated spatial scales. The model equations describe the interplay between macroscopic and microscopic pressures in an unsaturated heterogeneous medium with distributed microstructures as they often arise in modeling reactive flow in cementitious-based materials. Besides ensuring the well-posedness of our two-scale model, we design two-scale convergent numerical approximations and prove a priori error estimates for the semidiscrete case. We complement our analysis with simulation results illustrating the expected behaviour of the system.

coupled partial differential equations (PDEs) that explicitly encode two-scale interactions via transmission boundary conditions as well as production terms; see e.g. the PDE structures entering double or dual porosity models, models with distributed microstructures, fissured-media equations, as well as general two-scale models. Such models arise as descriptions of reactive flow through geometrically-structured porous media.
If the geometry of the porous media has a dual porosity structure, and hence, characteristic scales can possibly be separated, then PDE models with distributed microstructures are in theory able to describe the relevant multiscale spatial interactions like those occurring in gas-liquid mixtures. Now, the challenge shifts from the multiscale modeling to the computer implementation of multiscale models. Consequently, in this work we concern ourselves with the two-scale computability issue-complex systems of evolution equations acting on two spatial scales are notoriously hard to compute, especially if moving boundaries or stochastic dynamics are involved within e.g. the distributed microstructures. Combined with the so-called curse of dimensionality, this results in a computational problem of very high complexity.
In this paper, we discuss the case of an elliptic-parabolic coupling. We consider a coupled system of partial differential equations connected to multiscale descriptions of the evolution of the pressure arising in a compressible air-liquid mixture that distributes over two spatial scales (one called macroscopic, and one microscopic). This situation arises, for instance, in cementitious materials within concrete memberstypical composite porous materials where the amount of displaceable liquid is low and is practically trapped in the internal structure of the porous medium. The derivation of our particular model originates from applying a formal two-scale homogenization to a particular scaling of the level set equation coupled with Stokes equations for fluid flow (see [23] for details). Highlights of the more non-standard challenges of this system of equations are the two-scale coupling, the mismatch in structure between the two equations, i.e. the presence of a time derivative in the microscopic equation and its absence in the macroscopic one, as well as the nonlinear right hand side in the macroscopic equation. In order to tackle these challenges we perform most of our analysis on the finite element level, using techniques from e.g. [6,14,24].
If we assume the interface between air and liquid to remain fixed for a reasonable time span, then using homogenization techniques for locally periodic microstructures (compare with e.g. [5]) leads in suitable scaling regimes to a so-called two-pressure evolution systems. This system can be expressed as coupled elliptic-parabolic equations that describe the joint evolution in time t ∈ (0, T ) (T < +∞) of a parameter-dependent microscopic pressure Rρ(t, x, y) (where R represents the universal gas constant) evolving with respect to y ∈ Y ⊂ R d for any given macroscopic spatial position x ∈ Ω and a macroscopic pressure π(t, x) with x ∈ Ω for any given t. An illustration of the two-scale geometry we have in mind is depicted in Fig. 1.
In the final part of this manuscript, we discuss a numerical implementation of this scheme in the finite element library deal.ii [2]. Inspired by the Heterogeneous Multiscale Method framework (cf. e.g. [8]), we propose an implementation strategy that resolves the scale separation inherent in the two-scale structure of our problem. Fig. 1 The macroscopic domain Ω and microscopic pore Y at We consider the following problem, posed on two spatial scales Ω ⊂ R d 1 and Y ⊂ R d 2 with d 1 , d 2 ∈ {1, 2, 3} in the time interval t ∈ S := (0, T ) for some T > 0. Find the two pressures π : S × Ω → R and ρ : S × Ω × Y → R that satisfy: D∇ y ρ · n y = k(π + p F − Rρ) D∇ y ρ · n y = 0 i n S × Ω × Γ N , where Γ R ∪ Γ N = ∂Y , Γ R ∩ Γ N = ∅ and f , g are functions discussed in more detail below. We refer to (1)-(6) as (P 1 ). Note that (P 1 ) describes the interaction between a compressible viscous fluid (with density ρ) in a porous domain Ω, where the pores are partly filled with a gas that exerts an average (macroscopic) pressure π . The interaction between the fluid and the gas is determined by the right hand side of (1) and the microscopic boundary condition in (3), through the fluid-gas interface represented by Γ R . The mathematical problem stated in (P 1 ), contains a number of dimensional constant parameters: A (gas permeability), D (diffusion coefficient for the gaseous species), p F (atmospheric pressure) and ρ F (gas density). In addition, we need the dimensional functions k (Robin coefficient) and ρ I (initial liquid density). Except for the Robin coefficient k, all the model parameters and functions are either known or can be accessed directly via measurements. Even if the boundary Γ R is not accessible for measurements of parameters such as κ, this can be compensated by measuring on the boundary Γ N := ∂Y \Γ R . See e.g. [16]. We would like to point out that although we choose A and D constant, the analysis would be analogous for a system with A = A(x) and D = D(y), provided they satisfy certain suitable assumptions.
In this context, we prove existence and uniqueness of a discrete-in-space, continuous-in-time finite element element approximation and prove its convergence to the true solution of (P 1 ). The main results of this contribution are the well-posedness of the Galerkin approximation (Theorem 1), convergence rates for the approximation (Theorem 4), and the confirmation of the expected convergence rate with a numerical simulation (Sect. 5).
The choice of problem and approach is in line with other investigations running for two-scale systems, or systems with distributed microstructures, such as [15,17,21]. The reader is also referred to the FE 2 strategies developed by the engineering community to describe the evolution of mechanical deformations in structured heterogeneous materials; see e.g. [12] and references cited therein. Other classes of computationally challenging two-scale problems are mentioned, for instance, in [22], where the pore scale model has a priori unknown boundaries, and in [11] for a smoldering combustion scenario. This paper continues an investigation started in related works. In [16], we study the solvability issue and derive inverse Robin estimates for a variant of this model problem. Two-scale Galerkin approximations have been derived previously for related problem settings; see e.g. our previous investigations [4,15,18,19].
The rest of this paper is structured as follows. In Sect. 2, we discuss the technical concepts and requirements we need before starting our analysis. Then, in Sect. 3, we show the Galerkin approximation is well-posed and converges to the weak solution of the original system. In Sect. 4, we prove a priori convergence rates for the Galerkin approximation. Next, in Sect. 5, we propose a fully discrete scheme, an implementation of said scheme, and approximation errors of the finite element solutions on subsequently refined grids. Finally, in Sect. 6, we conclude this paper and provide an outlook into future research.

Weak solutions
We look for solutions to (P 1 ) in the weak sense. This is motivated by the fact that the underlying structured media can be of composite type, allowing for discontinuities in the model parameters. However, already at this stage it is worth mentioning that the solutions to (P 1 ) are actually more regular than stated, i.e. with minimal adaptations of the working assumptions the regularity of the solutions can be lifted so that they turn out to be strong or even classical. We will lift their regularity only when needed.

Assumptions
We introduce a set of assumptions that allows us to ensure the weak solvability and approximation of (P 1 ).
(A 1 ) The domains Ω and Y are convex polygons. where C Ω denotes the constant in Poincaré's inequality for H 1 0 (Ω) [see (18) (1) satisfies the following conditions: and where θ is small enough to satisfy all of the following: Note that this implies θ < 1.
An example of a functional g satisfying (A 6 ) is where L is a linear map. The fact that g defined as above satisfies (A 6 ) is a consequence of the interpolation-trace inequality [see (17) below].

Technical preliminaries
The rest of this section introduces the notation of the functional spaces and norms used in the paper. Let Tr 1 : H 1 (Ω) → L 2 (∂Ω) denote the (macroscopic) trace operator defined as and let Tr 2 : L 2 (Ω; H 1 (Y )) → L 2 (Ω × ∂Y ) denote the microscopic trace operator defined as Let f , g : D → R. Then the Lebesgue and Sobolev norms are defined as follows: with ∂ α f denoting derivatives in the weak sense. Furthermore, for L 2 (D) and H k (D) we have the following inner products.
Moreover, we use H 1 0 (D) to denote the following function space: Let B be a Banach space with norm · B . Then u belongs to Bochner space L 2 (S; B) if it has a finite L 2 (S; B) norm, defined as follows: An introduction to the concepts of Lebesgue and Bochner integration as well as on inner products and norms can be found in many functional analysis textbooks (e.g. [1]).

Auxiliary results
For the benefit of the reader, we collect a number of well-known results that we will use in the paper.
Lemma 1 (Young's inequality) Let E ⊆ R d be a measurable set and u, v ∈ L 2 (E). For any > 0 there holds It is well-known that Tr 2 defined between the function spaces specified above is a bounded linear operator. Thus, quantities of the type u L 2 (Ω×∂Y ) are well defined for u ∈ L 2 (Ω; H 1 (Y )).
Lemma 3 (Poincaré's inequality) There exists a constant C Ω depending only on Ω such that for all u ∈ H 1 0 (Ω).
Then the embedding of W into L 2 (S; B) is compact.
We refer the reader to [3] for the original proof of the statement.

Well-posedness
In this section we prove that (P 1 ) has a weak solution by approximating it with a Galerkin projection. We show the projection exists and is unique, and proceed by proving it converges to the weak solution of (P 1 ). First, we introduce the necessary tools for defining the Galerkin approximation. We use one mesh partition for each of the two spatial scales. Let P k be the space of polynomials up to degree k. Let B H be a mesh partition for Ω consisting of simplices. We denote the diameter of an element B ∈ B H with H B , and the global mesh size with H := max B∈B H H B . We introduce a similar mesh partition K h for Y with global mesh size h := max K ∈K h h K .
Our macroscopic and microscopic finite element spaces V H and W h are, respectively: We note that this approach is easily extensible to x-dependent microscopic domains. Let N 1 and N 2 denote the sets of degrees of freedom in B H and K h , respectively. Let span (ξ i ) = V H and span (η k ) = W h , and let α i , β ik : S → R denote the Galerkin projection coefficient for the ith and ikth degree of freedom, respectively. We introduce the following finite-dimensional Galerkin approximations of the functions π and ρ: Reducing the space of test functions to V H and W h , we obtain the following discrete weak formulation: find a solution pair and for any ϕ ∈ V H and ψ ∈ V H × W h and almost every t ∈ S. Furthermore, π H (0, x) = π H I (x) (see (9) and is a Galerkin approximation of (π I , ρ I ). These concepts lead us to the first theorem.
Theorem 1 (Existence and uniqueness of the Galerkin approximation) There exists a unique solution (π H , ρ H ,h ) to the system in (21) and (22).

Proof
The proof is divided in three steps. In step 1, the local existence in time is proven. In step 2, global existence in time is proven. Step 3 is concerned with proving the uniqueness of the system.
Since the system of ordinary differential equations in (25) is linear, we are able to explicitly formulate the solution representation for β ik with respect to α i . Let α i be given, and let Q and E denote matrices given by: Then β ik can be expressed as Substituting (31) in (25) results in the expression: ik be two function pairs that satisfy (31). Then it holds that with c α , c β defined as Now, we derive a time-dependent continuity estimate for sufficiently small t. Again picking two pairs (α * (t), β * (t)) and (α * * (t), β * * (t)) (not necessarily the same as in (33)): Here, the size of valid t is independent of the initial data. Using (35) we obtain a Lipschitz bound on all F i in the interval [0, τ ] for any choice of τ < t: Choosing τ small enough to satisfy c α +c β Cτ < 1 makes F a contraction on [0, τ ]. By Banach's fixed point theorem, it follows that the equation F(α(t), β(t)) = α(t) has a solution for α in L 2 (S). Substitution of α(t) into (31) leads to the corresponding β. Existence of π H and ρ H ,h follows directly.
Step 2: global existence of solutions to (21) and (22) Given that interval n has local solution β ([(n − 1)τ, nτ ]), we can obtain values β(nτ ), β (nτ ), α(nτ ) as initial values to the local system on interval n + 1, and show existence of a solution. This way, we are able to construct a solution satisfying (21) and (22) everywhere on S.
Step 3: uniqueness of solutions to (21) and (22): We decouple the system and use a fixed point argument to show that this system has a globally unique solution in time.
Note that showing the stability of the finite element approximation with respect to data and initial conditions follows an analogous argument. The proof is omitted here. The remaining part of this section is devoted to proving that the system in (21) and (22) converges to the solution of the Galerkin projection converges to the weak solution of (P 1 ). Our first aim is to derive standard energy estimates for the discrete solution (π H , ρ H ,h ).
Lemma 5 (Standard energy estimates) Let (π H , ρ H ,h ) be a solution to (21) and (22). Then we have the following energy estimates and and where C is independent on h and H , while it depends on the model parameters and the geometry of the domains.
We proceed with (42). Using Cauchy-Schwarz' inequality and (17) to the righthand side of (42) Then, we add to both sides of (42) a term D ρ H ,h 2 L 2 (Ω×Y ) to get After applying Young's inequality with the small parameter ε > 0, we get By applying Grönwall's inequality we obtain the desired estimates: with Finally, testing (22) where C is independent of H , h and ρ H ,h .
We shall need a bound for ∂ t π H and also an estimate for ∂ t ρ H ,h that is sharper than (40).

Remark 4
We will show below that the pair (π, ρ) provided by Proposition 1 is a weak solution to (P 1 ). However, the convergence statements of Proposition 1 are not strong enough to allow us to pass to limit in (21) and (22), due to the nonlinear term f (π H , g(ρ H ,h )). The next two lemmas provide additional regularity that will help us strengthen the convergence, and also be useful in Sect. 4. (π H , ρ H ,h ) be the solution to (21) and (22) for H , h > 0. Then

Lemma 7 Let
and Proof To prove (55) We introduce the directional finite difference for i ∈ {1, . . . , d 1 }. Let 0 < λ < δ and test (22) with which gives us: Because of the properties of the support of ζ , it holds that for any Applying the property in (59) to (58) yields leading to Using Young's inequality combined with the inequality, we estimate the third term of (61) as follows: Now, combining (62) with (61), we obtain Using Grönwall's inequality, we conclude that D λ i ρ H ,h ∈ L 2 (Ω × Y ), and by letting λ → 0, we obtain With the results obtained above, we are ready to state and prove the first two main result of this paper.
Proof Let (π H j , ρ H j ,h j ) and (π, ρ) be provided by Proposition 1. (21) and (22) and using the convergence statements of Proposition 1, it follows that (π, ρ) solves (7) and (8) if we can show that for all ϕ ∈ H 1 0 (Ω). Denote by By the Rellich-Kondrachov theorem, we have the compact embeddings and Hence, Aubin-Lions' lemma (Lemma 4) gives that }) that converge strongly in the spaces L 2 (S × Ω) and L 2 (S × Ω × Y ) respectively, and these strong limits must coincide with π and ρ. By continuity of f and g, (64) follows, and the theorem is proved.

Theorem 3 (Uniqueness of the weak solution)
The weak solution to problem (P 1 ) is unique.
Since θ < 2C Ω /A < 1 it also holds that θ < A/(2C Ω ) and thus Whence, By using Young's inequality with parameter (to be determined below) and (67) at the right-hand side of (66), we obtain Applying the interpolation-trace inequality with parameter 2 , we get Hence, Now take 0 such that then we get by Grönwall's inequality we obtain in other words ρ 1 = ρ 2 a.e., and it follows from (67) that π 1 = π 2 a.e.

Convergence rates for semidiscrete Galerkin approximations
In this section, we obtain convergence rates of the numerical approximations (21) and (22). The following argument is largely based on standard arguments from e.g. [14], adapted to multiscale systems.
Let R H : H 1 (Ω) → V H and R h : L 2 (Ω; H 1 (Y )) → L 2 (Ω; W h ) be the microscopic and macroscopic Ritz projection operator respectively. The projections R H r (x) and R h s(x, y) are defined such that: s(x, y)) · ∇ y w = 0, for all v ∈ V H and w ∈ L 2 (Ω; H 1 (Y )).

Lemma 8 (Projection error estimates)
Then there exists strictly positive constants γ l (with l ∈ {1, 2, 3, 4}), independent of h and H , such that projections R h π and R H ρ that satisfy Proof (68) and (69) are standard Ritz projection error estimates. For details on the proof, see for instance [14,24]. Specific to this context, (70) is a two-scale estimate which accounts for the presence of the microscopic Robin boundary condition (3) and therefore requires some tuning. See e.g. [18] for similar estimates. Here, we only present the proof of (70).
Applying the Ritz projection estimates (68) and (69), we obtain the following bound: Using Friedrich's inequality ϕ L 2 (Ω;H 2 (Y )) ≤C Δϕ L 2 (Ω;L 2 (Y )) =C ω L 2 (Ω;L 2 (Y )) for some C and choosing ε < c we obtain (73) yields: Finally, we can derive (70) as follows: By applying Lemmas 2 and 8, we can finally obtain the desired convergence rates. Let us denote the errors of the Galerkin projection as e π := π − π H , Theorem 4 (Convergence rates) Let (π H , ρ H ,h ) be a solution to (7) and (8) Proof By testing (7) and (21) with φ H ∈ V H and subtracting the equations, we obtain the identity for any φ H ∈ V H . Let ϕ H ∈ V H be arbitrary. By (78) applied with φ H = ϕ H − π H , we obtain By the inequality above, Cauchy-Schwarz' inequality and the triangle inequality, we get Applying Cauchy-Schwarz inequality to the last two terms above and denoting , we obtain Rearranging the above inequality gives Using (18) and (A 3 ) (i.e. 2C Ω /A < 1), we obtain As in the proof of Theorem 3, the assumptions (A 5 ) and (A 6 ) yield the estimate Using (18) again, we get By the assumptions, i.e. 1 + 4θ 2 < 3/2 and C Ω /A < 1/2, we get By (69), we have π − ϕ H Using a standard duality argument as in [14], we get e π L 2 (Ω) ≤ C H e π H 1 0 (Ω) . We proceed to demonstrate that e ρ L 2 (Ω×Y ) = O(H 2 + h 2 ). Write We bound ψ by using Lemma 8: and bound θ from (80) using the formulation: for all ϕ ∈ V h we have that Substituting ϕ = θ in (82) yields: Dividing the left and right hand side of (83) by θ , we obtain: Because of (A 4 ), the Galerkin projection error of the initial condition satisfies: Combining (81) and (84) proves the desired estimate in (77).

Implementation
In this section, we discuss a time-discretized version of (21) and (22) i.e., a fully discretized system, and provide details and performance results of the implementation of this system.

Setup
We implement the finite element formulation of this problem using deal.ii [2], a C++ library that computes numerical approximations to finite element problems on quadrilateral meshes.
To account for the scale separated multiscale structure, we implement this system using a heterogeneous multiscale method (see e.g. [8] for an introduction). Alternative multiscale finite element structures are the FEM 2 method [12] and the multiscale finite element method (MsFEM [7]). Both of these frameworks generally require a formulation in which the size relation ε between the macroscopic scale and the microscopic scale must be resolved. Since our problem is completely scale-separated, we opt for a heterogeneous multiscale method.
We build the microscopic systems by assigning a microscopic grid for every degree of freedom on the macroscopic grid. Since we use nodal basis functions, every degree of freedom corresponds to a single physical location in the macroscopic grid. By allowing the microscopic systems to correspond with degrees of freedom, by integrating the microscopic finite element functions on the finite element domain, we obtain a finite element function on the macroscopic domain, for which we can use classical finite element techniques.
deal.ii has no specific support for multiscale problems. However, we can build upon its structure to create new components that can deal with objects like multiscale functions and multiscale solutions.
For our multiscale implementation, we need not use a separate instance for each macroscopic degree of freedom. Specifically, we assume the same microscopic grid and triangulation for each microscopic instance. This allows us to reuse and share microscopic data structures throughout the simulation.

Manufactured system
We test the quality and correctness of the scheme and its implementation by simulating a more general problem, for which we can manufacture solutions. We compute the solution of this problem on subsequently finer meshes, and check if convergence rates are according to expectations. The manufactured problem (P M ) is defined as follows: where, in our test scenario, p, q, r , s, u are chosen such that they lead to a π and ρ of which we know the explicit form. Note that if we let p = q = r = s = u = 0, (P M ) reduces to (P 1 ). Additionally, f is defined in accordance to the first example provided in Sect. 2: f (π, g(ρ)) := θ min |π |, |π | α min (1, |g(ρ)|) , where θ is defined in Table 1, α = 0.5 and g is defined as  [20] and has a value of C Ω = 2 √ 2 π , where π in this expression represents the mathematical constant.

Time discretization
We discretize the microscopic equation in time with an implicit Euler scheme, while we use a Picard-like iteration scheme for the macroscopic equation. The Picard-like iterations avoid approximating the nonlinear term f via e.g. Newton's method, and the implicit Euler scheme ensures more stability in the microscopic equation without significant complications, since this equation is linear.
We discretize time domain S with time steps 0 < t 1 < · · · < t N T . Let τ n be the time step size in time step n. Then, given data from time step n − 1, we compute the approximations in time step n by solving: p represents the subsequent observed order of convergence of the finite element error, q represents the subsequent observed order of convergence of the error of its gradient The microscopic degrees of freedom (mDoFs) are summed over all microscopic grids. p represents the subsequent observed order of convergence of the finite element error, q represents the subsequent observed order of convergence of the error of its gradient

Results
We fix the parameters according to the values presented in Table 1.
We approximate the solutions to (P M ) with piece-wise linear basis functions, and compute the cell contributions with a third order Gaussian quadrature rule. We test the implementation by solving (P M ), choosing p, q, r , s such that solutions π, ρ become: Running this simulation for increasingly smaller τ n (to account for the time discretization error), h and H , yields the errors presented in Tables 2 and 3. The values p i and q i for i = 1, 2 represent the (subsequent) observed order of convergence of the finite element error and the gradient error, respectively. For example: refining the macroscopic grid results in an observed error of size e π L 2 (Ω) = O (H p 1 ).
Ignoring the time discretization as an error source, we observe the finite element error behaves according to the theory: e π L 2 (Ω) + e ρ L 2 (Ω;L 2 (Y )) = C h 2 + H 2 + h.o.t., thereby confirming Theorem 4. Noteworthy is the fact that the microscopic error accounts for a large part in the total error. This might be due to the fact that the macroscopic equatio lacks a time derivative. For t = 0.5, a representation of the macroscopic solution and a graphical representation of the microscopic solutions are represented in Figs. 2 and 3, respectively.

Conclusion and future work
We constructed a semidiscrete Galerkin approximation of our semi-linear ellipticparabolic two scale system and showed that this approximation is well-posed. Furthermore, the obtained sequence of Galerkin approximants based on finite elements converges in suitable spaces to the weak solution to the continuous system. Under additional regularity assumptions, we derived a priori rates of convergence of our approximation to the target weak solution. Finally, we implemented the fully discrete system in deal.ii, tested the convergence rates in practice, and observed the behaviour of the system for a certain set of parameters. We found convergence behaviour according to our analysis. Using this setup, we are now able to deal with truly scale-separated, two-scale problems, using a method fitting in the HMM framework.
As a natural next step, in a forthcoming work a numerical analysis study will address the quality of the fully discrete two-scale Galerkin approximation as well as improved numerical implementations of the method so that the proven convergence rates can be confirmed and variants of macroscopic refinement strategies can be tested, possibly a two-scale elliptic-elliptic PDE system (obtained by letting t → ∞ in our elliptic-parabolic formulation).
Acknowledgements Open access funding provided by Karlstad University. The authors acknowledge fruitful discussions with Prof. M. Asadzadeh (Chalmers University of Technology, Gothenburg, Sweden). OR thanks Dr. D. Tagami (Kyushu University, Japan) for valuable feedback and acknowledges partial support from Kungl. Vetenskapsakademien, Sweden. ML, AM and OR thank Dr. O. Lakkis and Dr. C. Venkataraman (both with University of Sussex, UK) for the intensive interactions during the Hausdorff Trimester Program "Multiscale Problems: Algorithms, Numerical Analysis and Computation" (Bonn, January 2017) and discussions at the University of Sussex. The authors would also like to thank the reviewers for their valuable comments that were a great help to improve the paper.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.