Weak and Classical Solutions to Multispecies Advection–Dispersion Equations in Multilayer Porous Media

Abstract

The basic model motivating this work is that of contaminant transport in the Earth’s subsurface, which contains layers in which analytical and semi-analytical solutions of the corresponding advection–dispersion equations could be derived. Then, using the interface relations between adjacent layers, one can streamline the study of the model to the solution to the initial boundary value problem for a coupled parabolic system on partitioned domains. For IBVPs, we set up weak formulations and prove the existence and uniqueness of solutions to appropriate Sobolev-like spaces. A priori estimates at different levels of input data smoothness were obtained. The nonnegativity preservation over time of the solution is discussed. We numerically demonstrate how to solve the reduced truncated problem instead of the original multispecies one with a large number of layers.

1. Introduction

Multispecies energy or mass transfer, along with contaminant transport within the Earth’s subsurface, are fundamental for many biological, chemical, environmental, and industrial heat-mass transfer processes; see, e.g., Calabro and Zunin [1], Carr [2,3], Carslaw [4], Chen et al. [5], Datta [6], Kacur et al. [7], Leij et al. [8], Straughan [9]. In Calabro and Zunin’s work [1], a problem is studied which is governed by linear parabolic PDEs, set in two disjoint domains, and coupled by Neumann–Dirichlet conditions. The theoretical results obtained were applied to investigate the behavior of a biological model for the transfer of chemicals through thin membranes.

In [10], Johnston and Simpsonare derived exact solutions for density profiles, survival probabilities, and splitting probabilities for diffusive processes on multiple growing domains. The exact solutions were obtained via a combination of transformation, convolution, and superposition techniques.

Advection–dispersion equation interface models are often used to predict such phenomena. The corresponding flux has two components: a diffusive one and a convective one. Here, we concentrate on diffusion in a one-dimensional domain with layers.

In [11], Govoli considered a PDE-represented mechanical system with artificial boundaries or an internal interface. With the exact representation of the solution of some layers, the original problem is reduced to a simpler one, in which nonlocal interfaces connect the solutions of partitioned domains.

In [12], Fitzgibbon et al. established the global existence for reaction–diffusion systems on multiple domains and discussed applications to problems of mathematical biology.

In this paper, we use the main ideas of some other papers, see, e.g., Leij et al. [8], Givoli [11], Jovanovic and Vulkov [13], Soliman and Fahmy [14], to replace the layers in which processes are well-known by a fictitious interface, namely a point, line, or surface (with zero thickness) in a 1D, 2D, or 3D analysis, respectively. Special jump conditions, instead of an interface to preserve the mathematical character of the layer, were imposed. In this present work, we impose nonlocal Robin-type jump conditions relating to the values of the solutions and their derivatives at the ends of two separate intervals. Our aim is to extend this idea to models that are similar to those outlined in recent papers by Carr [2,3] and Carslaw [4].

Regarding the multispecies model introduced in Section 2, some (or all of the) layers of the process are described by a system of differential equations; in contrast, in related papers, such as by Leij et al. [8], Givoli [11], Jovanovic and Vulkov [13], Koleva and Vulkov [15], Soliman and Fahmy [14], on each layer, only one equation is considered. These newly established results generalize or complement some earlier results studied in the papers, [8,11,13,14,15].

In [16], Givoli proposed an exact representation of the solution on the interface of the mechanical system. The advantages of this representation in comparison with the approximate approach were discussed as well.

Analytical solutions of the one-dimensional advection–dispersion equation for transport in a two-layered soil profile were discussed in Leij et al. [8]. For semi-infinite layers, solutions were obtained at the inlet boundary and the interface of the two layers.

Jovanovic and Vulkov established the well-posedness of the initial boundary value problem for parabolic equations on disjoint intervals in [13]. The energy stability of the weak solution was also discussed.

In [14], Soliman and Fahmy applied the finite element method to obtain approximate boundary conditions at the interface between a free fluid and a porous medium.

This paper is outlined as follows. In Section 2, we discuss the multispecies, multilayer model, as described by Carr [2]. In Section 3, we use a two-species three-layered interface problem to explain the motivation behind the new disjoint problem formulation. We study the well-posedness of the new one-dimensional problem in Section 4; the nonnegativity of the solution is reviewed in Section 5. The numerical solution to the new two-layered problem is discussed in Section 6. Finally, the paper concludes with some final remarks.

2. A Motivated Multispecies Model

As our model problem, we focus on the reaction–diffusion system [2,3]

$$\frac{\partial {\mathbf{C}}_i}{\partial t}=D_i\frac{{\partial }^2{\mathbf{C}}_i}{\partial x^2}+{\mathbf{M}}^i{\mathbf{C}}_i+{\mathbf{f}}_i\mathrm{for}x\in (l_{i-1},l_i),i=1,\cdots ,I$$

(1)

with initial conditions

$${\mathbf{C}}_i(x,0)={\mathsf{\ae }}_i\left(x\right)\mathrm{for}x\in (l_{i-1},l_i),i=1,\cdots ,I,$$

(2)

interface conditions

$${\mathbf{C}}_i(l_i,t)={\mathbf{C}}_{i+1}(l_i,t),i=1,\cdots ,I-1,$$

(3)

$${\theta }_i{D}_i\frac{\partial {\mathbf{C}}_i}{\partial x}(l_i,t)={\Theta }_{i+1}{D}_{i+1}\mathbf{C}(l_i,t),i=1,\cdots ,I-1,$$

(4)

and boundary conditions

$$-a_0\frac{\partial {\mathbf{C}}_1}{\partial x}(0,t)+b_0{\mathbf{C}}_1(0,t)={\mathbf{g}}_0\left(t\right),$$

(5)

$$a_L\frac{\partial {\mathbf{C}}_I}{\partial x}(L,t)+b_L(L,t){\mathbf{C}}_I={\mathbf{g}}_I\left(t\right),$$

(6)

Here the vectors and matrices are defined as follows:

$${\mathbf{C}}_i(x,t)=\left[\begin{array}{c}{C}_{i1}(x,t)\hfill \\ C_{i2}(x,t)\hfill \\ ⋮\hfill \\ C_{iN}(x,t)\hfill \end{array}\right],{\mathbf{M}}^i(x,t)=\left[\begin{array}{cccc}{\mu }_{11}^i& {\mu }_{12}^i& \cdots & {\mu }_{1N}^i\\ {\mu }_{21}^i& {\mu }_{22}& \cdots & {\mu }_{2N}^i\\ ⋮& ⋮& \ddots & \\ {\mu }_{N1}^i& {\mu }_{N2}^i& \cdots & {\mu }_{NN}^i\end{array}\right],$$

$${\mathbf{f}}_i(x,t)=\left[\begin{array}{c}{f}_{i1}\hfill \\ f_{i2}\hfill \\ ⋮\hfill \\ f_{iN}\hfill \end{array}\right],{\mathsf{\ae }}_i\left(x\right)=\left[\begin{array}{c}{\rho }_{i1}\hfill \\ {\rho }_{i2}\hfill \\ ⋮\hfill \\ {\rho }_{in}\hfill \end{array}\right],{\mathbf{g}}_0\left(t\right)=\left[\begin{array}{c}{g}_{01}\hfill \\ g_{02}\hfill \\ ⋮\hfill \\ g_{0N}\hfill \end{array}\right],{\mathbf{g}}_I\left(t\right)=\left[\begin{array}{c}{g}_{I1}\hfill \\ g_{I2}\hfill \\ ⋮\hfill \\ g_{NN}\hfill \end{array}\right],$$

In [2,3], the unknown function ${\mathbf{C}}_{in}(x,t)$ is the concentration $\left[M{L}^{-3}\right]$ of species n in the i-th layer, where $x\in (l_{i-1},l_i)$ and $t\ge 0$ are the spatial and time variables.

It is assumed that the medium is layered, as follows:

$$0=l_0<l_1<\cdots <l_{I-1}<l_I=L.$$

(7)

and let N be the number of species. Here, $R_i=[R_{i1},R_{i2},\cdots ,R_{iN}]$, $R_{ij}>0$ is the retardation function in layer i for species j. Next, $D_i$ is the dispersion coefficient $\left[L^2{T}^{-1}\right]$ in layer i, ${\theta }_i$ is the volume water content in layer i, inherent to Equation (4), and is the assumption of flow continuity ${\theta }_i{v}_I={\theta }_{i+1}{v}_{i+1}$ for all $i=1,\cdots ,m-1$ [8], $f_{in}$ is the rate of zero-order production $\left[M{L}^{-1}{T}^{-1}\right]$ of species n in layer i [2,3]. E.J. Carr [2] constructed a solution for arbitrary numbers (N and I) of species and layers, respectively, to Problem (1)–(6). The reaction between species is governed by the value ${\mu }_{nk}\left[T^{-1}\right]$.

The theoretical studying of parabolic systems (even linear) requires, in general, the additional assumption that the system has a (quasi-) monotonicity property, also called cooperativeness (this term comes from biology, where models from the population dynamics of species that cooperate with each other lead to quasi-monotone systems; see, e.g., [17]).

In our case, for every $i=1,\cdots ,I$ on the layer $(l_{i-1},l_i)$, system (1) is cooperative; see, e.g., [18].

The condition of being cooperative for all indices $i,j\in \{1,2,\cdots ,N\}$, $i\ne j$ implies that ${\mu }_{i,j}\ge 0$, a.e., in $(l_{i-1},l_i)$ (7). This condition is rather restrictive, but many important systems in applied sciences do satisfy it (for instance, in biology) [19].

