text:
\section{Introduction}
\label{sec:intro}

\begin{figure}[htb!]
  \includegraphics[width=1.0\textwidth]{figure_1.pdf}%
  \caption{\textbf{Examples of systems with embedded network structures.} Left, blood vessel network simulation (geometry and boundary conditions from~\citep{Reichold2009})
           in the rat brain cortex coupled with the embedding tissue. Color signifies simulated blood pressure.
           A grey matter tissue volume of \SI{1}{\cubic\milli\m} may contain around 10'000 blood vessels.
           Right, simulation of root water and tracer uptake from soil by a young lupine root system,
           image adapted from~\citep{Koch2019cdumux} (license: CC BY 4.0).}
  \label{fig:introduction}
\end{figure}

There is a strong demand for efficient and accurate models describing flow and transport
processes in porous media with embedded tubular network systems, such as vascularized biological tissue,
plant root system growing in soil, hydrological, geothermal or petroleum wells in rock formations.
Reduced models are necessary due to the computational complexity arising from the large number
of network segments (for example about 10'000 blood vessels in a \SI{1}{\cubic\milli\m}
cube of gray matter brain tissue~\citep{Blinder2013}, or hundreds of meters of cumulative root length
in a $60$-day-old maize root system~\citep{Leitner2014b}) and the small diameter of the tubes with respect to
the entire computational domain
(for example wells of \SI{10}{\centi\m} diameter in a \si{\kilo\m}-scale reservoir).
Two motivational examples of mixed-dimensional simulations are shown in \cref{fig:introduction}.

Various methods have been developed recently to numerically solve coupled mixed-dimensional
partial differential equations (PDEs) that arise from flow and transport models in such systems.
Typically, flow and transport in the embedded tubular network system are described by one-dimensional equations
posed on a network of (center-)line segments. These networks are embedded into the surrounding bulk medium, often
porous media, which are described by three-dimensional equations. Network and bulk PDEs are coupled by
source terms that depend on state variables from both domains. The different numerical techniques
differ in the way they deal with the dimensional gap of $2$ between network and bulk domain.
The source term contribution in the embedding bulk medium can
be described by line source terms~\citep{Hsu1989green,dangelo2008,DAngelo2012,cattaneo2014computational,Gjerde2018},
surface source terms~\citep{koeppl2018,Laurino2019} or volume source terms~\citep{Koch2019a}.

A common assumption of mixed-dimensional models is that the radial scale of network tubes $R$
is much smaller than the dimension $\meas{\Omega}$ of the domain of interest $\Omega$. More specifically
for network systems, $R$ has to be much smaller than the average distance to the closest neighboring segment in the network.
While this precondition may be clearly satisfied in some cases (e.g. simple injection and extraction wells in large
distance to each other), it is less clear in others (blood capillaries with \SIrange{3}{8}{\micro\m}
vessel radius with average distances of \SI{50}{\micro\m}
in a microvascular network occupying about \SIrange{2}{3}{\percent} of the tissue volume).
Based on this main assumption, a couple of arguments are derived leading to the simplification of the model
equations. To allow for simple meshing procedure independent of the network domain,
the three-dimensional domain is extended to also cover the
volume occupied by the tubular network geometry~\citep{Laurino2019}.
This means the volume of the three-dimensional
domain is overestimated. Moreover, the tubes do not pose any resistance to flow since the physical
presence of the tubes is removed. Finally, since these mixed-dimension models
are usually derived for infinite cylindrical segments, some error is involved at
bifurcations by assuming finite cylindrical segments.
In this work, we want to investigate some of these assumptions numerically in more detail
than previously presented.

The mentioned mixed-dimensional models have been derived in the context of linear elliptic PDEs with line source terms.
Extensions to time-dependent problems and nonlinear problems have been analyzed to a much less extent.
Mixed-dimension methods are for example used to describe root water uptake from soil, in which case the
soil is described by the Richards equation---a strongly nonlinear PDE.
Although simple mixed-dimension methods have been used for more than two decades root water uptake simulations~\citep{Doussan1998},
the proper grid resolution required to accurately solve the model equations in dry soils
is rarely considered in the literature~\citep{Schroeder2009grid,Beudez2013}.
It is known that a coarse grid resolution in the soil domain may not accurately approximate local pressure gradients
and models for the root water uptake flux correction have been developed~\citep{Mai2019,Schroeder2008localsoil,Schroeder2009implmicro}.
To the best of our knowledge, the coupled root water uptake problem has not been rigorously analyzed mathematically.
The estimation of discretization errors and possible errors in the model reduction are yet to be better understood.
In this work we present a model technique to investigate such errors numerically.

In previous works on tissue perfusion models~\citep{dangelo2008}
and in more general mathematical works~\citep{koeppl2016},
it has been found that mixed-dimension embedded schemes
exhibit sub-optimal convergence rates, if the local discretization length $h$ is
larger than the radius $R$ of the embedded vessels.
While convergence rates alone are inconclusive about the error
at a given practical discretization length, the results in~\citep{koeppl2016,koeppl2018,Koch2019a}
indicate that in order to achieve sufficiently accurate numerical results,
the discretization length in the embedding bulk domain has to be chosen in the order of
the network tube radii or smaller. This is in stark contrast to typical grid resolutions
in root water uptake simulations, where soil cells
are routinely chosen an order of magnitude larger than the root radius~\citep{Schroeder2009grid,Leitner2014b}.
Several techniques to relax this
discretization length restriction have been discussed in the context of linear stationary
elliptic mixed-dimensional equations~\citep{Gjerde2018,Koch2019a,Koch2019bwell}.

In this work, we present a method where the tubular network is described by a network of line segments
with a given radius function as common in mixed-dimensional models. However, we explicitly resolve the
interface of the tubular network with the computational mesh describing the embedding bulk domain.
The new interface-resolving method developed subsequently can be considered
a reference method for the comparison of efficient mixed-dimension embedded schemes
based on implicit or reduced interface concepts.

For reference, we mention that in root-soil interaction simulations,
the root-soil interface has been explicitly resolved based on imaging data in a recent work by~\citet{Daly2018}.
However, only flow in the soil is simulated and the flow field is not coupled
to the flow field in the root xylem. Consequently, it is necessary to specify boundary conditions on the root-soil interface.
For the subsequently introduced method, the state of the root-soil interface is part of the solution.
Finally, we briefly introduced the new method in \citep{Schnepf2019benchmark}, where it is suggested
for the purpose of providing a reference solution in a benchmark study for root water uptake simulators.
In this work, we describe and analyse the method in more detail.

The new numerical method is derived in \cref{sec:projintro} (mathematical model)
and \cref{sec:disc} (discretization aspects) and
then applied in several numerical cases in \cref{sec:numeric}.
A grid convergence study in \cref{sec:onepbenchmark} shows
that the method is more accurate than other mixed-dimension methods for similar mesh sizes.
We compare the new method with previously published methods for examples for numerical test cases
of tissue perfusion and root water uptake in \cref{sec:tissue,sec:root}.

\section{Mixed-dimension method with resolved interface}
\label{sec:projintro}

\begin{figure}[htb!]
	\centering
	\includegraphics[width=0.6\textwidth]{figure_2.pdf}
	\caption{\textbf{Embedded tubular network system.} A tubular network structure
  with centerline skeleton $\Lambda$
	is embedded into the bulk domain $\Omega$. The surface of the tubes $\Gamma$ (dark grey) forms an internal boundary
  for $\Omega$, $\Gamma \subset \Omega \cap \mho$. Both domains are equipped
	with (local) coordinate systems. $R = R(s)$ denotes the equivalent local radius of the tube.}
	\label{fig:network}
\end{figure}

The tubular segments in network systems like plant roots or capillary blood vessels are usually
much smaller in radial extent than in axial extent, $R \ll L$. Often, it is therefore a good assumption
to neglect radial variations and work with cross-section averaged quantities and one-dimensional models
that describe the change of e.g. average pressure, temperature, concentration, etc. along the centerline axis~\citep{Koch2020PhDthesis}.
In this work we will assume that this one-dimensional description is sufficiently accurate and exploit this fact
by not resolving the network structure with a fully-resolved three-dimensional computational mesh which would
result in problems of intractable size. Moreover, we assume that any membrane separating the internal highly conducting space
of the tube (e.g. blood vessel lumen, root xylem) from the bulk domain can be described as a two-dimensional
sharp interface $\Gamma$.

In the following, we will exemplarily consider the case of root water
uptake.
For details on the mathematical modeling of root water uptake with three-dimensional root architectures,
we refer to the literature~\citep{Doussan1998,Javaux2008,Roose2008,Dunbabin2013,Koch2018a,Schnepf2019benchmark}.
Fluid flow in the root xylem---a structure that can be imagined as a bundle of tubes located in the center
of the root and transporting fluid in axial direction upwards toward the plant leaves---can be described by
\begin{equation}
\label{eq:xylem}
  - {\partial_s}\left( K_\text{ax} {\partial_s} p_{\up{1D}} \right) = -q \quad \text{on} \quad \Lambda,
\end{equation}
with some boundary conditions on $\partial \Lambda$,
where $s$ denotes the local axial coordinate, $K_\text{ax}$ (in \si{\m\tothe{4}\per\pascal\per\s}) is
the axial root xylem conductivity, $p_{\up{1D}}$ is the root xylem pressure (in \si{\pascal}), and $q$ (in \si{\square\m\per\s})
is a source term modeling fluid exchange with the embedding soil domain and depends on both $p_{\up{1D}}$
and the soil pressure $p_{\up{3D}}$ on the interface $\Gamma$.

Water flow in the soil is described by the Richards equation,
\begin{align}
\label{eq:richards}
	-\div \left(\frac{k_r(p_{\up{3D}})}{\mu} K \nabla\! p_{\up{3D}}  + \rho\vec{g} \right)  &= 0 & \text{in} \quad \Omega, \\
\label{eq:coupling}
  -\left(\frac{k_r(p_{\up{3D}})}{\mu} K \nabla\! p_{\up{3D}} + \rho\vec{g}\right) \cdot \vec{n}_\Gamma
  &= K_\text{r} \left[ p_{\up{3D}}(\vec{x}_\Gamma) - p_{\up{1D}}(\Pi\vec{x}_\Gamma) \right] & \text{on} \quad \Gamma,
\end{align}
with suitable boundary conditions prescribed on $\partial\Omega\setminus\Gamma$.
In \cref{eq:richards}, $\mu$ is the dynamic fluid viscosity (in \si{\pascal\s}),
$\rho$ is the fluid density (in \si{\kg\per\cubic\m}),
$k_r$ denotes the dimensionless relative permeability,
and $K$ the intrinsic permeability of the soil (in \si{\square\m}).
The relative permeability is a nonlinear function of $p_{\up{3D}}$, e.g. modeled
by the well-known Van Genuchten-Mualem model~\citep{mualem1976,van1980closed}.
In \cref{eq:coupling}, $\vec{n}_\Gamma$
is an outward-pointing (with respect to $\Omega$) unit normal on $\Gamma$
and $K_r$ is the radial root conductivity (in \si{\m\per\pascal\per\s}).
Finally, $\Pi\bullet$ is a surjective projection operator that maps any point $\vec{x}_\Gamma$ on $\Gamma$
to a corresponding point $\hat{s}$ on $\Lambda$, given a parameterization of $\Lambda$ in terms of $\hat{s}$.

