# Well-posedness, discretization and preconditioners for a class of models for mixed-dimensional problems with high dimensional gap

In this work, we illustrate the underlying mathematical structure of mixed-dimensional models arising from the composition of graphs and continuous domains. Such models are becoming popular in applications, in particular, to model the human vasculature. We first discuss the model equations in the strong form which describes the conservation of mass and Darcy's law in the continuum and network as well as the coupling between them. By introducing proper scaling, we propose a weak form that avoids degeneracy. Well-posedness of the weak form is shown through standard Babuška-Brezzi theory. We also develop the mixed formulation finite-element method and prove its well-posedness. A mass-lumping technique is introduced to derive the two-point flux approximation type discretization as well, due to its importance in applications. Based on the Babuška-Brezzi theory, error estimates can be obtained for both the finite-element scheme and the TPFA scheme. We also discuss efficient linear solvers for discrete problems. Finally, we present some numerical examples to verify the theoretical results and demonstrate the robustness of our proposed discretization schemes.

## Authors

• 1 publication
• 16 publications
• 8 publications
09/05/2020

### Fully discrete finite element approximation for a family of degenerate parabolic mixed equations

The aim of this work is to show an abstract framework to analyze the num...
04/27/2021

### A mass-, kinetic energy- and helicity-conserving mimetic dual-field discretization for three-dimensional incompressible Navier-Stokes equations, part I: Periodic domains

We introduce a mimetic dual-field discretization which conserves mass, k...
09/18/2019

### A macroelement stabilization for multiphase poromechanics

Strong coupling between geomechanical deformation and multiphase fluid f...
09/29/2021

### A Mass Conserving Mixed hp-FEM Scheme for Stokes Flow. Part III: Implementation and Preconditioning

This is the third part in a series on a mass conserving, high order, mix...
12/20/2017

### Implementation of mixed-dimensional models for flow in fractured porous media

Models that involve coupled dynamics in a mixed-dimensional geometry are...
06/27/2021

### Radial complex scaling for anisotropic scalar resonance problems

The complex scaling/perfectly matched layer method is a widely spread te...
07/05/2019

### A note on optimal H^1-error estimates for Crank-Nicolson approximations to the nonlinear Schrödinger equation

In this paper we consider a mass- and energy--conserving Crank-Nicolson ...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

Coupled fluid flow in networks and porous domains arise in various applications, including blood flow in the human body as well as wells in geological applications. Such models are referred to as mixed-dimensional when the network flow is simplified to a family of 1D domains along with the network edges111Some authors refer to this problem as multiscale (see e.g. [DangeloQuarteroni2008, KopplVidottoWohlmuthZunino2018]), however, we prefer the nomenclature mixed-dimensional to avoid confusion with equidimensional multiscale methods such as are encountered in (numerical) homogenization problems.. Moreover, when the coupling between the network and the domain exceeds two topological dimensions, the model is referred to as having a high dimensional gap [Nordbotten2019, LaurinoZunino2019]. A high dimensional gap thus arises when the flow in the network is connected to a domain of dimension through its leaf nodes, or when the flow in the network is connected to a domain of dimension through its edges.

In this paper, we consider the problem composed of flow in one or more trees, coupled with a (porous) domain. This setting is motivated by blood flow in the brain, wherein the networks are the arterial and venous trees, and the domain is the sub-resolution capillary bed. Recognizing that the leaf nodes in the tree (referred to as ”terminals” hereafter) are in applications an artifact of limited imaging resolution, we consider in our equations a mesoscale model wherein fluid is distributed into the porous domain in a support region near the terminals. Such models have recently been introduced in [HodnelandHansonSaevareidNaevdalLundervoldSolteszovaMunthe-KaasDeistungReichenbachNordbotten2019] and also considered in [koch2020modeling], and are attractive also from a mathematical perspective, as they avoid the singularities which otherwise characterize the coupled equations. In this work, we will not adopt the precise models used in [HodnelandHansonSaevareidNaevdalLundervoldSolteszovaMunthe-KaasDeistungReichenbachNordbotten2019, koch2020modeling] directly, as they consider an explicitly given structure of fluid distribution between the network and the porous domain. In contrast, we will use a more canonical formulation, where the flow resistance is given, and the fluid distribution from the terminal is calculated.

