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.
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.
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.
The coefficient functions , , represent redistribution in a small region around the terminal node , thus, can be assumed to have compact support in some domain .
In practice the characteristic length scale of is comparable to the distance to the nearest neighbor, i.e.
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
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.,
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
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,
Similarly, we define the mixed-dimensional gradient as
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
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
Thus, , and we can introduce the rescaled divergence and gradients as and , respectively. The rescaled conservation equations are then summarized as
The rescaled conservation equations are summarized as
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 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,
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,
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,
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
By a direct calculation (using the re-scaled 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 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
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).
And the following norm on ,
Now we show the ellipticity of the inner product on the kernel of the mixed-dimensional divergence operator in the following lemma.
Next, we discuss the inf-sup condition of the bilinear form in the following lemma.
There exists a constant such that, for any given function ,
Here, the inf-sup constant depends on , the maximal number of overlaps between , structure of the trees , the domain , and the constants and .
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
Here, for trees , we set , , and for , we choose such that . The choice is not unique, and here we choose
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 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 , .
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
and, due the choice (33), the last term on the right-hand-side can be bounded by
with . Similarly, by Cauchy-Schwarz inequality, the second term on the right-hand-side can be bounded as follows,
Next we construct from and so that (18) is satisfied exactly, i.e., we define, for each terminal nodes ,
From the construction, we have
Finally, we consider the following mixed-formulation Laplacian problem
with boundary condition on . This problem is well-posed because