A model for the interaction of phytoplankton aggregates and the environment: approximation and parameter estimation

Abstract

We present a model for the dynamics of phytoplankton aggregates and their interaction with the environment in a well-mixed reactor or surface layer of a water column. The biophysical problem yields a nonlinear partial differential equation coupled with a general system of ordinary differential equations. We first develop a finite difference scheme for approximating the solution of this model. Then convergence results are established for the finite difference method and existence-uniqueness of solutions are obtained. A least-squares parameter estimation method is presented to illustrate the performance of the model against a real dataset, and to highlight the effects of aggregate growth on algal bloom dynamics and nutrient consumption rates. The model is well-suited to complement microcosm studies and can accommodate the effects of a wide range of environmental conditions.

Keywords:

• Phytoplankton
• aggregation
• environment
• finite difference approximation
• parameter estimation

AMS Subject Classifications:

• 35Q92
• 65M06

1. Introduction

Ample research has been undertaken to understand the role of phytoplankton in a variety of applied settings [1–8]. Phytoplankton represent an important layer in the aquatic food web of marine and freshwater ecosystems throughout the world [9]. Given the vast scope of natural phytoplankton systems (e.g. rivers, reservoirs, and oceans), it is challenging to address what-if type questions in a realistic setting [10,11]. Recent demonstrations show the potential of in situ algal stimulation by nutrient seeding in the open ocean at small scales [12,13]. However, seeding activities are challenged by logistics and other practical concerns associated with bloom stimulation in the natural environment. On the smaller scale of chemical engineering the development of commercially-viable algae-based clean energy alternatives and other green product technologies are also hindered by the difficulty of extrapolating laboratory results to practical scales of production [14–21].

Mathematical models can address some of the challenges associated with phytoplankton population dynamics, since they can include many of the variables which are too complex to investigate in practical settings. Phytoplankton populations can be modeled as a collection of irregularly-shaped particle aggregates bound together by some organic glue. The particles increase in size either by biological cell division or the physical process of coagulation when they encounter other suspended particles. Since phytoplankton growth depends on photosynthesis, cell death may occur at the core of large aggregates while cell division may persist on active portions of the particle surface. So, in general, only a portion of the aggregate experiences growth and uptake of nutrients from the water column and this phenomenon is challenging to model accurately.

Past work for the development of size-structured phytoplankton aggregation models and/or their numerical approximations includes [22–29]. Researchers have focused on establishing existence-uniqueness of solutions to the classical coagulation-fragmentation problem and an investigation of their long-term behavior [30–33]. Understanding what controls the distribution of phytoplankton and nutrient dyamics is a major challenge facing contemporary biological oceanographers, and such investigations could benefit from models capable of accomodating a wide variety of factors [34–37]. Previous studies have implemented the size-structured approach to model the behavior of flocculating bacterial and microbial aggregates in suspension systems [38–40]. However, to the best of our knowledge no size-structured model exists which includes a general array of environmental variables and numerically investigates the sensitivity of algal bloom formation and nutrient depletion in response to the interaction of coagulation, aggregate growth, and reproduction rates in a physically realistic context.

Our aim is to introduce a structured population model for phytoplankton which includes a general vector of environmental variables, growth of particle aggregates, reproduction, and coagulation in which the model parameters depend on the generalized environment. In addition to presenting the model, the objectives of this paper are (1) to present a first-order finite difference numerical approximation and to develop new innovative techniques to treat the nonlinearity associated with the coagulation term and to establish the necessary estimates for convergence of the first-order finite difference scheme, (2) to establish existence and uniqueness of the solution to the model, (3) to accurately model growth within aggregates and (4) to study the sensitivity of the model to its parameters and to utilize a least-squares approach to fit the model to a tank experiment data set.

This paper is organized as follows: in Section 2 we introduce the general nonlinear model and in Section 3 we introduce the notation and assumptions used in the forthcoming analysis. In Section 4 we setup a finite difference approximation and establish a series of technical lemmas using standard compactness arguments for scalar conservation laws in the spirit of [41] (Chapter 16) to show that the finite difference scheme converges to a weak bounded variation (BV) solution to the model. In Section 5 we show that the scheme converges to the unique BV solution. In Section 6 we develop sub-models for the vital rates and in Section 7 we setup a least-squares problem to compare the model to a dataset. This least-squares problem involves finding optimal parameter values for aggregate growth, and the fraction of live vs. dead material. In Section 8 we provide some concluding remarks. The details of the technical proofs are provided in the Appendix 1.

2. The model

Consider the well-mixed surface layer of a geographic parcel of water (commonly referred to as the water column in environmental engineering) and a population ofphytoplankton particles having number density per unit size u. The following model describes the size-structured dynamics of the system(2.1) $$\begin{array}{c}\hfill \begin{array}{c}\frac{\mathit{\partial }u}{\mathit{\partial }t}+\frac{\mathit{\partial }}{\mathit{\partial }x}\left(g(x,t,\mathit{\varphi },\overrightarrow{y})u\right)+R(x,t,\mathit{\varphi },\overrightarrow{y})u=A(u),(x,t)\in (0,x_{max}]\times (0,T],\hfill \\ \frac{\mathrm{d}{y}_{\ell }}{\mathrm{d}t}=f_{\ell }(t,{\mathit{\sigma }}_{\ell },\overrightarrow{y})-h_{\ell }(t,{\mathit{\sigma }}_{\ell },\overrightarrow{y})y_{\ell },t\in (0,T],\ell =1,\dots ,m,\hfill \\ g(0,t,\mathit{\varphi },\overrightarrow{y})u(0,t)=F(t,\overrightarrow{y})+{\int }_0^{x_{\text{max}}}\mathit{\beta }(x,t,\mathit{\varphi },\overrightarrow{y})u(x,t)\mathrm{d}x,t\in (0,T],\hfill \\ u(x,0)=u^0(x),x\in [0,x_{max}],y_{\ell }(0)=y_{\ell ,0}\ge 0,\ell =1,\dots ,m,\hfill \end{array}\end{array}$$(2.1)

where(2.2) $$\begin{array}{c}\hfill \begin{array}{c}A(u)=\frac{1}{2}{\int }_0^x\mathit{\eta }(x-s,s,t,\mathit{\varphi },\overrightarrow{y})u(x-s,t)u(s,t)\mathrm{d}s\hfill \\ -{\int }_0^{x_{max}}\mathit{\eta }(x,s,t,\mathit{\varphi },\overrightarrow{y})u(x,t)u(s,t)\mathrm{d}s,\hfill \end{array}\end{array}$$(2.2)

and(2.3) $$\begin{array}{c}\hfill \mathit{\varphi }(t)={\int }_0^{x_{max}}\mathit{\omega }(x)u(x,t)\mathrm{d}x,{\mathit{\sigma }}_{\ell }(t)={\int }_0^{x_{max}}{\mathit{\gamma }}_{\ell }(x)u(x,t)\mathrm{d}x,t\in (0,T].\end{array}$$(2.3)

The term A(u) is used to model the particle coalescence process [42,43]. The function $\mathit{\eta }(x,s,t,·)$ represents the rate at which particles of size x coagulate with particles of size s to yield particles of size $x+s$ at a given moment in time t. The parameter g is the size-specific growth rate of a particle aggregate, and R is the individual mortality, or rate of disappearance of aggregates of size x from the surface layer. The function F measures the addition of particles of smallest size from an outside source and $\mathit{\beta }$ is the individual reproduction rate of newborn cells. The term $\overrightarrow{y}$ is a vector containing the environmental variables (e.g. temperature, water chemistry, grazing, light intensity, nutrients, etc.) all of which interact with each other and the particle population. The components of the vector $\overrightarrow{y}$ are $y_{\ell }$ and the functions $f_{\ell }$ and $h_{\ell }$ are nonlinear environmental source and mortality terms, respectively which also depend on the particle population. The evolution of the environment is governed by an m-dimensional system of nonlinear ordinary differential equations. The variable $\mathit{\varphi }$ represents some weighted average of the population with weight function $\mathit{\omega }$. The variables ${\mathit{\sigma }}_{\ell }$ and ${\mathit{\gamma }}_{\ell }$ are weighted averages of the population and weight functions, respectively which interact with the environmental variables $\ell =1,\dots ,m$. In Section 6 and Table 1 we provide an example featuring specific sub-models, descriptions, and a summary of units for model variables and parameters in the SIGMA tank phytoplankton experiment [23].

The model (2.1(2.1) $$\begin{array}{c}\hfill \begin{array}{c}\frac{\mathit{\partial }u}{\mathit{\partial }t}+\frac{\mathit{\partial }}{\mathit{\partial }x}\left(g(x,t,\mathit{\varphi },\overrightarrow{y})u\right)+R(x,t,\mathit{\varphi },\overrightarrow{y})u=A(u),(x,t)\in (0,x_{max}]\times (0,T],\hfill \\ \frac{\mathrm{d}{y}_{\ell }}{\mathrm{d}t}=f_{\ell }(t,{\mathit{\sigma }}_{\ell },\overrightarrow{y})-h_{\ell }(t,{\mathit{\sigma }}_{\ell },\overrightarrow{y})y_{\ell },t\in (0,T],\ell =1,\dots ,m,\hfill \\ g(0,t,\mathit{\varphi },\overrightarrow{y})u(0,t)=F(t,\overrightarrow{y})+{\int }_0^{x_{\text{max}}}\mathit{\beta }(x,t,\mathit{\varphi },\overrightarrow{y})u(x,t)\mathrm{d}x,t\in (0,T],\hfill \\ u(x,0)=u^0(x),x\in [0,x_{max}],y_{\ell }(0)=y_{\ell ,0}\ge 0,\ell =1,\dots ,m,\hfill \end{array}\end{array}$$(2.1) ) studied in this paper utilizes concepts developed by Smoluchowski to investigate particle dynamics in the presence of coagulation [44]. It is a combination of the continuous version of the coagulation model used to study aggregation in well-mixed particle systems [45] and the classical Sinko-Streifer model of size-structured populations subject to growth, mortality, and reproduction [46]. This combination of factors takes into account the biological and physical variables which determine the fate of an individual phytoplankton aggregate. The model also generalizes the concept of environment by including an arbitrary number of variables ($\ell =1,\dots ,m$) which are coupled with the vital rates of the particle population. This yields a quasilinear hyperbolic partial differential equation coupled with a system of nonlinear ordinary differential equations containing nonlocal terms. The terms $\mathit{\varphi }$ and ${\mathit{\sigma }}_{\ell }$ are nonlocal in the sense that they depend on the entire particle population rather than the (local) population density at a given point. Theoretically, this model belongs to the class of 1-dimensional scalar conservation laws in which the conserved variable is the population density of cell aggregates. It is well known that solutions to the 1-dimensional particle aggregation problem can exhibit gelation (i.e. emergence of a particle having infinite size) in finite time for certain specific forms of the coagulation function. Gelation is not known to occur in real phytoplankton systems and we will establish conditions which guarantee the absence of gel-forming solutions to the model (2.1(2.1) $$\begin{array}{c}\hfill \begin{array}{c}\frac{\mathit{\partial }u}{\mathit{\partial }t}+\frac{\mathit{\partial }}{\mathit{\partial }x}\left(g(x,t,\mathit{\varphi },\overrightarrow{y})u\right)+R(x,t,\mathit{\varphi },\overrightarrow{y})u=A(u),(x,t)\in (0,x_{max}]\times (0,T],\hfill \\ \frac{\mathrm{d}{y}_{\ell }}{\mathrm{d}t}=f_{\ell }(t,{\mathit{\sigma }}_{\ell },\overrightarrow{y})-h_{\ell }(t,{\mathit{\sigma }}_{\ell },\overrightarrow{y})y_{\ell },t\in (0,T],\ell =1,\dots ,m,\hfill \\ g(0,t,\mathit{\varphi },\overrightarrow{y})u(0,t)=F(t,\overrightarrow{y})+{\int }_0^{x_{\text{max}}}\mathit{\beta }(x,t,\mathit{\varphi },\overrightarrow{y})u(x,t)\mathrm{d}x,t\in (0,T],\hfill \\ u(x,0)=u^0(x),x\in [0,x_{max}],y_{\ell }(0)=y_{\ell ,0}\ge 0,\ell =1,\dots ,m,\hfill \end{array}\end{array}$$(2.1) ). Also, in most applications one should expect that the initial data $u_{\text{0}}$ will contain gaps and noise (i.e. discontinuities), and these considerations motivate us to pursue weak solutions in the space of functions having bounded total variation, BV.

3. Assumptions and weak solution definition

Let $\mathrm{\Omega }\subset {\mathbb{R}}^d$, $d\in \mathbb{N}$. We denote by $C(\mathrm{\Omega })$ the space of continuous functions on $\mathrm{\Omega }$, by $C^1(\mathrm{\Omega })$ the space of continuously differentiable functions on $\mathrm{\Omega }$, and by $BV(\mathrm{\Omega })$ the space of functions with bounded total variation on $\mathrm{\Omega }$. Recall that a function f over $\mathrm{\Omega }$ is said to be locally Lipschitz continuous if for every $x\in \mathrm{\Omega }$ there exists a neighborhood U of x such that f restricted to U is Lipschitz continuous. We denote by ${\text{Lip}}_{loc}(\mathrm{\Omega })$ the space of locally Lipschitz continuous functions on $\mathrm{\Omega }$ and by $\text{Lip}(\mathrm{\Omega })$ the space of Lipschitz continuous functions on $\mathrm{\Omega }$. Throughout the paper we will use the shorthand notation $L_f$ to represent the Lipschitz constant for a Lipchitz continuous function f. For example, $L_{{\mathit{\gamma }}_{\ell }}$ will denote the Lipschitz constant for the function ${\mathit{\gamma }}_{\ell }$ and we omit the index $\ell$ when denoting the maximal Lipschitz constant for a given function over all environmental factors, i.e. $L_{\mathit{\gamma }}\equiv max\{L_{{\mathit{\gamma }}_1},\dots ,L_{{\mathit{\gamma }}_m}\}$.

Let $c_0,c_1,$ and ${\mathit{\omega }}_1$ be sufficiently large positive constants. Denote by ${\mathbb{R}}_{+}^d$ the positive cone of ${\mathbb{R}}^d$, $d\in \mathbb{N}$ and let $\mathbb{D}=[0,x_{max}]\times [0,T]\times {\mathbb{R}}_{+}^{m+1}$, ${\mathbb{D}}_x\equiv [0,x_{max}]\times [0,T]$, ${\mathbb{D}}_t=[0,T]\times {\mathbb{R}}_{+}^{m+2}$, and $|\overrightarrow{y}{|}_1={\mathrm{\Sigma }}_{\ell =1}^m|y_{\ell }|$. Throughout the discussion we impose the following conditions on our model parameters in (2.1(2.1) $$\begin{array}{c}\hfill \begin{array}{c}\frac{\mathit{\partial }u}{\mathit{\partial }t}+\frac{\mathit{\partial }}{\mathit{\partial }x}\left(g(x,t,\mathit{\varphi },\overrightarrow{y})u\right)+R(x,t,\mathit{\varphi },\overrightarrow{y})u=A(u),(x,t)\in (0,x_{max}]\times (0,T],\hfill \\ \frac{\mathrm{d}{y}_{\ell }}{\mathrm{d}t}=f_{\ell }(t,{\mathit{\sigma }}_{\ell },\overrightarrow{y})-h_{\ell }(t,{\mathit{\sigma }}_{\ell },\overrightarrow{y})y_{\ell },t\in (0,T],\ell =1,\dots ,m,\hfill \\ g(0,t,\mathit{\varphi },\overrightarrow{y})u(0,t)=F(t,\overrightarrow{y})+{\int }_0^{x_{\text{max}}}\mathit{\beta }(x,t,\mathit{\varphi },\overrightarrow{y})u(x,t)\mathrm{d}x,t\in (0,T],\hfill \\ u(x,0)=u^0(x),x\in [0,x_{max}],y_{\ell }(0)=y_{\ell ,0}\ge 0,\ell =1,\dots ,m,\hfill \end{array}\end{array}$$(2.1) ), for each $\ell =1,2,\dots ,m$.

(B1)	The function $g\in C^1(\mathbb{D})$ is nonnegative, and ${(g)}_x\in {\text{Lip}}_{\text{loc}}(\mathbb{D})$. Moreover, there exists a constant $\mathit{\zeta }>0$ such that $g(0,t,\mathit{\varphi },\overrightarrow{y})\ge \mathit{\zeta }>0$, and $g(x_{max},t,\mathit{\varphi },\overrightarrow{y})=0$.
(B2)	$R\in {\text{Lip}}_{\text{loc}}(\mathbb{D})$ is nonnegative.
(B3)	$f_{\ell }\in {\text{Lip}}_{\text{loc}}({\mathbb{D}}_t)$ is nonnegative, and for each $t\ge 0$ we have $f_{\ell }\le c_1+c_1|{\mathit{\sigma }}_{\ell }|+c_1{|\overrightarrow{y}|}_1$.
(B4)	$h_{\ell }\in {\text{Lip}}_{\text{loc}}({\mathbb{D}}_t)$ is nonnegative, and for each $t\ge 0$ $h_{\ell }\le c_1+c_1|{\mathit{\sigma }}_{\ell }|+c_1{|\overrightarrow{y}|}_1$.
(B5)	For any sufficiently small $\mathit{\delta }>0$ $$\begin{array}{c}\hfill \begin{array}{c}\underset{(x,t,\mathit{\varphi },\overrightarrow{y})\in \mathbb{D}}{sup}\{\left|\frac{g(x+\mathit{\delta },t,\mathit{\varphi },\overrightarrow{y})-g(x,t,\mathit{\varphi },\overrightarrow{y})}{\mathit{\delta }}\right|+R(x,t,\mathit{\varphi },\overrightarrow{y})\}\le {\mathit{\omega }}_1.\hfill \end{array}\end{array}$$
(B6)	$\mathit{\beta }\in {\text{Lip}}_{\text{loc}}(\mathbb{D})$ is nonnegative. Furthermore, $$\underset{(x,t,\mathit{\varphi },\overrightarrow{y})\in \mathbb{D}}{sup}\mathit{\beta }(x,t,\mathit{\varphi },\overrightarrow{y})\le c_0.$$
(B7)	$F\in {\text{Lip}}_{\text{loc}}({\mathbb{R}}^{m+1})$ is nonnegative. For each $t\ge 0$ we have $F(t,\overrightarrow{y})\le c_1+c_1{|\overrightarrow{y}|}_1$.
(B8)	The nonnegative weight functions $\mathit{\omega }\in \text{Lip}(\mathbb{R}),{\mathit{\gamma }}_{\ell }\in \text{Lip}(\mathbb{R})$ for all $\ell =1,\dots ,m$ with ${max}_{\ell =1,\dots ,m}\{\mathit{\omega },{\mathit{\gamma }}_{\ell }\}\le c_0$.
(B9)	$u^0$ and $y_{\ell }^0$ are nonnegative, and $u^0\in BV([0,x_{max}])$.
(B10)	$\mathit{\eta }\in {\text{Lip}}_{\text{loc}}({\mathbb{R}}^{4+m})$, is symmetric with respect to size, i.e. $\mathit{\eta }(x,s,t,\mathit{\varphi },\overrightarrow{y})=\mathit{\eta }(s,x,t,\mathit{\varphi },\overrightarrow{y})$, nonnegative, and uniformly bounded (i.e. $\mathit{\eta }(x,s,t,\mathit{\varphi },\overrightarrow{y})\le c_0$). Furthermore, $$\begin{array}{c}\hfill \begin{array}{c}\mathit{\eta }(x,s,t,\mathit{\varphi },\overrightarrow{y})\ge 0,x+s\le x_{max}\hfill \\ \mathit{\eta }(x,s,t,\mathit{\varphi },\overrightarrow{y})=0,x+s>x_{max}.\hfill \end{array}\end{array}$$

Next, we give the following definition of a weak solution to problem (2.1(2.1) $$\begin{array}{c}\hfill \begin{array}{c}\frac{\mathit{\partial }u}{\mathit{\partial }t}+\frac{\mathit{\partial }}{\mathit{\partial }x}\left(g(x,t,\mathit{\varphi },\overrightarrow{y})u\right)+R(x,t,\mathit{\varphi },\overrightarrow{y})u=A(u),(x,t)\in (0,x_{max}]\times (0,T],\hfill \\ \frac{\mathrm{d}{y}_{\ell }}{\mathrm{d}t}=f_{\ell }(t,{\mathit{\sigma }}_{\ell },\overrightarrow{y})-h_{\ell }(t,{\mathit{\sigma }}_{\ell },\overrightarrow{y})y_{\ell },t\in (0,T],\ell =1,\dots ,m,\hfill \\ g(0,t,\mathit{\varphi },\overrightarrow{y})u(0,t)=F(t,\overrightarrow{y})+{\int }_0^{x_{\text{max}}}\mathit{\beta }(x,t,\mathit{\varphi },\overrightarrow{y})u(x,t)\mathrm{d}x,t\in (0,T],\hfill \\ u(x,0)=u^0(x),x\in [0,x_{max}],y_{\ell }(0)=y_{\ell ,0}\ge 0,\ell =1,\dots ,m,\hfill \end{array}\end{array}$$(2.1) ).

Definition 3.1:

An $(m+1)$-tuple of functions $(u,\overrightarrow{y})\in L_{\infty }({\mathbb{D}}_x)\times {\mathrm{\Pi }}_{\ell =1}^{m}C[0,T]$ is called a weak solution to problem (2.1(2.1) $$\begin{array}{c}\hfill \begin{array}{c}\frac{\mathit{\partial }u}{\mathit{\partial }t}+\frac{\mathit{\partial }}{\mathit{\partial }x}\left(g(x,t,\mathit{\varphi },\overrightarrow{y})u\right)+R(x,t,\mathit{\varphi },\overrightarrow{y})u=A(u),(x,t)\in (0,x_{max}]\times (0,T],\hfill \\ \frac{\mathrm{d}{y}_{\ell }}{\mathrm{d}t}=f_{\ell }(t,{\mathit{\sigma }}_{\ell },\overrightarrow{y})-h_{\ell }(t,{\mathit{\sigma }}_{\ell },\overrightarrow{y})y_{\ell },t\in (0,T],\ell =1,\dots ,m,\hfill \\ g(0,t,\mathit{\varphi },\overrightarrow{y})u(0,t)=F(t,\overrightarrow{y})+{\int }_0^{x_{\text{max}}}\mathit{\beta }(x,t,\mathit{\varphi },\overrightarrow{y})u(x,t)\mathrm{d}x,t\in (0,T],\hfill \\ u(x,0)=u^0(x),x\in [0,x_{max}],y_{\ell }(0)=y_{\ell ,0}\ge 0,\ell =1,\dots ,m,\hfill \end{array}\end{array}$$(2.1) ) if for $t\in [0,T]$, $u(·,t)\in BV([0,x_{max}])$ and the following holds for every test function $\mathit{\rho }\in C^1((0,x_{max})\times (0,T))$ for all $\ell =1,2,\dots ,m$:$$\begin{array}{c}\hfill \begin{array}{cc}& {\int }_0^{x_{max}}(\mathit{\rho }(x,t)u(x,t)-\mathit{\rho }(x,0)u(x,0))\mathrm{d}x-{\int }_0^t{\int }_0^{x_{max}}\frac{\mathit{\partial }\mathit{\rho }}{\mathit{\partial }t}u\mathrm{d}x\mathrm{d}\mathit{\tau }\hfill \\ \hfill & -{\int }_0^t\mathit{\rho }(0,\mathit{\tau })(F+{\int }_0^{x_{max}}\mathit{\beta }u\mathrm{d}x)\mathrm{d}\mathit{\tau }-{\int }_0^t{\int }_0^{x_{max}}\frac{\mathit{\partial }\mathit{\rho }}{\mathit{\partial }x}gu\mathrm{d}x\mathrm{d}\mathit{\tau }+{\int }_0^t{\int }_0^{x_{max}}\mathit{\rho }Ru\mathrm{d}x\mathrm{d}\mathit{\tau }\hfill \\ \hfill & ={\int }_0^t{\int }_0^{x_{max}}{\int }_0^x\mathit{\rho }(x,\mathit{\tau })A(u)\mathrm{d}x\mathrm{d}\mathit{\tau },\hfill \end{array}\end{array}$$

and$$\begin{array}{c}\hfill \begin{array}{c}{y}_{\ell }(t)=y_{\ell }(0)+{\int }_0^t\left[f_{\ell }(\mathit{\tau },{\mathit{\sigma }}_{\ell }(\mathit{\tau }),y_{\ell }(\mathit{\tau }))-h_{\ell }(\mathit{\tau },{\mathit{\sigma }}_{\ell }(\mathit{\tau }),y_{\ell }(\mathit{\tau }))y_{\ell }(\mathit{\tau })\right]\mathrm{d}\mathit{\tau }.\hfill \end{array}\end{array}$$

4. Difference approximations

Suppose that the intervals $[0,x_{max}]$ and [0, T] are divided into N and L subintervals respectively. The following notation will be used throughout this section: $\mathrm{\Delta }x=x_{max}/N$, and $\mathrm{\Delta }t=T/L$ denote the spatial and time mesh sizes respectively. The mesh points are given by $x_j=j\mathrm{\Delta }x,j=0,1,2,\dots ,N$, and $t^k=k\mathrm{\Delta }t,k=0,1,2,\dots ,L$. For convenience of notation we apply equal mesh sizes in the x dimension and a uniform time step. The proof of convergence can be established in a similar manner using a more general non-uniform mesh. We denote by $u_j^k$, $y_{\ell }^k$, ${\mathit{\varphi }}^k$, and ${\mathit{\sigma }}_{\ell }^k$, the finite difference approximations of $u(x_j,t^k)$, $y_{\ell }(t^k)$, $\mathit{\varphi }(t^k)$, and ${\mathit{\sigma }}_{\ell }(t^k)$, respectively. The values of the model parameters evaluated at the mesh points are given below,$$\begin{array}{c}\hfill \begin{array}{c}{\mathit{\beta }}_j^k=\mathit{\beta }(x_j,t^k,{\mathit{\varphi }}^k,{\overrightarrow{y}}^k),g_j^k=g(x_j,t^k,{\mathit{\varphi }}^k,{\overrightarrow{y}}^k),R_j^k=R(x_j,t^k,{\mathit{\varphi }}^k,{\overrightarrow{y}}^k),\hfill \\ F^k=F(t^k,{\overrightarrow{y}}^k),f_{\ell }^k=f_{\ell }(t^k,{\mathit{\sigma }}_{\ell }^k,{\overrightarrow{y}}^k),h_{\ell }^k=h_{\ell }(t^k,{\mathit{\sigma }}_{\ell }^k,{\overrightarrow{y}}^k),{\mathit{\eta }}_{j,s}^k={\mathit{\eta }}_p(x_j,x_s,t^k,{\mathit{\varphi }}^k,{\overrightarrow{y}}^k),\hfill \\ {\mathit{\omega }}_j=\mathit{\omega }(x_j),{\mathit{\gamma }}_{\ell ,j}={\mathit{\gamma }}_{\ell }(x_j),\hfill \end{array}\end{array}$$

