Moment Methods for Advection on Networks and an Application to Forest Pest Life Cycle Models

We describe the network of coupled computational domains via the following mathematical framework. Let $\mathcal{G}=(\mathcal{V},\mathcal{E})$ be a finite directed graph with vertices $\mathcal{V}$ and edges $\mathcal{E}\subset\mathcal{V}\times\mathcal{V}$ such that if $(u,v)\in\mathcal{E}$, then $(v,u)\not\in\mathcal{E}$, i.e., between any pair of vertices, flow is allowed in one direction only. On each graph vertex $v\in\mathcal{V}$, we consider a density function $\rho^{v}(a,t)$, defined on the domain $\Omega:=(0,1)\times(0,\infty)$, that satisfies the governing linear advection-diffusion-reaction equation

subject to the boundary conditions

Here $F_{\text{in}}(t)=\nu(t)\rho(0,t)-\xi(t)\frac{\partial\rho}{\partial a}(0,t)$ is the influx into the domain at $a=0$ and it is assigned to equal $f_{\text{in}}(t)$, where $f_{\text{in}}$ is a given non-negative function of time. Note that in the diffusion-free case ($\xi^{v}=0$), the boundary condition at $a=1$ automatically becomes inactive.

For the first variable in $\rho^{v}(a,t)$ the letter $a$ is used, because in the SLF application (see §2.3) it represents the developmental age. Without loss of generality it is scaled to be $a\in[0,1]$ (in the application, $a$ represents the fraction of development achieved in a given life stage). The variable $t$ is time. In (1), the coefficients $\nu^{v}$, $\xi^{v}$, and $\mu^{v}$ are assumed non-negative at all times, and we further assume that on each domain, $\xi^{v}\leq\mathcal{O}(\nu^{v})$, i.e., as the advective speed $\nu$ approaches zero, the diffusion does so at the same rate (or faster). In the SLF application, $\nu$ is the developmental speed, $\xi$ is a spreading of age distributions (see §2.3), and $\mu$ is mortality. Suitable initial conditions $\rho^{v}(\cdot,0)$ are prescribed on all domains $v\in\mathcal{V}$.

With the governing equations and boundary conditions defined on each domain, the domains/vertices are now coupled together (along the edges $\mathcal{E}$) via fluxes. The outflux out of a domain (through $a=1$) is

where $(\xi\rho_{a})$-term vanishes due the outflow boundary condition in (2). The coupling conditions are defined as follows. For a given vertex/domain $v\in\mathcal{V}$, let

denote the set of vertices that lead into $v$, respectively the set of vertices that flow out of $v$. If $I_{v}$ is empty, then $v$ is a source node, and an influx function $f_{\text{in}}$ is to be prescribed. For every edge $(v,w)\in\mathcal{E}$, a “split ratio” or “flux fraction” $\alpha_{v,w}>0$ must be defined. If fluxes across domains exactly conserve mass, then one has for all $v\in\mathcal{V}$ that $\sum_{w\in O_{v}}\alpha_{v,w}=1$. However, it is also possible that not all the mass gets transferred across domains (e.g., mortality during life stage transitions), or that the mass multiplies (e.g., flux from egg-layers to eggs). Based on the split rations, the inflow into the domain $v$ is then given by

A critical aspect is that the boundary/coupling conditions are formulated to preserve advective causality, even if there exists some diffusive “smearing” in a domain. This is motivated by the application described in §2.3: in the age development model, there is no backwards flow of information, that is: developmental age never decreases. The diffusive term models the empirical phenomenon that age distributions are observed to widen as development progresses, see §2.3.

The above formulation avoids diffusive fluxes ($\xi\frac{\partial}{\partial a}\rho$) between domains. Moreover, because $\xi\leq\mathcal{O}(\nu)$, there is no outflux if there is no aging. One should also note that the model, as first described in [30], was formulated in terms of fractional steps: first, an advective sub-step is conducted, with coupled inflows and outflows, followed by a sub-step only pure diffusion and decay with zero flux at both boundaries. In the limit of the time step going to zero, the fractional steps and the model here are equivalent. Finally, the model in [30] also allows for an influx into a domain to result from the solution in another domain via a kernel, see below.