Previous mathematical analysis of models with high dimensional gap has to a large extent been focused on how to handle the singularities arising when the coupling is ”point-wise” between the network and the domain (see e.g. [DangeloQuarteroni2008, KopplVidottoWohlmuth2016, GjerdeKumarNordbottenWohlmuth2019]). In contrast, the model discussed herein has to our knowledge not been subjected to mathematical analysis before. In the absence of singularities, we exploit in this paper the framework recently developed for problems with small dimensional gap [BoonNordbottenVatne2019], and define mixed-dimensional variables and operators for the coupled problem. Together with appropriately defined integration and inner products, we then observe that we have available tools such as a mixed-dimensional Stokes’ theorem, integration by parts, and Hilbert spaces. This forms the building blocks for our well-posedness results and numerical analysis.

The main results of the paper are thus as follows:

• A general, non-singular model for a class of problems with a large dimensional gap.

• Well-posedness theory for both the continuous and finite-dimensional problem.

• Convergence results for mixed finite-element approximation and a finite volume variant.

• Efficient linear solvers for different discretizations.

• Numerical validation and application to a high-resolution data-set of a real human brain.

We structure the paper as follows. In Section 2 we present the model equations in both strong and weak forms and show well-posedness. In Section 3 and 4 we state and analyze the finite-element and finite volume approximations, respectively. Efficient solvers for the coupled problem are proposed in Section 5. The theoretical results are validated in Section 6. Finally, we give some conclusions in Section 7.

## 2 Model Equations

In this section, we discuss the basic model equation for modeling the coupled network-Darcy flow in brain. We will first discuss its strong form and then derive the weak form by introducing proper spaces.

### 2.1 Strong From

We are concerned with a domain  (which models the capillaries). In addition, we are concerned with a finite collection of rooted trees  with node (vertex) set  and edge set  (which model resolved arteries and veins). Those trees are disjoint and, therefore, form a forest with node set and edge set . The node set can be subdivided into three disjoint subsets, the first and last of which are assumed to be non-empty: root nodes , interior nodes , and terminal nodes . Note that and we use , , and to denote the root nodes, interior nodes, and the terminal nodes of a given tree , respectively. Naturally, we also have . We further divided the root nodes into two disjoint sets , which consists of the Dirichlet root nodes, and , which consists of the Neumann root nodes. The Dirichlet root nodes will be treat explicitly as Dirichlet boundary conditions and the Neumann root nodes will be implicitly handled through the right-hand side of the conservation laws on the graph. Following the same convention, and denotes the Dirichlet or Neumann root nodes of a given tree , respectively. Note that each tree can only have one root. Therefore, we can subdivide the forest into two disjoint sub-forests, i.e., Dirichlet rooted forest , which contains all the Dirichlet rooted trees , and Neumann rooted forest , which contains all the Neumann rooted trees . Naturally, and . Furthermore, we define, , , , , , and . We denote the set of the neighbors of the node as and the set of all the edges meeting at  as . Note that, and , if . These concepts are illustrated for in Figure 1.

As primary variables we choose the domain pressure potential  and the node pressure potentials . Furthermore, we consider the fluid mass fluxes denoted in the domain as , fluid mass flow from node  to   denoted  and fluid mass flow transferring from terminal node  to point  denoted . This last variable models the flow in unresolved arteries and veins, and is a novel component our this work.

