Analysis of a Coupled Darcy Multiple Scale Flow Model Under Geometric Perturbations of the Interface

For a Darcy flow coupled system modeling fluid exchange between two regions with fluid resistance at different scale, we address the question of continuity of the solution with respect to geometric perturbations of the interface. The original interface and its perturbation define three flow regions exchanging fluid, these are modeled as a coupled system in mixed variational formulation. For numerical purposes a-posteriori perturbation/error estimates are given.


Introduction
Layered media is a common structure in the analysis of fluid flow through geological systems, as in oil extraction or subsurface water preferential flow.Modeling this phenomenon presents naturally a multiple scale flow problem with regions of slow velocity of order O(1) in a rock matrix and fast flow of order O(1/ ) in the fissures.Here, the small number > 0 is introduced to account the order of scale at which the physical phenomenon is taking place in agreement with field data.A common approach consists in introducing a lower dimensional interface manifold Γ between regions, model each region with a suitable law and boundary conditions and finally state fluid transmission conditions between them across the interface, namely: normal stress and normal flux balance.Therefore, the placement and description of such interface becomes a delicate issue since it couples two physical phenomena occurring at different scale, hence a perturbation of the interface geometry can introduce a substantial error which needs to be addressed see Figure 2. On one hand such perturbation comes from the natural limitations when collecting data field of geological strata, on the other hand the numerical discretization of the models can not describe exactly a curved interface which is not a piecewise polynomial surface.
In the present work the stationary problem is modeled by a coupled system of partial differential equations with Darcy flow in both regions in a well-posed mixed variational formulation.The analysis will be centered at the weak variational formulation of the equations without stepping in critical point theory.Using coupled partial differential equations systems in analyzing this phenomenon has been largely studied with notorious success.From the analytic point of view see Arbogast and Lehr (2006) for a Darcy-Stokes coupled system model and Morales andShowalter (2010, 2012).for a Darcy-Darcy coupled model, see Chen, Gunzburger, and Wang (2009), Cao, Gunzburger, Hua, and Wang (2010) for discussion of the interface conditions, see Babuska and Gatica (2010), Gatica, Meddahi, and Oyarzúa (2009), Martín, Scheid, and Smaranda (2012) for the study of finite element discretization of the systems and Arbogast and Brunson (2007), Chen, Gunzburger, Hua, andWang (2011), Martin, Jaffré, andRoberts (2005) for modeling these phenomena from the numerical point of view.
A problem such as this has two sources of singularity: the geometric one coming from the width of the cracks which is much smaller than the porous medium and the physical one coming from the resistance to the flow much lower in the fractures than in the pores.However, most of the literature of theoretical approach in the field is focused on appropriate scaling the governing laws as in Lévy (1983), removing the singularity of the system by techniques of homogenization see Hornung (1997, Ed.) or studying the interface conditions as in Saffman (1971) which provide well-posed problems; in particular the mixed variational formulations and their associated finite element schemes have become a vigorous research area see Layton, Scheiweck, and Yotov (2003).However, the continuity of the solution with respect to the geometry of the interface is yet to be studied.
We introduce now equations, boundary and interface conditions modeling the problem in general, the unknowns are the fluid velocity u and pressure p.The constitutive law for Darcy's flow together with the mass conservation law Here Ω is the domain of interest, g(x) denotes the gravity force, F(x) the sources and R(x) is the flow resistance i.e. the fluid viscosity times the inverse of the permeability of the medium, to be scaled consistently with the fast and slow flow regions of the medium.Drained and non-flux boundary conditions on different parts of the domain boundary will be specified to set a boundary value problem.The fluid exchange across the interface separating the regions, see Figure 1 are given by p 1 − p 2 = α u 1 and (1c) The coefficient α indicates the fluid entry resistance of the rock matrix.

