Nonlinear mixed-dimension model for embedded tubular networks with application to root water uptake

We present a numerical scheme for the solution of nonlinear mixed-dimensional PDEs describing coupled processes in embedded tubular network system in exchange with a bulk domain. Such problems arise in various biological and technical applications such as in the modeling of root-water uptake, heat exchangers, or geothermal wells. The nonlinearity appears in form of solution-dependent parameters such as pressure-dependent permeability or temperature-dependent thermal conductivity. We derive and analyse a numerical scheme based on distributing the bulk-network coupling source term by a smoothing kernel with local support. By the use of local analytical solutions, interface unknowns and fluxes at the bulk-network interface can be accurately reconstructed from coarsely resolved numerical solutions in the bulk domain. Numerical examples give confidence in the robustness of the method and show the results in comparison to previously published methods. The new method outperforms these existing methods in accuracy and efficiency. In a root water uptake scenario, we accurately estimate the transpiration rate using only a few thousand 3D mesh cells and a structured cube grid whereas other state-of-the-art numerical schemes require millions of cells and local grid refinement to reach comparable accuracy.


Introduction
Nonlinear elliptic equations arise in the description of fluid flow in porous media where permeability depends on the water pressure (e.g. Richards' equation) or the description of heat conduction where the thermal conductivity depends on temperature. In this work, we discuss a numerical scheme to solve such equations in the presence of an embedded thin tubular transport system exchanging mass or energy with the embedding bulk domain. This exchange is modeled by local source terms and results in coupled systems of mixed-dimensional partial differential equations. A motivating example is the simulation of root water uptake from soil when considering complex three-dimensional root network architectures explicitly [1], which allows to predict the complex water distribution in the soil around roots.
Consider a domain Ω ⊂ R 3 with an embedded tubular network. The centerline is a network of curves connected at branching points and is denoted by Λ, as shown in Fig. 1. In this work, we want to discuss the stationary mixeddimensional nonlinear equation systems of the form where s is a local coordinate on Λ that has a unique mapping to some position x ∈ R 3 . The operators ∇· and ∇ are the spatial divergence and gradient operators in Ω, and ∂ s denotes the derivative in direction of s. The scalar unknowns in Ω and Λ are denoted byû b (b for "bulk") and u e (e for "embeddded"). The diffusion coefficients D b > 0 and D e > 0 are positive and continuous and D b is a (possibly nonlinear) function of u b . In the source term, γ is a diffusive permeability and P = P (s) is the cross-sectional tube perimeter at s. Moreover, Π P • is an average operator such that Π Pûb denotes the average ofû b on the perimeter P for a given s, (1.4) where ζ is some suitable parameterization of P in R 3 . Furthermore, we introducê This mixed-dimensional model based on line sources, Eqs. (1.1) to (1.3), can be derived from a corresponding model with the three-dimensional tubular network structure cut out from Ω, by assuming that the tube radii R are much smaller than average distance between tubes [2].
The water distribution around a three-dimensional root network taking up water from the embedding soil can be modeled by Eq. (1.1) [3,4,5,1]. In this case, the unknowns are hydraulic pressures in root and soil, and D b corresponds to the hydraulic conductivity. Soil can be viewed as a three-phasic porous medium composed of the solid matrix and two fluid phases, air and water. With decreasing water content the soil's hydraulic conductivity decreases drastically, and at the same time, capillary forces increase which attract water to the hydrophilic solid matrix [6]. Under the assumptions of local mechanical equilibrium a direct and nonlinear relationship between the local water pressure and the hydraulic conductivity can be found [7,8]. At low water saturation, high pressure gradients are necessary to move water in dry soil. According to the cohesion-tension theory [9,10], transpiration at the plant leaves causes a high suction potential in plant roots and plants can maintain water uptake even in relatively dry soils. The root water uptake rate (q) is proportional to the root-soil pressure difference [4]. This can cause large local pressure gradients around roots [10,11] which are difficult to approximate with standard numerical schemes. We will use root water uptake as the motivational application and in the numerical examples in this work.
Equations (1.1) to (1.3) also arises for heat conduction problems, for example when modeling the temperature distribution around geothermal well systems [12,13]. Here, D b corresponds to the thermal conductivity, which generally depends on temperature [14]. When modeling the flow field around a well in a confined aquifer under high injection rates the surrounding rock undergoes deformations. The hydraulic conductivity of the rock depends on the pore pressure.
The paper is structured as follows. Motivated by the results of [15], in Section 2 we propose a distributed source model to replace Eqs. (1.1) to (1.3). In Section 3, we design a numerical scheme to accurately approximate the interface unknowns and hence Eq. (1.3), locally for each tubular segment. The method is based on local analytical solution obtained by means of Kirchhoff transformation. In Section 4, we show and discuss numerical results. In particular, in Sections 4.3.1 and 4.3.2 we investigate the reconstruction scheme numerically for a series of carefully constructed verification scenarios. Finally, in Section 4.4, we simulate root water uptake with a realistic root network obtained from MRI measurement and compare our results against a numerical reference solution obtained with state-of-the-art methods.

The distributed source model
The introduced model formulation, Eqs. (1.1) to (1.3), leads to solutionsû b which exhibit singularities on Λ. It is therefore difficult to construct efficient and accurate numerical schemes for solving Eqs. (1.1) to (1.3) [16,17,15]. However, for a precise description of the source term coupling the network and bulk domain, it is crucial to accurately approximate the solution in the neighborhood of the network. Koch et al. [15] suggest to solve a modified problem where Φ Λ denotes a set of non-negative kernel functions Φ Λ,i that distribute q around a vessel segment i over a small radially-symmetric tubular support region with radius (s), where Φ Λ,i = 0 outside the support region (compact support), cf. Fig. 1.

Remark 1.
We choose kernel functions Φ Λ,i (s) along each segment i in the form [15] where the function ϕ is a positive symmetric mollifier [18,19] defined on a disc perpendicular to the vessel segment, cf. Fig. 1. An example for such kernel functions is given in Eq. (2.5). The kernel functions Φ Λ have dimension L −2 (L: length) and bridge the dimensional gap of 2 between Eq. (2.1) and Eq. (2.2). In this work we use a uniform distribution, For ≤ R,û b can be approximated by Π P u b . However, for > R, this approximation is poor as the distribution kernel leads to a locally regularized solution u b =û b (visualized in Fig. 2). Effect of the regularization kernel on the local soil pressure solution. Symbols: u b,0 , regularized unknown at the tube centerline;û b , (average) bulk interface quantity (e.g. soil pressure); ue, physical quantity inside the tube (e.g. root pressure); R, tube radius; , kernel support radius.
When simulating systems with large networks of thousands of tubes (e.g. root systems) resolving the local solution around each network segment requires fine local computational meshes. The kernel distributes the source or sink term in a local neighborhood around network segments. This leads to a smooth u b which is easy to approximate by standard numerical schemes, see Fig. 2 for schematic representation at the example of a root segment cross-section. Unfortunately, the value of u b at the tube-bulk interface does not correspond to the interface value of the line source model Eqs. (1.1) to (1.3) anymore, i.e. Π P u b =û b . In the following, we exploit the fact that we know how the introduced kernel function modifies the local solution around a single isolated tube segment. The design of a method to reconstructû b accurately from the smooth solution u b -in the presence of the nonlinearity introduced by the diffusion coefficient D b -is the main contribution of this work.
Finally, we note that the linear case, i.e. D b = const. is discussed in [15]. For the particular case of root water uptake, local corrections based on the analytical solutions of the Richards equation have been proposed in [20]. However, the scheme is only presented in the discrete setting for a single voxel. The authors of [11] suggest to solve local radial-symmetric problems at every root segment, introducing additional unknowns. Again, the method is only presented in the discrete setting. In this work, we follow [15] and present a model formulated on the continuous setting which allows generalization to any suitable discretization method and simplifies the analysis of possible sources of errors. Moreover, the source distribution kernel allows to control the accuracy of the interface reconstruction in interplay with the mesh size. Apart from getting rid of singularities, the problem formulation (2.1) to (2.3) has the advantage that it does not require that the computational grid resolves the length scale R to yield accurate ap-proximations of q [16], but relaxes this requirement to the grid being required to resolve the length scale , a selectable model parameter. It has been shown [15] that this formulation allows to significantly reduce the error in q, u e , and u b for coarse grids.

Local reconstruction of interface unknown and flux
In this section we describe a method to accurately reconstruct the interface unknownû b for a given tube segment from the evaluation of the smoothed solution u b on the centerline Λ. The Kirchhoff transformation, well-known from the solution of heat conduction problems [21,Eq. (10)] allows us to transform Eq. (2.1) such that the nonlinearity only appears in the source term. We then derive a local analytical solution for u b depending on tube and kernel radius. From this analytical solution, we deduct a nonlinear equation to computeû b from a point evaluation of u b . We conclude by discussing the validity of the approach for the case of tubular networks.

Kirchhoff transformation
Equations (2.1) to (2.3) is a nonlinear equation system, if the diffusion coefficient D b depends on u b . Let us introduce the following Kirchhoff transformation [22] T : and we defineψ analogously in terms ofû b . The chain rule yields which allows us to rewrite the left-hand-side of Eq. (2.1), in terms of the transformed variable ψ. Note that the left-hand-side is now a linear operator, and the nonlinearity is contained in the source term (3.4) in terms of the inverse Kirchhoff transformation. This is an essential step, since Eq. (3.3) can be solved analytically in a simple radially-symmetric setting, cf. Section 3.2. For example, for the exponential function D b = D 0 exp{k(u b − 1)} and k ∈ R + , a function that fulfills our preconditions on D b and that we will use subsequently in numerical verification tests, this is indeed the case. There, ψ c corresponds to u b = −∞ and is well-defined. Nevertheless, T −1 is typically ill-conditioned around ψ c (unbounded derivatives). This is evident in Fig. 3 for k = 5 where it becomes clear that the derivative can get arbitrarily large.  Proposition. Fortunately, this singularity disappears in the non-degenerate case [22] D b (u b ) ≥ D min for some D min > 0. (3.5) In this case both T and T −1 are defined on all of R. For instance, this property can be achieved for the exponential function D b by using D min as its lower bound, i.e. D b := max(D 0 exp{k(u b − 1)}, D min ).
The functions shown in Fig. 3 are regularized in this way using D min = 1 × 10 −6 (see k = 5).

Local cylinder model and interface reconstruction
Firstly, let us consider an infinitely long straight cylinder of radius R embedded in an infinite domain Ω ∞ ⊂ R 3 with u e being a given constant. Solving Eqs. (2.1) to (2.3) then reduces to finding radially symmetric solutions by solving on a cross-sectional plane in local cylinder coordinates. By using the Kirchhoff transformation, Eq. (3.1), and the definition ofû b , Eq. (1.5), we obtain As discussed above, q is now a nonlinear function ofψ (which solves Eq. (3.7) for → 0). Using the uniform distribution kernel, Eq. (2.5), which is only non-zero in the local support region defined by the kernel radius ≥ R, we obtain an analytical solution, Continuing in this setting and evaluating Eq. (3.8) at r = δ, where 0 ≤ δ < denotes a distance to some point close to the tube centerline (r = 0), yieldŝ where we introduced the symbol ψ δ := ψ(r = δ). We note that we evaluated ψ in the region regularized by the source distribution kernel. Moreover, recall that the coupling source term is given by (3.11) and with the Kirchhoff transformation and radial symmetry, we know that Result 1. If we can estimate u b,δ (for example by the discrete cell value in a finite volume discretization), we can find the interface pressureû b by solving a nonlinear equation composed of Eqs. (3.10) to (3.12), which then allows to compute the source term q, given by Eq. (3.11). This means that an accurate approximation of u b,δ and u e is sufficient to compute a good approximation of q(û b , u e ).
This also suggests that we can indeed solve Eqs.  4. In the case that D b is a constant, Eq. (3.13) can be explicitly solved forû b as shown in [15].

Remark 5.
To see that Eq. (3.13) has a unique solution, it is sufficient to show that the expression on the left side is monotone with respect toû b . In general, this only holds under some conditions. First, we assume that the kernel radius is chosen large enough such that ln R ≥ 0.5. With this assumption, the third term is a monotonically decreasing (linear) function ofû b . Due to the assumption of positive diffusion coefficients, the functional −T is also monotonically decreasing. It follows the monotonicity of the entire left hand side expression in Eq. (3.13) (constant terms do not influence the monotonicity), and Eq. (3.13) therefore has a unique solution.
For the general three-dimensional case, the second equality in Eq. (3.12) is only an approximation. The equality holds here due to the radial symmetry of the solution. This condition is violated in the presence of multiple arbitrarilyoriented tubes as they readily occur in tubular network structures (microvasculature, root systems, fibre networks). We motivate in the next section why this error is expected to be small.

Multiple interacting parallel tubes
Let us consider the case of many parallel tubes where the solution u b only varies linearly along the tubes. Hence, we can consider a two-dimensional crosssectional plane and denote with x ∈ R 2 a position on this plane. Due to the linearity of the Laplace operator in Eq. (3.3), we can obtain a general solution for ψ by superposition where H is some harmonic function (for example chosen such that some boundary conditions are satisfied), We remark thatû b,j in q j depends on contributions from all partial solutions ψ j . To simplify notation, we assume i ≤ R i in the following (û b,i = Π Piûb = Π Pi u b ). Applying the linear averaging operator (defined in Eq. (1.4)) on both sides of Eq. (3.14) and assuming that the tubes do not overlap yields where we used the mean value theorem for harmonic functions, i.e. Π Pi ψ j = ψ j (x i ). From this, it follows that yields the reconstruction equation which is equivalent to Eq. (3.13) with δ = 0. Theũ bi ≈û b,i resulting from solving Eq. (3.20) is an approximation due to Eq. (3.19). This extends our result from the single vessel case to multiple parallel vessels, but with the introduction of some approximation error due to the nonlinearity of D b . Since Eq. (3.20) only requires point evaluations of u b , the result extends to i > R i as long as the kernel support regions do not overlap.

Error estimate for the approximation of the average operator
In the following, we estimate the error associated with approximation (3.19) in the reconstruction of the interface unknown. Let us assume that ψ(x i ) and u e (x i ) are given and denote with u b the exact solution. Again, to simplify notation, we assume i ≤ R i . Subtracting Eq. (3.18) from Eq. (3.20) yields for some U ⊂ R such that u b (Ω) ⊂ U , and using Taylor's Theorem results in for any u ∈ U . Inserting the exact solution and applying the averaging operator on both sides give and inserting this expression into Eq. (3.21) results in (3.24) The above equation can be equivalently written as (3.25) By assuming that there exists some constantC i , independent of R i , such that If the solution u b is sufficiently smooth (e.g. C 1 or Lipschitz continuous), then we deduce from the inequality above that the error introduced by the . Furthermore, it also shows that there is no error if D b = const., i.e. Eq. (2.1) is a linear diffusion equation, or if u b is constant on P i (i.e. radial symmetric solution).
Remark 6. In a numerical scheme, the exact solution ψ is approximated by some discrete solution ψ h , for which it holds that |ψ( (h denotes the discretization length, see Section 4), when using a second order scheme. Therefore, the approximation This suggests that on grids where h > R i , Π Pi ψ ≈ T (û b,i ) yields a good approximation, without being the main source of error. This result is supported by the numerical results in Section 4. Since the goal of the kernel method is to allow coarser grids by choosing kernel support radii > R i while maintaining accuracy [15], grid resolutions with h > R i correspond to the typical use case.

Multiple arbitrarily-oriented tubes
The numerical method introduced above for single or parallel tubes can also be applied for the general three-dimensional case with arbitrarily-oriented tubes. However, for arbitrarily-oriented tubes an additional error is introduced because the mean value property of harmonic functions, used for ψ j to derive Eq. (3.18), is no longer valid.
The additional error depends on 16). Assuming that the kernel support regions are non-overlapping, we apply Taylor's Theorem to deduce the following estimate Assuming that a contribution of another tube has the shape of a line source, the constant C can be computed, yielding where B Ri (x i ) denotes a ball of radius R i centered at x i , and E ⊥ j orthogonally projects x onto the centerline of tube j. The derivation of the estimate is given in Section B. We note that most other tubes in a large network system embedded in Ω are far away from the segment i, and therefore the error is small. Close tubes may cause a signification error. However, average distances to the closest neighbor tube in relevant applications are often 10R i and larger. With these estimates and by using similar arguments as in Section 3.3.1, we conclude that also for the case of arbitrarily-oriented tubes, the error introduced by reconstructingũ b,i from Eq. (3.20) using approximation (3.19) is not dominant on grids where h > R i . This conclusion is supported by numerical experiments, e.g. [23].
Furthermore, the line segments in practical networks are finite and kinks and bifurcations may introduce additional errors since there the kernel support regions of two connected vessels may overlap and the assumption of discrete cylinders as segments may result in errors in the estimation of the network surface area. In [23], the author proposes a numerical method that does not suffer from these errors since the tube interface is explicitly resolved by the three-dimensional computational mesh. It is shown that the approximation of kinks and bifurcations with cylinder segments do not introduce significant errors and that these errors are irrelevant in practical simulations of tissue perfusion or root water uptake.

Numerical results and discussion
The three-dimensional bulk domain Ω and the network domain Λ are spatially decomposed into the meshes Ω h and Λ h consisting of control volumes (cells) K Ω ∈ Ω h and K Λ ∈ Λ h , respectively. The discretization length, i.e. the maximum cell diameter, is denoted as h. To discretize the nonlinear diffusion problem Eqs. (2.1) to (2.3) in space, a cell-centered finite volume method with a two-point flux approximation is employed [15]. The resulting nonlinear system of equations is solved with Newton's method. To solve the learnized system of equations within each Newton iteration, we use the same linear solver as described in [15], that is, a stabilized bi-conjugate gradient method with a block-diagonal incomplete LU-factorization. All presented methods and simulations are implemented using the open-source software framework DuMu x [24] with the network grid implementation dune-foamgrid [25] for representing the embedded network domain.
For general nonlinear constitutive models (e.g. the Van Genuchten curves in Section 4.4), an analytical expression for the Kirchhoff transformation Eq. (3.3) and its inverse cannot be found such that it has to be calculated numerically. Here, we use numerical integration based on the double exponential transformation [26] and Brent's method for the inverse transformation. For fast evaluation of the inverse transformation, the functional T −1 (ψ) is replaced by a lookup table with local linear interpolation and a high sampling rate. The nonlinear interface reconstruction, Eq. (3.13) (single tube) or Eq. (3.20) (multiple tubes), is solved using Brent's method.

Analytical solutions for multiple parallel tubes
In this section, we derive for the verification of the introduced method, two types of analytical solutions based on the superposition of point source solutions. We consider an infinite two-dimensional domain Ω ⊂ R 2 that cuts through N parallel non-overlapping circular tubes of different radii R i . The flow resistance that the tubes pose to flow through Ω and its volume are neglected so that Ω can be extended inside the tube radius and comprises the entire plane without circular cut-outs.
The error involved with this assumption has been numerically analyzed in [23] and found to be small. The assumption is also commonly used in 1d-3d models [27,28,2,29,30,16,31,15] based on the underlying assumption that tube radii are small. As noted in the beginning of this section, an analytical solution for ψ can be obtained by superposition, see Eq. (3.14). The analytical solution for u b is then found by numerical or exact inversion of the Kirchhoff transformation.
In the following, we compute the coefficients of such solutions numerically for given tube center positions x i , tube radii R i , and u e,i , fixed γ i , given D b (u b ), and i ≥ R i . For a fully determined analytical solution we require N average interface unknownsû b,i , and some constant C ψ (corresponds to the choice H ≡ C ψ in Eq. (3.14)). Furthermore, we have N equationŝ where x i,k ∈ R 2 , w i,k ∈ R + are K ip integration points and weights. We choose x i,k to be uniformly distributed on P i . To compute the solution numerically, Note that due to the steep gradients in the case k = 5 the minimum u b is found to be −6 × 10 4 . We therefore consider this case as a hard-to-solve test case. We note that by increasing the kernel radius , u b is much smoother, cf. Fig. 16. The interface reconstruction algorithm needs to reconstruct the average u b on the black line (tube-bulk interface) shown in this plot, regardless of .
we chooseû b,1 and then solve Eq. (4.1) with a Newton method, where in every step the dense linear system is solved, where u i is short notation forû b,i and are the nonlinear residuals. The partial derivatives are approximated by numerical differentiation. The resulting solution is denoted as U b . In a variation of the above algorithm we usẽ 4) and the resulting solution is denoted asŨ b . This variation exactly corresponds to the approximation (3.19). We will show in the subsequent numerical experiments that a distributed source scheme with an interface reconstruction based on Eq.

