MODELING, WELL-POSEDNESS AND DISCRETIZATION 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.

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 edges 1 . 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 [18,17]. A high dimensional gap thus arises when the flow in the network is connected to a domain of dimension d ≥ 2 through its leaf nodes, or when the flow in the network is connected to a domain of dimension d ≥ 3 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. This context is shown in Figure 1, which we will return to in the numerical results. 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 [13] and also considered in [14,21], and are attractive also from a mathematical perspective, as they avoid the singularities which otherwise characterize the coupled equations. In this work, Fig. 1. Illustration of a characteristic mixed-dimensional geometry associated with blood flow in the brain. This illustration is based on the the data-set used in the full-brain simulation study in section 5.3. The arterial tree is indicated in red, and the venous tree in blue. Note the complex geometry of the outer boundary of the brain (i.e. the domain Ω).
we will not adopt the precise models used in [13,14] 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. [10,15,11]). 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 [5], 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.
• 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. The theoretical results are validated in Section 5.
Finally, we give some conclusions in Section 6.

Model Equations.
In this section, we discuss the basic geometric setup and model equations for coupled network-Darcy flow in brain. We will both introduce the strong form and then derive the weak form by introducing proper spaces.
2.1. Geometry. We are concerned with a domain Ω ⊂ R n (which models the capillaries). In addition, we are concerned with a finite collection of rooted trees T with node (vertex) set N T and edge set E T (which model resolved arteries and veins). The arterial and venous trees are considered disjoint and, therefore, form a forest F with node set N = ∪ T ∈F N T and edge set E = ∪ T ∈F E T . We will refer to the composite (mixed-dimensional) problem domain of both Ω and F as the disjoint union B = Ω F.
We further distinguish the nodes of the forest as follows. The node set N can be subdivided into three disjoint subsets, the first and last of which are assumed to be non-empty: root nodes N R , interior nodes N I , and terminal nodes N T . Note that N = N R ∪ N I ∪ N T and we use N T ,R = N T ∩ N R , N T ,I = N T ∩ N I , and N T ,T = N T ∩ N T to denote the root nodes, interior nodes, and the terminal nodes of a given tree T , respectively. Naturally, we also have N T = N T ,R ∪ N T ,I ∪ N T ,T . We further divided the root nodes N R into two disjoint sets N D , which consists of the Dirichlet root nodes, and N N , 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, N T ,D and N T ,N denotes the Dirichlet or Neumann root nodes of a given tree T , 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 F D , which contains all the Dirichlet rooted trees T D , and Neumann rooted forest F N , which contains all the Neumann rooted trees T N . Naturally, We denote the set of the neighbors of the node i as N i and the set of all the edges meeting at i ∈ N as E i . Note that, |N i | = 1 and |E i | = 1, if i ∈ N R ∪ N T . These concepts are illustrated for n = 2 in Figure 2.
2.2. Strong From. As primary variables we choose the domain pressure potential p D (x) : Ω → R and the node pressure potentials p N : N → R. Furthermore, we consider the fluid mass fluxes denoted in the domain as q D (x) : Ω → R n , fluid mass flow from node i to j denoted q N i,j : N ×N → R and fluid mass flow transferring from terminal node i to point x denoted q T i (x) : Ω → R. 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) Here, the signs in (2.1)-(2.3) are chosen such that the right-hand-side terms represent sources added Fig. 2. Schematic illustrating the mixed-dimensional geometry B, and its subdivision into a continuous domain and a forest. Also shown is the coupling between the trees and the model domain.
to the system. Moreover, although both q N i,j and q N j,i are used in (2.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 (2.1)-(2.3) as follows, Where the last step is also known as the Graph-Stokes' Theorem, which is the counterpart of the Stokes' Theorem on graphs. We now propose constitutive laws for the flow. As our exposition primarily is concerned with the geometric complexity, we herein only consider linear constitutive laws, although it is reasonable that non-linear extensions may be required in applications (see e.g. [22]). We therefore introduce material parameters, all of which are assumed to be non-negative (precise bounds are given later). For each edge e(i, j) ∈ E, we assign a conductivity k N e(i,j) , which can be considered as the edge weights in certain sense. In the domain, for each x ∈ Ω we assign a permeability tensor k D (x) : Ω → R n×n . For each terminal node i ∈ N T , we assign connectivity function k T i (x) : Ω → R. Now based on the assumption that the potential flow is linear, we have the following constitutive laws.
(Potential flow in brain (Darcy)) q The coefficient functions k T i (x), i ∈ N T , represent redistribution in a small region around the terminal node i, thus, can be assumed to have compact support in some domain B i ⊆ Ω.
Remark 2.1. In practice the characteristic length scale of B i 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 h is the mesh size. The constants hidden in the O notation in (2.7) and (2.8) are usually between 2 to 10 in practical applications that we are interested in. Consequentially, q T i (x) also is compactly supported in B i . 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 N D , i.e., (2.9) q D · n = 0 on ∂Ω and p N i = 0, i ∈ N D . 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 q N i,j = −q N j,i . Therefore, although the total number of q N i,j is 2|E|, we only use half of them as the unknowns, i.e., one unknown, q N i,j or q N j,i , for each edge e(i, j) ∈ E. 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 e(i, j) ∈ E (i.e., if we choose q N i,j , which means the fluid mass flows from node i to node j, the edge e(i, j) is oriented such that it starts at node i and ends at node j). This allows us to define the following signed incidence matrix G ∈ R |E|×|N | , such that if flow on edge starts at node i −1, if flow on edge ends at node i 0, otherwise 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.3. Mixed-dimensional formulation and scaling. The model equations given above contain essentially three expressions of fluxes (q D , q T and q N ), and two expressions of potentials (p D and p N ). It will simplify the following exposition and analysis considerably to treat these as mixeddimensional variables on B, on which we define mixed-dimensional operators.
Therefore, let the mixed-dimensional pressure be denoted p : B → R, and defined as the doublet of pressures p := [p D , p N ]. Equivalently, the mixed-dimensional flux is defined as the triplet q := [q D , q T , q N ]. Now, we define the mixed-dimensional divergence operator D· as follows, Similarly, we define the mixed-dimensional gradient D as In addition, we introduce the function K which contains all the material functions k D , k T i (x), i ∈ N T , and k N e(i,j) , e(i, j) ∈ E, in (2.4) to (2.6), such that where k N = diag(k N e(i,j) ). It is now straight forward to verify that with these definitions, the conservation laws (2.1)-(2.3) can be summarized as where r ≡ [r D , r N ]. Furthermore, the constitutive laws (2.4) to (2.6) can be summarized as While the physical model formulation is satisfactory for non-degenerate k T i (x), i ∈ N T , it will be beneficial to rescale the coupling flux to avoid considering a degenerate inner product when k T i (x) → 0 for some points or region in B i . To that aim, we introduce the square-root of the respectively. In this setting, we allow for degeneracy of the coupling term in the sense that we allow k S i (x) → 0. However, we require that k S i is bounded from above, i.e., k S i (x) ≤ C k S for i ∈ N T and x ∈ Ω. Furthermore, for all i, we require it to hold that where c k S is a generic constant. We note that a similar scaling has been applied previously to handle degeneracies occurring in mantle dynamics [2] and flows in fractured porous media [6]. Equivalently, We denote the scaled mixed-dimensional flux on B as q S ≡ [q D , q S , q N ], and the scaling S such that Thus, q S = Sq, and we can introduce the rescaled divergence and gradients as D S · := D · S and D S := SD, respectively. The rescaled conservation equations are then summarized as The rescaled conservation equations are summarized as where K S = S −1 KS −1 , and thus Note in particular that (K S ) −1 applied to q S has unit weight, and therefore does not degenerate.

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 (2.15) and (2.16), as well as the re-scaled model, equations (2.20) and (2.21). Thus we will omit the superscript S 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 q S , and consider the coefficient k S to appear in the differential operator as opposed to the material law.
We first introduce proper function spaces on B. We begin by defining a mixed-dimensional square-integrable space for pressure as follows, where L 2 (Ω) is the standard L 2 space defined on domain Ω and l 2 (N \N D ) is the standard l 2 space defined on the node set N \N D . For flux, we consider a space with bounded mixed-dimensional divergence as follows, where H(div, Ω) is the space defined on Ω such that the functions and their divergence are both square-integrable. In addition, L 2 (B i ) are standard L 2 space defined on B i , i ∈ N T , and l 2 (E) is the standard l 2 space defined on the edge set E.
We associate the mixed-dimensional space L 2 (B) with the following inner product, Similarly, we introduce the following inner product on H(div, B), 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 2.2 (Integration by parts). For any q ∈ H(div, B) and p ∈ L 2 (B), we have Proof. 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 q D · n = 0 on ∂Ω, therefore, we define the following functions space with boundary conditions, where H 0 (div, Ω) := {q D ∈ H(div, Ω) | q D · n = 0, on ∂Ω}. In addition, with any material function K, we introduce a weighted inner product on H(div, B) as follows, Using the above function spaces and notation, together with the mixed-dimensional integration by parts formula (2.23) and the homogeneous Dirichlet boundary conditions (2.9) on N D , i.e., p N i = 0, i ∈ N D , we have the following weak form for the conservation laws (2.20), and constitutive laws (2.21): Find q ∈ H 0 (div, B) and p ∈ L 2 (B), 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 (2.24).
2.5. Well-posedness. In this subsection, we focus on the well-posedness of the weak formulation (2.24)-(2.25). As in the previous subsection, it is understood that we are considering the re-scaled formulation, even though the superscript S is suppressed. We first introduce the following norm on L 2 (B), And the following norm on H(div, B), We emphasize that the weights in this norm do not degenerate for the re-scaled equations sine the unite weight is applied to q S , see equation (2.22). The next lemma shows that the bilinear forms in the weak formulation (2.24)-(2.25) are continuous.
Lemma 2.3 (Continuity of (2.24)-(2.25)). For any q, v ∈ H(div, B) and w ∈ L 2 (B), we have Proof. The continuity of both bilinear forms follow directly from the Cauchy-Schwarz inequality and the definition of the norms (2.27) and (2.26). Now we show the ellipticity of the inner product (·, ·) K −1 on the kernel of the mixed-dimensional divergence operator D· in the following lemma.
Proof. Since D · q ∈ L 2 (B), from (2.29), we have Therefore, (2.30) follows directly from the above identity and the definition of the norm (2.27).
Next, we discuss the inf-sup condition of the bilinear form (r, D · q) in the following lemma.
Lemma 2.5 (Inf-sup condition of (2.24)-(2.25)). There exists a constant β > 0 such that, for any given function r ∈ L 2 (B), Here, the inf-sup constant β depends on |B i | = measure(B i ), the maximal number of overlaps between B i , structure of the trees T ∈ F, the domain Ω, and the constants c k S and C k S .
First step is to construct q N based on the forest F. Based on the signed incidence matrix G (2.10), we omit those columns that correspond to the Dirichlet root nodes to obtain the signed incidence matrix with boundary conditions G F . Then, we consider the following mixed-formulation graph Laplacian problem: The choice is not unique, and here we choose For trees T ∈ F D , we set (r F ) i = r N i , i ∈ N T ,T ∪ N T ,I , and, for i ∈ N T ,T , we set The reason of such a choice will be made clear later in the proof when we construct q D . Note that, With this choice of r F , the mixed-formulation graph Laplacian problem (2.32)-(2.33) is well-posed in the sense that ψ F is unique (up to a constant on the trees T ∈ F N ) and q F is uniquely defined.
From the mixed-formulation (2.32)-(2.33), we have the following estimates, where λ F min is the smallest non-zero eigenvalue of the weighted graph Laplaican of the forest F, i.e, L F = G T F KG F . We comment that λ F min is bounded below by the so-called Cheeger constant of the graph, so depends on the structure of the trees T in the forest F. Note that and, due the choice (2.34), the last term on the right-hand-side can be bounded by Similarly, by Cauchy-Schwarz inequality, the second term on the right-hand-side can be bounded as follows, Therefore, combining the estimates (2.35), (2.36), (2.37), and the definitions of G F and q N , the following estimate holds, . Next we construct q S from q N and r N so that (2.18) is satisfied exactly, i.e., we define, for each terminal nodes i ∈ N T , From the construction, we have Here we use the fact that q N Ni,i = −(r F ) i for i ∈ N T by our construction of q N , and the estimates (2.36) and (2.37) in the last step. Similarly, we also have . Finally, we consider the following mixed-formulation Laplacian problem with boundary condition q D · n = 0 on ∂Ω. This problem is well-posed because which verifies the consistency of the data with respect to the pure Nuemann boundary condition q D · n = 0 on ∂Ω. Furthermore, the following estimate holds, has been constructed based on [r D , r N ] and it satisfies and we have Now, based on (2.38), (2.40), and (2.45), we can derive that Remark 2.6. The inf-sup proof shows the importance of the using the scaled equations (2.20) and (2.21) in the case where k T goes to zero. Indeed, for the non-scaled equations, a similar approach would lead to an inf-sup constant depending on the pointwise lower bound on inf x∈Bi (k T i (x)), which may not be positive. In contrast, as seen in the proof above, for the scaled equations, inf-sup constant depends on the much less restrictive integrated bound c k S i . We now have the following well-posedness results. Proof. The result follows directly from the standard theory for saddle point problems, see, e.g. [4], and Lemmas 2.3, 2.4, and 2.5.
3. Finite-element Approximation. In this section, we propose the finite-element approximation for solving the weak formulation (2.24)-(2.25). The coupling between the graph and the porous domain, as well as the heterogeneous nature of the parameters found in applications, suggests that it is natural to consider low-order approximations. As a consequence, we only consider the lowest-order approximation here, recognizing that higher-order spaces can be introduced in the mixed formulation.
3.1. Mixed Finite-Element Method. Given a mesh M of the domain Ω, e.g., triangles/quadrilaterals in 2D and tetrahedrons/cuboids in 3D, we consider the standard RT0/P0 finite element for approximating the fluid flux q D and pressure P D in the domain and denote them by H h (div, M) and P 0 (M), respectively. For node pressure potentials p N , we use vertex degrees of freedom (DOFs) of the graph. For fluid flux on the tree edges, we use edge DOFs of the graph. For the fluid flux transferring from terminal i to point x, it appears natural to consider the piecewise constant finite element on M i (denoted as P 0 (M i )), which is the restriction of M to B i , i.e. M i = M ∩ B i . In summary, we consider the following conforming finite-element spaces its corresponding finite-element space with boundary conditions, Using the finite-element spaces introduced above, the mixed finite-element approximation of (2.24)- Remark 3.1. By considering a test function w h which is constant on a B i , we verify from equations 2.17 and 2.18 that the physical flux q T = k S i q S i is conserved. We note that the lowest-order mixed finite element approximation is locally conservative even when applied to scaled variables, in contrast to the situation observed when similar scalings are applied in the physical dimensions of Ω (see e.g. [2]).

Well-posedness.
In this subsection, we consider the well-posedness of the mixed finiteelement approximation (3.1)-(3.2). It is essentially the same as the well-posedness analysis for the weak formulation in Section 2.5, and our presentation will therefore be brief.
Since we use conforming finite-element spaces, the continuity results (Lemma 2.3) holds naturally on the discrete level.
For the ellipticity (Lemma 2.4), using the fact that the finite dimensional spaces are conforming in the sense that for q h ∈ H h,0 (div, B), then it holds that D · q h ∈ L 2 h , then the continuous results hold on the discrete level.
Moreover, the inf-sup condition (Lemma 2.5) can be derived in a similar fashion on the discrete level as well.  2)). There exists a constant β > 0 such that, for any given function r h ∈ L 2 h (B), Here, the inf-sup constant β depends on |M i | = measure(M i ) = O(h n ), the maximal number of overlaps between B i , structure of the trees T ∈ F, the domain Ω, and the constants c k S and C k S . B) is similar to the construction presented in the proof of Lemma 2.5. q N h can be constructed exactly the same as the construction of q N . Then q S h can be defined as (2.39) as well since such construction also makes sure that q S h ∈ i∈N T P 0 (M i ). The construction of q D h should be obtained by solving (2.42)-(2.43) with a mixed finite-element method using H h,0 (div, M) and P 0 (M). Such construction also makes sure that  3.3. Convergence. Based on Lemma 3.2, 3.3, and 3.4 and applying the general theory of Galerkin methods, see [7,4], we immediately gives a quasi-optimality error estimate.
As usual, to obtain the final convergence result, we use interpolations to bound the right-handside of the above error estimate (3.6). Here, we choose v D h = π div q D , where π div : H 1 (Ω) → H h (div, M) is the standard interpolation given by the H h (div, M) degrees of freedom, v S h = π 0 q S , where π 0 denotes the standard piecewice constant interpolation, and v N h = q N . With those choices and the classical error estimates for interpolations, together with Cauchy-Schwarz inequality, we naturally have Similarly, by choosing w D h = π 0 p D and w N h = p N , we have . Therefore, we have the overall convergence result for the finite-element method (3.1)-(3.2) as follows.
Corollary 3.7. Suppose that q ∈ H 0 (div, B) and p ∈ L 2 (B) satisfy the weak formulation (2.24)-(2.25), then the finite-element solution q h ∈ H h,0 (div, B) and p h ∈ L 2 h of the mixed finite-element approximation (3.1)-(3.2) satisfy that where the constant c depends on β and the quasi-uniformity of the mesh M.
Remark 3.8. In Corollary 3.7, we require ∇ · q D ∈ H 1 (Ω) because the convergence analysis is derived by following the standard Babuška-Brezzi theory. As it is well-known for the error analysis of the fixed-dimensional mixed finite-element method for second-order elliptic problem, this regularity requirement can be relaxed in the mixed-dimensional setting as well, i.e., we have the following error estimates Due to space constraints, we omit the derivation here but comment that it is essentially the same as the derivation for the fixed-dimensional case as shown in [4].