3. The Two-Species, Three-Layer Model Problem and Its Reduction

In this section, we show how, in a natural way, a new specific interface problem on a disjoint interval from a model of types (1)–(6) can arise. For the sake of clarity, we consider a three-layer advection–dispersion reaction model with two species of each layer.

New analytical solutions to advection–dispersion systems of Equations (1)–(6), which model solute transport, were derived by E.J. Carr in [2]. The transport coefficient (dispersion, coefficient, retardation factor, etc.) in each layer are functions of space and time. When following up on this work (and other ones where the coefficients are space-dependent, for example, [3,10,20,21,22]), we will consider the next representative to model (1)–(6). Namely, if we let $u_1=(u_{11},u_{12})$, $u_3=(u_{31},u_{32})$, then we have

$$\begin{array}{ccc}& & L_{11}{u}_1:=\frac{\partial u_{11}}{\partial t}-\frac{\partial }{\partial x}\left(p_{11}\left(x\right)\frac{\partial u_{11}}{\partial x}\right)-{\mu }_{11}^1{u}_{11}-{\mu }_{12}^1{u}_{12}=f_{11}(x,t),\hfill \end{array}$$

(8)

$$\begin{array}{ccc}& & L_{12}{u}_1:=\frac{\partial u_{12}}{\partial t}-\frac{\partial }{\partial x}\left(p_{12}\left(x\right)\frac{\partial u_{12}}{\partial x}\right)-{\mu }_{21}^1{u}_{11}-{\mu }_{22}^1{u}_{12}=f_{12}(x,t),\hfill \end{array}$$

(9)

$$\begin{array}{ccc}& & (x,t)\in Q_{1T}={\Omega }_1\times (0,T),{\Omega }_1=(l_0,l_1);\hfill \end{array}$$

$$\begin{array}{ccc}& & R_{21}\frac{\partial u_{21}}{\partial t}=\frac{\partial }{\partial x}\left(p_{21}\left(x\right)\frac{\partial u_{21}}{\partial x}\right)+{\mu }_1^2{u}_{21}+f_{21}(x,t),\hfill \end{array}$$

(10)

$$\begin{array}{ccc}& & L_{22}{u}_{22}:=R_{22}\frac{\partial u_{22}}{\partial t}-\frac{\partial }{\partial x}\left(p_{22}\left(x\right)\frac{\partial u_{22}}{\partial x}\right)-{\mu }_2^2{u}_{22}=f_{22}(x,t),\hfill \end{array}$$

(11)

$$\begin{array}{ccc}& & (x,t)\in Q_{2T}={\Omega }_2\times (0,T),{\Omega }_2=(l_1,l_2);\hfill \end{array}$$

$$\begin{array}{ccc}& & L_{31}{u}_3:=\frac{\partial u_{31}}{\partial t}-\frac{\partial }{\partial x}\left(p_{31}\left(x\right)\frac{\partial u_{31}}{\partial x}\right)-{\mu }_{11}^3{u}_{31}-{\mu }_{12}^3{u}_{32}=f_{31}(x,t),\hfill \end{array}$$

(12)

$$\begin{array}{ccc}& & L_{32}{u}_3:=\frac{\partial u_{32}}{\partial t}-\frac{\partial }{\partial x}\left(p_{32}\left(x\right)\frac{\partial u_{32}}{\partial x}\right)-{\mu }_{21}^3{u}_{31}-{\mu }_{22}^3{u}_{32}=f_{32}(x,t),\hfill \end{array}$$

(13)

$$\begin{array}{ccc}& & (x,t)\in Q_{3T}={\Omega }_3\times (0,T),{\Omega }_3=(l_2,l_3).\hfill \end{array}$$

Following the interface relations (3) and (4) of the leading models (1)–(6) of [2], we impose system (8)–(13), and the following interface conditions:

$$\begin{array}{ccc}& & u_{11}(l_1,t)=u_{21}(l_1,t),p_{11}\left(l_1\right)\frac{\partial u_{11}}{\partial x}(l_1,t)=p_{21}\left(l_1\right)\frac{\partial u_{21}}{\partial x}(l_1,t),\hfill \end{array}$$

(14)

$$\begin{array}{ccc}& & u_{12}(l_1,t)=u_{22}(l_1,t),p_{12}\left(l_1\right)\frac{\partial u_{12}}{\partial x}(l_1,t)=p_{22}\left(l_1\right)\frac{\partial u_{22}}{\partial x}(l_1,t)\hfill \end{array}$$

(15)

and

$$\begin{array}{ccc}& & u_{21}(l_2,t)=u_{31}(l_2,t),p_{21}\left(l_2\right)\frac{\partial u_{21}}{\partial x}(l_2,t)=p_{22}\left(l_2\right)\frac{\partial u_{31}}{\partial x}(l_2,t),\hfill \end{array}$$

(16)

$$\begin{array}{ccc}& & u_{22}(l_2,t)=u_{32}(l_2,t),p_{22}\left(l_2\right)\frac{\partial u_{22}}{\partial x}(l_2,t)=p_{32}\left(l_2\right)\frac{\partial u_{32}}{\partial x}(l_2,t).\hfill \end{array}$$

(17)

A schematic representation and physical interpretation of Problems (8)–(17) are given in Figure 1.

Here, we work with classical solutions, which need to be at least twice continuously differentiable in space and once in time, and they satisfy the differential at each point of $Q_T$.

First, we suppose that retardation $R_{21}\cong \epsilon >0$ is very small. Then, Equation (10) reduces to a second-order ordinary differential equation. In layer $(l_1,l_2)$, we solve (10) analytically. Its general solution is in the following form:

$$u_{21}(x,t)={\phi }_1\left(t\right)v_1\left(x\right)+{\phi }_2\left(t\right)v_2\left(x\right)+w(x,t),$$

(18)

where $v_1$, $v_2$ are known functions, and ${\phi }_1\left(t\right)$, ${\phi }_2\left(t\right)$ are unknown functions. Taking $x=l_1$ and $x=l_2$ in (17), we obtain the following relations:

$$\left[\begin{array}{cc}{v}_1\left(l_1\right)& v_2\left(l_1\right)\\ v_1\left(l_2\right)& v_2\left(l_2\right)\end{array}\right]\left[\begin{array}{c}{\phi }_1\left(t\right)\\ {\phi }_2\left(t\right)\end{array}\right]=\left[\begin{array}{c}{u}_{21}(l_1,t)-w(l_1,t)\\ u_{21}(l_2,t)-w(l_2,t)\end{array}\right].$$

(19)

By solving the system of equations, we express ${\phi }_1\left(t\right)$ and ${\phi }_2\left(t\right)$ in terms of $u_{21}(l_1,t)$, $u_{21}(l_2,t)$, and $w(l_1,t)$, $w(l_2,t)$. Let us denote

$$\begin{array}{c}\hfill \begin{array}{cc}& {\alpha }_1=\frac{p_{21}\left(l_1\right){▵}_1(l_2,l_1)}{▵(l_1,l_2)},{\beta }_1=\frac{p_{21}\left(l_1\right){▵}_1(l_1,l_1)}{▵(l_1,l_2)},\hfill \\ & {\alpha }_3=\frac{p_{21}\left(l_2\right){▵}_1(l_1,l_2)}{▵(l_1,l_2)},{\beta }_3=\frac{p_{21}\left(l_2\right){▵}_1(l_2,l_2)}{▵(l_1,l_2)},\hfill \\ & {\gamma }_1\left(t\right)=\frac{p_{21}\left(l_1\right)}{▵(l_1,l_2)}\left({▵}_1(l_2,l_1)w(l_1,t)-{▵}_1(l_1,l_1)w(l_2,t)+▵(l_1,l_2)\frac{\partial w}{\partial x}(l_1,t)\right),\hfill \\ & {\gamma }_3\left(t\right)=\frac{p_{21}\left(l_2\right)}{▵(l_1,l_2)}\left({▵}_1(l_2,l_2)w(l_1,t)-{▵}_1(l_1,l_2)w(l_2,t)+▵(l_1,l_2)\frac{\partial w}{\partial x}(l_2,t)\right),\hfill \end{array}\end{array}$$

(20)

where

$$\begin{array}{ccc}& & ▵(l_1,l_2)=\left|\begin{array}{cc}{v}_1\left(l_1\right)& v_2\left(l_1\right)\\ v_1\left(l_2\right)& v_2\left(l_2\right)\end{array}\right|=v_1\left(l_1\right)v_2\left(l_2\right)-v_1\left(l_2\right)v_2\left(l_1\right),\hfill \\ & & {▵}_1(r,s)=\left|\begin{array}{cc}{v}_1\left(r\right)& v_2\left(r\right)\\ {\displaystyle \frac{d{v}_1}{dx}\left(s\right)}& {\displaystyle \frac{d{v}_2}{dx}\left(s\right)}\end{array}\right|=v_1\left(r\right)\frac{d{v}_2}{dx}\left(s\right)-v_2\left(r\right)\frac{d{v}_1}{dx}\left(s\right).\hfill \end{array}$$

Now, we are in a position to prove the following assertion:

Theorem 1.

Let

$$p_{21}\left(x\right)\ge p_{21}^0=const.>0,\mathrm{a}.\mathrm{e}.in{\Omega }_2\mathit{and}{p}_{21}\in L_{\infty }\left({\Omega }_2\right).$$

(21)

Then,

$${\alpha }_i>0,{\beta }_i>0,i=1,3$$

(22)

and the nonlocal Robins–Dirichlet interface jump conditions are valid