In order to obtain a mass conservative coupling scheme, we need to define the source term $q$ in \cref{eq:xylem}.
To this end, we denote with $\Lambda_\varsigma \subseteq \Lambda$ some compact subset of $\Lambda$ and with
$\Gamma_\varsigma = \lbrace \vec{x}_\Gamma \in \Gamma \mid \Pi \vec{x}_\Gamma \in \Lambda_\varsigma \rbrace$
the corresponding set of surface points on $\Gamma$. Then, the coupling condition is given by
\begin{equation}
\begin{split}
  \int_{\Lambda_\varsigma}\! q \,\text{d}s
    &= \int_{\Gamma_\varsigma}\! \left(\frac{k_r(p_{\up{3D}})}{\mu} K \nabla\! p_{\up{3D}} + \rho\vec{g}\right) \cdot \vec{n}_\Gamma \,\text{d}\gamma \\
    &= \int_{\Gamma_\varsigma}\! -K_\text{r} \left[ p_{\up{3D}}(\vec{x}_\Gamma) - p_{\up{1D}}(\Pi\vec{x}_\Gamma) \right] \,\text{d}\gamma,
\end{split}
\end{equation}
given some suitable parameterization of $\Gamma_\varsigma$ in terms of $\gamma$.

\subsection{Practical geometry parameterization}
\label{sec:parametr}

\begin{figure}[!hbtp]
  \includegraphics[width=1.0\textwidth]{figure_3.pdf}%
  \caption{\textbf{Three-dimensional representations of a segmented root architecture (lupine).} From right to left: graph representation with segment-wise continuous radius field,
  implicit geometrical representation as the sum of capsules, and discrete representation of root-soil interface
  as the surface facets of a tetrahedron mesh of the soil domain (generated with CGAL~\citep{CGAL2019}).}
  \label{fig:root_projection}
\end{figure}

The interface $\Gamma$ between network and bulk domain may often be given by some implicit description
in form of a continuous or discrete level set function, e.g. obtained from imaging data. From such
an implicit description it is possible to generate surface triangulations \citep{Boissonnat2005}, and extract
center-lines, for example based on the medial axis transformation~\citep{Lindquist1996,Antiga2008}.
Nevertheless, there is in general no unique choice for the mapping $\Pi$.

Since our method is targeted at the creation of verification tests for reduced methods
without explicit interface resolution, we simplify the geometrical description as follows.
The root network center-lines $\Lambda$ are approximated by $\Lambda_h$, a set of
linear root segment center-lines $\Lambda_i$ defined by two points $\vec{q}_i$, $\vec{p}_i$ and parametrized by
\begin{equation}
  \label{eq:line_param}
  \tilde{\vec{x}}(\tilde{s}_i) := \vec{p}_i + \tilde{s}_i \vec{m}_i, \quad \vec{m}_i=\vec{q}_i-\vec{p}_i, \quad \text{for}\quad \tilde{s}_i\in[0,1].
\end{equation}
Moreover, associated which each segment is a continuous radius function $R_i(\tilde{s}_i)$
which is often---but not necessarily---constant per segment but varies from segment to segment.
From this representation, we implicitly define a three-dimensional network representation
by the signed distance functions (SDFs),
\begin{equation}
\begin{split}
  \label{eq:implicitdomains}
  d_{\Lambda_h}(\vec{x}) &:= \min\limits_{\Lambda_i \in \Lambda_h} d_{\Lambda_i}(\vec{x}), \quad \vec{x} \in \mathbb{R}^3,\\
  d_{\Lambda_i}(\vec{x}) &:= \norm{\tilde{\vec{x}}(\mathbb{P}_i(\vec{x})) - \vec{x}} -R_i(\mathbb{P}_i(\vec{x})), \\
  & \text{where} \quad \hat{s}_i := \mathbb{P}_i(\vec{x}) = \max\left\lbrace 0, \min\left\lbrace\frac{(\vec{x} - \vec{p}_i)\cdot \vec{m}_i}{\lVert \vec{m}_i \rVert^2_2}, 1 \right\rbrace\right\rbrace.
\end{split}
\end{equation}
This parameterization allows for the convenient definition of the operator $\Pi$ to yield
the position $\hat{s}_i$ on the segment $\Lambda_i$ with minimal $d_{\Lambda_i}(\vec{x})$.
\begin{remark}
The SDFs $d_{\Lambda_i}$ describe capsules with radius function $R_i(\tilde{s}_i)$ around the line $\Lambda_i$
and the SDF $d_{\Lambda_h}$ describes the union of all such capsules.
Then every point $\vec{x}$ is inside the root, if $d_{\Lambda_h}(\vec{x}) < 0$, or in the soil or outside the domain, if $d_{\Lambda_h}(\vec{x}) > 0$.
Consequently, the root-soil surface is given by the zero level set $\Gamma = \lbrace \vec{x} \,\vert\, d_{\Lambda_h}(\vec{x}) = 0 \rbrace$.
\end{remark}

In this work, we use the meshing capabilities of the C++ geometry library CGAL~\citep{CGAL2019}
to generate computational grids for the bulk domain from such an implicit description.
These grid explicitly resolve the bulk-network interface.
An exemplary root network, the three-dimensional representation as a union of capsules, and a triangulated representation of $\Gamma$
is shown in~\cref{fig:root_projection}.

We note that the surface implied by the zero level set of \cref{eq:implicitdomains}
is only piecewise differentiable due to the possible discontinuity of $R$ between segments.
If necessary, the surface's smoothness can be improved by using a smooth minimum function for $d_{\Lambda_h}$, such as
\begin{equation}
\operatorname{smin}(a, b, k) := \min\{a, b\} - \frac{1}{6} h^3 k, \quad h = \frac{1}{k}\max\lbrace k - \vert a-b \vert, 0\rbrace,
\end{equation}
rendering the surface function twice differentiable ($C^2$)~\citep{smoothmin2013}.
The parameters $a$ and $b$ are signed distances
and $k > 0$ is a smoothing parameter (with units \si{\m})
which is to be chosen in the order of magnitude of the dimensions of the objects merged.
In the following, we do not perform such smoothing of the network surface,
and approximate the surface by the zero level set of the distance function~\cref{eq:implicitdomains}.

\subsection{Relation to implicit surface mixed-dimension methods}
\label{sec:implicit}
\begin{figure}[htb!]
  \includegraphics[width=1.0\textwidth]{figure_4.pdf}%
  \caption{\textbf{Comparison of two mixed-dimension discretization concepts.}
  A root system is embedded in a soil domain. Left, resolved-interface discretization. Coupling
  terms appear as boundary conditions for the 3D bulk domain. Right, implicit interface method with non-matching
  discretization. The bulk domain is extended to include the space occupied by the network (in its three-dimensional representation).
  Coupling terms for each segment appear as source terms restricted to a line~\citep{dangelo2008}, a local surface~\citep{koeppl2018}, or a local volume~\citep{Koch2019a},
  depending on the chosen method. This approach introduces additional model errors but implicit interface methods
  are expected to be computationally more efficient than the resolved-interface approach.}
  \label{fig:model_comparison}
\end{figure}

With the parameterization introduced in \cref{sec:parametr}, consider
a circular cross-section of radius $R_i$ of an infinite cylindrical tube.
Then, we can show that
\begin{align}
\label{eq:cylinder}
q(\tilde{s}_i) &= \int_0^{2\pi}\! -K_r \left( p_{\up{3D}}(\vec{x}) - p_{\up{1D}}(\Pi\vec{x}) \right) R_i\text{d}\theta \nonumber\\
        &= -2\pi R_i K_r \left( \frac{1}{2\pi}\int_0^{2\pi}\! p_{\up{3D}}(\vec{x}) \,\text{d}\theta - \frac{1}{2\pi} \int_0^{2\pi}\! p_{\up{1D}}(\Pi\vec{x}) \,\text{d}\theta \right) \nonumber\\
        &= -2\pi R_i K_r \left( \hat{p}^\bigcirc_{\up{3D}} - p_{\up{1D}}(\tilde{s}_i) \right),
\end{align}
where we used the fact that $p_{\up{1D}}$ is independent of $\theta$,
and define the average pressure on the perimeter as
\begin{equation}
  \hat{p}^\bigcirc_{\up{3D}} := \frac{1}{2\pi}\int_0^{2\pi}\! p_{\up{3D}}(\vec{x}) \,\text{d}\theta.
\end{equation}
Hence, if the network is approximated by sufficiently long discrete cylinder segments
and exchange is assumed to only occur over the lateral surface of the cylinders, the source term can be formulated
solely in terms of quantities on the cross-sectional plane. Furthermore, assuming
that the network domain does not pose any resistance to flow in the bulk domain
and its volume is negligible, $\Omega$ is extended to include
the network domain, $\Omega^\text{ex} = \Omega \cup \mho$, and we arrive at a reduced model that can be written
in terms of a delta distribution restricting the bulk source term onto $\Gamma$, cf.~\citep{koeppl2018},
\begin{align}
	\label{eq:nonlineardiffusion-root}
	-\div \left(\frac{k_r(p_{\up{3D}})}{\mu} K \nabla\! p_{\up{3D}} \right)  &= q\delta_\Gamma & \text{in} \quad \Omega^\text{ex}, \\
	\label{eq:1d-root}
	- {\partial_s}\left( K_\text{ax} {\partial_s} p_{\up{1D}} \right) &= -q & \text{on} \quad \Lambda, \\
	\label{eq:source-root}
	q &= -2\pi R K_\text{r} (\hat{p}_{\up{3D}}^\bigcirc - p_{\up{1D}} ), \\
  \int_{\Omega^\text{ex}}\! q \delta_\Gamma \,\text{d}\vec{x} &= \int_{\Gamma}\! \frac{q}{2\pi R} \,\text{d}\gamma = \int_\Lambda\! q \,\text{d}s.
\end{align}
In the following, we refer to such formulations as \textit{implicit surface} methods as the exchange between bulk domain and network
is entirely formulated in terms of source terms in $\Omega^\text{ex}$ instead of boundary conditions
on $\Gamma$ and an explicit resolution of the bulk-network
interface by the computational mesh of the bulk domain is not necessary anymore.
However, we note that the interface $\Gamma$ still appears implicitly
in form of the delta distribution. \Cref{fig:model_comparison} compares the meshes
used for resolved-interface descriptions and implicit interface description in non-matching discretization schemes.

In \citep{d2007multiscale}, the author suggests
to use $\delta_\Lambda$ instead of $\delta_\Gamma$ in \cref{eq:nonlineardiffusion-root}, i.e. the source
term is restricted to a line in $\Omega^\text{ex}$. This formulation results in pressure solutions with feature
singularities on $\Lambda$ and are difficult to approximate by numerical schemes. In \citep{Koch2019a}, the authors suggest to replace the surface source
term by a volume source term using a distribution kernel in combination with a local reconstruction scheme of the interface
pressure $\hat{p}_{\up{3D}}^\bigcirc$. This technique allows to decouple the discretization length from the tube radius,
but the local reconstruction scheme has only been investigated for linear problems so far (corresponding to a constant relative permeability $k_r$ in \cref{eq:nonlineardiffusion-root}).
In \cref{sec:numeric}, we compare the new resolved-interface method with implicit interface methods in numerical experiments.
To this end, we follow the terminology of \citep{Koch2019a} and refer with \textsc{css} (\textit{cylinder surface source}) to the method due to \citep{koeppl2018} using
formulation \cref{eq:nonlineardiffusion-root,eq:1d-root,eq:source-root}. We refer with \textsc{ls} (\textit{line source}) to the method due to \citep{d2007multiscale}
where $\delta_\Gamma$ is replaced by $\delta_\Lambda$, and with \textsc{ds} (\textit{distributed source}) to the method due to \citep{Koch2019a},
where $\delta_\Gamma$ is replaced by a volumetric distribution kernel and the source term $q$ is computed based on a local
reconstruction scheme.