Mass Lumping and Two-Point Flux Approximation Scheme.
In practice, when the triangulation of the domain Ω is uniform, it is possible to simply the discretization scheme and use two-point flux approximation (TPFA) to discretize the PDE system given by the conservation laws (2.17), (2.2), (2.18) and the constitutive laws (2.4), (2.5), (2.19). This is particularly relevant for medical applications, where the data is frequently specified on voxels (i.e. regular Cartesian grids in 3D).
In this section, we, therefore, discuss the TPFA scheme for our coupled Network-Darcy model through its relationship with the mixed finite-element approximation (3.1) and (3.2) discussed in Section 3.
4.1. TPFA Scheme. On a given mesh M, similar to standard diffusion problems, the TPFA scheme can obtained by applying mass lumping to the mixed finite-element scheme (3.1)-(3.2) and then eliminating the flux q h . To this end, we define the following inner product on the finite element spaces H h (div, B), for q h and v h ∈ H h (div, B), being the average of k D on the element τ ∈ M and d f being the distance between the face f ∈ ∂τ and the cell center of τ . Now we define the mass lumping finite-element scheme as follows: Find q h ∈ H h,0 (div, B) and p h ∈ L 2 h (B), such that, Based on the inner product (4.1), we define a discrete gradient D h : L 2 h (B) → H h (div, B) via integration by part (Lemma 2.2), for any v h ∈ H h (div, B) and p h ∈ L 2 h (B), such that, Note that, due to the boundary conditions, v D h · n = 0 on ∂Ω and (p N Then the mass lumping mixed-formulation (4.2) and (4.3) can be written as, find q h ∈ H h,0 (div, B) and p h ∈ L 2 h (B), such that, The above formulation allows us to eliminate q h and obtain the TPFA scheme as follows, find Next we will explain the TPFA scheme (4.4) using matrix notation. The matrix form of the mass lumping finite-element scheme (4.2)-(4.3) can be written as  Since D D , D s , and D N are diagonal matrices, we can eliminate them by block Gaussian elimination and end up with a linear system only involves solving for p D h and p N h as follows, which is exactly the matrix form of the TPFA scheme (4.4).

Well-posedness.
Next we consider the well-posedness of the TPFA scheme (4.4). As we showed in the previous section, the TPFA scheme (4.4) is obtained from the mass lumpping mixed-formulation (4.2)-(4.3) by block Gaussian elimination. Therefore, we first show the wellposedness of the mass lumpping mixed-formulation (4.2)-(4.3) and then the well-posedness of the TPFA scheme (4.4) follows directly.
Since the only difference between the mixed-formulation (3.1)-(3.2) and the mass lumpping mixed-formulation (4.2)-(4.3) is the inner product used for H h (div, B), we first introduce the norm induced by the inner product (4.1) as follows, and show it is spectrally equivalent to the norm (2.28) in the following lemma.
where c 1 > 0 and c 2 > 0 are constants only depending on the shape regularity of the mesh M.
Proof. Based on the standard result, e.g., [12], we havē where the positive constantsc 1 andc 2 depend only the shape regularity of the mesh M. Then the spectral equivalence (4.5) follows directly from the definitions of the norms.

Define
We have the following lemmas concerning the continuity, ellipticity, and inf-sup condition for the mass lumping mixed-formulation (4.2)-(4.3).
Lemma 4.2 (Continuity of (4.2)-(4.3)). For any q h , v h ∈ H h (div, B) and w h ∈ L 2 h , we have For the ellipticity, again using the fact that, for q h ∈ H h (div, B), D · q h ∈ L 2 h , we have Lemma 4.3 (Ellipticity of (4.2)-(4.3)). If q h ∈ H h (div, B) satisfies Moreover, the inf-sup condition can be derived from the inf-sup condition (Lemma 3.4) and the spectral equivalence lemma (Lemma 4.1) Lemma 4.4 (Inf-sup condition of (4.2)-(4.3)). There exists a constant β > 0 such that, for any given function r h ∈ L 2 (B), Here, the inf-sup constant β depends on |M i | = measure(M i ) = O(h n ), the maximal number of overlaps between B i , structure of the trees T ∈ F, the domain Ω, the constants c k S and C k S , and the shape regularity of the mesh M.
Proof. The inf-sup condition (4.7) can be derived from the inf-sup condition (3.5) and the spectral equivalence result (4.5).

Convergence.
Regarding the convergence result of the TPFA scheme, since we use mass-lumping technique to derive it, existing theoretical tools developed in [3,8] can be adopted here. For the sake of the simplicity, in this subsection, we assume that k D is constant on each element τ ∈ M and the mesh M is uniform (e.g., rectangle/equilateral triangle in 2D, rectangular cuboid/regular tetrahedra in 3D). Under those conditions, as shown in [3], for τ ∈ M, f ∈∂τ ω f q D · n f v D · n f used in the definition (4.1) provides a numerical integration formula of τ (k D ) −1 q D v D dx and such a numerical integration is exact for constant functions on each element τ . Moreover, the following perturbation result holds for q D , v D ∈ H h,0 (div, M), Based on the above result, we can easily verify that, for q h , v h ∈ H h (div, B), Now, we can use the theory developed in [19] and conclude the convergence result of the TPFA scheme in the following theorem.
Theorem 4.7. Suppose that q ∈ H 0 (div, B) and p ∈ L 2 (B) satisfy the weak formulation (2.24)-(2.25), then the finite-element solution q h ∈ H h,0 (div, B) and p h ∈ L 2 h of the mass lumping mixed finite-element approximation (4.2)-(4.3) satisfy that where the constant c depends only on β, k D , the maximal number of the overlap between M i , max i {|M i |}, and quasi-uniformity of the mesh M.
Consequentially, this also implies the convergence result of the TPFA scheme because of the equivalence between the TPFA scheme (4.4) and the mass lumpping mixed-formulation (4.2) and (4.3).
Remark 4.8. As pointed out in Remark 3.8, the regularity requirement ∇ · q D ∈ H 1 (Ω) can be relaxed here as well and similar convergence analysis still holds. Remark 4.9. As shown in [3,8], similar results hold for some more general meshes. For example, the perturbation result (4.8) hold for general triangles in 2D with order h instead of order h 2 . However, this still leads to the error estimate (4.9) based on the same procedure. For general triangulation in 3D, convergence analysis for standard mixed-formulation Poisson problem with mass lumping was derived based on a different approach in [8]. We can also adopt a similar approach to derive the convergence result for our mass lumping mixed finite-element scheme as well to obtain the error estimate (4.9) for general triangulation as well.
5. Numerical Results. In this section, we include three numerical results to validate and explore the discretization and solver presented above. In particular, the first case contains the simplest possible geometry in 2D, on which we compare the discretization to a series solution (Bessel functions). In the second case, we have a more complex geometry embedded in 4D, which can be seen as a prototype of the geometries relevant for applications. In both the first and second cases, we perform convergence studies both for the discretization and multigrid solver. Finally, in the third case, we apply the methodology to a real dataset, based on the human brain.
The error is measured in the norms proposed in the analysis, in particular we measure the L 2 norm of pressure and the k −1/2 -weighted norm of flux. As is common for finite volume and mixed finite-element methods, we use cell-centered quadrature when evaluating the L 2 norm in the domain, which allows us to exhibit the usual super-convergence behavior for these methods on smooth problems.
Due to the prevalence of image data for the applications of interest, all the numerical experiments are conducted on uniform Cartesian grids and the TPFA scheme is used. To solve the resulting linear system, we use algebraic multigrid (AMG) preconditioned flexible GMRes (FGM-Res) method, as detailed in the Supplementary Materials ??. Here, an unsmoothed aggregation AMG method is used as the preconditioner. More precisely, one step of V-cycle AMG method is applied with one step of Gauss-Seidel method for both pre-and post-smoothing. The FGMRes method is terminated when the 2 -norm of the initial residual is reduced by a factor of 10 −6 . The solver performance for all three cases below is also reported in the Supporting Information. The implementations are in Matlab, and code is available from the authors on request. All runs are conducted on a Linux workstation using 40 Intel Xeon CPU processors (E5-2698 v4) at 2.20GHz clock speed, with 256 Gb RAM.

Case 1: Comparison to Convergent Series Solution.
Our first case is constructed such that a series solution (in terms of well-known Bessel functions) is available. The full derivation of the series solution is available in the Appendix, an illustration of the geometry, and the series solution is provided in Figure 3. Throughout this subsection, we consider the series solution as the exact solution of the equations, since arbitrary precision can be obtained using well-established implementations of table values [1].
The main features of the solution is a simple two-node tree, where node 0 is a Dirichlet boundary node, and node 1 is a terminal node. Correspondingly, there is a single edge in the network, which contains the network flux. The solution is constructed with a transfer function k T that has compact support on a disc of radius r 1 from the origin. We consider two variants of the case, case 1A has a smoothly degenerating transfer function such that (in terms of radial coordinates) k T (r) → 0 as r → r 1 , while case 1B has a constant k T within the disc (and zero outside), thus k T ∼ H(r 1 − r), where H denotes the Heaviside function. To drive the system, a quadratic source term is provided in the region r 2 < r ≤ r 3 .
We conduct numerical experiments with unit values, such that the domain Ω is the unit square centered at the origin, the domain and network permeabilities are unit valued, and the scaling of source term r D = 1. The transfer function k T has a unit maximum value at the origin, for both case A and B, thus in the notation of the appendix k T 0 = 1. As stated, we consider two versions of the case. For the case 1A, we consider a degenerating transfer function k T , with r 0 = 0.1, r 1 = 0.2, r 2 = 0.3, r 3 = 0.4. For case 1B, we let the transfer function abruptly go to zero by keeping all radii as in case 1A, except for r 0 = 0.2.
An important aspect of the implementation is the accuracy with which the right-hand-side and the inner products involving k S are evaluated. In the results reported here, we have used a fourth-order accurate numerical quadrature.
The convergence results of cases 1A and 1B are presented in Table 1 and 2. We show the convergence history separated into components similar to the analysis, i.e. Domain, Scaled terminal flux, and Network.
First note that for this example, since the network contains a single throat and the domain has Neumann boundary conditions, global conservation of mass implies that q N h will be exact up to the quadrature error in the evaluation of r D , and similarly for p N h . Thus the fourth-order convergence of these variables is expected.
As for the remaining variables, we observe in both Case 1A and Case 1B optimal second-order convergence of p D h and first-order convergence of q D h . In this example, the scaled terminal flux q S h is essentially just the weighted difference between p D h and p N h , and thus it inherits the (slower) convergence rate of the two, i.e. second-order. By comparing the two cases, we see that there is no influence of the degeneracy of k S .  5.2. Case 2: A Prototypical 4 Dimensional Case. Our second example is chosen to illustrate a typical case encountered in the modeling of tissue. The physical domain is 3-dimensional, however, due to the biomedical properties involved, the physical domain represents two or more continua (biomedically speaking, this corresponds to arterial and venal compartments, etc.). The continua are ordered, and communication between the compartments is only allowed between neighbors in the ordering. As such, the continua represent a discretization of an elliptic equation in a fourth dimension. The mathematical structure of the resulting system is thus one of a 4D elliptic equation, coupled to networks, and is naturally covered by the methods proposed analyzed in this paper.
To explore this concept, and validate the performance of our methods, we consider the following concrete problem, as illustrated schematically in Figure 4. Let the model domain be the unit 4cube. We consider Neumann boundary conditions on all faces of the domain. Furthermore, we consider two trees, which are named as "arterial tree" and "venous tree", respectively, to conform with applications and the next subsection. Each consists of four nodes connected in the shape of a "Y", wherein each tree, one node is a Dirichlet boundary node (p N D = 1 and p N D = 0 in arterial and venous Dirichlet nodes, respectively), while two nodes are terminal nodes. The arterial terminal nodes i are associated with transfer functions is the distance in the first three coordinates from the 3-points y i ,  We discretize the domain with an anisotropic Cartesian grid in the sense that the first three dimensions are discretized by a regular isotropic Cartesian grid. The fourth dimension is discretized by only two grid cells. This resulting system is equivalent to the common two-compartment model, where the cells in the fourth dimension with x 4 < 0.5 correspond to the arterial compartment, and the remaining cells the venous compartment. In accordance with the practice in applications, we will emphasize grid refinement over model refinement, and only consider refinement of the first three dimensions. Moreover, we will in accordance with the applications decompose the domain flux into two parts q D → [q D , q P ], where the flux in the fourth dimension q P is referred to as "perfusion". Model parameters are otherwise set to unity, k D = k P = k N = 1 where k P is the permeability constant of the flux in the fourth direction.
The convergence results for this case are presented in Table 3. All errors are reported relative to a numerical solution calculated with a resolution of h = 256 −1 , and convergence rates are therefore reported for grids up to a resolution of h = 128 −1 . As expected, we observe quasioptimal convergence rates in all variables. In contrast to case 1, we no longer have the artificial exact solutions in the network, where we observe second order convergence, as inherited from the interaction between the terminal nodes and the second-order accurate pressure in the domain.

Case 3:
Full-brain simulation study. As a final test case, we consider the application to a real data set, associated with blood flow in the human brain. As a modeling concept, we use the same general structure as illustrated in Figure 4. The data-set and parameterization is described in detail in [13], and is illustrated in Figure 1 of the introduction. Here we summarize the main features: the data contains two trees, corresponding to a segmentation of the arterial and venous systems, containing 355 and 1222 nodes, respectively. For the finest simulations, we consider full resolution MRI acquisitions, which after co-registration to the finest resolution image is a Cartesian grid with 346×448×319 grid cells, representing a brick-shaped field of view of 177×224×160mm 3 . The actual domain Ω is a 4D extrusion of the 3D subset of the field of view from a T1-weighted MR acquisition which contains segmentation of the brain acquired with the human brain segmentation software FreeSurfer [9]. Thus the mathematical formulation is a 4D model in the sense of the previous sub-section, and after discretizing the fourth dimension by two cells, the full model contains 17.5 million grid cells. The domain Ω is furthermore divided into two subdomains by the FreeSurfer segmentation (anatomically: white matter Ω W M and gray matter Ω GM ), with permeability in the three physical dimensions set to an isotropic value of k D = 10 −11 m 2 . The permeability k D acting in the 4th dimension (anatomically: the perfusion coefficient), is anisotropic relative to the physical dimensions, and is in the white matter set to k P = 10 −6 m · s · kg −1 , x ∈ Ω W M , and in grey matter is set to k P = 1.6 · 10 −6 m · s · kg −1 , x ∈ Ω GM . The transfer permeability is set according to equation (??), with r 1 = 30mm, r 0 = r 1 /2, and k T 0 = 10 −4 . The arterial and venous vessel trees are extracted down to voxel resolution from time-of-flight (TOF) and quantitative susceptibility mapping (QSM), respectively. Within both these MR acquisitions, a crude segmentation of the vessels is obtained by local adaptive thresholding, leading to a large number of disconnected structures. These binary satellites are connected with the main structure by repeatedly solving a boundary value problem around the main structure S for each satellite. Hence, the solution of the Eikonal equation |∇T | = f (x) −1 , T (x ∈ S) = 0 for the arrival time T (x) provides a geodesic distance map from x to the main structure. The Eikonal equation was solved using the fast marching method [20]. The function f (x) is user-supplied and is known as the speed of the arrival time field. Within TOF we use the image itself as speed function and for QSM the inverted image due to low contrast within vessels. The speed is set to zero outside the brain, possibly leading to curved geodesic trajectories, which is the reason why the signed distance function is not used. The arrival time itself is not of interest here, but rather the backtracing in the arrival time field from the satellite to the main structure providing a most probable path connecting these two structures with each other. The current approach favors probable paths to be aligned with dark-or bright-contrast structures that partly disappear within the images due to noise in the data. While backtracing, visited points are added to the main structure with a suitable vessel radius. Finally, the process of solving the Eikonal equation is repeated for each satellite, ultimately providing a connected structure, i.e. the arterial or venous vessel tree. For a more comprehensive description of how satellites are connected to the main structure, we refer the readers to [13].
The now connected binary trees are converted into abstract graphs using built-in Matlab routines for skeletonization, leaf (terminals and roots), and node detection. Vessel length is the geodesic distance along the edge between two connecting nodes, and the average vessel diameter is fitted by a Euclidean distance function around the centerline. The edge flow permeability k N is assigned individually for each edge based on Hagen-Poiseulle's law, using local estimates of vessel diameter and vessel length measured in the binary vessel trees. Both arterial and venous trees are modeled with Dirichlet root nodes as the main arterial inlets and main venous outlets. The only properties of the vessel trees that are needed for the simulation experiments are the edge flow permeability The full brain data contains many important qualitative properties, including connectivity of the trees after preprocessing of the initially disjoint trees, and connectivity of the brain geometry. These properties ensuring well-posedness, as well as the connected representation of grey and white matter, are not trivially preserved when coarsening the data. Thus instead of reporting relative results on a grid sequence for this case (which due to the above would have limited real value), we summarize the calculated solution on the image resolution in Figure 5 (the subfigures of this figure are shown in full size in section SM3 of the supporting information). While the quantitative aspects of the calculated results depend on parameters that are at present not fully justified by clinical measurements, our calculations verify that the proposed methods allow for efficient simulations at imaging resolution, preserving the qualitative properties of the solution corresponding to biomedical expectations.
6. Conclusions. We have proposed a mixed-dimensional mathematical model, closely related to models used for modeling fluid flow in human vasculature. We show the well-posedness of this model on the continuous level and develop suitable numerical discretizations, of both mixed finiteelement and finite volume types. These are shown to be stable and convergent.
Our theoretical results are complemented by numerical examples, which demonstrate superconvergence of the method in terms of the pressure variable on smooth solutions, and also verifies the stability and applicability of the method to large scale real-world data sets. In this section, we discuss how to solve the linear systems obtained by the TPFA scheme from the main part of the paper. The TPFA scheme solves the unknown p h first, after which the unknowns q h can be recovered. Therefore, we discuss multigrid methods for the TPFA scheme. In this section, we will solely be dealing with the discrete system, and thus omit the subscripts on q h and p h .
The TPFA scheme (??), following the discussion in Section ??, can be obtained by block Gaussion elimination of the mass lumping mixed-formulation, i.e., its compact matrix form is As we can see, this is a discretization of a Poisson type problem for the unknown p in the mixeddimensional setting. Since multigrid (MG) methods are designed for solving diffusion type problems efficiently [SM8,SM10], in our work, we use an algebraic variant of the MG methods to solve the resulting linear system. More precisely, we use the aggregation-based AMG method, which is a suitable choice for solving both Poisson-type problems, see [SM2,SM6], and graph Laplacian problems, see [SM3,SM7,SM4].
To validate the efficiency of multi-grid for this problem, we include the solver performance on the three test cases from the manuscript in the tables below.
The performance of the iterative solver for Case 1A is summarized in Table SM1. The results for Case 1B are essentially identical and are omitted (the deviation is less than 5% in all quantities) . This test case, due to the simplicity of the network, is from a linear algebra perspective almost identical to the 2D homogeneous Poisson problem. Therefore, the aggregation-based AMG preconditioner performs as expected. Here, the grid complexity (GridComp.) is defined as the ratio between the total size of the matrices on different levels and the size of the matrix on the finest level. The operator complexity (Oper.Comp.) is defined as the ratio between the total number of nonzeros of the matrices on different levels and the number of nonzeros of the matrix on the finest level. As we can see, both Grid Complexity and Operator Complexity are bounded by 2, which means that the aggregation-based AMG method roughly preserves the sparsity on the coarse levels and the overall computational cost of one V-cycle AMG method is optimal (i.e. linearly proportional to the cost of one step of Gauss-Seidel method). Since we use a simple unsmoothed aggregation AMG method here, the number of iterations grows with respect to the size of the linear system as expected. Overall, the computational cost of the multilevel solver is sub-optimal. We want to point out that it is possible to obtain optimal computational complexity when a more sophisticated AMG  cycle, such as the K-cycle is used [SM9,SM5]. This is a subject of our ongoing and future research since developing special tailored multilevel solvers for mixed-dimensional problems is a challenging and important task itself and beyond the scope of this work.
The solver results for case 2 are presented in Table SM2, including the performance on the reference solution, which has more than 33 million degrees of freedom. Again, we observe that both grid complexity and operator complexity are bounded while the number of iterations grows sublinearly with respect to the size of the linear systems. This means the overall computational complexity is sub-optimal, which is expected since we are using a simple unsmoothed aggregation AMG preconditioner.
The computational performance on Case 3, which consists of real data, is summarized in Table  SM3. All results are reported relative to the discretization level of the data, which we label h 0 . The results confirm that even for this large data-set, with data-driven and physuiologically realistic parameters, the solver performs efficiently. The number of iteration grows sublinearly and the CPU time also indicates that the current solver is sub-optimal. However, we do not see the degradation of the AMG solver comparing with the previous two cases with simpler geometry, which demonstrates that our discretization is stable and solver-friendly. Moreover, the efficiency of the discretization and solver combinations proposed above allow for analysis of full-resolution medical data.
SM2. Derivation of reference solution for section 5.1. The solution is valid for any domain Ω containing the origin and with zero Neumann boundary conditions, and for a tree consisting of two nodes: Node 0 is a Dirichlet node and node 1 is a terminal node. Please refer to Figure ?? for a visual illustration.
Setup. Consider polar coordinates (r, θ) around the origin, included in Ω. In addition to Neumann boundary conditions on Ω, the boundary pressure p N 0 in the Dirichlet node is taken as a given, and we consider for simplicity zero source terms in the terminal node, r N 1 = 0. Then we construct up to four concentric domains, bounded by 0 < r 0 ≤ r 1 ≤ r 2 < r 3 ≤ diam (Ω). The inner domain will be a disc, the remaining domains will be doughnuts. We consider homogeneous k D (x) = k D , while the remaing parameters in Ω are strictly functions of r. In

Table SM3
Solver report whole brain simulation. Note that we set h 0 as the resolution of the finest grid, since this is the resolution of the data. Also, the DOF column reflects that due to the nontrivial shape of the brain, the number of active cells in the computation does not scale exactly with the cube of the grid resolution. particular, we choose a transfer permeability which is non-zero only in the two inner-most regions It then follows that the pressure p N 1 at the terminal node is given by The combination of equations (SM2.4) and (SM2.5) have analytical solutions within each region of the domain. These are as follows (these are reported in numerous text books, we will use Chapter 9 of Abramovich and Stegun as a concrete reference for the use and manipulation of Bessel functions [?]).
Solution S1 : For r ≤ r 0 : By substitution of equation ( Constants. It remains to determine the constants c I , c J and c Y from solutions S1 and S2. The solutions S1 and S2 must satisfy three criteria: 1) The pressure must be continuous at r 0 , 2) the flux must be continuous at r 0 , and 3) The flux must be continuous at r 1 . This leads to three linear constraints on the constants.
These criteria can be made precise using the solutions derived above as follows.
For the special case where r 0 = r 1 , two of the constants, c J and c Y are superfluous, and the remaining constant is obtained by combining equations (SM2.28) and (SM2.29) as: (SM2.31) c I k D κ 0 I 1 (κ 0 r 1 ) = − q N 2π 1 r 1 SM3. Full-size figures for section 5.2. This section contains enlarged versions of the subfigures in Figure 5 of the article, which contains the colorbars for the figures in this section.