Population dynamics in river networks

TL;DR

This paper develops a mathematical framework to analyze how the structure of river networks influences population persistence, incorporating models, eigenvalue analysis, and numerical simulations to identify key factors affecting species survival.

Contribution

It introduces a process-oriented model for population dynamics in river networks, establishing fundamental theories and deriving thresholds for persistence based on eigenvalue analysis.

Findings

River network structure significantly impacts population persistence.

Hydrological, physical, and biological factors influence species survival.

Eigenvalue-based thresholds can predict population extinction or persistence.

Abstract

Natural rivers connect to each other to form networks. The geometric structure of a river network can significantly influence spatial dynamics of populations in the system. We consider a process-oriented model to describe population dynamics in river networks of trees, establish the fundamental theories of the corresponding parabolic problems and elliptic problems, derive the persistence threshold by using the principal eigenvalue of the eigenvalue problem, and define the net reproductive rate to describe population persistence or extinction. By virtue of numerical simulations, we investigate the effects of hydrological, physical, and biological factors, especially the structure of the river network, on population persistence.

Tables

Table 1. Table 5.1: The unit of parameters.

Parameter	$L$ & $l_j$	$D_j$	$v_j$	$A_j$	$r_j$	$m_j$	$Q_j$	$n_j$	$B_j$	$S_{0j}$	$y_j$
Unit	m	m${}_{}{}^2/$s	m/s	m2	1/s	1/s	m3/s	s/m1/3	m	m/m	m

Equations

$$E=\{e_i:1\le i\le N\},$$

$$K=\{k_j:j\in I_{N-1}\},n\ge 2,$$

$$I_{N-1}=\{1,2,\cdots ,N-1\}.$$