Discrete error measures
To quantify the discretization errors, we define the following relative discrete L 2 -errors for the unknown u b , its transformed variable ψ = T (u b ), and the source term q as where ψ ref is chosen as 0.1 (unless otherwise indicated), and where q KΛ and Q KΛ are the numerical and the exact source for the tube segment K Λ , defined as the integral of q in Eq. (2.1) over Finally, we analogously define relative discrete L 2 -errors with respect to the modified analytical solutionŨ (see Section 4.1), and accordingly denote them asẼ u b ,Ẽ ψ ,Ẽ q . We note that all errors are reported against the regularized solution u b rather than the line source solutionû b . The exact solutions for u e and q are equivalent in both settings.

Grid convergence tests
In the following, we present grid convergence tests against the analytical solution for a single tube and against the analytical solutions of Section 4.1 for multiple parallel tubes. The diffusion coefficient D b (u b ) is chosen as exponential function shown in Fig. 3 with D 0 = 0.5. Analytical expressions for the Kirchhoff transformation and its inverse are given in Section A. The diffusive wall permeability γ is chosen as 1.
As motivated in Section 4.1, we expect that the numerical solution converges to the modified solutionŨ b . However, U b andŨ b are expected to be very similar. They are identical for a single infinite tube, where the approximation (3.19) is exact. Let us consider Eq. (2.1) for an infinite straight tube with constant u e embedded in an infinite bulk domain. The problem reduces to solving Eq. (3.6). The analytical solution is given by Eq. (3.8). Note that for this radially symmetric case Eq. (3.19) is exact. We therefore only report errors with respect to U b since U b andŨ b (as described in Section 4.1) are identical for the single tube case. We choose R = 0.01, = 0.05,û b = 0.5, u e = 0.1, and numerically solve Eq. (3.6) on the unit interval with the analytical solution prescribed as Dirichlet boundary conditions at r = 1. The source term is computed based on u e using the proposed interface reconstruction, Eq. (3.13).