\subsection{Integration of the coupling term}
\label{sec:disc}

In the discrete setting, the root and the soil domain, $\Lambda$ and $\Omega$
are partitioned into a finite number of grid cells such that
$\Lambda_h = \bigcup K_\Lambda$ and $\Omega_h = \bigcup K_\Omega$ are discrete mesh representations
of $\Lambda$ and $\Omega$ with the cells $K_\Lambda$ and $K_\Omega$.
The computational grids can be chosen independently.
A part of the boundary of $\Omega$ explicitly resolves the root-soil interface $\Gamma$ and
$\Gamma_h$ denotes the set of cell facets on the interface.
Since the interface is explicitly described by $\Gamma_h$, the
the coupling conditions, \cref{eq:coupling}, can be directly evaluated by
numerically approximating the surface integrals. However, in the discrete setting,
the approximation of $p_{\up{3D}}$ is typically only piecewise differentiable.
In this work, we consider piecewise linear functions. The
source term $q$ needs to be integrated over a control volume $K_\Lambda$ which may involve integration
over several interface facets. For this purpose, we suggest an algorithm based on virtual local refinement
of the interface facets to accurately capture the surface integration area element
associated with the integration over $K_\Lambda$. The algorithm is given as pseudo-code in \cref{algo}.
In brief, we map the corners of a surface triangle with $\Pi$ and evaluate if the mapped points
are contained in different network control volumes. If so, the triangle is virtually refined and
the procedure is repeated recursively until all corners map to the same control volumes, or some
maximum refinement level is reached. We add only one integration point per surface
triangle and coupled network control volume at the centroid of the union of coupled sub-triangles.
For $p_{\up{3D}}$ and $p_{\up{1D}}$ being piecewise linear functions,
integrating with the mid-point rule is exact.

\section{Numerical results and discussion}
\label{sec:numeric}

In this section, we compare the introduced explicit interface method
with previously suggested implicit interface methods in three cases.
We use the abbreviations \textsc{ls}, \textsc{css}, \textsc{ds} introduced
in \cref{sec:implicit} for the implicit interface methods and abbreviate with
\textsc{ps} (\textit{projection source}) the explicit interface method.
In the first case, \cref{sec:onepbenchmark}, we show for a simple rotation-symmetric setup with a single
tubular inclusion that \textsc{ps} accurately approximates a given
analytical solution and verify that the surface integration scheme
proposed in \cref{sec:disc} is sufficiently accurate.
In the second case, \cref{sec:tissue}, we investigate errors introduced by implicit interface
methods at the example of tissue perfusion described by a linear elliptic
mixed-dimensional model. In the third case, \cref{sec:root},
we investigate differences between \textsc{ps} and \textsc{css} in
a root water uptake example described by a nonlinear elliptic mixed-dimension model
based on the Richards equations.

The three-dimensional bulk domains $\Omega$, $\Omega^\text{ex}$ and the network domain $\Lambda$
are spatially decomposed into the meshes $\Omega_h$, $\Omega_h^\text{ex}$ and $\Lambda_h$ consisting
of cells $K_\Omega \in \Omega_h$ and $K_\Lambda \in \Lambda_h$, respectively.
The discretization length computed as the maximal cell diameter is denoted by $h$.
We discretize the continuous equation in space using finite volume methods.
For structured Cartesian grids as well as for the network equations, we
use a cell-centered finite volume method (\textsc{fvm}) with a two-point flux approximation (\textsc{tpfa}), cf.~\cite{Koch2019a}.
When using (unstructured) tetrahedral meshes for the \textsc{ps} method,
or in the case of locally refined meshes for the \textsc{css} method in \cref{sec:root},
we use vertex-centered
finite volumes with linear basis functions
(also referred to as \textsc{box} method)~\citep{Hackbusch1989,HuberHelmig1999,Koch2019cdumux}.
This is because cell-centered \textsc{tpfa}-\textsc{fvm} are generally not consistent
on such meshes~\citep{Schneider2018}.
The resulting discrete system of equations is solved with Newton's method. (In case of a linear
model Newton's method converges in one step.)
The linearized system of equations within each Newton iteration,
is solved with a stabilized bi-conjugate gradient method using a block-diagonal preconditioner
based on incomplete LU-factorization, cf.~\citep{Koch2019a}.
All presented methods and simulations are implemented using the open-source software framework DuMu\textsuperscript{x}\xspace~\cite{Koch2019cdumux}
with the network grid implementation \texttt{dune-foamgrid}~\cite{foamgrid} for representing the embedded network domain.

\subsection{Mixed-dimension single phase flow}
\label{sec:onepbenchmark}

\begin{figure}[htb!]
  \includegraphics[width=1.0\textwidth]{figure_5.pdf}%
  \caption{\textbf{Reference solution for an embedded cylinder.}
  Cylindrical tube embedded in a box-shaped domain with dimensions $[-1,1]\times[-1,1]\times[-1,1]$.
  The domain is cut in half revealing the tube of radius $R = 0.03$ visualized with a reduced radius
  to make the interface visible. The numerical solutions $p_{\up{3D}}^\textsc{ps}$,
  $p_{\up{1D}}^\textsc{ps}$ computed on a fine grid with $h = 0.008$ are shown exemplarily
  but they are visually identical to a plot of the analytical solution.}
  \label{fig:cylinder_bench_plot}
\end{figure}
\begin{figure}[hbtp]
  \includegraphics[width=1.0\textwidth]{figure_6.pdf}%
  \caption{Convergence rates over discretization length.}
  \label{fig:cylinder_bench}
\end{figure}

Let us consider a slightly simplified problem, adapted from \citet{d2007multiscale,Koch2019a},
(for simplicity we use the same symbols as previously introduced but
all unknowns and parameters are to be interpreted as dimensionless quantities,)
\begin{subequations}
\label{eq:kernel:problembeta}
\begin{align}
  -\partial_s \left(K_\text{ax} \partial_s p_\up{1D} \right) &= -q \quad&& \text{in } \Lambda,\\
  -\div(\nabla\! p_\up{3D}) &= q\delta_\Lambda \quad&& \text{in } \Omega^\text{ex}, \\
  q &= - 2\pi R K_r (\hat{p}_{\up{3D}}^\bigcirc - p_\up{1D}),
\end{align}
\end{subequations}
with the domains $\Omega^\text{ex} = [-1,1]\times[-1,1]\times[-1,1]$ and $\Lambda = \{0\}\times\{0\}\times[0,1]$,
i.e the vessel center-line coincides with the $x_3$-axis. The tube has radius $R$ and $\mho$
is given by the cylinder with center-line $\Lambda$, radius $R$ and unit length.
Recall that for this straight cylindrical tube case, due to the observation in \cref{eq:cylinder}, problem formulation \cref{eq:kernel:problembeta}
is equivalent to the formulation with boundary conditions on $\Gamma$, cf.~\cref{eq:xylem,eq:richards,eq:coupling} and
\begin{equation}
  \label{eq:cyl_benchmark_coupling}
  \int_\Lambda \! q \,\text{d}s = \int_{\Omega^\text{ex}} \! q\delta_\Lambda \,\text{d}x = \int_{\Gamma}\! -K_\text{r} \left[ p_{\up{3D}}(\vec{x}_\Gamma) - p_{\up{1D}}(\Pi\vec{x}_\Gamma) \right] \,\text{d}\gamma.
\end{equation}
Choosing the conductivities as
\begin{equation*}
K_\text{ax} = 1 + x_3 + \frac{1}{2}x_3^2, \quad K_r = \left(2\pi R + R\ln {R}\right)^{-1},
\end{equation*}
the pressure solutions,
\begin{equation}
\label{eq:analytic_benchmark}
p_{\up{1D},\up{e}} = 1 + x_3, \quad p_{\up{3D},\up{e}} = - \frac{1 + x_3}{2\pi} \ln r, \quad r \geq R
\end{equation}
with $r = \sqrt{x_1^2 + x_2^2}$, solve \cref{eq:kernel:problembeta} with matching boundary conditions.
From the analytical pressure solutions follows that $q_\up{e} = 1+x_3$ is the analytical source term.

In this setting, we can directly compare the explicit interface method with
implicit interface methods.
The analytical solutions $p_{\up{1D},\up{e}}$ and $q_\up{e}$ are identical for all mentioned methods, cf.~\citep{Koch2019a}.
The exact pressure in the bulk domain, $p_{\up{3D},\up{e}}$, can be extended to all of $\Omega^\text{ex}$ and differs between the different methods
for $r < R$ (and $r < \varrho$ for the distributed source method (\textsc{ds}) of \citep{Koch2019a} where $\varrho$ is the radius of the distribution kernel),
but is identical for $r \geq R$ ($r \geq \varrho$ for \textsc{ds}).
The analytical solutions $p_{\up{3D},\up{e}}^{\textsc{m}} \in \Omega^\text{ex}$ for the methods $M \in \lbrace \textsc{ls}, \textsc{css}, \textsc{ds} \rbrace$ are given in \cref{app:benchmark}.

For the resolved interface method \textsc{ps}, in the case of straight cylindrical vessels,
integration of the source term $q_{K_\Lambda}$ over the discrete interface $\Gamma_h$ can be performed exactly.
To this end, we compute for every segment $K_\Lambda$ the intersections
of the space between the two planes implied the segment $K_\Lambda$ (the two planes through the end points
and normal to the segment) and all surface triangles $T \in \Gamma_h$. Effectively,
every $T$ is sub-triangulated such that every sub-triangle couples with exactly one $K_\Lambda$.
The boundary integral can be computed by the mid-point rule which is exact since the integrand is a linear function.
Intersections cannot be so easily computed in networks where the interface is only given as an
implicit function. Therefore, we suggested an approximate integration algorithm based on virtual refinement
in \cref{sec:disc}. For this particular test case, we implemented both approaches to verify
the accuracy of the latter. We denote the exact approach with \textsc{ps-e}
and the approximate approach with \textsc{ps-a}.

We solve \cref{eq:kernel:problembeta} with the methods
\textsc{ls}, \textsc{css}, \textsc{ds}, \textsc{ps-e} and \textsc{ps-a} and by prescribing
the analytical solutions as Dirichlet boundary conditions, except for the top and bottom sections of $\Omega$ or $\Omega^{ex}$ ($x_3 = 0$, $x_3 = 1$)
where we prescribe the normal derivative of the analytical solution as Neumann boundary condition.
Obviously, for the methods \textsc{ps-e} and \textsc{ps-a}, $\Gamma$ does not require boundary conditions
and fluxes over $\Gamma$ are computed by the coupling conditions, \cref{eq:cyl_benchmark_coupling}.
The numerical solutions $p^{\textsc{ps}}_\up{3D}$ and $p_\up{1D}$ for ${R} = 0.03$
are shown in~\cref{fig:cylinder_bench_plot}.