$$d_{ij}=\left\{\begin{array}[]{ll}1&\mbox{ if }\pi_{j}(l_{j})=e_{i}\ \mbox{(i.e., $e_{i}$ is the downstream end of the edge $k_{j}$)},\\ -1&\mbox{ if }\pi_{j}(0)=e_{i}\ \mbox{(i.e., $e_{i}$ is the upstream end of the edge $k_{j}$)},\\ 0&\mbox{ otherwise}\ \mbox{(i.e., $e_{i}$ is not a vertex on the edge $k_{j}$)}.\end{array}\right.$$

$$\begin{array}[]{l}E_{r}=\{e_{i}\in E:\gamma_{i}>1\}\mbox{ (ramification (or interior junction) vertices)},\\ E_{b}=\{e_{i}\in E:\gamma_{i}=1\}\mbox{ (boundary vertices)},\\ E_{u}=\{e_{i}\in E:\gamma_{i}=1,e_{i}\mbox{ is an upstream boundary vertex}\},\\ E_{d}=\{e_{i}\in E:\gamma_{i}=1,e_{i}\mbox{ is a downstream boundary vertex}\},\end{array}$$

$$\begin{array}[]{l}\Omega=G\times[0,T],\quad\Omega_{j}=[0,l_{j}]\times[0,T],\end{array}$$

$$\begin{array}[]{l}\Omega_{p}=(G\setminus E_{b})\times(0,T],\quad\omega_{p}=(G\times\{0\})\cup(E_{b}\times(0,T]).\end{array}$$

$$u_j(e_i,t)=u_j({\pi }_j^{-1}(e_i),t),u_{x_j}(e_i)=\frac{\partial }{\partial x_j}{u}_j({\pi }_j^{-1}(e_i),t),u_{x_j{x}_j}(e_i)=\frac{{\partial }^2}{\partial x_j^2}{u}_j({\pi }_j^{-1}(e_i),t).$$

$$C(G)=\{g:g_j\in C([0,l_j],\mathbb{R}),j\in I_{N-1}\}$$

$${\Vert g\Vert }_{C(G)}=\underset{j\in I_{N-1}}{\max }\underset{x\in [0,l_j]}{\max }|g_j(x)|.$$

$${\Vert g\Vert }_{C^m(G)}=\sum _{\beta =1}^m\Vert g^{(\beta )}{\Vert }_{C(G)}+{\Vert g\Vert }_{C(G)},$$

$$\|g\|_{L^{p}(G)}=\sum_{j=1}^{N-1}\Big{(}\int_{0}^{l_{j}}|g_{j}|^{p}\Big{)}^{1/p}.$$

$$C^{2+\alpha ,1+\frac{\alpha }{2}}(\Omega )=\{u\in C(\Omega ):u_j\in C^{2+\alpha ,1+\frac{\alpha }{2}}({\Omega }_j),\forall j\in I_{N-1}\}$$

$$\Vert u{\Vert }_{C^{2+\alpha ,1+\frac{\alpha }{2}}(\Omega )}=\sum _{j=1}^{N-1}\Vert u_j{\Vert }_{C^{2+\alpha ,1+\frac{\alpha }{2}}({\Omega }_j)}.$$

$$\frac{\partial u}{\partial t}=D\frac{{\partial }^{2}u}{\partial x^2}-v\frac{\partial u}{\partial x}+f(x,u)u,$$

$$\frac{\partial u_j}{\partial t}=D_j\frac{{\partial }^2{u}_j}{\partial x_j^2}-v_j\frac{\partial u_j}{\partial x_j}+f_j(x_j,u_j)u_j,x_j\in (0,l_j),j\in I_{N-1},t>0,$$

$$u_j(x_j,0)=u_j^0(x_j),x_j\in [0,l_j],j\in I_{N-1}.$$

$${\alpha }_{j,1}{u}_j(e_i,t)-{\beta }_{j,1}\frac{\partial u_j}{\partial x_j}(e_i,t)=0\text{mbox}with{\alpha }_{j,1}\ge 0,{\beta }_{j,1}\ge 0,{\alpha }_{j,1}+{\beta }_{j,1}>0,$$

$$\text{mbox}thezero-fluxboundarycondition:\left(D_j\frac{\partial u_j}{\partial x_j}-v_j{u}_j\right)(e_i,t)=0.$$

$${\alpha }_{j,2}{u}_j(e_i,t)+{\beta }_{j,2}\frac{\partial u_j}{\partial x_j}(e_i,t)=0\text{mbox}with{\alpha }_{j,2}\ge 0,{\beta }_{j,2}\ge 0,{\alpha }_{j,2}+{\beta }_{j,2}>0,$$

$$\text{mbox}thefreeflow(orNeumann)condition:\frac{\partial u_j}{\partial x_j}(e_i,t)=0\text{mbox}or$$

$$\text{mbox}thehostile(orDirichlet)condition:u_j(e_i,t)=0.$$

$$u_{i_1}(e_i,t)=u_{i_2}(e_i,t)=\cdots =u_{i_m}(e_i,t),$$

$$\sum _{j=i_1}^{i_m}{d}_{ij}{A}_j{D}_j\frac{\partial u_j}{\partial x_j}(e_i,t)=0,$$

$$\sum _{j=i_1}^{i_m}{d}_{ij}{A}_j{v}_j=0,$$

$$\sum _{j=i_1}^{i_m}{d}_{ij}{A}_j\left(D_j\frac{\partial u_j}{\partial x_j}-v_j{u}_j\right)(e_i,t)=0.$$

$$\begin{array}[]{l}E_{0}=\{e_{i}\in E_{b}:\mbox{ vertices with hostile (or Dirichlet) condition}\},\\ E_{b}\setminus E_{0}:\mbox{ vertices not assigned with hostile condition}.\end{array}$$

$$\text{noindent}(\mathbf{IBVP})(\text{ref}\mathbf{model}),(\text{ref}\mathbf{modelic}),(\text{ref}\mathbf{bcupstreamzf}),(\text{ref}\mathbf{bcdownstreamzf})\text{mbox}\mathbf{and}(\text{ref}\mathbf{interconds}).$$

$$X=\{\phi \in C^1(G):\phi \text{mbox}satisfies(\text{ref}bcupstreamzf)and(\text{ref}bcdownstreamzf)\},$$

$$X_{+}=\{\phi \in X:\phi \ge 0,\phi \not\equiv 0\}$$

$$X^o=\{\phi \in X:\phi >0\text{mbox}onG\setminus E_0,\text{mbox}and{d}_{ij}{\phi }_{x_j}(e_i)<0\text{mbox}if{e}_i\in E_0\}.$$

Full text

Keywords: River network, population persistence, eigenvalue problems, net reproductive rate

1 Introduction

River and stream ecosystems is a key component of the global environmental ecosystems, and it has gained increasing attention of ecologists and environmental scientists. The organisms living in river systems are subjected to the biased flow in the downstream direction. Instream flow needs (IFNs) and “drift paradox” are two important problems for stream ecologists and river managers. The former asks how much stream flow can be changed while still maintaining an intact stream ecology [3, 40], and the latter asks how stream-dwelling organisms can persist without being washed out when continuously subjected to a unidirectional water flow [42, 43, 44, 53]. The problems are challenging due to the complex and dynamic nature of interactions between the stream environment and the biological community. The study of population models in rivers or streams reveals the dependence of spatial population dynamics on environmental and biological factors in rivers or streams, hence it has become an important explanation of the Instream flow needs and “drift paradox”; see e.g., [1, 2, 9, 17, 34, 35, 36, 38, 40, 44, 53, 55, 56, 21, 20, 32, 33, 37, 51, 19].

In mathematical models, rivers and streams have been traditionally treated as a finite or infinite length one-dimensional interval on the real line (see e.g., [53, 40, 21, 20, 37]). When homogeneous river intervals (see e.g., [53, 17]) are recognized as oversimplification of real river systems, more realistic rivers have been approximated by alternating good and bad patches or pool-and-riffle channels (see e.g, [34, 22]), drift and benthic (or storage) zones (see e.g., [18, 44, 24]), or meandering rivers consisting of a main channel and point bars (see e.g., [23, 22]). These later generalizations are still in one-dimensional space or considered as one channel.

Nevertheless, natural river systems are often in a spatial network structure such as dendritic trees. The network topology (i.e., the topological structure of a network), together with other physical and hydrological features in a river network, can greatly influence the spatial distribution of the flow profile (including the flow velocity, water depth, etc.) in branches of the network. Moreover, the population dispersal vectors may be constrained by the network configuration and the flow profile, and species life history traits may depend on varying habitat conditions in the network. As a result, population distribution and persistence in river systems can be significantly affected by the network topology or structure, see e.g., [8, 11, 12, 14, 15, 16, 46]. Then there arise interesting questions such as whether a population can persist in the desert streams of the southwestern United State while the streams are experiencing substantial natural drying trends [11], or whether dendritic geometry enhances dynamic stability of ecological systems [15] etc. Furthermore, other related dynamics in the network, such as the dynamics of water-born infectious diseases like cholera may also be greatly affected by the river network geometry (see e.g., [6]).

Branches in a river network have been modeled as point nodes in a network of habitats in individual based models [11, 16] and matrix population models for stage-structured populations [14, 39]. However this oversimplifies the spatial heterogeneity of river networks. In a real river ecosystem, organisms mainly live in the branches of the network and the connections between branches (e.g., the network nodes) are mainly for population transitions from one branch to another. To take into account this realistic situation, in recent works [48, 49, 50], integro-differential equations and reaction-diffusion-advection equations were used to model population dynamics in river networks where the network branches, instead of the network nodes, are the main habitat for organisms. Here the river networks are modeled under the framework of metric graphs (or metric networks). A metric graph is a graph $G=(V,E)$ with a set $V$ of vertices and a set $E$ of edges, such that each edge $e\in E$ is associated with either a closed bounded interval. Mathematic notion of metric graphs was first introduced in the context of wave propagation on thin graph-like domains [5, 25], and they are also called quantum graphs.

The theories of parabolic and elliptic equations as well as the corresponding eigenvalue problems on metric networks are important in establishing population dynamics of biological species in river networks. The existence and uniqueness of solutions of linear parabolic equations and nonlinear parabolic equations have been established in [57] and [60], respectively. A maximum principle for semilinear parabolic network equations was obtained in [59], and the eigenvalue problems associated with parabolic equations in networks were studied in [58, 61]. Stability of steady states of parabolic equations were studied in [62, 65]. More studies of diffusion equations in networks can be found in e.g., [4, 63, 30, 31]. In these work, the model parameters are allowed to be time and/or space dependent; the so-called Kirchhoff laws or an excitatoric Kirchhoff condition (or dynamical node condition) are assumed at the interior connecting points.

The goal of this paper is to establish a mathematically rigorous foundation of reaction-diffusion-advection equations defined on a metric tree network, which models population dynamics of a biological species on a river network. The population model consists of reaction-diffusion-advection equations describing population dynamics on network branches and zero population flux at interior connecting points in the network, allowing variations of diffusion rates, advection rates, and growth rates throughout the network. We will rigorously derive the theories for the time-dependent parabolic equations, the corresponding elliptic equations for the steady states, and the associated eigenvalue problems for linearized equations, to establish population persistence conditions in terms of the principal eigenvalue and/or the net reproductive rate of the system. We will also prove the existence and uniqueness of a stable positive steady state when population persists under the logistic type growth rate. The theory of infinite-dimensional dynamical systems and existing theories of parabolic and elliptic equations as well as eigenvalue problems on a real line and on metric networks will be applied. We will also study how different factors influence population dynamics, especially persistence and the distribution of the stable positive steady state (if exists) in the whole network.

The population persistence in a spatial population model has been described by uniform persistence, the (in)stability of the trivial extinction solution, and the critical domain size (minimal length of the habitat such that a species can persist); see e.g., [20, 18, 27, 28, 29, 37, 40, 44]. For a single species population in one-dimensional rivers, the persistence theory was established in a homogeneous environment in [53, 55, 37, 32], in temporally periodically varying environments in [20], in temporally randomly varying environments in [19], and in spatially heterogeneous environment in [33]. For a benthic-drift population consisting of individuals drifting in water and individuals staying on the benthos, the critical domain size was studied in a spatially homogeneous river in [44] and in a river with alternating good and bad channels in [34]. In particular, persistence metrics (fundamental niche, source/sink metric, and the net reproductive rate) have been established for a single stage population in [40] and for a benthic-drift population in [18], respectively. Population persistence for a single species in river networks has also been studied in [48, 49, 50]. Integro-differential equations were used to describe population dynamics in river networks in [48], where the diffusion coefficients, advection rates and growth rates were assumed to be the same in all branches, and the population persistence was determined by the stability of the extinction state. Reaction-diffusion-advection equations were used in [49, 50], where again constant diffusion coefficients, advection rates, and growth rates were assumed throughout the network but zero-flux interior junction conditions were not assumed. The principal eigenvalue of the corresponding eigenvalue problem was used to determine the stability of the extinction state and also the population persistence. Most of the analyses and results about persistence conditions were restricted to radial trees, in which all branches on the same level are essentially the same habitats, and hence population dynamics in such networks is essentially equivalent to that in a one-dimensional river.

This paper is organized as follows. In section 2, we introduce the notion of river network of a general tree and the initial boundary value problem for population dynamics on the network. In section 3, we establish the existence of a principal eigenvalue of the corresponding eigenvalue problem and we show that it determines whether the population persists (the extinction solution is unstable) or becomes extinct (the extinction solution is globally asymptotically stable). Moreover, we obtain the existence of a globally asymptotically stable positive steady state when the population persists. In section 4, we define the next generation operator and the net reproductive rate $R_{0}$ of the population living in the river network, and we prove that $R_{0}=1$ can be used as a persistence threshold for the population. We also provide a method to calculate $R_{0}$. In section 5, by virtue of numerical simulations, we study the influences of hydrological, physical, and biological factors on the net reproductive rate as well as the positive steady state. In Appendix A, we provide the derivation of the theories for the parabolic and elliptic problems on networks, including the maximal principle, the comparison principle, and the existence, uniqueness and estimations of the solutions.

2 Model

2.1 The river network - a metric tree

In this work, we assume that the river network is a finite metric tree, i.e., a connected finite metric graph with no cycles, or equivalently, a finite metric graph on which any two vertices can be connected by a unique simple path.

We first introduce the mathematical definition of a river network (a finite tree) and notations on it (see e.g., [60]). Let $G$ be a $C^{\kappa}$-network for $\kappa\geq 2$ with the set of vertices

[TABLE]

the set of edges

[TABLE]

and arc length parameterization $\pi_{j}\in C^{\kappa}([0,l_{j}],\mathbb{R}^{2})$ on edge $k_{j}$, where $N$ and $N-1$ are the numbers of vertices and edges, respectively, and

[TABLE]

The edge $k_{j}$ is isomorphic to the interval $[0,l_{j}]$ with length $l_{j}$ and spatial variable $x_{j}$ on it, where $x_{j}=0$ and $x_{j}=l_{j}$ represent the upstream end and the downstream end of $k_{j}$, respectively. The topological graph $\Lambda=(E,K)$ embedded in $G$ is assumed to be simple and connected. Thus, $\Lambda$ admits the following properties: each $k_{j}$ has its endpoints in $E$, any two vertices in $E$ can be connected by a unique simple path with arcs in $K$, and any two distinct edges $k_{j}$ and $k_{h}$ intersect at no more than one point in $E$. See Figure 2.1 for an example of a river network.

Endowed with the above graph topology and metric defined on each edge, $G$ is a connected and compact subset of $\mathbb{R}^{2}$. The orientation of $G$ is given by the incidence matrix $(d_{ij})_{N\times(N-1)}$ with

[TABLE]

We distinguish the set $E$ of vertices as follows:

[TABLE]

where $\gamma_{i}=\gamma(e_{i})$ is the valency of $e_{i}$ that represents the number of edges that connect to $e_{i}$, and $E_{b}=E_{u}\cup E_{d}$.

Let $t$ be the time variable and for $T>0$, denote

[TABLE]

[TABLE]

For a function $u:\Omega\rightarrow\mathbb{R}$, we define $u_{j}=u\circ(\pi_{j},id):\Omega_{j}\rightarrow\mathbb{R}$. Differentiation is carried out on each edge $k_{j}$ with respect to the arc length parameter $x_{j}$. A function is differentiable on $G$ means that it is differentiable at all points $x\in G\setminus E$. We use the following notations for functions and derivatives at a vertex

[TABLE]

Any function $u$ on $\Omega$ satisfies $u_{j}(e_{i},t)=u_{h}(e_{i},t)\mbox{ if }k_{j}\cap k_{h}=\{e_{i}\}.$

We now introduce function spaces on $G$. Let

[TABLE]

with the norm:

[TABLE]

The Banach space $C^{m}(G)$ consists of all functions that are $m$ times continuously differentiable over $G$ with norm given by

[TABLE]

where $g^{(\beta)}$ is the $\beta$-th derivative of $g$. Similarly $L^{p}(G)$ is the Banach space of all real-valued functions defined on $G$ that are measurable and $p$-summable with respect to $G$ with $p\geq 1$. The norm in this space is defined by

[TABLE]

For $\alpha\in[0,1)$, define

[TABLE]

where $C^{2+\alpha,1+\frac{\alpha}{2}}(\Omega_{j})$ with the usual norm $\|\cdot\|_{C^{2+\alpha,1+\frac{\alpha}{2}}(\Omega_{j})}$ denotes the Banach space of functions $u$ on $\Omega_{j}$ having continuous derivatives $\frac{\partial^{r+s}u}{\partial t^{r}\partial x_{j}^{s}}$ for $2r+s\leq 2$ and finite Hölder constraints of the indicated exponents in the case of $\alpha>0$. Then $C^{2+\alpha,1+\frac{\alpha}{2}}(\Omega)$ is a Banach space endowed with the norm

[TABLE]

Similarly we can define $C^{2+\alpha}(G)$, $W^{2}_{p}(G)$ and $W^{2,1}_{p}(\Omega)$ for any fixed $\alpha\in[0,1)$ and $p\geq 1$.

2.2 The population model in the river network

Since Speirs and Gurney’s work [53], the dynamics of a population living in a one-dimensional river has been described by the following reaction-diffusion-advection equation:

[TABLE]

where $u(x,t)$ is the population density at location $x$ and time $t$, $D$ is the diffusion coefficient, $v$ is the flow velocity, and $f$ is the per capita growth rate.

We adapt model (2.2) to a population living in a river network $G$. The dynamics of the population can be described by

[TABLE]

where $u_{j}$ is the population density on the edge $k_{j}$, $D_{j}$ is the diffusion coefficient on $k_{j}$, $v_{j}$ is the flow velocity on $k_{j}$, and $f_{j}$ is the per capita growth rate on $k_{j}$. The initial population distribution in $G$ is $u^{0}$, that is,

[TABLE]

There are three types of vertices in the river network $G$: upstream boundary ends, downstream boundary ends, and interior junction vertices. Correspondingly, boundary or interface conditions are imposed at each vertex of $E$.

• •

At an upstream boundary point $e_{i}\in E_{u}$ that only connects to the edge $k_{j}$, the boundary condition can be assumed as

[TABLE]

for instance,

[TABLE]

• •

At a downstream boundary point $e_{i}\in E_{d}$ that only connects to the edge $k_{j}$, the boundary condition can be assumed as

[TABLE]

for instance,

[TABLE]

[TABLE]

• •

At an interior junction point $e_{i}\in E_{r}$, the population density is continuous and the total population flux in and out is zero. Hence, the connection conditions are the continuity conditions and Kirchhoff laws:

[TABLE]

where $e_{i}\in E_{r}$ connects to edges $k_{i_{1}}$, $k_{i_{2}}$, $\cdots$, and $k_{i_{m}}$, $d_{ij}$ is defined in (2.1), $A_{j}$ is the wetted cross-sectional area of the edge $k_{j}$, and (2.10b) is the result of substituting the continuity condition (2.10a) and the conservation of the flow at $e_{i}$

[TABLE]

into the zero-flux condition at $e_{i}$

[TABLE]

According to different ecological conditions at boundary vertices on $G$, we further use the following notations throughout the paper:

[TABLE]

We finally define an initial boundary value problem for a population in a river network:

[TABLE]

Furthermore, we impose the following assumptions in different parts of the paper.

1. [H1]

For each $j\in I_{N-1}$, $D_{j}>0$, $v_{j}\geq 0$, $A_{j}>0$. 2. [H2]

For each $j\in I_{N-1}$, $f_{j}:[0,l_{j}]\times[0,\infty)\rightarrow\mathbb{R}$ is continuous and there exists a constant $M_{j}>0$ such that $f_{j}(x,u)\leq 0$ for any $x\in[0,l_{j}]$ and $u\geq M_{j}$, and $f_{j}(\cdot,u_{j})u_{j}$ is Lipschitz continuous in $u_{j}$ with Lipschitz constant $L_{j}>0$. 3. [H3]

For each $j\in I_{N-1}$, $f_{j}(\cdot,u_{j})$ is monotonically decreasing in $u_{j}$.

By adapting theories for parabolic and elliptic equations on intervals and/or networks [41, 4, 13, 47, 57, 59, 60, 26, 52, 45, 66], we develop the fundamental theories of parabolic and elliptic problems on networks corresponding to (IBVP); see Appendix A. In particular, for linear parabolic problems, we establish the strong maximum principle (in Lemma A.1), Hopf boundary lemma for networks (in Lemma A.2), comparison principle (in Lemma A.3), and the existence, uniqueness, $L^{p}$ and Schauder estimates of solutions (in Theorem A.4, via writing the differential operator into a self-adjoint operator on the network); for the nonlinear problem (IBVP), we develop the theory of the existence, uniqueness and positivity of solutions (in Theorem A.7, by using the upper and lower solutions) and prove the monotonicity and strict subhomogeneity of the solution map (in Lemmas A.8 and A.9, respectively); for the corresponding elliptic problems, we also obtain the strong maximum principle (in Lemma A.10), Hopf boundary lemma (in Lemma A.11), comparison principle (in Lemma A.12), and the existence, uniqueness, $L^{p}$ and Schauder estimates of solutions (in Theorems A.13 and A.15). These mathematical preparations enable us to establish the extinction/persistence criteria for system (IBVP).

3 The eigenvalue problem and population persistence

In this section, we consider the eigenvalue problem corresponding to the linearized system of (IBVP) at the trivial solution, obtain the existence of the principal eigenvalue, and then use the principal eigenvalue as a threshold for population persistence and extinction. We also obtain the existence, uniqueness and stability of a positive steady state when the population persists. Assumptions [H1]-[H3] are all imposed throughout the rest of the paper.

3.1 The eigenvalue problem and its principal eigenvalue

We first introduce some Banach spaces which will be used frequently later. Denote

[TABLE]

and let

[TABLE]

be the positive cone in $X$. The interior of $X_{+}$ is

[TABLE]

Then $X_{+}$ is a solid cone of $X$ with nonempty interior $X^{o}$. We also write $\varphi_{1}\gg\varphi_{2}\ \ \mbox{if}\ \varphi_{1}-\varphi_{2}\in X^{o}.$

The linearization of (IBVP) at the trivial solution $u=0$ is

[TABLE]

Substituting $u_{j}(x_{j},t)=e^{\lambda t}\psi_{j}(x_{j})$ into (3.4), we obtain the corresponding eigenvalue problem

[TABLE]

For simplicity, denote $\mathcal{L}$ to be the operator such that $\mathcal{L}|_{k_{j}}=\mathcal{L}_{j}$, where

[TABLE]

The following result indicates that (3.5) admits a simple eigenvalue associated with a positive eigenfunction. The proof is given in Appendix B.

Proposition 3.1**.**

The eigenvalue problem (3.5) admits a simple eigenvalue $\lambda^{\ast}$ associated with a positive eigenfunction $\psi^{\ast}\in X^{o}$. None of the other eigenvalues of (3.5) corresponds to a positive eigenfunction; and if $\lambda\neq\lambda^{\ast}$ is an eigenvalue of (3.5), then $\mathrm{Re}(\lambda)\leq\lambda^{\ast}.$

Let $\lambda^{\ast}$ be the eigenvalue of the eigenvalue problem (3.5) with a corresponding positive eigenfunction $\psi^{\ast}\in X^{o}$. We call $\lambda^{\ast}$ the *principal eigenvalue * of (3.5).

We say that $\mathcal{L}$ has the strong maximum principle property if $u\in C^{2}(G)$ satisfying

[TABLE]

implies that $u>0$ in $G\backslash E_{0}$ unless $u\equiv 0$. We also say $u\in C^{2}(G)$ is an upper solution of $\mathcal{L}$ if (3.7) holds, and such $u$ is called as a strict upper solution of $\mathcal{L}$ if it is an upper solution but is not a solution. Then the analysis of [10, Theorem 2.4] can be easily adapted to conclude the following result.

Proposition 3.2**.**

The following statements are equivalent.

1. (i)

$\mathcal{L}$* has the strong maximum principle property;* 2. (ii)

$\mathcal{L}$* has a strict upper solution which is positive in $G\backslash E_{0}$;* 3. (iii)

$\lambda^{\ast}<0$.

3.2 Persistence and extinction

We now use the sign of the principal eigenvalue of (3.5) to determine the population persistence or extinction as well as the existence of a positive steady state for (IBVP), which satisfies the following elliptic equations:

[TABLE]

The following result shows that the principal eigenvalue $\lambda^{\ast}$ is the key threshold of extinction/persistence for (IBVP). The proof is given in Appendix C.

Theorem 3.3**.**

Let $\lambda^{\ast}$ be the principal eigenvalue of the eigenvalue problem (3.5) with corresponding eigenfunction $\psi^{\ast}\in X^{o}$. Then

1. (i)

If $\lambda^{\ast}\leq 0$, then $u\equiv 0$ is globally attractive for (IBVP) for all initial values in $X_{+}$. 2. (ii)

If $\lambda^{\ast}>0$, then (IBVP) admits a unique positive steady state $u^{\ast}\in X^{o}$ which is globally attractive for all initial values in $X_{+}\setminus\{0\}$.

The above theorem indicates that $\lambda^{\ast}=0$ is a threshold to determine population persistence on a river network. The idea has been used in [48, 49, 50] but this is the first time that it is rigorously proved for population models on river networks.

4 The net reproductive rate $\mathcal{R}_{0}$

The net reproductive rate has been defined and proved to be a threshold quantity for population persistence in a single river channel [40, 18]. In this section, we will define the next generation operator and the net reproductive rate $\mathcal{R}_{0}$ for (IBVP) and then use $\mathcal{R}_{0}$ to determine the population persistence and extinction. Moreover, we will provide a numerical method to calculate $\mathcal{R}_{0}$.

4.1 Definition of the net reproductive rate $\mathcal{R}_{0}$

Assume the growth rate of the population on edge $k_{j}$ satisfies $f_{j}(x_{j},u_{j})u_{j}=\tilde{f}_{j}(x_{j},u_{j})u_{j}-m_{j}(x_{j})u_{j}$, where $\tilde{f}_{j}$ is the recruitment rate and $m_{j}(x_{j})$ is the mortality rate. Let $r_{j}(x_{j})=\tilde{f}_{j}(x_{j},0)$ and assume $r,m\in C(G)$. Then

[TABLE]

For $\phi^{0}\in X$, assume that $\phi$ satisfies

[TABLE]

Define $\Gamma:X\rightarrow X$ by

[TABLE]

where $\phi$ is the solution of (4.1) with initial condition $\phi^{0}$. That is, $\Gamma$ is a linear operator mapping an initial distribution of the population to its offspring distribution. Hence, we call $\Gamma$ the next generation operator. Let

[TABLE]

where $r(\Gamma)$ is the spectral radius of the linear operator $\Gamma$ on $X$. Then $\mathcal{R}_{0}$ represents the average number of offsprings that an individual produces during its lifetime and we call $\mathcal{R}_{0}$ the net reproductive rate.

Let $B:\mathcal{D}\rightarrow\mathcal{D}$ with

[TABLE]

be defined by

[TABLE]

Let $\Gamma_{1}:X\rightarrow X$ be such that

[TABLE]

where $\phi$ is the solution of (4.1) with initial condition $\phi^{0}$. Similarly as in Proposition 2.10 of [40], we can prove that $-\Gamma_{1}$ is the inverse operator of $B$, i.e., $B^{-1}=-\Gamma_{1}$. Hence, $\Gamma(\phi)=-QB^{-1}(\phi)$, where the operator $Q$ is defined as

[TABLE]

Then $\mathcal{L}=B+Q$. By Proposition 3.1 and the above analysis, noting that $(\lambda I-B)^{-1}$ is defined such that

[TABLE]

where $\phi$ is the solution of (4.1) with initial condition $\phi^{0}$, we know that both $\mathcal{L}$ and $B$ are resolvent-positive operators in $X$. It follows from Propositions 3.1 and 3.2 that the spectral bound $s(B)$ of $B$ is the principal eigenvalue of $B$ and $s(B)<0$. We then obtain the following result by using [54, Theorem 3.5].

Lemma 4.1**.**

$\mathcal{R}_{0}-1$* and $\lambda^{\ast}$ have the same sign, where $\lambda^{\ast}$ is the principal eigenvalue of the eigenvalue problem (3.5).*

Theorem 3.3 and Lemma 4.1 imply the following result.

Corollary 4.2**.**

If $\mathcal{R}_{0}\leq 1$, then $u\equiv 0$ is globally attractive for (IBVP); if $\mathcal{R}_{0}>1$, then (IBVP) admits a unique positive steady state $u^{\ast}\in X^{o}$, which is globally attractive for all initial values in $X_{+}\setminus\{0\}$.

Therefore, $\mathcal{R}_{0}=1$ is the threshold for population persistence and extinction. The population will be extinct if $\mathcal{R}_{0}\leq 1$ and it is persistent if $\mathcal{R}_{0}>1$.

4.2 Calculation of $\mathcal{R}_{0}$

Let

[TABLE]

Integrating (4.1) with respect to $t$ from [math] to $\infty$ yields

[TABLE]

Note that $\phi_{j}(\cdot,t)\rightarrow 0$ as $t\rightarrow\infty$. The above equation implies that

[TABLE]

Therefore, $\Gamma_{1}(\phi^{0})$ is the solution of

[TABLE]

We define

[TABLE]

where $u$ is the solution of (4.4). Then $T_{1}$ is a compact and strongly positive operator on $X$. By the definition of $T_{1}$ and equations (4.3) and (4.4), we know that $\Gamma_{1}=T_{1}$ on $X$. Hence, $\Gamma=Q\Gamma_{1}$ is also a compact and strongly positive operator on $X$. It follows from [10, Theorem 1.2] that $\mathcal{R}_{0}=r(\Gamma)>0$ is a simple eigenvalue of $\Gamma$ with an eigenfunction $\phi^{\ast}\in X^{o}$, i.e.,

[TABLE]

and there is no other eigenvalues of $\Gamma$ associated with positive eigenfunctions.

By following the idea in the proof of [64, Theorem 3.2], we can obtain $\mathcal{R}_{0}$ via the principal eigenvalue of another eigenvalue problem.

Theorem 4.3**.**

If the eigenvalue problem

[TABLE]

admits a unique positive eigenvalue $\mu^{0}$ with a positive eigenfunction, then $\mathcal{R}_{0}=1/\mu^{0}$.

To numerically calculate $\mathcal{R}_{0}$, we use the finite difference method to discretize (4.5) and approximate (4.5) by

[TABLE]

where $\hat{r}$ is a diagonal matrix containing values of $r(x_{j})$ on the main diagonal and $\hat{A}$ is the discretization of the operator on the right-hand side of (4.5). The matrix $\hat{A}^{-1}\hat{r}$ is a non-negative and irreducible matrix. The Perron-Frobenius Theorem implies that $\hat{A}^{-1}\hat{r}$ admits a principal eigenvalue $\varsigma^{\ast}$, which is the unique simple eigenvalue of $\hat{A}^{-1}\hat{r}$ associated with a positive eigenvector $\psi^{*}$, that is,

[TABLE]

Then we approximate $1/\mu^{0}$ by $\varsigma^{\ast}$, i.e.,

[TABLE]

by using Theorem 4.3 and above approximating scheme.

It follows from (4.6) that $\hat{r}\psi^{\ast}$ is the eigenvector of $\hat{r}\hat{A}^{-1}$ corresponding to $\varsigma^{\ast}$, i.e., $\hat{r}\hat{A}^{-1}(\hat{r}\psi^{\ast})=\varsigma^{\ast}(\hat{r}\psi^{\ast})$. Note that the next generation operator $\Gamma$ can be approximated by $\hat{r}\hat{A}^{-1}$ and that $\mathcal{R}_{0}=1/\mu^{0}$. The eigenvalue problem $\Gamma\phi^{\ast}=\mathcal{R}_{0}\phi^{\ast}$ can be approximated by

[TABLE]

Then we obtain

[TABLE]

That is, $\hat{r}\psi^{\ast}$ can be used to approximate the eigenfunction of $\Gamma$ associated with the eigenvalue $\mathcal{R}_{0}$. We call $\phi^{\ast}$ (or $\hat{r}\psi^{\ast}$ as the approximation) the next generation distribution of the population.

5 The influences of factors on population persistence

The results in Sections 3 and 4 show that both the principal eigenvalue $\lambda^{\ast}$ of the eigenvalue problem (3.5) and the net reproductive rate $\mathcal{R}_{0}$ can be used to determine the population persistence. For the biological significance of $\mathcal{R}_{0}$, now we apply the theory in Section 4 to investigate the influences of biotic and abiotic factors on population dynamics (in particular, persistence or extinction) of (IBVP) via numerical studies of $\mathcal{R}_{0}$ and the stable positive steady state (if exists).

Real river networks are complex and the quantitative influence of a factor on the population persistence highly depends on the structure and scales of the network. While the choice of a general river network is random, we consider a few simple but typical river networks of trees with one, three, four, five, and seven branches, representing different types of network topologies, merging from the upstream or splitting into the downstream; see Figure 5.1. In particular, river networks (3-a), (3-b), (7-a) and (7-b) are radial trees, which are rooted trees with all tree features, including edge lengths, sectional areas, and boundary conditions depending only on the distance to the root (see Section of [50]).

Three sets of boundary conditions for (IBVP) are considered:

1. (ZF-FF)

Zero-flux condition (2.6) at the upstream and free flow (i.e., Neumann) condition (2.8) at the downstream; 2. (ZF-H)

Zero-flux condition (2.6) at the upstream and hostile condition (2.9) at the downstream; 3. (H-H)

Hostile condition (2.9) at both the upstream and the downstream.

We adopt the baseline parameters in [53] and vary their values to see the influences of different factors. The units of parameters are given in Table 5.1. Note that for simplicity, we choose a constant growth rate $r_{j}$ on each edge $k_{j}$, which may result in discontinuity of the growth rate $r(x)$ at the interior vertices of the network. As this assumption can be considered as an approximation of a continuous growth rate in the network, it does not change the essence of our results.

5.1 The influence of the river network structure on population persistence

Natural rivers are rarely in form of single branches but in various types of networks. The structure (or topology) of the river network influences hydrodynamics in the network as well as the intrinsic ecosystem dynamics.

To see how the network structure influences population persistence, we compare the values of $\mathcal{R}_{0}$ in a single branch river and in all river networks in Figure 5.1. Suppose that the growth rate and the diffusion rate do not change throughout each network. We vary the flow conditions by fixing the advection rates but varying the cross-sectional areas or fixing the cross-sectional areas but varying the advection rates, according to the conservation relation (2.11) of the flow at the interior junctions.

5.1.1 The influence of the total length of the network

In a specific network structure, assume that the lengths of all branches are the same and vary the total length of the network. Figure 5.2 shows that when the total length increases, the net reproductive rate $\mathcal{R}_{0}$ increases for all the networks in consideration. This coincides with the well-known result in one-dimensional river (see e.g., [20, 37]): given the same habitat conditions, increasing the total habitat size helps population persistence.

Figure 5.2 also shows that in radial trees ((3-a), (3-b), (7-a), and (7-b) in Figure 5.1), when the total length of the network is fixed, increasing the number of levels reduces the value of $\mathcal{R}_{0}$. It has been shown in [50] that population dynamics on a radial tree is equivalent to that in a one-dimensional river formed from one upstream to one downstream of the tree. Our observation coincides with this result in [50] since a radial tree of $3$ branches with total length $L$ is equivalent to a one-dimensional river of length $2L/3$, while a radial tree of $7$ branches with total length $L$ is equivalent to a one-dimensional river of length $3L/7$. However, one cannot conclude a general result from this that in river networks, when the total river length is fixed, the networks with more branches have smaller net reproductive rates. In the non-radial networks (4-a,b) and (5-a,b), when the flow velocity is fixed, the net reproductive rates of networks (5-a) and (5-b) are larger than those of networks (4-a) and (4-b); see Figure 5.2 (a,c). Networks with $5$ branches also have larger $\mathcal{R}_{0}$ than those with $4$ branches when the cross-sectional areas are fixed and the total river lengths are small; see Figure 5.2 (b,d).

5.1.2 The influence of boundary conditions

Figure 5.2 shows the relationship between the total river network length $L$ and the net reproductive rate $\mathcal{R}_{0}$ under different boundary conditions. When (ZF-FF) or (H-H) boundary conditions are applied, the same number of branches merging from the upstream or splitting into the downstream lead to the same net reproductive rate in radial trees (3-a,b) and (7-a,b) as well as in (5-a,b). That is, when both boundary conditions are bad (hostile) or not bad (not hostile), more branches in the upstream or in the downstream does not change the net reproductive rate, provided that two branches merge into one or one branch splits into two at each junction point. Nevertheless, when (ZF-H) boundary conditions are applied, if the total branch numbers are the same, the net reproductive rate of the network merging from the upstream is not less than the one of the network splitting into the downstream, in (3-a,b), (5-a,b) or (7-a,b). That is, when the downstream boundary condition is bad (hostile), more branches in the downstream cannot result in a larger net reproductive rate, provided that two branches merge into or are split from one branch at each junction point.

Among all the networks, the $4$-branch networks are exceptional. When (ZF-FF) or (H-H) boundary conditions are applied, $4$ branches merging from the upstream result in a larger $\mathcal{R}_{0}$ than splitting into the downstream. Nevertheless, when (ZF-H) boundary conditions are applied and the total river length is small, the network with $4$ branches merging from the upstream may have a smaller net reproductive rate than the one with $4$ branches splitting into the downstream if the advection rate is fixed.

Overall, in all these types of river networks with equal branch length, having more upstream branches helps the population persistence or at least does not accelerate the population extinction, provided that the upstream ends are not the only boundaries that are subjected to hostile conditions and that the total network length is sufficiently large.

To see the effect of boundary conditions on the net reproductive rate more closely, we focus on the $R_{0}$ values in networks with $3$ branches in Figure 5.3. It shows that $\mathcal{R}_{0}$ under zero-flux condition at the upstream and free flow condition at the downstream is the largest and the one under hostile condition at both ends is the smallest. When the total length of the network is small, $\mathcal{R}_{0}$ under hostile boundary conditions is much lower since individuals are close to the hostile boundary and are subjected heavy stress of being removed from the system. When the total length of the network is large, the hostile condition at the downstream does not make much difference on $\mathcal{R}_{0}$ compared to the free flow condition at the downstream if more branches are at the upstream, and the upstream conditions do not influence $\mathcal{R}_{0}$ very much if more branches are at the downstream with hostile condition. Therefore, both the boundary conditions and the type of river network affect the net reproductive rate.

5.2 Population persistence in uniform flows

Hydrological and physical factors in a real river are closely related to each other (see e.g., [7, 22]). To better see how the population persistence is influenced by physical, hydrological and biological conditions, we incorporate the explicit relation between hydrological and physical factors into the population model and investigate their effect on the net reproductive rate and the positive steady state (if exists).

Assume that in the network each branch $k_{j}$ has a constant bottom slope $S_{0j}$, a constant bottom Manning roughness $n_{j}$, and rectangular cross sections with a constant width $B_{j}$. We further assume that the water flow is at the steady state with flow discharge $Q_{j}$ in $k_{j}$, hence, there is a uniform flow in $k_{j}$ [7]. As a result, the water depth $y_{j}$ in $k_{j}$ can be estimated by the normal depth defined in (D.2) in Appendix D, that is,

[TABLE]

which yields the flow velocity in branch $k_{j}$ as

[TABLE]

We then substitute (5.2) into (2.3) to study how parameters influence population persistence in uniform flows in river networks where two or three branches merge into one branch (i.e., (3-a) and (4-a) in Figure 5.1).

5.2.1 The influence of the flow discharges on $\mathcal{R}_{0}$

Firstly we consider two isolated river branches, a longer and wider one with length $l_{1}=1600$ and width $B_{1}=20$, a shorter and narrower one with length $l_{2}=800$ and width $B_{2}=4$, both with the same bottom slope $S_{0j}=0.000001$ and bottom roughness $n_{j}=0.2$. Figure 5.4(a) shows the relationship between the net reproductive rate $\mathcal{R}_{0}$ and the flow discharge in each river, when the diffusion rate, birth and death rates are the same in both rivers and (ZF-FF) boundary conditions are applied. $\mathcal{R}_{0}$ decreases in both rivers as the upstream flow discharge $Q$ increases, but it is larger in the longer and wider river than in the shorter and narrower river.

Now we consider a river network of type (3-a), which is the result of the above small river merging into the large river at the midpoint of the large river. We use subscript 1 for parameters in the large river branch and subscript 2 for parameters in the small river branch. Figure 5.4(b) shows the dependence of $\mathcal{R}_{0}$ on the flow discharges $Q_{1}$ and $Q_{2}$ at the upstream of both river branches. Clearly, $\mathcal{R}_{0}$ is large when the upstream flow discharges $Q_{1}$ and $Q_{2}$ are both small and $\mathcal{R}_{0}$ is small when both $Q_{1}$ and $Q_{2}$ are large. Hence, in a river network, it is still true that low upstream flow discharges help the population persistence and high upstream flow discharges accelerate the population extinction. Moreover, for the parameters we choose, varying $Q_{1}$ changes $\mathcal{R}_{0}$ more than varying $Q_{2}$. That is, given the same habitat and demography conditions, in a merging river network, the upstream flow discharge in the large river influences the global population persistence/extinction more than the upstream flow discharge in the small river. Note that in our case, the population will be extinct in the small river if isolated (see the red curve in Figure 5.4(a)), but merging the rivers into a network helps the population to persist in the whole network, provided that the population can persist in the isolated large river.

5.2.2 The influence of the flow discharges and the widths on $\mathcal{R}_{0}$

We continue to consider river networks as the result of merging large and small rivers. In networks of type (3-a), assume that the large river is given and that the conditions in the small upstream branch vary. Figure 5.5 shows the dependence of $\mathcal{R}_{0}$ on the upstream flow discharge $Q_{2}$ and the width $B_{2}$ of the small river branch, in the cases where the population can persist in the isolated large river (Figure 5.5(a)) and where the population will be extinct in the isolated large river (Figure 5.5(b)). Both figures show the same phenomenon: to increase $\mathcal{R}_{0}$ or help the population persistence in the whole network, it is necessary to have small upstream flow discharge and large width in the small river, which essentially means that the flow discharge per unit width $Q_{2}/B_{2}$ should be low in the small river. If the persistence conditions in the large river become worse (e.g., from Figure 5.5(a) to Figure 5.5(b)), then lower $Q_{2}/B_{2}$ in the small branch is required to help the population persistence, i.e., in order for $\mathcal{R}_{0}>1$, in the whole network.

In the case where two small rivers with identical conditions merge into a large river and result in a network of type (4-a), the dependence of $\mathcal{R}_{0}$ and the flow discharge and the width of the small rivers is shown in Figure 5.6. Similarly, low upstream discharge and large width in the small rivers help the population persistence in the whole network. Comparing Figure 5.5(a) and Figure 5.6, we see that when more small rivers merge into a large river, lower flow discharge per unit width $Q_{2}/B_{2}$ is needed to help the population persistence, i.e., for $\mathcal{R}_{0}$ to go beyond $1$. However, once the conditions are good enough for population to persist, $\mathcal{R}_{0}$ is larger in the larger network (i.e., (4-a)) than that in the smaller network (i.e., (3-a)). In general, this indicates that the threshold conditions for population persistence in a large network may be stronger than those in a small network, but population grows better in the large network once it persists.

5.2.3 The influence of the flow discharge and biological conditions on $\mathcal{R}_{0}$

We continue to consider the networks of (3-a) and (4-a), where one or two identical small rivers merge into a large river. We vary the flow discharge $Q_{2}$, the diffusion rate $D_{2}$, and the birth rate $f_{2}$ in the small rivers to see the co-influence of flow conditions and biological conditions on population persistence and extinction. Figure 5.7 shows the regions for population persistence or extinction in the $D_{2}/D_{1}$-$f_{2}/f_{1}$ plane under different flow conditions. When the conditions in the large river are given, the larger the diffusion rate $D_{2}$ or the birth rate $f_{2}$ is, or the smaller discharge $Q_{2}$ is in the small river branch, the easier it is for the population to persist ($\mathcal{R}_{0}>1$) in the whole network. Comparing the two panels of Figure 5.7, we see that the parameter region for $\mathcal{R}_{0}>1$ is smaller in the network (4-a) than that in the network (3-a), under the same flow conditions. This confirms our previous observation that merging more small rivers with the same conditions into a large river leads to stronger threshold conditions for population persistence in the whole river network.

5.2.4 Good or bad regions in a river network

To see which part of a river network is good or bad for population persistence, we plot the spatial profile of the stable positive steady state (when $\mathcal{R}_{0}>1$) (see Figure 5.8) for the network (3-a), under different hydrological and biological conditions and three sets of boundary conditions. Our simulations show that when $\mathcal{R}_{0}$ in the whole network is large, the positive steady state is large in the upstream branch with better conditions (for persistence) and in the downstream branch. These branches are hence overall good regions for population growth and persistence. When $\mathcal{R}_{0}$ is less than $1$ or larger than but close to $1$, the positive steady state (if exists) is mainly distributed in the downstream branch, and hence only the downstream branch is the relatively better region in the whole network and the upstream (large and small) branches are considered to be bad regions for population growth or persistence. Figure 5.8 shows an example that the population can persist in the large river if isolated and the diffusion rate and the birth rate in the small river branch are larger than those in the large river branches. When the flow discharge in the small branch is much lower than that in the large branch, the resulted $\mathcal{R}_{0}$ is large and good regions are distributed in the middle part of the small branch and the downstream branch; when the flow discharge in the small river is not very low, the resulted $\mathcal{R}_{0}$ is less than $1$ or larger but close to $1$ and the downstream branch is relatively good. If hostile condition is applied at a boundary, then the best region (where the steady state takes the maximum value) shifts from the region near the boundary to the middle of the branch. In all the cases, it seems that the downstream branch is at least not the worst region in a merging river network. We can also plot the next generation distribution in these cases and the distributions are similar to the steady states in Figure 5.8.

6 Discussion

Real river systems are represented by networks. Previously, integro-differential equations and reaction-diffusion-advection equations were used to describe population dynamics in river networks [48, 49, 50], but parameters were assumed to be constants throughout the network and persistence theories in [49, 50] were mainly for radially symmetric trees.

We consider reaction-diffusion-advection equations in networks of general trees and establish the fundamental theories for the parabolic and elliptic problems in the networks (such as maximum principle, comparison principle, the existence, uniqueness and estimation of solutions, etc.). We also obtain the existence of the principal eigenvalue of the corresponding eigenvalue problem and establish the theory for population persistence and existence of a stable positive steady state when population persistence. Furthermore, we define the net reproductive rate $\mathcal{R}_{0}$, which biologically represents the average number of offsprings that a single individual produces during its lifetime, and prove that the population persists in the whole river network if $\mathcal{R}_{0}\geq 1$ but will be extinct if $\mathcal{R}_{0}<1$. The theories were rigorously developed and proved. We convert the calculation of the net reproductive rate $\mathcal{R}_{0}$ to that of the principal eigenvalue of a generalized eigenvalue problem. Then we use the finite difference method to discretize the eigenvalue problem and eventually use the principal eigenvalue of the eigenvalue problem for a matrix to approximate $\mathcal{R}_{0}$. We then use a vector related to the eigenvector associated with the matrix eigenvalue problem to approximate the eigenfunction of the next generation operator, which we call the next generation distribution. We also use the finite difference method to discretize the corresponding elliptic equations to obtain the positive steady state when $\mathcal{R}_{0}>1$.

We conduct numerical simulations in a few simple networks and investigate the influences of hydrological, physical, and biological factors on population persistence in terms of the net reproductive rate and on the spatial profile of the positive steady state (when exists). Our results coincide with existing ones in one-dimensional rivers: in a given type of river network, the population persistence becomes easier if the total river length is larger, the flow discharge/advection is lower, or the diffusion rate or the growth rate is higher. In particular, we compare population dynamics in a smaller network and in a larger network. It turns out that stronger (or better) conditions are required for a population to persist in a larger river network. However, if the population can persist in both networks, then the net reproductive rate is larger in the larger network than in the smaller network. This special case can be considered as a simplification of the real problem of adding or removing one branch from the network, which may happen when one upstream branch in the dessert suddenly disappears because of drought or when human beings plan to construct a new upstream or downstream branch to a network for some economic or other reasons. Our results answer the following questions: how will such phenomenon or activity influence population dynamics in the network, will the loss of a branch in the dessert cause the extinction of a species, and will a new branch attached to the network help the persistence of a species.

By using the interior connecting condition (2.10b), the differential operator in (2.3) in a tree can be proved to be self-adjoint by rearranging the nodes in the tree and the eigenvalue problem can be written into a Sturm-Liouville problem. Hence the theory for the principal eigenvalue and associated eigenfunction can be established. If the network is not a tree, then we cannot establish the principal eigenvalue with this method and the persistence theory still remains open in this situation. The dynamics of interacting species in river networks is also an interesting future problem.

Acknowledgments

Y.J. gratefully acknowledges NSF grant DMS-1411703. R.P. is supported by the NSF of China (No. 11671175, 11571200), the Priority Academic Program Development of Jiangsu Higher Education Institutions, Top-notch Academic Programs Project of Jiangsu Higher Education Institutions (No. PPZY2015A013) and Qing Lan Project of Jiangsu Province. J.S. is partially supported by NSF grant DMS-1715651.

Appendix A Theory of parabolic and elliptic equations on networks

A.1 The strong maximum principle and comparison principle

In this subsection, we establish the strong maximum principle and comparison principle for parabolic equations on metric graphs, which are fundamental in the investigation of existence, uniqueness and positivity of solutions to the nonlinear problem (IBVP).

The following strong maximum principle is an analogue of the classical one for open subsets in Euclidian spaces.

Lemma A.1**.**

Assume that $c(x,t)\geq 0$ and is bounded from above on $\Omega$. Let $u\in C(\Omega)\cap C^{2,1}(\Omega_{p})$ satisfy

[TABLE]

and

[TABLE]

Suppose that $u\leq M\ (u\geq m)$ on $\Omega$ and $u(x_{0},t_{0})=M\ (u(x_{0},t_{0})=m)$ at some point $(x_{0},t_{0})\in\Omega_{p}$. If $c(x,t)\not\equiv 0$, suppose that $M\geq 0\ (m\leq 0)$. Then

[TABLE]

Proof.

Suppose that $u(x_{0},t_{0})=M$ at some point $(x_{0},t_{0})\in\Omega_{p}$. We distinguish two cases: (i) $x_{0}\not\in E_{r}$; (ii) $x_{0}\in E_{r}$.

In Case (i), clearly $x_{0}$ is an interior point of some edge $k_{j}$. The direct application of the strong maximum principle for Euclidean domains (see, for example, [47, Theorem 4, Chapter 3]) gives that

[TABLE]

Let $k_{h}$ be an arbitrary edge such that $k_{h}\cap k_{j}=\{e_{i}\}$, where $e_{i}$ is an end point of $k_{j}$. If there exists an interior point $y_{0}$ of $k_{h}$ such that $u_{h}(y_{0},t_{0})=M$, then

[TABLE]

due to the strong maximum principle for domains. If such interior maximum point does not exist, then it is necessary that

[TABLE]

Thus, we can claim that for some $0<\hat{t}<t_{0}$, it holds

[TABLE]

Suppose that such a claim is false. Then we can find a sequence $\{(\hat{x}_{n},\hat{t}_{n})\}_{n=1}^{\infty}$ with $\hat{x}_{n}\in(0,l_{h})$, $\hat{t}_{n}<\hat{t}_{n+1}$ for all $n\geq 1$ and $\hat{t}_{n}\to t_{0}$ as $n\to\infty$ such that $u_{h}(\hat{x}_{n},\hat{t}_{n})=M$ for all $n\geq 1$. By the strong maximum principle for domains again, one has $u_{h}(x_{h},\hat{t}_{n})=M$ for all $x_{h}\in[0,l_{h}]$. Since $\hat{t}_{n}\to t_{0}$ as $n\to\infty$ and $u_{h}$ is continuous on $[0,l_{h}]\times[0,t_{0}]$, it easily follows that (A.1) holds, which contradicts with (A.2). Hence, the claim (A.3) is proved. In view of (A.3) and the fact that $u_{h}(x_{h},t)$ attains its maximum $M$ at the boundary point $(e_{i},t_{0})$ of the region $(0,l_{h})\times(\hat{t},t_{0})$, one then applies the classical Hopf boundary lemma for Euclidean domains (see [47, Theorem 3, Chapter 3]) to conclude that $d_{hi}\frac{\partial u_{h}}{\partial x_{h}}(e_{i},t_{0})>0$. Recall that $M$ is the maximum value of $u$ on $\Omega$. So for any $j=i_{1},\cdots,i_{m}$ with $j\not=h$ such that $k_{j}\cap k_{h}\neq\varnothing$, we have $d_{ij}\frac{\partial u_{j}}{\partial x_{j}}(e_{i},t_{0})\geq 0$. Therefore, it holds

[TABLE]

which is impossible due to our assumption. This contradiction shows that (A.1) must hold. As $G$ is connected, we can assert that $u=M$ in $G\times[0,t_{0}]$.

We now consider Case (ii). Take $k_{j}$ to be an arbitrary edge such that $x_{0}$ is its endpoint and $x_{0}\in E_{r}$. By what was proved in Case (i), we can suppose that $u_{j}(x_{j},t_{0})<M,\ \forall x_{j}\in(0,l_{j})$. However, the same argument as in the last two paragraphs in the proof of Case (i) leads to a contradiction. Thus, $u_{j}(x_{j},t)=M$ for all $(x_{j},t)\in[0,l_{j}]\times[0,t_{0}]$ and in turn $u=M$ in $G\times[0,t_{0}]$ by the arbitrariness of $k_{j}$. ∎

We remark that [59, Theorem 1] covers the special case of Lemma A.1 that $c(x,t)\equiv 0$ and

[TABLE]

Hence our result is more general and it is useful when dealing with upper-lower solutions. As a direct application of Lemma A.1 as the classical Hopf lemma for Euclidean domains (see, for example, [47, Theorem 3, Chapter 3]), we have the following Hopf boundary lemma for networks.

Lemma A.2**.**

Assume that $c(x,t)\geq 0$ and is bounded from above on $\Omega$. Let $u\in C(\Omega)\cap C^{2,1}(\Omega_{p})$ satisfy

[TABLE]

Suppose that $u$ is continuously differentiable at some point $(e_{i},t)\in E_{b}\times(0,T]$, $u(e_{i},t)=M$, and $u(x,t)<M$ ($>m$) for all $(x,t)\in\Omega_{p}$. If $c\not\equiv 0$, assume that $M\geq 0\ (m\leq 0)$. Then $d_{ij}u_{x_{j}}(e_{i},t)>0\ (<0)$.

As a corollary of Lemmas A.1 and A.2, we immediately obtain the following comparison principle.

Lemma A.3**.**

Assume that $c(x,t)$ is bounded on $\Omega$. Let $u\in C(\Omega)\cap C^{2,1}(\Omega_{p})$ satisfy

[TABLE]

and assume that $\frac{\partial{u}_{j}}{\partial x_{j}}(e_{i},t)$ exists for $t\in(0,T],\,e_{i}\in E_{b}$ if $\beta_{j,s}\not=0$ for some $j\in I_{N-1},\,s\in\{1,2\}$. Then $u(x,t)\geq 0$ for all $(x,t)\in\Omega$. If $u(x,0)\not\equiv 0$, then $u(x,t)>0$ for all $(x,t)\in G\backslash E_{0}\times(0,T]$.

Proof.

Denote $v(x,t)=e^{-\ell t}u(x,t)$ and take the constant $\ell>0$ to be large so that $\ell+c>0$ on $\Omega$. Elementary computation gives $v\in C(\Omega)\cap C^{2,1}(\Omega_{p})$ satisfy

[TABLE]

Hence, it follows from Lemmas A.1 and A.2 that $\min\limits_{\Omega}v(x,t)=m\geq 0$ and so $u(x,t)\geq 0$ on $\Omega$. When $u(x,0)\not\equiv 0$ (equivalently, $v(x,0)\not\equiv 0$), suppose that $u_{j}(x_{*},t_{*})=0$ for some $(x_{*},t_{*})\in\Omega_{p}$. Then $v_{j}(x_{*},t_{*})=0=\min_{\Omega}v(x,t)$, which implies $v\equiv 0$ on $\Omega$ by Lemma A.1, a contradiction. Thus, $u(x,t)>0$ for all $(x,t)\in\Omega_{p}$. Additionally, Lemma A.2 implies that $u(x,t)>0$ for all $(x,t)\in E_{b}\backslash E_{0}\times(0,T]$. ∎

We remark that [60, Theorem 1] states another type of comparison principle for the parabolic problem on graphs.

A.2 Linear parabolic problem

In this subsection, we aim to establish the existence, uniqueness, $L^{p}$ and Schauder estimates of solutions to the following linear parabolic problem:

[TABLE]

associated with the initial condition, boundary and interior connection conditions:

[TABLE]

where $T>0$ is a fixed number.

We would like to mention that when the initial data $u^{0}$ is smooth and satisfies compatibility conditions, von Below [57] already studied the solvability of (A.6). Here we want to discuss the same issue for the initial data $u^{0}\in L^{p}(G)\,(p>1)$ by appealing to the semigroup theory used in [41, 4, 13]. To the end we need to make a transformation so that problem (A.6) can be written in the form that the framework of [41] can apply.

Let

[TABLE]

where $\eta_{j}$ is a constant to be determined on edge $k_{j}$. Then (A.6) can be written as

[TABLE]

where

[TABLE]

Choose one upstream vertex and reorder the vertices and edges such that the chosen vertex is $e_{1}$, the edge connecting to $e_{1}$ is $k_{1}$, and the endpoint $e_{2}$ of $k_{1}$ connects to edges $k_{2}$, $k_{3}$, $\cdots$, and $k_{m}$. Then at this interior vertex $e_{2}$, the boundary condition is $\sum\limits_{j=1}^{m}d_{2j}A_{j}\left(D_{j}\frac{\partial u_{j}}{\partial x_{j}}\right)(e_{2})=0$. Define $\eta_{1}=1$ on the edge $k_{1}$. Choose suitable $\eta_{2}$, $\cdots$, $\eta_{m}$ such that

[TABLE]

where $p_{j}(e_{2})=p_{j}(0)$ or $=p_{j}(l_{j})$ depending on whether $e_{2}$ is the starting point or ending point of $k_{j}$. Then the interface boundary condition at $e_{2}$ is equivalent to $\sum\limits_{j=1}^{m}d_{2j}p_{j}(e_{2})\frac{\partial u_{j}}{\partial x_{j}}(e_{2})=0$. Since $G$ is a tree, we can similarly choose the values for other $\eta_{j}$’s and rewrite the interface boundary conditions (2.10b) at interior vertices as

[TABLE]

Introduce the inner product for functions $\psi,\,\phi\in L^{2}(G)$ as

[TABLE]

Then the differential operator $\mathcal{A}$ with the domain given as in [41, Lemma 4.11] is self-adjoint with respect to $\langle\cdot,\cdot\rangle$, and the similar analysis to that of [41] shows that $\mathcal{A}$ generates a compact, contractive and positive strongly continuous semigroup.

Denote

[TABLE]

We now state the following solvability result, and $L^{p}$ and Schauder estimates for problem (A.6)-(A.7).

Theorem A.4**.**

Assume that $c\in L^{\infty}(\Omega)$, $g\in L^{p}(\Omega)$ for fixed $p>1$. Then the initial boundary value problem (A.6)-(A.7) is well-posed on $L^{p}(G)$, i.e., for any initial data $u^{0}\in L^{p}(G)$, (A.6)-(A.7) admits a unique strong solution $u\in W^{2,1}_{p}(\Omega)$ for $t>0$ that continuously depends on the initial data. Moreover, the following estimates hold.

1. (i)

If $u^{0}\in W^{2}_{p}(G)$, then the unique solution $u$ satisfies

[TABLE]

for some constant $C>0$ independent of $u,\,u^{0},\,g$. 2. (ii)

If $c\in C^{\alpha,\alpha/2}(\Omega)$, $g\in C^{\alpha}(\Omega),\,u^{0}\in C^{2+\alpha}(G)$ for some $\alpha\in(0,1)$ and $u^{0}$ solves (A.6)-(A.7) at $E_{b}\times\{0\}$ and (2.10b) at $E_{r}\times\{0\}$, and $[D_{j}\frac{\partial^{2}u^{0}_{j}}{\partial x_{j}^{2}}-v_{j}\frac{\partial u^{0}_{j}}{\partial x_{j}}+c_{j}(\cdot,0)u^{0}_{j}+g_{j}(\cdot,0)](e_{i})=[D_{h}\frac{\partial^{2}u^{0}_{h}}{\partial x_{h}^{2}}-v_{h}\frac{\partial u^{0}_{h}}{\partial x_{h}}+c_{j}(\cdot,0)u^{0}_{h}+g_{h}(\cdot,0)](e_{i})$ when $k_{j}\cap k_{h}=\{e_{i}\}$, then the unique solution $u\in C^{2+\alpha,1+\frac{\alpha}{2}}(\Omega)$ and satisfies

[TABLE]

for some constant $C>0$ independent of $u,\,u^{0},\,g$.

Proof.

Note that (A.6) with interface condition (2.10b) can be written as (A.9) with condition (A.11), and the differential operator $\mathcal{A}$ is self-adjoint and generates a compact, contractive and positive strongly continuous semigroup. Thus, by adjusting the scalar product to (A.12) and defining corresponding functions and operators based on the boundary conditions in (A.7), the analysis in [41] (see also [4, 13]) can be borrowed to show that problem (A.6)-(A.7) with the initial data $u^{0}\in L^{p}(G)$ has a unique classical solution for $t>0$ that continuously depends on the initial data.

The Schauder estimates in the assertion (ii) have been derived by [57]. The $L^{p}$-estimates in the assertion (i) follow similarly as in the proof of the Theorem in [57]. The proof consists mainly of showing that the initial boundary value problem (A.6)-(A.7) is equivalent to a well-stated initial boundary value problem for a parabolic system, where the $L^{p}$-estimate results of [26, 52] for such a parabolic system can be applied. The details are omitted here. ∎

A.3 Nonlinear problem (IBVP)

This subsection is devoted to the existence, uniqueness and positivity of solutions to the nonlinear problem (IBVP). Assumptions [H1] and [H2] are assumed throughout this section.

We begin by introducing the definition of upper and lower solution associated with problem (IBVP) as follows.

Definition A.5**.**

1. A function $\overline{u}\in C^{2,1}(\Omega)$ is an upper solution of problem (IBVP) if $\overline{u}$ satisfies the following conditions

[TABLE] 2. A function $\underline{u}(x,t)\in C^{2,1}(\Omega)$ is a lower solution of problem (IBVP) if $\underline{u}$ satisfies the following conditions:

[TABLE]

According to the definition of upper and lower solutions, one can easily see from Lemmas A.2 and A.3 that

Lemma A.6**.**

Assume that $\overline{u}$ and $\underline{u}$ is a pair of upper and lower solutions of problem (IBVP) and $\overline{u}\geq\underline{u}$ on $G\times\{0\}$. Then $\overline{u}\geq\underline{u}$ on $\Omega$. If additionally $\overline{u}\geq,\not\equiv\underline{u}$ on $G\times\{0\}$, then $\overline{u}>\underline{u}$ on $(G\backslash E_{0})\times(0,T]$.

With the aid of the above preliminaries, we are now able to state the main result of this subsection.

Theorem A.7**.**

For any $u^{0}\in L^{p}(G)\,(p>1)$ with $u^{0}\geq 0$ a.e. in $G$, problem (IBVP) admits a unique classical solution $u$ for all $t>0$, which satisfies $u\geq 0$ in $\Omega$. If additionally, $u^{0}\not\equiv 0$, then $u(x,t)>0$ for all $t>0$ and $x\in G\setminus E_{0}$.

Proof.

Note that [math] and $M^{\ast}=\max\limits_{j\in I_{N-1}}\{M_{j}\}$ forms a pair of upper and lower solutions to (IBVP). Thus, in light of Theorem A.4 and Lemma A.6, the existence of strong solution follows from the standard iteration of lower and upper solutions; the obtained solution is classical due to Theorem A.4 again. We omit the details of the proof here and refer interesting readers to [45, 66]. The uniqueness and positivity of solutions are obvious consequences of Lemma A.6. The proof is thus complete. ∎

From now on, given $u^{0}\in L^{p}(G)$ for some $p>1$, denote by $u(x,t,u^{0})$ the unique solution to problem (IBVP). Clearly, we have

Lemma A.8**.**

For any $\psi_{1},\psi_{2}\in L^{p}(G)$ with $\psi_{1}\geq,\not\equiv\psi_{2}$ on $G$, $u(x,t,\psi_{1})>u(x,t,\psi_{2})$ for all $x\in G\setminus E_{0}$ and $t>0$.

We also have the following observation.

Lemma A.9**.**

If [H3] is also satisfied, for any $u^{0}\in L^{p}$ with $u^{0}\geq 0$ on $G$ and $\lambda\in(0,1)$, $u(x,t,\lambda u^{0})\geq\lambda u(x,t,u^{0})$ on $G$ for all $t>0$.

Proof.

For any $\lambda\in(0,1)$, clearly $\lambda u(x,t,u^{0})$ satisfies

[TABLE]

where we used assumption [H2]. It follows from Lemma A.6 that $\lambda u(x,t,u^{0})\leq u(x,t,\lambda u^{0})$ on $G$ for all $t>0$. ∎

A.4 Theory of elliptic equations

It is clear that Lemmas A.1 and A.2 imply the following strong maximum principle for elliptic equations and Hopf type boundary lemma.

Lemma A.10**.**

Assume that $c(x)\geq 0$ and is bounded from above on $G$. Let $u\in C(G)\cap C^{2}(G\backslash E_{b})$ satisfy

[TABLE]

and

[TABLE]

Suppose that $u\leq M\ (u\geq m)$ on $G$ and $u(x_{0})=M\ (u(x_{0})=M)$ at some point $x_{0}\in G\backslash E_{b}$. If $c(x)\not\equiv 0$, suppose that $M\geq 0\ (m\leq 0)$. Then

[TABLE]

Lemma A.11**.**

Assume that $c(x)\geq 0$ and is bounded from above on $G$. Let $u\in C(G)\cap C^{2}(G\backslash E_{b})$ satisfy

[TABLE]

Suppose that $u$ is continuously differentiable at some point $e_{i}\in E_{b}$, $u(e_{i})=M$, and $u(x)<M$ ($>M$) for all $x\in G$. If $c\not\equiv 0$, assume that $M\geq 0\ (M\leq 0)$. Then $d_{ij}u_{x_{j}}(e_{i})>0\ (<0)$.

The following comparison principle immediately follows from Lemmas A.10 and A.11.

Lemma A.12**.**

Assume that $c(x)\geq 0$, is bounded from above on $\Omega$ and $c_{j}\beta_{j,i}\not\equiv 0$ for some $j\in I_{N-1},\,i\in\{1,2\}$. Let $u\in C(G)\cap C^{2}(G\backslash E_{b})$ satisfy

[TABLE]

and assume that $\frac{\partial{u}_{j}}{\partial x_{j}}(e_{i})$ exists for $e_{i}\in E_{b}$ if $\beta_{j,s}\not=0$ for some $j\in I_{N-1},\,s\in\{1,2\}$. Then $u(x)\geq 0$ for all $x\in G$. If $u(x)\not\equiv 0$, then $u(x)>0$ for all $x\in G\backslash E_{0}$.

In what follows, we will establish the existence, uniqueness, $L^{p}$ and Schauder estimates of solutions to the following linear elliptic problem:

[TABLE]

Indeed, by using the similar idea to that of [57], we can write the boundary value problem on $G$ in (A.16) into an equivalent boundary value problem for an elliptic system and then obtain the following result about the existence and a priori estimates of solutions of (A.16).

Theorem A.13**.**

The following assertions hold.

1. (i)

Assume that $c$ is bounded on $G$ with $c(x)\geq,\not\equiv 0$ and $g\in L^{p}(G)$ ($p>1$), then (A.16) admits a unique strong solution $u\in W_{p}^{2}(G)$ and

[TABLE]

where the constant $C$ does not depend on $u,\,g$. 2. (ii)

Assume that $c\in C^{\alpha}(G)$ with $c(x)\geq,\not\equiv 0$ and $g\in C^{\alpha}(G)$, then (A.16) admits a unique solution $u\in C^{2+\alpha}(G)$ and

[TABLE]

where the constant $C$ does not depend on $u,\,g$.

Next, we develop the theory of upper and lower solutions to establish the existence and uniqueness of solution to the following nonlinear elliptic equations

[TABLE]

Definition A.14**.**

1. A function $\overline{u}\in C^{2}(G)$ is an upper solution of (A.17) if $\overline{u}$ satisfies

[TABLE] 2. A function $\underline{u}(x,t)\in C^{2}(G)$ is a lower solution of (A.17) if $\underline{u}$ satisfies

[TABLE]

Based on Lemma A.12 and Theorem A.13, one can use the standard iteration of lower and upper solutions (see, for instance [45, 66]) to conclude that

Theorem A.15**.**

Let $\overline{u}$ and $\underline{u}$ be a pair of upper and lower solutions of (A.17) satisfying $\overline{u}\geq\underline{u}$ on $G$ and $m=\min_{G}\underline{u}<M=\max_{G}\overline{u}$, and

[TABLE]

for some constants $K>0$ and $\alpha\in(0,1)$. Then (A.17) admits a solution $u\in C^{2+\alpha}(G)$ which satisfies $\underline{u}\leq u\leq\overline{u}$. Moreover, (A.17) admits a minimal $\tilde{w}$ and a maximal solution $\tilde{u}$ in $[\underline{u},\overline{u}]$ in the sense that for any solution $w$ of (A.17) satisfying $\underline{u}\leq w\leq\overline{u}$, we have $\tilde{w}\leq w\leq\tilde{u}$.

Appendix B Proof of Proposition 3.1

Choose $\xi>0$ large enough so that $f_{j}(\cdot,0)-\xi<0$ for all $j\in I_{N-1}$. For any $g\in X$, Theorem A.13 guarantees that the problem

[TABLE]

has a unique solution $u$ satisfying

[TABLE]

for some constants $C>0$ independent of $u$ and $g$.

Define the operator

[TABLE]

Then $T$ is a linear and continuous operator that maps a bounded set in $X$ into a bounded set in $C^{2+\alpha}(G)$. Note that a bounded set in $C^{2+\alpha}(G)$ is a sequentially compact set in $X$. This implies that $T$ maps a bounded set in $X$ into a sequentially compact set in $X$. Hence, $T$ is a compact operator on $X$. Moreover, by Lemmas A.2 and A.3, $Tg\geq 0$ if $g\in X_{+}$, and $u=Tg\in X^{o}$. Therefore, $T$ is strongly positive. Let $r(T)$ be the spectral radius of $T$. It follows from the well-known Krein-Rutman Theorem (see, for example, [10, Theorem 1.2]) that $r(T)>0$ is a simple eigenvalue of $T$ with an eigenfunction $g^{\ast}\in X^{0}$, i.e., $Tg^{\ast}=r(T)g^{\ast}$, and there is no other eigenvalue of $T$ associated with positive eigenfunctions. Thus, $\psi^{\ast}=Tg^{\ast}$ satisfies $-\mathcal{L}\psi^{\ast}+\xi\psi^{\ast}=(1/r(T))\psi^{\ast}$ in $G$, and hence, $\lambda^{\ast}=\xi-1/r(T)$ is a simple eigenvalue of (3.5) with positive eigenfunction $\psi^{\ast}\in X^{o}$ and no other eigenvalues of (3.5) correspond to positive eigenfunctions. Similarly as in the proof of [10, Theorem 1.4], we can obtain that if $\lambda\neq\lambda^{\ast}$ is an eigenvalue of (3.5), then $Re(\lambda)\leq\lambda^{\ast}.$

Appendix C Proof of Theorem 3.3

We first prove (i). Assume that $\lambda^{\ast}<0$. Let $\phi\in X_{+}$. Since $\psi^{\ast}\in X^{o}$, clearly there exists $\sigma>0$ such that $0\leq\phi\leq\sigma\psi^{\ast}$ on $G$. Let $u(\cdot,t,\phi)$ be the solution of (IBVP) with initial condition $\phi$ and $\overline{u}(\cdot,t,\sigma\psi^{\ast})=e^{\lambda^{\ast}t}\sigma\psi^{\ast}$ be the solution of (3.4) with initial condition $\sigma\psi^{\ast}$. By [H1], we have

[TABLE]

for any $x_{j}\in(0,l_{j})$ and $t>0$. Then

[TABLE]

on $(0,l_{j})$ for any $t>0$. It follows from Lemma A.3 and the fact $0\leq\phi\leq\sigma\psi^{\ast}$ that

[TABLE]

for any $t\geq 0$. Therefore,

[TABLE]

which implies that $u\equiv 0$ is globally attractive for all initial conditions in $X_{+}$.

Let $\lambda^{\ast}=0$. For any $\phi\in X^{o}$, there exists some $\sigma_{0}>0$ such that $0\leq\phi\leq\overline{u}:=\sigma_{0}\psi^{\ast}$ and that $\overline{u}$ is an upper solution of (3.8). Let $u^{(1)}$ and $u^{(2)}$ be solutions of (IBVP) with initial conditions $\phi$ and $\overline{u}$, respectively. It follows from Theorem A.3 that $0\leq u^{(1)}(x,t)\leq u^{(2)}(x,t)\leq\overline{u}(x)$. Note that $u^{(2)}$ is bounded and monotonically decreasing in $t$ by following the same proof of [66, Lemma 3.2.4]). Therefore, the similar proof of [66, Lemma 3.2.5]) shows that $\lim\limits_{t\rightarrow\infty}u^{(2)}(\cdot,t)=V\geq 0$ and $V$ is a classical solution of (3.8). If $V(x)>0$ at some $x\in G$, clearly $V\in X^{o}$. Then, by Proposition 3.2, we can easily obtain $\lambda^{\ast}>0$, which gives rise to a contradiction. Hence, $V(x)\equiv 0$ on $G$. Therefore, $\lim\limits_{t\rightarrow\infty}u^{(1)}(x,t)=0$. So we conclude that when $\lambda^{\ast}=0$, $u\equiv 0$ is  globally attractive for (IBVP) with respect to all initial conditions in $X_{+}$. Thus, (i) is verified.