Figure 1. Curved interface perturbation
Next we describe the geometry of the problem.The original and perturbed interfaces define three regions: two disconnected O 1 , O 3 where the permeability of the original an perturbed systems agree and a third region O 2 trapped between the interfaces, where the permeability of the original and perturbed systems disagrees, as Figure 1 depicts.The identification of these three regions suggests formulating the problem with a triple coupled system in mixed variational formulation in order to introduce the necessary degrees of freedom on the test functions and to make both solutions, the original and the perturbed one, comparable i.e. belonging to the same space.We analyze the particular geometric setting in which the interface Γ is a horizontal flat surface and the perturbation Γ ζ is a piecewise C 1 surface always above the original interface as in Figure 2. It is clear that such case is not useful in practice since flat surfaces do not induce curved perturbations in a natural way, nevertheless it gives significant simplicity to the notation and calculations without being conceptually far from the practical scenario.It is important to stress that the analysis should include continuous perturbations Γ ζ which are not differentiable everywhere as the red line in Figures 1 and 2, because this is the natural type of perturbation introduced by numerical discretization.
We state sufficient conditions on the forcing terms as well as the geometric perturbation in order to conclude continuity statements with respect to the interface perturbation for fixed data.The structure of the continuity estimates hints the necessity of introducing a-posteriori perturbation/error estimates for practical purposes i.e. in terms of the approximate solution obtained from numerical computation; these will be presented too.Vectors are denoted by boldface letters, as are vector-valued functions and corresponding function spaces.We use x to indicate a vector in 4 N−1 ; if x ∈ 4 N then the 4 N−1 × {0} projection is identified with x def = (x 1 , x 2 , . . ., x N−1 ) so that x = ( x, x N ).The symbol ∇ is the gradient in the first N − 1 derivatives.Given a function f : 4 N → 4 then A f , M f dS are the notations for its volume integral in the set A ⊆ 4 N and its surface integral on the N − 1 dimensional submanifold M ⊆ 4 N respectively.For any given set A we write A for its indicator function.λ N , λ N−1 will denote the Lebesgue measure in 4 N and 4 N−1 respectively.We use êN for the unitary vector in the N-th direction and ν for the outwards normal vector.Finally, n denotes the normal upwards vector i.e. n • êN ≥ 0.

Geometric Setting
We are to analyze the behavior of the solution under small perturbations of the interface, the essentials of this phenomenon can be captured in a very simple geometric setting.Let Γ be a connected set in 4 N−1 × {0} whose projection is relatively open in the 4 N−1 trace topology; in the following we make no distinction between these two domains.Let Ω 1 , Ω 2 be smooth open regions in 4 N separated by Γ, i.e. ∂Ω 1 ∩ ∂Ω 2 = Γ.
In order to represent perturbations of the interface Γ we introduce the following definition Definition 2.1 We say the set T (Γ, G) of piecewise C 1 perturbations of the interface Γ contained in Ω is given by Where Γ indicates the closure of Γ.The interface associated to ζ ∈ T (Γ, Ω) is given by the set The domains associated to ζ ∈ T (Γ, Ω) are defined by the sets

