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 mixeddimensional when the network flow is simplified to a family of 1D domains along with the network edges^{1}^{1}1Some authors refer to this problem as multiscale (see e.g. [DangeloQuarteroni2008, KopplVidottoWohlmuthZunino2018]), however, we prefer the nomenclature mixeddimensional 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 subresolution 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 [HodnelandHansonSaevareidNaevdalLundervoldSolteszovaMuntheKaasDeistungReichenbachNordbotten2019] 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 [HodnelandHansonSaevareidNaevdalLundervoldSolteszovaMuntheKaasDeistungReichenbachNordbotten2019, 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 ”pointwise” 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 mixeddimensional 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 mixeddimensional Stokes’ theorem, integration by parts, and Hilbert spaces. This forms the building blocks for our wellposedness results and numerical analysis.
The main results of the paper are thus as follows:

A general, nonsingular model for a class of problems with a large dimensional gap.

Wellposedness theory for both the continuous and finitedimensional problem.

Convergence results for mixed finiteelement approximation and a finite volume variant.

Efficient linear solvers for different discretizations.

Numerical validation and application to a highresolution dataset 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 wellposedness. In Section 3 and 4 we state and analyze the finiteelement 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 networkDarcy 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 nonempty: 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 righthand 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 subforests, 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.
(1)  
(2)  
(3) 
Here, the signs in (1)(3) are chosen such that the righthandside 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,
Where the last step is also known as the GraphStokes’ 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 nonnegative (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.(4)  
(5)  
(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.
(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
(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.,
(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 infsup 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
(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 Mixeddimensional 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 mixeddimensional variables, on which we define mixeddimensional operators.
Therefore, let the mixeddimensional flux be denoted , and defined as the triplet of fluxes . Equivalently, the mixeddimensional pressure is denoted and defined as . Now, we define the mixeddimensional divergence operator as follows,
(11) 
where
(12) 
Similarly, we define the mixeddimensional gradient as
(13) 
where
(14) 
In addition, we introduce the function which contains all the material functions , , , and , , in (4) to (6), such that
where . It is now straight forward to verify that with these definitions, the conservation laws (1)(3) can be summarized as
(15) 
where . Furthermore, the constitutive laws (4) to (6) can be summarized as
(16) 
While the physical model formulation is satisfactory for nondegenerate , , 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 squareroot of the transfer coefficient , , and the scaled transfer mass flux , , is defined as . Thus, we replace (1) (3), and (6) with
(17)  
(18)  
(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 mixeddimensional flux as , and the scaling such that
Thus, , and we can introduce the rescaled divergence and gradients as and , respectively. The rescaled conservation equations are then summarized as
(20) 
The rescaled conservation equations are summarized as
(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
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 rescaled model, equations (20) and (21). Thus we will omit the superscript on the mixeddimensional 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 mixeddimensional squareintegrable space for pressure as follows,
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 mixeddimensional divergence as follows,
where is the space defined on such that the functions and their divergence are both squareintegrable. In addition, are standard space defined on , , and is the standard space defined on the edge set .
We associate the mixeddimensional space with the following inner product,
Similarly, we introduce the following inner product on ,
It is important to note that the inner products are defined such that integrationbyparts holds in for the mixeddimensional operators (both original and rescaled cases).
Lemma 1 (Integration by parts)
For any and , we have
(22) 
Proof
By a direct calculation (using the rescaled operators and variables, the derivation for the original case is the same), we have that
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,
where . In addition, with the material function , we introduce a weighted inner product on as follows,
Using the above function spaces and notation, together with the mixeddimensional 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
(23)  
(24) 
Note that due to the integration by parts formula, if nonhomogeneous boundary data is considered, this would appear as extra righthand side terms in equation (23).
2.4 Wellposedness
In this subsection, we focus on the wellposedness of the weak formulation (23)(24). We first introduce the following norm on ,
(25) 
And the following norm on ,
(26) 
where
(27) 
Proof
Now we show the ellipticity of the inner product on the kernel of the mixeddimensional divergence operator in the following lemma.
Proof
Next, we discuss the infsup condition of the bilinear form in the following lemma.
Lemma 4 (Infsup condition of (23)(24))
There exists a constant such that, for any given function ,
(30) 
Here, the infsup 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 mixedformulation graph Laplacian problem: Find and
(31)  
(32) 
Here, for trees , we set , , and for , we choose such that . The choice is not unique, and here we choose
(33) 
For trees , we set , , and, for , we set
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 mixedformulation graph Laplacian problem (31)(32) is wellposed 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 mixedformulation (31)(32), we have the following estimates,
(34) 
where
is the smallest nonzero eigenvalue of the weighted graph Laplaican of the forest
, i.e, . We comment that is bounded below by the socalled Cheeger constant of the graph, so depends on the structure of the trees in the forest . Note thatand, due the choice (33), the last term on the righthandside can be bounded by
(35) 
with . Similarly, by CauchySchwarz inequality, the second term on the righthandside can be bounded as follows,
(36) 
where .
Therefore, combining the estimates (34), (35), (36), and the definitions of and , the following estimate holds,
(37) 
where .
Next we construct from and so that (18) is satisfied exactly, i.e., we define, for each terminal nodes ,
(38) 
From the construction, we have
(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
(40) 
with .
Finally, we consider the following mixedformulation Laplacian problem
(41)  
(42) 
with boundary condition on . This problem is wellposed because
Comments
There are no comments yet.