A P ] 1 4 A ug 2 01 8 A semidiscrete Galerkin scheme for a two-scale coupled elliptic-parabolic system : well-posedness and convergence approximation rates

In this paper, we study the evolution of a gas-liquid mixture via 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. Besides ensuring the well-posedness of our two-scale model, we design two-scale convergent numerical approximations and prove a priori error estimates and propose an a posteriori error estimator. Finally, we propose a macroscopic mesh refinement strategy that ensures a redistribution of the local macroscopic errors until an overall error reduction is achieved.


Introduction
This work is concerned with the design and approximation of evolution equations able to describe multiscale spatial interactions in gas-liquid mixtures. The long term goal and ultimate target is to set the foundation for a rigorous mathematical justification of Richards-like equations. Upscaled equations for the motion of flow in unsaturated porous media are chosen in a rather ad hoc manner by various engineering communities. The main issue is that one lacks a rigorous derivation of Darcy's law for such flows altering between compressibility and incompressibility. Therefore, the reliability of macroscopically imposed laws which are exclusively based on first principles and perhaps on single-scale fitting strategies is limited. One the other hand, the fully saturated case is clear. We refer the reader to Chapter 1 of [14] for a derivation of the Darcy's law in the saturated case using periodic homogenization arguments. Regarding the context of Darcy's law for unsaturated flows, related work is reported, for instance, in [26], [20] and [8].
If the geometry of the porous media has a dual porosity structure, and hence, characteristic scales can possible be separated, then PDE models with distributed microstructures are in theory able to describe the relevant multiscale spatial interactions occurring in gas-liquid mixtures. Now, the challenge shifts from multiscale modeling to the implementation of multiscale models. Consequently, we are concerned with the two-scale computability issue -complex systems of evolution equations acting on two spatial scales are notoriously hard to approximate especially if moving boundaries or stochastic dynamics are involved within e.g. the distributed microstructures.
In this paper, we consider a coupled system of partial differential equations describing the evolution of the pressure of a compressible air-liquid mixture on two spatial scales, when the amount of liquid is low and trapped in the internal structure of a 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 [28]).
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 e.g. [6]) leads in suitable scaling regimes to a so-called two-pressure evolution systems. This system can be expressed as a coupled elliptic-parabolic systems 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 Figure 1.
Y Ω x Figure 1: The macroscopic domain Ω and microscopic pore Y at x ∈ Ω.
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 (1)-(6) (referred to as (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.
If the Neumann part of the boundary Γ N := ∂Y \ Γ R is accessible via measurements, then the inaccessibility of the boundary Γ R can be compensated for in such a way that parameters like k entering two-scale transmission conditions can be identified (compare [19]).
In this context, we prove the existence and uniqueness of a discrete-in-space, continuous-in-time finite element element approximation and prove convergence of this approximation of (P 1 ). The main results of this contribution are the wellposedness of the Galerkin approximation (Theorem 1), convergence rates for the approximation (Theorem 2), and controlled error estimators which can then be used to refine the grid (Theorem 3).
The choice of problem and approach is in line with other investigations running for two-scale systems, or systems with distributed microstructures, like [18,21,25]. The reader is also referred to the FEM 2 strategies developed by the engineering community to describe the evolution of mechanical deformations in structured heterogeneous materials; see e.g. [16] and references cited therein. Other classes of computationally challenging two-scale problems are mentioned, for instance, in [27], where the pore scale model has a priori unknown boundaries, and in [15] for a smoldering combustion scenario. This paper continues an investigation started in related works. In [19], we study the solvability issue and derive inverse Robin estimates for a variant of this model problem. Twoscale Galerkin approximations have been derived previously for related problem settings; see e.g. our previous attempts [23], [22], [5], and [18]. Ref. [18] stands out since it is for the first time that the issue of feedback estimates is put in the context of computational efficiency of PDEs posed on multiple scales. Unfortunately, the obtained theoretical estimates in loc. cit. are not computable. This aspect is addressed here in Theorem 3.
The rest of this paper is structured as follows. In Section 2, we discuss the technical concepts and requirements we need before starting our analysis. Then, in Section 3, we show the Galerkin approximation is well-posed and converges to the weak solution of the original system. In Section 4, we prove a priori convergence rates for the Galerkin approximation. Next, in Section 5, we design an error estimator, prove upper and lower bounds with respect to the true error, and propose a macroscopic mesh refinement strategy. Finally, in Section 6, we conclude this paper and provide an outlook into future research.
2 Concept of weak solution, assumptions and technical preliminaries

Weak solutions
Essentially, 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 adaptions 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.
(A 2 ) All model parameters are positive; in particular D, R, p F , κ.
The value of C f and C π is given in Section 3. (1) satisfies the following structural conditions: (i) f is once continuously differentiable in π and ρ.
(ii) f (s, r) is a contraction in s for all r.
(iv) There exists a θ > 0 such that f (s, r) = 0 for all s > θ. Before moving on, we remark that the formal shape of the right hand side in (1) must be of the form f (π, g(ρ)) for some g : S × Ω × Y → S × Ω.

Technical preliminaries
The rest of the section introduces the notation of the functional spaces and norms used in the paper.
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: and H −1 (D) to denote the dual space of (13), equipped with the norm Let B be a Banach space with norm || · || B . Then u belongs to Bochner space L 2 (S; B) if its norm is finite, 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 any functional analysis textbook (e.g. [1]).
By π I ∈ H 1 0 (Ω) we denote the solution of −A∆ x π = f (π, ρ I ) in Ω, (16) is a stationary elliptic equation giving access to the value of π from (1) at time t = 0. Finally, we introduce several constants: c i refers to constants from the interpolation-trace theorem (see Lemma 4) and c p to constants arising from Poincaré's inequality. Moreover, we define the following two constants:

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 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∈BH 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: Let ξ B BH := span(V H ) and η K K h := span(W h ), and let α B , β BK : S → R denote the Galerkin projection coefficient for a patch B and B ×K, respectively. We introduce the following finite-dimensional Galerkin approximations of the functions π and ρ: where we clamp α B (t) = 0 for all B ∈ B H with ∂B ∩ Ω = ∅ to represent the macroscopic Dirichlet boundary condition.
Reducing the space of test functions to V H and W h , we obtain the following discrete weak formulation: find a solution pair (π H (t, x), ρ H,h (t, x, y)) ∈ L 2 (S; and for any ϕ ∈ V H and ψ ∈ V H × W h and almost every t ∈ S. These concepts lead us to the first proposition.
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.
We introduce an integer index for α B (t) and β BK (t) to increase the legibility of arguments in this proof. Let N 1 := {1, . . . , |B H |} and N 2 := {1, . . . , |K h |}. We introduce bijective mappings n 1 : N 1 → B H and n 2 : N 2 → K h , so that each index j ∈ N 1 corresponds to an element B ∈ B H and each index k ∈ N 2 corresponds to a K ∈ K h .
Since the system of ordinary differential equations in (21) 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 (26) in (21) results in the expression: (26) given some α 1 (t), α 2 (t). Then it holds that: with c α , c β defined as Now, we derive a time-dependent continuity estimate for sufficiently small t. Again picking β 1 (t) and β 2 (t) (not necessarily the same as in (28)): Using (30) 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 (26) leads to the corresponding β. Existence of π H and ρ H,h follows directly.
Step 2: global existence of solutions to (18) - (19): We cover time interval S into N intervals of length at most τ such that S ⊆ n ((n − 1)τ, nτ ]. From the arguments in the previous paragraph it is clear a solution exists on the first interval [0, τ ]. This allows us to provide an induction argument for the existence of a solution on interval n: 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 (18) -(19) everywhere on S.
Step 3: uniqueness of solutions to (18) - (19): 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 (18)- (19) converges to the solution of the Galerkin projection converges to the weak solution of (P 1 ) (as stated in Proposition 2). To this end, we first formulate the lemmata that help us prove this statement.
Then the embedding of W into L 2 (S; B) is compact.
We refer the reader to [2] for the original proof of the statement.
Proof. We use a weak maximum principle according to Stampacchia ([29]). Consequently, we test the weak formulation with ϕ = (π H − M 1 ) The left hand side of (34) can be manipulated as: which can be bounded with the right hand side of (34): Proceeding similarly with (19): Add (36) and (37) Lemma 3 (Regularity lift). Let (π H , ρ H,h ) be a solution to (18)- (19). Then it must hold that Proof. Testing (18) with ϕ = π H and (19) with ψ = ρ H,h yields identities and We recall the embedding which implies that there exists a c E such that for all u ∈ L 2 (Ω; H 1 (Y )). Using Cauchy-Schwarz' inequality and (43), we can bound the right hand side of (41) as Then, we add to both sides of (41) a term D||ρ H,h || 2 L 2 (Ω×Y ) to get 1 2 After applying Young's inequality with the small parameter ε > 0, we get By applying Grönwall's inequality we obtain the desired estimates: with To obtain a bound on π H , we test (18) with π H and then use (A 3 ) and (A 5 ): This yields the upper bound With these lemmata we are ready to state and prove a first convergence result.
Then, by choosing we satisfy the requirements of Aubin-Lions' lemma to get convergence of π H in L 2 (Ω).
Concerning ρ H,h , again we use Aubin-Lions' lemma to prove convergence, this time for the following spaces: Note that ρ H,h ∈ L 2 (S; B). To conclude the argument, what remains to show are the following steps: and ∇ x ρ H,h ∈ L 2 ((S; L 2 (Ω; H 1 (Y ))).
We start with the estimates that provide us (52) and (55) and then we handle (56).
Let f π and f ρ denote the partial derivatives of f to respectively π and ρ. We introduce u H := ∂ t π H and v h := ∂ t ρ H,h − ρ I . We differentiate (P 1 ) with respect to t and obtain the following system By multiplying (58) with v h and integrating the result over Ω × Y , we obtain an equation that can be bounded similarly to (44). This gives: From (63), by integration in time we obtain By multiplying (57) with v h and integrating the result over Ω, we obtain which, by combining this inequality with Poincaré's inequality for u H , yields Thus, we obtain an upper bound that holds in the interior of Ω, say Ω δ .
To handle the integral estimates on the boundary layer Ω δ \ Ω, we recall (19). For any Ω δ ⊂ Ω, by testing with ψ = ∂ t ρ H,h and integrating over the time domain, we obtain the identity Conveniently rearranging the terms of (67) yields: Now we can extend the bound in (66) to hold on the entire domain Ω, i.e.: (71) Finally, to obtain (56) we adapt an interior regularity argument from [11], Chapter 6. We let Ω δ ⊂⊂ W ⊂⊂ Ω and define a smooth cutoff function ζ : Ω → [0, 1] satisfying We introduce the directional finite difference for i ∈ {1, ..., d 2 }. We let λ be small and we test (19) with which gives us: Because of the properties of the support of ζ, it holds that for any f ∈ Ω Applying the property in (74) to (73) yields leading to Using Young's inequality combined with the inequality, we estimate the third term of (76) as follows: Now, combining (77) with (76), we obtain the required estimate for (56): Using Grönwall's inequality, we conclude that D λ i ρ H,h ∈ L 2 (Ω × Y ), and by letting λ → 0, we obtain With the newly found estimates (71), (79) and (64), we are able to apply Lemma 1 and we obtain that which proves: for h, H → 0.
The preliminary work allows us to state the first main result of this paper.
Proof. The proof of this theorem is a direct result of Proposition 1 and Proposition 2.

Convergence rates for semidiscrete Galerkin approximations
In this section, we obtain convergence rates of the numerical approximations (18) - (19). The following argument is largely based on standard arguments from [17], adapted to multiscale systems.
Proof. We omit the proof and refer to [13].
Proof. The proof follows from applying Young's inequality with a small parameter ε to the standard trace inequality.
Let R h and R H be the microscopic and macroscopic Ritz projection operator respectively.
Proof. (81) and (82) are standard Ritz projection error estimates. For details on the proof, see for instance [30] and [17]. Specific to this context, (83) 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. [22] for similar estimates. Here, we only present the proof of (83).
Applying the Ritz projection estimates (81) and (82), 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 (86) yields: Finally, we can derive (83) as follows: By applying Lemma 4 and Lemma 5, we can finally obtain the desired convergence rates. Let us denote the errors of the Galerkin projection as e π := π − π H , Theorem 2 (Convergence rates). Let (π H , ρ H,h ) be a solution to (7)- (8). Then the following statement holds: there exists constants M 1 , M 2 > 0 independent of h and H, such that Proof. By testing (7) with π and π H and testing (18) with π H , we obtain the following identities: +A which, summed up, results in the following identity: Let ϕ H ∈ V H be arbitrary. By applying the identity and using (95), we obtain the following estimate (omitting L 2 norm indications in their respective spaces for clarity): Moving the first and last term of (97) to the left hand side, we obtain the following inequality, which can be further bounded as follows: Finally, by compensating the small terms and using the finite element approximation property: applied to (98) we obtain the inequality with C a generic constant independent of H. Continuing, from (19), we get We bound ψ by using Lemma 5: and bound θ from (101) using the formulation: for all ϕ ∈ V h we have that ∂ t θ, ϕ L 2 (Ω;L 2 (Y )) + D ∇θ, ∇ϕ L 2 (Ω;L 2 (Y )) = − R h ∂ t ρ, ϕ L 2 (Ω;L 2 (Y )) − D ∇ρ, ∇ϕ L 2 (Ω;L 2 (Y )) , Substituting ϕ = θ in (103) yields: Dividing the left and right hand side of (104) by ||θ||, we obtain: Because of (A 4 ), the Galerkin projection error of the initial condition satisfies: Combining (102) and (105) proves the desired estimate in (90).

A posteriori refinement strategy
In this section we develop a computable error estimator which we will use to refine the finite element grid and obtain a lower overall error for the macroscopic equation (18). In this strategy, we aim for reliable error estimators, i.e. estimators which provide an upper and lower bound on the error.
There is a difference in usability between (a priori) error bounds and (a posteriori) error estimators. Where error bounds guarantee the upper bound of the error, often they are not sharp, not computable, and mainly useful for proving well-posedness of the numerical approximation. Error estimators, on the other hand, should be computable quantities approximating the true value of the error, and preferably provide both an upper bound and an lower bound on the error. An upper bound is required to guarantee the maximum error satisfies a certain tolerance. A lower bound makes sure that the error estimator does not overestimate the true error too much, ensuring the efficiency of the refinement strategy.
For a review on the different strategies in error control, we refer to e.g. [4] and [12]. The line of arguments we present, is based on [32].
In this section, we only describe how to obtain error estimators for the macroscopic equation. Strategies to obtain error estimators for parabolic equations can be found in e.g. [3], [24], [10] and [9].

Mesh-related notation
We assume the mesh partition B H with diameter H as defined in Section 3. For each element B ∈ B H , we denote the set of vertices with V(B) and the set of edges with E(B). The complete set of edges is denoted with E. Where necessary, we differ between vertices (edges) in Ω and on ∂Ω by denoting them as B H,Ω (V H,Ω ) and B H,∂Ω (V H,∂Ω ), respectively. To denote patches in Ω with certain structures, we use the following symbols: • ω B denotes the union of all elements that share an edge with B.
•ω B denotes the union of all elements that share a point with B.
• ω E denotes the union of all elements adjacent to E.
•ω E denotes the union of all elements that share a point with B.
• ω x denotes the union of all elements that have x as a vertex.
Furthermore, for any E ∈ E, n E denotes the unit vector orthogonal to E and J E (v) denotes the jump across E of some piece-wise continuous function v in the direction of n E .
Finally, for legibility we use the following notation to refer to norms on elements or edges.

Auxiliary results
For any x ∈ V(B), we denote the piece-wise linear basis function that takes value 1 in x and 0 in the other nodes by λ x . This allows us to define the following cutoff functions for any B ∈ B H : and for any E ∈ E Here, the coefficients α E , α B are chosen such that: With the cutoff functions defined in (108) and (109), we can obtain the following bounds for any v ∈ H 1 (B) and ϕ ∈ H 1 (E): for any element B ∈ B, edge E ∈ E.
In order to obtain error estimates on both elements and edges, we introduce a quasi-interpolation operator I H defined as This allows for the following estimates: We omit the derivation of (111) and (113) and refer the reader to [7] and Section 3.1 of [31], respectively.

Macroscopic error estimator
The error estimator is composed from the jump discontinuities on each edge and the residuals in each patch. Formally defined, for any B ∈ B H and E ∈ E, we define which can be combined into a complete residual operator R : V H → R, defined implicitly in the following identity: (115) must hold for all ϕ ∈ V H .
Recall e π = π − π H denotes the finite element error of π. There is an equivalence between the norm of residual R and the norm of true error. Subtracting (18) from (7), we obtain Picking a suitable c * and c 1 , by applying Poincaré's inequality and Cauchy-Schwarz' inequality shows that: On the other hand, recalling (14) and picking ϕ such that we obtain showing equivalence of the norms. Although (115) is a reliable estimator, it is not a computable quantity. Therefore, we introduce a new quantity η R , which is based only on computable values: Here, Theorem 3 (Reliable error estimation). Assume (A 1 ) and (A 5 ) For every t ∈ S, the error norm ||e π || can be approximated by the error estimator η R . The following bounds hold: Proof. For any ϕ ∈ V H with ||φ|| L 2 (Ω) = 1, the following estimate holds: This provides us with the first inequality in (120). The second inequality requires some auxiliary definitions. Let f H denote the Galerkin projection of f . We useR B andR E to denote the residuals where f is replaced with f H . Additionally, we introduce to conveniently manipulate the residual norm. Finally, let C f be a constant defined as Here, C is a constant independent of π, ρ, h and H.
Next, we bound each of the terms of η R to obtain a lower bound for the error.
Dividing (124) once by its common factor and rearranging terms results in which, after applying the triangle inequality results in To provide an upper bound on the second term of (120), we use the equivalence of the error and residual norms: Dividing both sides of (127) by its common factor results in: Combining (127) with (128) results in which yields the following lower bound on the error: We remark that the presence of H and h 2 error terms is a typical feature of two-scale models.

Macroscopic mesh refinement strategy
The strong separation of space scales in our setting (microscopic vs. macroscopic) allows us to propose a macroscopic mesh refinement strategy only weakly biased by the distribution of microscopic errors, based on local error indicator η R,B and global error indicator η R . Inspired by the popular and intuitive approach presented in e.g. [33], the goal of this refinement is to reduce the global approximation error and keep the error locally below a prescribed tolerance, i.e. select a refinement strategy satisfying where we denote the desired tolerance withη.
Our refinement strategy relies on the double sided estimate (120) stated in Theorem 3. This inequality gives us a satisfied upper bound on the global error. We solve (18) on B H for some t, compute the error estimator, and evaluate if refinement is necessary. If so, we repeat this process until the error estimator has reduced to a satisfactory level.
The set of triangles to be refined on each iteration follows directly: as well as boundary triangles Each element B ∈ Q B is partitioned into 2 d1 new elements (with d 1 the dimension of Ω), while each element B ′ ∈ Q ′ B is refined 2 d1−1 to ensure no vertices collide with edges. An illustration of this process in two dimensions is given in Figure 2 and   The convergence estimates from Theorem 3 ensure that this procedure will indeed halt for any fixedη.

Conclusion
We constructed a semidiscrete Galerkin approximation of our elliptic-parabolic two scale system (P 1 ) and showed that this approximation is well-posed and that the obtained sequence of Galerkin approximants converges in suitable spaces to the weak solution to the continuous system. Furthermore, we derived a priori rates of convergence and proposed an a posteriori grid refinement strategy at the macroscopic scale.
As natural next steps, future work will address the fully discrete two-scale Galerkin approximation as well as the numerical implementation of the method so that the proven convergence rates can be confirmed and the macroscopic refinement strategy can be tested. Additionally, a mesh refinement strategy on the microscopic level could also be considered for either this problem setting, or for its elliptic-elliptic variant (obtained by letting t → ∞ in (P 1 )).
At this stage, we would like to remark that the interaction between H B and h 2 in the error structure gives an indication on how to choose the mesh size h based on the error estimators in H B . Without going into details, it is worth mentioning that in principle, one can choose the microscopic mesh size to correspond to the macroscopic grid. This way one can ensure that the macroscopic and microscopic errors are roughly of the same order.