defined for all $j,s=0,\dots ,N$, $k=0,\dots ,L$, and $\ell =1,\dots ,m$.

We define the ${\ell }^1$ and ${\ell }^{\infty }$ norms and the TV (total variation) seminorm of the grid functions $u^k$ and ${\overrightarrow{y}}^k$ by$$\begin{array}{c}\hfill \begin{array}{c}\Vert u^k{\Vert }_1\equiv \sum _{j=1}^N\left|u_j^k\right|\mathrm{\Delta }x,\Vert u^k{\Vert }_{\infty }\equiv {max}_{0\le j\le N}\left|u_j^k\right|,|{\overrightarrow{y}}^k{|}_{\infty }\equiv {max}_{1\le \ell \le m}\{|y_{\ell }^k|\}\hfill \\ |{\overrightarrow{y}}^k{|}_1\equiv \sum _{\ell =1}^m|y_{\ell }^k|,TV(u^{k+1})\equiv {\sum }_{j=0}^{N-1}\left|u_{j+1}^{k+1}-u_j^{k+1}\right|\hfill \end{array}\end{array}$$

for $k=0,1,\dots ,L$. We also adopt the following nonstandard first-order finite difference scheme for $j=1,\dots ,N$:(4.1) $$\begin{array}{cc}& \frac{u_j^{k+1}-u_j^k}{\mathrm{\Delta }t}+\frac{g_j^k{u}_j^{k+1}-g_{j-1}^k{u}_{j-1}^{k+1}}{\mathrm{\Delta }x}+u_j^{k+1}\sum _{m=1}^N{\mathit{\eta }}_{m,j}^k{u}_m^k\mathrm{\Delta }x+R_j^k{u}_j^{k+1}\hfill \\ \hfill & =\frac{1}{2}\sum _{s=1}^{j-1}{\mathit{\eta }}_{j-s,s}^k{u}_{j-s}^k{u}_s^{k+1}\mathrm{\Delta }x,\hfill \\ \hfill & g_0^k{u}_0^{k+1}=F^k+\sum _{j=1}^N{\mathit{\beta }}_j^k{u}_j^k\mathrm{\Delta }x,u_j^0=u_0(x_j),\hfill \\ \hfill & \frac{y_{\ell }^{k+1}-y_{\ell }^k}{\mathrm{\Delta }t}=f_{\ell }^k-h_{\ell }^k{y}_{\ell }^{k+1},y_{\ell }^0=y_{\ell ,0},\ell =1,\dots ,m.\hfill \end{array}$$(4.1)

First observe that the first equation in (4.1(4.1) $$\begin{array}{cc}& \frac{u_j^{k+1}-u_j^k}{\mathrm{\Delta }t}+\frac{g_j^k{u}_j^{k+1}-g_{j-1}^k{u}_{j-1}^{k+1}}{\mathrm{\Delta }x}+u_j^{k+1}\sum _{m=1}^N{\mathit{\eta }}_{m,j}^k{u}_m^k\mathrm{\Delta }x+R_j^k{u}_j^{k+1}\hfill \\ \hfill & =\frac{1}{2}\sum _{s=1}^{j-1}{\mathit{\eta }}_{j-s,s}^k{u}_{j-s}^k{u}_s^{k+1}\mathrm{\Delta }x,\hfill \\ \hfill & g_0^k{u}_0^{k+1}=F^k+\sum _{j=1}^N{\mathit{\beta }}_j^k{u}_j^k\mathrm{\Delta }x,u_j^0=u_0(x_j),\hfill \\ \hfill & \frac{y_{\ell }^{k+1}-y_{\ell }^k}{\mathrm{\Delta }t}=f_{\ell }^k-h_{\ell }^k{y}_{\ell }^{k+1},y_{\ell }^0=y_{\ell ,0},\ell =1,\dots ,m.\hfill \end{array}$$(4.1) ) has a unique nonnegative algebraic solution and this result satisfies a biological requirement. This can be seen by multiplying the first row of (4.1(4.1) $$\begin{array}{cc}& \frac{u_j^{k+1}-u_j^k}{\mathrm{\Delta }t}+\frac{g_j^k{u}_j^{k+1}-g_{j-1}^k{u}_{j-1}^{k+1}}{\mathrm{\Delta }x}+u_j^{k+1}\sum _{m=1}^N{\mathit{\eta }}_{m,j}^k{u}_m^k\mathrm{\Delta }x+R_j^k{u}_j^{k+1}\hfill \\ \hfill & =\frac{1}{2}\sum _{s=1}^{j-1}{\mathit{\eta }}_{j-s,s}^k{u}_{j-s}^k{u}_s^{k+1}\mathrm{\Delta }x,\hfill \\ \hfill & g_0^k{u}_0^{k+1}=F^k+\sum _{j=1}^N{\mathit{\beta }}_j^k{u}_j^k\mathrm{\Delta }x,u_j^0=u_0(x_j),\hfill \\ \hfill & \frac{y_{\ell }^{k+1}-y_{\ell }^k}{\mathrm{\Delta }t}=f_{\ell }^k-h_{\ell }^k{y}_{\ell }^{k+1},y_{\ell }^0=y_{\ell ,0},\ell =1,\dots ,m.\hfill \end{array}$$(4.1) ) by $\mathrm{\Delta }t$ and rewriting equation to obtain(4.2) $$\begin{array}{c}\hfill \begin{array}{c}{u}_j^{k+1}=\frac{u_j^k+\frac{\mathrm{\Delta }t}{\mathrm{\Delta }x}{g}_{j-1}^k{u}_{j-1}^{k+1}+\frac{\mathrm{\Delta }t}{2}{\sum }_{s=1}^{j-1}{\mathit{\eta }}_{j-s,s}^k{u}_{j-s}^k{u}_s^{k+1}\mathrm{\Delta }x}{1+\frac{\mathrm{\Delta }t}{\mathrm{\Delta }x}{g}_j^k+\mathrm{\Delta }t{\sum }_{m=1}^N{\mathit{\eta }}_{m,j}^k{u}_m^k\mathrm{\Delta }x+\mathrm{\Delta }t{R}_j^k}.\hfill \end{array}\end{array}$$(4.2)

All terms on the right hand side of (4.2(4.2) $$\begin{array}{c}\hfill \begin{array}{c}{u}_j^{k+1}=\frac{u_j^k+\frac{\mathrm{\Delta }t}{\mathrm{\Delta }x}{g}_{j-1}^k{u}_{j-1}^{k+1}+\frac{\mathrm{\Delta }t}{2}{\sum }_{s=1}^{j-1}{\mathit{\eta }}_{j-s,s}^k{u}_{j-s}^k{u}_s^{k+1}\mathrm{\Delta }x}{1+\frac{\mathrm{\Delta }t}{\mathrm{\Delta }x}{g}_j^k+\mathrm{\Delta }t{\sum }_{m=1}^N{\mathit{\eta }}_{m,j}^k{u}_m^k\mathrm{\Delta }x+\mathrm{\Delta }t{R}_j^k}.\hfill \end{array}\end{array}$$(4.2) ) are nonnegative by the nonnegativity assumption (B9) on the initial data $u^0$. Furthermore, the solution to (4.2(4.2) $$\begin{array}{c}\hfill \begin{array}{c}{u}_j^{k+1}=\frac{u_j^k+\frac{\mathrm{\Delta }t}{\mathrm{\Delta }x}{g}_{j-1}^k{u}_{j-1}^{k+1}+\frac{\mathrm{\Delta }t}{2}{\sum }_{s=1}^{j-1}{\mathit{\eta }}_{j-s,s}^k{u}_{j-s}^k{u}_s^{k+1}\mathrm{\Delta }x}{1+\frac{\mathrm{\Delta }t}{\mathrm{\Delta }x}{g}_j^k+\mathrm{\Delta }t{\sum }_{m=1}^N{\mathit{\eta }}_{m,j}^k{u}_m^k\mathrm{\Delta }x+\mathrm{\Delta }t{R}_j^k}.\hfill \end{array}\end{array}$$(4.2) ) is obtained by a left-to-right substitution process. From assumptions (B3), (B4), and (B9), the algebraic solution to the system of ordinary differential equations given by the third equation in (4.1(4.1) $$\begin{array}{cc}& \frac{u_j^{k+1}-u_j^k}{\mathrm{\Delta }t}+\frac{g_j^k{u}_j^{k+1}-g_{j-1}^k{u}_{j-1}^{k+1}}{\mathrm{\Delta }x}+u_j^{k+1}\sum _{m=1}^N{\mathit{\eta }}_{m,j}^k{u}_m^k\mathrm{\Delta }x+R_j^k{u}_j^{k+1}\hfill \\ \hfill & =\frac{1}{2}\sum _{s=1}^{j-1}{\mathit{\eta }}_{j-s,s}^k{u}_{j-s}^k{u}_s^{k+1}\mathrm{\Delta }x,\hfill \\ \hfill & g_0^k{u}_0^{k+1}=F^k+\sum _{j=1}^N{\mathit{\beta }}_j^k{u}_j^k\mathrm{\Delta }x,u_j^0=u_0(x_j),\hfill \\ \hfill & \frac{y_{\ell }^{k+1}-y_{\ell }^k}{\mathrm{\Delta }t}=f_{\ell }^k-h_{\ell }^k{y}_{\ell }^{k+1},y_{\ell }^0=y_{\ell ,0},\ell =1,\dots ,m.\hfill \end{array}$$(4.1) ) is also nonnegative since this equation can be rewritten as(4.3) $$\begin{array}{c}\hfill \begin{array}{c}{y}_{\ell }^{k+1}=\frac{y_{\ell }^k+f_{\ell }^k\mathrm{\Delta }t}{1+h_{\ell }^k\mathrm{\Delta }t}\ge 0.\hfill \end{array}\end{array}$$(4.3)

4.1. Estimates for the finite difference approximations

The particle coagulation problem (2.1(2.1) $$\begin{array}{c}\hfill \begin{array}{c}\frac{\mathit{\partial }u}{\mathit{\partial }t}+\frac{\mathit{\partial }}{\mathit{\partial }x}\left(g(x,t,\mathit{\varphi },\overrightarrow{y})u\right)+R(x,t,\mathit{\varphi },\overrightarrow{y})u=A(u),(x,t)\in (0,x_{max}]\times (0,T],\hfill \\ \frac{\mathrm{d}{y}_{\ell }}{\mathrm{d}t}=f_{\ell }(t,{\mathit{\sigma }}_{\ell },\overrightarrow{y})-h_{\ell }(t,{\mathit{\sigma }}_{\ell },\overrightarrow{y})y_{\ell },t\in (0,T],\ell =1,\dots ,m,\hfill \\ g(0,t,\mathit{\varphi },\overrightarrow{y})u(0,t)=F(t,\overrightarrow{y})+{\int }_0^{x_{\text{max}}}\mathit{\beta }(x,t,\mathit{\varphi },\overrightarrow{y})u(x,t)\mathrm{d}x,t\in (0,T],\hfill \\ u(x,0)=u^0(x),x\in [0,x_{max}],y_{\ell }(0)=y_{\ell ,0}\ge 0,\ell =1,\dots ,m,\hfill \end{array}\end{array}$$(2.1) ) presents mathematical difficulties which are not present in the classical nonlinear Sinko-Streifer or McKendrick von-Forester type structured models. In the Sinko-Streifer size-structured model, individuals may advance in size only through the (continuous) process of biological growth. However, in a particle system subject to aggregation, coagulation causes individual particles to spontaneously increase in size via coalescence and this generally requires an unbounded region to ensure that no particles escape the domain. In numerical simulations and the problem considered here, we truncate the problem to a bounded domain $[0,x_{max}]\times [0,T]$ instead of the classical unbounded domain $[0,\infty )\times [0,T]$ used to study general particle aggregation problems. Since the coagulation process manufactures particles of increasing size, in order to prevent the unphysical emergence of a particle of infinite size (i.e. super-particle) in the context of phytoplankton aggregation, we require that $\mathit{\eta }(x,s)=0$ if $x+s>x_{max}$. The next challenge is to ensure the conservation of particles. In the continuous case, a straightforward calculation establishes the following result,$$\begin{array}{c}\hfill \begin{array}{c}{\int }_0^{x_{max}}{\int }_0^x\mathit{\eta }(x-s,s)u(x-s,t)u(s,t)\mathrm{d}s\mathrm{d}x={\int }_0^{x_{max}}u(x,t){\int }_0^{x_{max}}\mathit{\eta }(x,s)u(s,t)\mathrm{d}s\mathrm{d}x,\hfill \end{array}\end{array}$$

and hence in the absence of reproduction, mortality, and immigration (i.e. $\mathit{\beta }=0,R=0,F=0$), the total number of particles in the system satisfies(4.4) $$\begin{array}{c}\hfill \begin{array}{c}\frac{\mathrm{d}}{\mathrm{d}t}{\int }_0^{x_{max}}u(x,t)\mathrm{d}x={\int }_0^{x_{max}}A(u)\mathrm{d}x=-\frac{1}{2}{\int }_0^{x_{max}}u(x,t){\int }_0^{x_{max}}\mathit{\eta }(x,s)u(s,t)\mathrm{d}s\mathrm{d}x,\hfill \end{array}\end{array}$$(4.4)

which is an expression which measures the rate of change of the number of particles in the system. Physically this result states that in the absence of reproduction, the coagulation process reduces the total number of particles over time. This is an important result and the following lemma shows that it also holds in the discrete case as a result of the implicit differencing applied to the quadratic term on the right hand side in the first row of 4.1(4.1) $$\begin{array}{cc}& \frac{u_j^{k+1}-u_j^k}{\mathrm{\Delta }t}+\frac{g_j^k{u}_j^{k+1}-g_{j-1}^k{u}_{j-1}^{k+1}}{\mathrm{\Delta }x}+u_j^{k+1}\sum _{m=1}^N{\mathit{\eta }}_{m,j}^k{u}_m^k\mathrm{\Delta }x+R_j^k{u}_j^{k+1}\hfill \\ \hfill & =\frac{1}{2}\sum _{s=1}^{j-1}{\mathit{\eta }}_{j-s,s}^k{u}_{j-s}^k{u}_s^{k+1}\mathrm{\Delta }x,\hfill \\ \hfill & g_0^k{u}_0^{k+1}=F^k+\sum _{j=1}^N{\mathit{\beta }}_j^k{u}_j^k\mathrm{\Delta }x,u_j^0=u_0(x_j),\hfill \\ \hfill & \frac{y_{\ell }^{k+1}-y_{\ell }^k}{\mathrm{\Delta }t}=f_{\ell }^k-h_{\ell }^k{y}_{\ell }^{k+1},y_{\ell }^0=y_{\ell ,0},\ell =1,\dots ,m.\hfill \end{array}$$(4.1) . The result is established below.

Lemma 4.1:

Assume that $u_j^k$ is nonnegative. Then we have the following:(4.5) $$\begin{array}{c}\hfill \begin{array}{cc}& \frac{1}{2}{\sum }_{j=1}^N{\sum }_{s=1}^{j-1}{\mathit{\eta }}_{j-s,s}^k{u}_{j-s}^k{u}_s^{k+1}\mathrm{\Delta }x\mathrm{\Delta }x-{\sum }_{j=1}^N{u}_j^{k+1}{\sum }_{m=1}^N{\mathit{\eta }}_{m,j}^k{u}_m^k\mathrm{\Delta }x\mathrm{\Delta }x\hfill \\ \hfill & =-\frac{1}{2}{\sum }_{j=1}^N{u}_j^{k+1}{\sum }_{m=1}^N{\mathit{\eta }}_{m,j}^k{u}_m^k\mathrm{\Delta }x\mathrm{\Delta }x,\hfill \end{array}\end{array}$$(4.5)

for any $N\ge 1$.

Proof See Appendix 1.

Next we will establish an ${\ell }^1$ estimate on the finite difference scheme. One notable consequence of the ${\ell }^1$ estimate is that it ensures the absence of gelation (e.g. [47]). This is a biologically realisitic requirement for numerical solutions of the phytoplankton problem.

Lemma 4.2:

The scheme is bounded in ${\ell }^1$, i.e. there exists a constant $B_2$ such that for all $k=0,\dots ,L-1$,$$\begin{array}{c}\hfill \begin{array}{c}\Vert u^{k+1}{\Vert }_1+{|{\overrightarrow{y}}^{k+1}|}_1\le B_2\hfill \end{array}\end{array}$$

Proof See Appendix 1.

Note that from the inequality established in Lemma 4.2 we have,$$\begin{array}{c}\hfill \begin{array}{c}{\mathit{\varphi }}^{k+1}\le c_0{B}_2=B_{1,1}\hfill \end{array}\end{array}$$

and similarly ${\mathit{\sigma }}_{\ell }\le B_{1,1}$. Thus, using (B1)-(B10), we obtain a uniform upper bound on the parameters $g,{(g)}_x,R,F$ over the compact set $[0,x_{max}]\times [0,T]\times [0,B_{1,1}]\times {\mathrm{\Pi }}_{s=1}^m[0,B_2]$. For convenience we denote this bound by $B_{1,1}$ as well. Next we establish an ${\ell }^{\infty }$ stability bound on the numerical approximations.

Lemma 4.3:

There exists a constant $B_3$ such that if $\mathrm{\Delta }t$ satisfies $1-{\mathit{\omega }}_1\mathrm{\Delta }t>0$,$$\begin{array}{c}\hfill \begin{array}{c}\Vert u^{k+1}{\Vert }_{\infty }\le B_3.\hfill \end{array}\end{array}$$

Proof See Appendix 1.

Lemma 4.4:

There exists positive constants $B_{4,1}$ and $B_{4,2}$, independent of $\mathrm{\Delta }x$ and $\mathrm{\Delta }t$, such that

• For all $j,s=0,\dots ,N-1$ and $k=0,\dots ,L-1$, $$\begin{array}{c}\hfill \begin{array}{c}|R_{j+1}^k-R_j^k|\le B_{4,1}\mathrm{\Delta }x,\left|{\mathit{\eta }}_{j+1,s}^k-{\mathit{\eta }}_{j,s}^k\right|\le B_{4,1}\mathrm{\Delta }x,\left|g_{j+1}^k-g_j^k\right|\le B_{4,1}\mathrm{\Delta }x.\hfill \end{array}\end{array}$$

• For all $j=1,\dots ,N-1$ and $k=0,\dots ,L-1$, $$\begin{array}{c}\hfill \begin{array}{c}\left|g_{j+1}^k-2g_j^k+g_{j-1}^k\right|\le B_{4,1}\mathrm{\Delta }{x}^2.\hfill \end{array}\end{array}$$

• For all $j,s=1,\dots ,N$ and $k=0,\dots ,L-1$, $$\begin{array}{c}\hfill \begin{array}{c}{\sum }_{\ell =1}^m\left|y_{\ell }^{k+1}-y_{\ell }^k\right|\le B_{4,1}\mathrm{\Delta }t,\left|{\mathit{\varphi }}^{k+1}-{\mathit{\varphi }}^k\right|\le B_{4,1}\mathrm{\Delta }t,\left|g_j^{k+1}-g_j^k\right|\le B_{4,1}\mathrm{\Delta }t,\hfill \\ \left|{\mathit{\eta }}_{j,s}^{k+1}-{\mathit{\eta }}_{j,s}^k\right|\le B_{4,1}\mathrm{\Delta }t,\left|{\mathit{\beta }}_j^{k+1}-{\mathit{\beta }}_j^k\right|\le B_{4,1}\mathrm{\Delta }t,\left|F^{k+1}-F^k\right|\le B_{4,1}\mathrm{\Delta }t.\hfill \end{array}\end{array}$$

• For all $k=0,\dots ,L-1$, $$\begin{array}{c}\hfill \begin{array}{c}(1+\frac{\mathrm{\Delta }t}{\mathrm{\Delta }x}{g}_0^k)|u_1^{k+1}-u_0^{k+1}|\le |u_1^k-u_0^k|+B_3(c_0{B}_2+B_{4,1}+B_{1,1})\mathrm{\Delta }t\hfill \\ +B_{4,2}(1+TV(u^k))\mathrm{\Delta }t.\hfill \end{array}\end{array}$$

Proof See Appendix 1.

Next we establish an important property of the numerical scheme. In particular, we show that the scheme has bounded total variation.

Lemma 4.5:

If $\mathrm{\Delta }t$ satisfies ${\mathit{\omega }}_1\mathrm{\Delta }t<1$, there exists a constant $B_5$ such that$$\begin{array}{c}\hfill \begin{array}{c}TV(u^{k+1})\le B_5,\hfill \end{array}\end{array}$$

for every $k\ge 0$.

Proof See Appendix 1.

We now show that our difference approximations $u_j^k$ and $y_{\ell }^k$ are ${\ell }^1$-Lipschitz in t.

Lemma 4.6:

There exists a constant $B_6>0$, independent of $\mathrm{\Delta }x$ and $\mathrm{\Delta }t$ such that for $r>q$,$$\begin{array}{c}\hfill \begin{array}{c}{\sum }_{j=1}^N\left|\frac{u_j^r-u_j^q}{\mathrm{\Delta }t}\right|\mathrm{\Delta }x\le B_6(r-q),\hfill \\ {\sum }_{\ell =1}^n\left|\frac{y_{\ell }^r-y_{\ell }^q}{\mathrm{\Delta }t}\right|\le B_6(r-q).\hfill \end{array}\end{array}$$

Proof See Appendix 1.

5. Existence and uniqueness of the weak solution

Following [41] (Chapter 16, p.275–279) we define a family of functions ${\{U_{\mathrm{\Delta }x,\mathrm{\Delta }t}(x,t)\}}_{(x,t)\in {\mathbb{D}}_x}$ by extending the finite difference approximations $u_j^k$ in (4.1(4.1) $$\begin{array}{cc}& \frac{u_j^{k+1}-u_j^k}{\mathrm{\Delta }t}+\frac{g_j^k{u}_j^{k+1}-g_{j-1}^k{u}_{j-1}^{k+1}}{\mathrm{\Delta }x}+u_j^{k+1}\sum _{m=1}^N{\mathit{\eta }}_{m,j}^k{u}_m^k\mathrm{\Delta }x+R_j^k{u}_j^{k+1}\hfill \\ \hfill & =\frac{1}{2}\sum _{s=1}^{j-1}{\mathit{\eta }}_{j-s,s}^k{u}_{j-s}^k{u}_s^{k+1}\mathrm{\Delta }x,\hfill \\ \hfill & g_0^k{u}_0^{k+1}=F^k+\sum _{j=1}^N{\mathit{\beta }}_j^k{u}_j^k\mathrm{\Delta }x,u_j^0=u_0(x_j),\hfill \\ \hfill & \frac{y_{\ell }^{k+1}-y_{\ell }^k}{\mathrm{\Delta }t}=f_{\ell }^k-h_{\ell }^k{y}_{\ell }^{k+1},y_{\ell }^0=y_{\ell ,0},\ell =1,\dots ,m.\hfill \end{array}$$(4.1) ) to piecewise constant functions and a family of functions ${\{Y_{\ell ,\mathrm{\Delta }t}(t)\}}_{t\in [0,T]}$ by extending the finite difference approximations $y_{\ell }^k$ in (4.1(4.1) $$\begin{array}{cc}& \frac{u_j^{k+1}-u_j^k}{\mathrm{\Delta }t}+\frac{g_j^k{u}_j^{k+1}-g_{j-1}^k{u}_{j-1}^{k+1}}{\mathrm{\Delta }x}+u_j^{k+1}\sum _{m=1}^N{\mathit{\eta }}_{m,j}^k{u}_m^k\mathrm{\Delta }x+R_j^k{u}_j^{k+1}\hfill \\ \hfill & =\frac{1}{2}\sum _{s=1}^{j-1}{\mathit{\eta }}_{j-s,s}^k{u}_{j-s}^k{u}_s^{k+1}\mathrm{\Delta }x,\hfill \\ \hfill & g_0^k{u}_0^{k+1}=F^k+\sum _{j=1}^N{\mathit{\beta }}_j^k{u}_j^k\mathrm{\Delta }x,u_j^0=u_0(x_j),\hfill \\ \hfill & \frac{y_{\ell }^{k+1}-y_{\ell }^k}{\mathrm{\Delta }t}=f_{\ell }^k-h_{\ell }^k{y}_{\ell }^{k+1},y_{\ell }^0=y_{\ell ,0},\ell =1,\dots ,m.\hfill \end{array}$$(4.1) ) to piecewise linear functions as follows:(5.1) $$\begin{array}{c}\hfill \begin{array}{c}{U}_{\mathrm{\Delta }x,\mathrm{\Delta }t}\left(x,t\right)=u_j^k\text{for}x\in \left[x_{j-1},x_j\right),t\in [t^{k-1},t^k),\hfill \\ Y_{\ell ,\mathrm{\Delta }t}(t)=y_{\ell }^{k-1}+\frac{y_{\ell }^k-y_{\ell }^{k-1}}{\mathrm{\Delta }t}(t-t^{k-1}),\text{for}t\in [t^{k-1},t^k),\hfill \end{array}\end{array}$$(5.1) $j=1,\dots ,N,k=1,2,\dots ,L.$ Then the set of functions $\{U_{\mathrm{\Delta }x,\mathrm{\Delta }t}(x,t)\}$ is (sequentially) compact in the topology of $L^1({\mathbb{D}}_x)$ and the set of functions $Y_{\ell ,\mathrm{\Delta }t}(t)$ is (sequentially) compact in the topology of C([0, T]) and we have the following Theorem.