A key building block for the network flows defined above is a single transition between two domains connected via a single edge $(v,w)\in\mathcal{E}$, i.e., $O_{v}=\{w\}$ and $I_{w}=\{v\}$. We also assume $\alpha_{v,w}=1$. For this scenario, let us call $v$ the sending domain and $w$ the receiving domain. Here the receiving influx equals the sending outflux, $F^{w}_{\text{in}}=F_{\text{out}}^{v}$. Furthermore, we consider the sending influx to be some prescribed function $f^{v}_{\text{in}}$, and the receiving outflux $F_{\text{out}}^{w}$ is irrelevant.

To understand the key mechanism caused by flux couplings at varying advection speeds, let us also consider the case of no diffusion, $\xi^{v}=0=\xi^{w}$, and no decay, $\mu^{v}=0=\mu^{w}$. Intuitively, this situation is like two conveyor belts transporting sand, with the first belt dropping sand onto the second, and sand grains tending to pile perfectly on top of each other. If the sending speed $\nu^{v}$ equals the receiving speed $\nu^{w}$, the density profile passes across the domain boundary unmodified. If the sending speed is slower than the receiving speed, i.e., $\nu^{v}<\nu^{w}$, the profile widens upon domain transition. Conversely, if the sending speed is larger, i.e., $\nu^{v}>\nu^{w}$, the profile sharpens, i.e., it compresses in the $a$-direction and increases in the $\rho$-direction.

In the limiting situation of zero receiving speed, i.e., $\nu^{w}=0$ while $f^{w}_{\text{in}}>0$, the influx boundary condition (2), which here would read as $F^{w}_{\text{in}}=\nu^{w}\rho^{w}(0,t)=f^{w}_{\text{in}}$, cannot be satisfied in the sense of functions. However, in many applications such a situation makes perfect practical sense, for instance: eggs keep getting produced, but no egg development occurs; or: incoming sand keeps piling up at the start of a stopped receiving conveyor belt. It is therefore reasonable to allow for solutions in the sense of distributions, i.e., allow for the creation of a Dirac delta distribution at $a=0$. Clearly, the accurate representation of such solutions is numerically challenging; and one key advantage of moment methods is that they face no fundamental limitations in capturing this Dirac delta limit. In this work, we conduct a regularization by replacing Dirac peaks by narrow Gaussians (see §3.2). However, we do allow these Gaussians to be extremely narrow. Hence, for applications in which there is a lack of widening diffusive effects, moment methods that are asymptotic preserving in this described regime are of critical importance, as described in §4.2.

For many species, in the modeling of its population dynamics it is advantageous to break the population into life-stages where the populations’ dynamics (such as aging or death) are uniform within each life stage. For many insects such life stages come naturally due to their life cycle in instars, separated by molting events. Based on these life stages one can formulate a model for how a population may develop over time with a stage-structured system of age-structured PDE. This is a well-known methodology for modeling forest pest as well as other biological systems [19]. These systems often assume a form similar to equation (1) but with varying flux coupling conditions between stages. Models of this kind have been calibrated for many forest pests [11, 12, 33]. In this work we focus on a specific model, calibrated for the spotted lanternfly [30]. In §6, this model is used as a test problem for the moment method developed in §4.1, to assess the method’s utility in simulating the life cycle of forest pests efficiently.

In the model from [30], the SLF life cycle is organized into four ($v\in\{1,2,3,4\}$) life-stage domains: three egg stages and one motile stage. The reason for this division is that in each domain, individuals share the same response to temperature which drives their development and mortality, i.e., the parameters $\nu$, $\xi$, and $\mu$ in (1) are independent of the age variable $a$. However, the parameters do depend on time because they are functions of the ambient temperature $T$, i.e., for each domain there are functions $\nu^{v}(T)$, $\xi^{v}(T)$, and $\mu^{v}(T)$ (see [30]), and if a time-varying ambient temperature profile $T(t)$ is given, this effectively induces time-dependent coefficients in (1). The four stages, depicted in Fig. 1, are:

1.  $\bullet$

   Domain 1: “Non-diapause eggs”. Eggs develop normally, with the advection speed $\nu$ increasing with temperature (within certain limits [30]), and have normal resistance to extreme temperatures (expressed via the mortality function $\mu(T)$).

2.  $\bullet$

   Domain 2: “Diapause eggs”. Diapause is a state of slowed or paused development in exchange for improved cold resistance [24] in which warmth suppresses the advection speed $\nu$.

3.  $\bullet$

   Domain 3: “Post-diapause eggs”. Eggs’ development is driven by temperature as in domain 1. However, an increased cold resistance from the diapause state is retained.

4.  $\bullet$

   Domain 4: “Motiles”. This combines all non-egg SLF: nymphs, instars, and adults, including egg-laying adults. At the motile stage, the history of the (non-)diapause pathway is lost, and the response to temperature is uniform.

The transitions between life stages is modeled via the flux coupling conditions, depicted in in Fig. 1 via arrows. All depicted nonzero $\alpha_{v,w}=1$, except for the arrows emanating from the motile motile domain: the outflow (3) from domain 4 goes nowhere (death due to senescence); instead, eggs are produced via an egg laying kernel assigned to egg domains in a time-dependent fashion, see below.

On each domain $v\in\{1,2,3,4\}$, the function $\rho^{v}(a,t)$ represents the age density of females in that life stage. The stage-structured PDE assumes the form of (1) with standard flux coupling boundary conditions for the fluxes out of domains 1, 2, and 3:

The flux from domain 4 into domains 1 or 2 slightly deviates from the structure of §2.1 as the production of eggs not be determined by the outflux of domain 4, but rather by an egg-laying kernel as follows.

Here $I(t)$ is an indicator function that is $I(t)=1$ over the duration of the year when eggs are laid into non-diapause and $I(t)=0$ when eggs are laid into diapause. In [30], $I(t)$ is modeled based on the photo-period of the day when the eggs was laid, with $I=0$ after summer solstice but before winter solstice, and $I=1$ otherwise. The kernel $k(a)$ describes the rate of egg production of a female at a given motile developmental age $a$. Its integral $\int_{0}^{1}k(a)\,\mathrm{d}a$ is the total number of eggs produced by a female if it lives until death by senescence. The constant $\beta\in[0,1]$ captures that there is some background mortality $1\!-\!\beta$ not related to temperature.

The assumption $\xi^{v}\leq\mathcal{O}(\nu^{v})$ introduced in §2.1, which justifies the “causal” flux couplings across domains, is satisfied in the SLF model as follows. The ratio $\xi^{v}/\nu^{v}$ may differ across domains, but it is always independent of temperature. This expresses the fact that age distributions of SLF observed in the field widen in a diffusive fashion, plausibly caused by the real-world fact that at different times, different individuals develop/age at different rates, even when exposed to the same circumstances. Here, it is important to emphasize that for cold-blooded species, developmental age is not identical to the time spent alive.

A key feature of the ecological model studied here is the presence of multiple pathways that individuals can take through the life stage network. The diapause pathway equips eggs with an increased resistance to cold, which can be crucial to enable enough eggs to survive the cold season in temperate climates. However, there is an additional ecological effect of diapause: it tends to synchronize the generational cycle with the annual temperature. Without diapause, eggs deposited in late summer will start developing and, given a warm enough fall, may hatch just went winter comes which will kill all hatched motiles. Conversely, diapause prevents the possibility of multiple generations per year, so in extremely warm climates, the non-diapause pathway may in turn be advantageous. Model simulations, as conducted in §6, are therefore critical towards understanding the life cycle mechanics and predicting the establishment potential of invasive species like SLF.