Single tube convergence rates
Grid convergence results are shown in Fig. 5. The grid is uniformly refined, starting from h = 20R. Initially, we observe a pre-asymptotic range since the kernel support is not resolved by the grid yet. After the 3rd refinement, where h = , all errors decay quadratically with uniform grid refinement for all quantities, the primary variable u b , the transformed variable ψ, and the numerical source term q. As observed by [15], the onset of second-order convergence is determined by the kernel radius rather than the tube radius. This allows for good control of the error even for simulations where fine grids are not feasible.

Multiple parallel tubes
Next, we consider a two-dimensional domain Ω = [−1, 1] × [−1, 1] that perpendicularly cuts three tubes with radii (R max , 3 4 R max , 1 2 R max ) centered at  Fig. 4 for i = R i . The same solution for i = 2R i is shown in Section C.
In the following, we investigate the influence of the exponential rate parameter k of D b (u b ) and the influences of the tube radii R max on the discretization and the error involved in approximation Eq. (3.19). We compute grid convergence both against U b andŨ b . The mesh Ω h is uniformly refined starting with 4 × 4 cells. On boundaries, we enforce the respective analytical solution as Dirichlet boundary condition for u b .
In the first case, we set k = 1 and R max = 0.2. The errors E u b , E ψ and E q are shown in Fig. 6 with uniform grid refinement. As motivated in Section 4.1,  we see convergence to the modified analytical solutionŨ b . Second-order convergence is observed for all relevant quantities. For the convergence test against U b , we observe a non-reducible error. Due to the way the analytical solution is constructed, we can identify this error as the model error caused by approximation (3.19). However, most interestingly, this error is very small (less than 0.1 % forẼ q ) in this case and shows only an influence in the convergence plot for grid discretization length below the tube radius, i.e. h < R max . In the second case, we investigate more challenging cases by varying k, cf. Fig. 3. The source error for different k is shown in Fig. 7. As convergence of all three fields u b , ψ, q are confirmed in Figs. 5 and 6, we only report errors for q in the following, allowing for a more concise presentation.
For k = 0 the exponential function reduces to a constant. Therefore, this case corresponds to the linear stationary diffusion equation analyzed in [15]. Due to the mean value theorem for harmonic functions, approximation (3.19) is exact and U b =Ũ b . The larger k the stronger the changes in D b . In particular, as shown in Fig. 4, the gradient of u b at the tube interface is large for larger k and the value of u b varies considerably along the tube perimeter. Nevertheless, the largest observed error for k = 5 is E q ≈ 1 %, and only dominates the discretization error for h < R max .
In a third case, we fix k = 1 and vary the tube radii R max = 0.2, 0.1, 0.05 and 0.02. The analytical solution is computed such that the source term of the largest tube is equivalent for all cases. The source errors are presented in Fig. 8. It is evident that the error due to approximation (3.19) decreases with smaller tube radius. This is in good agreement with the error estimate in Section 3.3.1. In particular, we can see again that the error is only relevant in comparison to the discretization error for h < R max .