Theorem 5.1:

For $\ell =1,\dots ,m$, there exists a subsequence $\{U_{\mathrm{\Delta }{x}_i,\mathrm{\Delta }{t}_i}\}\subset \{U_{\mathrm{\Delta }x,\mathrm{\Delta }t}(x,t)\}$ which converges to a function of bounded total variation $u\in L_{\infty }({\mathbb{D}}_x)$, and a subsequence $\{Y_{\ell ,\mathrm{\Delta }{t}_j}{\}}_{j=1}^{\infty }\subset \{Y_{\ell ,\mathrm{\Delta }t}\}$ which converges to a function $y_{\ell }\in C[0,T]$ in the sense that for all $t>0$ $$\begin{array}{c}\hfill \begin{array}{c}{\int }_0^{x_{max}}|U_{\mathrm{\Delta }{x}_i,\mathrm{\Delta }{t}_i}(x,t)-u(x,t)|\mathrm{d}x\to 0,\hfill \\ {\int }_0^T{\int }_0^{x_{max}}|U_{\mathrm{\Delta }{x}_i,\mathrm{\Delta }{t}_i}(x,t)-u(x,t)|\mathrm{d}x\mathrm{d}t\to 0,\hfill \end{array}\end{array}$$

and$$\begin{array}{c}\hfill \begin{array}{c}\Vert Y_{\ell ,\mathrm{\Delta }{t}_j}-y_{\ell }{\Vert }_{C[0,T]}\to 0,\hfill \end{array}\end{array}$$

as $i,j\to \infty$.

Proof See Appendix 1.

Theorem 5.2:

The limit functions u and $\overrightarrow{y}$ constructed by the finite difference scheme in Theorem () define a weak solution of the model (2.1(2.1) $$\begin{array}{c}\hfill \begin{array}{c}\frac{\mathit{\partial }u}{\mathit{\partial }t}+\frac{\mathit{\partial }}{\mathit{\partial }x}\left(g(x,t,\mathit{\varphi },\overrightarrow{y})u\right)+R(x,t,\mathit{\varphi },\overrightarrow{y})u=A(u),(x,t)\in (0,x_{max}]\times (0,T],\hfill \\ \frac{\mathrm{d}{y}_{\ell }}{\mathrm{d}t}=f_{\ell }(t,{\mathit{\sigma }}_{\ell },\overrightarrow{y})-h_{\ell }(t,{\mathit{\sigma }}_{\ell },\overrightarrow{y})y_{\ell },t\in (0,T],\ell =1,\dots ,m,\hfill \\ g(0,t,\mathit{\varphi },\overrightarrow{y})u(0,t)=F(t,\overrightarrow{y})+{\int }_0^{x_{\text{max}}}\mathit{\beta }(x,t,\mathit{\varphi },\overrightarrow{y})u(x,t)\mathrm{d}x,t\in (0,T],\hfill \\ u(x,0)=u^0(x),x\in [0,x_{max}],y_{\ell }(0)=y_{\ell ,0}\ge 0,\ell =1,\dots ,m,\hfill \end{array}\end{array}$$(2.1) ) and satisfy(5.2) $$\begin{array}{c}\hfill \begin{array}{c}{\Vert u(·,t)\Vert }_{L^1([0,x_{max}])}+{\sum }_{\ell =1}^m|y_{\ell }(t)|\le B_2,\hfill \\ {\Vert u(·,t)\Vert }_{L^{\infty }([0,x_{max}])}\le B_3,\hfill \\ {\Vert u(·,t)\Vert }_{BV([0,x_{max}])}\le B_7,\hfill \end{array}\end{array}$$(5.2)

for some constant $B_7$ and each $t\in [0,T]$.

Proof See Appendix 1.

Towards demonstrating that the weak solution is unique, we first establish the following lemma which guarantees the continuous dependence of the solutions $\{u_j^k\}$ and $\{{\overrightarrow{y}}^k\}$ to the coupled PDE system (2.1(2.1) $$\begin{array}{c}\hfill \begin{array}{c}\frac{\mathit{\partial }u}{\mathit{\partial }t}+\frac{\mathit{\partial }}{\mathit{\partial }x}\left(g(x,t,\mathit{\varphi },\overrightarrow{y})u\right)+R(x,t,\mathit{\varphi },\overrightarrow{y})u=A(u),(x,t)\in (0,x_{max}]\times (0,T],\hfill \\ \frac{\mathrm{d}{y}_{\ell }}{\mathrm{d}t}=f_{\ell }(t,{\mathit{\sigma }}_{\ell },\overrightarrow{y})-h_{\ell }(t,{\mathit{\sigma }}_{\ell },\overrightarrow{y})y_{\ell },t\in (0,T],\ell =1,\dots ,m,\hfill \\ g(0,t,\mathit{\varphi },\overrightarrow{y})u(0,t)=F(t,\overrightarrow{y})+{\int }_0^{x_{\text{max}}}\mathit{\beta }(x,t,\mathit{\varphi },\overrightarrow{y})u(x,t)\mathrm{d}x,t\in (0,T],\hfill \\ u(x,0)=u^0(x),x\in [0,x_{max}],y_{\ell }(0)=y_{\ell ,0}\ge 0,\ell =1,\dots ,m,\hfill \end{array}\end{array}$$(2.1) ) with respect to the initial conditions $u^0,{\overrightarrow{y}}^0$.

Lemma 5.3:

Let $\{u_j^k\}$, $\{{\widehat{u}}_j^k\}$, $\{{\overrightarrow{y}}^k\}$ and $\{{\overrightarrow{\widehat{y}}}^k\}$ be the solutions of (2.1(2.1) $$\begin{array}{c}\hfill \begin{array}{c}\frac{\mathit{\partial }u}{\mathit{\partial }t}+\frac{\mathit{\partial }}{\mathit{\partial }x}\left(g(x,t,\mathit{\varphi },\overrightarrow{y})u\right)+R(x,t,\mathit{\varphi },\overrightarrow{y})u=A(u),(x,t)\in (0,x_{max}]\times (0,T],\hfill \\ \frac{\mathrm{d}{y}_{\ell }}{\mathrm{d}t}=f_{\ell }(t,{\mathit{\sigma }}_{\ell },\overrightarrow{y})-h_{\ell }(t,{\mathit{\sigma }}_{\ell },\overrightarrow{y})y_{\ell },t\in (0,T],\ell =1,\dots ,m,\hfill \\ g(0,t,\mathit{\varphi },\overrightarrow{y})u(0,t)=F(t,\overrightarrow{y})+{\int }_0^{x_{\text{max}}}\mathit{\beta }(x,t,\mathit{\varphi },\overrightarrow{y})u(x,t)\mathrm{d}x,t\in (0,T],\hfill \\ u(x,0)=u^0(x),x\in [0,x_{max}],y_{\ell }(0)=y_{\ell ,0}\ge 0,\ell =1,\dots ,m,\hfill \end{array}\end{array}$$(2.1) ) corresponding to the initial conditions $u_{\text{i}}^0$, ${\widehat{u}}_{\text{i}}^0$, ${\overrightarrow{y}}^0$, and ${\overrightarrow{\widehat{y}}}^0$ respectively. Assume that $\mathrm{\Delta }t$ is chosen so that $({\mathit{\omega }}_1+c_0{B}_2/2)\mathrm{\Delta }t<1$. Then there exists a constant $B_8>0$ such that$$\begin{array}{c}\hfill \begin{array}{cc}(\Vert v^{k+1}{\Vert }_1+{\sum }_{\ell =1}^m|y_{\ell }^{k+1}-{\widehat{y}}_{\ell }^{k+1}|)\hfill & \le e^{B_{8}T}(\Vert v^0{\Vert }_1+{\sum }_{\ell =1}^m|y_{\ell }^0-{\widehat{y}}_{\ell }^0|).\hfill \end{array}\end{array}$$

Proof See Appendix 1.

Lemma 5.3 is then applied to prove that the weak solution defined in Theorems () and () is unique. Using techniques similar to those presented in [48], and the results established in Lemma 5.3, we arrive at the following result for some constant $\tilde{C}$,(5.3) $$\begin{array}{c}\hfill \begin{array}{cc}& \Vert u(·,t)-\widehat{u}{(·,t)\Vert }_1+\sum _{\ell =1}^m|y_{\ell }(·,t)-{\widehat{y}}_{\ell }(·,t)|\hfill \\ \hfill & \le e^{2CT}[\Vert u(·,0)-\widehat{u}{(·,0)\Vert }_1+\sum _{\ell =1}^m|y_{\ell }(0)-{\widehat{y}}_{\ell }(0)|\hfill \\ \hfill & +\tilde{C}{\int }_0^t(\Vert u(·,s)-\widehat{u}{(·,s)\Vert }_1+\sum _{\ell =1}^m|y_{\ell }(s)-{\widehat{y}}_{\ell }(s)|)\mathrm{d}s].\hfill \end{array}\end{array}$$(5.3)

Theorem 5.4:

Suppose u, $\overrightarrow{y}$ and $\widehat{u}$, $\widehat{\overrightarrow{y}}$ are two bounded variation weak solutions of (2.1(2.1) $$\begin{array}{c}\hfill \begin{array}{c}\frac{\mathit{\partial }u}{\mathit{\partial }t}+\frac{\mathit{\partial }}{\mathit{\partial }x}\left(g(x,t,\mathit{\varphi },\overrightarrow{y})u\right)+R(x,t,\mathit{\varphi },\overrightarrow{y})u=A(u),(x,t)\in (0,x_{max}]\times (0,T],\hfill \\ \frac{\mathrm{d}{y}_{\ell }}{\mathrm{d}t}=f_{\ell }(t,{\mathit{\sigma }}_{\ell },\overrightarrow{y})-h_{\ell }(t,{\mathit{\sigma }}_{\ell },\overrightarrow{y})y_{\ell },t\in (0,T],\ell =1,\dots ,m,\hfill \\ g(0,t,\mathit{\varphi },\overrightarrow{y})u(0,t)=F(t,\overrightarrow{y})+{\int }_0^{x_{\text{max}}}\mathit{\beta }(x,t,\mathit{\varphi },\overrightarrow{y})u(x,t)\mathrm{d}x,t\in (0,T],\hfill \\ u(x,0)=u^0(x),x\in [0,x_{max}],y_{\ell }(0)=y_{\ell ,0}\ge 0,\ell =1,\dots ,m,\hfill \end{array}\end{array}$$(2.1) ) corresponding to the initial conditions $u^0$,${\overrightarrow{y}}^0$ and ${\widehat{u}}^0$,${\widehat{\overrightarrow{y}}}^0$ respectively. Then there exists a constant $B_9$ such that(5.4) $$\begin{array}{c}\hfill \begin{array}{cc}& \Vert u(·,t)-\widehat{u}{(·,t)\Vert }_1+\sum _{\ell =1}^m|y_{\ell }(t)-{\widehat{y}}_{\ell }(t)|\le B_9({\Vert u(·,0)-\widehat{u}(·,0)\Vert }_1+\sum _{\ell =1}^m|y_{\ell }(0)-{\widehat{y}}_{\ell }(0)|).\hfill \end{array}\end{array}$$(5.4)

Proof The result follows immediately by applying the Gronwall inequality to Equation (5.3(5.3) $$\begin{array}{c}\hfill \begin{array}{cc}& \Vert u(·,t)-\widehat{u}{(·,t)\Vert }_1+\sum _{\ell =1}^m|y_{\ell }(·,t)-{\widehat{y}}_{\ell }(·,t)|\hfill \\ \hfill & \le e^{2CT}[\Vert u(·,0)-\widehat{u}{(·,0)\Vert }_1+\sum _{\ell =1}^m|y_{\ell }(0)-{\widehat{y}}_{\ell }(0)|\hfill \\ \hfill & +\tilde{C}{\int }_0^t(\Vert u(·,s)-\widehat{u}{(·,s)\Vert }_1+\sum _{\ell =1}^m|y_{\ell }(s)-{\widehat{y}}_{\ell }(s)|)\mathrm{d}s].\hfill \end{array}\end{array}$$(5.3) ) and letting $B_9=e^{2CT+\tilde{C}T{e}^{2CT}}$.

6. Development of sub-models for the vital rates

In this section we parameterize the model (2.1(2.1) $$\begin{array}{c}\hfill \begin{array}{c}\frac{\mathit{\partial }u}{\mathit{\partial }t}+\frac{\mathit{\partial }}{\mathit{\partial }x}\left(g(x,t,\mathit{\varphi },\overrightarrow{y})u\right)+R(x,t,\mathit{\varphi },\overrightarrow{y})u=A(u),(x,t)\in (0,x_{max}]\times (0,T],\hfill \\ \frac{\mathrm{d}{y}_{\ell }}{\mathrm{d}t}=f_{\ell }(t,{\mathit{\sigma }}_{\ell },\overrightarrow{y})-h_{\ell }(t,{\mathit{\sigma }}_{\ell },\overrightarrow{y})y_{\ell },t\in (0,T],\ell =1,\dots ,m,\hfill \\ g(0,t,\mathit{\varphi },\overrightarrow{y})u(0,t)=F(t,\overrightarrow{y})+{\int }_0^{x_{\text{max}}}\mathit{\beta }(x,t,\mathit{\varphi },\overrightarrow{y})u(x,t)\mathrm{d}x,t\in (0,T],\hfill \\ u(x,0)=u^0(x),x\in [0,x_{max}],y_{\ell }(0)=y_{\ell ,0}\ge 0,\ell =1,\dots ,m,\hfill \end{array}\end{array}$$(2.1) ) to numerically investigate the role of the phytoplankton aggregate growth rate process on algae bloom development and nutrient dynamics. The basis of this biogeochemical application is the 1993 SIGMA reactor experiment [23,49]. The researchers employed a version of the model (2.1(2.1) $$\begin{array}{c}\hfill \begin{array}{c}\frac{\mathit{\partial }u}{\mathit{\partial }t}+\frac{\mathit{\partial }}{\mathit{\partial }x}\left(g(x,t,\mathit{\varphi },\overrightarrow{y})u\right)+R(x,t,\mathit{\varphi },\overrightarrow{y})u=A(u),(x,t)\in (0,x_{max}]\times (0,T],\hfill \\ \frac{\mathrm{d}{y}_{\ell }}{\mathrm{d}t}=f_{\ell }(t,{\mathit{\sigma }}_{\ell },\overrightarrow{y})-h_{\ell }(t,{\mathit{\sigma }}_{\ell },\overrightarrow{y})y_{\ell },t\in (0,T],\ell =1,\dots ,m,\hfill \\ g(0,t,\mathit{\varphi },\overrightarrow{y})u(0,t)=F(t,\overrightarrow{y})+{\int }_0^{x_{\text{max}}}\mathit{\beta }(x,t,\mathit{\varphi },\overrightarrow{y})u(x,t)\mathrm{d}x,t\in (0,T],\hfill \\ u(x,0)=u^0(x),x\in [0,x_{max}],y_{\ell }(0)=y_{\ell ,0}\ge 0,\ell =1,\dots ,m,\hfill \end{array}\end{array}$$(2.1) ) which was used to investigate the effects of coagulation on the population dynamics of a phytoplankton community (mostly diatoms, i.e. Chaetoceros sp., Thalassiosira sp.) subject to nutrient depletion. The model results matched the observations well, and size-structured model parameters were made available as a result of this experiment. The key parameters under investigation in the SIGMA experiment were the stickiness parameter and the contact efficiency power relationship. An individual based approach using only single cells was used to track nutrient depletion, and nutrient uptake was assumed to be a function of the single cells only. Our goal is to enhance the SIGMA reactor aggregation model [23] to include aggregate growth and reproduction and to demonstrate the effect of including these terms on the coupled particle and nutrient dynamics.

In order to extend the model presented in [23] we assume that cell aggregates in general can be partitioned into living material on the particle surface and non-living material at the core. The analysis of a size-structured model which includes growth and aggregation has been investigated using the finite element method (FEM) (see, [25]). The numerical results established in [25] suggest that growth and aggregation are processes that enhance each other and the inclusion of growth increases the formation of large aggregates compared to models which include aggregation only. The following analysis extends the previous work by investigating the role of multi-celled aggregate growth, coagulation and nutrient availability in a realistic setting. We begin by deriving concrete expressions for the vital rate parameters which are necessary to perform the numerical trials.

Note that in this application, $R=0$ since the reactor system is closed and well-mixed throughout the experiment (i.e. negligible sedimentation effects), and the effects of grazing were negligible. Also, $F=0$ throughout the simulation since the tank was not artificially seeded with single cells after the beginning of the experiment. We also point out that we understand the slight abuse of notation in this section because some parameters described herein are not as general as those described in (2.1(2.1) $$\begin{array}{c}\hfill \begin{array}{c}\frac{\mathit{\partial }u}{\mathit{\partial }t}+\frac{\mathit{\partial }}{\mathit{\partial }x}\left(g(x,t,\mathit{\varphi },\overrightarrow{y})u\right)+R(x,t,\mathit{\varphi },\overrightarrow{y})u=A(u),(x,t)\in (0,x_{max}]\times (0,T],\hfill \\ \frac{\mathrm{d}{y}_{\ell }}{\mathrm{d}t}=f_{\ell }(t,{\mathit{\sigma }}_{\ell },\overrightarrow{y})-h_{\ell }(t,{\mathit{\sigma }}_{\ell },\overrightarrow{y})y_{\ell },t\in (0,T],\ell =1,\dots ,m,\hfill \\ g(0,t,\mathit{\varphi },\overrightarrow{y})u(0,t)=F(t,\overrightarrow{y})+{\int }_0^{x_{\text{max}}}\mathit{\beta }(x,t,\mathit{\varphi },\overrightarrow{y})u(x,t)\mathrm{d}x,t\in (0,T],\hfill \\ u(x,0)=u^0(x),x\in [0,x_{max}],y_{\ell }(0)=y_{\ell ,0}\ge 0,\ell =1,\dots ,m,\hfill \end{array}\end{array}$$(2.1) ).

Consider an individual phytoplankton cell aggregate of volume x subject to cell growth, coagulation, and encapsulation. Encapsulation is assumed to occur when cell aggregates become sufficiently large in size to the point where cell division is arrested in portions of the aggregate (Figure 1). As indicated in [25], the maximum specific growth rate of single cells based on experimental data is approximately 1.2 per day where the specific growth rate is the ratio of the growth rate to the total size of the aggregate. Encapsulation would then be reflected in specific growth rates smaller than unity based on the assumption that a necrotic portion (i.e. core) of the aggregate does not undergo photosynthesis (e.g. as a result of inadequate light exposure or nutrients). This assumption has been applied to the numerical investigation of aggregating bacterial systems in [39]. The following derivation is in the spirit of the arguments presented in [22,23]. The forthcoming assumptions are made in an effort to develop biologically relevant expressions of the size structured vital rate parameters:

(N1)	The total volume of a cell aggregate can be partitioned into core volume and active volume, i.e. $x=x_c+x_a$.
(N2)	The core volume $x_c$ does not undergo photosynthesis or cell division.
(N3)	The active volume, $x_a$ consists of the non-core volume. The active volume is proportional to the total size of the aggregate and the proportionality factor is termed the active fraction $K_a$ which is assumed to depend on size, i.e. $x_a=K_a(x)x$ where $0\le K_a(x)\le 1$. We also assume that the active fraction is a twice continuously differentiable monotone function satisfying $K_a(0)=1$ and $K_a(x_{max})=0$.
(N4)	Nutrient uptake and cell division occurs only in the active portion of the aggregate. The nutrient uptake rate is proportional to the total active biomass in the system.
(N5)	The total aggregate growth rate, $g_b$ is defined as the rate of increase in size as a result of cell division. Aggregate growth is linearly related to the active volume of the aggregate.
(N6)	The total aggregate growth rate is the sum of attached growth $g_a$ and offspring growth $g_o$, i.e. $g_b=g_a+g_o$. The offspring growth rate measures the production of newborn aggregates which fall off the parent aggregate. Offspring contributes to the reproduction term and consists entirely of living cells.

Based on (N5) we assume the total aggregate growth rate (uninhibited by nutrient limitation) has the following form$$\begin{array}{c}\hfill \begin{array}{c}{g}_b(x)=K_a(x)(g_{\text{min}}+{\mathit{\kappa }}_{0}x),\hfill \end{array}\end{array}$$

where $K_a(x)=1-{(x/x_{max})}^{P_a}$ and the exponent $P_a\ge 1$ is a parameter which determines the influence of encapsulation with larger values corresponding to greater resistance to the effects of encapsulation. The parameter $g_{\text{min}}$ represents the growth rate of the smallest sized aggregate considered in the analysis. The constant ${\mathit{\kappa }}_0$ is a parameter which also reflects the effects of encapsulation with larger values corresponding to a greater resistance to the effects of encapsulation. The attached growth rate is defined by $g_a=\mathit{\lambda }(x)g_b(x)$ where $\mathit{\lambda }(x)=1-x/x_{max}$. Conversely, the reproduction rate has the following form$$\begin{array}{c}\hfill \begin{array}{c}\mathit{\beta }(x)={\mathit{\beta }}_f(1-\mathit{\lambda }(x))g_b(x),\hfill \end{array}\end{array}$$

where ${\mathit{\beta }}_f$ is a reproduction factor.

Nutrient concentration was the only environmental variable considered in the trial and hence $m=1$ in (2.1(2.1) $$\begin{array}{c}\hfill \begin{array}{c}\frac{\mathit{\partial }u}{\mathit{\partial }t}+\frac{\mathit{\partial }}{\mathit{\partial }x}\left(g(x,t,\mathit{\varphi },\overrightarrow{y})u\right)+R(x,t,\mathit{\varphi },\overrightarrow{y})u=A(u),(x,t)\in (0,x_{max}]\times (0,T],\hfill \\ \frac{\mathrm{d}{y}_{\ell }}{\mathrm{d}t}=f_{\ell }(t,{\mathit{\sigma }}_{\ell },\overrightarrow{y})-h_{\ell }(t,{\mathit{\sigma }}_{\ell },\overrightarrow{y})y_{\ell },t\in (0,T],\ell =1,\dots ,m,\hfill \\ g(0,t,\mathit{\varphi },\overrightarrow{y})u(0,t)=F(t,\overrightarrow{y})+{\int }_0^{x_{\text{max}}}\mathit{\beta }(x,t,\mathit{\varphi },\overrightarrow{y})u(x,t)\mathrm{d}x,t\in (0,T],\hfill \\ u(x,0)=u^0(x),x\in [0,x_{max}],y_{\ell }(0)=y_{\ell ,0}\ge 0,\ell =1,\dots ,m,\hfill \end{array}\end{array}$$(2.1) ). From hereinafter, the symbol $y_1(t)=y(t)$ is used to model the average nutrient concentration in the tank at time t. The nutrient reduction rate in the reactor at time t is assumed to be proportional to the active biomass in the following manner,(6.1) $$\begin{array}{c}\hfill \begin{array}{c}\frac{\mathrm{d}y}{\mathrm{d}t}=-N_f{\int }_0^{x_{max}}{K}_a(x)xu(x,t)\mathrm{d}x.\hfill \end{array}\end{array}$$(6.1)

Here $N_f$ is a nutrient consumption parameter. We then modify (6.1(6.1) $$\begin{array}{c}\hfill \begin{array}{c}\frac{\mathrm{d}y}{\mathrm{d}t}=-N_f{\int }_0^{x_{max}}{K}_a(x)xu(x,t)\mathrm{d}x.\hfill \end{array}\end{array}$$(6.1) ) by imposing a Monod-type function to model the self-limiting effects of nutrient consumption during bloom formation to obtain(6.2) $$\begin{array}{c}\hfill \begin{array}{cc}& \frac{\mathrm{d}y}{\mathrm{d}t}=-{\mathrm{\Gamma }}_1(y)N_f{\int }_0^{x_{max}}{K}_a(x)xu(x,t)\mathrm{d}x,{\mathrm{\Gamma }}_1(y)=\frac{y}{k_y+y},\hfill \end{array}\end{array}$$(6.2)

where ${\mathrm{\Gamma }}_1(y)$ is the nutrient consumption rate limiter. We also apply a growth rate limiter ${\mathrm{\Gamma }}_2(y;\mathit{\epsilon })={\mathrm{\Gamma }}_1(y)+{\mathit{\epsilon }}_y$ to the growth function where ${\mathit{\epsilon }}_y$ is a small constant which reflects the presence of alternative nutrients once nitrates are depleted (e.g. Figure 5, [49]). Using (6.2(6.2) $$\begin{array}{c}\hfill \begin{array}{cc}& \frac{\mathrm{d}y}{\mathrm{d}t}=-{\mathrm{\Gamma }}_1(y)N_f{\int }_0^{x_{max}}{K}_a(x)xu(x,t)\mathrm{d}x,{\mathrm{\Gamma }}_1(y)=\frac{y}{k_y+y},\hfill \end{array}\end{array}$$(6.2) ) we are able to simulate the nutrient dynamics coupled with the size-structured phytoplankton dynamics in the reactor. The particle density is modeled using the approach presented in [22] and is coupled with (6.2(6.2) $$\begin{array}{c}\hfill \begin{array}{cc}& \frac{\mathrm{d}y}{\mathrm{d}t}=-{\mathrm{\Gamma }}_1(y)N_f{\int }_0^{x_{max}}{K}_a(x)xu(x,t)\mathrm{d}x,{\mathrm{\Gamma }}_1(y)=\frac{y}{k_y+y},\hfill \end{array}\end{array}$$(6.2) ) through the nutrient-dependent growth and birth rate functions. The nutrient dependent attached growth rate (Figure 1) is given by(6.3) $$\begin{array}{c}\hfill \begin{array}{c}g(x,y)={\mathrm{\Gamma }}_2(y)g_a(x)0\le x\le x_{max},\hfill \end{array}\end{array}$$(6.3)