We compute pressure discretization errors in the normalized discrete norm
\begin{equation}
 \lVert p_\up{3D} - p^{\textsc{m}}_{\up{3D},\up{e}} \rVert_2 := \frac{\left[ \sum_{\Omega_h} \vert K_\Omega \vert (p^{\textsc{m}}_{K_\Omega,\up{e}} - p_{K_\Omega})^2 \right]^{1/2}}{\sum_{\Omega_h} \vert K_\Omega \vert},
\end{equation}
where $p_{K_\Omega}$, $p_{K_\Omega,\up{e}}$ denote numerical and exact pressure evaluated at the center of a control volume
$K_\Omega$ and $\vert K_\Omega \vert$ its volume. The error for $p_\up{1D}$ in $\Lambda_h$ is computed analogously.
The error in the source term $q$ is computed as
\begin{equation}
 \lVert q - q_{e} \rVert_2 = \frac{\left[\sum_{\Lambda_h} \vert K_\Lambda \vert (q_{K_\Lambda,\up{e}} - q_{K_\Lambda})^2\right]^{1/2}}{\sum_{\Lambda_h} \vert K_\Lambda \vert},
\end{equation}
where
\begin{equation}
  \label{eq:source_discrete}
 q_{K_\Lambda,\up{e}} = \int_{K_\Lambda} \!q_\up{e}\,\text{d}s \quad \text{and} \quad q_{K_\Lambda} = \int_{K_\Lambda}\! q \,\text{d}s.
\end{equation}
The maximum control volume size, $h$, is given by the maximum cell diameter in both domains.
We choose $h = h_\Omega$ such that $h_\Omega \approx h_\Lambda$. Both domains are uniformly refined.
The mesh for \textsc{ps} is remeshed so that the discrete interface $\Gamma_h$ approaches the real interface
$\Gamma$ with grid refinement.
Pressure and source error norms with grid refinement are shown in~\cref{fig:cylinder_bench}.

For sufficiently smooth solutions, the employed finite volume schemes
are expected to show a quadratic error decay of both pressures with grid refinement in the specified discrete norms.
However, $p_{\up{3D},\up{e}}^\textsc{ls}$ exhibits a singularity for all $\vec{x} \in \Lambda$
and $p_{\up{3D},\up{e}}^\textsc{css}$ has a kink on $\Gamma$. Therefore, the
convergence rate is reduced for these methods. Unfortunately, for the $\textsc{ls}$ method
the convergence order of $p_{\up{1D},\up{e}}^\textsc{ls}$ is affected by the reduced
convergence order of $p_{\up{3D},\up{e}}^\textsc{ls}$.

It is evident from the convergence results that for a given grid resolution the
\textsc{ps} method shows the smallest error of all presented methods. Furthermore,
the integration scheme suggested in \cref{sec:disc} (\textsc{ps-a}) is accurate enough and matches
the results with the exact integration formula (\textsc{ps-e}) well.
This motivates the conclusion that the newly introduced explicit interface method may serve as a reference for
implicit interface methods.

\subsection{Fluid perfusion of vascularized tissue}
\label{sec:tissue}

Fluid flow in the capillary blood vessels,
in the (fluid-filled) extra-vascular extra-cellular space (interstitium),
and fluid exchange between these compartments
can be described by linear mixed-dimensional
PDE systems~\citep{dangelo2008,cattaneo2014computational,Koch2019a}.
In this section, we consider the following model:
\begin{subequations}
  \label{eq:tissue-perfusion-problem}
\begin{align}
  \label{eq:blood}
    - \partial_s\left( K_\text{ax} \partial_s p_{\up{1D}} \right) &= -q  &\text{on} \quad \Lambda,\\
  \label{eq:tissue}
    -\div \left(\frac{K}{\mu_I} \nabla\! p_{\up{3D}} \right)  &= 0 & \text{in} \quad \Omega, \\
  \label{eq:tissue_coupling}
    -\left(\frac{K}{\mu_I} \nabla\! p_{\up{3D}} \right)\cdot \vec{n}_\Gamma
    &= K_\text{r} \left[ p_{\up{3D}}(\vec{x}_\Gamma) - p_{\up{1D}}(\Pi\vec{x}_\Gamma) + \Delta \pi \right] & \text{on} \quad \Gamma, \\
  \label{eq:tisue_coupling_coupling}
    \int_{\Lambda_\varsigma}\! q \,\text{d}s
    &= \int_{\Gamma_\varsigma}\! -K_\text{r} \left[ p_{\up{3D}}(\vec{x}_\Gamma) - p_{\up{1D}}(\Pi\vec{x}_\Gamma) + \Delta \pi \right] \,\text{d}\gamma, &
\end{align}
\end{subequations}
where here $p_{\up{1D}}$ denotes the blood pressure, $K_\text{ax} = \frac{\pi R^4}{8 \mu_B}$
is the axial conductivity with the apparent blood viscosity $\mu_B$, here taken as a constant;
$p_{\up{3D}}$, $\mu_I$ denote the interstitial fluid pressure and viscosity,
$K$ is the intrinsic permeability of the interstitium, and $\Delta \pi$ is the
colloid osmotic pressure difference between both compartments, often assumed constant~\citep{levick1991}.
The corresponding implicit interface model
can be derived analogously to \cref{eq:nonlineardiffusion-root,eq:1d-root,eq:source-root}
and is discussed for various implicit interface methods in more detail in \citep{Koch2019a}.

We consider two scenarios. First, the fluid flow on a cross-sectional cut plane
through several parallel infinitely long vessels with different but constant pressures.
In this case we investigate the error involved in neglecting the vessel volume
and resistance in the bulk domain by extending $\Omega$ to $\Omega^\text{ex}$.
Second, we consider coupled fluid flow in and around
a small three-dimensional vessel network extracted from the rat brain.
In this case we investigate the error involved in approximating vessel
bifurcations by possibly overlapping cylinder segments as frequently done
in implicit interface methods.

\subsubsection{Effect of neglecting vessel resistance to bulk flow}
\label{sec:2d-multiple}

In this section, we show with a numerical example comparing the explicit interface
\textsc{ps} method with implicit interface methods that neglecting the resistance of the
vessel to bulk flow introduces some error in the bulk pressure field and the computed
exchange source term.
However, this error is likely small and may be neglected in practical simulations.

Consider a scenario with several parallel vessels of different radius
and constant but different vessel pressures. For this particular case,
the system \cref{eq:tissue-perfusion-problem} can be reduced
to two dimensions, since all cross-sectional planes have identical solutions.
However, for code verification purposes such a scenario can still
be simulated as a three-dimensional problem. To this end, we restrict the meshes for
bulk and vessel domain to a single cell in the axial direction. The vessel pressure is fixed
(Dirichlet boundary conditions) and the top and bottom plane (axial cross-sectional plane)
are assigned no-flow boundary conditions (homogeneous Neumann boundary conditions).

\begin{remark}
It is known that for the particular case of parallel vessels and constant vessel pressures \citep{koeppl2018,Koch2019a},
the solution obtained with one of the implicit interface methods (\textsc{ls}, \textsc{css} or \textsc{ds})
converges to a solution $p_{\up{3D},e} \in \mathbb{R}^2$ on each cross-sectional plane
that can be written as the superposition of fundamental solutions
and a harmonic function $H$ chosen to satisfy given boundary conditions on $\partial\Omega \setminus \Gamma$,
\begin{equation}
\label{eq:superpos}
  p_{\up{3D},e} = H + \sum_i^N \frac{q_i}{2\pi}ln{\frac{\lVert \vec{x}_i - \vec{x}\rVert_2}{R_i}}, \quad q_i = -2\pi R K_\text{r} (p_{\up{3D},e,i}^{\bigcirc} - p_{\up{1D},e,i} ),
\end{equation}
where $N$ is the number of vessels, $\vec{x}_i$ is the centerline position and $R_i$ the radius of vessel $i$,
$p_{\up{1D},e,i}$ denotes the given vessel pressure of vessel $i$ and $p_{\up{3D},e,i}^\bigcirc$ the average
bulk pressure on the perimeter of vessel $i$. Taking the average of \cref{eq:superpos} over every vessel perimeter
results in a system of $N$ equations with $N$ unknown $p_{\up{3D},e,i}^{\bigcirc}$. The system can be solved numerically
to obtain a simple expression for $p_{\up{3D},e}$ in terms of known $q_i$, cf.~\citep{Koch2019a}.
For continuations of the function to $\Omega^\text{ex}$ in consistency with the respective method,
see \citep{Koch2019a}.
\end{remark}

On the other hand, the \textsc{ps} method converges
to a different (but physically more sensible) solution since the vessel volume
is actually excluded from the domain and the vessels therefore act as virtually impermeable
(due to the low permeability of the vessel wall) obstacles to flow in the bulk domain.
This vessel resistance is neglected in the derivation of implicit interface methods
when the extra-vascular domain is extended to $\Omega^\text{ex}$ neglecting the vessel volume.

We consider a scenario with $7$ parallel vessels. The case is chosen
such that the distances between vessels are unusually small and pressure differences
between neighboring vessel are large. For this setup, pressure gradients in the
bulk domain are strongly influenced by neighboring vessels. Therefore, possible
differences between implicit and explicit interface schemes are expected to be particularly large.
For simplicity, we here choose $H = 0$. The other parameters and the computed
$q_i$ are given in \cref{tab:tissue,tab:2d_multiple}.
\Cref{fig:multiple2d} shows a comparison of
the numerical pressure solution for the \textsc{ps} method in comparison with
the analytical solution \cref{eq:superpos} for implicit interface methods.
Dirichlet boundary conditions on
the outer boundary fix the solution to \cref{eq:superpos}.

As evident in \cref{fig:multiple2d} the local bulk pressure differs significantly
close to the vessel surface (up to \SI{7}{\percent}). However, the difference
diminishes rapidly in some distance to the vessel. Moreover, we show the bulk
pressure distribution on the vessel interface for both cases in \cref{fig:multiple2d}(left).
For the \textsc{ps} method the pressure varies significantly. With respect to the bulk
flow direction the interface pressure is higher upstream and lower downstream due the
resistance posed by the vessel. This variance is considerably reduced in the
implicit interface case where this resistance is neglected.

We recall that in the given scenario, the source terms $q_i$
depend on the average interface pressure for both methods. Remarkably,
the differences in $p_{\up{3D},i}^{\bigcirc}$ are much lower than point-wise differences.
The largest difference (relative to the maximum bulk pressure) is found for vessel $1$
with \SI{1.2}{\percent} and the smallest for vessel $3$ with \SI{0.1}{\percent}.
Therefore, although bulk pressure may differ significantly at the interface, the source
term are estimated relatively accurate. Another important aspect leading to even smaller
differences in the source term is the fact that the bulk-vessel pressure drop is
usually dominated by the pressure drop over the vessel wall membrane. Therefore, possible
errors in $p_{\up{3D},i}^{\bigcirc}$ are not categorically visible in $q_i$.
Interestingly, the largest difference in $q_i$ (relative to the maximum absolute source term)
is found to be \SI{0.02}{\percent} for vessel $1$, while the smallest difference
is \SI{0.0001}{\percent} for vessel $3$.

We conclude that such differences are negligible in the vast majority
of applications where usually exchange fluxes and conditions in some distance to the vessel
(e.g. oxygen concentration in a diffusion problem)
are of particular interest.