$$\begin{array}{ccc}& & l_1(u_{11},u_{31}):=p_{11}\left(l_1\right)\frac{\partial u_{11}}{\partial x}(l_1,t)+{\alpha }_1{u}_{11}(l_1,t)-{\beta }_1{u}_{31}(l_2,t)={\gamma }_1\left(t\right),\hfill \end{array}$$

(23)

$$\begin{array}{ccc}& & l_2(u_{11},u_{31}):=-p_{31}\left(l_2\right)\frac{\partial u_{31}}{\partial x}(l_2,t)+{\alpha }_3{u}_{31}(l_2,t)-{\beta }_3{u}_{11}(l_1,t)={\gamma }_3\left(t\right).\hfill \end{array}$$

(24)

Proof. 

Functions $v_1\left(t\right)$ and $v_2\left(t\right)$ are two linearly independent solutions of the corresponding (10) homogeneous ordinary differential equation:

$$\mathcal{L}v=-\frac{d}{dx}\left(p_{21}\left(x\right)\frac{dv}{dx}\right)+{\mu }_1^{2}v=0,x\in {\Omega }_2.$$

(25)

Indeed, we can take $v_1\left(t\right)$ and $v_2\left(t\right)$ as the solutions of the Cauchy problems:

$$\begin{array}{ccc}& & \mathcal{L}{v}_1=0,x\in {\Omega }_2,v_1\left(l_1\right)=0,p_{21}\left(l_1\right)\frac{d{v}_1}{dx}\left(l_1\right)=1,\hfill \end{array}$$

(26)

$$\begin{array}{ccc}& & \mathcal{L}{v}_2=0,x\in {\Omega }_2,v_2\left(l_2\right)=0,p_{21}\left(l_2\right)\frac{d{v}_2}{dx}\left(l_2\right)=-1.\hfill \end{array}$$

(27)

A maximum principle application [18] results in the following:

$$v_1\left(l_2\right)>0,\frac{d{v}_1}{dx}\left(l_1\right)>0,\frac{d{v}_1}{dx}\left(l_2\right)>0,$$

and

$$v_2\left(l_1\right)>0,\frac{d{v}_2}{dx}\left(l_1\right)<0,\frac{d{v}_2}{dx}\left(l_2\right)<0,$$

which provide (22).

Next, we place Formula (18) in (14) and (16), and after a suitable rearrangement, using the above-obtained formulas for ${\alpha }_i$, ${\beta }_i$, ${\gamma }_i$, $i=1,3$ and $▵$, ${▵}_1$, we obtain the interface conditions (23) and (24). □

Let us suppose that boundary conditions of types (5) and (6) are imposed for the unknown functions $u_{11}(x,t)$, $u_{12}(x,t)$, and $u_{31}(x,t)$, $u_{32}(x,t)$:

$$\begin{array}{ccc}& & l_{0i}\left(u_{1i}\right):=-a_{0i}\frac{\partial u_{1i}}{\partial x}(l_0,t)+b_{0i}{u}_{1i}(l_0,t)=g_{0i}\left(t\right),i=1,2,t\in (0,T),\hfill \end{array}$$

(28)

$$\begin{array}{ccc}& & l_{3i}\left(u_{3i}\right):=a_{3i}\frac{\partial u_{3i}}{\partial x}(l_3,t)+b_{3i}{u}_{3i}(l_3,t)=g_{3i}\left(t\right),i=1,2,t\in (0,T).\hfill \end{array}$$

(29)

Also, assume that the initial conditions are imposed

$$u_{1i}(x,0)={\psi }_{1i}\left(x\right),u_{3i}(x,0)={\psi }_{3i}\left(x\right),i=1,2,u_{22}(x,0)={\psi }_{22}\left(x\right),$$

(30)

Now, we are in a position to formulate the reduced truncated problem, obtained by using the possibility to construct the exact solution of Equation (10). Namely, we seek the unknown functions $u_{1i}(x,t)$, $i=1,2$, $u_{22}(x,t)$, and the pair$u_{3i}(x,t)$, $i=1,2$, satisfied the parabolic Equations (9) and (11)–(13), the interface conditions (14)–(17), initial and boundary conditions (28)–(30).

The demonstrated technique of the representation of the solution $u_{21}(x,t)$ can be implemented for important problems and defined on a large number of layers of applied mechanics. Thus, the new computational domain can be significantly smaller than the original one.

In Figure 2, we schematically illustrate the reduced truncated problem.

Remark 1.

In a similar way, we can consider the case where one can find the exact solution to Equation (11) and obtain a new IBVP for the unknown functions $u_{1i}$, $i=1,2$, $u_{21}$, and $u_{3i}$, $i=1,2$ with appropriate initial, boundary, and interface conditions, together with new nonlocal Robins–Dirichlet conditions of types (23) and (24) for functions $u_{22}$ and $u_{31}$.

It is also possible to find the exact solutions of equations for $u_{2i}$, $i=1,2$, even for nonlinear cases [17]. Then, the reduced IBVP finds the unknown functions $u_{1i}$, $u_{1i}$, and $i=1,2$ with the appropriate initial, boundary, and interface conditions, and new nonlocal Robins–Dirichlet conditions of types (23) and (24) for functions $u_{12}$ and $u_{31}$.

Remark 2.

The approach developed could be implemented for problems with other types of interface conditions instead of (3) and (4). For example,

(i)

Type 1:

$$\begin{array}{ccc}& & {\mathbf{C}}_i(l_i,t)={\theta }_i{\mathbf{C}}_{i+1}(l_i,t),{\theta }_i=const.>0,i=1,2,\cdots ,I-1,\hfill \\ & & D_i\frac{\partial {\mathbf{C}}_i}{\partial x}(l_1,t)=D_{i+1}{\mathbf{C}}_{i+1}(l_i,t),i=1,2,\cdots ,I-1.\hfill \end{array}$$

(ii)

Type 2:

$$\begin{array}{ccc}& & D_i\frac{\partial {\mathbf{C}}_i}{\partial x}(l_1,t)={\theta }_i({\mathbf{C}}_{i+1}(l_i,t)-{\mathbf{C}}_i(l_i,t)),{\theta }_i=const.>0,i=1,2,\cdots ,I-1,\hfill \\ & & D_{i+1}\frac{\partial {\mathbf{C}}_i}{\partial x}(l_1,t)={\theta }_i({\mathbf{C}}_{i+1}(l_i,t)-{\mathbf{C}}_i(l_i,t)),i=1,2,\cdots ,I-1,\hfill \end{array}$$

4. Existence and Uniqueness of Solutions

A weak solution needs to be only weakly differentiable while the classical solution must be at least continuously differentiable. The intermediate solutions concern the so-called strong solutions.

Now, we can study the well-posedness of the problem formulated in the previous section. Thus, we consider the following: the existence and uniqueness of weak, strong, and classical solutions; the smoothness and positivity of the solutions; energy stability, etc. In this section, we prove the existence and uniqueness of the weak solution in Sobolev spaces. For the concept of strong solutions, the readers can read, for example, the references [18,23]. Strong and classical solutions are not the same. For parabolic problems, a strong solution is two times weakly differentiable in space, one in time, and satisfies the differential problem almost everywhere. The use of strong solutions instead of weak solutions can increase the rate of convergence of the difference scheme; see, e.g., [24].

Existence and Uniqueness of Weak Solutions

We introduce the product space

$$\begin{array}{ccc}\hfill L& =& L_2\left({\Omega }_1\right)\times L_2\left({\Omega }_1\right)\times L_2\left({\Omega }_2\right)\times L_2\left({\Omega }_3\right)\times L_2\left({\Omega }_3\right)\hfill \\ & \equiv & \{v=(v_1,v_2,v_3,v_4,v_5)|v_i\in L_2\left({\Omega }_i\right),i=1,2,v_3\in L_2\left({\Omega }_2\right),v_4,v_5\in L_2\left({\Omega }_3\right)\},\hfill \end{array}$$

endowed with the inner product and the associated norm

$$\begin{array}{ccc}\hfill {(u,v)}_L& =& {\beta }_3{(u_{11},v_1)}_{L_2\left({\Omega }_1\right)}+{(u_{12},v_2)}_{L_2\left({\Omega }_1\right)}+{(u_{22},v_3)}_{L_2\left({\Omega }_2\right)}\hfill \\ & & +{\beta }_1{(u_{31},v_4)}_{L_2\left({\Omega }_3\right)}+{(u_{32},v_5)}_{L_2\left({\Omega }_3\right)},\hfill \end{array}$$

where

$${(u_{11},v_1)}_{L_2\left({\Omega }_1\right)}={\int }_{{\Omega }_1}{u}_{11}{v}_{1}dx,\mathrm{etc.},u\equiv (u_{11},u_{12},u_{22},u_{31},u_{32}).$$

We define the space

$$\begin{array}{ccc}\hfill H^1& =& \{v=(v_1,v_2,v_3,v_4,v_5)|v_1\in H^1\left({\Omega }_1\right)\mathrm{and}{v}_1\left(l_0\right)=0,v_2\in H^1\left({\Omega }_1\right),\hfill \\ & & v_3\in H^1\left({\Omega }_2\right),v_4\in H^1\left({\Omega }_3\right),v_5\in H^1\left({\Omega }_3\right)\mathrm{and}{v}_5\left(l_3\right)=0\},\hfill \end{array}$$

endowed with the inner product and norm

