Mathematics in Applied Sciences and Engineering https://ojs.lib.uwo.ca/mase Volume 2, Number 4, December 2021, pp.235-272 https://doi.org/10.5206/mase/14112 PREVAILING WINDS AND SPRUCE BUDWORM OUTBREAKS: A REACTION-DIFFUSION-ADVECTION MODEL ABBY ANDERSON AND OLGA VASILYEVA Abstract. We extend the classical reaction-diffusion model for spatial population dynamics of spruce budworm on a finite domain with hostile boundary conditions by including an advection term repre- senting biased unidirectional movement of individuals due to a prevailing wind. We use phase plane techniques to establish existence of a critical value of advection speed that prevents outbreak solutions on any finite domain while possibly allowing an endemic solution. We obtain lower and upper bounds for this critical advection value in terms of biological parameters involved in the reaction term. We also perform numerical simulations to illustrate the effect of advection on the dependence of the domain size on the maximal population density of a steady state solution and on critical domain sizes for endemic and outbreak solutions. The results are also applicable to other ecological settings (rivers, climate change) where a logistically growing population is subject to predation by a generalist, diffusion and biased movement. 1. Introduction This paper is inspired by classical 1979 work [11] of D. Ludwig, D.G. Aronson and H.F. Weinberger where the authors introduced and analyzed a reaction-diffusion model for the spatial dynamics of spruce budworm (SBW) population in which the insects were assumed to disperse by diffusive movement only. On page 251 in [11], the authors state the following: “The limitation that dispersal is modelled by pure diffusion is a serious one. There are prevailing wind patterns over New Brunswick which produce systematic, convective motions of the budworm. The addition of a convective term to the present model presents no particular difficulties for the treatment of the differential equation. The boundary conditions introduce serious complications.” In this paper, we extend the reaction-diffusion model from [11] to include the convective/advective motion in the form of an advection term. The boundary conditions remain hostile (lethal), as they were in [11]. The main difficulty in generalizing the results from [11] is the presence of the advection term which requires new techniques since as the first integral method and comparison methods used in [11] are not applicable anymore. We use the geometric method based on the phase plane analysis of a boundary value problem given by a system of two first order ODEs obtained from the second order ODE describing the profile of a steady state solution of the original reaction-diffusion-advection (RDA) problem. This approach was used in [1, 22, 21, 2] to study steady states or traveling wave solutions arising in nonlinear RDA models. Received by the editors 18 July 2021; 24 October 2021; published online 3 November 2021. 2010 Mathematics Subject Classification. Primary 35K57; Secondary 34B15. Key words and phrases. reaction-diffusion-advection models, steady state, critical advection, stable and unstable manifolds, spruce budworm modeling. A. Anderson was supported in part by an Undergraduate Student Research Award from the Natural Sciences and Engineering Research Council of Canada. O. Vasilyeva was supported in part by a Discovery Grant from the Natural Sciences and Engineering Research Council of Canada (RGPIN-2017-04376). 235 https://ojs.lib.uwo.ca/mase https://dx.doi.org/https://doi.org/10.5206/mase/14112 236 A. ANDERSON AND O. VASILYEVA Historically, RDA models have been used to adrress the well-known ”drift paradox” in river ecology: the ability of population to persist in finite river segments despite being constantly washed downstream [20, 13, 15, 4]. The RDA approach had also been used to model the dynamics of sinking phytoplankton and populations experiencing a shift in habitat boundaries due to climate change [7, 17, 1]. A recent survey on a variety of RDA models is given in [9]. The development of RDA models progressed from the linear reaction term case in [20], to monostable logistic term (e.g. in [22, 10]), and to the bistable Allee effect term [3, 23]. In some RDA models, additional complexity was introduced by incorporating spatial dependence of the reaction term as in [23] or by considering river networks instead of the single river [18, 19, 21, 6, 2]. The mathematical motivation of our paper is in extending the RDA techniques to the case of a reaction term that combines the logistic growth with predation by a generalist. The resulting nonlinear reaction term has up to four zeros. In the non-spatial setting considered in [12], zero state is an unstable equilibrium, and there are stable endemic and outbreak equilibria, with the latter two separated by an unstable threshold. As the complexity of the reaction term increases, it is natural to expect multiple steady state solutions, and indeed, in [11], the main focus was on the analysis of endemic and outbreak steady state solutions. In Section 2 (Preliminaries), we give an overview the classical nonspatial and spatial models intro- duced in [12] and [11]. First, we outline how by scaling the original nonspatial model one obtains a model involving two biological parameters Q and R [12]. Then we focus on the “cusped” region in the QR-plane for which the nonspatial model has exactly three positive equilibria (we denote this region by Ω) . Then we describe the spatial (reaction-diffusion) model [11] focusing on the subregion Ω∗ of Ω for which the spatial model admits both endemic and outbreak solutions. We also recall the behavior of the spatial model when the parameters (Q,R) are taken outside of Ω∗, as described in [11]. In Section 3, we set up an extension of the spatial model from [11] obtained by adding an advection term. Using our geometric approach, we study steady state solutions of the resulting boundary value problem by reducing it to a nonlinear ODE system and analyzing its orbits in the phase plane. We identify the equilibria of the ODE system and their nature depending on the advection velocity q. We introduce notation for stable and unstable manifolds of the two saddle point equilibria of the system that will play the crucial role in our analysis and find the slopes of the stable and unstable manifolds at these saddle points. In Section 4, we focus on the behavior of unstable manifold of the first saddle point as we increase advection and prove technical lemmas that allow us to establish the existence of a critical value of advection q∗cr beyond which the outbreak solution is not possible on any finite domain. We also obtain an upper bound on q∗cr. We perform further geometric analysis of the vector field associated with the ODE system and construct a piecewise linear approximation of the orbit connecting two saddle points that allows us to obtain a lower bound on q∗cr (the details are given in Appendix A). In Section 5, we perform a detailed analysis of behavior of various orbits in the phase plane using trapping region method as one of the main tools. This leads to our main result that describes the global picture of orbits in the phase plane. Depending on the advection speed, we classify all possible scenarios in terms of orbits representing positive steady state solutions. Endemic and outbreak orbits are classified according to the maximal density or the upstream or downstream density gradients. We also discuss the dependence of habitat length on the maximal density of a steady state solution. Numerical simulations are performed in Section 6. We test the estimates for critical advection q∗cr for outbreak solutions, and explore the effect of advection on the relation between the habitat length and the maximal population density, and how the increase of advection speed changes the nature of steady states. Section 7 contains a discussion of our results and provides their biological interpretation. PREVAILING WINDS AND SPRUCE BUDWORM OUTBREAKS: AN RDA MODEL 237 2. Preliminaries: an overview of the classical non-spatial and spatial spruce budworm models In this section, we recall the classical non-spatial [12] and spatial [11] model for spruce budworm (SBW) dynamics, summarizing and illustrating the relevant results from these two papers. For the sake of completeness and for the reader’s convenience, we expand some of the arguments appearing in [12, 11] by adding details and providing plots and diagrams. 2.1. The non-spatial Ludwig-Jones-Holling model. The non-spatial SBW model introduced in [12], is given by dB dT = rB ( 1 − B K′S ) −β B2 (α′S)2 + B2 , (2.1) where B is the budworm density (larvae per acre) at time T (measured in years). The first term on the right gives the logistic population growth in the absense of predators. Here r is the per capita growth rate at low density, K′ is the carrying capacity per branch, and S is the number of (standard size) branches per acre. The second term on the right is the predation term. It uses a Type 3 functional response which takes into account the predation pattern between the spruce budworm and generalist predators (birds). Namely, when the budworm density is low, the predation rate is very small as the birds tend to ignore the budworms and feed on other sources. When the budworm density increases, the birds start noticing their presence and start actively feeding on them. Finally, if the density increases beyond a certain point, the birds become saturated and as the result the predation level stabilizes at a constant value β, the limiting consumption rate (larvae per year), or saturation rate. Here α′S is the density of budworm (larvae per acre) that results in the predators consuming at half saturation rate (β 2 ). Using u = B α′S , t = rT, R = rα ′S β and Q = K ′ α′ (2.1) is reduced to the non-dimensional version du dt = φ(u; Q,R), (2.2) where φ(u; Q,R) = u− u2 Q − 1 R u2 1 + u2 = u ( 1 − u Q ) − 1 R u2 1 + u2 . (2.3) The number of equilibria of (2.2) varies between 2 and 4. Indeed, they are the roots of the equation u ( 1 Q u3 −u2 + ( 1 Q + 1 R ) u− 1 ) = 0. We will denote the zero solution of the above equation by u0. As in [11], in this paper we will focus on the case when the above equation has three distinct positive roots. It was observed in [11] that the set of pairs (Q,R) for which (2.2) has three positive equilibria corresponds to a region in the QR-plane bounded by two curves R = R1(Q) and R = R2(Q) where Q > 3 √ 3 and R2(Q) < R1(Q). Here, the lines passing through (Q, 0) and (0,Ri(Q)), i = 1, 2, are tangent to the curve v = u 1 + u2 in the uv-plane. Fixing Q > 3 √ 3 and varying R from R2(Q) to R1(Q), one can observe that the line v = R ( 1 − u Q ) 238 A. ANDERSON AND O. VASILYEVA Figure 2.1. Curve in the QR-plane bounding the region Ω corresponding to (2.2) having three positive equilibria. transitions from being tangent to the curve v = u 1 + u2 at a point ( ū2, ū2 1+ū22 ) to having three intersection points with the curve, and to being tangent to the curve at a point ( ū1, ū1 1+ū21 ) , where 1 < ū1 < ū2 < Q. Let Ω ⊂ R2 be the region in the QR-plane defined by Ω = {(Q,R) : φ(u; Q,R) has three positive zeros}. Note that Ω is an open set.To obtain the boundary of Ω, one can simply trace a parametric curve (Q(s),R(s)) where s > 1 and the tangent line to the graph of v = u u2+1 at u = s has the u-intercept Q(s) and the v-intercept R(s). In fact, it is given by the following parametric equation: (Q(s),R(s)) = ( s + s(s2 + 1) s2 − 1 , 2s3 (s2 + 1)2 ) , s > 1. A plot of (Q(s),R(s)) is given in Figure 2.1. From now on, we will be working with the case when (2.2) has three positive equilibria. For simplicity, we will denote φ(u) = φ(u; Q,R). Thus, we assume that φ(u) has four distinct roots u0 = 0 < u1 < u2 < u3, the equilibria of system (2.2). In this case, φ ′(u0),φ ′(u2) > 0 and φ ′(u1),φ ′(u3) < 0 (see the graph in Figure 2.2 for Q = 6.5 and R = 0.6). It follows that the equilibria u0,u2 of (2.2) are unstable and u1,u3 are stable. As in [11], we will refer to the equilibria u0 = 0,u1,u2,u3 as the extinction, endemic, threshold and outbreak equilibria, respectively. Indeed, u ≡ 0 corresponds to the complete absence of individuals, the lower of the two stable states corresponds to the endemic situation, while the higher stable state describes the population outbreak. The unstable steady state separating the endemic and the outbreak equilibria serves as a threshold. PREVAILING WINDS AND SPRUCE BUDWORM OUTBREAKS: AN RDA MODEL 239 Figure 2.2. Graph of φ(u) = u− 1 Q u2 − 1 R u2 1+u2 for Q = 6.5,R = 0.6. The following lemma describing the behavior of the function φ(·; Q,R) will be useful in our further analysis. Lemma 2.1. For (Q,R) ∈ Ω, the function φ(·) = φ(·; Q,R) has exactly four inflection points ±û1,±û2 such that 0 < û1 < 1 < û2. Proof. Differentiating (2.3) three times we get: φ′(u) = 1 − 2 Q u− 1 R 2u (u2 + 1)2 , φ′′(u) = − 2 Q − 2 R · 1 − 3u2 (u2 + 1)3 , and φ(3)(u) = 1 R · 24u(1 −u2) (u2 + 1)4 . Note that since φ(·) has four zeros 0 = u0 < u1 < u2 < u3 in each of which it changes sign, it has at least two positive inflection points û1 < û2. Since φ ′′(·) is an even function, −û1 and −û2 are also inflection points of φ(·). Note also that φ(3)(u) has three zeros, u = −1, 0, 1. Clearly, φ′′(·) has local maxima at −1 and 1 and local minimum at 0, and therefore it has at most four zeros. Thus, φ(·) has exactly four inflection points given by ±û1 and ±û2 where 0 < û1 < 1 < û2. � 2.2. The spatial Ludwig-Aronson-Weinberger model. The spatial version of the SBW model described below was introduced in [11]. In this model, the habitat is represented by an infinite strip in the XY -plane of the form [ −L 2 , L 2 ] × R. The population density is assumed to depend only on X, and thus the authors reduce the setting to that of a finite interval [ −L 2 , L 2 ] of the real line (where the spatial variable X is measured in kilometers). Thus, the quantity of interest becomes the linear density of SBW (individuals per kilometer) given by the function u(X,T) at location X ∈ [ −L 2 , L 2 ] at time T ≥ 0 (years). In addition to the birth and death dynamics described by the reaction term discussed 240 A. ANDERSON AND O. VASILYEVA earlier in the non-spatial setting ([12]), the individuals are assumed to move randomly. The dispersal of individuals from a fixed location over a short time interval δT is described by a normal distribution with mean 0 and variance σ2δT. Thus, the corresponding reaction-diffusion equation takes the form: ∂B ∂T = σ2 2 ∂2B ∂X2 + rB ( 1 − B K′S ) −β B2 (α′S)2 + B2 , (2.4) The boundary conditions considered in [11] are of the hostile type: u ( − L 2 ,T ) = u ( L 2 ,T ) = 0, for T > 0. Introducing the scaled variables t = rT, x = 2r σ X, u = B α′S , one obtains the non-dimensionalized reaction-diffusion equation as follows: ∂u ∂t = ∂2u ∂x2 + φ(u), (2.5) where the reaction term is given by (2.3). The boundary conditions for the non-dimensional case will be u ( − l 2 , t ) = u ( l 2 , t ) = 0, (2.6) for t > 0, where l = 2r σ L. A steady-state solution u(x) of (2.5, 2.6) satisfies the following boundary value problem on [ − l 2 , l 2 ] : u′′ + φ(u) = 0, u ( − l 2 ) = u ( l 2 ) = 0. (2.7) For our further analysis we will use the following fact about the spatial dynamics of SBW model given by (2.5, 2.6) as observed in [11]. Fact 2.2. ([11]) Let U : [ − l 2 , l 2 ] → R be a positive steady state solution of (2.5, 2.6) (solution of (2.7)). Also let µ = max− l 2 ≤x≤ l 2 U(x). Let F(u) = ∫ u 0 φ(w)dw = 1 2 u2 − 1 3Q u3 − u R + 1 R arctan u. Note that since U′′(x) + φ(U(x)) ≡ 0, the expression 1 2 (U′(x))2 + F(U(x)) in constant for − l 2 < x < l 2 . Thus, 1 2 (U′(x))2 + F(U(x)) = F(µ) for any − l 2 < x < l 2 . In particular, for any 0 < w < µ we have F(w) < F(µ). It follows from Fact 2.2 that for any positive steady-state solution U(−) of (2.5, 2.6), or equivalently, a positive solution U(−) of (2.7), the maximal value of U(−), µ = max− l 2 ≤x≤ l 2 u(x) satisfies 0 < µ < u1 or u∗ < µ < u3, where u2 < u ∗ < u3 is such that F(u ∗) = F(u1) (see Figure 2.3). Note that such u ∗ exists only when F(u3) > F(u1). As we will see later in this paper, the value u ∗ will play a crucial role in our analysis. When such u∗ does not exist, we have 0 < µ < u1. Moreover, it was shown in [11], as a consequence of Fact 2.2, that the habitat length l can be expressed explicitly in terms of the maximal density µ of the positive steady state solution, as follows: l(µ) = √ 2 ∫ µ 0 dw√ F(µ) −F(w) . Furthermore, it was shown in [11] that here exists a function R̃(Q) such that R2(Q) < R̃(Q) < R1(Q) for any Q > 3 √ 3, so that for any (Q,R) ∈ Ω, we have F(u3) > F(u1) exactly when R̃(Q) < R < R1(Q). PREVAILING WINDS AND SPRUCE BUDWORM OUTBREAKS: AN RDA MODEL 241 Figure 2.3. Graphs of F(u) = 1 2 u2 − 1 3Q u3 − u R + 1 R arctan u (solid) and φ(u) = u − 1 Q u2 − 1 R u2 1+u2 (dashed) for Q = 8.5,R = 0.5. Here u2 < u ∗ < u3 is such that F(u∗) = F(u1). Figure 2.4. Six subregions of QR-plane. In order to summarize the behavior of l = l(µ) and the nature of steady states as described in [11], we will consider subregions of the first quadrant of the QR-plane with boundaries given by Q = 3 √ 3, R = R1(Q), R = R2(Q) and R = R̃(Q), as well as R = R̄(Q) (to be introduced below). We will enumerate the regions as I-VI. See Figure 2.4. (I) {(Q,R) : 0 < Q ≤ 3 √ 3,R > 0} (this region was not studied in [11], we include it for complete- ness): here φ(u) has a unique positive zero u1; l(µ) is defined for 0 < µ < u1 and is increasing from its left limit of π to ∞ at u1. Thus, there exists a critical domain size lc1 = π such that for any l > lc1 there is a unique positive steady state solution U1(x) of (2.5, 2.6), it is bounded above by u1. For l ≤ lc1 there is no positive steady state solution. (II) {(Q,R) : Q > 3 √ 3, 0 < R ≤ R2(Q)}: When R < R2(Q) φ(u) has only one positive zero u1, for R = R2(Q) a new positive (double) zero u2 = u3 appears. The behavior of l(µ) and the nature of steady states are identical to that in region (I). (III) {(Q,R) : Q > 3 √ 3,R2(Q) < R ≤ R̃(Q)}: Now φ(u) has three positive zeros u1 < u2 < u3 and F(u1) ≥ F(u3) (with strict inequality for R < R̃(Q)); l(µ) behaves as in (I) and (II). We still 242 A. ANDERSON AND O. VASILYEVA Figure 2.5. Dependence of habitat length l on maximal denisty µ for different values of R (for a fixed Q). have only one positive steady state U1(x) bounded above by u1 for l > l c 1 = π. For l ≤ lc1 there is no positive steady state solution. (IV) {(Q,R) : Q > 3 √ 3, R̃(Q) < R < R1(Q)}: In this case φ(u) still has three positive zeros u1 < u2 < u3 and F(u1) < F(u3) (so, u ∗ exists); l(µ) has two branches, the left branch behaves as in (I-III), the right branch is defined for u∗ < µ < u3, and it has a minimal value of l c 2 and approaches infinity at u∗ and u3. Thus, there exists another critical habitat length l c 2 > l c 1 = π such that for l > lc2, there exist three distinct positive steady-state solutions U1(x),U2(x),U3(x) of (2.5, 2.6) , where U1(x) < U2(x) < U3(x) for all x ∈ ( − l 2 , l 2 ) . Denote the maximal value of Ui(x) (maximal density) on [ − l 2 , l 2 ] by µi, i = 1, 2, 3. The smallest of the solutions, U1(x), satisfies µ1 < u1 and is referred to as an endemic solution. The largest of them, U3(x), satisfies u ∗ < µ3 < u3 and is called an outbreak solution. Both U1(x) and U3(x) are asymptotically stable. The intermediate solution U2(x) satisfies u∗ < µ2 < µ3 and is unstable. If the initial density of organisms is below U2(x) (throughout the habitat) then in the long term the density approaches the endemic state U1(x) (when l > l c 1). If the initial density is above U2(x) then the density approaches the outbreak state U3(x) (provided l > lc2). For l c 1 < l ≤ lc2, there is only one (endemic) positive steady state solution bounded above by the constant u1. For l ≤ lc1 there is no positive steady state solution. (V) {(Q,R) : Q > 3 √ 3,R1(Q) ≤ R < R̄(Q)}: When R reaches R1(Q), the first two positive zeros of φ, u1 and u2, coincide, then they disappear for R > R1(Q). Let u1 denote the only positive zero that φ(u) has in this case. Any positive steady state will have its maximal density bounded by u1. As we keep increasing R until a certain value that we call R̄(Q), we observe that l(µ) is defined for 0 < µ < u1, has limit π at 0, followed by local maximum (with value l c 3), local minimum (with value lc2) and goes to infinity at u1. Thus, we still observe three distinct positive steady state solutions (stable endemic and outbreak, and unstable threshold) for l in the finite PREVAILING WINDS AND SPRUCE BUDWORM OUTBREAKS: AN RDA MODEL 243 interval lc2 < l < l c 3, while for π = l c 1 < l ≤ lc2 there is only one (endemic) positive solution. For l ≥ lc3, only a positive outbreak solution is left. (VI) {(Q,R) : Q > 3 √ 3,R ≥ R̄(Q)}: As we increase R further, φ(u) still has only one positive zero u1; now l(µ) is increasing on the interval 0 < µ < u1, has limit π at 0, and goes to infinity at u1. Thus, we can see a single positive steady state solution for all l > l c 1 = π. Therefore, R = R̄(Q) separates the three steady state case from the single steady state. The scenarios above are illustrated by l vs. µ plots in Figure 2.5. Note that for R = R1(Q), the two vertical asymptotes µ = u1 and µ = u ∗ collapse into one. One of our main goals in Sections 5 and 6 will be to explore the effect of advection on the above scenarios. Most of our analysis will be carried out in the case when φ(u) has three positive zeros and the system (2.5, 2.6) has both endemic and outbreak solutions for sufficiently large domains, that is, we will consider (Q,R) in Region (IV), which we will denote Ω∗: (Q,R) ∈ Ω∗ = {(Q,R) : Q > 3 √ 3, R̃(Q) < R < R1(Q)}. In this case, there is a clear distinction between endemic and outbreak solutions since their maximal densities are separated by the zeros of φ(u). In the region (V), the difference is more subtle and the phase plane technique that we use in this paper is not as effective. However, we do perform numerical simulations in all six regions of the QR-plane. From now on, unless specified otherwise, we will assume that (Q,R) ∈ Ω∗. 3. Model set up In our setting, all individuals from the population are subject to a unidirectional movement bias given in the form of wind transporting the larval population in one direction. This phenomenon is represented by an advection term −q∂u ∂x (where q ≥ 0). Therefore, we add the advection term to the non-dimensionalized version of a reaction-diffusion equation (2.5) to get the following PDE boundary value problem:{ ∂u ∂t = ∂ 2u ∂x2 −q∂u ∂x + φ(u), u ( − l 2 , t ) = u ( l 2 , t ) = 0. (3.1) where, as before, φ(u) = u− 1 Q u2 − 1 R u2 1+u2 . The steady state solutions of (3.1) can be viewed as solutions of the following second order ODE boundary value problem: { u′′ −qu′ + φ(u) = 0 u ( − l 2 ) = u ( l 2 ) = 0. (3.2) By introducing v(x) = u′(x), we rewrite (3.2) as the following system of first order ODE:{ du dx = v dv dx = qv −φ(u) (3.3) subject to the boundary conditions u ( − l 2 ) = u ( l 2 ) = 0. (3.4) For any l > 0, a solution of (3.3, 3.4) is represented by an orbit in the uv-plane that starts at a point on the positive v-semiaxis (when x = − l 2 ) and ends at a point on the negative v-semiaxis (when x = l 2 ). Observe that for any solution (u(x),v(x)), we have du dx > 0 (< 0) for v > 0 (v < 0), while for 244 A. ANDERSON AND O. VASILYEVA Figure 3.1. Phase plane portrait for system (3.3) when α = 0.154,β = 1.67 and q = 0. The four equilibria of (3.3) are easily identified as two centers and two saddle points. The portions of the closed orbits in the first and fourth quadrants form endemic and outbreak solutions of (3.3) corresponding to different habitat lengths. v = 0 we have dv dx = −φ(u) < 0 exactly when 0 < u < u1 or u2 < u < u3. Thus, for a fixed l > 0, the positive solutions of (3.3, 3.4) are of two kinds: those that cross the u-axis between 0 and u1 (endemic solutions/orbits) and those that cross the u-axis between u2 and u3 (outbreak solutions/orbits). Note that the u-coordinate of the intersection of an orbit (u(x),v(x)) with the u-axis is the maximum value of u(x) on the interval − l 2 ≤ x ≤ l 2 (maximum density in the habitat). Thus, we can characterize the endemic and outbreak solutions by max − l 2 ≤x≤ l 2 u(x) < u1 and max − l 2 ≤x≤ l 2 u(x) > u2, respectively. We will use the word “orbits” instead of “solutions” when the actual value of l > 0 is not specified. As in the non-advective case, the equilibria of (3.3) are given by zeros of φ(u), i.e. (u(x),v(x)) = (ui, 0) where 0 ≤ i ≤ 3. We will now determine their nature. Consider the Jacobian matrix of system (3.3): J(u,v) = [ 0 1 −φ′(u) q ] . Its trace and determinant are given by tr(J(u,v)) = q ≥ 0 and det(J(u,v)) = φ′(u). We observe the following regarding each equilibrium point of (3.3): • (u0, 0) = (0, 0): When q = 0, the the second order differential equation in (3.2) has a first integral (see Fact 2.2) and the phase plane portrait is symmetric with respect to the u-axis. In fact, when q = 0, the orbits of the system (3.3) are level curves of 1 2 v2 +F(u), where F(u) is given in the Fact 2.2. Thus, orbits are given by v = ±1 2 √ −F(u) + c, where c are constants. Note that PREVAILING WINDS AND SPRUCE BUDWORM OUTBREAKS: AN RDA MODEL 245 Figure 3.2. Phase plane portraits for system (3.3) when α = 0.154,β = 1.67, for q = 0.1, 0.2, 0.27, 0.33. The nature of saddle points remains the same, while both centers become unstable spirals. −F(·) has a local maximum at u = 0 (see Figure 2.3). It follows that in a neighborhood of (0, 0) the phase portrait of the system (3.3) consists of closed curves. Therefore, this equilibrium point is a center (see p. 296 in [8] for a similar argument). Note that φ′(0) = 1. When 0 < q < 2, we have 4 det(J(0, 0)) − tr2(J(0, 0)) = 4 − q2 < 0 while tr(J(0, 0)) = q > 0, so the equilibrium is an unstable spiral. When q ≥ 2, we have an unstable node. • (u1, 0): Since det(J(u1, 0)) = φ′(u1) < 0, this equilibrium point is a saddle for any value of q. • (u2, 0): Similarly to (0, 0), when q = 0, this equilibrium point is also a center (note that −F(·) has a local maximum at u = u2 as well). Note that φ ′(u2) > 0. When 0 < q < 2 √ φ′(u2), we have 4 det(J(u2, 0)) − tr2(J(u2, 0)) = 4(φ′(u2))2 − q2 < 0 while tr(J(u2, 0)) = q > 0, so the equilibrium is an unstable spiral. When q ≥ 2 √ φ′(u2), we have an unstable node. • (u3, 0): As in the case of (u1, 0), we have det(J(u3, 0)) = φ′(u3) < 0, and therefore this equilibrium point is a saddle point for any value of q. The changes in the phase plane portrait of system (3.3) as we increase the advection are illustrated in Figure 3.2. As noted above, for any value of q, the system (3.3) has exactly two saddle points: (u1, 0) and (u3, 0). We will denote the upper and lower branches of the stable and unstable manifolds of (ui, 0) (where i = 1 246 A. ANDERSON AND O. VASILYEVA or 3) by Sst+i ,S st− i ,S unst+ i and S unst− i , respectively. Here “+” corresponds to “upper” (first quadrant) and “-” corresponds to “lower” (fourth quadrant). Thus, Sst+i and S st− i have (ui, 0) as their ω-limit, while Sunst+i and S unst− i have the same point as their α-limit. Since the right hand side of (3.3) is given by C∞ functions of u and v, by the result of Guysinsky et al. [5], the homeomorphisms involved in the Grobman-Hartman theorem can be chosen to be differentiable at the corresponding equilibrium point and having Jacobian matrix equal to the identity matrix at that point. It follows that the lines through each saddle point spanned by the eigenvectors of the Jacobian matrix J(ui, 0) are tangent to their stable and unstable manifolds (see also Lemma 6.2 in [21]). Alternatively, the slopes of the stable and unstable manifolds of these equilibria can be found by viewing them as solution curves of an ordinary differential equation involving v as a function of u. Namely, since every solution (u(x),v(x)) of system (3.3) satisfies u′(x) = v(x), the portions of the solution curves of (3.3) located either in the upper half of the uv-plane {(u,v) : v > 0} or in the lower half {(u,v) : v < 0} satisfy the vertical line test. Such a curve located in either half of the uv-plane can be viewed as the graph of a function v = v(u) defined on a certain interval. Such a function will be a solution of the first order differential equation dv du = q − φ(u) v . (3.5) Lemma 3.1. The slopes of the tangent lines to the stable and unstable manifolds at the saddle point (ui, 0) (i = 1, 3) of (3.3) are given by m = q ± √ q2 − 4φ′(ui) 2 . Furthermore, the slope of Sst+i and S st− i is given by m = q− √ q2−4φ′(ui) 2 < 0 and the slope of Sunst+i and Sunst−i is given by m = q+ √ q2−4φ′(ui) 2 > 0. Proof. Let v = v(u) be the function whose graph is the upper portion of the stable or unstable manifold of the saddle point (ui, 0). Thus, v(u) is defined on an interval of the form (ui − ε,ui) or (ui,ui + ε), v(u) > 0 for all such u, limu→ui v(u) = 0, and v(u) satisfies (3.5) on its domain. Our goal is to find m = limu→ui v ′(u). Noticing that φ(ui) = 0 and limu→ui v(u) = 0 and using L’Hospital’s rule, we get: m = lim u→ui v′(u) = q − lim u→ui φ(u) v(u) = q − lim u→ui φ′(u) v′(u) = q − φ′(ui) m . The above equation can be written as m2 −qm + φ′(ui) = 0. Hence, the slope of the tangent line to the stable or unstable manifold at (ui, 0) is given by m = q± √ q2−4φ′(ui) 2 . The second statement of the lemma follows by noticing that du dx = v > 0 in the first quadrant and du dx = v < 0 in the fourth quadrant, and φ′(ui) > 0. � 4. The critical advection for outbreak orbits In this section we will study conditions for the existence of outbreak orbits in terms of the advection speed q. We establish existence of a critical value of advection beyond which the outbreak solutions do not exist. We obtain upper and lower bounds for this threshold value. PREVAILING WINDS AND SPRUCE BUDWORM OUTBREAKS: AN RDA MODEL 247 Figure 4.1. Behaviour of Sunst+1 : three possible cases. 4.1. Existence and upper bound for critical advection for outbreaks. We will begin by de- scribing the effect of increasing q on the behavior of the unstable manifold of the saddle point (u1, 0). Due to the fact that the vector field on the u-axis is given by ( du dx , dv dx ) = (0,−φ(u)), and φ(u) < 0 for u1 < u < u2 and u > u3 and φ(u) > 0 for u2 < u < u3, there are three possible types of behavior of the upper portion of the unstable manifold Sunst+1 of (u1, 0): • Case 1: Sunst+1 crosses the u-axis between u2 and u3; • Case 2: Sunst+1 becomes a saddle-to-saddle conection with the ω-limit at (u3, 0); • Case 3: Sunst+1 stays in the first quadrant without approaching any equilibrium. We will denote the point where Sunst+1 intersects u-axis (if any) as (u ∗, 0). To indicate the dependence of u∗ on q we will use the notation u∗ = u∗(q). We will write u∗(q) = u3 in Case 2, and u ∗(q) = ∞ in Case 3. We will also write u∗(q) < ∞ in Cases 1 and 2. See Figure 4.1 for an illustration of the three cases. Recall that the portion of the curve located in the first quadrant is the graph of a smooth function of u, which we denote v = v(u; q). It is defined for u1 ≤ u ≤ u∗(q) (or u ≥ u1 in Case (3)) and is positive and satisfies the differential equation (3.5) for u1 < u < u ∗(q) (or for u > u1 in Case (3)), and we have v(u1; q) = v(u ∗; q) = 0 (or just v(u1; q) = 0 in Case (3)). The following lemma shows monotonicity of v(u; q) and u∗(q) with respect to q. Lemma 4.1. Suppose 0 ≤ q1 < q2 and u∗(q2) < ∞. Then also u∗(q1) < ∞, u∗(q1) < u∗(q2), and for any u1 < u < u ∗(q1) we have v(u; q1) < v(u; q2). Proof. By Lemma 3.1, the slope of v(u; q) at u = u1 equals m(q) = q+ √ q2−φ′(u1) 2 , which is a strictly increasing funciton of q. Thus, m(q1) < m(q2). Hence, lim u→u+1 v(u; q1) −v(u1; q1) u−u1 = lim u→u+1 v(u; q1) u−u1 = m(q1) < m(q2) = lim u→u+1 v(u; q2) −v(u1; q2) u−u1 = lim u→u+1 v(u; q2) u−u1 . Then for some ε > 0, we have v(u; q1) < v(u; q2) for all u1 < u < u1 + ε. Next, we will show that this inequality holds for all values u > u1 for which both curves are in the first quadrant of the uv-plane, i.e. for all u1 < u < min(u ∗(q1),u ∗(q2)). Suppose that for some u1 < ũ < min(u ∗(q1),u ∗(q2)) we have v(ũ; q1) = v(ũ; q2) (> 0). We may assume that ũ is the smallest such value. Then for all u1 < u < ũ we have 0 < v(u; q1) < v(u; q2), 248 A. ANDERSON AND O. VASILYEVA which implies that d du v(ũ; q1) ≥ dduv(ũ; q2). But d du v(ũ; q1) = q1 − φ(ũ) v(ũ; q1) = q1 − φ(ũ) v(ũ; q2) < q2 − φ(ũ) v(ũ; q2) = d du v(ũ; q2), a contradiction. Thus, for all u1 < u < min(u ∗(q1),u ∗(q2)) we have (0 <)v(u; q1) < v(u; q2). Then we must have u∗(q1) ≤ u∗(q2). We claim that the inequality is strict. Indeed, suppose u∗(q1) = u∗(q2). Then for some u2 < u < u ∗(q1) close to u ∗(q1) we have d du v(u; q1) > d du v(u; q2). Since u > u2, we have φ(u) > 0. Then d du v(u; q1) = q1 − φ(u) v(u; q1) < q2 − φ(u) v(u; q1) < q2 − φ(u) v(u; q2) = d du v(u; q2), a contradiction. Thus, u∗(q1) < u ∗(q2), and for any u1 < u < u ∗(q1) we have v(u; q1) < v(u; q2), as needed. � Lemma 4.2. For any u1 < u < u2, v(u; q) is a continuous function of q ≥ 0. Proof. First, recall that φ(u) < 0 for all u1 < u < u2 Let 0 ≤ q1 < q2 and u1 < u < u2. By Lemma 4.1, 0 < v(u; q1) < v(u; q2). Thus, d du (v(u; q2) −v(u; q1)) = q2 −q1 + φ(u) ( 1 v(u; q1) − 1 v(u; q2) ) < q2 −q1. By the Mean Value Theorem applied to the function f(u) = v(u; q2) − v(u; q1) on the interval [u1,u], we have: 0 < v(u; q2) −v(u; q1) < (q2 −q1)(u−u1). It follows that for any q0 ≥ 0 we have lim q→q0 v(u; q) = v(u; q0), as needed. � The next lemma shows that for sufficiently large values of q, Case 3 takes place (i.e. u∗(q) = ∞), and moreover, v(u; q) is an increasing function of u ≥ u1. Lemma 4.3. Let φ∗ = maxu2≤u≤u3 φ(u) and q̄ = √ φ∗ u2−u1 . Then (i) for any q ≥ 0, v(u2; q) > q(u2 −u1); (ii) for any q ≥ q̄, v(u; q) is a nondecreasing function of u ≥ u1; (iii) for any q ≥ q̄ we have u∗(q) = ∞, and the system (3.3) has no outbreak orbits. Proof. (i) First, note that for any q ≥ 0 and u1 < u < u2 we have φ(u) < 0 and v(u; q) > 0, and therefore d du v(u; q) = q − φ(u) v > q. Let us define v(u1; q) = limu→u+1 v(u; q) = 0. Then by the Mean Value Theorem, for some u1 < ξ < u2, we have v(u2; q) = v(u2; q) −v(u1; q) = d du v(ξ; q)(u2 −u1) > q(u2 −u1). (ii) Suppose q ≥ q̄. Let g(u,v) = q − φ(u) v . Then g(u,v) ≥ 0 for any (u,v) in the quadrant u ≥ u2, v ≥ q(u2 −u1). PREVAILING WINDS AND SPRUCE BUDWORM OUTBREAKS: AN RDA MODEL 249 Indeed, for any u ≥ u2, we have φ(u) ≤ φ∗ (note that φ∗ > 0 and φ(u) < 0 for u > u3). Thus, for u ≥ u2, v ≥ q(u2 −u1) we have g(u,v) = q − φ(u) v ≥ q − φ∗ v ≥ q − φ∗ q(u2 −u1) = 1 q ( q2 − q̄2 ) ≥ 0, as needed. By (i), v(u; q) > q(u2 − u1) for u = u2. Next, we wil show that this inequality holds for all u ≥ u2. If not, let û > u2 be the smallest u > u2 such that v(û; q) = q(u2 − u1). Then, by the Mean Value Theorem, there exists u2 < ξ < û such that d du v(ξ; q) < 0. But (ξ,v(ξ; q)) belongs to the quadrant u ≥ u2, v ≥ q(u2−u1), and, therefore, dduv(ξ; q) = g(ξ,v(ξ; q)) ≥ 0, a contradiction. Thus, for any u ≥ u2, we have v(u; q) > q(u2 −u1) and hence dduv(u; q) = g(u,v(u; q)) ≥ 0 (in fact, ≥ 1 q (q2 − q̄2)). Combining the above with the fact that d du v(u; q) > q > 0 for all u1 < u < u2, as established int he proof of (i), we conclude that the function v(u; q) is nondecreasing for all u > u1. (iii) Clearly, v(u; q) > 0 for any q ≥ q̄ and u > u1, i.e. u∗(q) = ∞. Therefore, the boundary value problem (3.3, 3.4) has no solution whose orbit crosses the u-axis to the right from u1. Hence, there is no outbreak orbits for q ≥ q̄. � Let Γ = {q ≥ 0 : u∗(q) < ∞}. Thus, Γ is the set of all values of q for which either Case 1 or Case 2 occurs. By Lemma 4.1, Γ is an interval of the real line containing zero as its left endpoint. Lemma 4.3 shows that q̄ 6∈ Γ and q̄ therefore serves as an upper bound for Γ. Let us define q∗cr = sup Γ. Clearly, q∗cr ≤ q̄ = √ φ∗ u2−u1 . The following Proposition clarifies the role of q∗cr as a critical value of q for existence of outbreak orbits. Proposition 4.4. (i) For any 0 ≤ q < q∗cr we have u∗(q) < u3. (ii) For any q > q∗cr we have u ∗(q) = ∞. (iii) u∗(q∗cr) = u3 (i.e. for q = q ∗ cr, S unst+ 1 becomes a heteroclinic connection between (u1, 0) and (u3, 0) and coincides with S st+3). (iv) limq→q∗−cr u ∗(q) = sup u∗(Γ) = u3. Proof. (i) Suppose 0 ≤ q < q∗cr . Let q̃ = q+q∗cr 2 . Then q < q̃ < q∗cr and q, q̃ ∈ Γ. By Lemma 4.1, u∗(q) < u∗(q̃) ≤ u3, as needed. (ii) Follows by noticing that q > q∗cr = sup Γ implies q 6∈ Γ. (iii) Let ū be any value between u1 and u2, e.g. ū = u1+u2 2 . For any q ≥ 0, let (U(x; µ,ν,q),V (x; µ,ν,q)), x ≥ 0, be the solution (u(x),v(x)) of the system (3.3) such that (u(0),v(0)) = (µ,ν). Then by global con- tinuity with respect to the initial conditions and parameters, for any fixed x > 0, (U(x; µ,ν,q),V (x; µ,ν,q)) is continuous with respect to (µ,ν,q). Let (uq(x),vq(x)) = (U(x; ū,v(ū; q),q),V (x; ū,v(ū; q),q)), which defines the portion of Sunst+1 starting at (ū,v(ū; q)) for the given value of q. By Lemma 4.2, v(ū; q) is a continuous function of q. Then for any fixed x = L > 0, (uq(L),vq(L)) is continuous with respect to q. To show u∗(q∗cr) = u3 we need to eliminate the two other cases: u ∗(q∗cr) < u3 and u ∗(q∗cr) = ∞. Case (1): Suppose u∗(q∗cr) < u3. Then for q = q ∗ cr, the orbit S unst+ 1 crosses into the fourth quadrant, and therefore, for some L > 0 we will have vq∗cr (L) < 0. Let ε = −vq∗cr (L) > 0. For any q > q ∗ cr, we have q 6∈ Γ, and thus, the orbit Sunst+1 stays in the first quadrant. Therefore, for any x > 0 we will 250 A. ANDERSON AND O. VASILYEVA have vq(x) > 0. In particular, for any q > q ∗ cr, we have ‖(uq(L),vq(L)) − (uq∗cr (L),vq∗cr (L))‖≥ ε, contradicting the continuity of (uq(L),vq(L)) at q = q ∗ cr. Case (2): Suppose u∗(q∗cr) = ∞. Then for q = q∗cr, the orbit S unst+ 1 stays in the first quadrant and eventually crosses the vertical line u = u3. Thus, for some L > 0, we have uq∗cr (L) > u3. Let ε = uq∗cr (L) −u3 > 0. For any q < q ∗ cr, we have q ∈ Γ, and thus, the orbit S unst+ 1 stays to the left from the line u = u3. Therefore, for any x > 0 we will have uq(x) < u3. In particular, for any q < q ∗ cr, we have ‖(uq(L),vq(L)) − (uq∗cr (L),vq∗cr (L))‖≥ ε, contradicting the continuity of (uq(L),vq(L)) at q = q ∗ cr. (iv) In this proof we will use the notation (uq(x),vq(x)) as defined in the proof of (iii). Note that by monotonicity of u∗(q) (Lemma 4.1) and the fact that u∗(q) ≤ u3 for all q ∈ Γ, the set u∗(Γ) is bounded above by u3 and limq→q∗−cr u ∗(q) = sup u∗(Γ) ≤ u3. Suppose that sup u∗(Γ) < u3. Let ε = u3−sup u∗(Γ) 2 . Note that ε > 0. By (iii), we have uq∗cr (x) → u3 as x →∞. Thus, for some L > 0, we have uq∗cr (L) > u3 −ε. In particular, for any q < q ∗ cr, we have ‖(uq(L),vq(L)) − (uq∗cr (L),vq∗cr (L))‖≥ ε, and as in (iii), it contradicts the continuity of (uq(L),vq(L)) at q = q ∗ cr. Thus, sup u ∗(Γ) = u3, as needed. � 4.2. Lower bound for critical advection for outbreaks. We will now work towards obtaining a lower bound for q∗cr. Recall that v(u; q) satisfies the differential equation (3.5) for u1 < u < u ∗(q) and v(u1; q) = v(u ∗(q); q) = 0 (when q ≤ q∗cr). Multiplying both sides of (3.5) by v(u; q) and integrating from u1 to u ∗(q) gives: ∫ u∗(q) u1 (qv(u; q) −φ(u))du = ∫ u∗(q) u1 v(u; q) dv du du = ∫ v(u∗(q);q) v(u1;q) vdv = (v(u∗(q); q))2 2 − (v(u1; q)) 2 2 = 0 − 0 = 0. Therfore, we obtain the following: q ∫ u∗(q) u1 v(u; q)du = ∫ u∗(q) u1 φ(u)du. Since, ∫u∗(q) u1 v(u; q)du > 0, we get q = ∫u∗(q) u1 φ(u)du∫u∗(q) u1 v(u; q)du . (4.1) Note that φ(u) < 0 for u1 < u < u2 and φ(u) > 0 for u2 < u < u3. Let A1 = − ∫u2 u1 φ(u)du, A∗ =∫u∗(q) u2 φ(u)du and A2 = ∫u3 u2 φ(u)du. Then A1,A ∗,A2 > 0 and A ∗ ≤ A2. Notice that ∫u∗(q) u1 φ(u)du = A∗ − A1 and ∫u3 u1 φ(u)du = A2 − A1. Recall that for q = 0, outbreak orbits exist exactly when u∗(0) < u3. This is also equivalent to existence of outbreak orbits for some q ≥ 0, or for all sufficiently small q ≥ 0. Since the existence of an outbreak orbit for q = 0 is equivalent to F(u1) < F(u3), and F(u3)−F(u1) = A2 −A1, an outbreak orbit exists for q = 0 (or for sufficiently small q ≥ 0) if and only if A1 < A2. Below we give an alternative proof of this fact based on the phase plane argument. PREVAILING WINDS AND SPRUCE BUDWORM OUTBREAKS: AN RDA MODEL 251 Proposition 4.5. Suppose φ(u) has three positive zeros. Then outbreak orbits exist for sufficiently small q ≥ 0 if and only if A2 > A1. Proof. Suppose that for some q ≥ 0, an outbreak solution exists (Case 1 holds). Then u2 < u∗(q) < u3, and by (4.1) 0 ≤ q = A∗ −A1∫u∗(q) u1 v(u; q)du , and thus, A∗ ≥ A1. Note that since u∗(q) < u3, we have A2 > A∗. Therefore A2 > A1. Suppose that for q = 0, an outbreak solution does not exist (Case 2 or Case 3 holds), but A2 > A1. Then u∗(0) = u3 or u ∗(0) = ∞. We have: 0 > A1 −A2 = − ∫ u3 u1 φ(u)du = ∫ u3 u1 v dv du du = (v(u3)) 2 2 − (v(u1)) 2 2 = (v(u3)) 2 2 − 0 = (v(u3)) 2 2 , a contradiction. Thus, an outbreak solution exists. � Figure 4.2 shows the phase plane of system (3.3) for Q = 6.5 and q = 0 and three different values of β: (a) R = 0.588,A1 < A2, outbreak solutions exist (here R̃(Q) < R < R1(Q), i.e. we are in Ω ∗ or Region IV of the QR-plane); (b) R = 0.574,A1 = A2, no outbreak solution is possible (here R = R̃(Q), which corresponds to the lower boundary of Ω∗); (c) R = 0.571,A1 > A2, no outbreak solution is possible (here R2(Q) < R < R̃(Q), which corresponds to Region III). Recall that limq→q∗cr− u ∗(q) = u3. Passing to the limit as q → q∗cr − in (4.1), we get q∗cr = ∫u3 u1 φ(u)du∫u3 u1 v(u; q∗cr)du . (4.2) Note that the denominator of the fraction in (4.2) represents the area of the region in the uv-plane bounded above by the curve v = v(u; q∗cr), which is the orbit S unst+ 1 = S st+ 3 , where u1 ≤ u ≤ u3, and the u-axis. We will denote this orbit by S13. Our next objective is to find an upper bound for this area. Let umin ∈ [u1,u2] be such that φ′(umin) is the minimal value of φ′ on [u1,u2]. If φ is concave up on [u1,u2], then umin = u1, otherwise it is the inflection point 0 < û1 < 1 of φ(·) (see Lemma 2.1). The proof of the following proposition amounts to constructing a triangular “trapping region” enclosing the orbit S13 and having two of its vertices at the saddle points (u1, 0) and (u3, 0). Proposition 4.6. For q ≥ 0, let A(q) = ( √ q2 − 4φ′(umin) + q)( √ q2 − 4φ′(u3) −q)(u3 −u1)2 4( √ q2 − 4φ′(umin) + √ q2 − 4φ′(u3)) . Then ∫u3 u1 v(u; q∗cr)du < A(q ∗ cr). Proof. See Appendix A. � Using the above proposition and the fact that ∫u3 u1 φ(u)du > 0, we obtain: q∗cr = ∫u3 u1 φ(u)du∫u3 u1 v(u,q∗cr)du > ∫u3 u1 φ(u)du A(q∗cr) . 252 A. ANDERSON AND O. VASILYEVA Figure 4.2. The phase plane of system (3.3) for Q = 6.5,q = 0 and R = 0.588 (R̃(Q) < R < R1(Q)),R = 0.574 (R = R̃(Q)), R = 0.571 (R2(Q) < R < R̃(Q)). Graphs of φ(u) with the areas A1 and A2 indicated are shown on the right. Note that A(q) is well-defined, continuous and positive for all q ≥ 0. For q ≥ 0, let G(q) = q − ∫u3 u1 φ(u)du A(q) , PREVAILING WINDS AND SPRUCE BUDWORM OUTBREAKS: AN RDA MODEL 253 which is a continuous function for q ≥ 0. Then G(0) = − ∫u3 u1 φ(u)du A(0) < 0, and G(q∗cr) > 0. Let q > 0 be the smallest positive zero of G(q). Then 0 < q < q ∗ cr, that is, q serves as a lower bound for q∗cr. 5. Geometric analysis of endemic and outbreak orbits and steady state solutions In the first part of this section, we will use phase plane analysis to study conditions for existence of endemic and outbreak orbits. Then we will discuss the dependence of habitat length on the maximal density of a steady state solution. 5.1. Conditions for existence of endemic and outbreak orbits. First, we will show that for q ≥ 2 there is no positive solution (i.e. satisfying u(x) > 0 for − l 2 < x < l 2 ) to the boundary value problem (3.3, 3.4) for any l > 0, i.e. (3.3) has no endemic or outbreak orbits. Proposition 5.1. Suppose q ≥ 2 and l > 0. Then (3.3, 3.4) has no positive solution. In particular, (3.3) has no endemic or outbreak orbits. Proof. Note that a positive solution of (3.3,3.4) satisfies v ( − l 2 ) > 0 and v (x∗) = 0 for some − l 2 < x∗ < l 2 . Consider the half-line v = u, u > 0 in the uv-plane. To establish that there is no positive solution of (3.3,3.4) it suffices to show that the slope of vector field of the system (3.3) when restricted to this half-line is greater than 1. That is, qu−φ(u) u > 1 for all u > 0. We have φ(u) u = 1 − u Q − 1 R · u u2 + 1 < 1 ≤ q − 1, for all u > 0, hence qu−φ(u) u = q − φ(u) u > q −q + 1 = 1, as needed. � The next lemma demonstrates that there is a horizontal line such that once an orbit of (3.3) passes above this line it will stay above it. Recall that φ(u) < 0 outside the interval [0,u3], and thus φ(u) reaches its maximal value on R at some 0 ≤ u ≤ u3. Lemma 5.2. Let q > 0 and M = 1 q max0≤u≤u3 φ(u). Then for any solution (u(x),v(x)) of (3.3), if for some x∗ we have v(x∗) > M, then v(x) > M for all x > x∗. Proof. Suppose for some x̂ > x∗ we have v(x̂) ≤ M. By the Intermediate Value Theorem, there is x∗ < x̃ ≤ x̂ such that v(x̃) = M. We may assume that for all x∗ < x < x̃, v(x) > M. By the Mean Value Theorem, there exists x∗ < ξ < x̃ such that v′(ξ) < 0. But v(ξ) > M, and thus v′(ξ) = qv(ξ) −φ(u(ξ)) > qM − max 0≤u≤u3 φ(u) = 0, a contradiction. � We will now show the existence of an orbit located in the first quadrant of the uv-plane, connecting a point on the positive v-semiaxis and the saddle point (u1, 0) (in fact, a portion of S st+ 1 ). 254 A. ANDERSON AND O. VASILYEVA Figure 5.1. Region W used in the proof of Lemma 5.3. Lemma 5.3. Let 0 < q < 2. Then there exists a solution (u(x),v(x)), x ≥ 0, of (3.3) such that lim x→∞ (u(x),v(x)) = (u1, 0), (u(0),v(0)) = (0,v1) where 0 < v1 ≤ M, and (u(x),v(x)) ∈ (0,u1) × (0,M) for all x > 0. Proof. Let (u(x),v(x)), x ∈ R, be the solution of (3.3) such that lim x→∞ (u(x),v(x)) = (u1, 0) and v(x) > 0 as x → ∞. Thus, it represents the upper portion of the stable manifold of the saddle point (u1, 0). Note that u(x) is a non-decreasing function of x as long as v(x) ≥ 0. Consider the region W = [0,u1] × [0,M] of the uv-plane (see Figure 5.1). Note that as x → ∞, (u(x),v(x)) ∈ (0,u1) × (0,M) (the interior of W). Since the only other equilibrium point in W is (0, 0), as we decrease x, the only two possibilities are: (1) (u(x),v(x)) remains in W and approaches (0, 0) as x →−∞; (2) (u(x),v(x)) leaves W through (2a) the upper boundary of W (v = M, 0 < u < u1), (2b) the lower boundary of W (v = 0, 0 < u < u1), (2c) or the left boundary of W (u = 0, 0 < v ≤ M). The case (1) is not possible since (0, 0) is an unstable spiral, so there is no solution that approaches (0, 0) as x →−∞ entirely in the first quadrant. If (2a) holds then there is x∗ ∈ R such that v(x∗) = M and for some ε > 0, v(x) > M for all x∗ − ε < x < x∗, which contradicts Lemma 5.2. The case (2b) contradicts the fact that dv dx = −φ(u(x)) < 0 on the lower boundary. Thus, (2c) holds, i.e. there exists x∗ ∈ R such that (u(x∗),v(x∗)) = (0,v1) where 0 < v1 ≤ M, and for all x > x∗ we have (u(x),v(x)) ∈ W . Shifting by −x∗ gives us x∗ = 0, as needed. � Lemma 5.4. Let 0 < q < 2 and suppose u∗(q) < u3 (Case 1). Then there exists a solution (u(x),v(x)), x ≥ 0, of (3.3) such that lim x→∞ (u(x),v(x)) = (u3, 0), (u(0),v(0)) = (0,v3) where v1 < v3 ≤ M, and (u(x),v(x)) ∈ (0,u3) × (0,M) for all x > 0. Proof. Consider the portions of the stable and unstable manifolds of the saddle point (u1, 0) of (3.3) located in the first quadrant (i.e. the first quadrant portions of the curves Sst+1 and S usnt+ 1 ). Note that both curves can be viewed as graphs of continuous functions of u, namely, v = vl(u) defined for PREVAILING WINDS AND SPRUCE BUDWORM OUTBREAKS: AN RDA MODEL 255 Figure 5.2. Region W ′ = W1 ∪W2 ∪W3 used in the proof of Lemma 5.4. 0 ≤ u ≤ u1 (by Lemma 5.3) and v = vr(u) defined for u1 ≤ u ≤ u∗(q). By Lemma 5.3, vl(u) > 0 for all 0 ≤ u < u1, and vl(u1) = 0. We also have vr(u1) = vr(u∗(q)) = 0 and vr(u) > 0 for u1 < u < u∗(q). By Lemma 5.2, vl(u),vr(u) ≤ M throughout their respective domains. Consider the following region of the uv-plane: W ′ = W1 ∪W2 ∪W3, where W1 = {(u,v) : 0 ≤ u < u1,vl(u) ≤ v ≤ M}, W2 = {(u,v) : u1 ≤ u < u∗(q),vr(u) ≤ v ≤ M} and W3 = [u ∗(q),u3] × [0,M]. Note that W ′ is a closed simply connected plane region (see Figure 5.2). Let (u(x),v(x)), x ∈ R, be the solution of (3.3) such that limx→∞(u(x),v(x)) = (u3, 0) and v(x) > 0 as x → ∞ (it represents part of the upper portion of the stable manifold Sst+3 of the saddle point (u3, 0)). Recall that u(x) is a non-decreasing function of x as long as v(x) ≥ 0. Note that as x → ∞, (u(x),v(x)) ∈ (u∗(q),u3) × (0,M) (the interior of W3). Since the only other equilibrium point in W ′ is (u1, 0), as we decrease x, the only two possibilities are: (1) (u(x),v(x)) remains in W ′ and approaches (u1, 0) as x →−∞; (2) (u(x),v(x)) leaves W through (2a) the upper boundary of W ′ (v = M, 0 < u < u3), (2b) the lower boundary of W ′ (v = vl(u), 0 ≤ u < u1, or v = vr(u), u1 ≤ u ≤ u∗(q), or v = 0, u∗(q) < u < u3), (2c) or the left boundary of W ′ (u = 0,v1 < v ≤ M). The case (1) is not possible since (u1, 0) is a saddle point having v = v r(u) as its unstable manifold in W ′ (thus, no other orbit in W ′ can approach it as x →−∞). The case (2a) impossible for the same reason as in the proof of Lemma 5.3. The case (2b) is impossible due to the uniqueness of the solution of an initial value problem (for lower boundaries of W1 or W2), or due to the fact that dv dx = −φ(u(x)) < 0 on the lower boundary of W3 (u ∗(q) < u < u3,v = 0). Thus, (2c) holds, i.e. there exists x ∗ ∈ R such that (u(x∗),v(x∗)) = (0,v3) where v1 < v3 ≤ M, and for all x > x∗ we have (u(x),v(x)) ∈ W ′. Shifting by −x∗, we may assume that x∗ = 0, as needed. 256 A. ANDERSON AND O. VASILYEVA Figure 5.3. Region Tim used to “capture” the curve S unst− i . � In order to establish the existence of endemic and outbreak orbits, we will now turn our attention to the fourth quadrant of the uv-plane. Let i = 1 or 3. Recall that (ui, 0) is a saddle point, and by Lemma 3.1 the slope of Sunst−i (the lower branch of its unstable manifold) is given by q+ √ q2−4φ′(ui) 2 . Our goal is to find a closed triangular region Tim = {(u,v) : 0 ≤ u ≤ ui, m(u−ui) ≤ v ≤ 0} bounded by positive u-semiaxis, negative v-semiaxis and the line through (ui, 0) with a positive slope m, such that the portion of Sunst−i in the fourth quadrant will lie entirely in that region (see Figure 5.3). To guarantee that the small portion of this curve near the saddle point is located inside this region, we will need m > q+ √ q2−4φ′(ui) 2 . To guarantee that Sunst−i does not cross the hypothenuse of T i m, we will require that the vector field of the system (3.3), when restricted to the line segment v = m(u−ui), 0 < u < ui, forms an acute angle with the normal vector (−m, 1) pointing towards the origin. Thus, the condition is (−m, 1) ·(v,qv−φ(u)) > 0 where v = m(u−ui) and 0 < u < ui. Equivalently, we need m(u−ui)(q −m) −φ(u) > 0, or m2 −mq > − φ(u) u−ui . This inequality has to hold for all 0 < u < ui. For i = 1 or 3, define a function hi : [0,ui] → R as follows hi(u) = { − φ(u) u−ui , u < ui −φ′(ui), u = ui. PREVAILING WINDS AND SPRUCE BUDWORM OUTBREAKS: AN RDA MODEL 257 Clearly, hi is continuous on [0,ui], and thus, it attains its maximal value, h i max(> 0) in [0,ui]. Note that by the definition of himax, we have − φ(u) u−ui ≤ himax for all 0 < u < ui. Hence, it suffices to require m2 −mq −himax > 0. Let mi be the positive root of the quadratic equation m 2 −mq −himax = 0, i.e. mi = q + √ q2 + 4himax 2 . Clearly, m2 −mq −himax > 0 holds for m > mi. We are now ready to analyze the behavior of the curves Sunst−1 , S unst+ 1 and S unst− 3 in the fourth quadrant of the uv-plane. Lemma 5.5. Suppose 0 ≤ q < 2. The fourth quadrant portion of the curve Sunst−1 is contained in the region T 1m1 . It connects its α-limit (u1, 0) with a point (0,v − 1 ) where −m1u1 < v − 1 < 0. Proof. Since T 1m1 = ⋂ m>m1 T 1m, it suffices to prove the first statement for all m > m1. We have m > m1 ≥ q+ √ q2−4φ′(u1) 2 . Hence, near the point (u1, 0), the curve S unst− 1 stays above the line v = m(u−u1) (as its tangent at (u1, 0) has slope < m). Since m > m1, by the vector field arguments preceding this lemma, the curve cannot cross this line at any point (u,m(u−u1)) where 0 < u < u1. Therefore the curve is contained in T 1m for any m > m1, and thus, also in T 1 m1 . The second statement follows easily by noticing that Sunst−1 cannot approach the origin (an unstable equilibrium) and cannot cross the u-axis (as the vector field is directed downward on the u-axis between u = 0 and u = u1). � For q ≥ 0, let Sst+1 (q) and S st− 1 (q) denote the curves S st+ 1 and S st− 1 for this particular value of q. Lemma 5.6. Suppose 0 < q < 2. The fourth quadrant portion of the curve Sunst+1 (q) is contained in the region T 3m3\D, where D is the region bounded by the curve S unst− 1 (q), the positive u-semiaxis and negative v-semiaxis. The curve Sunst+1 connects (u ∗(q), 0) with the point (0,v∗) where −m3u3 < v∗ < v−1 . Proof. As in the previous lemma, it suffices to prove the first statement for T 3m for any m > m3. Let (u(x),v(x)),x ∈ R be a solution of (3.3) representing Sunst+1 (q). Thus, limx→−∞(u(x),v(x)) = (u1, 0) and for some x∗, we have (u(x∗),v(x∗)) = (u∗(q), 0). As we increase x > x∗, u(x) will decrease and, by Poincare-Bendixson theorem, the point (u(x),v(x)) either leaves the region T 3m or approaches an equilibrium point in T 3m (note that by the Bendixson-Dulac criterion, the system (3.3) has no closed orbits for q > 0, as the vector field has positive divergence q). Recall that the vector field of the system (3.3) points downward on the u-axis for 0 < u < u1 and u2 < u < u3. Note that S st− 1 (0) is the lower portion of the right-side homoclinic orbit of (u1, 0) (for q = 0), it connects the point (u∗(0), 0) to (u1, 0) in the fourth quadrant, and is a graph of a function of u. Recall that for any q ≥ 0 and point (u,v) where v 6= 0, the slope of the solution curve of (3.3) passing through (u,v) is given by dv du = q − φ(u) v . Then for any q > 0, the vector field of (3.3) on the curve Sst−1 (0) will point outside the region above the curve. In addition, by Lemma 3.1, the slope of the curve Sst−1 (q) at (u1, 0) is given by m(q) = q − √ q2 −φ′(u1) 2 , and for q > 0, m(q) > m(0). 258 A. ANDERSON AND O. VASILYEVA Figure 5.4. Regions D and T 3m3\D from the proof of Lemma 5.6. Thus, we conclude that (u(x),v(x)) cannot cross the u-axis again for any x > x∗, and cannot approach (u1, 0) as x → ∞. Clearly, it cannot approach (u2, 0) either. Note that Sunst+1 (q) cannot cross the curve Sunst−1 (0) or the segment v = m(u−u3), 0 < u < u3 (by the vector field argument). Since there are no more equilibria in the region T 3m\D, the only possibility for Sunst+1 (q) is to leave the region T 3 m\D (and thus, the fourth quadrant) through the v-axis between v = v−1 and −m3u3. This proves the lemma (see Figure 5.4 for an illustration). � Lemma 5.7. Suppose 0 < q < 2. The fourth quadrant portion of the curve Sunst−3 is contained in the region T 3m3\E, where E is the region bounded by the curve S unst+ 1 , the positive u-semiaxis and negative v-semiaxis. It connects its α-limit (u3, 0) with a point (0,v − 3 ) where −m3u3 < v − 3 < v ∗. Proof. The proof is similar to the proof of the previous two lemmas. � See Figure 5.5 for an illustration of Lemma 5.7. Remark 5.8. Note that when q = 0, we have Sunst+1 = S st− 1 forming a loop symmetric with respect to the u-axis. Note also that Lemma 5.5 includes the case q = 0, and Lemma 5.7 can be modified to include the case q = 0. Moreover, because of the symmetry of the phase plane about the u-axis, this gives us the existence of orbits described in Lemmas 5.3 and 5.4 for the non-advective case, which agrees with the analysis of the steady state solutions in the non-advective case done in [11]. We are now ready to state our main results. Theorem 5.9. Suppose (Q,R) ∈ Ω∗. Let u1 < u2 < u3 be the three positive zeros of φ(u). (i) Suppose 0 < q < min(2,q∗cr). Then there exist v − 3 < v ∗ < v−1 < 0 < v1 < v3 with the following properties: PREVAILING WINDS AND SPRUCE BUDWORM OUTBREAKS: AN RDA MODEL 259 Figure 5.5. Regions E and T 3m3\E from the proof of Lemma 5.7. (1) Sst+1 connects (0,v1) with (u1, 0) (in the first quadrant); (2) Sunst−1 connects (u1, 0) with (0,v − 1 ) (in the fourth quadrant); (3) Sunst+1 connects (u1, 0) with (0,v ∗), passing through the first and the fourth quadrants and crossing the u-axis at the point (u∗(q), 0) where u2 < u ∗(q) < u3; (4) Sst+3 connects (0,v3) with (u3, 0), passing above S st+ 1 and S unst+ 1 (in the first quadrant); (5) Sunst−3 connects (u3, 0) with (0,v − 3 ), passing below S unst+ 1 (in the fourth quadrant); (6) for any 0 < µ < u1 there exists l = l(µ; q) > 0 and a solution (u(−),v(−)) : [ − l 2 , l 2 ] → R2 of (3.3, 3.4) such that u(x) > 0 for all x ∈ ( − l 2 , l 2 ) , 0 < v ( − l 2 ) < v1 and v − 1 < v ( l 2 ) < 0, and (u(x∗),v(x∗)) = (µ, 0) for some x∗ ∈ ( − l 2 , l 2 ) (endemic orbits); (7) for any u∗(q) < µ < u3 there exists l = l(µ; q) > 0 and a solution (u(−),v(−)) : [ − l 2 , l 2 ] → R2 of (3.3, 3.4) such that u(x) > 0 for all x ∈ ( − l 2 , l 2 ) , v1 < v ( − l 2 ) < v3 and v − 3 < v ( l 2 ) < v∗, and (u(x∗),v(x∗)) = (µ, 0) for some x∗ ∈ ( − l 2 , l 2 ) (outbreak orbits); 260 A. ANDERSON AND O. VASILYEVA Figure 5.6. Phase plane in the case when both endemic and outbreak orbits exist. (8) there are no other positive solutions of (3.3, 3.4) for any l > 0; (9) endemic orbits provide monotone bijections between intervals (0,v1) and (v − 1 , 0) of the v-axis via the interval (0,u1) of the u-axis; (10) outbreak orbits provide monotone bijections between intervals (v1,v3) and (v − 3 ,v ∗) of the v-axis via the interval (u∗(q),u3) of the u-axis. (ii) Suppose q∗cr < 2, and q ∗ cr < q < 2. Then there exist v − 1 < 0 < v1 such that (i)(1,2,6,8,9) hold. In this case, Sunst+1 lies entirely in the first quadrant and passes above the point (u3, 0), and the system does not admit outbreak orbits. (iii) Suppose q ≥ 2. Then (3.3, 3.4) has no positive solutions for any l > 0 (no endemic or outbreak orbits). Proof. Parts (i)(1-5) follow from Lemmas 5.3, 5.5, 5.6 and 5.7. Parts (i)(6-10) follow by noticing the direction of the vector field of (3.3) on both u and v-axes, and the fact that the stable and unstable manifolds form “trapping regions” in the uv-plane. See Figure 5.6. The proof of part (ii) is similar to that of part (i). See Figure 5.7. Part (iii) is given by Proposition 5.1. � Remark 5.10. For (Q,R) is Region III (that is, Ω\Ω∗, see Figure 2.4), the function φ(u) still has three positive roots, and for q ≥ 0, the unstable manifold Sunst+1 remains in the first quadrant of the PREVAILING WINDS AND SPRUCE BUDWORM OUTBREAKS: AN RDA MODEL 261 Figure 5.7. Phase plane in the case (ii) of Theorem 5.9 (when only endemic orbits exist). Figure 5.8. Phase plane when (Q,R) 6∈ Ω and q < 2. uv-plane. Thus, there are no outbreak orbits in this case. The equilibrium (u1, 0) remains a saddle point, and the phase portrait is similar to the one in Figure 5.7. Therefore, the same argument as the one used in the proof of Theorem 5.9(ii) shows that for q < 2, there are endemic orbits with maximal density given by values 0 < µ < u1. For q ≥ 2, the conclusion of Theorem 5.9(iii) holds as well. When (Q,R) is outside of Ω (that is, in Regions I, II, V or VI ), the function φ(u) has only one positive zero, we will still denote it u1, and we have φ(u) > 0 for 0 < u < u1 and φ(u) < 0 for u > u1. In the uv-plane, the point (u1, 0) is still a saddle point (for any q ≥ 0). In this case, the situation is again completely analogous to Theorem 5.9(ii, iii). See Figure 5.8 for a typical orbit when q < 2. 5.2. Habitat length vs. maximal density. Now that we have established conditions for existence of orbits representing positive solutions of (3.3, 3.4) for some l > 0, we can focus on the function l = l(µ; q) 262 A. ANDERSON AND O. VASILYEVA that provides the length of the habitat corresponding to such an orbit passing through the point (µ, 0) on the u-axis; we can refer to l(µ; q) as the parametric length of the orbit. Recall, that the function l(µ; q) is defined for 0 ≤ q < 2 and 0 < µ < u1 or, in the case the conditions of Theorem 5.9(i) are satisfied (i.e. (Q,R) ∈ Ω∗ and q < q∗cr), u∗(q) < µ < u3. Recall from Section 2, that in the non-advective case, there is an explicit integral formula for l(µ; 0) derived in [11], however, the first integral method does not apply in the advective case. We can compute l(µ; q) as follows. Let (uµ(x),vµ(x)) be the solution of the initial value problem consisiting of the system (3.3) and the condition (uµ(0),vµ(0)) = (µ, 0), defined on R. Let l+(µ; q) be the smallest positive solution of uµ(−x) = 0 (the parametric length of the upper part of the orbit) and l−(µ; q) be the smallest positive solution of uµ(x) = 0 (the parametric length of the lower part of the orbit). Then l(µ; q) = l+(µ; q) +l−(µ; q). Note that (uµ(x),vµ(x)) is a continuous function of x and µ, and therefore, both l+(µ; q) and l−(µ; q) are continuous functions of µ (as solutions of uµ(−x) = 0 and uµ(x) = 0). Thus, l(µ; q) is also a continuous function of µ. In the non-advective case, the behavior of l(µ; 0) was described in [11], and is summarized in Section 2 according to the choice of (Q,R) is one of the five regions of the QR-plane (see Figure 2.4). As we increase q, we expect that the general shape of l(µ; q) as a function of µ for q > 0 will remain of the same four types outlined in Figure 2.5, but with some vertical and horizontal transformations and transitions from one type of shape to another. Of the main interest to us is the behavior of the critical domain lengths for the endemic and outbreak solutions, lc1(q) and l c 2(q), respectively, as we vary q. Recall that in all the cases, lc1 = limµ→0 l(µ; 0) = π. In the advective case, linearizing (3.1) at zero steady state, we can find lc1(q) = lim µ→0 l(µ; q) = 2π√ 4 −q2 . Note that lc1(q) increases with q and lim q→2 lc1(q) = ∞. It has been shown in [11] that lc2(0) > l c 1(0). We conjecture that the same inequality holds for 0 ≤ q < min(q∗cr, 2), and moreover, q ∗ cr < 2. In the next section, we will use numerics to investigate the effect of advection on the behavior of l(µ; q) and the nature of steady state solutions. 6. Numerical simulations In this section, we will use numerics to test our approximations of the critical advection q∗cr for existence of outbreak orbits. For different choices of (Q,R) ∈ Ω∗, we obtain the values of q∗cr and its lower and upper estimates q and q̄, respectively. We will also explore the dependence of the habitat length l(µ; q) on the maximal population density µ of a steady state solution for different value of advection speed q. We will perform these numerical simulations for (Q,R) in Regions II, IV and V (see Figure 2.4). 6.1. Estimating critical advection. We test our lower and upper estimates for q∗cr as follows. To obtain the lower estimate q, recall that q is the smallest positive root of the equation G(q) = 0, where G(q) = q − ∫ u3 u1 φ(u)du A(q) and A(q) is given by (A.1). Recall that he upper estimate q̄ is given by q̄ = √ φ∗ u2 −u1 PREVAILING WINDS AND SPRUCE BUDWORM OUTBREAKS: AN RDA MODEL 263 where φ∗ = maxu2≤u≤u3 φ(u). Finally, we can approximate the actual value of q ∗ cr by using Euler’s method to obtain the solution curve v = v(u; q) of the differential equation dv du = q − φ(u) v with the property lim u→u+1 v(u) = 0, which is the unstable manifold Sunst+1 . Since we know that the slope of S unst+ 1 at (u1, 0) is given by m = q+ √ q2−4φ′(u1) 2 (see Lemma 3.1), we can use this information to perform the first step of the Euler’s method to avoid division by zero. We determine whether v(u; q) intersects the u-axis (i.e. whether u∗(q) < ∞). Increasing the advection q by small increments, we iterate until we reach the situation when u∗(q) = ∞ (the stopping criterion is when v(u; q) takes a value above the threshold M = 1 q max 0 ≤ u ≤ u3φ(u)). The following table summarizes the numerical results for estimating critical advection. All values of (Q,R) are chosen in Region IV (Ω∗). Q R q q∗cr q̄ 6.5 0.6 0.22 0.325 0.5507 10 0.5556 0.48 0.705 1.6122 10 0.5128 0.3 0.455 0.7454 15 0.5128 0.53 0.855 1.6031 15 0.4 0.3 0.475 0.7687 Figures B.1 - B.5 illustrate the lower and upper estimates for q∗cr for different choices of parameters. We observe that q∗cr < 2 for each choice of (Q,R), and q < q ∗ cr < q̄. 6.2. Dependence of habitat length on maximal density and advection. To investigate the effect of advection on the habitat length l = l(µ; q) we choose (Q,R) = (15, 0.1) (Region II), (Q,R) = (15, 0.5128) (Region IV or Ω∗), and (Q,R) = (15, 0.7) (Region V). In each case we consider three values of q: q = 0, 0.35, 0.7 for Region IV (recall that in this case q∗cr = 0.855) and q = 0, 0.7, 1.4 for regions II and V. In each case, we observe that l = l(µ; q) is monotonously increasing with respect to q. Moreover, we make the following observations in each of the three examples. • (Q,R) = (15, 0.1) (Region II, see Figure B.6): As we increase q, the curve l = l(µ; q) (in the (µ,l)-plane) retains its monotonously increasing concave up shape, but moves upward. As a result, for a fixed l > π, we see that the maximal density µ of the unique positive steady state decreases, and for q > 2 √ 1 − π2 l2 only the the extinction steady state left. • (Q,R) = (15, 0.5128) (Region IV, see Figure B.7): As we increase q, the left component of the curve l = l(µ; q), defined for 0 > µ < u1, retains its monotonously increasing concave up shape, but moves upward. The right component of the curve l = l(µ; q), defined for u∗(q) < µ < u3, retains its concave up shape with two vertical asymptotes. The vertical asymptote µ = u∗(q) is moving to the right as we increase q (thus reducing the domain of l = l(−; q)). The right component moves upward and is compressed horizontally. The minimal value lc2(q) of l = l(µ; q) on the interval (u∗(q),u3) is increasing accordingly. Note that for q = q ∗ cr = 0.855, 264 A. ANDERSON AND O. VASILYEVA the two vertical asymptotes collapse into one, and the right component of the curve l = l(µ; q) disappears. As a result, for a fixed l > π we see that the density of the endemic equilibrium is decreasing, and if l > lc2(q), the density of the threshold steady state is increasing and the density of the outbreak state is decreasing. At some value of q < q∗cr, the threshold and outbreak states are lost, and only endemic and extinction states are left. For q > 2 √ 1 − π2 l2 only the extinction steady state is left. In this example, lc2(q) > l c 1(q), as expected. • (Q,R) = (15, 0.7) (Region V, see Figure B.8): We observe that for all values of q, the curve l = l(µ; q) is defined for 0 < µ < u1 and has a vertical asymptote at µ = u1. As in the previous cases, it moves upward as we increase q. For smaller values of q it still has two internal extrema: the local maximum value lc3(q) and the local minimum l c 2(q). However, for larger values of q, l(µ; q) becomes an increasing function of µ, similar to the situation in Region VI. As a result, for a fixed π < l < lc2(0) we observe only an endemic state with its maximal density decreasing with q. For relatively small values of l > lc2(0), we can have the following scenarios as we increase q: (1) outbreak (large density) only with decreasing maximal density; (2) three states (endemic, threshold, outbreak) with endemic and outbreak maximal densities decreasing and threshold maximal density increasing; (3) endemic (small density) only, with decreasing maximal density; moreover, depending on l, the following transitions are possible as we increase q: (2 → 3), (1 → 2 → 3) or (1 → 3). As before, for q > 2 √ 1 − π2 l2 only the extinction steady state is left. 7. Discussion In this paper we have considered an advective version of the classical non-dimensional reaction- diffusion model for population dynamics of spruce budworm introduced in [11]. Biologically, the goal was to understand the possible effect of biased movement caused by prevailing winds on the spatial distribution and occurrence of SBW outbreaks. As in [11], we have used hostile boundary conditions; an appropriate setting could be that of an island or a peninsula (e.g. Newfoundland, Cape Breton, Anticosti) that experiences relatively strong prevailing winds and is known to have occasional SBW outbreaks, such as the 1970’s outbreak on the island Newfoundland, as well as the current (2021) increase in SBW numbers observed on the island [14, 16]. Mathematically, we explored the new RDA setting where the reaction term φ(u) = u− 1 Q u2 − 1 R u2 1 + u2 combines the logistic growth and the loss due to predation, and results in richer dynamics than in the logistic RDA setting and allows multiple steady states. Here Q and R are non-dimensional biological parameters: Q represents the carrying capacity of the tree branches and R characterizes the generalist predation impact (increasing R leads to decrease in predation). The main focus of this study was on the effect that advection can have on outbreaks as opposed to the endemic states. The original nonspatial spruce budworm model introduced in [12] focused on the case of three positive equilibria 0 < u1 < u2 < u3 given by the zeros of the reaction term φ(u), where u1 was the stable endemic equilibrium, u2 was the unstable threshold, and u3 was the stable outbreak equilibrium. We identify the endemic and outbreak solutions of the corresponding RDA model given by ∂u ∂t = ∂2u ∂x2 −q ∂u ∂x + φ(u) PREVAILING WINDS AND SPRUCE BUDWORM OUTBREAKS: AN RDA MODEL 265 (where q is the advection speed) with hostile boundary conditions as positive steady state solutions, where the maximal density of the endemic solution is bounded above by u1, and the maximal density of the outbreak solution is bounded below by u2. As in [22] and [21], here we employed the geometric approach to the study of steady state solutions of the RDA model, by representing them as orbits in the phase plane a system of first order ODEs{ du dx = v dv dx = qv −φ(u) obtained by setting ∂u ∂t = 0 and letting v = ∂u ∂x . The orbits that originate on the positive v-semiaxis and terminate at the negative v-semiaxis are viewed as steady state solutions of the original RDA equation subject to the hostile (Dirichlet) conditions on some finite domain [0, l]. If the orbit intersects the u-axis at the point (µ, 0), we can identify µ as the maximal density of the corresponding steady state. One can then distinguish endemic and outbreak steady state solutions by noticing that the former are represented by orbits with 0 < µ < u1 while the latter correspond to u2 < µ < u3. The resulting system has four equilibria: (0, 0), (u1, 0), (u2, 0), (u3, 0). In the non-advective case considered in [11], the phase portrait is symmetric about the u-axis, equilibria (0, 0) and (u2, 0) are centers, while (u1, 0) and (u3, 0) are saddle points. In this case, there is an explicit formula for habitat length l(µ) as a function of maximal density µ. Once we introduce advection, the phase plane portrait changes: it is no longer symmetric, and the two centres become unstable spirals for relatively small values of q, and become unstable nodes with the further increase of q. The other two equilibria remain saddle points. The first integral method that lead to the explicit formula for l(µ) in the non-advective case is no longer applicable for q > 0. The key element of our analysis was the fact that when the non-advective model admits an outbreak solution, the increase of advection forces the u-intercept (u∗, 0) of the upper branch of the unstable manifold of (u1, 0) to move towards (u3, 0), and at a certain critical value of advection speed q ∗ cr, we observe the appearance of a heteroclinic orbit connecting (u1, 0) and (u3, 0). As a consequence, outbreak orbits can only exist for q < q∗cr, as they have to cross the u-axis between u ∗ and u3. After establishing the existence of this critical advection value for outbreaks, we obtained upper and lower bounds for q∗cr, expressed in terms of the non-dimensional biological parameters Q and R. We also observed that for q ≥ 2, no positive steady states exist on any finite domain. Based on numerical evidence, we conjecture that q∗cr < 2 for any choice of the biological parameters Q and R for which φ(u) has three positive zeros and non-advective model admits an outbreak solution. For 0 ≤ q < min(2,q∗cr), we showed, using phase plane techniques such as trapping regions, that the RDA model has both endemic and outbreak steady state solutions for some finite domains, while for q satisfying min(2,q∗cr) < q < 2 (if any) only endemic solutions remain. Our main objective in this paper was to establish and estimate the critical value of advection for existence of outbreak solutions on some finite domains. Further exploration of the advective model for spruce budworm will require an analysis of the behaviour of the habitat length l(µ; q) as the function of both the maximal density µ and the advection speed q. In this paper, we made the first step in this direction by exploring l(µ; q) numerically and making some observations regarding its behaviour. In particular, we observe that l(µ; q) is increasing with respect to q and the critical domain size for outbreak solutions is always greater than that of the endemic solution (the latter given by lc1(q) = 2π√ 4−q2 ). Most of our analysis focused on the case when in the non-advective setting we have three positive steady states (endemic, threshold and outbreak). In that case, the numerics suggest that when we fix the habitat size and let the advection increase, the outbreak and threshold solutions always disappear first, followed by disappearance of the endemic state. This agrees well with our analytical results indicating 266 A. ANDERSON AND O. VASILYEVA that increasing advection beyond q∗cr prevents existence of outbreak solutions on any finite domain while still admitting endemic (low denisity) states, provided q∗cr < 2. We have also observed that in the case when the non-advective setting only allows an outbreak (high-density) solution, the increase of advection may lead to appearance of three states, replaced by a single endemic (low-density) state. Exploring this case analytically is an interesting direction for future research. While our focus in this paper was on the model for SBW in a “windy island” habitat, there are many other advective settings where the reaction term that we considered would be appropriate. The same model can be applied to other defoliating insect species (e.g. gypsy moth), but more importantly, it is also applicable to logistically growing aquatic organisms subject to generalist predation. In this case, different boundary conditions may have to be imposed. A typical upstream boundary condition will be the zero flux condition which, in our non-dimensionized setting, can be visualized by the half line v = qu in the first quadrant of the uv-plane. The typical downstream conditions considered in such models are either hostile (u = 0) or outflow (v = 0). In either case, we can still talk about outbreak orbits versus endemic orbits, but the situation with existence of outbreak orbits is more complicated: even if the unstable manifold of (u1, 0) stays in the first quadrant, there still might be orbits connecting the half line v = qu and the half line u = 0 or v = 0 representing the downstream boundary condition. This setting will require further exploration and additional phase-plane analysis. Other settings that our analysis can be applicable to are the effect of shifting habitat boundaries due to climate change on population dynamics as in [17, 1] and sinking phytoplankton as in [7], if one assumes the presence of a generalist predator. Acknowledgements. The second author would like to thank Andre Arsenault, Jean-Noel Can- dau, Frithjof Lutscher and Sebastien Portalier for valuable discussions. The authors also thank the anonymous referee for valuable comments and suggestions. Appendix A. Proof of Proposition 4.6 Consider the line passing through (u1, 0) with the slope m > 0. Thus, its equation is given by v = m(u − u1). In order to guarantee that the line is above the curve S13 given by v = v(u; q∗cr) for u1 ≤ u ≤ u3, it suffices to require that the slope of the line is greater than the slope of the direction field of the differential equation (3.5) at every point on the line for u1 ≤ u ≤ u3, i.e. m ≥ dv du |v=m(u−u1) = ( q∗cr − φ(u) v ) |v=m(u−u1) = q ∗ cr − φ(u) m(u−u1) , for every u1 ≤ u ≤ u3. Note that by Lemma 3.1, the slope of the tangent line to v = v(u; q∗cr) at u = u1 is given by m1 = q∗cr+ √ (q∗cr) 2−4φ′(u1) 2 . Clearly, we must have m ≥ m1. Recall that φ(u) ≥ 0 for u2 ≤ u ≤ u3. Thus, since m1 > q ∗ cr, for any u2 ≤ u < u3 and any m ≥ m1, we have m ≥ q∗cr ≥ q∗cr − φ(u) m(u−u1) . It remains to consider the case when u1 < u < u2. By the Mean Value Theorem, for any u1 < u < u2 there exists u1 < ũ < u such that φ(u) = φ(u) −φ(u1) = φ′(ũ)(u−u1). Thus, we need q − 1 m φ′(ũ) ≤ m, for all u1 < u < u2. Recall that umin ∈ [u1,u2] is such that φ′(umin) is the minimal value of φ′ on [u1,u2]. Then for any u1 < ũ < u2 and m > 0 we have q − 1 m φ′(ũ) ≤ q∗cr − 1 m φ′(umin). PREVAILING WINDS AND SPRUCE BUDWORM OUTBREAKS: AN RDA MODEL 267 Thus, it suffices to choose the smallest m > 0 such that q∗cr − 1 m φ′(umin) ≤ m. Such m is given by m = q∗cr + √ (q∗cr) 2 − 4φ′(umin) 2 . Note that m ≥ m1 (and m = m1 if umin = u1). Therefore, the line v = q∗cr + √ (q∗cr) 2 − 4φ′(umin) 2 (u−u1) lies above the curve S13. We will now follow similar argument to construct another side of a triangular “trapping region” for S13 with the base given by the segment [u1,u3] of the u-axis. Namely, we are looking for a line v = k(u−u3) that stays above S13. Similar to the previous case, it suffices to ensure that the slope of the line is less than the slope of the direction field of the differential equation (3.5) at every point on the line where u1 ≤ u ≤ u3, i.e. k ≤ dv du |v=k(u−u3) = ( q∗cr − φ(u) v ) |v=k(u−u3) = q ∗ cr − φ(u) k(u−u3) , for all u1 ≤ u ≤ u3. Note that by Lemma 3.1, the slope of the tangent line to v = v(u) at u = u3 is given by m3 = q∗cr− √ (q∗cr) 2−4φ′(u3) 2 . Note that we need k ≤ m3. Recall that φ(u) ≤ 0 for u1 ≤ u ≤ u2. Also note that m3 < 0 and u−u3 < 0. Thus, for any u1 ≤ u ≤ u2 and any k ≤ m3, we have k ≤ 0 ≤ q − φ(u) k(u−u3) . It remains to consider the case when u2 < u < u3. By the Mean Value Theorem, for any u2 < u < u3 there exists u < ũ < u3 such that φ(u) = φ(u) −φ(u3) = φ′(ũ)(u−u3). Thus, we need q − 1 k φ′(ũ) ≥ k, for all u2 < u < u3. By Lemma 2.1, φ(·) has only two positive inflection points: 0 < û1 < û2. Clearly, û1 < u2 and φ′(û2) > 0. Therefore, the mimimum of φ ′(u) on [u2,u3] occurs at u = u3. Then for any u2 < ũ < u3 and k < 0 we have q − 1 k φ′(ũ) ≥ q∗cr − 1 k φ′(u3). Thus, it suffices to choose the largest k < 0 such that q∗cr − 1 k φ′(u3) ≥ k. Such k is given by k = q∗cr − √ (q∗cr) 2 − 4φ′(u3) 2 = m3. Therefore, the line v = q∗cr+ √ (q∗cr) 2−4φ′(u3) 2 (u−u3) lies above the curve S13 (given by v = v(u), u1 ≤ u ≤ u3). Hence, we have constructed a triangular “trapping region” for the curve S13 with the base given by the segment connecting (u1, 0) and (u3, 0) and the sides given by v = m(u−u1) and v = m3(u−u3) (the latter is the tangent line to S13 at (u3, 0)). Note that if S13 is concave down (which seems to be the case according to the plots), or if û1 < u1 (e.g., if u1 ≥ 1), m is the slope of the tangent line to S13 at (u1, 0) as well. Let θ1 = arctan m1 and θ2 = arctan m3. Then we have tan θ1 = m1 = q∗cr + √ q∗cr 2 − 4φ′(umin) 2 268 A. ANDERSON AND O. VASILYEVA Figure A.1. Estimating ∫u3 u1 v(u)du using a triangular trapping region. and tan θ2 = m3 = −q∗cr + √ q∗cr 2 − 4φ′(u3) 2 . Thus, the area of the triangular region, is given by A = 1 2 (b1 + b2)h = 1 2 (u3 −u1)h. To find h, we note that h = b1 tan θ1 = b2 tan θ2 = (u3 −u1 − b1) tan θ2. Solving for b1, we get b1 = (u3 −u1) tan θ2 tan θ1 + tan θ2 , and, therefore, h = (u3 −u1) tan θ1 tan θ2 tan θ1 + tan θ2 . Hence, A = 1 2 · (u3 −u1)2 tan θ1 tan θ2 tan θ1 + tan θ2 = A(q∗cr), where A(q) = ( √ q2 − 4φ′(umin) + q)( √ q2 − 4φ′(u3) −q)(u3 −u1)2 4( √ q2 − 4φ′(umin) + √ q2 − 4φ′(u3)) (A.1) is viewed as a function of q ≥ 0. By the construction of the triangular region, we conclude that∫ u3 u1 v(u; q∗cr)du < A(q ∗ cr). Appendix B. Figures for Section 6 PREVAILING WINDS AND SPRUCE BUDWORM OUTBREAKS: AN RDA MODEL 269 Figure B.1. Biological parameter values: Q = 6.5,R = 0.6. Left: the graph of G(q) with the values q = 0.22 (circle), q∗cr = 0.325 (star) and q̄ = 0.5507 (square). Right: orbits Sunst+1 for q ≥ q with increments of ∆q = 0.01. Roots u0 = 0 < u1 < u2 < u3 of φ(u) are marked with stars. Figure B.2. Biological parameter values: Q = 10,R = 0.5556. Left: the graph of G(q) with the values q = 0.48 (circle), q∗cr = 0.705 (star) and q̄ = 1.6122 (square). Right: orbits Sunst+1 for q ≥ q with increments of ∆q = 0.01. Roots u0 = 0 < u1 < u2 < u3 of φ(u) are marked with stars. Figure B.3. Biological parameter values: Q = 10,R = 0.5128. Left: the graph of G(q) with the values q = 0.3 (circle), q∗cr = 0.455 (star) and q̄ = 0.7454 (square). Right: orbits Sunst+1 for q ≥ q with increments of ∆q = 0.01. Roots u0 = 0 < u1 < u2 < u3 of φ(u) are marked with stars. 270 A. ANDERSON AND O. VASILYEVA Figure B.4. Biological parameter values: Q = 15,R = 0.5128. Left: the graph of G(q) with the values q = 0.53 (circle), q∗cr = 0.855 (star) and q̄ = 1.6031 (square). Right: orbits Sunst+1 for q ≥ q with increments of ∆q = 0.05. Roots u0 = 0 < u1 < u2 < u3 of φ(u) are marked with stars. Figure B.5. Biological parameter values: Q = 15,R = 0.4. Left: the graph of G(q) with the values q = 0.3 (circle), q∗cr = 0.475 (star) and q̄ = 0.7687 (square). Right: orbits Sunst+1 for q ≥ q with increments of ∆q = 0.05. Roots u0 = 0 < u1 < u2 < u3 of φ(u) are marked with stars. Figure B.6. Habitat length l vs. maximal density µ for Q = 15,R = 0.1 and q = 0, 0.7, 1.4. PREVAILING WINDS AND SPRUCE BUDWORM OUTBREAKS: AN RDA MODEL 271 Figure B.7. Habitat length l vs. maximal density µ for Q = 15,R = 0.5128 and q = 0, 0.35, 0.7. The positive zeros u1 < u2 < u3 of φ(u) are marked by stars. Figure B.8. Habitat length l vs. maximal density µ for Q = 15,R = 0.7 and q = 0, 0.7, 1.4. 272 A. ANDERSON AND O. VASILYEVA References [1] H. Berestycki, O. Diekmann, C.J. Nagelkerke, P.A. Zegeling, Can a species keep pace with a shifting climate? Bul. Math. Biol. 71 (2009), 399–429. [2] Y. Du, B. Lou, R. Peng, M. Zhou, The Fisher-KPP equation over simple graphs: varied persistence states in river networks, J. Math. Biol. 80 (2020), 1559–1616. [3] H. Gu, B. Lou, Spreading in advective environment modelled by a reaction diffusion equation with free boundaries, J. Diff. Equations 260 (2016), no.5, 3991– 4015. [4] W. Hao, K.-Y. Lam, Y. Lou, Ecological and evolutionary dynamics in advective environments: Critical domain size and boundary conditions, Disc. & Cont. Dyn. Systems - B 26 (2021), no.1, 367–400. [5] M. Guysinsky, B. Hasselblatt, V. Rayskin, Differentiability of the Hartman–Grobman linearization, Disc. Contin. Dyn. Syst. 9 (2003), no.4, 979–984. [6] Y. Jin, R. Peng, J. Shi, Population Dynamics in River Networks, J Nonlinear Sci. 29 (2019), 2501–2545. [7] T. Kolokolnikov, C. Ou, Y. Yuan, Phytoplankton depth profile and their transitions near the critical sinking velocity, J. Math. Biol. 59 (2009), no.1, 105–122. [8] M. Kot, Elements of Mathematical Ecology, Cambridge University Press, Cambridge, 2001. [9] K.-Y. Lam, S. Liu, Y. Lou, Selected Topics on Reaction-Diffusion-Advection Models from Spatial Ecology, Math. Appl. Sci. Eng. 1 (2020), no.2, 91–206. [10] K.-Y. Lam, Y. Lou, F. Lutscher, The Emergence of Range Limits in Advective Environments, SIAM J. Appl. Math. 76(2016), no.2, 641–662. [11] D. Ludwig, D.G. Aronson, H.F. Weinberger, Spatial patterning of the bruce budworm, J. Math. Biol. 8 (1979), no.3, 217–258. [12] D. Ludwig, D. Jones, C. Holling, Qualitative Analysis of Insect Outbreak Systems: The Spruce Budworm and Forest, J. Animal Ecology 47 (1978), no.1, 315–332. [13] F. Lutscher, E. Pachepsky, M. Lewis, The effect of dispersal patterns on stream populations, Appl. Math. 47 (2005), no.5, 749–772. [14] I.S. Otvos, B.H. Moody, The spruce budworm in Newfoundland: history, status and control, Fisheries and Environ- ment Canada, Canadian Forest Service, Newfoundland Forest Research Centre, St. Johns, NF, Information Report N-X-150 (1978). [15] E. Pachepsky, F. Lutscher, R. Nisbet, M. Lewis, Persistence, spread and the drift paradox, Theor. Popul. Biol. 67 (2005), 61–73. [16] Parks Canada – Gros Morne National Park – Spruce Budworm: https://www.pc.gc.ca/en/pn-np/nl/grosmorne/decouvrir-discover/sb, accessed July 12, 2021. [17] A.B. Potapov, M.A. Lewis, Climate and competition: the effect of moving range boundaries on habitat invasibility, Bull. Math. Biol. 66 (2004), no.5, 975–1008. [18] J. M. Ramirez, Population persistence under advection-diffusion in river networks, J. Math. Biol. 65 (2012), 919–942. [19] J. Sarhad, R. Carlson, K. E. Anderson, Population persistence in river networks, J. Math. Biol. 69 (2014) , 401–448. [20] D.C. Speirs, W.S.C. Gurney, Population persistence in rivers and estuaries, Ecology 82 (2001), no.5, 1219–1237. [21] O. Vasilyeva, Population dynamics in river networks: analysis of steady states, J. Math. Biol. 79 (2019), 63–100. [22] O. Vasilyeva, F. Lutscher, Population dynamics in rivers: analysis of steady states, Can. Appl. Math. Q. 18 (2010), no.4, 439–469. [23] Y. Wang, J. Shi, J. Wang, Persistence and extinction of population in reaction-diffusion-advection model with strong Allee effect growth, J. Math. Biol. (2019), 2093–2140. Department of Mathematics and Statistics, University of Guelph, 50 Stone Road East, Room 437 Mac- Naughton Building, Guelph, ON N1G 2W1, Canada E-mail address: aander20@uoguelph.ca Corresponding author, Memorial University of Newfoundland, Grenfell Campus, Corner Brook, NL A2H 5G4, Canada E-mail address: ovasilyeva@grenfell.mun.ca https://www.pc.gc.ca/en/pn-np/nl/grosmorne/decouvrir-discover/sb 1. Introduction 2. Preliminaries: an overview of the classical non-spatial and spatial spruce budworm models 2.1. The non-spatial Ludwig-Jones-Holling model. 2.2. The spatial Ludwig-Aronson-Weinberger model. 3. Model set up 4. The critical advection for outbreak orbits 4.1. Existence and upper bound for critical advection for outbreaks. 4.2. Lower bound for critical advection for outbreaks. 5. Geometric analysis of endemic and outbreak orbits and steady state solutions 5.1. Conditions for existence of endemic and outbreak orbits 5.2. Habitat length vs. maximal density 6. Numerical simulations 6.1. Estimating critical advection 6.2. Dependence of habitat length on maximal density and advection 7. Discussion Appendix A. Proof of Proposition 4.6 Appendix B. Figures for Section 6 References