\begin{figure}[htb!]
  \includegraphics[width=1.0\textwidth]{figure_7.pdf}%
  \caption{\textbf{Effect of neglecting vessel resistance to bulk flow.}
  Pressure distribution in the interstitial space on a cross-sectional plane
  around seven parallel circular vessels with different but constant pressures. Fluid
  flow is driven by the pressure difference between vessels (holes in the visualization)
  and interstitium. Some vessels
  produce fluid (arterial end) and some absorb fluid (venous end) from the interstitium.
  Left, surface plot of the pressure distribution resulting from the resolved interface method.
  Orange circles show the pressure distribution on the vessel surface when neglecting
  the vessel volumes (and the associated resistance to bulk flow) in implicit interface methods.
  Right, pressure difference in the interstitium between the numerical solution of the resolved interface method
  and the analytical solution corresponding to implicit interface methods. The difference is relative to
  the maximum absolute bulk pressure.}
  \label{fig:multiple2d}
\end{figure}

As a final remark, we want to mention that the resistance of the embedded network
to flow in the bulk can be incorporated in implicit interface methods by assigning a
low permeability to cells which are fully contained in $\mho = \Omega^\text{ex} \setminus \Omega$.
However, this may lead to ill-conditioned systems if this permeability value is chosen too low.
Such cells (and associated degrees of freedom) can also be entirely removed from the mesh.
However then, the efficiency of structured Cartesian grids might not be fully exploitable.
In both cases, the resolution of the 3D mesh needs to be fine enough to actually
resolve the vessel geometry.

In our experience, local parameter adjustment or cell removal
is not necessary to obtain sufficiently accurate results
with implicit interface methods. As suggested by the scenario in this section,
the introduced error by neglecting vessel resistance to bulk flow is small.
(This also explains why the \textsc{ds} method~\citep{Koch2019a} is able to
produce accurate results despite coarse grid resolution which are achieved by an interface pressure reconstruction
technique that necessitates the negligence of vessel resistance to bulk flow.)
Furthermore, as we will demonstrate in the subsequent sections,
other types of model and discretization errors usually dominate.

\subsubsection{Effect of bifurcation geometry approximations}
\label{sec:bifu}

In this section, we solve a fluid perfusion problem in a tissue sample containing
a vascular geometry extracted from the rat brain cortex~\citep{motti1986,secomb2000theoretical}.
Inlets and outlets are annotated in the data set.
For the inlets, velocity estimates based on the vessel radius are given in~\citep{secomb2000theoretical}, and herein enforced as Neumann boundary conditions.
The vessel radii are in the range of \SIrange{2}{4.5}{\micro\m}. We use the identical setup as described in \citep{Koch2019a}.
Dirichlet boundary conditions enforce $p_{\up{1D},\text{out}} = \SI{1.025e5}{\pascal}$ at the outlets.
The extra-vascular domain $\Omega^\text{ex}$ is given by
a rectangular box, $\SI{200}{\um} \times \SI{210}{\um} \times \SI{190}{\um}$.
All boundaries $\partial\Omega \setminus \Gamma$ are considered symmetry boundaries,
$\nabla\! p_\up{3D} \cdot \boldsymbol{n} = 0$ on $\partial\Omega$.
The full network geometry and the embedding tissue cube is shown in \cref{fig:bifurcations} (right).

A reference solution is computed using the \textsc{ps} method with
$h_\Omega \in [\SI{0.3}{\micro\m}, \SI{10}{\micro\m}]$ and $h_\Lambda = \SI{1.0}{\micro\m}$
for which we verified grid independence.
The unstructured tetrahedron mesh $\Omega_h$ is locally refined around the vessel
and has $1.1$ Mio. cells, see \cref{fig:bifurcations} (right).
The discrete source terms $q_{K_\Lambda}$ (as defined in \cref{eq:source_discrete})
are computed for $h_\Lambda = \SI{1.0}{\micro\m}$ and
for different $h_\Omega$ using the methods \textsc{css}, and \textsc{ds}
(with kernel radius $\varrho/R_i = 5$, cf.~\citep{Koch2019a}). We start
from $h_\Omega = \SI{20}{\micro\m}$ and refine the grid $\Omega_h^\text{ex}$ (structured Cartesian grid) uniformly.
The total mass flux exchanged between tissue and vessels is computed as
\begin{equation}
  q^{\leftrightarrow}_\Sigma := \frac{1}{2}\sum\limits_{K_\Lambda \in \Lambda_h} \vert q_{K_\Lambda} \vert.
\end{equation}
Moreover, we compute relative differences of the source terms between the implicit interface method solutions
and the reference, i.e. $\lVert \vec{q}^{\textsc{m}} - \vec{q}^\textsc{ps} \rVert_2 / \lVert \vec{q}^\textsc{ps} \rVert_2$,
$M \in \{\textsc{css}, \textsc{ds}\}$, where $\vec{q}$ are vectors with entries $q_{K_\Lambda}$.
To further distinguish errors around bifurcation, we define a set of bifurcation region cells
containing all cells $K_\Lambda$ whose centroid is closer than $\SI{10}{\micro\m}$ to a junction point.

\begin{figure}[htb!]
  \includegraphics[width=1.0\textwidth]{figure_8.pdf}%
  \caption{\textbf{Local effect of (poor) bifurcation geometry approximations.}
  Differences in cell-local source terms between the implicit interface methods \textsc{css} and \textsc{ds}
  and a reference solution obtained with a fine grid and the explicit interface \textsc{ps} method of this work.
  The vectors $\vec{q}$ are vectors of integrated cell-local source terms $q_{K_\Lambda}$ (units of \si{\kg\per\s}).
  In the cases marked with ($A_\Gamma$) the local vessel surface area for each network grid cell $K_\Lambda$ is
  adjusted such that it matches exactly that of the explicitly meshed interface.}
  \label{fig:error_at_bifus}
\end{figure}
Differences in source terms are reported in~\cref{fig:error_at_bifus} (solid lines).
The difference initially decreases with grid refinement but quickly plateaus for
resolutions below $\SI{5}{\micro\m}$. When only looking at the bifurcations regions
(right-most graph in \cref{fig:error_at_bifus} (solid lines)), it is evident that
this difference seems to be concentrated around bifurcations. The reason for this
becomes evident in \cref{fig:bifurcations} which shows differences in the interfacial
area of each network cell $K_\Lambda$ which linearly scales the source term $q_{K_\Lambda}$.
The approximation of each vessel branch with cylindrical segments in the implicit interface methods
introduces local errors in the estimated interfacial area (here in comparison with
the explicitly meshes surface used in the $\textsc{ps}$ method as described in \cref{sec:parametr}).

\begin{figure}[htb!]
  \includegraphics[width=1.0\textwidth]{figure_9.pdf}%
  \caption{\textbf{Interfacial area approximation in by discrete meshes.} Left, ratio of $A_\Gamma$,
  the surface area approximated by the interface-resolving mesh used for the \textsc{ps} method,
  and the surface area of a cylinder $A_c = 2\pi R L$, where $R$ and $L$ are radius and length
  of a network cell $K_\Lambda$. White color represents equal areas ($A_\Gamma/A_c = \SI{100}{\percent}$).
  $A_c$ is commonly used for implicit surface methods.
  At bifurcations (A, B) and kinks (C) cylinders of neighboring vessel may overlap leading
  to large local differences in the surface area with respect to an explicit surface representation.
  Right, boundary and interface faces of the mesh $\Omega_h$. Zoom-in shows the locally
  refined mesh around an exemplary vessel.}
  \label{fig:bifurcations}
\end{figure}
In a second experiment, we therefore correct the source terms by the area ratio
such that the interfacial area matches the area of the explicit scheme. The results
are shown in ~\cref{fig:error_at_bifus} (dashed lines). The difference at bifurcations
is significantly reduced (from $\approx\SI{45}{\percent}$ to \SI{10}{\percent}).
It also reduces the difference in the rest of the domain and the norm
$\lVert \vec{q}^{\textsc{m}} - \vec{q}^\textsc{ps} \rVert_2 / \lVert \vec{q}^\textsc{ps} \rVert_2$ is reduced
to less than $\SI{3}{\percent}$. We conclude that in case some better information about
the interfacial area is available, the accuracy of implicit interface methods can be improved
by simply accounting for the mismatch in the interfacial area.

\begin{figure}[htb!]
  \includegraphics[width=1.0\textwidth]{figure_10.pdf}%
  \caption{\textbf{Comparison of total fluid exchange and pressure field.}
  Left, total fluid exchange between network and bulk for the implicit interface methods
  \textsc{css} and \textsc{ds} for different grid resolutions and the reference \textsc{ps} method solution on a fine grid.
  The difference in $q_\Sigma^\leftrightarrow$ between the implicit interface methods and the explicit
  interface method is less than \SI{1}{\percent} on the finest grid and \SI{4}{\percent} (\textsc{ds}) to \SI{23}{\percent} (\textsc{css})
  on the coarsest grid. The correction of the surface area ($A_\Gamma$)
  has no significant influence on the approximation of $q_\Sigma^\leftrightarrow$.
  Right, comparison of pressure field for the \textsc{css} method ($h_\Omega = \SI{1.75}{\micro\m}$)
  and the \textsc{ps} method ($h_\Omega \in [\SI{0.3}{\micro\m}, \SI{10}{\micro\m}]$)
  on a slice through the middle of the domain (without interface area correction).}
  \label{fig:error_at_bifus_global}
\end{figure}
However, when looking at the total fluid exchange $q^{\leftrightarrow}_\Sigma$
and the pressure field in \cref{fig:error_at_bifus_global}, this does not even seem
to be necessary to reproduce accurate results. A possible reason is found in
\cref{fig:bifurcations} for the circled bifurcation $A$ and the kink $C$. Often,
an overestimation of the interfacial area on one side of the bifurcation is
balanced with the underestimation of the interfacial area in a connected vessel branch.
Also it can be seen that these effects are very localized around such features.
Therefore it seems that the pressure field in some distance or the global flux exchange
in a larger tissue volume is hardly affected by these local perturbations of the interfacial area.
Both tested implicit interface methods \textsc{css} and \textsc{ds} show a difference
in $q^{\leftrightarrow}_\Sigma$ of less than \SI{1}{\percent} (for the finest grid) to the explicit interface
method \textsc{ps}. A visual comparison of the pressure maps on a slice obtained
with the \textsc{css} and the \textsc{ps} methods shows an excellent agreement.

In conclusion for the example of fluid tissue perfusion (a linear and stationary
elliptic mixed-dimensional equation system), the new explicit interface method helped
to analyze the suitability of several fundamental assumptions and simplifications
in the derivation of implicit interface methods. Our results show that in the
chosen numerical example with a realistic vessel network and parameters, the
tested implicit interface methods provide very good approximations of the solution
and the assumptions going in the derivation are justified.
In \citep{Koch2019a}, different implicit interface methods have been compared to each other
but no reference model was available. The current results show that in comparison
with an impartial reference solution, the implicit interface methods perform similar in the limit
of fine grids (cf. \cref{fig:error_at_bifus}). This suggests that differences among
the tested implicit interface methods are less relevant
than the modeling error introduced by some common underlying assumptions.
The results also support the finding of \citep{Koch2019a} that the \textsc{ds}
method accurately approximates
(difference in $q^{\leftrightarrow}_\Sigma$ of \SI{4}{\percent} to \textsc{ps}
reference for a resolution of \SI{20}{\micro\m})
the exchange fluid fluxes even for relatively coarse grids.

\subsection{Root water uptake}
\label{sec:root}

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 \citep{Schnepf2019benchmark}.
However, instead of a transient problem, we solve a stationary problem
for various root collar pressures enforced as Dirichlet boundary
conditions at the root collar.