The Strong Problem and Its Interface Perturbation
Consider the solution of the strong problem (5) presented in Morales and Showalter (2012) and in a less general version in Morales and Showalter (2010).Such solution is to be compared with the solution of the system (6) corresponding to a geometric perturbation of the interface Here u, v are the velocities and p, q the pressures, g is the gravity force and F the fluid sources of the medium.It is assumed that the forcing terms of normal stress h and normal flux f are defined on Γ as well as on Γ ζ ; for a full exposition of their associated interface balance statements see Morales and Showalter (2012).Drained (5c), (6c) and null normal flux (5h), (6h) boundary conditions on the indicated parts of the boundary are a natural choice.
In both problems the scaling of the resistance coefficient R and R in different parts of the domain will impose the flow field velocity been slow and fast flow in the corresponding regions.Due to the physics of the problem it can be assumed that R 0 ≤ R(x) ≤ R 1 for positive constants R i , i = 0, 1.For notation clarity, the resistance coefficient will be omitted in the following.
Definition 2.2 Define the open set where the perturbation Γ ζ is above the original interface Γ and the open set where the perturbation Γ ζ is below the original interface Γ A complete analysis of the question would demand both of the above sets not to be trivial.However, in order to ease calculations we assume ζ( x) ≥ 0 for all x ∈ Γ.
Furthermore, due to these observations the following notation is introduced Definition 2.3 Let ζ ∈ T (Γ, Ω) a non-negative function, then define the following associated domains Remark 2.2 Observe the following identities Remark 2.3 For the treatment of the general case, when the perturbation ζ is allowed to change sign see part (vi) in section 5.
Finally, in order to make comparable the solutions (u, p), (v, q) of problems ( 5) and ( 6) we force them to be in the same space introducing an artificial interface in each system.Problems ( 5) and ( 6) are reformulated in the following way Remark 2.4 (i) Notice that the modeling technique used in Morales and Showalter (2012) shows that in order to have the well-posedness of the problems ( 10) and ( 6) the presence of the three regions is not necessary.This question can be addressed using only two regions as in problems ( 5) and ( 6) respectively.
(ii) The Equations ( 10f), (10g) as well as the statements of balance across the interface Γ ζ (10h) and (10i) ensure that the solutions of problems ( 6) and (10) coincide.Similarly, the Equations (11f), (11g) and the statements of balance across the interface Γ (11d) and (11e) are introduced to make the solutions of problems ( 6) and (11) agree.
(iii) Overall the introduction of artificial interfaces permits the treatment of the problems in a common geometric setting, consequently the solutions to both problems can be modeled in the same space of functions.Without the presence of artificial interfaces, common spaces could not be established and consequently it is impossible to compare the solutions.This will be seen more clearly in section 3.

The Weak Problem and Its Perturbation in Three-Region Mixed Formulation
Consider the spaces endowed with their corresponding norms With the spaces described above the problem (10) has the following mixed variational formulation.
Similarly, the mixed variational formulation of the perturbed problem ( 11) is given by.
The variational statements ( 13) and ( 14) are problems in mixed variational formulation, see section 8.6 in Atkinson and Han (2005) or chapter II in Brezzi and Fortin (1991).
Remark 2.5 (i) It is paramount to observe that the spaces defined in Equations ( 12) are fully decoupled, i.e. they do not require matching conditions of fluid exchange on the test functions.The properties of normal balance stress (10d), ( 10h), ( 11d), (11h) and mass conservation balance (10e), ( 10i), ( 11e), (11i) across the interfaces Γ and Γ ζ occur only on the solutions (u, p), (v, q) ∈ V × Q.Such freedom of the quantifiers w ∈ V and r ∈ Q allows to test the variational statements of problems ( 13) and ( 14) with great flexibility; this property is crucial in the following section and consequently for the whole technique.
(ii) Also notice that the condition ζ ∈ W 1,∞ (Γ) and ζ ∈ C(Γ) implies that the domains O i for i = 1, 2, 3 are smooth and therefore the integration by parts applies on the function spaces.Without these conditions the weak variational formulations ( 13) and ( 14) would not be possible.