Influence of the kernel radius
In the fourth test case, we investigate the influence of the kernel radius on the discretization error. Therefore, we choose a case where the approximation error due to Eq. (3.19) does not dominate the total model error. We choose R max = 0.05 and k = 1. The discretization length h is fixed at 0.125. All other parameters are the same as in Section 4.3.2. The kernel radius is increased from 2R i to 12R i . The resulting source errors E q with respect to the analytical solution U b are shown in Fig. 9. The error decays with increasing kernel radius.  Interestingly, the speed of this decay is comparable with the error decay by grid refinement and matches the observations in [15] for the linear diffusion equation. It can also be seen that as soon as the kernel regions of the largest two tubes start to overlap the source error increases again. The error decay is best explained by the better approximation of u b (x i ) used in the reconstruction algorithm, cf. Eq. (3.13) and Eq. (3.20). Since u b is increasingly regularized with increasing kernel support, it becomes easier to approximate the function numerically.
4.3.4. Discrete evaluations of u b,δ or u b,0 As noted in Section 3.2, the reconstruction allows to consider that the numerically measured quantity u b,KΩ in cell K Ω does not represent the value on the centerline but rather a value in some distance δ. This is particularly relevant for coarse discretizations where it may make a difference whether the root segment is located in the middle or the corner of a 3D cell.
As a fifth case, we present the results of the second case (where u b,KΩ is interpreted as u 0 ) in comparison with results for which we assumed in the reconstruction that u b,KΩ represents the u b,δ and δ is chosen as the mean minimum distance between tube segment and bulk cell [32], The resulting source errors are shown in Fig. 10. We observe that the error is  significantly reduced. As for the previous cases, the error curve flattens as soon as the error is dominated by the mean value approximation error.