$$\begin{array}{ccc}\hfill {(u,v)}_{H^1}& =& {\beta }_3{(u_{11},v_1)}_{H^1\left({\Omega }_1\right)}+{(u_{12},v_2)}_{H^1\left({\Omega }_1\right)}+{(u_{22},v_3)}_{H^1\left({\Omega }_2\right)}\hfill \\ & & +{\beta }_1{(u_{31},v_4)}_{H^1\left({\Omega }_3\right)}+{(u_{32},v_5)}_{H^1\left({\Omega }_3\right)}.\hfill \end{array}$$

Finally, with $u=(u_{11},u_{12},u_{22},u_{31},u_{32})$ and $v=(v_1,v_2,v_3,v_4,v_5)$ we define the following bilinear form

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill A(u,v)=& {\beta }_3{\int }_{{\Omega }_1}\left(p_{11}\frac{d{u}_{11}}{dx}\frac{d{v}_1}{dx}\right)dx+{\int }_{{\Omega }_2}{p}_{12}\frac{\partial u_{12}}{\partial x}.\frac{\partial v_2}{\partial x}dx+{\int }_{{\Omega }_2}{p}_{22}\frac{d{u}_{22}}{\partial x}\frac{d{v}_3}{dx}dx\hfill \\ \hfill & +{\beta }_1\int \left(p_{31}\frac{d{u}_{31}}{dx}\frac{d{v}_4}{dx}\right)dx+{\int }_{{\Omega }_3}{p}_{32}\frac{d{u}_{32}}{dx}\frac{d{v}_5}{dx}dx\hfill \\ \hfill & +{\beta }_3{\alpha }_1{u}_{11}\left(l_1\right)v_1\left(l_1\right)+{\beta }_1{\alpha }_3{u}_{31}\left(l_2\right)v_4\left(l_2\right)-{\beta }_1{\beta }_3[u_{11}\left(l_1\right)v_4\left(l_2\right)+u_{11}\left(l_1\right)v_4\left(l_2\right)]\hfill \\ \hfill & +I_1+I_2+I_3,\hfill \end{array}\end{array}$$

(31)

where

$$I_1={\int }_{{\Omega }_1}({\beta }_3{\mu }_{11}^1{u}_{11}^2+{\mu }_{22}^1{u}_{12}^2)dx+{\int }_{{\Omega }_1}({\beta }_3{\mu }_{12}^1+{\mu }_{21}^1)u_{11}{u}_{12}dx,$$

$$I_2={\int }_{{\Omega }_2}{\mu }_2^2{u}_{22}^{2}dx,I_3={\int }_{{\Omega }_3}({\beta }_1{\mu }_{12}^3+{\mu }_{21}^3)u_{31}{u}_{32}dx.$$

Next, using the well-known substitution $w(x,t):=e^{\rho t}w(x,t)$ where $\rho >0$ is a sufficiently large constant and w is the solution $u_{1j},j=1,2$ or $u_3,j=1,2$, we can assume that the constants ${\mu }_{ii}^1,i=1,2$ and ${\mu }_{ii}^3,i=1,2$ are negative and sufficiently large in the absolute value when compared to ${\mu }_{ij}^1,{\mu }_{ij}^3,i\ne j.$

Lemma 1.

Suppose that

$$p_{1i}\left(x\right)\in L_{\infty }\left({\Omega }_1\right),i=1,2;p_{22}\left(x\right)\in L_{\infty }\left({\Omega }_2\right),p_{3i}\left(x\right)\in L_{\infty }\left({\Omega }_3\right),i=1,2$$

(32)

and there exist positive constants $c_{1i},i=1,2$, $c_2$ and $c_{3i},i=1,2$ such that

$$0<c_{1i}\le p_{1i}\left(x\right),i=1,2,x\in {\Omega }_1;and{c}_{3i}\le p_{3i}\left(x\right),i=1,2,x\in {\Omega }_3$$

(33)

In addition to (22), we also assume

$${\beta }_1{\beta }_3\le {\alpha }_1{\alpha }_3.$$

(34)

Then, under conditions (22) and (32)–(34), the bilinear form A, defined by (31), is symmetric and bounded on $H^1\times H^1$. Moreover, this form is also coercive, i.e., there exists a constant $C_0>0$, such that

$$A(u,u)\ge C_0{\parallel u\parallel }_{H^1}^2\forall H^1.$$

Proof. 

We multiply Equation (8) by ${\beta }_3{u}_{11}$ and (9) by $u_{12}$; Equation (10) by $u_{21}$ and (11) by $u_{22}$; Equation (12) by ${\beta }_1{u}_{31}$, and (13) by $u_{32}$.

Then, we integrate by parts the results on the corresponding space intervals, taking into account the initial and boundary conditions, to obtain (31).

The symmetry of A is obvious, while the boundedness follows from (32) and the embeddings $H^1\left({\Omega }_i\right)\subset C\left({\overline{\Omega }}_i\right)$, $i=1,2,3$.

Under Conditions (22) and (34), we have

$${\beta }_3{\alpha }_1{u}_{11}^2\left(l_1\right)+{\beta }_1{\alpha }_3{u}_{31}^2\left(l_2\right)-2{\beta }_1{\beta }_3{u}_{11}\left(l_1\right)u_{31}\left(l_2\right)\ge 0$$

Next, for the estimation of the integrals $I_1$, $I_2$, and $I_3$, we can use the following trick. By the well-known substitution $w(x,t):=e^{\rho t}w(x,t)$, where $\rho >0$ is a sufficiently large constant and w is the solution $u_{ij},j=1,2$ or $u_{3j},j=1,2$, we can assume that ${\mu }_{ii}^1,i=1,2$ and ${\mu }_{ii}^3,i=1,2$ are negative and sufficiently large in absolute value in comparison with ${\mu }_{ij}^1>0,{\mu }_{ij}^3,i\ne j,i,j=1,2$. Then, with the simple inequality $2ab\le a^2+b^2$ and together with (33) and Poincare-type inequalities,

$${\int }_a^b{w}^2\left(x\right)dx\le \frac{b^2-a^2}{2}{\int }_a^b{\left(\frac{du}{dx}\right)}^{2}dx$$

provides the coercive of operator A. □

We use the following notations. Let $\Omega$ be a bounded domain in $R^n$ and $v\left(t\right)$ a function mapping into a Hilbert space $H>0$. In a standard manner, see, e.g., [23,25], we define the Sobolev space of vector-valued functions $H^k(\Omega ,H)$, endowed with the inner product

$${(v,w)}_{H^k(\Omega ,H)}={\int }_{\Omega }\sum _{\left|\alpha \right|\le k}{(D^{\alpha }v\left(t\right),D^{\alpha }w\left(t\right))}_{H}dt,k=0,1,2,dots$$

and with the usual modification for non-integer k.

We set $L_2(\Omega ,H)=H^0(\Omega )$ and define $L^2((0,T),H^1)\cap H^{1/2}((0,T),L)$.

Let $H^{-1}={\left(H^1\right)}^{*}$ be the dual space for $H^1$. The spaces $H^1,L$ and $H^{-1}$ form a Gelfand triple $H^1\subset L\subset H^{-1}$ (see [23,25]), with continuous and dense embedding. We also introduce the space

$$W(0,T)=\left\{v\right|v\in L_2((0,T),H^1),\frac{\partial v}{\partial t}\in L_2((0,T),H^1)\}$$

with the inner product

$${(v,w)}_{W(0,T)}={\int }_0^T\left[{(v(·,t),w(·,t))}_{H^1}dt+{\left(\frac{\partial v}{\partial t}(·,t),\frac{\partial w}{\partial t}(·,t)\right)}_{H^{-1}}\right]dt.$$

Now, we are in the position to introduce a weak solution of the truncated problem.

Assume for simplicity that

$${\gamma }_i\left(t\right)=0,i=1,3,a_{0i}\left(t\right)=a_{3i}\left(t\right)=g_{0i}\left(t\right)=g_{3i}\left(t\right)=0,i=1,2.$$

(35)

Multiplying Equation (8) by a function $v_{11}(x,t)$, such that $v_{11}(l_0,t)=0$ and integrating the result by parts on ${\Omega }_1$, using Condition (23), we have the following:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {\left(\frac{\partial u_{11}}{\partial t}(·,t),v_{11}(·,t)\right)}_{L_2\left({\Omega }_1\right)}+{\int }_{{\Omega }_1}{p}_{11}\frac{\partial u_{11}}{\partial x}\frac{\partial v_{11}}{\partial x}dx+{\int }_{{\Omega }_1}({\mu }_{11}^1{u}_{11}{v}_{11}+{\mu }_{12}^1{u}_{12}{v}_{12})dx\\ \hfill +{\alpha }_1{u}_{11}(l_1,t)v_{11}(l_1,t)-{\beta }_1{u}_{31}(l_2,t)v_{11}(l_1,t)=(f_{11}(·,t),v_1(·,t)).\end{array}\end{array}$$

(36)

In a similar way, multiplying Equation (12) by $v_{31}(x,t)$ such that $v_{31}(l_3,t)=0$ and integrating by parts on ${\Omega }_3$, taking into account Condition (24), we obtain

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {\left(\frac{\partial u_{31}}{\partial t}(·,t),v_{31}(·,t)\right)}_{L_2\left({\Omega }_3\right)}+{\int }_{{\Omega }_3}{p}_{31}\frac{\partial u_{31}}{\partial x}\frac{\partial v_{31}}{\partial x}dx+{\int }_{{\Omega }_3}({\mu }_{11}^3{u}_{31}{v}_{31}+{\mu }_{12}^3{u}_{32}{v}_{32})dx\\ \hfill +{\alpha }_3{u}_{31}(l_2,t)v_{31}(l_2,t)-{\beta }_3{u}_{11}(l_1,t)v_{32}(l_2,t)=(f_{31}(·,t),v_{32}(·,t)).\end{array}\end{array}$$