First, we consider the model equations for mass conservation and they are given as follows based on the above definitions and notation.

 (Conservation of mass in brain tissue)∇⋅qD−∑i∈NTqTi=rD,in Ω (1) (Conservation of mass at interior nodes)−∑j∈NiqNj,i=rNi,for all i∈NI∪NN (2) (Conservation of mass at terminal nodes)∫ΩqTi(x)dx−qNNi,i=rNifor% all i∈NT (3)

Here, the signs in (1)-(3) are chosen such that the right-hand-side terms represent sources added to the system. Moreover, although both and are used in (2) for notational convenience, they should be understood as one unknown with a sign difference.

Next we verify the global conservation of mass based on (1)-(3) as follows,

 (Stokes' Theorem)∫∂ΩqD⋅ndx =∫Ω∇⋅qDdx (By (1)) =∑i∈NT∫ΩqTi(x)dx+∫ΩrDdx (By (3)) =∑i∈NTqNNi,i+∑i∈NTrNi+∫ΩrDdx (By (2)) =∑i∈NDqNi,Ni+∑i∈NN∪NI∪NTrNi+∫ΩrDdx.

Where the last step is also known as the Graph-Stokes’ Theorem, which is the counterpart of the Stokes’ Theorem on graphs.

In order to present the constitutive laws properly, we need the introduce some material parameters, all of which are assumed to be non-negative (precise bounds are given later). For each edge , we assign a conductivity , which can be considered as the edge weights in certain sense. In the domain, for each

we assign a permeability tensor

. For each terminal node , we assign connectivity function . Now based on the assumption that the potential flow is linear, we have the following constitutive laws.

 (Potential flow in brain (Darcy))qD=−kD(x)∇pD,in Ω, (4) (Potential flow in network (Poiseuille))qNi,j=−kNe(i,j)(pNj−pNi),for e(i,j)∈E, (5) (Potential flow from network to brain)qTi(x)=−kTi(x)(pD−pNi),for i∈NT. (6)

The coefficient functions , , represent redistribution in a small region around the terminal node , thus, can be assumed to have compact support in some domain .

###### Remark 1

In practice the characteristic length scale of is comparable to the distance to the nearest neighbor, i.e.

 diam(Bi)=O(minj∈NT|xi−xj|). (7)

Moreover, the grid is frequently given by the voxel resolution of the image and the terminals are due to a finite resolution effect, and thus

 minj∈NT|xi−xj|=O(h), (8)

where is the mesh size. The constants hidden in the notation in (7) and (8) are usually between to in practical applications that we are interested in. Consequentially, also is compactly supported in . While these considerations could be applied to further refine some of the constants in the proofs below, we will not exploit these details in this paper.

In addition to the conservation laws and constitutive laws, we also need boundary conditions to close the system. For the sake of simplicity, we only consider the case of homogeneous Neumann data on  and the Dirichlet root nodes , i.e.,

 qD⋅n=0on ∂ΩandpNi=0,i∈ND. (9)

We want to point out that our results and analysis below also hold for other types of boundary conditions as well only at the cost of extra notation. The choice of Neumann data on is in a sense the most difficult case, as the inf-sup proofs can be simplified considerably in the case where there is a measurable subset of the boundary with Dirichlet data.

We close this subsection by the observation that by definition . Therefore, although the total number of is , we only use half of them as the unknowns, i.e., one unknown, or , for each edge . The choice is arbitrary. In this work, we choose the one follows the direction from the root node to the terminal nodes. This direction is also the assigned orientation of the corresponding edge (i.e., if we choose , which means the fluid mass flows from node to node , the edge is oriented such that it starts at node and ends at node ). This allows us to define the following signed incidence matrix , such that

 Gℓ,i=⎧⎨⎩1,if flow on edge ℓ starts at % node i−1,if flow on edge ℓ ends at node i0,otherwise (10)

We want to point out that the signed incidence matrix represents a discrete gradient on the graph and its transpose serves as a discrete divergence.

### 2.2 Mixed-dimensional formulation and scaling

The model equations given above contain essentially three expressions of fluxes ( and ), and two expressions of potentials ( and ). It will simplify the following exposition and analysis considerably to treat these as mixed-dimensional variables, on which we define mixed-dimensional operators.

Therefore, let the mixed-dimensional flux be denoted , and defined as the triplet of fluxes . Equivalently, the mixed-dimensional pressure is denoted and defined as . Now, we define the mixed-dimensional divergence operator as follows,

 D⋅q=D⋅[qD,qT,qN]:=[uD,uN], (11)

where

 uD:=∇⋅qD−∑i∈NTqTiand% uNi={−∑j∈NiqNj,i,i∈ NI∪NN,∫BiqTi(x)dx−qNNi,i,i∈NT. (12)

Similarly, we define the mixed-dimensional gradient as

 Dp=D[pD,pN]:=[vD,vT,vN], (13)

where

 vD:=∇pDandvTi(x):=pD(x)−pNi, i∈NT,andvN=GpN. (14)

In addition, we introduce the function which contains all the material functions , , , and , , in (4) to (6), such that

 K−1[qD,qT,qN]:=[(kD)−1qD,(kT)−1qT,(kN)−1qN],

where . It is now straight forward to verify that with these definitions, the conservation laws (1)-(3) can be summarized as

 D⋅q=r, (15)

where . Furthermore, the constitutive laws (4) to (6) can be summarized as

 q=−KDp. (16)

While the physical model formulation is satisfactory for non-degenerate , , it will be beneficial to rescale the coupling flux to avoid considering a degenerate mass matrix when for some points or region in . To that aim, we introduce the square-root of the transfer coefficient , , and the scaled transfer mass flux , , is defined as . Thus, we replace (1) (3), and (6) with

 (Conservation of mass in brain tissue)∇⋅qD−∑i∈NTkSiqSi=rD,in Ω, (17) (Conservation of mass at terminal nodes)∫BikSiqSi(x)dx−qNNi,i=rNifor all i∈NT, (18) (Potential flow from network to brain)qSi(x)=−kSi(x)(pD−pNi),for i∈NT, (19)

respectively. We note that a similar scaling has been applied previously to handle degeneracies occurring in mantle dynamics [arbogast2017mixed] and flows in fractured porous media [boon2018robust].

Equivalently, We denote the scaled mixed-dimensional flux as , and the scaling such that

 S−1[qD,qT,qN]:=[qD,(kS)−1qT,qN].

Thus, , and we can introduce the rescaled divergence and gradients as and , respectively. The rescaled conservation equations are then summarized as

 DS⋅qS=r. (20)

The rescaled conservation equations are summarized as

 qS=−KSDSp, (21)

where .

In this setting, we allow for degeneracy of the coupling term in the sense that we allow . However, we require that is bounded from above, i.e., for and . Furthermore, for all , we require it to hold that

 ∫BikSidx=ckSi≥ckS>0,

where is a generic constant.

### 2.3 Weak Form

In this subsection, we derive the weak formulation of the system. The development will be equally valid for both the original model, equations (15) and (16), as well as the re-scaled model, equations (20) and (21). Thus we will omit the superscript on the mixed-dimensional operators and variables to reduce notational overload. Nevertheless, in order to allow for degeneracies, we will always have the rescaled equations in mind, and thus when we need to specifically refer to , and consider the coefficient to appear in the differential operator as opposed to the material law.

We first introduce proper function spaces. We begin by defining a mixed-dimensional square-integrable space for pressure as follows,

 L2:=L2(Ω)×l2(N∖ND),

where is the standard space defined on domain and is the standard space defined on the node set . For flux, we consider a space with bounded mixed-dimensional divergence as follows,

 H(div):=H(div,Ω)×∏i∈NTL2(Bi)×l2(E)

where is the space defined on such that the functions and their divergence are both square-integrable. In addition, are standard space defined on , , and is the standard space defined on the edge set .

We associate the mixed-dimensional space with the following inner product,

 (p,w)=([pD,pN],[wD,wN]):=∫ΩpDwDdx+∑i∈N∖NDpNiwNi,∀p,w∈L2.

Similarly, we introduce the following inner product on ,

It is important to note that the inner products are defined such that integration-by-parts holds in for the mixed-dimensional operators (both original and re-scaled cases).

###### Lemma 1 (Integration by parts)

For any and , we have

 (Dp,q)+(p,D⋅q)=∫∂ΩpDqD⋅ndx+∑i∈NDpNiqNNi,i (22)

###### Proof

By a direct calculation (using the re-scaled operators and variables, the derivation for the original case is the same), we have that

 (p,D⋅q)= ∫ΩpD⎛⎝∇⋅qD−∑i∈NTkSiqSi⎞⎠dx−∑i∈NI∪NNpNi∑j∈NiqNj,i +∑i∈NTpNi(∫BikSiqSidx−qNNi,i) = −∫Ω∇pD⋅qDdx+∫∂ΩpDqD⋅ndx−∑e(i,j)∈E(pNi−pNj)qNi,j +∑i∈NDpNiqNNi,i−∑i∈NT∫BikSi(pD−pNi)qSidx = ∫∂ΩpDqD⋅ndx+∑i∈NDpNiqNNi,i−(Dp,q),

which completes the proof.

To derive the weak formulation, we need to incorporate the boundary conditions. Recall that we consider on , therefore, we define the following functions space with boundary conditions,

 H0(div):=H0(div,Ω)×∏i∈NTL2(Bi)×l2(E)⊂H(div),

where . In addition, with the material function , we introduce a weighted inner product on as follows,

 (q,v)K−1:=(K−1q,v).

Using the above function spaces and notation, together with the mixed-dimensional integration by parts formula (22) and the homogeneous Dirichlet boundary conditions (9) on , i.e., , , we have the following weak form for the conservation laws (17), (2), (18) and constitutive laws (4), (5), (19): Find  and , such that

 (q,v)K−1−(p,D⋅v)=0,∀ v∈H0(div), (23) −(D⋅q,w)=−(r,w),∀ w∈L2. (24)

Note that due to the integration by parts formula, if non-homogeneous boundary data is considered, this would appear as extra right-hand side terms in equation (23).

### 2.4 Well-posedness

In this subsection, we focus on the well-posedness of the weak formulation (23)-(24). We first introduce the following norm on ,

 ∥p∥2L2:=(p,p). (25)

And the following norm on ,

 ∥q∥2H(div):=∥q∥2K−1+∥D⋅q∥2L2, (26)

where

 ∥q∥2K−1:=(q,q)K−1. (27)

Next lemma shows that the bilinear forms in the weak formulation (23)-(24) are continuous.

###### Lemma 2 (Continuity of (23)-(24))

For any and , we have

 (q,v)K−1 ≤∥q∥H(div)∥v∥H(div), (D⋅q,w) ≤∥q∥H(div)∥w∥L2.

###### Proof

The continuity of both bilinear forms follow directly from the Cauchy-Schwarz inequality and the definition of the norms (26) and (25).

Now we show the ellipticity of the inner product on the kernel of the mixed-dimensional divergence operator in the following lemma.

###### Lemma 3 (Ellipticity of (23)-(24))

If satisfies

 (D⋅q,w)=0,∀w∈L2, (28)

then

 (q,q)K−1=∥q∥2H(div). (29)

###### Proof

Since , from (28), we have

 ∥D⋅q∥L2=0.

Therefore, (29) follows directly from the above identity and the definition of the norm (26).

Next, we discuss the inf-sup condition of the bilinear form in the following lemma.

###### Lemma 4 (Inf-sup condition of (23)-(24))

There exists a constant such that, for any given function ,

 supq∈H0(div)(r,D⋅q)∥q∥H(div)≥β∥r∥L2. (30)

Here, the inf-sup constant depends on , the maximal number of overlaps between , structure of the trees , the domain , and the constants and .

###### Proof

Assume given, we first aim to construct such that .

First step is to construct based on the forest . Based on the signed incidence matrix  (10), we omit those columns that correspond to the Dirichlet root nodes to obtain the signed incidence matrix with boundary conditions . Then, we consider the following mixed-formulation graph Laplacian problem: Find and

 K−1qF−GFψF =0, (31) GTFqF =rF. (32)

Here, for trees , we set , , and for , we choose such that . The choice is not unique, and here we choose

 (rF)i=rNi−∑i∈NTrNi|NT,T|,i∈NT,T,T∈FN. (33)

For trees , we set , , and, for , we set

 (rF)i=rNi+1|NFD,T|∫ΩrD dx+1|NFD,T|∑i∈NFNrNi,i∈NFD,T.

The reason of such a choice will be made clear later in the proof when we construct . Note that, since the degree of node is one, once is fixed, we natrually have . With this choice of , the mixed-formulation graph Laplacian problem (31)-(32) is well-posed in the sense that is unique (up to a constant on the trees ) and is uniquely defined. Once is obtained, we define by , .

From the mixed-formulation (31)-(32), we have the following estimates,

 ∥GTFqF∥2=∥rF∥2and(K−1qF,qF)≤(λFmin)−1∥rF∥2, (34)

where

is the smallest non-zero eigenvalue of the weighted graph Laplaican of the forest

, i.e, . We comment that is bounded below by the so-called Cheeger constant of the graph, so depends on the structure of the trees in the forest . Note that

 ∥rF∥2=∑i∈NI∪NN(rNi)2+∑i∈NFD,T((rF)i)2+∑i∈NFN,T((rF)i)2

and, due the choice (33), the last term on the right-hand-side can be bounded by

 ∑i∈NFN,T((rF)i)2≤CN∑i∈NFN(rNi)2 (35)

with . Similarly, by Cauchy-Schwarz inequality, the second term on the right-hand-side can be bounded as follows,

 ∑i∈NFD,T((rF)i)2≤CD⎡⎢⎣∑i∈NFD,T(rNi)2+∑i∈NFN(rNi)2+∫Ω(rD)2dx⎤⎥⎦, (36)

where .

Therefore, combining the estimates (34), (35), (36), and the definitions of and , the following estimate holds,

 ∑e(i,j)∈E(kNe(i,j))−1|qNi,j|2≤CqN⎡⎣∑i∈NI∪NN(rNi)2+∑i∈NT(rNi)2+∫Ω(rD)2dx⎤⎦, (37)

where .

Next we construct from and so that (18) is satisfied exactly, i.e., we define, for each terminal nodes ,

 qSi(x)=qNNi,i+rNicki,x∈Bi. (38)

From the construction, we have

 ∑i∈NT∫Bi|qSi(x)|2dx ≤2|Bi|c2kS⎡⎣∑i∈NT|rNi|2+∑i∈NT|qNNi,i|2⎤⎦ ≤C1qS⎡⎢⎣∑i∈NT|rNi|2+∑i∈NFN|rNi|2+∫Ω(rD)2dx⎤⎥⎦, (39)

where . Here we use the fact that for by our construction of , and the estimates (35) and (36) in the last step. Similarly, we also have

 ∑i∈NT∫Bi|kSiqSi|2dx≤C2qS⎡⎢⎣∑i∈NT|rNi|2+∑i∈NFN|rNi|2+∫Ω(rD)2dx⎤⎥⎦ (40)

with .

Finally, we consider the following mixed-formulation Laplacian problem

 (kD)−1qD+∇ψ =0 (41) ∇⋅qD =rD+∑i∈NTkSiqSi (42)

with boundary condition on . This problem is well-posed because

 ∫ΩrDdx+∑i∈NT∫BikSi(x)qSi(x)dx =∫ΩrDdx+∑i∈NT∫BikSi(x)qNNi,i+rNickidx =∫ΩrDd