We next prove (ii). Assume that $\lambda^{\ast}>0$. Let $\psi^{\ast}$ be the eigenfunction associated with $\lambda^{\ast}$. For sufficiently small $\epsilon>0$, we have $f_{j}(x_{j},\epsilon)\geq f_{j}(x_{j},0)-\lambda^{\ast}$ for all $x_{j}\in(0,l_{j})$, $j\in\{1,\cdots,N-1\}$. This implies that $w_{1}=\epsilon\psi^{\ast}$ is a lower solution of (3.8). Note that for any constant $K^{\ast}>\max\{M_{1},\cdots,M_{N-1}\}$, $w_{2}(x)=K^{\ast}$ is an upper solution of (3.8). Let $\epsilon>0$ be sufficiently small such that $w_{2}(x)\geq w_{1}(x)=\epsilon\psi^{\ast}(x)$ and $K^{\ast}>\min_{x\in G}\{\epsilon\psi^{\ast}(x)\}$. It follows from Theorem A.15 that (3.8) admits a positive solution $u^{\ast}\in X^{o}$.

Assume there are two distinct positive steady states of (IBVP): $u_{1}^{\ast}$ and $u_{2}^{\ast}$. Since $K^{*}$ can be arbitrarily large and $\epsilon$ can be arbitrarily small, without loss of generality, assume that $u_{1}^{\ast}$ is the maximal solution of (3.8) and $u_{2}^{\ast}$ is the minimal solution of (3.8) in $[\epsilon\psi^{\ast},K^{\ast}]$. Then $u_{2}^{\ast}\leq u_{1}^{\ast}$ on $G$. It suffices to show $u_{2}^{\ast}=u_{1}^{\ast}$. Suppose that $u_{2}^{\ast}\leq,\not=u_{1}^{\ast}$ on $G$. Recall that $u_{1}^{\ast},\,u_{2}^{\ast}\in X^{o}$. By defining $\tau^{0}=\sup\{\tau>0:\ u_{2}^{\ast}\geq\tau u_{1}^{\ast}\ \mbox{on}\ G\}$, we then have $\tau^{0}\in(0,1)$ and $u_{2}^{\ast}\geq\tau^{0}u_{1}^{\ast}$ on $G$. Due to [H3], we further observe that $u_{2}^{\ast}\geq,\not\equiv\tau^{0}u_{1}^{\ast}$. Thus, by Lemmas A.2, A.8 and A.9, we get