(37)

Further, we multiply Equation (9) by $v_{12}(x,t)$, such that $v_{12}(l_0,t)=0$ and integrating by parts on ${\Omega }_1$ to obtain

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {\left(\frac{\partial u_{12}}{\partial t}(·,t),v_{12}(·,t)\right)}_{L_2\left({\Omega }_1\right)}+{\int }_{{\Omega }_1}{p}_{12}\frac{\partial u_{12}}{\partial x}\frac{\partial v_{12}}{\partial x}dx{\int }_{{\Omega }_1}({\mu }_{21}^1{u}_{21}{v}_{12}+{\mu }_{22}^1{u}_{22}{v}_{12})dx\\ \hfill +p_{12}\left(l_1\right)\frac{\partial u_{12}}{\partial x}(l_1,t)v_{12}(l_1,t)={(f_{12}(·,t),v_{12}(·,t))}_{L_2\left({\Omega }_1\right)}.\end{array}\end{array}$$

(38)

In a similar way, we multiply Equation (11) by $v_{22}(x,t)$, such that $v_{22}(l_1,t)=v_{12}(l_1,t)$ and integrating by parts on ${\Omega }_2$, to derive

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {\left(\frac{\partial u_{22}}{\partial t}(·,t),v_{22}(·,t)\right)}_{L_2\left({\Omega }_2\right)}+{\int }_{{\Omega }_2}{p}_{22}\frac{\partial u_{22}}{\partial x}\frac{\partial v_{22}}{\partial x}dx+p_{22}\left(l_3\right)\frac{\partial u_{22}}{\partial x}(l_3,t)v_{22}(l_3,t)\\ \hfill -p_{22}\left(l_2\right)\frac{\partial u_{22}}{\partial x}(l_2,t)v_{22}(l_2,t)+{\int }_{{\Omega }_2}{\mu }_2^2{u}_{22}{v}_{22}dx={(f_{22}(·,t),v_{22}(·,t))}_{L_2\left({\Omega }_2\right)}.\end{array}\end{array}$$

(39)

Finally, multiplying Equation (15) by function $v_{32}(x,t)$, such that $v_{32}(l_2,t)=v_{22}(l,t)$, $v_{32}(l_3,t)=0$ and integrating by parts on ${\Omega }_3$, we have

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {\left(\frac{\partial u_{32}}{\partial t}(·,t),v_{32}(·,t)\right)}_{L_2\left({\Omega }_3\right)}+{\int }_{{\Omega }_3}{p}_{32}\frac{\partial u_{32}}{\partial x}\frac{\partial v_{32}}{\partial x}dx+p_{32}\left(l_3\right)\frac{\partial u_{32}}{\partial x}(l_3,t)v_{32}(l_3,t)\\ \hfill -p_{32}\left(l_2\right)\frac{\partial u_{32}}{\partial x}(l_2,t)v_{32}(l_2,t)+{\int }_{{\Omega }_3}(p_{21}^3{u}_{31}{v}_{32}+p_{22}^3{u}_{32}{v}_{32})dx\\ \hfill ={(f_{32}(·,t),v_{32}(·,t))}_{L_2\left({\Omega }_3\right)}.\end{array}\end{array}$$

(40)

Now, we multiply identity (36) by ${\beta }_3$ (37) by ${\beta }_1$ and we sum the results. Next, we sum up the identities (38)–(40), taking into account the interface conditions (16) and (17). Finally, we sum up the last two results to obtain the weak formulation of the truncated problems (8), (9), and (11)–(13):

$${\left(\frac{\partial u}{\partial t}(·,t),v(·,t)\right)}_L+A(u(·,t),v(·,t))={(f(·,t),v(·,t))}_L$$

(41)

An application of Theorem 26.1 from [14] to (41) provides the existence and uniqueness of the weak solution to the reduced (truncated) problem.

Theorem 2.

Let us assume that (21), (22), and (34) hold, and

$$\begin{array}{ccc}& & u_0\equiv u(x,0)=({\psi }_{11},{\psi }_{12},{\psi }_{22},{\psi }_{31},{\psi }_{32})\in L,\hfill \\ & & f=(f_{11},f_{12},f_{22},f_{31},f_{32})\in L_2((0,T),H^{-1})\hfill \end{array}$$

and the conditions (21), (22), and (35) are fulfilled. Then, the reduced (truncated) IBVP (8), (9), (11)–(13), (15)–(17), and (23)–(29) has a unique solution $u\in W(0,T)$; it continuously depends on f and $u_0$.

The a priori estimate holds:

$${\parallel u\parallel }_W^2\le C\left(T\right)\left(\parallel u_0{\parallel }_L^2+\underset{0}{\overset{T}{\int }}{\parallel f(·,t)\parallel }_{H^{-1}}^{2}dt\right),$$

where $C\left(T\right)>0$.

5. Positivity Preservation

The model described in Section 2 can be used to describe the transport of nonnegative concentration species of contaminants or the temperature of interacting bodies, etc. So, the concentration and temperature are nonnegative over time if the initial state and boundary conditions of the medium are nonnegative. This is an essential requirement for mathematical modeling and is called nonnegativity (or positivity) preservation.

Following reference [26], we introduce the following notations. For function $w_i\in H^1\left({\Omega }_i\right)$, $i=1,2,3$, we denote the negative and positive parts by $w_i^{-}$ and $w_i^{+}$, respectively. That is $w_i=w_i^{+}+w_i^{-}$, where $w_i^{+}=max(w_i,0)\ge 0$ and $w_i^{-}=min(w_i,0)\le 0$. By ([26], Lemma 7.6, p. 150), it follows that