Ecological stage- and age-structured PDE models have been explored computationally prior, e.g., in [14, 34, 18], however, only via standard numerical methods (which track many numerical degrees of freedom per domain) and not exploring the significant model reduction offered by a well-designed moment method. Another important novelty is that the specific model from [30] provides pathways of co-existence of diapause and non-diapause that have not been widely explored numerically. The structure of the model, particularly the transitions between domains involving diapause, provides precisely the mechanism that tends to produce Gaussian density profiles, as described in §2.2, at least in a wide range of relevant temperature profiles (see the results in §6).

Due to the structure of the governing equations described in §2.1, the considered models have a tendency to produce solutions that are close to Gaussians (that may decay and widen as they travel through the domains). Because a Gaussian is uniquely defined by three parameters, we choose moment methods that track exactly three moments on each domain in the network. Specifically, we track the lowest three monomial moments of the density function $\rho(a,t)$ on the domain $a\in(0,1)$, which are:

As the results for the SLF application in §6 show, the assumption of the density on each domain being close to a Gaussian is indeed well satisfied in the current establishment region in the United States (see Fig. 12), where the life cycle is dominated by diapause which synchronized species development and thus creates very peaked distribution functions $\rho$.

Closely associated with the monomial moments (6) are the derived key quantities: (i) the mass $M$ of species in the domain; (ii) the mean (developmental age) $E$ in the domain; and (iii) the variance (of developmental age) $V$ in the domain, which are given in terms of the monomial moments as:

Computationally, we choose to track the monomial moments $(m_{0},m_{1},m_{2})$, because their governing moment equations (see below) are easier. However, modulo the singularity at $M\to 0$, a description in terms of $(M,E,V)$ would be equivalently possible.

Numerical methodologies that represent the density function $\rho$ via pointwise values or cell averages (cf. §4.3) tend to require hundreds of degrees of freedom per domain. In comparison, tracking only three moments represents a significant reduction in dimensionality, which—if successful—can lead to substantial reductions in computational demand. In addition, there is a conceptual advantage of descriptions in terms of moments: the quantities total population $M$, mean developmental age $E$, and age variance $V$ are very intuitive and familiar concepts for ecologists and stakeholders working with these pests. While those moments could of course also be computed (post-hoc) for traditional methods like finite volumes, a reduced model and software that more directly track/update moments have the potential to be more accessible.

We now derive the equations that govern the rates of change of the moments (6). We differentiate (6) with respect to $t$, invoke the governing PDE (1), employ integration-by-parts, and/or cancel/substitute quantities via the boundary conditions (2). For notational efficiency here we suppress domain indices and the time variable. In that abridged notation, (2) reads as $(\nu\rho-\xi\rho_{a})(0)=f_{\text{in}}$ and $\xi\rho_{a}(1)=0$. With that, we obtain:

While some of the terms in the ODEs’ right hand sides can be written in terms of the moments themselves, other terms require the knowledge of the actual solution $\rho$ of the PDE (1). As the goal is to formulate a numerical method that tracks (on each domain) solely the moments $(m_{0},m_{1},m_{2})$, a reconstruction mapping is needed that assigns to a given vector of moments a reconstructed function. Such mappings are constructed in §3.2. Then, in §4 the components are combined into two versions of closed moment methods that approximate the original network PDE problems.

Careful attention must be paid to the reconstruction of the Gaussian from the moments. The reconstruction mapping is necessary to evaluate the fluxes, and for many applications, the reconstructed population function is the quantity of interest. The construction herein is a special case of the Hausdorff moment problem [31, 36], with a particular attention to the fact that exactly 3 moments are considered. We wish to design methods to represent a Gaussian on each domain, parameterized by $C$, $a_{0}$, and $\sigma$:

However, the moments $m_{0}$, $m_{1}$, and $m_{2}$ are being tracked. Hence, the mapping

and its inverse

must be understood and numerically approximated. In particular, the methods constructed in §4 will need the reconstruction mapping