Root-soil interaction scenario
In the following application scenario, we compute root water uptake with small root system architecture obtained from MRI measurements. The scenario is similar to benchmark scenario C1.2 presented in [33]. However, we solve a stationary problem for various root collar pressures enforced as Dirichlet boundary conditions at the root collar.  Figure 11: Root conductivities and radius for a lupin root system. Left and middle, age-dependent hydraulic root conductivities from [33]. The axial root conductivity Kax corresponds to De and the radial root conductivity Kr corresponds to γ in the nonlinear diffusion equation. Right, 8-day-old lupin root system reconstructed from MRI data (courtesy of M. Landl, FZ Jülich). Grid data available from [34]. The root segment radius is visualized to scale. The rooting depth is about 10 cm. Figure adapted from [35].
To make it easier for readers familiar with root-soil interaction, we introduce several symbols and relate them to the symbols in Eqs. (2.1) to (2.3). The soil pressure p s and root pressure p r (in Pa) correspond to the unknowns u b and u e . The term µ −1 Kk r corresponds to D b , where µ = 1 × 10 −3 Pa s is the viscosity of water, K (in m 2 ) is the intrinsic permeability of the solid matrix, and k r is the dimensionless relative permeability. Relative permeability is commonly modeled as a nonlinear function of water content which in turn can be described by a nonlinear function of soil water pressure. The axial root conductivity K ax (in m 4 Pa −1 s −1 ) corresponds to D e , and the radial root conductivity K r (in m Pa −1 s −1 ) corresponds to γ. We can then reformulate Eqs. (2.1) to (2.3) to obtain a stationary root-water uptake model neglecting gravity, wherep s denotes the average soil pressure at the root-soil interface. Let θ = S w θ s denote the water content, where S w is the water saturation and θ s is the water content at saturation (equal to the porosity of the soil), and where p c (in Pa) is called capillary pressure, S e = θ r θ −1 s and θ r , θ s , α, n, m = 1 − n −1 are material-dependent parameters. In the following, we use a parametrization corresponding to loam given in [33], see Fig. 12. The functions in Eq. (4.12) are plotted in Fig. 12. The axial and radial root conductivities vary along the roots dependent on the root age. These root conductivity values are plotted in Fig. 11. For tabularized values, we refer to [33].
The root system shown in Fig. 11 is embedded in a box-shaped domain with dimension 8 × 8 × 15 cm. The top of the box intersects with the root collar at x 3 = 0 cm. The bottom of the domain is located at x 3 = −15 cm. We prescribe a water saturation of S w = 0.4 (corresponding to p s = 0.78 × 10 5 Pa) at all sides except for the top boundary where we enforce a zero-flow Neumann boundary condition. In the root domain, we prescribe no-flow boundary conditions at root tips and a fixed pressure p r,c at the root collar. We solve the same scenario for p r,c = 0.0, −0.5 × 10 5 , −1.0 × 10 5 , −2.5 × 10 5 , and −5.0 × 10 5 Pa. With decreasing root pressure, the flow rate of water leaving the domain at the root collar (transpiration rate) increases and the root-soil interface dries out. Dry soil (low water saturation) corresponds to a strong decrease of the local hydraulic conductivity, cf. Fig. 12.
We compare the results obtained with the presented distribution kernelbased method (DS) for i = 3R i and δ chosen as in Section 4.3.4 with two previously published methods. In [23] the root-soil interface is fully resolved by a locally refined unstructured three-dimensional mesh. The roots are modeled with Eq. (4.10) on a network of line segments. The solution p r (or u e ) is projected onto the closest surface on the three-dimensional mesh to evaluate a source term similar to Eq. (4.11). Since the root-soil interface is resolved, this resolved-interface method (PS) does not suffer from the approximation errors described in the previous sections, however, this comes with a higher computational cost as we will demonstrate. Secondly, we compare with the method presented in [16] and adapted for root-soil interaction as described in [32]. All methods are implemented in DuMu x [24,32] and are therefore easily comparable. In fact they share most of the source code and mostly differ in the way the soil and root domain are coupled. The simulation result for p r,c = −1 × 10 5 Pa and both methods PS and DS is shown in Fig. 13. A close-up shows the locally refined grid necessary to accurately resolve the root-soil interface. Due to the regularizing effect of both the distribution kernel and the coarse grid the DS solution u b does not contain the low saturation values found on the interface in the PS solution. Hence, the question is if these values can be accurately reconstructed from u b . As a global measure of how accurate the source terms q are approximated, we compute the transpiration rate at the root collar. Due to mass conservation, the transpiration rate is given by (4.13)   Transpiration rates for all three methods and all root collar pressures are shown in Fig. 14. Firstly, it can be observed that the DS method approximates all transpiration rates with a maximum relative difference of 3 % to the high fidelity PS solution, notably using a quite coarse grid (h = 5 mm). Using the same grid resolution the CSS method shows a difference of 10 % for p r,c = 0 Pa and even 60 % for p r,c = −5 × 10 5 Pa. Since CSS is consistent the results improve with grid refinement. The difference to the reference solution is comparable to that of DS in terms of the transpiration rate using strong local grid refinement resulting in a mesh of 1.9M cells. Although the finest mesh used for the PS method in this work is locally refined with a total of 12M cells and well-resolved root-soil interface (cf. Fig. 13), the grid convergence results shown in the rightmost plot in Fig. 14 suggest that the transpiration rate is not fully converged yet. Interestingly, following the trend, the transpiration rate is expected to get even closer to the solution of the DS method which is already stable with grid refinement for rather coarse grids. Figure 15 shows the reconstructed interface soil pressuresp s for all root segments in the mesh over the domain depth for the case p r,c = −5 × 10 5 Pa. For PS,p s is computed as a numerical integral over all coupling surface facets; for CSS as numerical integral over the cylinder surface; for DS it is the reconstructed value using the proposed reconstruction algorithm Eq. (3.20). The PS solution is considered a reference, although the results in Fig. 14 suggest that the result is not fully converged, the results can be considered quite accurate with less than 3 % difference in the transpiration rate to the DS method and the CSS method with strong local grid refinement. It can be seen that for the CSS method and the coarse grid, the interface soil pressure is overestimated explaining the large overestimation of the transpiration rate. With local grid refinement there is a much closer match with the PS solution. However, we note that in particular for the larger roots (lowp s ), difference between PS and CSS are still clearly visible. On the other hand, the DS closely matches the PS solution even for a coarse grid discretization. The largest difference is observed towards the root tips (highest p s ). The difference is improved by grid refinement. Note that at h = 2.5 mm and given that i = 3R i the kernel support around most smaller root segments is hardly resolved by one mesh cell, cf. Fig. 11. Therefore, this close match between the novel DS method and the PS method using a fully-resolved rootsoil interface and more than a thousand times more mesh cells is remarkable. This shows that the root water uptake problem, in particular in drier soils is clearly dominated by large and very localized gradients at the root-soil interface that are difficult to approximate for standard numerical schemes.
Finally, we remark that there is a strong reduction in computational time associated with the advantage of using a coarse grid discretization and a structured cube grid for the three-dimensional domain. While the DS method used 1 s (wall-clock time) per simulation for the 8k grid (5 min for 4M cells), the interface resolving PS method used 2.7 h for 11M cells. The DS method also required on average less Newton iterations. This can be attributed to the fact that the nonlinearity in D b is shifted into the reconstruction algorithm, whereas the numerical solution is regularized. For the projection method the large interface gradients have to be approximated by the discrete solution.