$$\frac{d{w}_i^{+}}{dx}=\left\{\begin{array}{cc}{\displaystyle \frac{d{w}_i}{dx},}\hfill & \mathrm{if}{w}_i\ge 0,\hfill \\ 0,\hfill & \mathrm{otherwise}\hfill \end{array}\right.,\frac{d{w}_i^{-}}{dx}=\left\{\begin{array}{cc}{\displaystyle \frac{d{w}_i}{dx},}\hfill & \mathrm{if}{w}_i<0,\hfill \\ 0,\hfill & \mathrm{otherwise}\hfill \end{array}\right..$$

These inequalities imply

$$w_i^{+}{w}_i^{-}=\frac{d{w}_i^{+}}{dx}\frac{d{w}_i^{-}}{dx}=w_i^{-}\frac{d{w}_i^{+}}{dx}+w_i^{+}\frac{d{w}_i^{-}}{dx}=0,\mathrm{a}.\mathrm{e}.\mathrm{in}{\Omega }_i,i=1,2,3.$$

For some $T>0$, let us consider the following inequalities:

$$\begin{array}{ccc}& & {\mathcal{L}}_{1i}{u}_i\ge 0,\mathrm{in}{Q}_T,i=1,2,3,\hfill \\ & & {\mathcal{L}}_2{u}_{22}\ge 0,\mathrm{in}{Q}_{2T},\hfill \\ & & {\mathcal{L}}_{3i}{u}_3\ge 0,\mathrm{in}{Q}_{3T},i=1,2,\hfill \end{array}$$

(42)

$$\begin{array}{ccc}& & l_2(u_{11},u_{31})\ge 0,i=1,2,0<t<T,\hfill \end{array}$$

(43)

$$\begin{array}{ccc}& & l_{0i}\left(u_{1i}\right)\ge 0,l_{3i}\left(u_{3i}\right)\ge 0,i=1,2,0<t<T.\hfill \end{array}$$

Let V (see [27], [p. 6]) be the Banach space of all functions in $H^{1,0}$ with finite norm

$${\parallel w\parallel }_H=\underset{1\le i\le 3}{max}\left(\underset{0\le t\le T}{sup}{\parallel w(·,t)\parallel }_{L_2\left({\Omega }_i\right)}+{∥\frac{\partial w_i}{\partial x}∥}_{L_2\left(Q_{iT}\right)}\right),$$

where $w_i=u_{1i}$, $i=1,2$, $w_2=u_{22}$, $w_{3+i}=u_{3i}$, $i=1,2$.

A function $u=(u_{11},u_{12},u_{22},u_{31},u_{32})\in V$ is said to satisfy (42) weakly, if for any $w=(w_1,w_2,w_3,w_4,w_5)\in H^{1,1}$ if for $w_i\ge 0$, $i=1,2,\cdots ,5$:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & {\beta }_3\underset{0}{\overset{l_1}{\int }}{u}_{11}(x,t)w_1(x,t)dx+{\beta }_1\underset{l_2}{\overset{l_3}{\int }}{u}_{31}(x,t)w_4(x,t)dx-{\beta }_3\underset{0}{\overset{l_1}{\int }}{u}_{11}(x,0)w_1(x,0)dx\hfill \\ & -{\beta }_1\underset{l_2}{\overset{l_3}{\int }}{u}_{31}(x,0)w_4(x,0)dx-{\beta }_3\underset{0}{\overset{t}{\int }}\underset{0}{\overset{l_1}{\int }}{u}_{11}(x,\tau )\frac{\partial w_1}{\partial t}(x,\tau )dxd\tau \hfill \\ & -{\beta }_1\underset{0}{\overset{t}{\int }}\underset{l_2}{\overset{l_3}{\int }}{u}_{31}(x,\tau )\frac{\partial w_4}{\partial t}(x,\tau )dxd\tau +\underset{0}{\overset{l_1}{\int }}{u}_{12}(x,t)w_2(x,t)dx+\underset{l_1}{\overset{l_2}{\int }}{u}_{22}(x,t)w_2(x,t)dx\hfill \\ & \underset{l_2}{\overset{l_3}{\int }}{u}_{32}(x,t)w_5(x,t)dx-\underset{0}{\overset{l_1}{\int }}{u}_{12}(x,0)w_2(x,0)dx-\underset{l_1}{\overset{l_2}{\int }}{u}_{22}(x,0)w_3(x,0)dx\hfill \\ & -\underset{l_2}{\overset{l_3}{\int }}{u}_{32}(x,0)w_5(x,0)dx+\underset{0}{\overset{t}{\int }}\mathcal{A}\left(u(·,\tau ),w(·,\tau )\right)d\tau \ge 0,\hfill \end{array}\end{array}$$

(44)

for almost every $t\in (0,T)$. It is to show that if $u_{1i}\in C^{2,1}\left({\overline{Q}}_{1T}\right)$, $i=1,2$, $u_{22}\in C^{2,1}\left({\overline{Q}}_{2T}\right)$, $u_{3i}\in C^{2,1}\left({\overline{Q}}_{3T}\right)$, $i=1,2$, satisfy (42) and (43) pointwise, then (44) holds. Conversely, if $u\in V$ satisfy (44) and $u_{1i}$, $i=1,2$, $u_{22}$, $u_{3i}$, $i=1,2$ are sufficiently smooth, then they also satisfy (42) and (43) pointwise. We are now in a position to prove the following assertion.

Theorem 3.

Let $u=(u_{11},u_{12},u_{22},u_{31},u_{32})\in V$ satisfy (42) and (43) weakly, and let (21), (22), and (34) hold. If $u_{1i}(x,0)\ge 0$, $i=1,2$, $u_{22}(x,0)\ge 0$, $u_{3i}(x,0)\ge 0$, $i=1,2$, then $u_{1i}(x,t)\ge 0$, $i=1,2$ a.e. in $Q_{1T}$, $u_{22}(x,0)\ge 0$ a.e. in $Q_{2T}$, $u_{3i}(x,0)\ge 0$, $i=1,2$, a.e., in $Q_{3T}$.

Proof. 

Following the methodology of ([18], [Chapter 3]), using the Steklov averaging and passing to the limit, we can take $w=-u^{-}$ in (44), to obtain

$$\begin{array}{ccc}& & \frac{1}{2}{\beta }_3\underset{0}{\overset{l_1}{\int }}{\left(u_{11}^{-}(x,t)\right)}^{2}dx+\frac{1}{2}{\beta }_1\underset{l_2}{\overset{l_3}{\int }}{\left(u_{31}^{-}(x,t)\right)}^{2}dx-\frac{1}{2}{\beta }_3\underset{0}{\overset{l_1}{\int }}{\left(u_{11}^{-}(x,0)\right)}^{2}dx\hfill \\ & & -\frac{1}{2}{\beta }_1\underset{l_2}{\overset{l_3}{\int }}{\left(u_{31}^{-}(x,0)\right)}^{2}dx+\frac{1}{2}\underset{0}{\overset{l_1}{\int }}{\left(u_{12}^{-}(x,t)\right)}^{2}dx+\frac{1}{2}\underset{l_1}{\overset{l_2}{\int }}{\left(u_{22}^{-}(x,t)\right)}^{2}dx\hfill \\ & & +\frac{1}{2}\underset{l_2}{\overset{l_3}{\int }}{\left(u_{32}^{-}(x,t)\right)}^{2}dx-\frac{1}{2}\underset{0}{\overset{l_1}{\int }}{\left(u_{12}^{-}(x,0)\right)}^{2}dx-\underset{l_1}{\overset{l_2}{\int }}{\left(u_{22}^{-}(x,0)\right)}^{2}dx-\underset{l_2}{\overset{l_3}{\int }}{\left(u_{32}^{-}(x,0)\right)}^{2}dx\hfill \\ & & +{\beta }_3\underset{0}{\overset{l_1}{\int }}{\left(p_{11}\left(x\right)\frac{\partial u_{11}}{\partial x}\frac{\partial u_{11}^{-}}{\partial x}+{\mu }_{11}^1{u}_{11}{u}_{11}^{-}+{\mu }_{12}^1{u}_{12}{u}_{11}^{-}\right)}^{2}dx\hfill \\ & & +{\beta }_1\underset{l_2}{\overset{l_3}{\int }}{\left(p_{31}\left(x\right)\frac{\partial u_{31}}{\partial x}\frac{\partial u_{31}^{-}}{\partial x}+{\mu }_{11}^3{u}_{31}{u}_{31}^{-}+{\mu }_{12}^3{u}_{32}{u}_{31}^{-}\right)}^{2}dx\hfill \\ & & +\underset{0}{\overset{l_1}{\int }}{\left(p_{12}\left(x\right)\frac{\partial u_{12}}{\partial x}\frac{\partial u_{12}^{-}}{\partial x}+{\mu }_{21}^1{u}_{11}{u}_{12}^{-}+{\mu }_{22}^1{u}_{12}{u}_{12}^{-}\right)}^{2}dx\hfill \\ & & +\underset{l_1}{\overset{l_2}{\int }}{\left(p_{23}\left(x\right)\frac{\partial u_{22}}{\partial x}\frac{\partial u_{22}^{-}}{\partial x}+{\mu }_2^2{u}_{22}{u}_{22}^{-}\right)}^{2}dx\hfill \\ & & +\underset{l_2}{\overset{l_3}{\int }}{\left(p_{32}\left(x\right)\frac{\partial u_{32}}{\partial x}\frac{\partial u_{32}^{-}}{\partial x}+{\mu }_{21}^3{u}_{32}{u}_{31}^{-}+{\mu }_{22}^3{u}_{32}{u}_{32}^{-}\right)}^{2}dx\hfill \\ & & \le \underset{0}{\overset{t}{\int }}(-{\beta }_3{\alpha }_1{u}_{11}(l_1,\tau )u_{11}^{-}(l_1,\tau )-{\beta }_1{\alpha }_3{u}_{31}(l_2,\tau )u_{31}^{-}(l_2,\tau )\hfill \\ & & +{\beta }_1{\beta }_2[u_{11}(l_1,\tau )u_{31}^{-}(l_2,\tau )+u_{31}(l_2,\tau )u_{11}^{-}(l_1,\tau )])d\tau \hfill \end{array}$$

(45)

Next, we cultivate (45) using the relation

$$\frac{\partial w_i}{\partial x}\frac{\partial w_i^{-}}{\partial x}={\left(\frac{\partial u_i^{-}}{\partial x}\right)}^2,w_i{w}_i^{-}={\left(w_i^{-}\right)}^2,$$

$$\begin{array}{ccc}& & u_{11}(l_1,t)u_{31}^{-}(l_2,t)+w_{31}(l_2,t)u_{11}^{-}(l_1,t)\hfill \\ & & =2u_{11}^1(l_1,t)u_{31}^{-}(l_2,t)+u_{11}^{+}(l_1,t)u_{31}^{-}(l_2,t)+u_{11}^{-}(l_1,t)u_{31}^{-}(l_2,t)\hfill \\ & & \le 2u_{11}^{-}(l_1,t)u_{31}(l_2,t).\hfill \end{array}$$

Now, using assumptions $u_{1i}^{-}(x,0)=0$, $i=1,2$, $u_{22}^{-}(x,0)=0$, $u_{3i}^{-}(x,0)=0$, $i=1,2$, we derive

$$\begin{array}{ccc}& & \frac{1}{2}{\beta }_3\underset{0}{\overset{l_1}{\int }}{\left(u_{11}^{-}(x,t)\right)}^{2}dx+{\beta }_3\underset{0}{\overset{t}{\int }}\underset{0}{\overset{l_1}{\int }}{\left(p_{11}(x,t){\left(\frac{\partial u_{11}^{-}}{\partial x}\right)}^2+{\mu }_1{\left(u_{11}\right)}^2\right)}^{2}dxd\tau \hfill \\ & & +\frac{1}{2}{\beta }_1\underset{l_2}{\overset{l_3}{\int }}{\left(u_{31}^{-}(x,t)\right)}^{2}dx+{\beta }_1\underset{0}{\overset{t}{\int }}\underset{l_2}{\overset{l_3}{\int }}{\left(p_{31}(x,t){\left(\frac{\partial u_{31}^{-}}{\partial x}\right)}^2+{\mu }_3{\left(u_{31}\right)}^2\right)}^{2}dxd\tau \hfill \\ & & \le -\underset{0}{\overset{t}{\int }}[{\beta }_3{\alpha }_1{\left(u_{11}^{-}(l_1,t)\right)}^2+{\beta }_1{\alpha }_3{\left(u_{31}^{-}(l_2,t)\right)}^2-2{\beta }_1{\beta }_3{u}_{11}(l_1,\tau )u_{31}^{-}(l_2,\tau )]d\tau \le 0.\hfill \end{array}$$

(46)

Here, ${\mu }_1$ and ${\mu }_2$ are nonnegative constants, obtained by the trick used at the estimation of the integrals $I_1$, $I_2$, and $I_3$ of Lemma 1.

We also use inequalities (21), (22), and (34). Therefore, (46) implies $u_{11}^{-}(x,t)=u_{11}^{+}(x,t)$, a.e., in $Q_{1T}$ and $u_{31}^{-}(x,t)=u_{31}^{+}(x,t)\ge 0$, a.e., in $Q_{3T}$. □

6. Numerical Test Example

In this section, we discuss the numerical implementation of the proposed approach. The illustrative example is the case of a three-layer advection–dispersion reaction model with two species in the first layer and one species in both the second and third layers in the more general case of time-dependent diffusion coefficients. Specifically, we consider the system (8)–(10), (12), (14), and (16) with the following model parameters

$$\begin{array}{cc}& {\displaystyle {\mu }_{11}^1=-2,{\mu }_{12}^1=0,{\mu }_{21}^1=3,{\mu }_{22}^1=-2,{\mu }_1^2=1,{\mu }_{11}^3=-1,{\mu }_{12}^3=0,R_{21}=1,}\\ & {\displaystyle p_{11}=t(1-4x)+1,p_{12}=2t+1,p_{21}=1,p_{31}=(t^2+2x+1)(x-1)+1,}\\ & f_{21}(x,t)=e^t(e^x+e^{-x}),l_0=-1/4,l_1=1/4,l_2=1,l_3=2.\end{array}$$

(47)

It is easy to check that function $u_{21}=e^t(e^x+e^{-x})$ satisfies Equation (10). Further we determine $v_1$ and $v_2$ from (25)–(27) and ${\varphi }_1\left(t\right)$, and ${\varphi }_2\left(t\right)$ from (19) for $w(x,t)=e^{t+x}$, as follows:

$$v_1=\frac{e^{x-l_1}+e^{l_1-x}}{2},v_2=\frac{e^{l_2-x}+e^{x-l_2}}{2},{\phi }_1\left(t\right)=\frac{2e^{t-l_2}}{e^l-e^{-l}},{\phi }_2\left(t\right)=\frac{2e^{t-l_1}}{e^l-e^{-l}},$$

(48)

where $l=l_2-l_1$.

Note that, for the present solution $u_{21}$, we may obtain the representation (18) in a more straightforward manner, namely

$$v_1\left(x\right)=e^x,v_2\left(x\right)=e^{-x},{\phi }_1\left(t\right)={\phi }_2\left(t\right)=e^t,w(x,t)=0.$$

(49)

Here, for the computational tests, we apply a general approach by using (19) and (25)–(27), i.e., the representation (18) via (48). Of course, expressions (48) and (49) of (18) produce similar results.

Next, we determine ${\alpha }_i$, ${\beta }_i$, ${\gamma }_i$, $i=1,3$ from (20). Thus, the reduced truncated problem with simplified notations, i.e., $U_1:=u_{11}$, $U_2:=u_{12}$, $U_3=u_{31}$, $p_1:=p_{11}$, $p_2:=p_{12}$, $p_3:=U_{31}$, $a_1:={\mu }_{11}^1$, $a_2:={\mu }_{22}^1$, $b_2:={\mu }_{21}^1$, $a_3:={\mu }_{11}^3$, has the following representation

$$\begin{array}{ccc}& & \frac{\partial U_1}{\partial t}-\frac{\partial }{\partial x}\left(p_1(x,t)\frac{\partial U_1}{\partial x}\right)-a_1{U}_1=f_1(x,t),(x,t)\in Q_{1T},\hfill \end{array}$$

(50)

$$\begin{array}{ccc}& & \frac{\partial U_2}{\partial t}-\frac{\partial }{\partial x}\left(p_2(x,t)\frac{\partial U_2}{\partial x}\right)-a_2{U}_2-b_2{U}_1=f_2(x,t),(x,t)\in Q_{1T},\hfill \end{array}$$

(51)

$$\begin{array}{ccc}& & \frac{\partial U_3}{\partial t}-\frac{\partial }{\partial x}\left(p_3(x,t)\frac{\partial U_3}{\partial x}\right)-a_3{U}_3=f_3(x,t),(x,t)\in Q_{3T},\hfill \end{array}$$

(52)

$$\begin{array}{ccc}& & p_1(l_1,t)\frac{\partial U_1}{\partial x}(l_1,t)+{\alpha }_1{U}_1(l_1,t)-{\beta }_1{U}_3(l_2,t)={\gamma }_1\left(t\right),t\in [0,T],\hfill \end{array}$$

(53)

$$\begin{array}{ccc}& & -p_3(l_2,t)\frac{\partial U_3}{\partial x}(l_2,t)+{\alpha }_3{U}_3(l_2,t)-{\beta }_3{U}_1(l_1,t)={\gamma }_3\left(t\right),t\in [0,T],\hfill \end{array}$$

(54)

$$\begin{array}{ccc}& & U_1(l_0,t)={\lambda }_{10}\left(t\right),U_2(l_0,t)={\lambda }_{20}\left(t\right),U_2(l_1,t)={\lambda }_{21}\left(t\right),U_3(l_3,t)={\lambda }_{33}\left(t\right),t\in [0,T],\hfill \end{array}$$

(55)

$$\begin{array}{ccc}& & U_i(x,0)={\psi }_i\left(x\right),i=1,2,x\in {\Omega }_1,U_3(x,0)={\psi }_3\left(x\right),x\in {\Omega }_3.\hfill \end{array}$$

(56)

In the vein of the present work, we construct explicit–implicit discretization for solving Problems (50)–(56), which perform fast and save computational time.

We define a uniform temporal mesh (${\overline{\omega }}_{\tau }$) and spatial (${\overline{\omega }}_{h_i}$, $i=1,3$) meshes in ${\Omega }_1$ and ${\Omega }_3$, in order to solve Problems (50)–(56) numerically

$$\begin{array}{ccc}& & {\overline{\omega }}_{h_1}=\{x_{1,j_1}=l_0+j_1{h}_1,j_1=0,1,\cdots ,N_1,h_1=(l_1-l_0)/N_1\},\hfill \\ & & {\displaystyle {\overline{\omega }}_{h_3}=\{x_{3,j_3}=l_2+j_3{h}_3,j_3=0,1,\cdots ,N_3,h_3=(l_3-l_2)/N_3\},}\hfill \\ & & {\overline{\omega }}_{\tau }=\{t_n=n\tau ,n1=0,1,\cdots ,M,\tau =T/M\}.\hfill \end{array}$$

The numerical solution at the grid node $(x_{i,j_i},t_n)\in {\overline{\omega }}_{h_i}\times {\overline{\omega }}_{\tau }$ is denoted by $U_{i,j_i}^n$, i.e., $U_{i,j_i}^n=U_i(x_{i,j_i},t_n)$. Further, we introduce the notations

$$U_{i,x_{j_i}}^n=\frac{U_{i,j_i+1}^n-U_{i,j_i}^n}{h_i},U_{i,{\overline{x}}_{j_i}}^n=\frac{v_{i,j_i}^n-U_{i,j_i-1}^n}{h_i},p_{j,i_j\pm 1/2}=p_i\left(x_{i,j_i\pm 1/2}\right).$$

For the spatial approximation of (50)–(56), we apply the finite volume method. We consider dual meshes of cell-centered grid nodes $x_{i,j_i-1/2}=x_{i,j_i}-h_i/2$, $j_i=0,1,\cdots ,N_i+1$, where $x_{i,-1/2}=x_{i,0}$, $x_{i,N_i+1/2}=x_{i,N_i}$, $i=1,3$.

Integrating Equations (50) and (51) over the volumes $(x_{1,j_1-1/2}$, $x_{1,j_1+1/2})$, $j_1=1,2,\cdots ,$$N_1-1$, we obtain the discretization at inner grid nodes of the discrete domain of ${\overline{\omega }}_{h_1}\times [0,T]$. Similarly, integrating Equation (53) over the volumes $(x_{3,j_3-1/2}$, $x_{3,j_3+1/2})$, $j_3=1,2,\cdots ,N_3-1$, we obtain the discretization at the inner grid nodes of ${\overline{\omega }}_{h_3}\times [0,T]$.

Then we integrate Equation (50) over the volume $(x_{1,N_1-1/2},x_{1,N_1})$, and using the boundary condition (54), we obtain the approximation at boundary $x=l_1$. In the same manner, we derive the approximation at interface $x=l_2$, dealing with the boundary condition (53).

Further, applying implicit–explicit time stepping, we obtain

$$\begin{array}{c}\hfill \begin{array}{cc}& \frac{U_{1,j_1}^{n+1}}{\tau }-\frac{1}{h_1}\left(p_{1,j_1+1/2}^{n+1}{U}_{1,x_{j_1}}^{n+1}-p_{1,j_1-1/2}^{n+1}{U}_{1,{\overline{x}}_{j_1}}^{n+1}\right)-a_1{U}_{1,j_1}^{n+1}\hfill \\ & =f_{1,j_1}^{n+1}+\frac{U_{1,j_1}^n}{\tau },j_1=1,2,\cdots ,N_1-1,\hfill \\ & \frac{U_{2,j_1}^{n+1}}{\tau }-\frac{1}{h_1}\left(p_{2,j_1+1/2}^{n+1}{U}_{2,x_{j_1}}^{n+1}-p_{2,j_1-1/2}^{n+1}{U}_{2,{\overline{x}}_{j_1}}^{n+1}\right)-a_2{U}_{2,j_1}^{n+1}\hfill \\ & =b_2{U}_{1,j_1}^n+f_{2,j_1}^{n+1}+\frac{U_{1,j_1}^n}{\tau },j_1=1,2,\cdots ,N_1-1,\hfill \\ & \frac{U_{3,j_3}^{n+1}}{\tau }-\frac{1}{h_3}\left(p_{3,i_3+1/2}^{n+1}{U}_{3,x_{i_3}}^{n+1}-p_{3,j_3-1/2}^{n+1}{U}_{3,{\overline{x}}_{i_3}}^{n+1}\right)-a_3{U}_{3,j_3}^{n+1}\hfill \\ & =f_{3,j_3}^{n+1}+\frac{U_{3,j_3}^n}{\tau },j_3=1,2,\cdots ,N_3-1,\hfill \\ & \left(\frac{1}{\tau }+\frac{2{\alpha }_1}{h_1}\right)U_{1,N_1}^{n+1}+\frac{2}{h_1}{p}_{1,N_1-1/2}^{n+1}{U}_{1,{\overline{x}}_{N_1}}^{n+1}-a_1{U}_{1,N_1}^{n+1}\hfill \\ & =\frac{2{\beta }_1}{h_1}{U}_{3,0}^n+\frac{2}{h_1}{\gamma }_1^{n+1}+f_{1,N_1}^{n+1}+\frac{U_{1,N_1}^n}{\tau },\hfill \\ & \left(\frac{1}{\tau }+\frac{2{\alpha }_3}{h_2}\right)U_{3,0}^{n+1}-\frac{2}{h_2}{p}_{2,1/2}^{n+1}{U}_{3,x_0}^{n+1}-a_3{U}_{3,0}^{n+1}\hfill \\ & =b_2{U}_{1,N_1}^n+\frac{2{\beta }_3}{h_2}{U}_{1,N_1}^n+\frac{2}{h_2}{\gamma }_3^{n+1}+f_{3,0}^{n+1}+\frac{U_{3,0}^n}{\tau },\hfill \\ & U_{1,0}^{n+1}={\lambda }_{10}^{n+1},U_{2,0}^{n+1}={\lambda }_{20}^{n+1},U_{2,N_2}^{n+1}={\lambda }_{21}^{n+1},U_{3,N_2}^{n+1}={\lambda }_{33}^{n+1},\hfill \\ & U_{i,j_1}^0={\varphi }_{i,j_1},i=1,3,j_1=0,1,\cdots ,N_1,u_{3,j_3}^0={\varphi }_{3,j_3}{j}_3=0,1,\cdots ,N_3.\hfill \end{array}\end{array}$$

(57)

The performance of the numerical scheme (57) requires the inversion of two matrices of size $(N_1,N_1)$ and one of size $(N_3,N_3)$ at each time level, in contrast to the corresponding fully implicit scheme, which involves the matrix inversion of a large $(2N_1+N_3,2N_1+N_3)$ matrix at each time layer. Meanwhile, as it is well-known, implicit–explicit discretizations have better stability properties than explicit schemes [28].

We illustrate the accuracy and efficiency of the proposed discrete scheme (57).

For the numerical tests, we chose the residual functions $f_{11}(x,t)$, $f_{12}(x,t)$, $f_{31}(x,t)$, such that

$$\begin{array}{ccc}& & u_{11}(x,t)=e^t(e^x+e^{-x})+(2+t){(x-l_1)}^2,\hfill \\ & & u_{31}(x,t)=e^t(e^x+e^{-x})+t^2\left(sin\left(\pi x\right)-\pi (l_2-x)\right),\hfill \\ & & u_{12}(x,t)=e^{2t}cos(\pi x/4)+e^{5x},\hfill \end{array}$$

to be exact solutions of (8)–(10), (12), (14), (16), and (47), namely systems (50)–(56) and (47).

The errors in maximal norm and convergence rates are computed as follows

$$\begin{array}{cc}& {\displaystyle \mathcal{E}\left(U_i\right)=U_i(t_n,x_{i,j_i})-U_{i,j_i}^n,{\mathcal{E}}^{N_i,M}\left(U_i\right)=\underset{0\le n\le M}{max}\underset{0\le j_i\le N_i}{max}\left|\right\}\left\}E\left(U_i\right)\right|,}\\ & {\displaystyle C{R}_h\left(U_i\right)={log}_2\frac{{\mathcal{E}}^{2N_i,M}\left(U_i\right)}{{\mathcal{E}}^{N_i,M}\left(U_i\right)},C{R}_{\tau }={log}_2\frac{{\mathcal{E}}^{N_i,2M}\left(U_i\right)}{{\mathcal{E}}^{N_i,M}\left(U_i\right)}.}\end{array}$$