The domain of this mapping (and others presented in this section) is a non-trivial issue, and will be discussed in §3.2.1.

Due to the homogeneity of the mappings $L$ and $R$, it suffices to characterize the mappings between $(a_{0},\sigma)$ and the statistical moments $(E,V)$, where $C$ may be chosen as a convenient non-negative constant:

and its inverse

As neither the forward nor backwards map may be described with elementary functions, the mappings are implemented via a lookup table, moving the bulk of the computational complexity to a pre-computation that is conducted only once. Here we discuss some of the concerns that arise in the creation of the lookup table, and potential avenues to alleviate them.

We now combine the governing equations for the moments (7) with the (regularized) reconstruction mapping (9) to obtain the system of ODEs, in which for each domain $v\in\mathcal{V}$ one has

where we use the short notation for the reconstruction $\rho^{v}(a)=\mathcal{R}(m_{0}^{v},m_{1}^{v},m_{2}^{v})(a)$. Moreover, via (3) the outflux from domain $v$ is $F_{\text{out}}^{v}=\nu^{v}\rho^{v}(1)$, again using the reconstruction. The outfluxes from all domains, via the flux coupling conditions (4), determine the influxes $f_{\text{in}}^{w}$ for all domains $w\in\mathcal{V}$.

This yields a $3n$-dimensional ODE system, where $n$ is the number of domains. Within each domain, the three moments are coupled via the reconstruction mapping $\mathcal{R}$, and the equations across domains couple via (4) from outfluxes to influxes. An important fact is that the ODE system (17) is nonlinear, even though the underlying network PDE model is linear (see the discussion in §4.3). Due to the regularized reconstruction (16), the right hand side of (17) is always well-defined. However, it is not a-priori guaranteed that the solution of (17) remain close to the true moments of the solution of the PDE model. The numerical results in §5 and §6 demonstrate the success and possible modes of failure of the moment method.

It is also possible that a given initial condition for the PDE (1) is not Gaussian, or even not realizable. Hence, we project the moments of any initial condition via (15) before starting the computation.

A key structural advantage of the moment method approach presented herein is that it generates a generic ODE system which in principle can be advanced via any time stepping method of choice. In particular, existing methodologies and toolboxes for high-order-in-time, adaptive time step choice, dense output, or even implicit time stepping [38] could easily be employed. Moreover, the structure of (17) can even be amenable to positivity-preserving ODE solvers [3]. In turn, if no such specialized ODE solver is used, time-stepping errors could produce non-realizable moments. Of particular concern is the possibility that $m_{0}$ could become negative in a domain. In Runge-Kutta schemes, the constraint $m_{0}\geq 0$ can be ensured via an ad-hoc flux limiting applied in each stage: any outflux (3) is capped to $F_{\mathrm{out}}^{v}=\min(\nu^{v}(t)\rho^{v}(1,t),m_{0}^{v}(t)/\Delta t)$, where $\Delta t$ is the time step. This ensures that the mass that leaves a domain never exceeds the mass that is inside the domain.

In all numerical tests conducted in §5 and §6, the use of a simple RK4 with uniform time steps and the above flux limiting generated positive solutions; moreover, for reasonably small time steps, it even retained realizability, even though that property is not guaranteed. That is in contrast to automated methods like Matlab’s ode45.m without flux limiting. For sharply peaked Gaussians, negative $m_{0}$-values frequently emerged; and even when that did not happen, the adaptive time stepper sometimes terminated without success because the target tolerance could not be met. For the same peaked Gaussian problems, RK4 with flux limiting always generated reasonable solutions, albeit not always at a desirable accuracy, see the results in Fig. 8.

With only three degrees of freedom per domain, the dimensionality of the ODE system (17) is likely as low as a reasonable reduced model of the full network PDE model can be, and the accuracy of this description is enabled by the problem-specific mechanism described in §2.2. That being said, low-dimensionality does not automatically translate into low computational cost. One right hand side evaluation of (17) requires (at least) one reconstruction per domain, whose fast evaluation relies on an efficient lookup table, created in §3.2.