[TABLE]

where $Q_{t}$ is the solution map of (IBVP) defined as $Q_{t}(\psi)=u(x,t,\psi)$ for the solution $u(x,t,\psi)$ of (IBVP) with initial condition $\psi$. This implies that $u_{2}^{\ast}(x)-\tau^{0}u_{1}^{\ast}(x)\in X^{o}$, which in turn implies that $u_{2}^{\ast}(x)-\tau^{0}u_{1}^{\ast}(x)\geq\tau_{0}u_{1}^{\ast}(x)$ on $G$ for some small $\tau_{0}>0$. This is a contradiction to the definition of $\tau^{0}$. Therefore, there is only a unique positive steady state $u^{\ast}$ of (IBVP).

For any $u^{0}\in X_{+}\setminus\{0\}$, obviously the unique solution $u$ of (IBVP) satisfies that $u(\cdot,t)\in X^{o}$ for any $t>0$. Thus, we can assume that $u^{o}\in X^{o}$. So there exist some $\epsilon_{0}>0$ and $\sigma_{0}\geq 1$ such that $\underline{u}=\epsilon_{0}\psi^{\ast}$ and $\overline{u}=\sigma_{0}K^{\ast}$ are lower and upper solutions of (3.8), respectively, and

[TABLE]

Let $u_{1},\,u_{2}$ be solutions of (IBVP) with initial conditions $\underline{u}$ and $\overline{u}$, respectively. It follows from Lemma A.6 that $\underline{u}(x)\leq u_{1}(x,t)\leq u(x,t)\leq u_{2}(x,t)\leq\overline{u}(x)$. As before, it can be easily proved that $u_{1}$ and $u_{2}$ are monotonically increasing and decreasing in t, respectively ([66, Lemma 3.2.4]). Therefore, $u_{1}$ and $u_{2}$ are bounded and monotonic with respect to $t$. Additionally, we can claim $\lim\limits_{t\rightarrow\infty}u_{1}(\cdot,t)=U$ and $\lim\limits_{t\rightarrow\infty}u_{2}(\cdot,t)=V$. Furthermore, we can prove $U$ and $V$ are solutions of (3.8) (see Lemma 3.2.5 in [66]). Then $\underline{u}(x)\leq U(x)\leq V(x)\leq\overline{u}(x)$. Hence, $U(x)=V(x)=u^{\ast}$, and $\lim\limits_{t\rightarrow\infty}u(x,t)=u^{\ast}(x)$. Therefore, we have proved that $u^{\ast}$ is globally attractive with respect to any positive initial values in $X_{+}\setminus\{0\}$.

Appendix D The hydrological relation in a gradually varying flow

Recall that the governing equation for the gradually varied flow is given by

[TABLE]

(see (5-7) in [7]), where $x$ (unit: m) represents the longitudinal location along the river, $y(x)$ (unit: m) is the water depth at location $x$, $S_{0}(x)$ is the slope of the channel bed at location $x$, $S_{f}$ is the friction slope, i.e., the slope of the energy grade line, $F_{r}$ is the Froude number that is defined as the ratio between the flow velocity and the water wave propagation velocity. In the case where the river has a rectangular cross section with a constant width $B$ (unit: m) and a constant bed slope $S_{0}$, the water depth $y(x)$ is stabilized at the normal depth

[TABLE]

where $Q$ (unit: m3/s) is the flow discharge, $k=1$ is a dimensionless conversion factor, and $n$ (unit: s/m1/3) is Manning’s roughness coefficient, which represents the resistance to water flows in channels and depends on factors such as the bed roughness and sinuosity. The flow in such a river is called a uniform flow. See more details in [7] or Appendices C and D in [22].