We illustrate the spatial (

Table 1. Errors and spatial order of convergence of numerical solutions $U_i$, $i=1,2,3$.

${\mathit{N}}_1$	${\mathcal{E}}^{{\mathit{N}}_1,\mathit{M}}\left({\mathit{U}}_1\right)$	${\mathit{CR}}_{\mathit{h}}\left({\mathit{U}}_1\right)$	${\mathcal{E}}^{{\mathit{N}}_1,\mathit{M}}\left({\mathit{U}}_2\right)$	${\mathit{CR}}_{\mathit{h}}\left({\mathit{U}}_2\right)$	${\mathcal{E}}^{{\mathit{N}}_2,\mathit{M}}\left({\mathit{U}}_3\right)$	${\mathit{CR}}_{\mathit{h}}\left({\mathit{U}}_3\right)$
10	$3.709\times {10}^{-3}$		$1.377\times {10}^{-2}$		$4.737\times {10}^{-3}$
20	$9.325\times {10}^{-4}$	1.992	$3.446\times {10}^{-3}$	1.998	$1.187\times {10}^{-3}$	1.996
40	$2.335\times {10}^{-4}$	1.998	$8.619\times {10}^{-4}$	2.000	$2.972\times {10}^{-4}$	1.999
80	$5.838\times {10}^{-5}$	1.999	$2.155\times {10}^{-4}$	2.000	$7.427\times {10}^{-5}$	2.000
160	$1.460\times {10}^{-5}$	2.000	$5.387\times {10}^{-5}$	2.000	$1.857\times {10}^{-5}$	2.000

) and temporal (

Table 2. Errors and temporal order of convergence of numerical solutions $U_i$, $i=1,2,3$.

M	${\mathcal{E}}^{{\mathit{N}}_1,\mathit{M}}\left({\mathit{U}}_1\right)$	${\mathit{CR}}_{\mathit{h}}\left({\mathit{U}}_1\right)$	${\mathcal{E}}^{{\mathit{N}}_1,\mathit{M}}\left({\mathit{U}}_2\right)$	${\mathit{CR}}_{\mathit{h}}\left({\mathit{U}}_2\right)$	${\mathcal{E}}^{{\mathit{N}}_2,\mathit{M}}\left({\mathit{U}}_3\right)$	${\mathit{CR}}_{\mathit{h}}\left({\mathit{U}}_3\right)$
40	$2.578\times {10}^{-2}$		$1.906\times {10}^{-2}$		$1.413\times {10}^{-3}$
80	$1.300\times {10}^{-2}$	0.9871	$9.801\times {10}^{-3}$	0.960	$7.870\times {10}^{-4}$	0.8441
160	$6.531\times {10}^{-3}$	0.9936	$4.968\times {10}^{-3}$	0.980	$4.139\times {10}^{-4}$	0.9270
320	$3.273\times {10}^{-3}$	0.9968	$2.501\times {10}^{-3}$	0.990	$2.121\times {10}^{-4}$	0.9646
640	$1.638\times {10}^{-3}$	0.9984	$1.255\times {10}^{-3}$	0.995	$1.074\times {10}^{-4}$	0.9826
1280	$8.195\times {10}^{-4}$	0.9992	$6.283\times {10}^{-4}$	0.998	$5.399\times {10}^{-5}$	0.9913

) orders of convergence in the maximal norm. To this aim, in Table 1, computations are performed for $N_2=2N_1$ (i.e., $h_1=h_2=h$) and $\tau =h^2$, while in Table 2, for all runs, we take $N_1=$, $N_2=2N_1$, namely $h_1=h_2=\tau$. We observe that the accuracy of the discretization (57) is $O(\tau +h^2)$.