and the reproduction term (i.e. left boundary condition) is given by$$\begin{array}{c}\hfill \begin{array}{c}g(0,y)u(0,t)={\mathit{\beta }}_f{\mathrm{\Gamma }}_2(y){\int }_0^{x_{max}}(1-\mathit{\lambda }(\mathit{\xi }))\left(K_a(\mathit{\xi })(g_{\text{min}}+{\mathit{\kappa }}_0\mathit{\xi })\right)u(\mathit{\xi },t)\mathrm{d}\mathit{\xi }.\hfill \end{array}\end{array}$$

The coagulation function is based on the formulation given in [23],$$\begin{array}{c}\hfill \begin{array}{c}\mathit{\eta }(x,s)={\mathit{\eta }}_f(b_s(x,s)+b_d(x,s)).\hfill \end{array}\end{array}$$

Here, ${\mathit{\eta }}_f$ is the sticking probability, and $b_s$ and $b_d$ are the individual collision rates due to shear and differential settling of particles, respectively as defined below [23,24],$$\begin{array}{c}\hfill \begin{array}{cc}& b_s(x,s)=0.163{\left(\frac{\mathit{\epsilon }}{\mathit{\nu }}\right)}^{1/2}(d_x+d_s)E{C}_s(x,s),b_d(x,s)=\frac{\mathit{\pi }{(d_x+d_s)}^2}{4}|w_x-w_s|E{C}_{x,s}^d\hfill \\ \hfill & E{C}_s(x,s)=E{C}_{0}exp\left(\frac{-0.1}{d_{\text{min}}/d_{max}}\right),E{C}_0=\frac{1}{exp(-0.1)},\hfill \\ \hfill & E{C}_d(x,s)=\frac{{(d_{\text{min}}/d_{max})}^{1.3}}{2{(1+d_{\text{min}}/d_{max})}^{1.3}},d_x=2{\left(\frac{3x}{4\mathit{\pi }}\right)}^{1/3}\hfill \\ \hfill & d_{\text{min}}=min\{d_x,d_s\},d_{max}=max\{d_x,d_s\},w_x=f_{\mathit{\rho }}\mathit{\xi }{\left(\frac{d_x}{2}\right)}^{1.17}.\hfill \end{array}\end{array}$$

Figure 1. Modeled phytoplankton aggregate growth rates. (a) Particle growth parameters, (b) Particle doubling rates.

Figure 2. Coagulation rates, $\mathit{\eta }(x,s)$.

Here, $E{C}_s$ and $E{C}_d$ denote the contact efficiencies due to shear and differential settling respectively. The parameter $w_x$ represents the fall velocity of a particle of volume x and $d_x$ represents the diameter of a particle having volume x. See Figure 2 for a graph of the coagulation rates. The remaining constants and variables used in the model are summarized in Table 1.

Table 1. Variables and parameter values used in the model.