Perturbation Estimates and Geometric Aspects
In the present section we obtain estimates in norm for the difference of solutions, original and perturbed and derive sufficient conditions on the forcing terms as well as the geometric perturbations ζ which allow to conclude continuity and convenient estimates.Test (13a) with (u ( On the other hand, testing (13b) with (p The sum gives of expressions ( 15) and ( 16) gives (17) We repeat the same procedure for the perturbed problem (14).Test (14a) with (u (18) Finally, subtracting the added systems ( 17) and ( 18) we have The association the first and second summands in the second line of (19) simplifies the integral O 2 ∇q • (u − v).Recalling (10f) and reordering we have (20) Next the forcing terms of the interfaces need to be analyzed.

Regularity Conditions for the Interfaces Forcing Terms and Related Estimates
In this subsection we recall ν denotes the unitary vector outwards normal vector to the boundary of a given open smooth region G ⊆ 4 N .The classical duality relation below Tartar (2007) will be frequently used where w ∈ H div (G) and r ∈ H 1 (G) are arbitrary .

The Normal Flux Interface Forcing Term f
In order to assure continuity estimates for the normal flux forcing term f across the interfaces certain type of regularity is required.From now on we assume there exists G ⊆ Ω a smooth open region such that Γ ∪ Γ ζ ⊆ G and that there exists a function ), e.g.consider Φ a continuous differentiable function from G to 4 N .Hence, the first two summands of the second line in (20) transform as The equality comes from the normal trace duality relation ( 21) for elements Φ ∈ H div (O 2 ) and the last inequality is the mere application of Cauchy-Schwartz.The L 2 (O 2 )-norm of the pressures difference in the expression above needs to be estimated, due to the drained boundary conditions ( 10c), (11c) we know Here d > 0 denotes the diameter of Ω 1 ∪ Ω 2 ; the above holds since the Poincaré constant is always less or equal than the diameter of the domain divided over √ 2. Next, recall the Darcy-type relations (10a), (10f), ( 10a) and (10f), therefore the expression can be estimated by Where we used the fact that x 2 ≤ x 1 for all x ∈ 4 2 .

The Normal Stress Interface Forcing Term h
From now on assume there exists G ⊆ Ω a smooth open region such that Γ ∪ Γ ζ ⊆ G and that h ∈ H 1 (G); without loss of generality it can be assumed the set G is the same for both forcing terms.Recalling (10e), (10i), (11e) and (11i) the following identity holds Using the normal trace duality relation ( 21) for elements w ∈ H div (O 2 ) and Cauchy-Schwartz inequality we have Equations ( 10g) and (11g) yield Combining ( 22) and ( 23) with (20) yields

Expression (21) gives the inequality
. Finally, we reorder the expression as Notice the first summand on the right hand side of the inequality above is independent from the solution.In the final section the last two summands are estimated.

Flux Perturbation Estimates
In this section two possible estimates for the last two summands of the second line in (24) are given; the first one aiming to conclude continuous dependence with respect to the geometric perturbation and the second one seeking a-posteriori estimates.

Continuity Estimates
In order to attain continuity estimates we modify the last two summands of (24) in the following way Together with (24) yields

A-Posteriori Estimates
The last two summands of the second line in inequality (24) are estimated in terms of the approximate solution v. Vol. 5, No. 4;2013 We manipulate the summands and apply Cauchy-Schwartz to get Associating the above estimate with (24), reordering and rearranging the terms on both sides of the inequality furnishes (26)

Continuity and a-Posteriori Analysis of the Perturbation
This section begins recalling a well-known result Lemma 4.1 Let {ζ n : n ∈ N} be a sequence in C 1 (G) and A, B, x: for all n ∈ N then (i) The explicit estimate holds Proof.(i) The inequality (28) follows from completing squares on the inequality ( 27); (ii) is an immediate conclusion of (i).It is also important to stress that, by definition, the involved functions A, B and x always assume non-negative values.

Continuity Analysis of the Solution
Now we prove the continuity of the flux n) , q (n) ) denotes the solution of ( 14) for the perturbation 25), applying Cauchy-Schwartz on the last two summands of the second line in right hand side yields Where O n 1 , O n 2 and O n 3 are the domains defined by the perturbation ζ n .From the expression above we define the functions Finally, in order to get a-posteriori estimates for the pressure we repeat the technique of Theorem 4.3 and get Recalling inequality (35) and combining it with (37) we have