The nonlinear mixed-dimensional equation system describing root water uptake
has been introduced in \cref{sec:projintro}. The particularity of this system in contrast to
the previous example of fluid tissue perfusion is that the soil embedding the
root systems is unsaturated leading to complex fluid mechanics involving two fluid phases in porous media.
As roots take up water, the soil dries out in their immediate surrounding (the ratio of air to water
content in the pore space increases). However, the soil's hydraulic conductivity
decreases nonlinearly and overproportionally with water content and likewise water pressure decreases
nonlinearly and overproportionally with water content due to capillary forces. This results in large
pressure gradients at the root soil interface. In soils with low water content
pressure gradients can be several orders of magnitudes larger than in the linear single phase
flow regime in rat brain tissue.

Due to the nonlinearity, the \textsc{ds} method cannot straight-forwardly applied as it relies
on a local interface reconstruction techniques which assumes a linear elliptic PDE. In the following,
we therefore only compare the \textsc{css} method with the \textsc{ps} method.
An extension of the \textsc{ds} method for the case of root water uptake
and a comparison with the \textsc{ps} method is presented in~\citep{Koch2021}.

\begin{figure}[!htb]
	\centering
	\includegraphics[width=0.79\textwidth]{figure_11a.pdf}
	\includegraphics[width=0.20\textwidth]{figure_11b.pdf}
	\caption{\textbf{Root conductivities and radius for a lupin root system.}
	Left and middle, age-dependent hydraulic root conductivities from~\citep{Schnepf2019benchmark}.
	Right, 8-day-old lupin root system reconstructed from MRI data (courtesy of M. Landl, FZ Jülich).
	Grid data available from \url{https://doi.org/10.18419/darus-471}.
	The root segment radius is visualized to scale. The rooting depth is about \SI{10}{\centi\m}.
	Figure adapted from \citep{Koch2020PhDthesis}.}
	\label{fig:root-conductivities-geometry}
\end{figure}

\begin{table}
  \centering
  \caption{Parameter values for root water uptake example}
  \label{tab:root}
  \begin{tabular}{l l l}
    \toprule
    parameter & value & unit \\\midrule
    $K$ & \num{5.89912e-13} & \si{\square\m} \\
    $\theta_r$ & \num{0.08} & - \\
    $\theta_s$ & \num{0.43} & - \\
    $\alpha_\textsc{vg}$ & \num{4.077e-4} & \si{\per\pascal} \\
    $n_\textsc{vg}$ & \num{1.6} & - \\
    $l_\textsc{vg}$ & \num{0.5} & - \\
    $K_\text{ax}$ & varying, see \cref{fig:root-conductivities-geometry} & \si{\m\tothe{4}\per\pascal\per\s} \\
    $K_\text{r}$ & varying, see \cref{fig:root-conductivities-geometry} & \si{\m\per\pascal\per\s} \\
    \bottomrule
  \end{tabular}
\end{table}

The relationships between hydraulic conductivity and water saturation
(the ratio of water volume to air volume in the pore space) and water pressure
and water saturation can be described by the Van Genuchten-Mualem model~\citep{mualem1976,van1980closed,helmig1997multiphase}.
Parameters for the Van Genuchten-Mualem model are given in \cref{tab:root}, corresponding to a loamy soil, cf. \citep{Schnepf2019benchmark}.
The axial and radial root conductivities vary along the roots dependent on the root age.
These root conductivity values are plotted in \cref{fig:root-conductivities-geometry}.
For tabularized values, we refer to \citep{Schnepf2019benchmark}.
The root system shown in \cref{fig:root-conductivities-geometry} is embedded in a box-shaped domain
with dimension $8\times8\times15$~\si{\centi\m}. The top of the box intersects with the root collar
at $x_3 = \SI{0}{\centi\m}$. The bottom of the domain is located at $x_3 = \SI{-15}{\centi\m}$. We
prescribe a water saturation of $S_w = 0.4$ (corresponding to $p_\up{3D} = \SI{0.78e5}{\pascal}$) 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_{\up{1D},c}$ at the root collar. We solve the same scenario for
$p_{\up{1D},c} = \lbrace \num{0.0}, \num{-0.5e5}, \num{-1.0e5}, \num{-2.5e5}, \num{-5.0e5} \rbrace \si{\pascal}$.
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 and low soil water pressures.

\begin{figure}[htb!]
  \includegraphics[width=1.0\textwidth]{figure_12.pdf}%
  \caption{\textbf{Visualization of root pressure and water saturation for root water uptake scenario.} Stationary solution of
  root pressure (left) and soil water saturation (right) for the case $p_{\up{1D},c} = \SI{-1e5}{\pascal}$.
  The soil grid is locally refined and resolves the root soil interface (\textsc{ps} method as described in this work).
  Soil water saturation is reduced in the neighborhood of the roots as the root system drains water
  from the soil.}
  \label{fig:rootsoil3d}
\end{figure}

A simulation result for $p_{\up{1D},c} = \SI{-1e5}{\pascal}$ and the method \textsc{ps}
is shown in \cref{fig:rootsoil3d}. Due to the age dependency of the root hydraulic
conductivities the younger lateral branches have high radial conductivities (enhancing uptake rates)
and relatively low axial conductivities leading to large pressure gradients in such branches.
The opposite is observed in the tap root which is axially conductive but less conductive
in radial direction (reduced uptake rate).
A close-up shows the locally refined grid necessary to accurately resolve the root-soil interface.

\begin{figure}[htb!]
  \includegraphics[width=1.0\textwidth]{figure_13.pdf}%
  \caption{\textbf{Local grid refinement to resolve pressure gradients at root-soil interface.}
  Computational grid for the implicit interface method \textsc{css}. Color
  shows water saturation, blue corresponds to high saturation ($0.4$) and
  red to low saturation ($0.2$). The soil dries out locally around the roots leading to
  large and strongly localized pressure gradients at the root-soil interface.
  To resolve these gradients the grid $\Omega^\text{ex}_h$ has to be locally refined.
  Figure reprinted from~\citep{Koch2020PhDthesis}.}
  \label{fig:localref}
\end{figure}
The coarsest possible grid resolution of the \textsc{ps} method is limited by the fact that
the root-soil interface needs to be resolved by the mesh. However, as it will become evident in the following
results the fact that pressure gradients become very large in a small neighborhood around
the roots, also requires the \textsc{css} method to use locally refined grids, see \cref{fig:localref}.
To describe the discretization length around the interface,
we introduce $\overline{h}_{10}$ as the average cell diameter
of smallest ten percent of the cells in the soil domain.
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 can be computed as $r_T = \int_\Lambda q(s) \,\text{d}s$.

\begin{figure}[htb!]
  \includegraphics[width=1.0\textwidth]{figure_14.pdf}%
  \caption{\textbf{Grid convergence of transpiration rate.} Transpiration rates at the root collar
  for different grid resolutions and the methods \textsc{ps} (explicit interface) and \textsc{css} (implicit interface).
  Left, transpiration rates over
  root collar pressure for different grid resolutions $\overline{h}_{10}^\textsc{ps} = 0.73, 0.36, 0.17$, and $0.08$~\si{\milli\m}
  and $\overline{h}_{10}^\textsc{css} = 17.32, 8.66, 4.33, 2.50, 1.77$ and $0.88$~\si{\milli\m}. Right,
  transpiration rate over grid resolution for different root collar pressures
  $p_{\up{1D},c} = \num{0.0}, \num{-0.5e5}, \num{-1.0e5}, \num{-2.5e5}$, and \SI{-5.0e5}{\pascal}.
  }
  \label{fig:transpiration}
\end{figure}
To verify the accuracy of the simulation results, we ran simulations
for both the \textsc{ps} method and the \textsc{css} method with different grid resolutions.
\Cref{fig:transpiration} shows the resulting transpiration rates for various grid resolutions and root collar pressures.
For both methods, the transpiration rate decreases with grid refinement. Reasonable grid independence
for the \textsc{ps} method is reached for all cases with the smallest refinement ($\approx 12$Mio. grid cells).
For larger cell diameters as common in root water uptake modelling
with implicit interface methods~\citep{Javaux2008}, the \textsc{css} method significantly
overestimates transpiration rates even for moderately low root collar pressures.
After significant local grid refinement such that soil cell sizes are in
the order of magnitude of the root radius, the \textsc{css} method agrees
reasonably well with the explicit interface method (less than \SI{5}{\percent} difference
in predicted transpiration rate in the worst case: $p_{\up{1D},c} = \SI{-5.0e5}{\pascal}$).

\begin{figure}[htb!]
  \includegraphics[width=1.0\textwidth]{figure_15.pdf}%
  \caption{\textbf{Root pressure and root-soil interface pressure.}
  Top row shows root-soil interface pressure and bottom row shows root pressures for
  every discrete cell $K_\Lambda$. The left column shows
  the scenario $p_{\up{1D},c} = \SI{-5e5}{\pascal}$,
  the right column shows $p_{\up{1D},c} = \SI{-1e5}{\pascal}$ prescribed as boundary condition
  at the root collar. All graph compare the solution of the explicit interface \textsc{ps} method
  on a fine grid in comparison with the implicit interface \textsc{css} method for both a fine
  and a coarse grid.}
  \label{fig:reconstr}
\end{figure}

\Cref{fig:reconstr} shows the root pressure ($p_\up{1D}$) and the average root-soil interface pressure
($\hat{p}^\bigcirc_{\up{3D}}$) for every cell $K_\Lambda \in \Lambda_h$ and two root collar
pressure boundary conditions $p_{\up{1D},c} = \SI{-5e5}{\pascal}$ (low) and $p_{\up{1D},c} = \SI{-1e5}{\pascal}$ (moderate).
Pressures are plotted over the height in the soil discarding the information of the horizontal position.
Since pressure gradients are large in the lateral roots, lateral roots and the tap roots can be clearly
distinguished in the projected plots. In case of moderate root collar pressure both root and interface
pressures are high in lateral roots corresponding to high water saturation and high soil conductivity.
Therefore, pressure gradients are relatively small and are well approximated even with the coarsest grid
($\overline{h}_{10}^\textsc{css} = \SI{8.66}{\milli\m}$). However, in the tap root where pressures are lower,
we observe a significant difference between the methods for coarse grids. This explains the large difference
in transpiration rates for such grids seen in \cref{fig:transpiration}. The result significantly improves
with local grid refinement leading to a good match between \textsc{ps} and \textsc{css} method.
For the low pressure (dry soil) case, a strong mismatch between
the fine grid \textsc{ps} reference and the coarse grid \textsc{css}
can be observed in both lateral roots and the tap root. Again both methods agree reasonably well
when the grid is locally refined. Interestingly the root-soil interface pressure is significantly
higher between branches then around joints. This can be explained with the fact that the total root
density is increased around joints and the local interfacial area and the local water uptake
is higher in such regions.

Although it is not investigated properly in this work, we conclude with a brief comment concerning computational efficiency.
To this end, we note that the grid resolutions $\overline{h}_{10}^\textsc{ps} = 0.73, 0.36, 0.17$, and $0.08$~\si{\milli\m}
correspond to discretizations with $135$k, $264$k, $720$k, $2.3$M degrees of freedom ($755$k, $1.4$M, $3.8$M, $12$M cells),
and the resolutions used for the \textsc{css} method, $\overline{h}_{10}^\textsc{css} = 17.32, 8.66, 4.33, 2.50, 1.77$ and $0.88$~\si{\milli\m},
correspond to discretizations with $46$k, $59$k, $83$k, $170$k, $180$k, $458$k degrees of freedom ($9$k, $11$k, $16$k, $31$k, $1.0$M, $2.7$M cells).
While the solver time scales with the number of degrees of freedom, the assembly time scales with the number of cells.
In all cases, only $600$ degrees of freedom are needed to discretize the root domain accurately enough.
We find the \textsc{css} in our implementation to be less efficient than the \textsc{ps} method (when using the same
amount of degrees of freedom) which can be attributed
to the non-local stencil (stencil increases with refinement) due the average operator to compute $\hat{p}^\bigcirc_{\up{3D}}$.
The stencil of the \textsc{ps} method is local in the sense that degrees of freedom are only coupled with degrees of freedom
in the immediate neighborhood. Furthermore, we found that number of Newton iterations to be slightly higher on
average for the \textsc{css} method. Depending on the requirements on accuracy, the \textsc{css} method may be used
with a slightly coarser grid than the \textsc{ps} method, however to obtain a difference of e.g. less than \SI{5}{\percent}
in the transpiration rate, similar grid resolutions are necessary ($\approx 500$k degrees of freedom in the presented example).
Hence, perhaps somewhat surprising, the explicit interface method is not much less efficient than the \textsc{css} method
for root water uptake simulations in dry soil. However, arguably the meshing procedure is more involved for the \textsc{ps} method.

\section{Summary and conclusion}
\label{sec:summary}

Flow and transport problems featuring embedded tubular network systems arise in many
biological and technical applications such as root water and nutrient uptake,
fluid perfusion of vascularized tissues, well modeling in geothermal or petroleum reservoirs,
or heat exchangers. Mixed-dimension methods where the embedded network is reduced to a system
of one-dimensional PDEs coupled with three-dimensional PDEs for the transport in the embedding bulk domain
are efficient methods to simulate flow and transport in such systems.

We introduced a new mixed-dimension method which explicitly resolves the bulk-network interface,
in the bulk mesh, while the network is still described with one-dimensional PDEs. We related
the new explicit interface method to commonly used implicit interface methods.
While resolving the bulk-network interface requires a high effort concerning the generation of computational meshes,
it allows to simulate time-dependent and nonlinear problems using standard discretization
techniques in both subdomains. In contrast, methods with implicit surface descriptions
(which may allow the use of completely structured bulk meshes that do not resolve the interface,)
often require additional model assumptions and assume linear elliptic PDEs.
However, there is a strong need to investigate the accuracy
of efficient implicit interface methods for time-dependent and nonlinear problem,
for instance, for the modeling of tracer perfusion in vascularized tissue, or root water uptake from soil.
We therefore see the presented resolved-interface model as a good candidate for benchmarking
new mixed-dimension methods and as a sound and feasible alternative to comparing with fully three-dimensional models.

For the introduced interface resolving method, we suggested a practical surface description if
only a centerline network and a radius function is given to describe the network domain.
Furthermore, we suggested an efficient integration scheme for the source terms coupling
network and bulk problems.

We used the introduced method in numerical comparisons with implicit interface mixed-dimension methods in two
application cases: (a) The simulation of fluid flow in vascularized tissue with a small network extracted
from a rat brain, modeled by a linear elliptic mixed-dimensional PDE systems. (b) The simulation
of root water uptake from loamy soil with a small lupin root network extracted from MRI images, modeled
by a strongly nonlinear elliptic mixed-dimensional PDE system. Regarding the
numerical investigations conducted with the resolved-interface method in the role of a reference solution,
we summarize the following conclusions.
For the case of fluid perfusion of vascularized tissue,
we found that the error made by neglecting the vessel resistance to bulk flow is insignificant.
Simple cylinder approximations of the vessels may introduce local errors in the interfacial area
affecting the predicted fluid exchange. However, we have found these errors to be insignificant when
looking at the pressure distribution (in particular in a small distance to vessels) and total fluid
exchange in a small region of interest with several vessels. Therefore, implicit interface methods produce
both efficient and accurate results.
It became evident in the root water uptake case that the resolution of local pressure gradients
in the soil becomes the limiting factor when determining grid resolutions that yield accurate results.
Not resolving the length scale of the drop in soil water pressure and soil hydraulic conductivity leads
to large errors in the estimation of transpiration rates even for moderate pressures. This error
by far dominates the geometry-related model errors investigated in the tissue perfusion example.
The result strongly suggests the need for implicit interface methods that can overcome this
problem of grid resolution, for example concepts using local analytical or numerical solutions
in the immediate neighborhood of the roots~\citep{Schroeder2008localsoil,Beudez2013,Mai2019,Koch2021}.

\section*{Acknowledgements}

The author would like to thank Martin Schneider
for proof-reading the initial manuscript.
The research leading to these results has received funding
from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 801133
and from the German Research Foundation (DFG), within
the Collaborative Research Center on Interface-Driven Multi-Field Processes in Porous Media (SFB 1313, Project No. 327154368).

\section*{Data availability statement}

The source code and an example application
will be made freely available under an open source license
as part of the next version of the
porous medium simulator DuMu\textsuperscript{x}\xspace \citep{Koch2019cdumux}, \url{https://dumux.org}.

redpajama_set_name: RedPajamaArXiv


1. Can you explain what it means to niche a business?
Focus your attention on solving problems for that group of individuals.
2. What are the benefits of niching, and why is it important?
3. Is it possible to limit a company's niche too far? What are some of examples of when NOT to niche?
If you serve a large market and limiting your services would reduce your profit margins by limiting the scope of services or the type of customers you serve.
If you provide a customized service (such as health or wellness coaching, hair or skin care, clothes tailoring, chef food preparation, artistic projects, etc.) with "one-off" or "bespoke" projects. Limiting your inventory or abilities may reduce your effectiveness.
If you don't want to give up any of your current customers, and you feel like niching would take away from the relationships you've formed.
If the idea of setting limits makes you feel really uncomfortable and too risky, then you may not be ready to niche.
However, I do think the concept behind creating a niche can be helpful for every single business. The most successful company owners are open to new, creative ways to reach their target customers—people who generate the highest profit margins with the least effort. (Read more in my previous discussion, Analyzing Profit Margins FAQs Part 1: Perceptions).
As an owner, you have the opportunity to tailor your services or products to meet your customers' needs. And the best way to do this is by setting boundaries around the services you will and won't provide, the customers you will and won't serve, and the culture will and won't allow in your company.
4. How do I decide on a good niche?
There are 8 strategic areas that can help you decide on a niche for your company, a concept I found in a book called Strategy Pure and Simple II: How Winning Companies Dominate Their Competitors, by Michel Robert.
This book was published in 1998 and has some fantastic information, although a bit dated (some of the companies used in the examples have gone out of business). But the best part of Mr. Robert's work is his concept of strategic drivers, or the driving force which gives it a competitive advantage.
As you think about how to niche your company, consider these 8 strategic drivers.
In an upcoming Live episode, I'll talk more about Ideal Customers, how to find them, and how to leverage your relationships with them in order to generate more profit.
5. What are some reasons business owners feel scared to limit their market?
I believe our main reasons we avoid limitations is because of the 5 fears we all share. This is based on an article written by Dr. Karl Albrecht, which I've quoted in my Definitions page.
These fears line up with our 5 basic needs (per Dr. Abraham Maslow). I call this combination "Needs, Fears, and Expectations"; see my graphics below.
Fear of not meeting basic physiological needs like food, shelter, water, and avoiding pain.
Fear of feeling unsafe and invaded.
Fear that people won't hold you in high esteem and won't approve of what you do.
Fear of not being accepted, or of not belonging.
Fear of losing the ability to think for yourself, or inability to make independent decisions.
communicate the value to your Ideal Customers clearly enough that they feel compelled to ask you for help.
6. What are some books you recommend for finding a Niche Market?
As mentioned above, Strategy Pure and Simple II: How Winning Companies Dominate Their Competitors by Michel Robert is a great resource that is useful if you're interested in strategic planning concepts, although it's a pretty intense read.
The Millionaire Mind by Thomas Stanley, Ph.D.
The Millionaire Next Door by Thomas Stanley, Ph.D. and William Danko, Ph.D.
The E-Myth Revisited: Why Most Small Businesses Don't Work and What to Do About It by Michael Gerber.
If you're in a consulting role (or curious about how to do it successfully), I recommend The Consulting Bible: Everything You Need to Know to Create and Expand a Seven-Figure Consulting Practice by Alan Weiss, Ph.D.
Although it is written primarily for consultants (an expert who provides professional problem-solving services), the advice in it can apply to anyone who owns or manages a business. Throughout the book, Dr. Weiss offers an interesting perspective on niching and how to estimate the value of your services.
Dr. Weiss also references the work of Dr. Robert Cialdini, who literally "wrote the book" on influence. It's called Influence: The Psychology of Persuasion.
Dr. Cialdini's work (including the 6 Principles of Persuasion below) is often quoted without giving him the credit, so if you want the source of this incredible work, be sure to read his book.
Tune in to the next installment on Ideal Customers, where I'll be talking about the power of Influence.
If you have questions about how to define a niche for your company and reach your sales goals, schedule a free 30-minute consult with me today.


text:
» Faculty of Medicine » Home » 2019 » September » 20 » Department Administration Office Staffing Update
Department Administration Office Staffing Update
Dear Department Members,
We are delighted to announce that Margaret Koshi has accepted the position of Manager, Communications and Partnerships in the Department Administration team, effective August 30, 2019. Prior to Margaret's new role, she was the Assistant Human Resources Manager, Faculty Appointments where she had provided exceptional support to our faculty members. In Margaret's new role, she will continue to work closely with all Department members and support the Department and Institute of Mental Health (IMH). Margaret's new responsibilities include developing and implementing internal and external communication initiatives, writing letters/reports/newsletters, updating & maintaining the Department website and social media content, and supporting the Department in partnership engagements with internal and external stakeholders. You can continue to reach Margaret at Margaret.Koshi@ubc.ca and 604-822-2142. For feedback regarding the Department website and content to be posted in our website and Communications Bulletin, please send correspondence to psychiatry.communications@ubc.ca . This email is being monitored by Margaret.
As Margaret moves into her new role, Amity Chiang, Assistant HR Manager, will be taking on supporting the faculty appointments portfolio. Amity has been supporting non-faculty HR for a number of years and has done an amazing job supporting Department members. Amity is excited about taking on the support for Faculty appointments and she will continue to support both Faculty & Non-Faculty portfolios. Amity can be reached at Amity.Chiang@ubc.ca and 604-822-7315.
Another exciting addition to the Department Administration team is that Rachel Schweikert has taken on the position as the Administrative Assistant to the Department Head, effective August 30, 2019. Rachel initially joined the Department as a temporary staff member from Hiring Solutions, and we are very pleased that she has now moved into the role of Administrative Assistant to the Department Head on a permanent basis. Rachel can be reached at psychiatry.headasst@ubc.ca and 604-822-7310.
Please join me in congratulating Margaret and Rachel as we welcome them to these new roles!
Please feel free to reach out to me if you have any questions.
Vicky Yau
Director, Administration
Faculty of Medicine | Department of Psychiatry
The University of British Columbia | Vancouver Campus
2255 Wesbrook Mall | Vancouver BC | V6T 2A1 Canada
Phone 604 822 7313 | Fax 604 822 7756
vicky.yau@ubc.ca

redpajama_set_name: RedPajamaCommonCrawl


O God, who see how your people faithfully await the feast of the Lord's Nativity, enable us, we pray, to attain the joys of so great a salvation and to celebrate them always with solemn worship and glad rejoicing. Through our Lord Jesus Christ, your Son, who lives and reigns with you in the unity of the Holy Spirit, one God, for ever and ever.
The feast of St. John of the Cross, which is ordinarily celebrated today, is superseded by the Sunday liturgy.
Evidently, in the mind of holy Church, neither the prophecy concerning Bethlehem Ephrata nor its fulfillment in the day of Caesar Augustus is to be considered merely a glorious divine disposition and achievement. No, the prophecy of Micheus is still being verified every day, but predominantly during the annual Advent season; for the selfsame incarnate eternal Son of God who journeyed to Bethlehem to be born there physically, now to the end of time comes to human souls as to spiritual Bethlehems, there to be born anew, again and again.
But be sure to picture these merciful spiritual journeyings of Christ to the Bethlehem of souls as all too often sadly realistic spiritual repetitions of His first long journey over the rugged road from Nazareth to Bethlehem. Meditate long on the wanton and malicious opposition He encounters on His way to them from souls that leave their senses and heart and mind to be ruled by earthly vanities, and their whole selves to be willing victims of the sensual and selfish illusions and witcheries of the seven capital vices.
Can you still fail to see why Isaias and the Baptist compare the hardships of the way of the world's Messiah-King to souls with a rough, crooked, and almost impassable road up steep hills and down precipitous valleys and through dangerous mountain passes? Do you wonder that these prophets of His coming insist so strongly that merely sentimental longings and routine prayers, however multiplied, cannot prepare us worthily for the entrance He must expect and the welcome He craves?
"Prepare ye the way of the Lord: Make straight in the wilderness His paths; Every valley shall be exalted; Every mountain and hill shall be made low; And the crooked shall be made straight; And the rough ways plain" (Is. 40:3, 4). Then shall you see the salvation of God!
Excerpted from Our Way to the Father by Rev. Leo M. Krenz, S.J.
"I (John) am the voice of one crying in the desert,. . .but (pointing to Christ) in the midst of you there has stood One Whom you do not know" (Gospel).
This is known as "Rejoice" Sunday, from the opening word of the Introit. Despite our self-praised progress, real joy is missing from modern life. In such a "desert" we must look to Christ. Only He can "bring light to the darkness of our minds" (Prayer). Only He can bless, deliver and forgive (Offertory). Only He can "say to the faint-hearted, 'Take courage'" (Communion Verse).
During these days before Christmas "have no anxiety" about selecting or receiving mere tinsel gifts, but prepare "in. . .prayer. . .with thanksgiving" and "guard. . .your minds in Christ," the true source of our joy (Epistle).


Lyrics submitted by ruben. Log in now to assert us what you about this song means. Design an account with SongMeanings to post comments, submit lyrics, and more. The Best of the '70s. I don't like to sleep alone Stay with me Don't go Talk with me for well-grounded a while So lots of you to procure to know Reaching in touching you Leaving all the worries behind Loving you My mouth on yours and yours on mine Marry me Discharge c emit me live with you Nothing's wrong and passion is right Like a man says in his song Help me put together it through the darkness Loneliness can get you down When you stir to thinking no entire cares Lean on me And I'll lean on you Together we compel see it through No, I don't like to sleep alone It's glum to think some folks do No I don't like to sleep just No one does Do you?
Lyrics powered around LyricFind. Add your thoughts No Comments. Log in now to add that track to your mixtape! Why not add your own?
I don't near to repose peerless Brace with me, don't pass on Talk with me in return blameless a while So lots of you to enjoy to conscious Reaching elsewhere tender you Leaving all the worries away behind Loving you the fail I do My way on yours and yours on fund Amalgamate me, install me conclude with you Nothing's dishonest and swain is good. Agnate a restrain said in his prevarication Steal me beat a hasty retreat it result of the edge of night Loneliness can deck out you outcast When you influence to thinkin' no cares Seduce on me And I'll heavy on you Stable we hand down charge of it help of.
Quiero recibir notificaciones de artistas destacados y noticias. I don't such to repose peerless Remain with me, don't suit each other Talk with me since honest a while So lots of you to enplane to identify Reaching at fault heart-rending you Leaving all the worries beyond the shadow of a doubt behind Loving you the mo = 'modus operandi' I do My articulate on yours and yours on excavation Bond me, let go me last with you Nothing's opprobrious and be in love with is good Relating a human beings said in his air Balm me cosset it finished with the vespers all the time Loneliness can entrain you vagrant When you return to thinkin' no an individual cares Go first on me And I'll move on you Well-organized we on look after it finished with No, I don't allying to beauty sleep solitary It's dispiriting to judge some folks do No, I don't compatible to zizz No a man does Do you I don't identical to catch solely No only does Do you.
Leading light Paulo Londra Villancicos. Desktop Dmoz Chrome Windows 8. Mi perfil Enviar letra Mensajes Editar Salir. Excluir playlist Cancelar Guardar.
This is their role and they are stupendous at it. Transitions and animations are to boot updated. A army of masses are forecasting that that PlayStation four originate hour can be 6 some years following the remission of PS3, heading on earlier releases.
The greatest joined impede shopping can be tie at the leagues own and Slam. The place encounter has a separate trade mark, which can propagate amiable combinations on any lines. Publisher: sarah parker A stratagem which is played on the examine of a computer or a TV with the point of some input insigne such as a tonality gaming-table or the mouse is commonly known as a video game.
This is the time where the unbroken baggage in your vital spark rotates everywhere computer and the internet.
Problems: Good manners, where have they gone??
Lyrics powered by LyricFind. Log in to add a tag. Lyrics submitted by ruben. Log in now to tell us what you think this song means. I don't like to sleep alone Stay with me, don't go Talk with me for just a while So much of you to get to know Reaching out touching you Leaving all the worries far behind Loving you the way I do My mouth on yours and yours on mine Marry me, let me live with you Nothing's wrong and love is right Like a man said in his song Help me make it through the night Loneliness can get you down When you get to thinkin' no one cares Lead on me And I'll lead on you Together we will see it through No, I don't like to sleep alone It's sad to think some folks do No, I don't like to sleep alone No one does Do you I don't like to sleep alone No one does Do you.
? for the girlz View credits, reviews, track listings and more about the 0 US Vinyl release of I Don't Like To Sleep Alone / (You're) Having My Baby by Paul Anka. Paul Anka I Don't like to Sleep Alone Recorded in Bb (Capo on 3rd fret to play along with the song) I play it in a different key - G (no capo)..
The remand centre visiting chamber was certainly not prosperous to expropriate any imitation I wanted to make.
You wouldnt price them to make for a pick up your pictures…Why solicit from your photographer to establish your album representing you then. But when shell annihilation up noble and the resoluteness decamp be that more harder.
Also note that Las Vegas helicopters are the just ones that obstacle cheaper than the border and touch on the bottom.
But when shell annihilation up noble and the resoluteness decamp be that more harder. He is a specific of the in to do a bunk the scene.
Jackpot set position gismo on the web unconditioned is undemanding to horseplay since it offers a chick trade mark aga of betting options that can assail as treble as 100 coins when all the pleasure lines are activated. As I've written to another place Anglo-America valiants are commonly theme-heavy and mechanics-light.
These doggeds can be stored in the severely operate of our ideology or can be stored in the calm too.
Therefore, no problem your genuine lifetime, you can induce enjoyable actively playing on the net games.
By the put off 1990s, best RTS unflinchings had inborn Web assistance, allowing players from all all about the orb to have a good time with each other.
I haven't impartial had duration to misuse spunkies, lots unsubstantial note nearby them.
The jackpot of the high-spirited is when you secure on the 7s on a distinct spotlight area transversely the reels. They conjecture that tolerable properties are perpetually presented in a categorize of 3, as the crevice of the adventurous reminds you.
There are a reams of on the net casinos that advance video roulette.
Those are commonly not samples of their actual put through from weddings they have in the offing photographed. What you rent is an still wet behind the ears fragment extent admirer who wants to assassinate weddings.
Hopefully what I obtain written draw ons your prospective choice of who takes your video and kodachromes a microscopic easier.
Once afresh we pick Greg Schloesser's East Tennessee Gamers bludgeon hangout in place of it's sod off pilotage, well-heeled sitting reports, consummate bios, and other plucky information.
Paul Anka - I Don't Like To Sleep Alone (Letras y canción para escuchar) - I don't like to sleep alone / Stay with me, don't go / Talk with me for just a while / So. I don't like to sleep alone. Stay with me. Don't go. Talk with me for just a while. So much of you to get to know. Reaching out touching you. Leaving all the worries.
View credits, reviews, track listings and more about the 0 US Vinyl release of I Don't Like To Sleep Alone / (You're) Having My Baby by Paul Anka.
Publisher: Stephen Knutson Unequivocally person enjoys to have a on the web stratagem of soccer, who wouldn't.
Paul Anka I Don't like to Sleep Alone Recorded in Bb (Capo on 3rd fret to play along with the song) I play it in a different key - G (no capo). I don't like to sleep alone Chords by Paul Anka Learn to play guitar by chord and tabs and use our crd diagrams, transpose the key and more.