Symbol	Interpretation	Value (Unit)	Source
u(x, t)	Density of aggregates	$(\#/\text{liter})/{\text{mm}}^3$
$x_{max}$	Volume of largest aggregate	$1.0({\text{mm}}^3)$	[23]
$x_c(x)$	Volume of core material	$({\text{mm}}^3)$
$x_a(x)$	Volume of active material	$({\text{mm}}^3)$
$K_a(x)$	Active fraction	Dimensionless
$g_b(x)$	Total aggregate growth rate	$({\text{mm}}^3/\text{day})$
$g_a(x)$	Attached growth rate	$({\text{mm}}^3/\text{day})$
$g_o(x)$	Offspring growth rate	$({\text{mm}}^3/\text{day})$
g(x, y)	Nutrient-limited growth rate	$({\text{mm}}^3/\text{day})$
${\mathit{\kappa }}_0$	Encapsulation factor	(1/day)	Estimated
$P_a$	Encapsulation exponent	Dimensionless	Estimated
$g_{\text{min}}$	Growth rate of smallest aggregate	$({\text{mm}}^3/\text{day})$	Estimated
${\mathrm{\Gamma }}_1(y)$	Nutrient reduction rate limiter	Dimensionless	[23]
${\mathrm{\Gamma }}_2(y)$	Growth rate limiter	Dimensionless
y(t)	Nitrate concentration	$(\mathit{\mu }\text{M})$
$k_y$	Half saturation constant	0.5 $(\mathrm{\mu }\text{M})$	[23]
$\mathit{\epsilon }$	Dissipation rate	${10}^{-6}$ $({\text{m}}^2{s}^{-3})$	[23]
$N_f$	Biomass nutrient consumption constant	$(\mathit{\mu }\text{M}/((\#/\text{liter})·{\text{mm}}^3·\text{day}))$	Estimated
${\mathit{\beta }}_f$	Reproduction factor	$({\text{mm}}^{-3})$	Estimated
$\mathit{\nu }$	Kinematic viscosity	$1.4\times {10}^{-6}$ $({\text{m}}^2{s}^{-1})$	[23]
${\mathit{\epsilon }}_y$	Alternate nutrient source	0.01 $(\mathit{\mu }\text{M})$
$\mathit{\eta }(x,s)$	Coagulation rate	$1/((\#/\text{liter})·\text{day})$
${\mathit{\eta }}_f$	Sticking probability	Dimensionless	Estimated
$b_s(x,s)$	Shear collision rate	(1/day)	[23]
$b_d(x,s)$	Differential settling collision rate	(1/day)	[23]
$d_x$	Particle diameter	(mm)	[23]
$E{C}_s$	Shear contact efficiency	Dimensionless	[23]
$E{C}_d$	Differential settling contact efficiency	Dimensionless	[23]
$E{C}_0$	Shear contact efficiency constant	Dimensionless	[23]
$w_x$	Fall velocity of a particle	cm/s	[23]
$f_{\mathit{\rho }}$	Fall velocity scale parameter	2; Dimensionless	[24]
$\mathit{\xi }$	Curve fit for w parameter	1.24; $({\text{cm}}^{-0.17}{\text{s}}^{-1})$	[24]
p	Contact efficiency exponent	1.3; Dimensionless	[23]

7. Least-squares method and parameter estimation results

In this section, we apply the vital rates and nutrient consumption formulation outlined in Section 6 to the model (2.1(2.1) $$\begin{array}{c}\hfill \begin{array}{c}\frac{\mathit{\partial }u}{\mathit{\partial }t}+\frac{\mathit{\partial }}{\mathit{\partial }x}\left(g(x,t,\mathit{\varphi },\overrightarrow{y})u\right)+R(x,t,\mathit{\varphi },\overrightarrow{y})u=A(u),(x,t)\in (0,x_{max}]\times (0,T],\hfill \\ \frac{\mathrm{d}{y}_{\ell }}{\mathrm{d}t}=f_{\ell }(t,{\mathit{\sigma }}_{\ell },\overrightarrow{y})-h_{\ell }(t,{\mathit{\sigma }}_{\ell },\overrightarrow{y})y_{\ell },t\in (0,T],\ell =1,\dots ,m,\hfill \\ g(0,t,\mathit{\varphi },\overrightarrow{y})u(0,t)=F(t,\overrightarrow{y})+{\int }_0^{x_{\text{max}}}\mathit{\beta }(x,t,\mathit{\varphi },\overrightarrow{y})u(x,t)\mathrm{d}x,t\in (0,T],\hfill \\ u(x,0)=u^0(x),x\in [0,x_{max}],y_{\ell }(0)=y_{\ell ,0}\ge 0,\ell =1,\dots ,m,\hfill \end{array}\end{array}$$(2.1) ). We then utilize the finite difference scheme (4.1(4.1) $$\begin{array}{cc}& \frac{u_j^{k+1}-u_j^k}{\mathrm{\Delta }t}+\frac{g_j^k{u}_j^{k+1}-g_{j-1}^k{u}_{j-1}^{k+1}}{\mathrm{\Delta }x}+u_j^{k+1}\sum _{m=1}^N{\mathit{\eta }}_{m,j}^k{u}_m^k\mathrm{\Delta }x+R_j^k{u}_j^{k+1}\hfill \\ \hfill & =\frac{1}{2}\sum _{s=1}^{j-1}{\mathit{\eta }}_{j-s,s}^k{u}_{j-s}^k{u}_s^{k+1}\mathrm{\Delta }x,\hfill \\ \hfill & g_0^k{u}_0^{k+1}=F^k+\sum _{j=1}^N{\mathit{\beta }}_j^k{u}_j^k\mathrm{\Delta }x,u_j^0=u_0(x_j),\hfill \\ \hfill & \frac{y_{\ell }^{k+1}-y_{\ell }^k}{\mathrm{\Delta }t}=f_{\ell }^k-h_{\ell }^k{y}_{\ell }^{k+1},y_{\ell }^0=y_{\ell ,0},\ell =1,\dots ,m.\hfill \end{array}$$(4.1) ) to solve a least squares parameter estimation problem that compares the model output to real data. The objective is to obtain a suitable set of parameters from which we can gain insight on the interaction between aggregate growth, coagulation, and nutrient dynamics and evaluate the sensitivity of the model outputs near the optimal parameter values. Given that the model contains many parameters, such an exercise may be useful both in experimental design and also in identifying the most significant factors driving the overall dynamics of the system.

The focus of this section is the study of the following parameter identification problem: given observations $D_r^N$ and $D_{r,i}^Q$ which correspond to the measured daily average nutrient concentration and total particle number in the size interval $[{\mathit{\xi }}_i,{\mathit{\xi }}_{i+1}),i=1,2,\dots ,10,$ respectively at time ${\mathit{\tau }}_r,r=1,2,\dots ,11$ days from the beginning of the simulation, find a parameter $\mathit{\theta }=({\mathit{\theta }}_1,\dots ,{\mathit{\theta }}_6)$ which minimizes the least squares cost functional $f(\mathit{\theta })$ defined by(7.1) $$\begin{array}{c}\hfill \begin{array}{cc}f(\mathit{\theta })\hfill & ={\sum }_{r=1}^{11}{|D_r^N-y({\mathit{\tau }}_r;\mathit{\theta })|}^2\hfill \\ \hfill & +\mathit{\gamma }{\sum }_{r=1}^{11}{\sum }_{i=1}^{10}{|{log}_{10}(D_{r,i}^Q+1)-{log}_{10}({\int }_{{\mathit{\xi }}_{i-1}}^{{\mathit{\xi }}_i}u(s,{\mathit{\tau }}_r;\mathit{\theta })\mathrm{d}s+1)|}^2,\hfill \end{array}\end{array}$$(7.1)

where the subscript $i=1,2,\dots ,10$ denotes the population histogram bin corresponding to the particle size class ${\mathit{\xi }}_i\in [0,x_{max}]$. The vector of parameters to be estimated is $\mathit{\theta }=({\mathit{\theta }}_1,\cdots ,{\mathit{\theta }}_6)$, where ${\mathit{\theta }}_1=N_f,{\mathit{\theta }}_2={\mathit{\beta }}_f,{\mathit{\theta }}_3={\mathit{\eta }}_f,{\mathit{\theta }}_4={\mathit{\kappa }}_0,{\mathit{\theta }}_5=P_a$ and ${\mathit{\theta }}_6=g_{\text{min}}$. The factor $\mathit{\gamma }$ is introduced to balance the objective function since the quantities of the nutrients and phytoplankton have different units and relative magnitudes.

Figure 3. Nutrient model performance.

To compute minimizers for the above inverse problem we first approximate the objective function (7.1(7.1) $$\begin{array}{c}\hfill \begin{array}{cc}f(\mathit{\theta })\hfill & ={\sum }_{r=1}^{11}{|D_r^N-y({\mathit{\tau }}_r;\mathit{\theta })|}^2\hfill \\ \hfill & +\mathit{\gamma }{\sum }_{r=1}^{11}{\sum }_{i=1}^{10}{|{log}_{10}(D_{r,i}^Q+1)-{log}_{10}({\int }_{{\mathit{\xi }}_{i-1}}^{{\mathit{\xi }}_i}u(s,{\mathit{\tau }}_r;\mathit{\theta })\mathrm{d}s+1)|}^2,\hfill \end{array}\end{array}$$(7.1) ) in the following manner:$$\begin{array}{c}\hfill \begin{array}{cc}f{(\mathit{\theta })}_{\mathrm{\Delta }x,\mathrm{\Delta }t}\hfill & ={\sum }_{r=1}^{11}{|D_r^N-Y_{\mathrm{\Delta }t}({\mathit{\tau }}_r;\mathit{\theta })|}^2\hfill \\ \hfill & +\mathit{\gamma }{\sum }_{r=1}^{11}{|{log}_{10}(D_{r,i}^Q+1)-{log}_{10}({\int }_{{\mathit{\xi }}_{i-1}}^{{\mathit{\xi }}_i}{U}_{\mathrm{\Delta }x,\mathrm{\Delta }t}(s,{\mathit{\tau }}_r;\mathit{\theta })\mathrm{d}s+1)|}^2,\hfill \end{array}\end{array}$$

where the functions $Y_{\mathrm{\Delta }t}$ and $U_{\mathrm{\Delta }x,\mathrm{\Delta }t}({\mathit{\tau }}_r;\mathit{\theta })$ as defined in (5.1(5.1) $$\begin{array}{c}\hfill \begin{array}{c}{U}_{\mathrm{\Delta }x,\mathrm{\Delta }t}\left(x,t\right)=u_j^k\text{for}x\in \left[x_{j-1},x_j\right),t\in [t^{k-1},t^k),\hfill \\ Y_{\ell ,\mathrm{\Delta }t}(t)=y_{\ell }^{k-1}+\frac{y_{\ell }^k-y_{\ell }^{k-1}}{\mathrm{\Delta }t}(t-t^{k-1}),\text{for}t\in [t^{k-1},t^k),\hfill \end{array}\end{array}$$(5.1) ) are obtained from solving the difference equations (4.1(4.1) $$\begin{array}{cc}& \frac{u_j^{k+1}-u_j^k}{\mathrm{\Delta }t}+\frac{g_j^k{u}_j^{k+1}-g_{j-1}^k{u}_{j-1}^{k+1}}{\mathrm{\Delta }x}+u_j^{k+1}\sum _{m=1}^N{\mathit{\eta }}_{m,j}^k{u}_m^k\mathrm{\Delta }x+R_j^k{u}_j^{k+1}\hfill \\ \hfill & =\frac{1}{2}\sum _{s=1}^{j-1}{\mathit{\eta }}_{j-s,s}^k{u}_{j-s}^k{u}_s^{k+1}\mathrm{\Delta }x,\hfill \\ \hfill & g_0^k{u}_0^{k+1}=F^k+\sum _{j=1}^N{\mathit{\beta }}_j^k{u}_j^k\mathrm{\Delta }x,u_j^0=u_0(x_j),\hfill \\ \hfill & \frac{y_{\ell }^{k+1}-y_{\ell }^k}{\mathrm{\Delta }t}=f_{\ell }^k-h_{\ell }^k{y}_{\ell }^{k+1},y_{\ell }^0=y_{\ell ,0},\ell =1,\dots ,m.\hfill \end{array}$$(4.1) ). For the remainder of this section we will use the notation for simplicity,$$\begin{array}{c}\hfill \begin{array}{c}{Q}_i({\mathit{\tau }}_r;\mathit{\theta })={\int }_{{\mathit{\xi }}_{i-1}}^{{\mathit{\xi }}_i}u(s,{\mathit{\tau }}_r;\mathit{\theta })\mathrm{d}s,Q{}_i^{\mathrm{\Delta }x,\mathrm{\Delta }t}({\mathit{\tau }}_r;\mathit{\theta })={\int }_{{\mathit{\xi }}_{i-1}}^{{\mathit{\xi }}_i}{U}_{\mathrm{\Delta }x,\mathrm{\Delta }t}(s,{\mathit{\tau }}_r;\mathit{\theta })\mathrm{d}s.\hfill \end{array}\end{array}$$

For the computations presented below we divide the model domain $[0,x_{max}]\times [0,T]$ using uniformly spaced mesh points with $T=12$ days and $x_{max}=1.0$ ${\text{mm}}^3$. The computational point spacing along the x-axis was consistent with the numerical trials presented in [23], i.e. $N=100$, yielding $\mathrm{\Delta }x=0.01$ ${\text{mm}}^3$ and $x_j=j\mathrm{\Delta }x$ where $j=1,\dots ,N$. The size bins were defined as ${\mathit{\xi }}_i=x_{10i}$. The computational time step is selected as follows: $k=0,1,\dots ,L$ and $L=2400$ simulation steps yielding $\mathrm{\Delta }t=0.005,$ $t_k=0.005k,$ and ${\mathit{\tau }}_r=t_{200r}$ days respectively.

We point out that the raw experimental data was not directly available. As a result the observed nutrient data $D_r^N$ was digitized from the figures provided in [23]. Similarly, the observed particle histogram data $D_{r,i}^Q$ was inferred from the figures presented in [29] noting that the non-uniform measurements were converted to a uniform particle size bin in order to compare the measurements with the numerical model results. The initial conditions for both the nutrients and the particle densities were inferred from the digitized experimental data. A series of preliminary numerical trials were performed first to determine the weight factor $\mathit{\gamma }=13$. This was done for convenience, but in general the weight parameter $\mathit{\gamma }$ may be estimated from the data directly using the methods discussed in [50]. The least squares minimization problem was then numerically solved using the standard Fortran LMDIF1 minimization package.

The final squared residuals for the nutrients and the phytoplankton population were 36.4 and 22.4, respectively. This shows that our choice of $\mathit{\gamma }=13$ produces residuals which are within the same order of magnitude for both datasets. In Figures 3 and 4 we present the data fits for the nutrients and the phytoplankton aggregates that resulted from solving the above least-squares problem and in Table 2 we present the values of the estimated parameters.

The results suggest that the model performed reasonably well in capturing the overall behavior of the bloom, i.e. formation time for large aggregates (day 6) and order of magnitude for particles in the middle of the size spectrum (Figure 4). The nutrient model also performed reasonably well particularly in tracking the disappearance rate and the time at which nitrates fell below detection limits (day 9) inside the tank (Figure 3). However, the model overestimated the total particle numbers during the middle of the experiment but converged with the measurements towards the end of the simulation run (Figure 5). We also point out that the model tends to overestimate the total particles until day 8 at which time the model begins to underestimate the aggregates uniformly except for particles in the size range smaller than 0.2 ${\text{mm}}^3$. After this time the modeled histogram generally fell short of the measurements over the final days of the simulation. We point out that the total particle numbers converged at the end of the bloom which suggests the results may be improved in the future by including an additional process which has the ability to distribute particles uniformly across the size range. Such a process should theoretically operate throughout the course of the bloom even when growth and reproduction processes are restricted by a lack of nutrients. A more realistic scenario may be that cell aggregates fracture according to a distribution of sizes, particularly in mixed (i.e. turbulent) environments such as those maintained during the SIGMA experiment. Distributed breakage leads to a fragmentation term which should be investigated in future efforts to improve the model predictions.

The calculation of standard errors and sensitivity for the estimated parameters follows the methods outlined in [51]. In particular, to compute the variance in the estimated vector of parameters, $\mathit{\theta }$, we start by computing the sensitivity matrix given by$$Z(\mathit{\theta })=\left[\begin{array}{cccc}{y}_{{\mathit{\theta }}_1}({\mathit{\tau }}_1;\mathit{\theta })& y_{{\mathit{\theta }}_2}({\mathit{\tau }}_1;\mathit{\theta })& \cdots & y_{{\mathit{\theta }}_6}({\mathit{\tau }}_1;\mathit{\theta })\\ ⋮& ⋮& & ⋮\\ y_{{\mathit{\theta }}_1}({\mathit{\tau }}_{11};\mathit{\theta })& y_{{\mathit{\theta }}_2}({\mathit{\tau }}_{11};\mathit{\theta })& \cdots & y_{{\mathit{\theta }}_6}({\mathit{\tau }}_{11};\mathit{\theta })\\ \frac{\overline{c}{(Q_1)}_{{\mathit{\theta }}_1}({\mathit{\tau }}_1;\mathit{\theta })}{1+Q_1({\mathit{\tau }}_1;\mathit{\theta })}& \frac{\overline{c}{(Q_1)}_{{\mathit{\theta }}_2}({\mathit{\tau }}_1;\mathit{\theta })}{1+Q_1({\mathit{\tau }}_1;\mathit{\theta })}& \cdots & \frac{\overline{c}{(Q_1)}_{{\mathit{\theta }}_6}({\mathit{\tau }}_1;\mathit{\theta })}{1+Q_1({\mathit{\tau }}_1;\mathit{\theta })}\\ ⋮& ⋮& & ⋮\\ \frac{\overline{c}{(Q_{10})}_{{\mathit{\theta }}_1}({\mathit{\tau }}_1;\mathit{\theta })}{1+Q_{10}({\mathit{\tau }}_1;\mathit{\theta })}& \frac{\overline{c}{(Q_{10})}_{{\mathit{\theta }}_2}({\mathit{\tau }}_1;\mathit{\theta })}{1+Q_{10}({\mathit{\tau }}_1;\mathit{\theta })}& \cdots & \frac{\overline{c}{(Q_{10})}_{{\mathit{\theta }}_6}({\mathit{\tau }}_1;\mathit{\theta })}{1+Q_{10}({\mathit{\tau }}_1;\mathit{\theta })}\\ \frac{\overline{c}{(Q_1)}_{{\mathit{\theta }}_1}({\mathit{\tau }}_2;\mathit{\theta })}{1+Q_1({\mathit{\tau }}_2;\mathit{\theta })}& \frac{\overline{c}{(Q_1)}_{{\mathit{\theta }}_2}({\mathit{\tau }}_2;\mathit{\theta })}{1+Q_1({\mathit{\tau }}_2;\mathit{\theta })}& \cdots & \frac{\overline{c}{(Q_1)}_{{\mathit{\theta }}_6}({\mathit{\tau }}_2;\mathit{\theta })}{1+Q_1({\mathit{\tau }}_2;\mathit{\theta })}\\ ⋮& ⋮& & ⋮\\ \frac{\overline{c}{(Q_{10})}_{{\mathit{\theta }}_1}({\mathit{\tau }}_2;\mathit{\theta })}{1+Q_{10}({\mathit{\tau }}_2;\mathit{\theta })}& \frac{\overline{c}{(Q_{10})}_{{\mathit{\theta }}_2}({\mathit{\tau }}_2;\mathit{\theta })}{1+Q_{10}({\mathit{\tau }}_2;\mathit{\theta })}& \cdots & \frac{\overline{c}{(Q_{10})}_{{\mathit{\theta }}_6}({\mathit{\tau }}_2;\mathit{\theta })}{1+Q_{10}({\mathit{\tau }}_2;\mathit{\theta })}\\ ⋮& ⋮& & ⋮\\ \frac{\overline{c}{(Q_1)}_{{\mathit{\theta }}_1}({\mathit{\tau }}_{11};\mathit{\theta })}{1+Q_1({\mathit{\tau }}_{11};\mathit{\theta })}& \frac{\overline{c}{(Q_1)}_{{\mathit{\theta }}_2}({\mathit{\tau }}_{11};\mathit{\theta })}{1+Q_1({\mathit{\tau }}_{11};\mathit{\theta })}& \cdots & \frac{\overline{c}{(Q_1)}_{{\mathit{\theta }}_6}({\mathit{\tau }}_{11};\mathit{\theta })}{1+Q_1({\mathit{\tau }}_{11};\mathit{\theta })}\\ ⋮& ⋮& & ⋮\\ \frac{\overline{c}{(Q_{10})}_{{\mathit{\theta }}_1}({\mathit{\tau }}_{11};\mathit{\theta })}{1+Q_{10}({\mathit{\tau }}_{11};\mathit{\theta })}& \frac{\overline{c}{(Q_{10})}_{{\mathit{\theta }}_2}({\mathit{\tau }}_{11};\mathit{\theta })}{1+Q_{10}({\mathit{\tau }}_{11};\mathit{\theta })}& \cdots & \frac{\overline{c}{(Q_{10})}_{{\mathit{\theta }}_6}({\mathit{\tau }}_{11};\mathit{\theta })}{1+Q_{10}({\mathit{\tau }}_{11};\mathit{\theta })}\end{array}\right],$$

where $\overline{c}=\sqrt{\mathit{\gamma }}/ln(10)$. To compute ${(Q_i)}_{{\mathit{\theta }}_p}$ and $y_{{\mathit{\theta }}_p}$, ($i=1,\dots ,10$ and $p=1,\dots ,6$) we used the following approximation$$\begin{array}{c}\hfill \begin{array}{c}{(Q_i)}_{{\mathit{\theta }}_p}({\mathit{\tau }}_r;\mathit{\theta })\approx \frac{Q{}_i^{\mathrm{\Delta }x,\mathrm{\Delta }t}({\mathit{\tau }}_r;{\mathit{\theta }}_1,\dots ,{\mathit{\theta }}_p+\mathrm{\Delta }{\mathit{\theta }}_p,\dots ,{\mathit{\theta }}_6)-Q{}_i^{\mathrm{\Delta }x,\mathrm{\Delta }t}({\mathit{\tau }}_r;{\mathit{\theta }}_1,\dots ,{\mathit{\theta }}_6)}{\mathrm{\Delta }{\mathit{\theta }}_p},\hfill \\ y_{{\mathit{\theta }}_p}\approx \frac{Y_{\mathrm{\Delta }t}({\mathit{\tau }}_r;{\mathit{\theta }}_1,\dots ,{\mathit{\theta }}_p+\mathrm{\Delta }{\mathit{\theta }}_p,\dots ,{\mathit{\theta }}_6)-Y_{\mathrm{\Delta }t}({\mathit{\tau }}_r;{\mathit{\theta }}_1,\dots ,{\mathit{\theta }}_6)}{\mathrm{\Delta }{\mathit{\theta }}_p},\hfill \end{array}\end{array}$$

for $r=1,\dots ,11$ days from the beginning of the experiment. The increment $\mathrm{\Delta }{\mathit{\theta }}_p=0.01{\mathit{\theta }}_p,p=1,2,\dots ,6,$ was used in the approximation of the sensitivity matrix Z. Additional simulations were taken to investigate the impact of the step size on the sensitivity calculations. The results suggest that reducing $\mathrm{\Delta }{\mathit{\theta }}_p$ beyond $0.01{\mathit{\theta }}_p$ did not significantly alter the sensitivity metrics. An alternative method for computing derivative approximations which addresses the ‘step size dilemma’ is discussed in [52].

Under the standard assumptions of classical nonlinear regression theory, we know that if the errors ${\mathit{ϵ}}_r$ satisfy ${\mathit{ϵ}}_r\sim \mathcal{N}(0,{\mathit{\sigma }}^2)$ where ${\mathit{ϵ}}_r$ is the difference between the model and the observations at time ${\mathit{\tau }}_r$, then the least-squares estimate ${\mathit{\theta }}^{\ast }$ is expected to be asymptotically normally distributed. In particular, for large samples, we may assume(7.2) $$\begin{array}{c}\hfill \begin{array}{c}{\mathit{\theta }}^{\ast }\sim \mathcal{N}[{\mathit{\theta }}^0,{\mathit{\sigma }}^2{\{Z^T({\mathit{\theta }}^0)Z({\mathit{\theta }}^0)\}}^{-1}],\hfill \end{array}\end{array}$$(7.2)

where ${\mathit{\theta }}^0$ is the true vector of parameters and ${\mathit{\sigma }}^2{\{Z^T({\mathit{\theta }}^0)Z({\mathit{\theta }}^0)\}}^{-1}$ is the true covariance matrix. Since ${\mathit{\theta }}^0$ and ${\mathit{\sigma }}^2$ are not available, we follow a standard statistical practice in which we substitute the computed estimate ${\mathit{\theta }}^{\ast }$ for ${\mathit{\theta }}^0$ and approximate the error variance ${\mathit{\sigma }}^2\approx {\widehat{\mathit{\sigma }}}^2$ in (7.2(7.2) $$\begin{array}{c}\hfill \begin{array}{c}{\mathit{\theta }}^{\ast }\sim \mathcal{N}[{\mathit{\theta }}^0,{\mathit{\sigma }}^2{\{Z^T({\mathit{\theta }}^0)Z({\mathit{\theta }}^0)\}}^{-1}],\hfill \end{array}\end{array}$$(7.2) ) using the residuals given by$$\begin{array}{c}\hfill \begin{array}{cc}{\mathit{\sigma }}^2=\hfill & \frac{1}{{11}^2-6}\left({\sum }_{r=1}^{11}|D_r^N-y({\mathit{\tau }}_r;{\mathit{\theta }}^{\ast }){|}^2+\mathit{\gamma }{\sum }_{r=1}^{11}{\sum }_{i=1}^{10}|{log}_{10}(D_{r,i}^Q+1)\right.\hfill \\ \hfill & \left.-{log}_{10}(Q{}_i({\mathit{\tau }}_r;{\mathit{\theta }}^{\ast })+1){|}^2\right),\hfill \end{array}\end{array}$$

which are computed using the finite difference approximation as follows:$$\begin{array}{c}\hfill \begin{array}{c}{\widehat{\mathit{\sigma }}}^2\approx \frac{1}{{11}^2-6}\left({\sum }_{r=1}^{11}|D_r^N-Y_{\mathrm{\Delta }t}({\mathit{\tau }}_r;{\mathit{\theta }}^{\ast }){|}^2+\mathit{\gamma }{\sum }_{r=1}^{11}{\sum }_{i=1}^{10}|{log}_{10}(D_{r,i}^Q+1)\right.\hfill \\ \left.-{log}_{10}(Q{}_i^{\mathrm{\Delta }x,\mathrm{\Delta }t}({\mathit{\tau }}_r;{\mathit{\theta }}^{\ast })+1){|}^2\right).\hfill \end{array}\end{array}$$

The covariance matrix is then approximated by$$\begin{array}{c}\hfill \begin{array}{c}E\approx {\widehat{\mathit{\sigma }}}^2{\{Z^T({\mathit{\theta }}^{\ast })Z({\mathit{\theta }}^{\ast })\}}^{-1}\hfill \end{array}\end{array}$$

and the standard errors are approximated by taking the diagonal elements of the covariance matrix of the approximations, i.e. $SE({\mathit{\theta }}_p)\approx \sqrt{E_{p,p}}$ for each parameter ${\mathit{\theta }}_p,p=1,2,\dots ,6$ respectively. The coefficient of variation (CV) is defined as $C{V}_p=SE({\mathit{\theta }}_p)/{\mathit{\theta }}_p$, respectively. The CV is also provided in Table 2 for an assessment of the degree of uncertainty in the estimates. The standard errors suggest that the greatest degree of uncertainty was found in the estimate for the encapsulation exponent, $P_a$.

Figure 4. Phytoplankton bloom formation: (a) measurements adapted from Figure 2 in [23] and (b) model results.

Figure 5. Phytoplankton bloom formation: (a) error levels given by ${log}_{10}(D_{r,i}^Q+1)-{log}_{10}(Q_i({\mathit{\tau }}_r;{\mathit{\theta }}^{\ast })+1)$ and (b) incremental total particle comparison.

Table 2. Least-squares parameter values and their standard errors.

Parameter	Interpretation	Estimated Value	Standard Error, $\widehat{\mathit{\sigma }}$	CV(%)
${\mathit{\theta }}_1=N_f$	Nutrient consumption constant	$6.1230\times {10}^{-2}$	$2.7638\times {10}^{-2}$	45
${\mathit{\theta }}_2={\mathit{\beta }}_f$	Reproduction factor	$3.1260\times {10}^2$	$2.6119\times {10}^2$	84
${\mathit{\theta }}_3={\mathit{\eta }}_f$	Sticking probability	$2.1555\times {10}^{-3}$	$2.3742\times {10}^{-3}$	110
${\mathit{\theta }}_4={\mathit{\kappa }}_0$	Encapsulation factor	$1.4316\times {10}^{-1}$	$2.2724\times {10}^{-1}$	159
${\mathit{\theta }}_5=P_a$	Encapsulation exponent	$4.4075\times {10}^0$	$4.2532\times {10}^1$	965
${\mathit{\theta }}_6=g_{\text{min}}$	Growth rate of smallest aggregate	$1.5942\times {10}^{-2}$	$2.2502\times {10}^{-2}$	141

The sensitivity plots given in Figure 6 suggest that the particle histograms (0.5 ${\text{mm}}^3$ size class) were most sensitive to the growth rate of the smallest aggregate ${\mathit{\theta }}_6=g_{\text{min}}$, the sticking probability ${\mathit{\theta }}_3={\mathit{\eta }}_f$, and the aggregate encapsulation factor ${\mathit{\theta }}_4={\mathit{\kappa }}_0$ respectively. This outcome underscores the importance of the aggregate growth process in modeling phytoplankton bloom and nutrient dynamics in general. The histogram was relatively insensitive to perturbations in the reproduction factor ${\mathit{\beta }}_f$ or the encapsulation exponent $P_a$. The nutrient dynamics were sensitive to the aggregate growth parameters as well as the nutrient consumption parameter. The relative sensitivities of the nutrient dynamics with respect to all parameters approached zero after the nutrients fell below detection limits in the tank (day 9). We point out that the sensitivity of the histogram of mid-sized particles with respect to changes in the nutrient consumption rate ${\mathit{\theta }}_1=N_f$ increased from virtually zero after nutrients fell below the detection limit. This suggests a possible time delay between the nutrient dynamics and bloom formation which is physically reasonable given the lag which occurs between nutrient consumption, storage, and cell growth. The overall conclusion of the sensitivity analysis is that in addition to coagulation, the phytoplankton-nutrient dynamics were most sensitive to the growth and nutrient consumption processes, in particular the growth rate of small aggregates and the role of encapsulation as reflected in ${\mathit{\kappa }}_0$. These findings point to the significance of size-specific processes which should be investigated in experimental studies as well as the potential role of fragmentation in improving the bloom peak predictions in the future.

Figure 6. Relative sensitivity of model output to parameter changes: (a) phytoplankton sensitivity ($0.5{\text{mm}}^3$ particles) and (b) nutrient sensitivity.

7.1. Influence of physical processes on bloom formation

As noted earlier, the particle dynamics are driven primarily by physical processes (i.e. coagulation and fragmentation) once nutrients fall below the detection limit. Since the model (2.1(2.1) $$\begin{array}{c}\hfill \begin{array}{c}\frac{\mathit{\partial }u}{\mathit{\partial }t}+\frac{\mathit{\partial }}{\mathit{\partial }x}\left(g(x,t,\mathit{\varphi },\overrightarrow{y})u\right)+R(x,t,\mathit{\varphi },\overrightarrow{y})u=A(u),(x,t)\in (0,x_{max}]\times (0,T],\hfill \\ \frac{\mathrm{d}{y}_{\ell }}{\mathrm{d}t}=f_{\ell }(t,{\mathit{\sigma }}_{\ell },\overrightarrow{y})-h_{\ell }(t,{\mathit{\sigma }}_{\ell },\overrightarrow{y})y_{\ell },t\in (0,T],\ell =1,\dots ,m,\hfill \\ g(0,t,\mathit{\varphi },\overrightarrow{y})u(0,t)=F(t,\overrightarrow{y})+{\int }_0^{x_{\text{max}}}\mathit{\beta }(x,t,\mathit{\varphi },\overrightarrow{y})u(x,t)\mathrm{d}x,t\in (0,T],\hfill \\ u(x,0)=u^0(x),x\in [0,x_{max}],y_{\ell }(0)=y_{\ell ,0}\ge 0,\ell =1,\dots ,m,\hfill \end{array}\end{array}$$(2.1) ) does not include fragmentation, a least squares simulation was performed to investigate the impact of up-scaling the coagulation function by increasing the coagulation rate of small to intermediate sized particles. This optimization run used the least squares parameters provided in Table 2 to initialize the parameter vector. The coagulation (Figure 7) scale-up factor is given by$$\begin{array}{c}\hfill \begin{array}{c}\mathrm{\Pi }(x,s)={\left(2-\frac{x+s}{x_{max}}\right)}^{P_{\mathit{\eta }}},\hfill \end{array}\end{array}$$

where the exponent $P_{\mathit{\eta }}$ is an estimated parameter which controls the amplification rate. Larger values of $P_{\mathit{\eta }}$ correspond to an accelerated coagulation rate of small to intermediate sized particles. An initial value of $P_{\mathit{\eta }}=1$ was taken to complete the least squares initial parameter vector. The shape factor $\mathrm{\Pi }$ is scaled against the coagulation function to obtain the accelerated coagulation function. The assumption of this trial was to confirm whether altering the coagulation function in this manner could compensate for omitting the fragmentation process (see Figure 9).

Figure 7. Accelerated coagulation rates, $\widehat{\mathit{\eta }}(x,s)=\mathrm{\Pi }(x,s)\mathit{\eta }(x,s)$, $P_{\mathit{\eta }}=5.88$.

Figure 8. Phytoplankton bloom formation: (a) measurements adapted from Figure 2 in [23] and (b) model results with accelerated coagulation.

Figure 9. Nutrient model performance with accelerated coagulation.

Figure 10. Phytoplankton bloom formation with accelerated coagulation: (a) error levels given by ${log}_{10}(D_{r,i}^Q+1)-{log}_{10}(Q_i({\mathit{\tau }}_r;{\mathit{\theta }}^{\ast })+1)$ and (b) incremental total particle comparison.

The results (Figures 8 and 9) suggest that increasing the coagulation rate of small particles slightly improves the late-stage performance of the particle histogram model for particles between 0.4 ${\text{mm}}^3$ and 1.0 ${\text{mm}}^3$ without compromising the performance of the nutrient model (Figure 9). However, accelerated coagulation rates increase the reproduction factor in this model by an order of magnitude (${\mathit{\beta }}_f=1.97\times {10}^3$) which leads to a uniform overestimation of the total particles (Figure 10). These results further suggest that the reproduction process and hence the late stage particle dynamics could be modeled more accurately in the future by including the fragmentation process. Ample research exists investigating the role of the fragmentation process on size-structured aggregating particle systems, and parameter estimation within an inverse problem framework. Moreover, the techniques developed and insights presented by the researchers in [53–59] can be employed to include a fragmentation process explicitly which should improve the phytoplankton model predictions presented here.

8. Conclusions

This paper presents a new finite difference scheme (4.1(4.1) $$\begin{array}{cc}& \frac{u_j^{k+1}-u_j^k}{\mathrm{\Delta }t}+\frac{g_j^k{u}_j^{k+1}-g_{j-1}^k{u}_{j-1}^{k+1}}{\mathrm{\Delta }x}+u_j^{k+1}\sum _{m=1}^N{\mathit{\eta }}_{m,j}^k{u}_m^k\mathrm{\Delta }x+R_j^k{u}_j^{k+1}\hfill \\ \hfill & =\frac{1}{2}\sum _{s=1}^{j-1}{\mathit{\eta }}_{j-s,s}^k{u}_{j-s}^k{u}_s^{k+1}\mathrm{\Delta }x,\hfill \\ \hfill & g_0^k{u}_0^{k+1}=F^k+\sum _{j=1}^N{\mathit{\beta }}_j^k{u}_j^k\mathrm{\Delta }x,u_j^0=u_0(x_j),\hfill \\ \hfill & \frac{y_{\ell }^{k+1}-y_{\ell }^k}{\mathrm{\Delta }t}=f_{\ell }^k-h_{\ell }^k{y}_{\ell }^{k+1},y_{\ell }^0=y_{\ell ,0},\ell =1,\dots ,m.\hfill \end{array}$$(4.1) ) for the size-structured particle aggregation problem (2.1(2.1) $$\begin{array}{c}\hfill \begin{array}{c}\frac{\mathit{\partial }u}{\mathit{\partial }t}+\frac{\mathit{\partial }}{\mathit{\partial }x}\left(g(x,t,\mathit{\varphi },\overrightarrow{y})u\right)+R(x,t,\mathit{\varphi },\overrightarrow{y})u=A(u),(x,t)\in (0,x_{max}]\times (0,T],\hfill \\ \frac{\mathrm{d}{y}_{\ell }}{\mathrm{d}t}=f_{\ell }(t,{\mathit{\sigma }}_{\ell },\overrightarrow{y})-h_{\ell }(t,{\mathit{\sigma }}_{\ell },\overrightarrow{y})y_{\ell },t\in (0,T],\ell =1,\dots ,m,\hfill \\ g(0,t,\mathit{\varphi },\overrightarrow{y})u(0,t)=F(t,\overrightarrow{y})+{\int }_0^{x_{\text{max}}}\mathit{\beta }(x,t,\mathit{\varphi },\overrightarrow{y})u(x,t)\mathrm{d}x,t\in (0,T],\hfill \\ u(x,0)=u^0(x),x\in [0,x_{max}],y_{\ell }(0)=y_{\ell ,0}\ge 0,\ell =1,\dots ,m,\hfill \end{array}\end{array}$$(2.1) ) in the context of phytoplankton dynamics. The model includes growth and coagulation terms and a general system of differential equations which represent the environment. We prove that the numerical scheme (4.1(4.1) $$\begin{array}{cc}& \frac{u_j^{k+1}-u_j^k}{\mathrm{\Delta }t}+\frac{g_j^k{u}_j^{k+1}-g_{j-1}^k{u}_{j-1}^{k+1}}{\mathrm{\Delta }x}+u_j^{k+1}\sum _{m=1}^N{\mathit{\eta }}_{m,j}^k{u}_m^k\mathrm{\Delta }x+R_j^k{u}_j^{k+1}\hfill \\ \hfill & =\frac{1}{2}\sum _{s=1}^{j-1}{\mathit{\eta }}_{j-s,s}^k{u}_{j-s}^k{u}_s^{k+1}\mathrm{\Delta }x,\hfill \\ \hfill & g_0^k{u}_0^{k+1}=F^k+\sum _{j=1}^N{\mathit{\beta }}_j^k{u}_j^k\mathrm{\Delta }x,u_j^0=u_0(x_j),\hfill \\ \hfill & \frac{y_{\ell }^{k+1}-y_{\ell }^k}{\mathrm{\Delta }t}=f_{\ell }^k-h_{\ell }^k{y}_{\ell }^{k+1},y_{\ell }^0=y_{\ell ,0},\ell =1,\dots ,m.\hfill \end{array}$$(4.1) ) converges to the unique weak solution of the model (2.1(2.1) $$\begin{array}{c}\hfill \begin{array}{c}\frac{\mathit{\partial }u}{\mathit{\partial }t}+\frac{\mathit{\partial }}{\mathit{\partial }x}\left(g(x,t,\mathit{\varphi },\overrightarrow{y})u\right)+R(x,t,\mathit{\varphi },\overrightarrow{y})u=A(u),(x,t)\in (0,x_{max}]\times (0,T],\hfill \\ \frac{\mathrm{d}{y}_{\ell }}{\mathrm{d}t}=f_{\ell }(t,{\mathit{\sigma }}_{\ell },\overrightarrow{y})-h_{\ell }(t,{\mathit{\sigma }}_{\ell },\overrightarrow{y})y_{\ell },t\in (0,T],\ell =1,\dots ,m,\hfill \\ g(0,t,\mathit{\varphi },\overrightarrow{y})u(0,t)=F(t,\overrightarrow{y})+{\int }_0^{x_{\text{max}}}\mathit{\beta }(x,t,\mathit{\varphi },\overrightarrow{y})u(x,t)\mathrm{d}x,t\in (0,T],\hfill \\ u(x,0)=u^0(x),x\in [0,x_{max}],y_{\ell }(0)=y_{\ell ,0}\ge 0,\ell =1,\dots ,m,\hfill \end{array}\end{array}$$(2.1) ). We then use the model to investigate the role of the aggregate growth processes on the formation of an algal bloom based on digitized measurements collected in the 1993 SIGMA reactor experiment [23,49]. The model simulations suggest that the aggregate growth rate processes (i.e. reproduction, encapsulation, growth rate of smallest aggregates) have an appreciable effect on the timing of nutrient depletion and bloom peak predictions.

The determination of a biologically sound mathematical expression for cell aggregate growth rates is confounded by the general assumption that only a portion of the aggregate is alive (i.e. active fraction) and capable of growth and cell division at any given time [25]. We also note that the mathematical model of growth presented here and in many applications is based on the simplification that the aggregate particles can be represented using the concept of equivalent spherical diameter (ESD) determined by a volume balance. Real flocculated particles are irregular and more accurately characterized using fractal geometry techniques [60]. The fractal properties of aggregates could affect physical (e.g. settling velocity) and biological characteristics (e.g. active volume, reproduction). These challenges are magnified when coupled with the notion of a variable internal storage pool which controls nutrient depletion rates. Moreover, nutrient uptake by aggregates depends not only on the outward size of the particle but more directly on the internal storage pool which allows growth and cell reproduction to continue when nutrients fall below the limits of chemical detection [61]. The implication here is that in order to improve the predictions of algae bloom development and nutrient uptake, in addition to the environment and floc characteristics, there are at least three fundamental mechanisms requiring further investigation. We believe that the first mechanism is distributed reproduction, i.e. fragmentation which could serve to distribute particles more evenly across the size classes during the bloom peak. The second is the nutrient adsorption rate (e.g. as modeled by the Langmuir isotherm) which serves to pull nutrients into cell storage, and should depend on the active volume of the aggregate. The next (and perhaps more difficult) mechanism is the internal storage pool of cell aggregates and the rate at which nutrients are incorporated into the construction of new cell material (i.e. cell growth) [62]. The rigorous development and application of an internal pool model (e.g. [63]) including the adsorption and incorporation mechanisms in the more general context of cell aggregates should improve the performance of the model in biogeochemical applications.

Additional information

Funding

This work is supported in part by the National Science Foundation [grant number DMS-1312963].

Notes

No potential conflict of interest was reported by the authors.

Appendix 1

Proofs of lemmas and theorems

Proof of Lemma 4.1:

First observe that for any real numbers $a_{j,s}$ and positive integer N we have the interchange of summation ${\mathrm{\Sigma }}_{j=1}^N{\mathrm{\Sigma }}_{s=1}^j{a}_{j,s}={\mathrm{\Sigma }}_{s=1}^N{\mathrm{\Sigma }}_{j=s}^N{a}_{j,s},$ so by direct calculation and assumption (B10),$$\begin{array}{c}\hfill \begin{array}{cc}{\sum }_{j=1}^N{\sum }_{s=1}^{j-1}{\mathit{\eta }}_{j-s,s}^k{u}_{j-s}^k{u}_s^{k+1}\hfill & ={\sum }_{s=1}^N{\sum }_{j=s+1}^N{\mathit{\eta }}_{j-s,j}^k{u}_{j-s}^k{u}_s^{k+1}={\sum }_{s=1}^N{\sum }_{m=1}^{N-s}{\mathit{\eta }}_{m,s}^k{u}_m^k{u}_s^{k+1}\hfill \\ \hfill & ={\sum }_{j=1}^N{u}_j^{k+1}{\sum }_{m=1}^{N-j}{\mathit{\eta }}_{m,j}^k{u}_m^k={\sum }_{j=1}^N{u}_j^{k+1}{\sum }_{m=1}^N{\mathit{\eta }}_{m,j}^k{u}_m^k.\hfill \end{array}\end{array}$$

where the last step follows from assumption (B10), namely ${\mathit{\eta }}_{m,j}^k=0$ if $m+j>N$. Therefore,(A1) $$\begin{array}{c}\hfill \begin{array}{cc}& (\frac{1}{2}{\sum }_{j=1}^N{\sum }_{s=1}^{j-1}{\mathit{\eta }}_{j-s,s}^k{u}_{j-s}^k{u}_s^{k+1}-{\sum }_{j=1}^N{u}_j^{k+1}{\sum }_{m=1}^N{\mathit{\eta }}_{m,j}^k{u}_m^k)\mathrm{\Delta }x\mathrm{\Delta }x\hfill \\ \hfill & =(\frac{1}{2}{\sum }_{j=1}^N{u}_j^{k+1}{\sum }_{m=1}^N{\mathit{\eta }}_{m,j}^k{u}_m^k-{\sum }_{j=1}^N{u}_j^{k+1}{\sum }_{m=1}^N{\mathit{\eta }}_{m,j}^k{u}_m^k)\mathrm{\Delta }x\mathrm{\Delta }x\hfill \\ \hfill & =-\frac{1}{2}{\sum }_{j=1}^N{u}_j^{k+1}{\sum }_{m=1}^N{\mathit{\eta }}_{m,j}^k{u}_m^k\mathrm{\Delta }x\mathrm{\Delta }x.\hfill \end{array}\end{array}$$(A1)

This completes the proof.$\square$

Proof of Lemma 4.2:

First multiply the first row of (4.1(4.1) $$\begin{array}{cc}& \frac{u_j^{k+1}-u_j^k}{\mathrm{\Delta }t}+\frac{g_j^k{u}_j^{k+1}-g_{j-1}^k{u}_{j-1}^{k+1}}{\mathrm{\Delta }x}+u_j^{k+1}\sum _{m=1}^N{\mathit{\eta }}_{m,j}^k{u}_m^k\mathrm{\Delta }x+R_j^k{u}_j^{k+1}\hfill \\ \hfill & =\frac{1}{2}\sum _{s=1}^{j-1}{\mathit{\eta }}_{j-s,s}^k{u}_{j-s}^k{u}_s^{k+1}\mathrm{\Delta }x,\hfill \\ \hfill & g_0^k{u}_0^{k+1}=F^k+\sum _{j=1}^N{\mathit{\beta }}_j^k{u}_j^k\mathrm{\Delta }x,u_j^0=u_0(x_j),\hfill \\ \hfill & \frac{y_{\ell }^{k+1}-y_{\ell }^k}{\mathrm{\Delta }t}=f_{\ell }^k-h_{\ell }^k{y}_{\ell }^{k+1},y_{\ell }^0=y_{\ell ,0},\ell =1,\dots ,m.\hfill \end{array}$$(4.1) ) by $\mathrm{\Delta }t\mathrm{\Delta }x$ and sum from $j=1$ to N to get(A2) $$\begin{array}{c}\hfill \begin{array}{cc}& {\sum }_{j=1}^N(1+\mathrm{\Delta }t{R}_j^k)u_j^{k+1}\mathrm{\Delta }x\hfill \\ \hfill & ={\sum }_{j=1}^N{u}_j^k\mathrm{\Delta }x-\mathrm{\Delta }t{\sum }_{j=1}^N(g_j^k{u}_j^{k+1}-g_{j-1}^k{u}_{j-1}^{k+1})\hfill \\ \hfill & +\frac{\mathrm{\Delta }t}{2}{\sum }_{j=1}^N{\sum }_{s=1}^{j-1}{\mathit{\eta }}_{j-s,s}^k{u}_{j-s}^k{u}_s^{k+1}\mathrm{\Delta }x\mathrm{\Delta }x-\mathrm{\Delta }t{\sum }_{j=1}^N{u}_j^{k+1}{\sum }_{m=1}^N{\mathit{\eta }}_{m,j}^k{u}_m^k\mathrm{\Delta }x\mathrm{\Delta }x\hfill \\ \hfill & ={\sum }_{j=1}^N{u}_j^k\mathrm{\Delta }x+\mathrm{\Delta }t(F^k+{\sum }_{j=1}^N{\mathit{\beta }}_j^k{u}_j^k\mathrm{\Delta }x)-\frac{\mathrm{\Delta }t}{2}{\sum }_{j=1}^N{u}_j^{k+1}{\sum }_{m=1}^N{\mathit{\eta }}_{m,j}^k{u}_m^k\mathrm{\Delta }x\mathrm{\Delta }x\hfill \\ \hfill & \le {\sum }_{j=1}^N{u}_j^k\mathrm{\Delta }x+\mathrm{\Delta }t(F^k+{\sum }_{j=1}^N{\mathit{\beta }}_j^k{u}_j^k\mathrm{\Delta }x),\hfill \end{array}\end{array}$$(A2)

which follows from assumption (B1), Lemma 4.1, and nonnegativity of the numerical scheme, respectively. Since the scheme is nonnegative, we have ${\mathrm{\Sigma }}_{j=1}^N{u}_j^{k+1}\mathrm{\Delta }x={\Vert u^{k+1}\Vert }_1$ and hence by (A1(A1) $$\begin{array}{c}\hfill \begin{array}{cc}& (\frac{1}{2}{\sum }_{j=1}^N{\sum }_{s=1}^{j-1}{\mathit{\eta }}_{j-s,s}^k{u}_{j-s}^k{u}_s^{k+1}-{\sum }_{j=1}^N{u}_j^{k+1}{\sum }_{m=1}^N{\mathit{\eta }}_{m,j}^k{u}_m^k)\mathrm{\Delta }x\mathrm{\Delta }x\hfill \\ \hfill & =(\frac{1}{2}{\sum }_{j=1}^N{u}_j^{k+1}{\sum }_{m=1}^N{\mathit{\eta }}_{m,j}^k{u}_m^k-{\sum }_{j=1}^N{u}_j^{k+1}{\sum }_{m=1}^N{\mathit{\eta }}_{m,j}^k{u}_m^k)\mathrm{\Delta }x\mathrm{\Delta }x\hfill \\ \hfill & =-\frac{1}{2}{\sum }_{j=1}^N{u}_j^{k+1}{\sum }_{m=1}^N{\mathit{\eta }}_{m,j}^k{u}_m^k\mathrm{\Delta }x\mathrm{\Delta }x.\hfill \end{array}\end{array}$$(A1) ), noting that $R_j^k\ge 0$ the relationship (A2(A2) $$\begin{array}{c}\hfill \begin{array}{cc}& {\sum }_{j=1}^N(1+\mathrm{\Delta }t{R}_j^k)u_j^{k+1}\mathrm{\Delta }x\hfill \\ \hfill & ={\sum }_{j=1}^N{u}_j^k\mathrm{\Delta }x-\mathrm{\Delta }t{\sum }_{j=1}^N(g_j^k{u}_j^{k+1}-g_{j-1}^k{u}_{j-1}^{k+1})\hfill \\ \hfill & +\frac{\mathrm{\Delta }t}{2}{\sum }_{j=1}^N{\sum }_{s=1}^{j-1}{\mathit{\eta }}_{j-s,s}^k{u}_{j-s}^k{u}_s^{k+1}\mathrm{\Delta }x\mathrm{\Delta }x-\mathrm{\Delta }t{\sum }_{j=1}^N{u}_j^{k+1}{\sum }_{m=1}^N{\mathit{\eta }}_{m,j}^k{u}_m^k\mathrm{\Delta }x\mathrm{\Delta }x\hfill \\ \hfill & ={\sum }_{j=1}^N{u}_j^k\mathrm{\Delta }x+\mathrm{\Delta }t(F^k+{\sum }_{j=1}^N{\mathit{\beta }}_j^k{u}_j^k\mathrm{\Delta }x)-\frac{\mathrm{\Delta }t}{2}{\sum }_{j=1}^N{u}_j^{k+1}{\sum }_{m=1}^N{\mathit{\eta }}_{m,j}^k{u}_m^k\mathrm{\Delta }x\mathrm{\Delta }x\hfill \\ \hfill & \le {\sum }_{j=1}^N{u}_j^k\mathrm{\Delta }x+\mathrm{\Delta }t(F^k+{\sum }_{j=1}^N{\mathit{\beta }}_j^k{u}_j^k\mathrm{\Delta }x),\hfill \end{array}\end{array}$$(A2) ) satisfies$$\begin{array}{c}\hfill \begin{array}{cc}\Vert u^{k+1}{\Vert }_1\hfill & \le \Vert u^k{\Vert }_1+\mathrm{\Delta }t(c_1+c_2|y^k{|}_1+{\Vert \mathit{\beta }\Vert }_{\infty }\Vert u^k{\Vert }_1)\hfill \\ \hfill & ={(1+\Vert \mathit{\beta }\Vert }_{\infty }\mathrm{\Delta }t)\Vert u^k{\Vert }_1+(c_1+c_2|y^k{|}_1)\mathrm{\Delta }t.\hfill \end{array}\end{array}$$

Next observe that by nonnegativity of the scheme and (B3) and (B4),$$\begin{array}{c}\hfill \begin{array}{c}{y}_{\ell }^{k+1}=\frac{y_{\ell }^k+f_{\ell }^k\mathrm{\Delta }t}{1+h_{\ell }^k\mathrm{\Delta }t}\le (1+c_3\mathrm{\Delta }t)y_{\ell }^k+(1+c_2\mathrm{\Delta }t)|{\mathit{\sigma }}_{\ell }^k|+c_1\mathrm{\Delta }t.\hfill \end{array}\end{array}$$

By nonnegativity of ${\overrightarrow{y}}^k$ we have ${\mathrm{\Sigma }}_{\ell =1}^m{y}_{\ell }^k={|y^k|}_1$ and since $|{\mathit{\sigma }}_{\ell }^k{|\le \Vert \mathit{\gamma }\Vert }_{\infty }{\Vert u^k\Vert }_1$, we get$$\begin{array}{c}\hfill \begin{array}{c}|y^{k+1}{|}_1\le (1+c_3\mathrm{\Delta }t)|y^k{|}_1+m{c}_2{\Vert \mathit{\gamma }\Vert }_{\infty }{\Vert u^k\Vert }_1\mathrm{\Delta }t+m{c}_1\mathrm{\Delta }t.\hfill \end{array}\end{array}$$

Thus,$$\begin{array}{c}\hfill \begin{array}{cc}\Vert u^{k+1}{\Vert }_1+{|y^{k+1}|}_1\hfill & \le {(1+(\Vert \mathit{\beta }\Vert }_{\infty }+c_2{\Vert \mathit{\gamma }\Vert }_{\infty })\mathrm{\Delta }t)\Vert u^k{\Vert }_1+(1+c_3\mathrm{\Delta }t){|y^k|}_1+2c_1\mathrm{\Delta }t.\hfill \end{array}\end{array}$$

Let $C_{1,1}={\Vert \mathit{\beta }\Vert }_{\infty }+c_2{\Vert \mathit{\gamma }\Vert }_{\infty }$, $C_{1,2}=2c_1$, and $C_{1,3}=max\{C_{1,1},c_3\}$.(A3) $$\begin{array}{c}\hfill \begin{array}{cc}\Vert u^{k+1}{\Vert }_1+{|y^{k+1}|}_1\hfill & \le (1+C_{1,3}\mathrm{\Delta }t)(\Vert u^k{\Vert }_1+|y^k{|}_1)+C_{1,2}\mathrm{\Delta }t.\hfill \end{array}\end{array}$$(A3)

The result now follows from inductive application of (A3(A3) $$\begin{array}{c}\hfill \begin{array}{cc}\Vert u^{k+1}{\Vert }_1+{|y^{k+1}|}_1\hfill & \le (1+C_{1,3}\mathrm{\Delta }t)(\Vert u^k{\Vert }_1+|y^k{|}_1)+C_{1,2}\mathrm{\Delta }t.\hfill \end{array}\end{array}$$(A3) ) and elementary calculations.$\square$

Proof of Lemma 4.3:

Define$$\begin{array}{c}\hfill \begin{array}{c}{u}_{j_0}^{k+1}={max}_{i=0,1,\dots ,N}\left(u_i^{k+1}\right)={\Vert u^{k+1}\Vert }_{\infty }\hfill \end{array}\end{array}$$

If $j_0=0$, substituting the boundary conditions at the corner, and appealing to assumptions (B6), (B7), and Lemma 4.2,$$\begin{array}{c}\hfill \begin{array}{c}\Vert u^{k+1}{\Vert }_{\infty }\le \frac{c_2+(c_3+{\Vert \mathit{\beta }\Vert }_{\infty })B_2}{\mathit{\zeta }}.\hfill \end{array}\end{array}$$

Otherwise, fix $j_0$, and multiply the numerical scheme by $\mathrm{\Delta }t$,(A4) $$\begin{array}{c}\hfill \begin{array}{cc}(1+\mathrm{\Delta }t{R}_{j_0}^k)u_{j_0}^{k+1}\hfill & =u_{j_0}^k-\frac{\mathrm{\Delta }t}{\mathrm{\Delta }x}((g_{j_0}^k-g_{j_0-1}^k)u_{j_0}^{k+1}+(u_{j_0}^{k+1}-u_{j_0-1}^{k+1})g_{j_0-1}^k)\hfill \\ \hfill & +\frac{\mathrm{\Delta }t}{2}{\sum }_{s=1}^{j_0-1}{\mathit{\eta }}_{j_0-s,s}^k{u}_{j_0-s}^k{u}_s^{k+1}\mathrm{\Delta }x-\mathrm{\Delta }t{u}_{j_0}^{k+1}{\sum }_{m=1}^N{\mathit{\eta }}_{j_0,m}^k{u}_{j_0}^k\mathrm{\Delta }x.\hfill \end{array}\end{array}$$(A4)

Using $u_{j_0}^{k+1}\ge u_{j_0-1}^{k+1}$, nonnegativity of the scheme, and $g_{j_0-1}^k\ge 0$ we have$$\begin{array}{c}\hfill \begin{array}{cc}\left(1+\mathrm{\Delta }t(R_{j_0}^k+\frac{g_{j_0}^k-g_{j_0-1}^k}{\mathrm{\Delta }x})\right)u_{j_0}^{k+1}\hfill & \le u_{j_0}^k+\frac{\mathrm{\Delta }t}{2}\Vert u^k{\Vert }_{\infty }{\Vert \mathit{\eta }\Vert }_{\infty }{\Vert u^{k+1}\Vert }_1\hfill \\ \hfill & \le (1+\frac{\mathrm{\Delta }t}{2}{\Vert \mathit{\eta }\Vert }_{\infty }\Vert u^{k+1}{\Vert }_1)\Vert u^k{\Vert }_{\infty }.\hfill \end{array}\end{array}$$

Let $C_{\infty }=B_2{\Vert \mathit{\eta }\Vert }_{\infty }/2$, and by assumption (B5) we arrive at$$\begin{array}{c}\hfill \begin{array}{c}(1-{\mathit{\omega }}_1\mathrm{\Delta }t)\Vert u^{k+1}{\Vert }_{\infty }\le (1+C_{\infty }\mathrm{\Delta }t){\Vert u^k\Vert }_{\infty }.\hfill \end{array}\end{array}$$

By induction on k and assumption (B9) we establish the upper bound $B_3$.$\square$

Proof of Lemma 4.4:

We omit the details for the constant $B_{4,1}$, since the results follow directly from Lemma 4.2 and the arguments presented in Lemma 5.3 in [48]. Observe that by (4.1(4.1) $$\begin{array}{cc}& \frac{u_j^{k+1}-u_j^k}{\mathrm{\Delta }t}+\frac{g_j^k{u}_j^{k+1}-g_{j-1}^k{u}_{j-1}^{k+1}}{\mathrm{\Delta }x}+u_j^{k+1}\sum _{m=1}^N{\mathit{\eta }}_{m,j}^k{u}_m^k\mathrm{\Delta }x+R_j^k{u}_j^{k+1}\hfill \\ \hfill & =\frac{1}{2}\sum _{s=1}^{j-1}{\mathit{\eta }}_{j-s,s}^k{u}_{j-s}^k{u}_s^{k+1}\mathrm{\Delta }x,\hfill \\ \hfill & g_0^k{u}_0^{k+1}=F^k+\sum _{j=1}^N{\mathit{\beta }}_j^k{u}_j^k\mathrm{\Delta }x,u_j^0=u_0(x_j),\hfill \\ \hfill & \frac{y_{\ell }^{k+1}-y_{\ell }^k}{\mathrm{\Delta }t}=f_{\ell }^k-h_{\ell }^k{y}_{\ell }^{k+1},y_{\ell }^0=y_{\ell ,0},\ell =1,\dots ,m.\hfill \end{array}$$(4.1) ),$$\begin{array}{c}\hfill \begin{array}{c}{u}_1^{k+1}-u_0^{k+1}=(u_1^{k+1}-u_1^k)+(u_1^k-u_0^k)+(u_0^k-u_0^{k+1})\hfill \end{array}\end{array}$$

and,$$\begin{array}{c}\hfill \begin{array}{cc}{u}_1^{k+1}-u_1^k\hfill & =-\mathrm{\Delta }t{u}_1^{k+1}{\sum }_{m=1}^N{\mathit{\eta }}_{m,1}^k{u}_m^k\mathrm{\Delta }x-\frac{\mathrm{\Delta }t}{\mathrm{\Delta }x}((g_1^k-g_0^k)u_1^{k+1}+(u_1^{k+1}-u_0^{k+1})g_0^k)\hfill \\ \hfill & -\mathrm{\Delta }t{R}_1^k{u}_1^{k+1},\hfill \end{array}\end{array}$$

and hence by assumption (B2) and Lemmas 4.2 and 4.3,$$\begin{array}{c}\hfill \begin{array}{cc}(1+\frac{\mathrm{\Delta }t}{\mathrm{\Delta }x}{g}_0^k)|u_1^{k+1}-u_0^{k+1}|\hfill & \le \mathrm{\Delta }t\Vert u^{k+1}{\Vert }_{\infty }{(\Vert \mathit{\eta }\Vert }_{\infty }\Vert u^k{\Vert }_1+B_{4,1}+{\Vert R\Vert }_{\infty })\hfill \\ \hfill & +|u_1^k-u_0^k|+|u_0^k-u_0^{k+1}|\hfill \\ \hfill & \le |u_1^k-u_0^k|+B_3(c_0{B}_2+B_{4,1}+B_{1,1})\mathrm{\Delta }t+|u_0^k-u_0^{k+1}|\hfill \end{array}\end{array}$$

Next observe that by the boundary condition in (4.1(4.1) $$\begin{array}{cc}& \frac{u_j^{k+1}-u_j^k}{\mathrm{\Delta }t}+\frac{g_j^k{u}_j^{k+1}-g_{j-1}^k{u}_{j-1}^{k+1}}{\mathrm{\Delta }x}+u_j^{k+1}\sum _{m=1}^N{\mathit{\eta }}_{m,j}^k{u}_m^k\mathrm{\Delta }x+R_j^k{u}_j^{k+1}\hfill \\ \hfill & =\frac{1}{2}\sum _{s=1}^{j-1}{\mathit{\eta }}_{j-s,s}^k{u}_{j-s}^k{u}_s^{k+1}\mathrm{\Delta }x,\hfill \\ \hfill & g_0^k{u}_0^{k+1}=F^k+\sum _{j=1}^N{\mathit{\beta }}_j^k{u}_j^k\mathrm{\Delta }x,u_j^0=u_0(x_j),\hfill \\ \hfill & \frac{y_{\ell }^{k+1}-y_{\ell }^k}{\mathrm{\Delta }t}=f_{\ell }^k-h_{\ell }^k{y}_{\ell }^{k+1},y_{\ell }^0=y_{\ell ,0},\ell =1,\dots ,m.\hfill \end{array}$$(4.1) ), we have(A5) $$\begin{array}{c}\hfill \begin{array}{cc}|u_0^{k+1}-u_0^k|\hfill & =\frac{1}{g_0^{k-1}}|F^k-F^{k-1}|+\frac{1}{g_0^{k-1}}{\sum }_{j=1}^N|{\mathit{\beta }}_j^k-{\mathit{\beta }}_j^{k-1}|u_j^k\mathrm{\Delta }x\hfill \\ \hfill & +\frac{1}{g_0^{k-1}}|{\sum }_{j=1}^N{\mathit{\beta }}_j^{k-1}(u_j^k-u_j^{k-1})\mathrm{\Delta }x|+\frac{1}{g_0^{k-1}}|g_0^k-g_0^{k-1}|u_0^{k+1}.\hfill \end{array}\end{array}$$(A5)

Next observe that$$\begin{array}{c}\hfill \begin{array}{cc}& {\sum }_{j=1}^N{\mathit{\beta }}_j^{k-1}(u_j^k-u_j^{k-1})\mathrm{\Delta }x=-\frac{\mathrm{\Delta }t}{\mathrm{\Delta }x}{\sum }_{j=1}^N{\mathit{\beta }}_j^{k-1}((g_j^{k-1}-g_{j-1}^{k-1})u_{j-1}^k+(u_j^k-u_{j-1}^k)g_{j-1}^{k-1})\mathrm{\Delta }x\hfill \\ \hfill & -\mathrm{\Delta }t{\sum }_{j=1}^N{\mathit{\beta }}_j^{k-1}{u}_j^k{\sum }_{m=1}^N{\mathit{\eta }}_{m,j}^{k-1}{u}_m^{k-1}\mathrm{\Delta }x\mathrm{\Delta }x+\frac{\mathrm{\Delta }t}{2}{\sum }_{j=1}^N{\mathit{\beta }}_j^{k-1}{\sum }_{s=1}^{j-1}{\mathit{\eta }}_{j-s,s}^{k-1}{u}_{j-s}^k{u}_s^{k-1}\mathrm{\Delta }x\mathrm{\Delta }x\hfill \\ \hfill & -\mathrm{\Delta }t{\sum }_{j=1}^N{\mathit{\beta }}_j^{k-1}{R}_j^{k-1}{u}_j^k\mathrm{\Delta }x,\hfill \end{array}\end{array}$$

from which we obtain,$$\begin{array}{c}\hfill \begin{array}{cc}& |{\sum }_{j=1}^N{\mathit{\beta }}_j^{k-1}(u_j^k-u_j^{k-1}){\mathrm{\Delta }x|\le \Vert \mathit{\beta }\Vert }_{\infty }\frac{\mathrm{\Delta }t}{\mathrm{\Delta }x}{\sum }_{j=1}^N(|g_j^{k-1}-g_{j-1}^{k-1}|u_{j-1}^k+|u_j^k-u_{j-1}^k|g_{j-1}^{k-1})\mathrm{\Delta }x\hfill \\ \hfill & +{\mathrm{\Delta }t\Vert \mathit{\beta }\Vert }_{\infty }{\sum }_{j=1}^N{u}_j^k{\sum }_{m=1}^N{\mathit{\eta }}_{m,j}^{k-1}{u}_m^{k-1}\mathrm{\Delta }x\mathrm{\Delta }x\hfill \\ \hfill & +\frac{\mathrm{\Delta }t}{2}{\Vert \mathit{\beta }\Vert }_{\infty }{\sum }_{j=1}^N{\sum }_{s=1}^{j-1}{\mathit{\eta }}_{j-s,s}^{k-1}{u}_{j-s}^k{u}_s^{k-1}\mathrm{\Delta }x\mathrm{\Delta }x+{\mathrm{\Delta }t\Vert \mathit{\beta }\Vert }_{\infty }{\Vert R\Vert }_{\infty }{\Vert u^k\Vert }_1\hfill \\ \hfill & \le {\Vert \mathit{\beta }\Vert }_{\infty }{B}_{4,1}\Vert u^k{\Vert }_{\infty }{x}_{max}\mathrm{\Delta }t+c_0\mathrm{\Delta }t{\sum }_{j=1}^N|u_j^k-u_{j-1}^k|\hfill \\ \hfill & +\frac{3}{2}{\Vert \mathit{\beta }\Vert }_{\infty }{\Vert \mathit{\eta }\Vert }_{\infty }\Vert u^k{\Vert }_1\Vert u^{k-1}{\Vert }_1\mathrm{\Delta }t+{\mathrm{\Delta }t\Vert \mathit{\beta }\Vert }_{\infty }{\Vert R\Vert }_{\infty }{\Vert u^k\Vert }_1\hfill \\ \hfill & \le {\Vert \mathit{\beta }\Vert }_{\infty }{B}_{4,1}{B}_3{x}_{max}\mathrm{\Delta }t+c_0\mathrm{\Delta }tTV(u^k)+\frac{3}{2}{\Vert \mathit{\beta }\Vert }_{\infty }{\Vert \mathit{\eta }\Vert }_{\infty }{B}_2^2\mathrm{\Delta }t+{\mathrm{\Delta }t\Vert \mathit{\beta }\Vert }_{\infty }{\Vert R\Vert }_{\infty }{B}_2.\hfill \end{array}\end{array}$$

Let $C_{4,1}={\Vert \mathit{\beta }\Vert }_{\infty }{B}_{4,1}{B}_3{x}_{max}+{3\Vert \mathit{\beta }\Vert }_{\infty }{\Vert \mathit{\eta }\Vert }_{\infty }{B}_2^2/2+{\Vert \mathit{\beta }\Vert }_{\infty }{\Vert R\Vert }_{\infty }{B}_2$ and from Lemmas 4.1 through 4.3 and previous results established in this Lemma we obtain$$\begin{array}{c}\hfill \begin{array}{cc}|{\sum }_{j=1}^N{\mathit{\beta }}_j^{k-1}(u_j^k-u_j^{k-1})\mathrm{\Delta }x|\hfill & \le (C_{4,1}+c_{0}TV(u^k))\mathrm{\Delta }t.\hfill \end{array}\end{array}$$

By (A5(A5) $$\begin{array}{c}\hfill \begin{array}{cc}|u_0^{k+1}-u_0^k|\hfill & =\frac{1}{g_0^{k-1}}|F^k-F^{k-1}|+\frac{1}{g_0^{k-1}}{\sum }_{j=1}^N|{\mathit{\beta }}_j^k-{\mathit{\beta }}_j^{k-1}|u_j^k\mathrm{\Delta }x\hfill \\ \hfill & +\frac{1}{g_0^{k-1}}|{\sum }_{j=1}^N{\mathit{\beta }}_j^{k-1}(u_j^k-u_j^{k-1})\mathrm{\Delta }x|+\frac{1}{g_0^{k-1}}|g_0^k-g_0^{k-1}|u_0^{k+1}.\hfill \end{array}\end{array}$$(A5) ) we have$$\begin{array}{c}\hfill \begin{array}{cc}|u_0^{k+1}-u_0^k|\hfill & \le \frac{1}{\mathit{\zeta }}(B_{4,1}+B_{4,1}{B}_2+B_3{B}_{4,1}+C_{4,1}+c_{0}TV(u^k))\mathrm{\Delta }t.\hfill \end{array}\end{array}$$

This completes the proof.$\square$

Proof of Lemma 4.5:

From the first row of (4.1(4.1) $$\begin{array}{cc}& \frac{u_j^{k+1}-u_j^k}{\mathrm{\Delta }t}+\frac{g_j^k{u}_j^{k+1}-g_{j-1}^k{u}_{j-1}^{k+1}}{\mathrm{\Delta }x}+u_j^{k+1}\sum _{m=1}^N{\mathit{\eta }}_{m,j}^k{u}_m^k\mathrm{\Delta }x+R_j^k{u}_j^{k+1}\hfill \\ \hfill & =\frac{1}{2}\sum _{s=1}^{j-1}{\mathit{\eta }}_{j-s,s}^k{u}_{j-s}^k{u}_s^{k+1}\mathrm{\Delta }x,\hfill \\ \hfill & g_0^k{u}_0^{k+1}=F^k+\sum _{j=1}^N{\mathit{\beta }}_j^k{u}_j^k\mathrm{\Delta }x,u_j^0=u_0(x_j),\hfill \\ \hfill & \frac{y_{\ell }^{k+1}-y_{\ell }^k}{\mathrm{\Delta }t}=f_{\ell }^k-h_{\ell }^k{y}_{\ell }^{k+1},y_{\ell }^0=y_{\ell ,0},\ell =1,\dots ,m.\hfill \end{array}$$(4.1) ), we obtain$$\begin{array}{c}\hfill \begin{array}{cc}& (1+\frac{\mathrm{\Delta }t}{\mathrm{\Delta }x}{g}_j^k+\mathrm{\Delta }t{R}_j^k+\mathrm{\Delta }t\frac{g_j^k-g_{j-1}^k}{\mathrm{\Delta }x}+\mathrm{\Delta }t{\sum }_{m=1}^N{\mathit{\eta }}_{m,j+1}^k{u}_m^k\mathrm{\Delta }x)(u_{j+1}^{k+1}-u_j^{k+1})\hfill \\ \hfill & +\frac{\mathrm{\Delta }t}{\mathrm{\Delta }x}(g_{j+1}^k-2g_j^k+g_{j-1}^k)u_{j+1}^{k+1}-\frac{\mathrm{\Delta }t}{\mathrm{\Delta }x}(u_j^{k+1}-u_{j-1}^{k+1})g_{j-1}^k\hfill \\ \hfill & +\mathrm{\Delta }t{u}_j^{k+1}{\sum }_{m=1}^N({\mathit{\eta }}_{m,j+1}^k-{\mathit{\eta }}_{m,j}^k)u_m^k\mathrm{\Delta }x\hfill \\ \hfill & =(u_{j+1}^k-u_j^k)+\frac{\mathrm{\Delta }t}{2}{\mathit{\eta }}_{1,j}^k{u}_1^k{u}_j^{k+1}\mathrm{\Delta }x-\mathrm{\Delta }t(R_{j+1}^k-R_j^k)u_{j+1}^{k+1}\hfill \\ \hfill & +\frac{\mathrm{\Delta }t}{2}{\sum }_{s=1}^{j-1}\left(({\mathit{\eta }}_{j+1-s,s}^k-{\mathit{\eta }}_{j-s,s}^k)u_{j+1-s}^k+(u_{j+1-s}^k-u_{j-s}^k){\mathit{\eta }}_{j-s,s}^k\right)u_s^{k+1}\mathrm{\Delta }x.\hfill \end{array}\end{array}$$

Applying the triangle inequality, nonnegativity of the numerical scheme, and assumptions (B1) and (B2), we have$$\begin{array}{c}\hfill \begin{array}{c}(1-{\mathit{\omega }}_1\mathrm{\Delta }t+\frac{\mathrm{\Delta }t}{\mathrm{\Delta }x}{g}_j^k+\mathrm{\Delta }t{\sum }_{m=1}^N{\mathit{\eta }}_{m,j+1}^k{u}_m^k\mathrm{\Delta }x)|u_{j+1}^{k+1}-u_j^{k+1}|-\frac{\mathrm{\Delta }t}{\mathrm{\Delta }x}|u_j^{k+1}-u_{j-1}^{k+1}|g_{j-1}^k\hfill \\ \le |u_{j+1}^k-u_j^k|+\frac{\mathrm{\Delta }t}{2}{\mathit{\eta }}_{1,j}^k{u}_1^k{u}_j^{k+1}\mathrm{\Delta }x\hfill \\ +\frac{\mathrm{\Delta }t}{2}{\sum }_{s=1}^{j-1}|{\mathit{\eta }}_{j+1-s,s}^k-{\mathit{\eta }}_{j-s,s}^k|u_{j+1-s}^k{u}_s^{k+1}\mathrm{\Delta }x+\frac{\mathrm{\Delta }t}{2}{\sum }_{s=1}^{j-1}{\mathit{\eta }}_{j-s,s}^k|u_{j+1-s}^k-u_{j-s}^k|u_s^{k+1}\mathrm{\Delta }x\hfill \\ +\mathrm{\Delta }t{u}_j^{k+1}{\sum }_{m=1}^N|{\mathit{\eta }}_{m,j+1}^k-{\mathit{\eta }}_{m,j}^k|u_m^k\mathrm{\Delta }x+\frac{\mathrm{\Delta }t}{\mathrm{\Delta }x}|g_{j+1}^k-2g_j^k+g_{j-1}^k|u_{j+1}^{k+1}+\mathrm{\Delta }t|R_{j+1}^k-R_j^k|u_{j+1}^{k+1}.\hfill \end{array}\end{array}$$

Noting that ${\mathrm{\Sigma }}_{m=1}^N{\mathit{\eta }}_{m,j+1}^k{u}_m^k\mathrm{\Delta }x\ge 0$, then summing from $j=1$ to $N-1$,$$\begin{array}{c}\hfill \begin{array}{c}(1-{\mathit{\omega }}_1\mathrm{\Delta }t){\sum }_{j=1}^{N-1}|u_{j+1}^{k+1}-u_j^{k+1}|+\frac{\mathrm{\Delta }t}{\mathrm{\Delta }x}{g}_{N-1}^k|u_N^{k+1}-u_{N-1}^{k+1}|-\frac{\mathrm{\Delta }t}{\mathrm{\Delta }x}|u_1^{k+1}-u_0^{k+1}|g_0^k\hfill \\ \le {\sum }_{j=1}^{N-1}|u_{j+1}^k-u_j^k|+\frac{\mathrm{\Delta }t}{2}{\sum }_{j=1}^{N-1}{\mathit{\eta }}_{1,j}^k{u}_1^k{u}_j^{k+1}\mathrm{\Delta }x+\frac{\mathrm{\Delta }t}{2}{\sum }_{j=1}^{N-1}{\sum }_{s=1}^{j-1}|{\mathit{\eta }}_{j+1-s,s}^k-{\mathit{\eta }}_{j-s,s}^k|u_{j+1-s}^k{u}_s^{k+1}\mathrm{\Delta }x\hfill \\ +\frac{\mathrm{\Delta }t}{2}{\sum }_{j=1}^{N-1}{\sum }_{s=1}^{j-1}{\mathit{\eta }}_{j-s,s}^k|u_{j+1-s}^k-u_{j-s}^k|u_s^{k+1}\mathrm{\Delta }x+\mathrm{\Delta }t{\sum }_{j=1}^{N-1}{u}_j^{k+1}{\sum }_{m=1}^N|{\mathit{\eta }}_{m,j+1}^k-{\mathit{\eta }}_{m,j}^k|u_m^k\mathrm{\Delta }x\hfill \\ +\frac{\mathrm{\Delta }t}{\mathrm{\Delta }x}{\sum }_{j=1}^{N-1}|g_{j+1}^k-2g_j^k+g_{j-1}^k|u_{j+1}^{k+1}+\mathrm{\Delta }t{\sum }_{j=1}^{N-1}|R_{j+1}^k-R_j^k|u_{j+1}^{k+1}.\hfill \end{array}\end{array}$$

By Lemmas 4.1 and 4.4,$$\begin{array}{c}\hfill \begin{array}{cc}& (1-{\mathit{\omega }}_1\mathrm{\Delta }t){\sum }_{j=1}^{N-1}|u_{j+1}^{k+1}-u_j^{k+1}|\le {\sum }_{j=1}^{N-1}|u_{j+1}^k-u_j^k|+\frac{\mathrm{\Delta }t}{2}{\Vert \mathit{\eta }\Vert }_{\infty }\Vert u^k{\Vert }_{\infty }{\Vert u^{k+1}\Vert }_1\hfill \\ \hfill & +\frac{\mathrm{\Delta }t}{2}{B}_{4,1}\Vert u^k{\Vert }_1\Vert u^{k+1}{\Vert }_1+\frac{\mathrm{\Delta }t}{2}{\Vert \mathit{\eta }\Vert }_{\infty }\Vert u^{k+1}{\Vert }_1{\sum }_{j=1}^{N-1}|u_{j+1}^k-u_j^k|\hfill \\ \hfill & +\mathrm{\Delta }t{B}_{4,1}{\Vert \mathit{\eta }\Vert }_{\infty }\Vert u^{k+1}{\Vert }_1\Vert u^k{\Vert }_1+\mathrm{\Delta }t{B}_{4,1}\Vert u^{k+1}{\Vert }_1+\frac{\mathrm{\Delta }t}{\mathrm{\Delta }x}|u_1^{k+1}-u_0^{k+1}|g_0^k+\mathrm{\Delta }t{B}_{4,1}{\Vert u^{k+1}\Vert }_1.\hfill \end{array}\end{array}$$

By assumption (B2), Lemmas 4.2 and 4.3,(A6) $$\begin{array}{c}\hfill \begin{array}{cc}& (1-{\mathit{\omega }}_1\mathrm{\Delta }t){\sum }_{j=1}^{N-1}|u_{j+1}^{k+1}-u_j^{k+1}|\le {\sum }_{j=1}^{N-1}|u_{j+1}^k-u_j^k|+\frac{\mathrm{\Delta }t}{2}{c}_0{B}_2{B}_3\hfill \\ \hfill & +\frac{\mathrm{\Delta }t}{2}{B}_2{B}_3{B}_{4,1}+\frac{\mathrm{\Delta }t}{2}{c}_0{B}_2{\sum }_{j=1}^{N-1}|u_{j+1}^k-u_j^k|+\mathrm{\Delta }t{c}_0{B}_{4,1}{B}_2^2+\mathrm{\Delta }t{B}_2{B}_{4,1}\hfill \\ \hfill & +\frac{\mathrm{\Delta }t}{\mathrm{\Delta }x}|u_1^{k+1}-u_0^{k+1}|g_0^k+\mathrm{\Delta }t{B}_{4,1}{B}_2.\hfill \end{array}\end{array}$$(A6)

Adding $|u_1^{k+1}-u_0^{k+1}|$ to both sides of (A6(A6) $$\begin{array}{c}\hfill \begin{array}{cc}& (1-{\mathit{\omega }}_1\mathrm{\Delta }t){\sum }_{j=1}^{N-1}|u_{j+1}^{k+1}-u_j^{k+1}|\le {\sum }_{j=1}^{N-1}|u_{j+1}^k-u_j^k|+\frac{\mathrm{\Delta }t}{2}{c}_0{B}_2{B}_3\hfill \\ \hfill & +\frac{\mathrm{\Delta }t}{2}{B}_2{B}_3{B}_{4,1}+\frac{\mathrm{\Delta }t}{2}{c}_0{B}_2{\sum }_{j=1}^{N-1}|u_{j+1}^k-u_j^k|+\mathrm{\Delta }t{c}_0{B}_{4,1}{B}_2^2+\mathrm{\Delta }t{B}_2{B}_{4,1}\hfill \\ \hfill & +\frac{\mathrm{\Delta }t}{\mathrm{\Delta }x}|u_1^{k+1}-u_0^{k+1}|g_0^k+\mathrm{\Delta }t{B}_{4,1}{B}_2.\hfill \end{array}\end{array}$$(A6) ) we obtain$$\begin{array}{c}\hfill \begin{array}{cc}& (1-{\mathit{\omega }}_1\mathrm{\Delta }t)TV(u^{k+1})\le {\sum }_{j=1}^{N-1}|u_{j+1}^k-u_j^k|+\frac{\mathrm{\Delta }t}{2}{c}_0{B}_2{B}_3+\frac{\mathrm{\Delta }t}{2}{B}_2{B}_3{B}_{4,1}\hfill \\ \hfill & +\frac{\mathrm{\Delta }t}{2}{c}_0{B}_2{\sum }_{j=1}^{N-1}|u_{j+1}^k-u_j^k|+\mathrm{\Delta }t{c}_0{B}_{4,1}{B}_2^2+2B_2{B}_{4,1}\mathrm{\Delta }t+(1+\frac{\mathrm{\Delta }t}{\mathrm{\Delta }x}{g}_0^k)|u_1^{k+1}-u_0^{k+1}|.\hfill \end{array}\end{array}$$

By Lemma 4.4,$$\begin{array}{c}\hfill \begin{array}{cc}(1-{\mathit{\omega }}_1\mathrm{\Delta }t)TV(u^{k+1})\hfill & \le {\sum }_{j=1}^{N-1}|u_{j+1}^k-u_j^k|+\frac{\mathrm{\Delta }t}{2}{c}_0{B}_2{B}_3+\frac{\mathrm{\Delta }t}{2}{B}_2{B}_3{B}_{4,1}\hfill \\ \hfill & +\frac{\mathrm{\Delta }t}{2}{c}_0{B}_2{\sum }_{j=1}^{N-1}|u_{j+1}^k-u_j^k|+\mathrm{\Delta }t{c}_0{B}_{4,1}{B}_2^2+2B_2{B}_{4,1}\mathrm{\Delta }t\hfill \\ \hfill & +|u_1^k-u_0^k|+B_3(c_0{B}_2+B_{4,1}+B_{1,1})\mathrm{\Delta }t+B_{4,2}(1+TV(u^k))\mathrm{\Delta }t\hfill \\ \hfill & \le (1+(B_{4,2}+\frac{c_0{B}_2}{2})\mathrm{\Delta }t)TV(u^k)+C_{\text{TV1}}\mathrm{\Delta }t.\hfill \end{array}\end{array}$$

where $C_{\text{TV1}}=(c_0{B}_2{B}_3/2)+(B_2{B}_3{B}_{4,1}/2)+(c_0{B}_{4,1}{B}_2^2)+2B_2{B}_{4,1}+B_3(c_0{B}_2+B_{4,1}+B_{1,1})+B_{4,2}$. Let $C_{\text{TV2}}=B_{4,2}+c_0{B}_2/2$ to arrive at the following recurrence relationship,(A7) $$\begin{array}{c}\hfill \begin{array}{cc}(1-{\mathit{\omega }}_1\mathrm{\Delta }t)TV(u^{k+1})\hfill & \le (1+C_{\text{TV2}}\mathrm{\Delta }t)TV(u^k)+C_{\text{TV1}}\mathrm{\Delta }t.\hfill \end{array}\end{array}$$(A7)

The constant $B_5$ is established using a standard induction argument on (A7(A7) $$\begin{array}{c}\hfill \begin{array}{cc}(1-{\mathit{\omega }}_1\mathrm{\Delta }t)TV(u^{k+1})\hfill & \le (1+C_{\text{TV2}}\mathrm{\Delta }t)TV(u^k)+C_{\text{TV1}}\mathrm{\Delta }t.\hfill \end{array}\end{array}$$(A7) ).$\square$

Proof of Lemma 4.6:

Multiplying the first row in the numerical scheme (4.1(4.1) $$\begin{array}{cc}& \frac{u_j^{k+1}-u_j^k}{\mathrm{\Delta }t}+\frac{g_j^k{u}_j^{k+1}-g_{j-1}^k{u}_{j-1}^{k+1}}{\mathrm{\Delta }x}+u_j^{k+1}\sum _{m=1}^N{\mathit{\eta }}_{m,j}^k{u}_m^k\mathrm{\Delta }x+R_j^k{u}_j^{k+1}\hfill \\ \hfill & =\frac{1}{2}\sum _{s=1}^{j-1}{\mathit{\eta }}_{j-s,s}^k{u}_{j-s}^k{u}_s^{k+1}\mathrm{\Delta }x,\hfill \\ \hfill & g_0^k{u}_0^{k+1}=F^k+\sum _{j=1}^N{\mathit{\beta }}_j^k{u}_j^k\mathrm{\Delta }x,u_j^0=u_0(x_j),\hfill \\ \hfill & \frac{y_{\ell }^{k+1}-y_{\ell }^k}{\mathrm{\Delta }t}=f_{\ell }^k-h_{\ell }^k{y}_{\ell }^{k+1},y_{\ell }^0=y_{\ell ,0},\ell =1,\dots ,m.\hfill \end{array}$$(4.1) ) by $\mathrm{\Delta }x$ we obtain$$\begin{array}{c}\hfill \begin{array}{cc}\frac{u_j^{k+1}-u_j^k}{\mathrm{\Delta }t}\mathrm{\Delta }x\hfill & =-((g_j^k-g_{j-1}^k)u_{j-1}^{k+1}+(u_j^{k+1}-u_{j-1}^{k+1})g_{j-1}^k)-R_j^k{u}_j^{k+1}\mathrm{\Delta }x\hfill \\ \hfill & -u_j^{k+1}{\sum }_{m=1}^N{\mathit{\eta }}_{m,j}^k{u}_m^k\mathrm{\Delta }x\mathrm{\Delta }x+\frac{1}{2}{\sum }_{s=1}^{j-1}{\mathit{\eta }}_{j-s,s}^k{u}_{j-s}^k{u}_s^{k+1}\mathrm{\Delta }x\mathrm{\Delta }x\hfill \end{array}\end{array}$$

Applying the triangle inequality, summing from $j=1,\dots ,N$, and applying Lemma 4.4,(A8) $$\begin{array}{c}\hfill \begin{array}{cc}{\sum }_{j=1}^N\left|\frac{u_j^{k+1}-u_j^k}{\mathrm{\Delta }t}\right|\mathrm{\Delta }x\hfill & \le (B_{4,1}+{\Vert R\Vert }_{\infty })\Vert u^{k+1}{\Vert }_1+B_{1,1}TV(u^{k+1})\hfill \\ \hfill & +\frac{3}{2}{\Vert \mathit{\eta }\Vert }_{\infty }\Vert u^k{\Vert }_1{\Vert u^{k+1}\Vert }_1\le \overline{L}.\hfill \end{array}\end{array}$$(A8)

It clearly follows that$$\begin{array}{c}\hfill \begin{array}{c}{\sum }_{j=1}^N\left|\frac{u_j^r-u_j^q}{\mathrm{\Delta }t}\right|\mathrm{\Delta }x\le {\sum }_{j=1}^N{\sum }_{k=q}^{r-1}\left|\frac{u_j^{k+1}-u_j^k}{\mathrm{\Delta }t}\right|\mathrm{\Delta }x\le \overline{L}(r-q)\hfill \end{array}\end{array}$$

and by Lemma 4.4,$$\begin{array}{c}\hfill \begin{array}{c}{\sum }_{\ell =1}^m\left|\frac{y_{\ell }^r-y_{\ell }^q}{\mathrm{\Delta }t}\right|\le {\sum }_{\ell =1}^m{\sum }_{k=q}^{r-1}\left|\frac{y_{\ell }^{k+1}-y_{\ell }^k}{\mathrm{\Delta }t}\right|\le m{B}_{4,1}(r-q)\hfill \end{array}\end{array}$$

Let $B_6=max\{\overline{L},m{B}_{4,1}\}$ to complete the proof.$\square$

Proof of Theorem 5.1:

The result for the population density u can be established using similar techniques as in the proof of Lemma 16.10 (p.279) in [41]. By Lemma 4.6 each component ($\ell =1,\dots ,m$) in the sequence ${\overrightarrow{y}}^k$ approximating the environment satisfies a Lipschitz condition on [0, T] and by the assumptions on $f_{\ell }$ and $h_{\ell }$, each component of the approximating sequence defines a bounded equi-Lipschitz family. Thus the family $\{Y_{\ell ,\mathrm{\Delta }t}\}$ is equicontinuous and by the Arzela–Ascoli Theorem one can extract a uniformly convergent subsequence in C[0, T].$\square$

Proof of Theorem 5.2:

First observe that convergence of a sequence $y_{\ell }^k$ in C[0, T] is uniform convergence. Applying the Picard Iteration to the limit function $\overrightarrow{y}$ as outlined in [64], (Pgs 112-127) we can show that $\overrightarrow{y}$ is a weak solution to our ODE system. So it suffices to establish the weak solution result for u. Let $\mathit{\rho }\in C^1({\mathbb{D}}_x)$ and denote the values of $\mathit{\rho }$ on the numerical mesh $\mathit{\rho }(x_j,t^k)$ by ${\mathit{\rho }}_j^k$. Multiplying the first equation of the difference scheme (4.1(4.1) $$\begin{array}{cc}& \frac{u_j^{k+1}-u_j^k}{\mathrm{\Delta }t}+\frac{g_j^k{u}_j^{k+1}-g_{j-1}^k{u}_{j-1}^{k+1}}{\mathrm{\Delta }x}+u_j^{k+1}\sum _{m=1}^N{\mathit{\eta }}_{m,j}^k{u}_m^k\mathrm{\Delta }x+R_j^k{u}_j^{k+1}\hfill \\ \hfill & =\frac{1}{2}\sum _{s=1}^{j-1}{\mathit{\eta }}_{j-s,s}^k{u}_{j-s}^k{u}_s^{k+1}\mathrm{\Delta }x,\hfill \\ \hfill & g_0^k{u}_0^{k+1}=F^k+\sum _{j=1}^N{\mathit{\beta }}_j^k{u}_j^k\mathrm{\Delta }x,u_j^0=u_0(x_j),\hfill \\ \hfill & \frac{y_{\ell }^{k+1}-y_{\ell }^k}{\mathrm{\Delta }t}=f_{\ell }^k-h_{\ell }^k{y}_{\ell }^{k+1},y_{\ell }^0=y_{\ell ,0},\ell =1,\dots ,m.\hfill \end{array}$$(4.1) ) by ${\mathit{\rho }}_j^{k+1}$, and adding and subtracting several terms:$$\begin{array}{c}\hfill \begin{array}{cc}& \frac{u_j^{k+1}{\mathit{\rho }}_j^{k+1}-u_j^k{\mathit{\rho }}_j^k}{\mathrm{\Delta }t}+\frac{g_j^k{u}_j^{k+1}{\mathit{\rho }}_j^{k+1}-g_{j-1}^k{u}_{j-1}^{k+1}{\mathit{\rho }}_{j-1}^{k+1}}{\mathrm{\Delta }x}-\frac{{\mathit{\rho }}_j^{k+1}-{\mathit{\rho }}_j^k}{\mathrm{\Delta }t}{u}_j^k\hfill \\ \hfill & -\frac{{\mathit{\rho }}_j^{k+1}-{\mathit{\rho }}_{j-1}^{k+1}}{\mathrm{\Delta }x}{g}_{j-1}^k{u}_{j-1}^{k+1}+{\mathit{\rho }}_j^{k+1}{R}_j^k{u}_j^{k+1}+{\mathit{\rho }}_j^{k+1}{u}_j^{k+1}{\sum }_{m=1}^N{\mathit{\eta }}_{m,j}^k{u}_m^k\mathrm{\Delta }x\hfill \\ \hfill & =\frac{1}{2}{\mathit{\rho }}_j^{k+1}{\sum }_{s=1}^{j-1}{\mathit{\eta }}_{j-s,s}^k{u}_{j-s}^k{u}_s^{k+1}\mathrm{\Delta }x.\hfill \end{array}\end{array}$$

Multiplying by $\mathrm{\Delta }x$ and $\mathrm{\Delta }t$ and summing from $k=0,1,\dots ,L-1$, and $j=1,2,\dots ,N$, evaluating the telescoping sums, noting that $g_N^k=0$, and inserting the corner boundary condition we have(A9) $$\begin{array}{c}\hfill \begin{array}{cc}& {\sum }_{j=1}^N(u_j^L{\mathit{\rho }}_j^L-u_j^0{\mathit{\rho }}_j^0)\mathrm{\Delta }x-{\sum }_{k=0}^{L-1}{\sum }_{j=1}^N\frac{{\mathit{\rho }}_j^{k+1}-{\mathit{\rho }}_j^k}{\mathrm{\Delta }t}{u}_j^k\mathrm{\Delta }x\mathrm{\Delta }t\hfill \\ \hfill & -{\sum }_{k=0}^{L-1}{\mathit{\rho }}_0^{k+1}(F+{\sum }_{j=1}^N{\mathit{\beta }}_j^k{u}_j^k\mathrm{\Delta }x)\mathrm{\Delta }t-{\sum }_{k=0}^{L-1}{\sum }_{j=1}^N\frac{{\mathit{\rho }}_j^{k+1}-{\mathit{\rho }}_{j-1}^{k+1}}{\mathrm{\Delta }x}{g}_{j-1}^k{u}_{j-1}^{k+1}\mathrm{\Delta }x\mathrm{\Delta }t\hfill \\ \hfill & +{\sum }_{k=0}^{L-1}{\sum }_{j=1}^N{\mathit{\rho }}_j^{k+1}{R}_j^k{u}_j^{k+1}\mathrm{\Delta }x\mathrm{\Delta }t\hfill \\ \hfill & =\frac{1}{2}{\sum }_{k=0}^{L-1}{\sum }_{j=1}^N{\mathit{\rho }}_j^{k+1}{\sum }_{s=1}^{j-1}{\mathit{\eta }}_{p,j-s,s,l}^k{u}_{p,j-s,l}^k{u}_{p,s,l}^{k+1}\mathrm{\Delta }x\mathrm{\Delta }x\mathrm{\Delta }t\hfill \\ \hfill & -{\sum }_{k=0}^{L-1}{\sum }_{j=1}^N{\mathit{\rho }}_j^{k+1}{u}_j^{k+1}{\sum }_{m=1}^N{\mathit{\eta }}_{m,j}^k{u}_m^k\mathrm{\Delta }x\mathrm{\Delta }x\mathrm{\Delta }t\hfill \end{array}\end{array}$$(A9)

Using similar techniques as in the proof of Lemma 16.10 in [41] (p.279) we can show that the above Equation (A9(A9) $$\begin{array}{c}\hfill \begin{array}{cc}& {\sum }_{j=1}^N(u_j^L{\mathit{\rho }}_j^L-u_j^0{\mathit{\rho }}_j^0)\mathrm{\Delta }x-{\sum }_{k=0}^{L-1}{\sum }_{j=1}^N\frac{{\mathit{\rho }}_j^{k+1}-{\mathit{\rho }}_j^k}{\mathrm{\Delta }t}{u}_j^k\mathrm{\Delta }x\mathrm{\Delta }t\hfill \\ \hfill & -{\sum }_{k=0}^{L-1}{\mathit{\rho }}_0^{k+1}(F+{\sum }_{j=1}^N{\mathit{\beta }}_j^k{u}_j^k\mathrm{\Delta }x)\mathrm{\Delta }t-{\sum }_{k=0}^{L-1}{\sum }_{j=1}^N\frac{{\mathit{\rho }}_j^{k+1}-{\mathit{\rho }}_{j-1}^{k+1}}{\mathrm{\Delta }x}{g}_{j-1}^k{u}_{j-1}^{k+1}\mathrm{\Delta }x\mathrm{\Delta }t\hfill \\ \hfill & +{\sum }_{k=0}^{L-1}{\sum }_{j=1}^N{\mathit{\rho }}_j^{k+1}{R}_j^k{u}_j^{k+1}\mathrm{\Delta }x\mathrm{\Delta }t\hfill \\ \hfill & =\frac{1}{2}{\sum }_{k=0}^{L-1}{\sum }_{j=1}^N{\mathit{\rho }}_j^{k+1}{\sum }_{s=1}^{j-1}{\mathit{\eta }}_{p,j-s,s,l}^k{u}_{p,j-s,l}^k{u}_{p,s,l}^{k+1}\mathrm{\Delta }x\mathrm{\Delta }x\mathrm{\Delta }t\hfill \\ \hfill & -{\sum }_{k=0}^{L-1}{\sum }_{j=1}^N{\mathit{\rho }}_j^{k+1}{u}_j^{k+1}{\sum }_{m=1}^N{\mathit{\eta }}_{m,j}^k{u}_m^k\mathrm{\Delta }x\mathrm{\Delta }x\mathrm{\Delta }t\hfill \end{array}\end{array}$$(A9) ) converges (as $L\to \infty$) to Equation (3.1) for any convergent subsequence defined in Theorem 5.1. The $L^1$, $L^{\infty }$ and BV bounds on u follow by observing that $L^1$ convergence obtained in Theorem 5.1 implies pointwise convergence almost everywhere (along a subsequence) and then taking limits in the estimates established in Lemmas 4.2, 4.3 and 4.5.$\square$

Proof of Lemma 5.3:

We let $v_j^k=u_j^k-{\widehat{u}}_j^k$, and observe that $v_j^k$ satisfies(A10) $$\begin{array}{c}\hfill \begin{array}{cc}& \frac{v_j^{k+1}-v_j^k}{\mathrm{\Delta }t}+\frac{g_j^k{v}_j^{k+1}-g_{j-1}^k{v}_{j-1}^{k+1}}{\mathrm{\Delta }x}+\frac{(g_j^k-{\widehat{g}}_j^k){\widehat{u}}_j^{k+1}-(g_{j-1}^k-{\widehat{g}}_{j-1}^k){\widehat{u}}_{j-1}^{k+1}}{\mathrm{\Delta }x}\hfill \\ \hfill & =\frac{1}{2}{\sum }_{s=1}^{j-1}({\mathit{\eta }}_{j-s,s}^k-{\widehat{\mathit{\eta }}}_{j-s,s}^k)u_{j-s}^k{u}_s^{k+1}\mathrm{\Delta }x+\frac{1}{2}{\sum }_{s=1}^{j-1}{\widehat{\mathit{\eta }}}_{j-s,s}^k{v}_{j-s}^k{u}_s^{k+1}\mathrm{\Delta }x\hfill \\ \hfill & +\frac{1}{2}{\sum }_{s=1}^{j-1}{\widehat{\mathit{\eta }}}_{j-s,s}^k{\widehat{u}}_{j-s}^k{v}_s^{k+1}\mathrm{\Delta }x-{\widehat{R}}_j^k{v}_j^{k+1}-(R_j^k-{\widehat{R}}_j^k)u_j^{k+1}-v_j^{k+1}{\sum }_{m=1}^N{\widehat{\mathit{\eta }}}_{m,j}^k{\widehat{u}}_m^k\mathrm{\Delta }x\hfill \\ \hfill & -{\widehat{u}}_j^{k+1}{\sum }_{m=1}^N({\mathit{\eta }}_{m,j}^k-{\widehat{\mathit{\eta }}}_{m,j}^k)u_m^k\mathrm{\Delta }x-{\widehat{u}}_j^{k+1}{\sum }_{m=1}^N{\widehat{\mathit{\eta }}}_{m,j}^k{v}_m^k\mathrm{\Delta }x.\hfill \end{array}\end{array}$$(A10)

Multiplying both sides of (A10(A10) $$\begin{array}{c}\hfill \begin{array}{cc}& \frac{v_j^{k+1}-v_j^k}{\mathrm{\Delta }t}+\frac{g_j^k{v}_j^{k+1}-g_{j-1}^k{v}_{j-1}^{k+1}}{\mathrm{\Delta }x}+\frac{(g_j^k-{\widehat{g}}_j^k){\widehat{u}}_j^{k+1}-(g_{j-1}^k-{\widehat{g}}_{j-1}^k){\widehat{u}}_{j-1}^{k+1}}{\mathrm{\Delta }x}\hfill \\ \hfill & =\frac{1}{2}{\sum }_{s=1}^{j-1}({\mathit{\eta }}_{j-s,s}^k-{\widehat{\mathit{\eta }}}_{j-s,s}^k)u_{j-s}^k{u}_s^{k+1}\mathrm{\Delta }x+\frac{1}{2}{\sum }_{s=1}^{j-1}{\widehat{\mathit{\eta }}}_{j-s,s}^k{v}_{j-s}^k{u}_s^{k+1}\mathrm{\Delta }x\hfill \\ \hfill & +\frac{1}{2}{\sum }_{s=1}^{j-1}{\widehat{\mathit{\eta }}}_{j-s,s}^k{\widehat{u}}_{j-s}^k{v}_s^{k+1}\mathrm{\Delta }x-{\widehat{R}}_j^k{v}_j^{k+1}-(R_j^k-{\widehat{R}}_j^k)u_j^{k+1}-v_j^{k+1}{\sum }_{m=1}^N{\widehat{\mathit{\eta }}}_{m,j}^k{\widehat{u}}_m^k\mathrm{\Delta }x\hfill \\ \hfill & -{\widehat{u}}_j^{k+1}{\sum }_{m=1}^N({\mathit{\eta }}_{m,j}^k-{\widehat{\mathit{\eta }}}_{m,j}^k)u_m^k\mathrm{\Delta }x-{\widehat{u}}_j^{k+1}{\sum }_{m=1}^N{\widehat{\mathit{\eta }}}_{m,j}^k{v}_m^k\mathrm{\Delta }x.\hfill \end{array}\end{array}$$(A10) ) by $\mathrm{\Delta }x\mathrm{\Delta }t$ we obtain$$\begin{array}{c}\hfill \begin{array}{cc}& (1+\frac{\mathrm{\Delta }t}{\mathrm{\Delta }x}{g}_j^k+\mathrm{\Delta }t{\widehat{R}}_j^k+\mathrm{\Delta }t{\sum }_{m=1}^N{\widehat{\mathit{\eta }}}_{m,j}^k{\widehat{u}}_m^k\mathrm{\Delta }x)v_j^{k+1}\mathrm{\Delta }x-\frac{\mathrm{\Delta }t}{\mathrm{\Delta }x}{g}_{j-1}^k{v}_{j-1}^{k+1}\mathrm{\Delta }x\hfill \\ \hfill & =v_j^k\mathrm{\Delta }x-\mathrm{\Delta }t[(g_j^k-{\widehat{g}}_j^k)-(g_{j-1}^k-{\widehat{g}}_{j-1}^k)]{\widehat{u}}_j^{k+1}-\mathrm{\Delta }t(g_{j-1}^k-{\widehat{g}}_{j-1}^k)({\widehat{u}}_j^{k+1}-{\widehat{u}}_{j-1}^{k+1})\hfill \\ \hfill & +\frac{\mathrm{\Delta }t}{2}{\sum }_{s=1}^{j-1}({\mathit{\eta }}_{j-s,s}^k-{\widehat{\mathit{\eta }}}_{j-s,s}^k)u_{j-s}^k{u}_s^{k+1}\mathrm{\Delta }x\mathrm{\Delta }x+\frac{\mathrm{\Delta }t}{2}{\sum }_{s=1}^{j-1}{\widehat{\mathit{\eta }}}_{j-s,s}^k{v}_{j-s}^k{u}_s^{k+1}\mathrm{\Delta }x\mathrm{\Delta }x\hfill \\ \hfill & +\frac{\mathrm{\Delta }t}{2}{\sum }_{s=1}^{j-1}{\widehat{\mathit{\eta }}}_{j-s,s}^k{\widehat{u}}_{j-s}^k{v}_s^{k+1}\mathrm{\Delta }x\mathrm{\Delta }x-\mathrm{\Delta }t(R_j^k-{\widehat{R}}_j^k)u_j^{k+1}\mathrm{\Delta }x\hfill \\ \hfill & -\mathrm{\Delta }t{\widehat{u}}_j^{k+1}{\sum }_{m=1}^N({\mathit{\eta }}_{m,j}^k-{\widehat{\mathit{\eta }}}_{m,j}^k)u_m^k\mathrm{\Delta }x\mathrm{\Delta }x-\mathrm{\Delta }t{\widehat{u}}_j^{k+1}{\sum }_{m=1}^N{\widehat{\mathit{\eta }}}_{m,j}^k{v}_m^k\mathrm{\Delta }x\mathrm{\Delta }x.\hfill \end{array}\end{array}$$

Applying the triangle inequality and nonnegativity of the parameters ${\widehat{R}}_j^k$, ${\widehat{\mathit{\eta }}}_{m,j}^k$, and the numerical scheme ${\widehat{u}}_m^k$ we obtain,$$\begin{array}{c}\hfill \begin{array}{cc}& (1+\frac{\mathrm{\Delta }t}{\mathrm{\Delta }x}{g}_j^k+\mathrm{\Delta }t{\widehat{R}}_j^k+\mathrm{\Delta }t{\sum }_{m=1}^N{\widehat{\mathit{\eta }}}_{m,j}^k{\widehat{u}}_m^k\mathrm{\Delta }x)|v_j^{k+1}|\mathrm{\Delta }x-\frac{\mathrm{\Delta }t}{\mathrm{\Delta }x}{g}_{j-1}^k|v_{j-1}^{k+1}|\mathrm{\Delta }x\hfill \\ \hfill & \le |v_j^k|\mathrm{\Delta }x+\mathrm{\Delta }t|(g_j^k-{\widehat{g}}_j^k)-(g_{j-1}^k-{\widehat{g}}_{j-1}^k)|u_j^{k+1}+\mathrm{\Delta }t|g_{j-1}^k-{\widehat{g}}_{j-1}^k||{\widehat{u}}_j^{k+1}-{\widehat{u}}_{j-1}^{k+1}|\hfill \\ \hfill & +\frac{\mathrm{\Delta }t}{2}{\sum }_{s=1}^{j-1}|{\mathit{\eta }}_{j-s,s}^k-{\widehat{\mathit{\eta }}}_{j-s,s}^k|u_{j-s}^k{u}_s^{k+1}\mathrm{\Delta }x\mathrm{\Delta }x+\frac{\mathrm{\Delta }t}{2}{\sum }_{s=1}^{j-1}{\widehat{\mathit{\eta }}}_{j-s,s}^k|v_{j-s}^k|u_s^{k+1}\mathrm{\Delta }x\mathrm{\Delta }x\hfill \\ \hfill & +\frac{\mathrm{\Delta }t}{2}{\sum }_{s=1}^{j-1}{\widehat{\mathit{\eta }}}_{j-s,s}^k{\widehat{u}}_{j-s}^k|v_s^{k+1}|\mathrm{\Delta }x\mathrm{\Delta }x+\mathrm{\Delta }t|R_j^k-{\widehat{R}}_j^k|u_j^{k+1}\mathrm{\Delta }x\hfill \\ \hfill & +\mathrm{\Delta }t{\widehat{u}}_j^{k+1}{\sum }_{m=1}^N|{\mathit{\eta }}_{m,j}^k-{\widehat{\mathit{\eta }}}_{m,j}^k|u_m^k\mathrm{\Delta }x\mathrm{\Delta }x+\mathrm{\Delta }t{\widehat{u}}_j^{k+1}{\sum }_{m=1}^N{\widehat{\mathit{\eta }}}_{m,j}^k|v_m^k|\mathrm{\Delta }x\mathrm{\Delta }x.\hfill \end{array}\end{array}$$

Applying nonnegativity on the parameters and summing from $j=1$ to N,$$\begin{array}{c}\hfill \begin{array}{cc}& {\sum }_{j=1}^N|v_j^{k+1}|\mathrm{\Delta }x+\frac{\mathrm{\Delta }t}{\mathrm{\Delta }x}{\widehat{g}}_N^k|v_N^{k+1}|\mathrm{\Delta }x-\mathrm{\Delta }t{\widehat{g}}_0^k|v_0^{k+1}|\hfill \\ \hfill & \le {\sum }_{j=1}^N|v_j^k|\mathrm{\Delta }x+\mathrm{\Delta }t{\sum }_{j=1}^N|(g_j^k-{\widehat{g}}_j^k)-(g_{j-1}^k-{\widehat{g}}_{j-1}^k)|u_j^{k+1}\hfill \\ \hfill & +\mathrm{\Delta }t{\sum }_{j=1}^N|g_{j-1}^k-{\widehat{g}}_{j-1}^k||{\widehat{u}}_j^{k+1}-{\widehat{u}}_{j-1}^{k+1}|\hfill \\ \hfill & +\frac{\mathrm{\Delta }t}{2}{\sum }_{j=1}^N{\sum }_{s=1}^{j-1}|{\mathit{\eta }}_{j-s,s}^k-{\widehat{\mathit{\eta }}}_{j-s,s}^k|u_{j-s}^k{u}_s^{k+1}\mathrm{\Delta }x\mathrm{\Delta }x+\frac{\mathrm{\Delta }t}{2}{\sum }_{j=1}^N{\sum }_{s=1}^{j-1}{\widehat{\mathit{\eta }}}_{j-s,s}^k|v_{j-s}^k|u_s^{k+1}\mathrm{\Delta }x\mathrm{\Delta }x\hfill \\ \hfill & +\frac{\mathrm{\Delta }t}{2}{\sum }_{j=1}^N{\sum }_{s=1}^{j-1}{\widehat{\mathit{\eta }}}_{j-s,s}^k{\widehat{u}}_{j-s}^k|v_s^{k+1}|\mathrm{\Delta }x\mathrm{\Delta }x+\mathrm{\Delta }t{\sum }_{j=1}^N|R_j^k-{\widehat{R}}_j^k|u_j^{k+1}\mathrm{\Delta }x\hfill \\ \hfill & +\mathrm{\Delta }t{\sum }_{j=1}^N{\widehat{u}}_j^{k+1}{\sum }_{m=1}^N|{\mathit{\eta }}_{m,j}^k-{\widehat{\mathit{\eta }}}_{m,j}^k|u_m^k\mathrm{\Delta }x\mathrm{\Delta }x+\mathrm{\Delta }t{\sum }_{j=1}^N{\widehat{u}}_j^{k+1}{\sum }_{m=1}^N{\widehat{\mathit{\eta }}}_{m,j}^k|v_m^k|\mathrm{\Delta }x\mathrm{\Delta }x.\hfill \end{array}\end{array}$$

For convenience, let $G^k=|{\mathit{\varphi }}^k-{\widehat{\mathit{\varphi }}}^k|+{\mathrm{\Sigma }}_{\ell =1}^m|y_{\ell }^k-{\widehat{y}}_{\ell }^k|$. Applying the triangle inequality and the Lipschitz conditions on the model parameters we obtain,$$\begin{array}{c}\hfill \begin{array}{cc}& {\sum }_{j=1}^N|v_j^{k+1}|\mathrm{\Delta }x={\sum }_{j=1}^N|v_j^k|\mathrm{\Delta }x\hfill \\ \hfill & +L_{g_x}{G}^k\Vert u^{k+1}{\Vert }_1\mathrm{\Delta }t+L_g{G}^{k}TV({\widehat{u}}^{k+1})\mathrm{\Delta }t+\frac{L_{\mathit{\eta }}}{2}{G}^k\Vert u^k{\Vert }_1{\Vert u^{k+1}\Vert }_1\mathrm{\Delta }t\hfill \\ \hfill & +\frac{\mathrm{\Delta }t}{2}{\Vert \mathit{\eta }\Vert }_{\infty }\Vert u^{k+1}{\Vert }_1{\sum }_{j=1}^N|v_j^k|\mathrm{\Delta }x+\frac{\mathrm{\Delta }t}{2}{\Vert \mathit{\eta }\Vert }_{\infty }\Vert {\widehat{u}}^k{\Vert }_1{\sum }_{j=1}^N|v_j^{k+1}|\mathrm{\Delta }x+\mathrm{\Delta }t{L}_R{G}^k{\Vert u^{k+1}\Vert }_1\hfill \\ \hfill & +\mathrm{\Delta }t{L}_{\mathit{\eta }}{G}^k\Vert {\widehat{u}}^{k+1}{\Vert }_1\Vert u^k{\Vert }_1+{\mathrm{\Delta }t\Vert \mathit{\eta }\Vert }_{\infty }\Vert u^{k+1}{\Vert }_1{\sum }_{j=1}^N|v_j^k|\mathrm{\Delta }x+\mathrm{\Delta }t{\widehat{g}}_0^k|v_0^{k+1}|.\hfill \end{array}\end{array}$$

Evaluating the telescoping sum on the left hand side,(A11) $$\begin{array}{c}\hfill \begin{array}{cc}\Vert v^{k+1}{\Vert }_1\hfill & \le \Vert v^k{\Vert }_1+L_{g_x}{B}_2{G}^k\mathrm{\Delta }t+L_g{B}_5{G}^k\mathrm{\Delta }t+\frac{L_{\mathit{\eta }}}{2}{B}_2^2{G}^k\mathrm{\Delta }t+\frac{\mathrm{\Delta }t}{2}{c}_0{B}_2{\Vert v^k\Vert }_1\hfill \\ \hfill & +\frac{\mathrm{\Delta }t}{2}{c}_0{B}_2\Vert v^{k+1}\Vert +\mathrm{\Delta }t{L}_R{G}^k{B}_2\hfill \\ \hfill & +\mathrm{\Delta }t{L}_{\mathit{\eta }}{G}^k{B}_2^2+\mathrm{\Delta }t{c}_0{B}_2\Vert v^k{\Vert }_1+\mathrm{\Delta }t{\widehat{g}}_0^k|v_0^{k+1}|.\hfill \end{array}\end{array}$$(A11)

Note that,$$\begin{array}{c}\hfill \begin{array}{cc}{v}_0^{k+1}\hfill & =u_0^{k+1}-{\widehat{u}}_0^{k+1}=\left(\frac{{\widehat{g}}_0^k-g_0^k}{{\widehat{g}}_0^k{g}_0^k}\right)(F^k+{\sum }_{j=1}^N{\mathit{\beta }}_j^k{u}_j^k\mathrm{\Delta }x)\hfill \\ \hfill & +\frac{1}{{\widehat{g}}_0^k}\left((F^k-{\widehat{F}}^k)+{\sum }_{j=1}^N({\mathit{\beta }}_j^k-{\widehat{\mathit{\beta }}}_j^k)u_j^{k-1}\mathrm{\Delta }x+{\sum }_{j=1}^N{\widehat{\mathit{\beta }}}_j^k(u_j^k-{\widehat{u}}_j^k)\mathrm{\Delta }x\right)\hfill \end{array}\end{array}$$

and the triangle inequality yields,$$\begin{array}{c}\hfill \begin{array}{cc}{\widehat{g}}_0^k|v_0^{k+1}|\hfill & \le \frac{L_g}{\mathit{\zeta }}(c_1+c_2{B}_2+{\Vert \mathit{\beta }\Vert }_{\infty }{B}_2)G^k+(L_F{G}^k+L_{\mathit{\beta }}{B}_2{G}^k+\Vert \widehat{\mathit{\beta }}{\Vert }_{\infty }\Vert v^k{\Vert }_1).\hfill \end{array}\end{array}$$

Define the following constants: $C_{6,1}=L_g(c_1+c_2{B}_2+{\Vert \mathit{\beta }\Vert }_{\infty }{B}_2)/\mathit{\zeta }+L_F+L_{\mathit{\beta }}{B}_2+L_{g_x}{B}_2+L_g{B}_5+3L_{\mathit{\eta }}{B}_2^2/2+L_R{B}_2$, $c_5=3c_0{B}_2/2+{\Vert \mathit{\beta }\Vert }_{\infty }$, $c_4=c_0{B}_2/2$. By (A11(A11) $$\begin{array}{c}\hfill \begin{array}{cc}\Vert v^{k+1}{\Vert }_1\hfill & \le \Vert v^k{\Vert }_1+L_{g_x}{B}_2{G}^k\mathrm{\Delta }t+L_g{B}_5{G}^k\mathrm{\Delta }t+\frac{L_{\mathit{\eta }}}{2}{B}_2^2{G}^k\mathrm{\Delta }t+\frac{\mathrm{\Delta }t}{2}{c}_0{B}_2{\Vert v^k\Vert }_1\hfill \\ \hfill & +\frac{\mathrm{\Delta }t}{2}{c}_0{B}_2\Vert v^{k+1}\Vert +\mathrm{\Delta }t{L}_R{G}^k{B}_2\hfill \\ \hfill & +\mathrm{\Delta }t{L}_{\mathit{\eta }}{G}^k{B}_2^2+\mathrm{\Delta }t{c}_0{B}_2\Vert v^k{\Vert }_1+\mathrm{\Delta }t{\widehat{g}}_0^k|v_0^{k+1}|.\hfill \end{array}\end{array}$$(A11) ) and assumption (B1) we have(A12) $$\begin{array}{c}\hfill \begin{array}{c}(1-c_4\mathrm{\Delta }t)\Vert v^{k+1}{\Vert }_1\le (1+c_5\mathrm{\Delta }t){\Vert v^k\Vert }_1+C_{6,1}{G}^k\mathrm{\Delta }t\hfill \end{array}\end{array}$$(A12)

Observe also that by the fifth row in (4.1(4.1) $$\begin{array}{cc}& \frac{u_j^{k+1}-u_j^k}{\mathrm{\Delta }t}+\frac{g_j^k{u}_j^{k+1}-g_{j-1}^k{u}_{j-1}^{k+1}}{\mathrm{\Delta }x}+u_j^{k+1}\sum _{m=1}^N{\mathit{\eta }}_{m,j}^k{u}_m^k\mathrm{\Delta }x+R_j^k{u}_j^{k+1}\hfill \\ \hfill & =\frac{1}{2}\sum _{s=1}^{j-1}{\mathit{\eta }}_{j-s,s}^k{u}_{j-s}^k{u}_s^{k+1}\mathrm{\Delta }x,\hfill \\ \hfill & g_0^k{u}_0^{k+1}=F^k+\sum _{j=1}^N{\mathit{\beta }}_j^k{u}_j^k\mathrm{\Delta }x,u_j^0=u_0(x_j),\hfill \\ \hfill & \frac{y_{\ell }^{k+1}-y_{\ell }^k}{\mathrm{\Delta }t}=f_{\ell }^k-h_{\ell }^k{y}_{\ell }^{k+1},y_{\ell }^0=y_{\ell ,0},\ell =1,\dots ,m.\hfill \end{array}$$(4.1) ) we have,$$\begin{array}{c}\hfill \begin{array}{cc}|y_{\ell }^{k+1}-{\widehat{y}}_{\ell }^{k+1}|\hfill & \le |f_{\ell }^k-{\widehat{f}}_{\ell }^k|\mathrm{\Delta }t+|h_{\ell }^k-{\widehat{h}}_{\ell }^k{|\Vert y\Vert }_{\infty }\mathrm{\Delta }t+|y_{\ell }^k-{\widehat{y}}_{\ell }^k|\hfill \\ \hfill & \le (L_f+L_h{\Vert y\Vert }_{\infty })(|{\mathit{\sigma }}_{\ell }^k-{\widehat{\mathit{\sigma }}}_{\ell }^k|\mathrm{\Delta }t+{\sum }_{\ell =1}^m|y_{\ell }^k-{\widehat{y}}_{\ell }^k|\mathrm{\Delta }t)+|y_{\ell }^k-{\widehat{y}}_{\ell }^k|.\hfill \end{array}\end{array}$$

Letting $c_7=L_f+L_h{\Vert y\Vert }_{\infty }$, $c_8=m(L_f+L_h{\Vert y\Vert }_{\infty })$ and summing from $\ell =1,\dots ,m$,(A13) $$\begin{array}{c}\hfill \begin{array}{c}{\sum }_{\ell =1}^m|y_{\ell }^{k+1}-{\widehat{y}}_{\ell }^{k+1}|\le c_7{\sum }_{\ell =1}^m|{\mathit{\sigma }}_{\ell }^k-{\widehat{\mathit{\sigma }}}_{\ell }^k|\mathrm{\Delta }t+(1+c_8\mathrm{\Delta }t){\sum }_{\ell =1}^m|y_{\ell }^k-{\widehat{y}}_{\ell }^k|.\hfill \end{array}\end{array}$$(A13)

From (A12(A12) $$\begin{array}{c}\hfill \begin{array}{c}(1-c_4\mathrm{\Delta }t)\Vert v^{k+1}{\Vert }_1\le (1+c_5\mathrm{\Delta }t){\Vert v^k\Vert }_1+C_{6,1}{G}^k\mathrm{\Delta }t\hfill \end{array}\end{array}$$(A12) ) and (A13(A13) $$\begin{array}{c}\hfill \begin{array}{c}{\sum }_{\ell =1}^m|y_{\ell }^{k+1}-{\widehat{y}}_{\ell }^{k+1}|\le c_7{\sum }_{\ell =1}^m|{\mathit{\sigma }}_{\ell }^k-{\widehat{\mathit{\sigma }}}_{\ell }^k|\mathrm{\Delta }t+(1+c_8\mathrm{\Delta }t){\sum }_{\ell =1}^m|y_{\ell }^k-{\widehat{y}}_{\ell }^k|.\hfill \end{array}\end{array}$$(A13) ) we obtain,(A14) $$\begin{array}{c}\hfill \begin{array}{cc}& (1-c_4\mathrm{\Delta }t)(\Vert v^{k+1}{\Vert }_1+{\sum }_{\ell =1}^m|y_{\ell }^{k+1}-{\widehat{y}}_{\ell }^{k+1}|)\le (1+c_5\mathrm{\Delta }t){\Vert v^k\Vert }_1+\mathrm{\Delta }t{C}_{6,1}{G}^k\hfill \\ \hfill & +c_7{\sum }_{\ell =1}^m|{\mathit{\sigma }}_{\ell }^k-{\widehat{\mathit{\sigma }}}_{\ell }^k|\mathrm{\Delta }t+(1+c_8\mathrm{\Delta }t){\sum }_{\ell =1}^m|y_{\ell }^k-{\widehat{y}}_{\ell }^k|.\hfill \end{array}\end{array}$$(A14)

By the triangle inequality and assumption (B8) on $\mathit{\omega }$ and ${\mathit{\gamma }}_{\ell }$, we obtain $G^k\le {\mathrm{\Sigma }}_{\ell =1}^m|y_{\ell }^k-{\widehat{y}}_{\ell }^k{|+\Vert \mathit{\omega }\Vert }_{\infty }{\Vert v^k\Vert }_1$ and that ${\mathrm{\Sigma }}_{\ell =1}^m|{\mathit{\sigma }}_{\ell }^k-{\widehat{\mathit{\sigma }}}_{\ell }^k{|\le m\Vert \mathit{\gamma }\Vert }_{\infty }{\Vert v^k\Vert }_1$. As a result, there exists a constant $C^{†}$ such that(A15) $$\begin{array}{c}\hfill \begin{array}{cc}(1-c_4\mathrm{\Delta }t)(\Vert v^{k+1}{\Vert }_1+{\sum }_{\ell =1}^m|y_{\ell }^{k+1}-{\widehat{y}}_{\ell }^{k+1}|)\hfill & \le (1+C^{†}\mathrm{\Delta }t)(\Vert v^k{\Vert }_1+{\sum }_{\ell =1}^m|y_{\ell }^k-{\widehat{y}}_{\ell }^k|).\hfill \end{array}\end{array}$$(A15)

Applying a standard induction argument on k we arrive at$$\begin{array}{c}\hfill \begin{array}{cc}(\Vert v^{k+1}{\Vert }_1+{\sum }_{\ell =1}^m|y_{\ell }^{k+1}-{\widehat{y}}_{\ell }^{k+1}|)\hfill & \le e^{(c_4+C^{†})T}(\Vert v^0{\Vert }_1+{\sum }_{\ell =1}^m|y_{\ell }^0-{\widehat{y}}_{\ell }^0|).\hfill \end{array}\end{array}$$

Let $B_8=c_4+C^{†}$ to complete the proof.$\square$