Discussion and Concluding Remarks
The present work yields several conclusions.We start listing those related to the variational formulation of the problem (i) The geometric perturbation of the interface introduces regions of disagreement in the coefficients of resistance.In order to get the estimates it is necessary to test the variational statements with convenient functions which are not continuous across the interface as in section 3.Even if the interface forcing terms f , h were null this necessity remains.
(ii) A direct variational formulation of the problem does not allow testing with functions which are discontinuous across the interface since the functional setting for test functions is H 1 .
(iii) The mixed variational formulations L 2 − H 1 and H div − L 2 both demand coupling conditions of the spaces on the interfaces Γ, Γ ζ .In the first case the space of pressures must be continuous and in the second the normal traces of the flux must be continuous across the interface.Therefore, in both cases the possibility of testing with discontinuous functions across the interfaces is not allowed.Finally, the first formulation requires h = 0 and the second requires f = 0 i.e. in both cases some generality with respect to the formulation used in this work is lost.
(iv) The mixed-mixed variational formulation of section (2.3) introduced in Morales and Showalter (2012) is the unique variational setting introducing the necessary degrees of freedom in the function spaces i.e. permitting jump-discontinuities across the interfaces on the test functions.(vi) For the general type of perturbation, when ζ can take positive and negative values the strategy is identical to the one presented here.There will be need to consider both subdomains O ζ , U ζ given in Definition 2.2.The variational statements corresponding to ( 13) and ( 14) would involve more terms and in particular, the non-symmetric interface terms on ∂O ζ and ∂U ζ would have opposite sign.The function spaces V, Q given in Equation ( 12) are still adequate to address the question.Aside from the aforementioned observations and larger expressions the strategy and reasoning follows exactly the one presented here.
Finally we list the conclusions with regard to the non-linearity of the problem as well as the estimates order (i) It is important to observe in both theorems of continuity above (4.2),(4.3) the role played by the term u L 2 (O n 2 ) .Well-posedness of the problem (13) furnishes global estimates with respect to the forcing terms i.e. it can be claimed u L 2 (O n 2 ) ≤ C( F L 2 (Ω) + g L 2 (Ω) + h L 2 (Γ) + f L 2 (Γ) ) and certain constant C > 0, nevertheless we can not assure a local estimate such as u L 2 (O n 2 ) ≤ C( F L 2 (O n 2 ) + g L 2 (O n 2 ) + h L 2 (Γ) + f L 2 (Γ) ).Therefore the term u L 2 (O n 2 ) can not be excluded in the continuity statements.(ii) The problem of continuous dependence with respect to the geometric interface perturbation is nonlinear as it is evident by the nature of the continuity and a-posteriori estimates in Equations ( 25), (26) respectively.
(iii) The introduction of the mixed variational formulation in Equations ( 13), ( 14) solves the issue of a functional context consistent with the necessity of jump-discontinuities across the interface.However it does not linearize the problem itself as the implicit nature of the attained estimates show.
(iv) The nature of the continuity estimates suggests numerical experimentation as the most feasible approach to investigate the convergence rate problem.
(v) Observe that the flux continuity Theorem 4.3 involves building functions A(•), B(•), both of them are multiplied by 1/ .Though this amplification factor is fixed, it reveals that a geometric perturbation of the interface introduces an error of order O(1/ ) which is rather significant.The same observations hold for the a-posteriori estimates.(38) since the functions a(•), b(•) are both affected by the same amplification factor.
(vi) For geometric perturbations independent from time, the evolution problem can be analyzed by the same technique presented in Morales and Showalter (2012) using analytic semigroups.Since the perturbation occurs in space and not in time the essentials are captured by the stationary problem we have presented here.
(vii) Finally, the work presented here is a first step towards the understanding of flow in deformable porous media i.e. when the geometry changes with time.This is future work, however it is the author's opinion that conditions on the time deformation ratio of the perturbations will play a fundamental role in the analysis; as crucial as the introduction of an artificial interface has been here.

Figure 2 .
Figure 2. Flat interface perturbation (v) Stability and continuity statements demand extra hypothesis of regularity on the forcing terms discussed in sections (3.1.1)and (3.1.2) which are reasonable for sources in fluid flow problems.