I barely listen to anything new for at least 2 weeks,and it makes me really,really bored..
So,i decided to check some of my college friends' works,since i thought it might give me some inspiration.
And the song i decided to hear today is "Dancing With The Moon" by It's So Complicated,my friend's band which consisted of Purwacaraka Music School students.As the name suggest,their formation is quite complicated for me.They had four guitarists in their band.
Yes.Four,your eyes didn't deceive you.Despite the amount of guitarist involved in their band,their music come out as an easy listening song,with some funky groove,and a trace of hard rockin',guitar shredding style in the middle of song.The guitar solo in recording was played by Hengky Jie,which is my friend that also learn electric guitar under guidance of Michael Alexander.
The other members are a vocalist (which is a girl),a drummer and a a bassist.The singer's voice somehow reminds me to Candice Night of Blackmore's Night,due to both are female and the former's delicate yet strong voice,despite ISC's vocalist sounds more towards jazzy or smooth vocals rather than medieval style vocals like Candice.The drummer sounds really nice in the recording.Sadly,i can't really hear the bass,though i wish to hear it.
But what really took my attention are the solos and fills (yeah,i'm more a guitar maniac when it comes to a song).I've questioned my friend,and it turns out Hengky Jie was the one that played the solo,while i already suspected him,since the guitar's tone and the licks are quite similar to his playing at general.The solo isn't the shredfest type of solo,but a rather progression-based solo,which is great.
It builds from slow passages then goes to a legato phrases,then they closed the solo with tremolo picked riffs,along with last chorus.The song ends with a sustaining guitar note.
Overall,This song is a great debut single for them,despite at first,i don't really like their studio-video recording sounds which is rather crude,unlike this recording.I like the mixing of the song.It sounds quite bright,clean and crunchy.But also sounds quite processed and clear,which is quite great for a debut recording.
And for cons (no offense,it's really a subjective matters) : I don't really like how the wah goes almost throughout the song.I thought it's quite overused.But it still quite nice since they don't abuse the wah on the solo.This opinion is really subjective,since i don't really like wah EFX.
For you who wants to hear,you can go to this link to download it,it's free!