As described in §2.2, the network advection problems studied herein possess a natural mechanism to create strongly peaked distributions. These pose computational challenges for the ODE-based moment method from §4.1. Here we develop a numerical methodology that remedies some of these challenges. Since the problem is fully rooted in the advection part of (1), in this section we restrict the formulation to (1) without diffusion or decay, i.e., $\xi$ and $\mu$ are zero.

Consider the case of a narrow Gaussian of width $\sigma\ll 1$ traveling in a domain of speed $\nu=\mathcal{O}(1)$ to the domain’s outflow boundary. This induces the characteristic time scale $\tau=\sigma/\nu$, on which solution values change at fixed $a$-locations. In order to accurately resolve this time scale, non-specialized time-stepping methods for the ODE-based moment methods tend to require time step sizes $\Delta t\leq\tau$. Hence, when $\tau$ is very small, as is the case when $\sigma\ll 1$ and $\nu=\mathcal{O}(1)$, traditional time-stepping methods become computationally expensive. A desirable numerical scheme tracks very narrow Gaussians without incurring the above restriction. Such numerical schemes that also work well when $\Delta t\gg\tau$, but sufficiently small for any accuracy requirements, are denoted asymptotic preserving [20]. The asymptotic preserving property is critical particularly in the case when flow of mass in a network of domains where the advection speed $\nu^{v}(t)$ in each domain $v$ can vary, creating Dirac delta distributions that must be accurately tracked. The AP scheme here tracks the mass that is “pushed” from one domain to another during one time step, effectively shifting the temporal rates of change (as in the the ODE-based method in §4.1) to rates of change in the $a$-variable that can be resolved via direct formulas or simple quadratures.

To derive the method, let $\rho^{v}(a,t)$ be the density function in domain $v\in\mathcal{V}$ (considering at least two domains) for $a\in[0,1]$, time $t$, and traveling with speed $\nu^{v}(t)$. Given corresponding moments $(m_{0}^{v}(t),m_{1}^{v}(t),m_{2}^{v}(t))$ at time $t$, as defined in (6), we wish to determine the value of the moments at time $t+\Delta t$, accounting for the mass flow across domains. A key quantity is the length scale over which information is advected in domain $v$ from $t$ to $t+\Delta t$, which is:

For general unknown velocity functions $\nu^{v}(t)$, the integral needs to be approximated numerically. A simple way is to use $\Delta a^{v}\approx\Delta t\,\nu^{v}(t+\Delta t/2)$. However, more accurate quadratures are also possible. The mass outflux from domain $v$ is given by $\int_{1-\Delta a^{v}}^{1}\rho^{v}(a,t)\ \,\mathrm{d}a$ where

For computational accuracy and efficiency, we equivalently write the mass outflux as

Using the same shift and change of variables we additionally define the following convenient formulations for the first and second moment respectively:

The implementation of these integrals is crucial to the working of the numerical method. Note that this is the same problem addressed in §3.2.2 on computing the forward map which can be done via analytical Gaussian error function formulas or numerical quadrature. The computation of the integrals is pseudo-exact, incurring numerical errors only when numerical quadrature is used.

We now describe the changes due to the flow of mass between a sending domain $v$ and receiving domain $w$ that are part of an edge $(v,w)\in\mathcal{E}$. To obtain the correct behavior in domain $w$, when domain $v$ and domain $w$ have different time-dependent speeds, we define the ratio $\displaystyle\gamma^{v\to w}=\frac{\nu^{w}(t+\Delta t/2)}{\nu^{v}(t+\Delta t/2)}$ (which approximates $\displaystyle\frac{\nu^{w}}{\nu^{v}}$ on the interval $[t,t+\Delta t]$ with $\mathcal{O}(\Delta t^{2})$ when $\nu^{w}$ and $\nu^{v}$ are sufficiently smooth). Changes in the moments of the sending domain are