In Figure 3 and Figure 4, we plot the numerical solution and error throughout the entire computational domain for $\tau =h^2$, $N_1=21$, $N_2=41$.

7. Conclusions

In this paper, we began with a multispecies multilayer system of porous media; exploiting the fact that an exact solution could be constructed in some layers, we formulated a new problem with fewer layers and nonstandard interface relations. For the new problem, we set up a weak formulation and then proved the existence and uniqueness of the solution in appropriate Sobolev-like spaces. Since the unknown functions model concentrations or temperatures of physical processes, we also discussed their nonnegativity.

We illustrated the numerical realization of the proposed approach by using fast and accurate implicit–explicit-based finite volume discretization for the reduced truncated problem.

It is very natural to ask whether the results obtained for the linear parabolic systems can be extended to nonlinear systems. For biology models with a standard family of quadratic nonlinearities, see reference [19], and for nonlinear interface conditions, see [2].

The investigations in this paper are mainly analytical, providing a background for a numerical study. Since the applications for the numerical solutions of the multilayer multispecies problem are very important, we plan to further develop in this direction. In references [1,2,3,16], computational results for multispecies models of the type described in Section 2 were obtained. We plan to numerically solve such models on the basis of the reduced truncated problem, derived in the present paper. It would also be interesting to solve inverse coefficient or source problems for models of the type introduced in Section 2, applying the methods proposed in [15,20,29].

One future direction of our work is to consider this model problem in extreme cases, where some retardation effects are very large, e.g., at least $R_{ij}>{10}^4$. We plan to apply the methodology presented in [30,31] in order to handle the corresponding singularly perturbed equations.