Summary and final remarks
Mixed-dimensional methods are efficient methods for solving coupled mixeddimensional PDEs arising from flow and transport processes in systems with tubular network systems embedded in a bulk domain. The bulk is represented by a three-dimensional mesh, the tubes are given as a network of cylinder center-line segments, and the meshes are typically non-matching. If the diffusion coefficient in the bulk domain depends on the unknown concentration a coupled nonlinear diffusion equation has to be solved. Its solution may exhibit large local gradients at the tube-bulk interface. This is for example the case when modeling water transport in soils with embedded root systems. Roots take up water from the soil and transport it upwards toward the atmosphere. In particular in dry soils, water uptake causes a strong local drop in the hydraulic conductivity leading to large pressure gradients around root segments.
We introduced an efficient numerical method for the solution of such nonlinear mixed-dimensional PDEs. The method is based on source distribution in a finite neighborhood region of the network, in combination with a nonlinear reconstruction scheme for interface unknowns. The method is based on several approximations we have analyzed. We estimated the errors associated with the approximations and showed in series of numerical verification tests that these errors remain small in practical applications. We use the new method to simulate root water uptake using a realistic root network. In comparison with existing methods we showed that the novel method outperforms other methods in both accuracy and efficiency. While the numerical results clearly show that the method accurately solves stationary problems, time-dependent problems are yet to be investigated in future work. However, preliminary results with root water uptake and slowly varying conditions (e.g. diurnal cycles) suggest that the method remains accurate.
Due to the possibility to use coarse computational grids in comparison with the tube diameter of the embedded networks, the presented method allows to perform simulations with large networks with reasonable 3D grid resolutions. Coarse grid discretizations (h R) in mixed-dimensional methods effectively lead to a distribution of any exchange source term q of an embedded tube into a neighborhood of diameter h. However, since local variations of a bulk unknown u b cannot be resolved, the approximation of a source term depending on u b evaluated on the tube-bulk interface suffers from significant errors. By deliberately introducing a distribution kernel for the source term, u b is locally modified and deviates from the real solution. However, by construction, the behavior of u b in the vicinity of the tube-bulk interface is better understood and it is possible to develop an interface reconstruction scheme, see Eq. (3.13) and Section 3.2. We have shown that the reconstruction significantly reduces the discretization error. Another approach is the Peaceman well model known in reservoir engineering [27,36,37]. Peaceman devises a reconstruction method eliminating discretization errors for one specific discretization scheme, and with some assumptions on the structure of the mesh and the orientation of the well tube. Methods based on the discrete formulation have also been explored for root water uptake simulations [20,38,11]. In contrast, the distributed source approach used in this work is formulated in the continuous setting which allows the analysis of the problem from a different perspective, and is applicable for any discretization scheme. Furthermore, the model remains valid in the case that the discretization length is smaller than the tube radius, which may readily occur if the tube radii vary significantly in the network (e.g large root systems). The distribution kernel effectively relaxes the strong discretization length restriction on the tube radius present in methods that require a direct numerical approximation of the interface unknown [15]. Our observations extend the results of [15] to the nonlinear case and indicate that as soon as the discretization length is in the range of the kernel radius (h ≈ ) the error in the source term is already reasonably small for many application scenarios. In effect, this property allows to improve the solution accuracy for simulations where finer grid discretization are not feasible, e.g. because of high computational costs.
where D 0 is a constant diffusion coefficient, u c = 1 + 1 k ln{ Dmin D0 } such that D b is continuous, and D min = D 0 (where > 0 is a small constant) is a minimal diffusion coefficient, for the purpose of rendering the Kirchhoff transformation, T (u b ), and its inverse, T −1 (ψ), well-defined on all of R (see Section 3.1).
We obtain analytical expressions for T (u b ) and T −1 (ψ) by splitting the integral range depending on the sign of u c and depending which one of u b and u c is larger. The Kirchhoff transformation, as defined in Eq. (3.1), is given for u c ≤ 0 by and for u c > 0 by where D b,r (u b ) = exp{k(u b − 1)}. The inverse transformation, for the case that u c ≤ 0, is given by The functions D b (u b ), T (u b ), and T −1 (ψ) are plotted for D 0 = 1 and = 1 × 10 −6 in Fig. 3.

B. Error estimate for arbitrarily-oriented tubes in 3D
In the following we want to derive an estimate for |Π Pi f j − f j (x i )|. Using the multivariate Taylor expansion for f j for any x ∈ B Ri (x i ) yields with the second-order residual |R 2 (x − x i )| ≤ C 2 x − x i We then assume that the contributions from neighboring tubes contributions from all neighboring tubes j can be formulated in terms of line sources and that the tubes are non-overlapping. The influence of segment sources decays faster with the distance than for the line source, so we consider this a conservative assumption. The corresponding functions f j are given by where E ⊥ j orthogonally projects x onto the centerline of tube j. For such functions (assuming non-overlapping kernel support regions), the constant C can explicitly be estimated by calculating the second-order derivatives. We note that gradients of the assumed f j are aligned with the radial direction of a cylinder coordinate system (r, θ, s) implied by centerline j. Therefore, it is enough to analyse radial derivatives of f j : , With the definition of C, this results in the estimate where for the last equality we have used the assumption that E ⊥ j (x i ) ∈ B Ri (x i ).
C. Analytical solution for three tubes ( i = 2R i ) In Fig. 4, we show an analytical solution for three point sources plotted for i = R i to show the variations ofû b over the tube-bulk interface that determine the approximation error of the interface reconstruction scheme. However, the numerical tests are run for the analytical solution with i = 2R i which suffer from the same approximation error of the interface reconstruction scheme but u b is much smoother due to the distribution kernel. To show how much influence this actually has on very nonlinear D b we here show in Fig. 16 the solution for i = 2R i . As noted in the introduction, the distribution kernel turns a difficult to approximate functionû b into a much smoother u b better suited for numerical approximation. However, using the developed reconstruction algorithm, the correct source term q can be reconstructed from the regularized solution.