Approximate analytical solutions for multispecies convection-dispersion transport equation with variable parameters

1 Introduction

Pollution is caused by a number of sources, such as agricultural fertilizers, erosion, industries, energy generation, waste disposal, hospital wastage, vehicular transport, nuclear waste, and domestic waste. A significant portion of pollution generated by these sources finds its way either directly or indirectly into ground water (Sposito et al., 1979; Domenico 1987; Bear 1988; Clement et al., 1998; Batu 2005). Often pollutants are different species, and their chemical characteristics are different. In general, pollutant transport in groundwater is modeled using the conservation of mass of water, flux law for flow of water, conservation of mass of pollutant, and flux law for pollutant transport. Depending on the type, a pollutant can be physical, chemical, or biological. A chemical pollutant can be conservative or non-conservative for which production and decay functions need to be specified. When these equations are combined, the resulting equation is the convection-dispersion equation. This equation is for a specific species. This paper discusses the case of two species pollutant transport.

There are several works done for the analysis of multispecies pollutant transport problems. These works adopted either analytical or numerical methods. The analytical solutions were obtained using Fourier transforms, Laplace transform, general integral transforms, decomposition methods, series solutions, etc. (Lunn et al., 1996; Van Genuchten 1985; Fujikawa and Fukui 1990; Sun and Clement 1999; Sun et al., 1999; Chamkha 2005; Slodička and Balážová, 2008; Slodička and Balážová, 2010; Natarajan and Kumar 2010; Chen et al., 2012; Simpson and Ellery 2014). Numerical methods basically adopted the finite difference method (Arnold et al., 2017; Natarajan and Kumar 2018). However, the afore-mentioned works used constant transport parameter—a consideration that makes it easier to deal it mathematically. In reality, it is necessary to model the system using variable parameters in order to fully understand the pollutant dispersion behaviour. There are a few works available for the variable transport parameters; however, some are for single species or restricted formulation. Chaudhary and Singh (2020) used the homotopy analysis method (HAM) for two-species convection-dispersion transport equation, considering spatial-temporal varying dispersion coefficient and seepage velocity. They used the standard HAM to deal with the system of PDEs governing the phenomenon. However, since the inception of HAM, there have been several modifications done. The most popular variants of the method are homotopy perturbation method (HPM) and optimal homotopy asymptotic method (OHAM). On the other hand, it is worth mentioning that the transport of heavy metals in groundwater and the interactions between them and the soil medium are key concepts in environmental engineering. Such transport processes are useful for studying the damage to the environment caused by the smelting of metallic materials, the seepage of landfill leachates, the use of pesticides and chemical fertilizers, the treatment of municipal wastewater, etc. Furthermore, the solid particles (bacteria, silicon powders, etc.) present in the soil interact with the transport mechanism of the heavy metals by accelerating the migration rate of those. Thus, in order to choose a suspended matter of particular for removing the heavy metals from soils, it is crucial to understand the coupling mechanism. In fact, because of the suspended particles adsorb the heavy metals, they show an influence on the migration process of the metals in soil under seepage conditions. There are several factors affecting the coupling mechanism of heavy metals and suspended particles, such as the size and concentration of the particles, seepage velocity, types of heavy metals, porous media, coupling mechanism (solid-solid, solid-liquid, gas-liquid, etc.). The detail discussion can be found in Bai et al. (2021a), Bai et al. (2021b).

The core idea of homotopy-based methods is based on the concept of homotopy from topology. It creates a continuous mapping that deforms continuously to obtain one function from another. Liao (1992) used this concept to develop an analytical method for the solution of non-linear problems in terms of a series solution. Since then, it has been used extensively in different areas of science and engineering (Liao, 2003; Liao, 2012). On the other hand, He (1999) presented an analytical methodology named HPM. Recently, Marinca and Herişanu (2008) extended the concept of HAM using an approximation method to derive the so-called OHAM. These methods have their own advantages/disadvantages in terms of their applicability. Therefore, the present work presents three kinds of approximate series solutions using HAM, HPM, and OHAM. This paper considers the case of spatially varying transport parameters. Indeed, the methodologies can be adopted for other cases such as the ones provided in Chaudhary and Singh (2020).

2 Brief overview of homotopy-based methods

Here, first we describe the homotopy-based methods in a general framework considering a system of PDEs. Before doing so, it is pertinent to mention that all these methods are based on the mathematical concept called ‘homotopy’ from topology (Liao 1992). Two objects (mathematical) are homotopic if one can be continuously deformed into the other. Mathematically, a homotopy $\mathcal{H}\left(t;q\right)$ between two functions $f\left(t\right)$; $g\left(t\right)$, where $t$ is a dimension (space or time), is itself a continuous function, defined as $\mathcal{H}:T\times \left[\mathrm{0,1}\right]\to U$ and satisfies $\mathcal{H}\left(t;q\right)=f\left(t\right)$ when $q=0$ and $\mathcal{H}\left(t;q\right)=g\left(t\right)$ when $q=1$, where $T$ and $U$ are the topological spaces. This shows that as $q$ goes from 0 to 1, $\mathcal{H}\left(t;q\right)$ varies from $f\left(t\right)$ to $g\left(t\right)$. The functions $f\left(t\right)$ and $g\left(t\right)$ are called homotopic. For example, in 2D, a circle can be continuously deformed into an ellipse or a square; similarly, in 3D, a doughnut and a coffee cup are homotopic. Since algebraic or differential equations represent curves (or functions), the concept of homotopy can be employed for solving non-linear differential or algebraic equations. Using this concept, Liao (1992) proposed the so-called homotopy analysis method (HAM). Two other popular variants of this method are homotopy perturbation method (HPM) and optimal homotopy asymptotic method (OHAM). These three methods are employed in this paper.

2.1 Homotopy analysis method

Let a system of PDEs be written as:

${\mathcal{N}}_i\left(y_i\left(x,t\right)\right)=0(1)$

where ${\mathcal{N}}_i$ are the non-linear operators or the original operators of the system of equations; $y_i$ are the unknown variables for say $i=\mathrm{1,2},\dots ,n$; $x$ and $t$ are the independent variables. Now, the zeroth-order deformation equation can be constructed as follows (Liao 2003):

$\left(1-q\right){\mathcal{L}}_i\left[{\Phi }_i\left(x,t;q\right)-y_{i,0}\left(x,t\right)\right]=q{\mathrm{ħ}}_i{H}_i\left(x,t\right){\mathcal{N}}_i\left[{\Phi }_i\left(x,t;q\right)\right],i=\mathrm{1,2},\dots ,n(2)$

subject to the initial and boundary conditions:

$\mathcal{B}\left({\Phi }_i,\frac{\partial {\Phi }_i}{\partial x},\frac{\partial {\Phi }_i}{\partial t}\right)=0,{\rm I}\left({\Phi }_i,\frac{\partial {\Phi }_i}{\partial x},\frac{\partial {\Phi }_i}{\partial t}\right)=0(3)$

where $q$ is the embedding parameter, ${\Phi }_i\left(x,t;q\right)$ are the representations of solutions across $q$, $y_{i,0}\left(x,t\right)$ are the initial approximations, ${ħ}_i$ are the auxiliary parameters, $H_i\left(x,t\right)$ are the auxiliary functions, and ${\mathcal{L}}_i$ and ${\mathcal{N}}_i$ are the linear and non-linear operators, respectively. Eq. 2 shows that as $q=0$, ${\Phi }_i\left(x,t;0\right)=y_{i,0}\left(x,t\right)$, and as $q=1$, ${\Phi }_i\left(x,t;1\right)=y_i\left(x,t\right)$. This means as $q$ varies from 0 to 1, ${\Phi }_i\left(x,t;q\right)$ transforms from the initial approximation to the final solution. The higher-order terms can be calculated from the higher-order deformation equation given as follows (Liao 2003):

${\mathcal{L}}_i\left[y_{i,m}\left(x,t\right)-{\chi }_m{y}_{i,m-1}\left(x,t\right)\right]={\mathrm{ħ}}_i{H}_i\left(x,t\right)R_{i,m}\left({\overrightarrow{y}}_{i,m-1}\right),m=\mathrm{1,2,3},...(4)$

where

${\chi }_m=\left\{\begin{array}{c}0\mathrm{w}\mathrm{h}\mathrm{e}\mathrm{n}m=1,\\ 1\mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}\end{array}\right.(5)$

and

$R_{i,m}\left({\overrightarrow{y}}_{i,m-1}\right)={\left.\frac{1}{\left(m-1\right)!}\frac{{\partial }^{m-1}{\mathcal{N}}_i\left[{\Phi }_i\left(x,t;q\right)\right]}{{\partial q}^{m-1}}\right|}_{q=0}(6)$

where $y_{i,m}$ for $m\ge 1$ are the higher-order terms. Eq. 4 is obtained by differentiating the zeroth-order deformation Eq. 2 in succession. Also, it is interesting to note from Eq. 4 that the original governing equations are transformed into an infinite system of linear equations, which are easier to solve. This idea is a central concept for all perturbation-based approaches. Now, the final solutions can be obtained as follows:

$y_i\left(x,t\right)=y_{i,0}\left(x,t\right)+{\displaystyle \sum }_{m=1}^{\infty }{y}_{i,m}\left(x,t\right)(7)$

One can truncate the series Eq. 7 to obtain an approximate solution to a non-linear system. It can be seen that the method involves several operators and functions, which should be adequately chosen in order to have a convergent series. For that purpose, Liao (2003) proposed three generalized rules, namely rule of solution expressions, rule of coefficient ergodicity, and rule of solution existence. These rules will be discussed in a later section in relation to the problem under consideration.

2.2 Homotopy perturbation method

Here, we rewrite a system of PDEs as follows:

${\mathcal{N}}_i\left(y_i\left(x,t\right)\right)=f\left(x,t\right)(8)$

Now, we can construct a homotopy that satisfies (He 1999):

$\left(1-q\right)\left[{\mathcal{L}}_i\left({\Phi }_i\left(x,t;q\right)\right)-{\mathcal{L}}_i\left(y_{i,0}\left(x,t\right)\right)\right]+q\left[{\mathcal{N}}_i\left[{\Phi }_i\left(x,t;q\right)\right]-f\left(x,t\right)\right]=0,i=\mathrm{1,2},\dots ,n(9)$

where the symbols denote the same variables as mentioned in the previous section. Also, similar to HAM, Eq. 9 shows that at $q=0$, ${\Phi }_i\left(x,t;0\right)=y_{i,0}\left(x,t\right)$, and at $q=1$, ${\Phi }_i\left(x,t;1\right)=y_i\left(x,t\right)$. Let us now express ${\Phi }_i\left(x,t;q\right)$ as a series in terms of $q$ as follows:

${\Phi }_i\left(x,t;q\right)={\Phi }_{i,0}+q{\Phi }_{i,1}+q^2{\Phi }_{i,2}+q^3{\Phi }_{i,3}+q^4{\Phi }_{i,4}+...(10)$

where ${\Phi }_{i,m}$ for $m\ge 1$ are the higher-order terms. As $q\to 1$, Eq. 10 produces the final solutions as:

$y_i\left(x,t\right)=\underset{q\to 1}{\lim }{\Phi }_i\left(x,t;q\right)={\displaystyle \sum }_{k=0}^{\infty }{\Phi }_{i,k}(11)$

In this method, similar to the classical perturbation approach, one can substitute the series (10) in Eq. 9 and then equate the like powers of $q$ to obtain the terms of the series ${\Phi }_{i,k}$. The difference between HAM and HPM is evident from the construction of the zeroth-order deformation equation. HAM contains an additional auxiliary function and auxiliary parameters, which help obtain a rapid convergence series (Liao, 2003; Liao, 2012). Interestingly, subject to the same set of linear and non-linear operators and unit auxiliary function, HPM solution is a subcase of HAM for ${ħ}_i=-1$.

2.3 Optimal homotopy asymptotic method

In order to obtain an accurate approximate solution with just a few terms of the series, Marinca and Herişanu (2008) proposed a variant of homotopy-based method, known as the optimal homotopy asymptotic method (OHAM), using asymptotic expansions of the functions and operators. We consider a non-linear system of PDEs as:

${\mathcal{L}}_i\left[y_i\left(x,t\right)\right]+{\mathcal{N}}_i\left[y_i\left(x,t\right)\right]+h_i\left(x,t\right)=0,i=\mathrm{1,2},\dots ,n(12)$

subject to the initial and boundary conditions:

$\mathcal{B}\left(y_i,\frac{\partial y_i}{\partial x},\frac{\partial y_i}{\partial t}\right)=0,{\rm I}\left(y_i,\frac{\partial y_i}{\partial x},\frac{\partial y_i}{\partial t}\right)=0(13)$

where symbols denote the same variables as discussed in the previous section. Following HAM, one can construct a family of equations as follows:

$\left(1-q\right)\left[{\mathcal{L}}_i\left({\Phi }_i\left(\mathrm{x},\mathrm{t},C_{i,j};q\right)\right)+h_i\left(x,t\right)\right]=H_i\left(x,\mathrm{t},C_{i,j};q\right)\left[{\mathcal{L}}_i\left({\Phi }_i\left(\mathrm{x},\mathrm{t},C_{i,j};q\right)\right)+h_i\left(x,t\right)+{\mathcal{N}}_i\left({\Phi }_i\left(\mathrm{x},\mathrm{t},C_{i,j};q\right)\right)\right](14)$

where symbols have their usual meaning and $C_{i,j}$ are the unknown parameters. The auxiliary functions are defined as:

$H_i\left(x,\mathrm{t},C_{i,j};q\right)=0\mathrm{f}\mathrm{o}\mathrm{r}q=0\left.\ne 0\mathrm{f}\mathrm{o}\mathrm{r}q\in (\mathrm{0,1}\right](15)$

It is easy to verify that when $q=0$, ${\Phi }_i\left(\mathrm{x},\mathrm{t},C_{i,j};q\right)=y_{i,0}\left(x,t\right)$ and when $q=1$, ${\Phi }_i\left(\mathrm{x},\mathrm{t},C_{i,j};q\right)=y_i\left(x,t\right)$, which is the same as HAM and HPM, i.e., as $q$ goes from 0 to 1, we have the continuous deformation from the initial approximation to the final solution. Now, the initial approximations $y_{i,0}\left(x,t\right)$ can be found by solving the following set of equations, which is obtained after substituting $q=0$ in Eq. 14:

${\mathcal{L}}_i\left(y_{i,0}\left(x,t\right)\right)+h_i\left(x,t\right)=0(16)$

subject to the boundary conditions:

$\mathcal{B}\left(y_{i,0},\frac{\partial y_{i,0}}{\partial x},\frac{\partial y_{i,0}}{\partial t}\right)=0,{\rm I}\left(y_{i,0},\frac{\partial y_{i,0}}{\partial x},\frac{\partial y_{i,0}}{\partial t}\right)=0(17)$

In OHAM, expanding the auxiliary function with respect to the embedding parameter is the key step. The auxiliary function $H_i\left(x,\mathrm{t},C_{i,j};q\right)$ can be expressed as follows:

$H_i\left(x,\mathrm{t},C_{i,j};q\right)=q{H}_{i,1}\left(x,\mathrm{t},C_{i,j}\right)+q^2{H}_{i,2}\left(x,\mathrm{t},C_{i,j}\right)+q^3{H}_{i,3}\left(x,\mathrm{t},C_{i,j}\right)+...(18)$

where $H_i\left(x,\mathrm{t},C_{i,j}\right)$’s are the auxiliary functions that depend on the independent variables $x$ and $t$ and parameters $C_{i,j}$. It can be noted that the series Eq. 18 is in accordance with the property Eq. (15). Now, the solution of Eq. 12 can be assumed to be of the form:

${\Phi }_i\left(\mathrm{x},\mathrm{t},C_{i,j};q\right)=y_{i,0}\left(x,t\right)+{\displaystyle \sum }_{j=1}^{\infty }{y}_{i,j}\left(x,t,C_{i,j}\right)q^j(19)$

Substituting Eq. 19 into Eq. 14, and equating the like powers of $q$, the following equations are obtained [$q^0$ corresponds to Eqs. 16, 17:

${\mathcal{L}}_i\left(y_{i,1}\left(x,t,C_{i,j}\right)\right)=H_{i,1}\left(x,t,C_{i,j}\right){\mathcal{N}}_{i,0}\left(y_{i,0}\left(x,t\right)\right)(20)$

$\mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}$ the initial and boundary conditions:

$\mathcal{B}\left(y_{i,1},\frac{\partial y_{i,1}}{\partial x},\frac{\partial y_{i,1}}{\partial t}\right)=0,{\rm I}\left(y_{i,1},\frac{\partial y_{i,1}}{\partial x},\frac{\partial y_{i,1}}{\partial t}\right)=0(21)$

For $k=\mathrm{2,3,4},...$,

$\begin{array}{c}{\mathcal{L}}_i\left[y_{i,k}\left(x,t,C_{i,j}\right)-y_{i,k-1}\left(x,t,C_{i,j}\right)\right]=H_{i,k}\left(x,t,C_{i,j}\right){\mathcal{N}}_{i,0}\left(y_{i,0}\left(x,t\right)\right)\\ +{\displaystyle \sum }_{j=1}^{k-1}{H}_{i,j}\left({x,t,C}_{i,j}\right)[{\mathcal{L}}_i\left[y_{i,k-j}\left(x,t,C_{i,j}\right)\right]\\ +{\mathcal{N}}_{i,k-j}[y_{i,0}\left(x,t\right),y_{i,1}\left(x,t,C_{i,j}\right),\dots ,y_{i,k-j}\left(x,t,C_{i,j}\right)]]\end{array}(22)$

subject to the initial and boundary conditions:

$\mathcal{B}\left(y_{i,k},\frac{\partial y_{i,k}}{\partial x},\frac{\partial y_{i,k}}{\partial t}\right)=0,{\rm I}\left(y_{i,k},\frac{\partial y_{i,k}}{\partial x},\frac{\partial y_{i,k}}{\partial t}\right)=0(23)$

where the term ${\mathcal{N}}_{i,k-j}\left[y_{i,0}\left(x,t\right),y_{i,1}\left(x,t,C_{i,j}\right),\dots ,y_{i,k-j}\left(x,t,C_{i,j}\right)\right]$ is the coefficient of $q^m$, which is obtained by expanding ${\mathcal{N}}_i\left({\Phi }_i\left(\mathrm{x},\mathrm{t},C_{i,j};q\right)\right)$ as follows:

${\mathcal{N}}_i\left({\Phi }_i\left(\mathrm{x},\mathrm{t},C_{i,j};q\right)\right)={\mathcal{N}}_{i,0}\left(y_{i,0}\left(x,t\right)\right)+q{\mathcal{N}}_{i,1}\left(y_{i,0}\left(x,t\right),y_{i,1}\left(x,t,C_{i,j}\right)\right)+q^2{\mathcal{N}}_{i,2}\left(y_{i,0}\left(x,t\right),y_{i,1}\left(x,t,C_{i,j}\right),y_{i,2}\left(x,t,C_{i,j}\right)\right)+...(24)$

It can be observed from Eq. 22 that like HAM and HPM, OHAM also converts the original non-linear equation into an infinite set of linear sub-equations. Further, the method does not depend on the presence of a small parameter in the governing equation. The convergence of the series Eq. 19 depends on the choice of $H_{i,j}\left({x,t,C}_{i,j}\right)$, and there exist many ways to choose it. According to Marinca and Herişanu (2015), one can choose $H_{i,j}\left({x,t,C}_{i,j}\right)$ in such a way that the product $H_{i,j}\left({x,t,C}_{i,j}\right)\left[{\mathcal{L}}_i\left[y_{i,k-j}\left(x,t,C_{i,j}\right)\right]+{\mathcal{N}}_{i,k-j}\left[y_{i,0}\left(x,t\right),y_{i,1}\left(x,t,C_{i,j}\right),\dots ,y_{i,k-j}\left(x,t,C_{i,j}\right)\right]\right]$ of Eq. 22 is of the same form as that of $H_{i,j}\left({x,t,C}_{i,j}\right)$. The considered functions can be of any type, such as polynomial, exponential, trigonometric, etc. Now, based on the choice of auxiliary function, if the series Eq. 19 converges at $q=1$, then

$y_i\left(x,t,C_{i,j}\right)=y_{i,0}\left(x,t\right)+{\displaystyle \sum }_{j=1}^{\infty }{y}_{i,j}\left(x,t,C_{i,j}\right)(25)$

Finally, the approximate solution can be obtained by considering a finite number of terms of the series Eq. 25. The unknown parameters $C_{i,j}$ and the choice of auxiliary function will be discussed later in relation to the problem under consideration.

3 Governing equation and solution methodologies

In a one-dimensional flow field, the general form for a reactive transport system which describes the phenomenon of multi-species migration having spatially and temporally varying parameters can be modelled as follows (Chaudhary and Singh 2020):

$r_i\frac{\partial C_i}{\partial t}=\frac{\partial }{\partial x}\left(D\left(x,t\right)\frac{\partial C_i}{\partial x}-v\left(x,t\right)C_i\right)+k_{i-1}{C}_{i-1}-k_i{C}_i,i=1,2,...,n(26)$

Here, $k_0=0$ and $n$ represents the number of species; $x$ and $t$ denote the spatial and temporal coordinates, respectively; $r_i$ is the retardation factor for the $i$-th species; $C_i$ is the concentration strength for the $i$-th species; $k_i$ is the decay rate coefficient for the $i$-th species; and $D\left(x,t\right)$ and $v\left(x,t\right)$ are the dispersion coefficient and seepage velocity, respectively.

The present work considers two-species system (i.e., $i=1,2$) and includes the effect of spatially and temporally dependent transport parameters on pollutant migration. To that end, the model becomes (Chaudhary and Singh 2020):

$r_1\frac{\partial C_1}{\partial t}=\frac{\partial }{\partial x}\left(D\left(x,t\right)\frac{\partial C_1}{\partial x}-v\left(x,t\right)C_1\right)-k_1{C}_1(27)$

$r_2\frac{\partial C_2}{\partial t}=\frac{\partial }{\partial x}\left(D\left(x,t\right)\frac{\partial C_2}{\partial x}-v\left(x,t\right)C_2\right)-k_2{C}_2+k_1{C}_1(28)$

where $C_1$ and $C_2$ represent the pollutant concentration level for parent and daughter species, respectively.

The initial and boundary conditions for the model can be given as follows:

$C_1\left(x,0\right)=f_1\left(x\right),C_1\left(0,t\right)=f_1\left(0\right),\underset{x\to \infty }{\lim }{C}_1\left(x,t\right)=0(29)$

$C_2\left(x,0\right)=f_2\left(x\right),C_2\left(0,t\right)=f_2\left(0\right),\underset{x\to \infty }{\lim }{C}_2\left(x,t\right)=0(30)$

where $f_1\left(x\right)$ and $f_2\left(x\right)$ are considered as (Simpson and Ellery 2014):

$f_1\left(x\right)=a_{1}xexp\left(-a_{2}x\right),f_2\left(x\right)=0(31)$

Eqs 29, 30, 31 show that spatially dependent pollution exists initially for the parent species. However, for the daughter species, the initial concentration is zero.

Based on the forms of $D\left(x,t\right)$ and $v\left(x,t\right)$, several scenarios can be considered. Here, we consider a specific case, specifically the steady migration phenomenon in the case of steady groundwater flow. It may be noted that the solution methodologies reported here can be adopted for other cases such as the transient migration in steady and unsteady flows (Chaudhary and Singh 2020). However, as our objective is to show the applicability of different homotopy-based methods by extending the work of Chaudhary and Singh (2020), we restrict our analysis for the case of steady migration in steady flow. To that end, we assume:

$D\left(x,t\right)=D_0{\left({\alpha }_1+{\alpha }_{2}x\right)}^2,v\left(x,t\right)=v_0\left({\alpha }_1+{\alpha }_{2}x\right)(32)$

where $D_0$ and $v_0$ are the initial dispersion and velocity, respectively; ${\alpha }_1$ and ${\alpha }_2$ are the parameters, where ${\alpha }_2$ is knows as the heterogeneity parameter. Eq. 32 shows that dispersion is directly proportional to the square of velocity (Batu 2005). For ${\alpha }_1=1$ and ${\alpha }_2\to 0$, i.e., ${\alpha }_2$ is very small, the system converts into a homogeneous one, where the effect of space on dispersion and seepage velocity is absent. As ${\alpha }_2$ increases, the system becomes heterogeneous. It may be noted that the dispersion can be expressed differently but Eq. 32 is reasonable for purposes of this study. Further, seepage velocity given in Eq. 32 can also be expressed exponentially or in power form. However, for purposes of this study, Eq. 32 is reasonable, especially for homogeneous aquifers. Before applying the homotopy-based methods, for convenience, let us rewrite the governing Eqs 27, 28 in the following form:

$r_1\frac{\partial C_1}{\partial t}=D\frac{{\partial }^2{C}_1}{\partial x^2}+\left(\frac{\partial D}{\partial x}-v\right)\frac{\partial C_1}{\partial x}-\left(k_1+\frac{\partial v}{\partial x}\right)C_1(33)$

$r_2\frac{\partial C_2}{\partial t}=D\frac{{\partial }^2{C}_2}{\partial x^2}+\left(\frac{\partial D}{\partial x}-v\right)\frac{\partial C_2}{\partial x}-\left(k_2+\frac{\partial v}{\partial x}\right)C_2+k_1{C}_1(34)$

3.1 HAM-based solution

As discussed in Section 2.1, one can apply the HAM to the system of Eqs 33, 34 by constructing some operators and functions that are given below. The non-linear operators for the problem are selected as:

${\mathcal{N}}_1\left[{\Phi }_i\left(x,t;q\right)\right]=r_1\frac{\partial {\Phi }_1\left(x,t;q\right)}{\partial t}-D\frac{{\partial }^2{\Phi }_1\left(x,t;q\right)}{\partial x^2}-\left(\frac{\partial D}{\partial x}-v\right)\frac{\partial {\Phi }_1\left(x,t;q\right)}{\partial x}+\left(k_1+\frac{\partial v}{\partial x}\right){\Phi }_1\left(x,t;q\right)(35)$

${\mathcal{N}}_2\left[{\Phi }_i\left(x,t;q\right)\right]=r_2\frac{\partial {\Phi }_2\left(x,t;q\right)}{\partial t}-D\frac{{\partial }^2{\Phi }_2\left(x,t;q\right)}{\partial x^2}-\left(\frac{\partial D}{\partial x}-v\right)\frac{\partial {\Phi }_2\left(x,t;q\right)}{\partial x}+\left(k_2+\frac{\partial v}{\partial x}\right){\Phi }_2\left(x,t;q\right)-k_1{\Phi }_1\left(x,t;q\right)(36)$

The above equations are the original operators of the equations, which is always an easy and convenient consideration for the non-linear operators in the framework of HAM. Therefore, using Eqs 35, 36, terms $R_{i,m}$ can be calculated from Eq. 6 as follows:

$R_{1,m}\left({\overrightarrow{C}}_{i,m-1}\right)=r_1\frac{\partial C_{1,m-1}}{\partial t}-D\frac{{\partial }^2{C}_{1,m-1}}{\partial x^2}-\left(\frac{\partial D}{\partial x}-v\right)\frac{\partial C_{1,m-1}}{\partial x}+\left(k_1+\frac{\partial v}{\partial x}\right)C_{1,m-1}(37)$

$R_{2,m}\left({\overrightarrow{C}}_{i,m-1}\right)=r_2\frac{\partial C_{2,m-1}}{\partial t}-D\frac{{\partial }^2{C}_{2,m-1}}{\partial x^2}-\left(\frac{\partial D}{\partial x}-v\right)\frac{\partial C_{2,m-1}}{\partial x}+\left(k_2+\frac{\partial v}{\partial x}\right)C_{2,m-1}-k_1{C}_{1,m-1}(38)$

Now, as discussed earlier, we need to follow three fundamental rules provided by Liao (2003) to have a convergent series solution. First, we consider the following set of base functions to represent the solutions $C_i\left(x,t\right)$ of the present problem:

$\left\{x^m{t}^p\mathit{exp}\left(-nbx\right)|m,n,p=\mathrm{0,1,2},3,...\right\}(39)$

so that

$C_i\left(x,t\right)={\displaystyle \sum }_{p=0}^{\infty }{\displaystyle \sum }_{m=0}^{\infty }{\displaystyle \sum }_{n=1}^{\infty }{\beta }_{m,n,p}{x}^m{t}^p\mathit{exp}\left(-nbx\right)(40)$

where ${\beta }_{m,n,p}$ are the coefficients of the series. Eq. 40 provides the so-called rule of solution expression. Following the rule of solution expression, the linear operators and the initial approximations are chosen, respectively, as follows:

${\mathcal{L}}_i\left[{\Phi }_i\left(x,t;q\right)\right]=r_i\frac{\partial {\Phi }_i\left(x,t;q\right)}{\partial t}\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{t}\mathrm{y}{\mathcal{L}}_i\left[E_2\right]=0(41)$

$C_{\mathrm{1,0}}\left(x,t\right)=f_1\left(x\right)=a_{1}x\exp \left(-a_{2}x\right),C_{\mathrm{2,0}}\left(x,t\right)=f_2\left(x\right)=0(42)$

where $E_2$ is an integral constant. It may be noted that the choice for Eq. 41 is not unique. Using Eq. 41, the higher-order terms can be obtained from Eq. 4 as follows:

$C_{i,m}\left(x,t\right)={\chi }_m{C}_{i,m-1}\left(x,t\right)+\left(\frac{{\hslash }_i}{r_i}\right){\int }_0^t{H_i\left(x,t\right)R}_{i,m}\left({\overrightarrow{C}}_{i,m-1}\right)dt+E_2,m=\mathrm{1,2,3},...(43)$

where $R_{i,m}$ are given by Eqs 37, 38, and constant $E_2$ can be determined from the initial condition for the higher-order deformation equations, which simply yields $E_2=0$ for all $m\ge 1$.

Now, the auxiliary functions $H_i\left(x,t\right)$ can be determined from the rule of coefficient ergodicity. Based on the rule of solution expression, the general form of $H_i\left(x,t\right)$ should be:

$H_i\left(x,t\right)=x^{n_1}{t}^{n_3}\exp \left(-b{n}_{2}x\right)(44)$

where $n_1$, $n_2$, and $n_3$ are the integers. However, for simplicity, we can take $H_i\left(x,t\right)=1$ (Vajravelu and Van Gorder 2013). Finally, the approximate solutions can be obtained as follows:

$C_i\left(x,t\right)\approx C_{i,0}\left(x,t\right)+{\displaystyle \sum }_{m=1}^M{C}_{i,m}\left(x,t\right)(45)$

For the convenience of readers, a flowchart containing the steps of the method is provided in Figure 1.

FIGURE 1

FIGURE 1. Flowchart for the HAM solution (45).

3.2 HPM-based solution

In relation to the discussion in Section 2.2, for simplicity, we consider the linear operator as follows:

${\mathcal{L}}_i\left[{\Phi }_i\left(x,t;q\right)\right]=r_i\frac{\partial {\Phi }_i\left(x,t;q\right)}{\partial t},i=\mathrm{1,2}(46)$

The non-linear operators are selected as Eqs 35, 36. In the framework of HPM also, the selection of these operators is not unique; indeed, one can choose any other forms to check for a better solution. Using these expressions in Eq. 9 and then substituting the series Eq. 10, we obtain the following systems of differential equations after equating the like powers of $q$:

$r_1\left(\frac{\partial {\Phi }_{\mathrm{1,0}}}{\partial t}-\frac{\partial C_{\mathrm{1,0}}}{\partial t}\right)=0\mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}{\Phi }_{\mathrm{1,0}}\left(x,0\right)=f_1\left(x\right),{\Phi }_{\mathrm{1,0}}\left(0,t\right)=f_1\left(0\right),\underset{x\to \infty }{\lim }{\Phi }_{\mathrm{1,0}}\left(x,t\right)=0(47)$

$r_2\left(\frac{\partial {\Phi }_{\mathrm{2,0}}}{\partial t}-\frac{\partial C_{\mathrm{2,0}}}{\partial t}\right)=0\mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}{\Phi }_{\mathrm{2,0}}\left(x,0\right)=f_2\left(x\right),{\Phi }_{\mathrm{2,0}}\left(0,t\right)=f_2\left(0\right),\underset{x\to \infty }{\lim }{\Phi }_{\mathrm{2,0}}\left(x,t\right)=0(48)$

$\begin{array}{c}{r}_1\left(\frac{\partial {\Phi }_{\mathrm{1,1}}}{\partial t}-\left(\frac{\partial {\Phi }_{\mathrm{1,0}}}{\partial t}-\frac{\partial C_{\mathrm{1,0}}}{\partial t}\right)\right)+r_1\frac{\partial {\Phi }_{\mathrm{1,0}}}{\partial t}-D\frac{{\partial }^2{\Phi }_{\mathrm{1,0}}}{\partial x^2}-\left(\frac{\partial D}{\partial x}-v\right)\frac{\partial {\Phi }_{\mathrm{1,0}}}{\partial x}+\left(k_1+\frac{\partial v}{\partial x}\right){\Phi }_{\mathrm{1,0}}\\ =0\mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}{\Phi }_{\mathrm{1,1}}\left(x,0\right)=0,{\Phi }_{\mathrm{1,1}}\left(0,t\right)=0,\underset{x\to \infty }{\lim }{\Phi }_{\mathrm{1,1}}\left(x,t\right)=0\end{array}(49)$

$\begin{array}{c}{r}_2\left(\frac{\partial {\Phi }_{\mathrm{2,1}}}{\partial t}-\left(\frac{\partial {\Phi }_{\mathrm{2,0}}}{\partial t}-\frac{\partial C_{\mathrm{2,0}}}{\partial t}\right)\right)+r_2\frac{\partial {\Phi }_{\mathrm{2,0}}}{\partial t}-D\frac{{\partial }^2{\Phi }_{\mathrm{2,0}}}{\partial x^2}-\left(\frac{\partial D}{\partial x}-v\right)\frac{\partial {\Phi }_{\mathrm{2,0}}}{\partial x}+\left(k_2+\frac{\partial v}{\partial x}\right){\Phi }_{\mathrm{2,0}}-k_1{\Phi }_{\mathrm{1,0}}\\ =0\mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}{\Phi }_{\mathrm{2,1}}\left(x,0\right)=0,{\Phi }_{\mathrm{2,1}}\left(0,t\right)=0,\underset{x\to \infty }{\lim }{\Phi }_{\mathrm{2,1}}\left(x,t\right)=0\end{array}(50)$

$r_1\left(\frac{\partial {\Phi }_{\mathrm{1,2}}}{\partial t}-\frac{\partial {\Phi }_{\mathrm{1,1}}}{\partial t}\right)+r_1\frac{\partial {\Phi }_{\mathrm{1,1}}}{\partial t}-D\frac{{\partial }^2{\Phi }_{\mathrm{1,1}}}{\partial x^2}-\left(\frac{\partial D}{\partial x}-v\right)\frac{\partial {\Phi }_{\mathrm{1,1}}}{\partial x}+\left(k_1+\frac{\partial v}{\partial x}\right){\Phi }_{\mathrm{1,1}}=0\mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}{\Phi }_{\mathrm{1,2}}\left(x,0\right)=0,{\Phi }_{\mathrm{1,2}}\left(0,t\right)=0,\underset{x\to \infty }{\lim }{\Phi }_{\mathrm{1,2}}\left(x,t\right)=0(51)$

$\begin{array}{c}{r}_2\left(\frac{\partial {\Phi }_{\mathrm{2,2}}}{\partial t}-\frac{\partial {\Phi }_{\mathrm{2,1}}}{\partial t}\right)+r_2\frac{\partial {\Phi }_{\mathrm{2,1}}}{\partial t}-D\frac{{\partial }^2{\Phi }_{\mathrm{2,1}}}{\partial x^2}-\left(\frac{\partial D}{\partial x}-v\right)\frac{\partial {\Phi }_{\mathrm{2,1}}}{\partial x}+\left(k_2+\frac{\partial v}{\partial x}\right){\Phi }_{\mathrm{2,1}}-k_1{\Phi }_{\mathrm{1,1}}\\ =0\mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}{\Phi }_{\mathrm{2,2}}\left(x,0\right)=0,{\Phi }_{\mathrm{2,2}}\left(0,t\right)=0,\underset{x\to \infty }{\lim }{\Phi }_{\mathrm{2,2}}\left(x,t\right)=0\end{array}(52)$

Proceeding in a like manner, one can arrive at the following recurrence relation:

$\begin{array}{c}{r}_1\left(\frac{\partial {\Phi }_{1,m}}{\partial t}-\frac{\partial {\Phi }_{1,m-1}}{\partial t}\right)+r_1\frac{\partial {\Phi }_{1,m-1}}{\partial t}-D\frac{{\partial }^2{\Phi }_{1,m-1}}{\partial x^2}-\left(\frac{\partial D}{\partial x}-v\right)\frac{\partial {\Phi }_{1,m-1}}{\partial x}+\left(k_1+\frac{\partial v}{\partial x}\right){\Phi }_{1,m-1}\\ =0\mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}{\Phi }_{1,m}\left(x,0\right)=0,{\Phi }_{1,m}\left(0,t\right)=0,\underset{x\to \infty }{\lim }{\Phi }_{1,m}\left(x,t\right)=0\mathrm{f}\mathrm{o}\mathrm{r}m\ge 2\end{array}(53)$

$\begin{array}{c}{r}_2\left(\frac{\partial {\Phi }_{2,m}}{\partial t}-\frac{\partial {\Phi }_{2,m-1}}{\partial t}\right)+r_2\frac{\partial {\Phi }_{2,m-1}}{\partial t}-D\frac{{\partial }^2{\Phi }_{2,m-1}}{\partial x^2}-\left(\frac{\partial D}{\partial x}-v\right)\frac{\partial {\Phi }_{2,m-1}}{\partial x}+\left(k_2+\frac{\partial v}{\partial x}\right){\Phi }_{2,m-1}-k_1{\Phi }_{1,m-1}\\ =0\mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}{\Phi }_{2,m}\left(x,0\right)=0,{\Phi }_{2,m}\left(0,t\right)=0,\underset{x\to \infty }{\lim }{\Phi }_{2,m}\left(x,t\right)=0\mathrm{f}\mathrm{o}\mathrm{r}m\ge 2\end{array}(54)$

The initial approximation can be chosen as ${\Phi }_{\mathrm{1,0}}=f_1\left(x\right)=a_{1}xexp\left(-a_{2}x\right)$ and ${\Phi }_{\mathrm{2,0}}=f_2\left(x\right)=0$. Using this initial approximation, we can solve the equations iteratively using a symbolic software. Finally, the HPM-based solutions can be approximated as follows:

$C_i\left(x,t\right)\approx {\displaystyle \sum }_{k=0}^M{\Phi }_{i,k}(55)$

A flowchart containing the steps of the method is given in Figure 2.

FIGURE 2

FIGURE 2. Flowchart for the HPM solution (55).

3.3 OHAM-based solution

In relation to the discussion given in Section 2.3, we choose the following operators:

${\mathcal{L}}_i\left[C_i\left(x,t\right)\right]=r_i\frac{\partial C_i\left(x,t\right)}{\partial t}(56)$

${\mathcal{N}}_1\left[C_i\left(x,t\right)\right]=-D\frac{{\partial }^2{C}_1\left(x,t\right)}{\partial x^2}-\left(\frac{\partial D}{\partial x}-v\right)\frac{\partial C_1\left(x,t\right)}{\partial x}+\left(k_1+\frac{\partial v}{\partial x}\right)C_1\left(x,t\right)(57)$

${\mathcal{N}}_2\left[C_i\left(x,t\right)\right]=-D\frac{{\partial }^2{C}_2\left(x,t\right)}{\partial x^2}-\left(\frac{\partial D}{\partial x}-v\right)\frac{\partial C_2\left(x,t\right)}{\partial x}+\left(k_2+\frac{\partial v}{\partial x}\right)C_2\left(x,t\right)-k_1{C}_1\left(x,t\right)(58)$

and

$h_i\left(x,t\right)=0(59)$

With these considerations, we can solve the zeroth-order Eq. 16 to have the following solution:

$C_{i,0}\left(x,t\right)=f_i\left(x\right)(60)$

Now to have the higher-order equations, we need the expressions ${\mathcal{N}}_{i,0}$, ${\mathcal{N}}_{i,1}$, ${\mathcal{N}}_{i,2}$, etc. as can be seen from Eq. 20. For that, one can use Eq. 24 in relation with Eqs 57, 58. Using this we have the first-order equations:

$r_1\frac{\partial C_{\mathrm{1,1}}}{\partial t}=H_{\mathrm{1,1}}\left(x,t,D_{1,i}\right)\left[-D\frac{{\partial }^2{C}_{\mathrm{1,0}}}{\partial x^2}-\left(\frac{\partial D}{\partial x}-v\right)\frac{\partial C_{\mathrm{1,0}}}{\partial x}+\left(k_1+\frac{\partial v}{\partial x}\right)C_{\mathrm{1,0}}\right]\mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}{C}_{\mathrm{1,1}}\left(x,0\right)=0,C_{\mathrm{1,1}}\left(0,t\right)=0,\underset{x\to \infty }{\lim }{C}_{\mathrm{1,1}}\left(x,t\right)=0(61)$

$r_2\frac{\partial C_{\mathrm{2,1}}}{\partial t}=H_{\mathrm{2,1}}\left(x,t,D_{2,i}\right)\left[-D\frac{{\partial }^2{C}_{\mathrm{2,0}}}{\partial x^2}-\left(\frac{\partial D}{\partial x}-v\right)\frac{\partial C_{\mathrm{2,0}}}{\partial x}+\left(k_2+\frac{\partial v}{\partial x}\right)C_{\mathrm{2,0}}-k_1{C}_{\mathrm{1,0}}\right]\mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}{C}_{\mathrm{2,1}}\left(x,0\right)=0,C_{\mathrm{2,1}}\left(0,t\right)=0,\underset{x\to \infty }{\lim }{C}_{\mathrm{2,1}}\left(x,t\right)=0(62)$

The auxiliary functions can be chosen in many ways. Here, we select $H_{\mathrm{1,1}}\left(x,t,D_{1,i}\right)=D_{\mathrm{1,1}}$ and $H_{\mathrm{2,1}}\left(x,t,D_{2,i}\right)=D_{\mathrm{2,1}}.$ Putting $k=2$ in Eq. 22 and then using the expansion of $\mathcal{N}$, we have the following second-order equation:

$\begin{array}{c}{r}_1\frac{\partial C_{\mathrm{1,2}}}{\partial t}=r_1\frac{\partial C_{\mathrm{1,1}}}{\partial t}+D_{\mathrm{1,2}}\left[-D\frac{{\partial }^2{C}_{\mathrm{1,0}}}{\partial x^2}-\left(\frac{\partial D}{\partial x}-v\right)\frac{\partial C_{\mathrm{1,0}}}{\partial x}+\left(k_1+\frac{\partial v}{\partial x}\right)C_{\mathrm{1,0}}\right]\\ +D_{\mathrm{1,1}}\left[r_1\frac{\partial C_{\mathrm{1,1}}}{\partial t}-D\frac{{\partial }^2{C}_{\mathrm{1,1}}}{\partial x^2}-\left(\frac{\partial D}{\partial x}-v\right)\frac{\partial C_{\mathrm{1,1}}}{\partial x}+\left(k_1+\frac{\partial v}{\partial x}\right)C_{\mathrm{1,1}}\right]\mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}{C}_{\mathrm{1,2}}\left(x,0\right)\\ =0,C_{\mathrm{1,2}}\left(0,t\right)=0,\underset{x\to \infty }{\lim }{C}_{\mathrm{1,2}}\left(x,t\right)=0\end{array}(63)$

$\begin{array}{c}{r}_2\frac{\partial C_{\mathrm{2,2}}}{\partial t}=r_2\frac{\partial C_{\mathrm{2,1}}}{\partial t}+D_{\mathrm{2,2}}\left[-D\frac{{\partial }^2{C}_{\mathrm{2,0}}}{\partial x^2}-\left(\frac{\partial D}{\partial x}-v\right)\frac{\partial C_{\mathrm{2,0}}}{\partial x}+\left(k_2+\frac{\partial v}{\partial x}\right)C_{\mathrm{2,0}}-k_1{C}_{\mathrm{1,0}}\right]\\ +D_{\mathrm{2,1}}\left[r_2\frac{\partial C_{\mathrm{2,1}}}{\partial t}-D\frac{{\partial }^2{C}_{\mathrm{2,1}}}{\partial x^2}-\left(\frac{\partial D}{\partial x}-v\right)\frac{\partial C_{\mathrm{2,1}}}{\partial x}+\left(k_2+\frac{\partial v}{\partial x}\right)C_{\mathrm{2,1}}-k_1{C}_{\mathrm{1,1}}\right]\\ \mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}{C}_{\mathrm{1,2}}\left(x,0\right)=0,C_{\mathrm{1,2}}\left(0,t\right)=0,\underset{x\to \infty }{\lim }{C}_{\mathrm{1,2}}\left(x,t\right)=0\end{array}(64)$

Here, we choose $H_{\mathrm{1,2}}\left(x,t,D_{2,i}\right)=D_{\mathrm{1,2}}$ and $H_{\mathrm{2,2}}\left(x,t,D_{2,i}\right)=D_{\mathrm{2,2}}$. The terms of the series can be computed using the equations developed here. One can compute these terms without any difficulty using a symbolic computation software, such as MATLAB. Further, following Eq. 22, the higher-order terms can be computed in a similar manner. However, our aim is to produce an accurate solution with just two-three terms of OHAM-based series. Therefore, we restrict our calculation up to $k=2$. Finally, the approximate solution can be found as:

$C_i\left(x,t\right)\approx C_{i,0}\left(x,t\right)+C_{i,1}\left(x,t,D_{i,i}\right)+C_{\mathrm{1,2}}\left(x,t,D_{i,i}\right)(65)$

where the terms are given by Eqs. 60, 61, 62, 63, 64. A flowchart containing the steps of the method is provided in Figure 3.

FIGURE 3

FIGURE 3. Flowchart for the OHAM solution (65).

In Section 2, we discussed the homotopy-based methods from a general framework of system of PDEs. Indeed, the methodologies are equally applicable to single equations, integral equations, etc. Then, in Section 3, the methodologies were applied to a system of equations describing two-species pollutant transport. It can be seen from the application of these methods that they are different from the construction of different functions, operators, and parameters, and therefore may have advantages or disadvantages while dealing with a particular problem. From a theoretical perspective, the convergence theorems for HAM- and OHAM-based solutions are provided in Appendix A.

4 Results and discussion

First, we discuss the expressions and the values of parameters needed for the assessment of the developed solutions. Then, the numerical convergence of the HAM-based analytical solution is established for a specific test case, and then the solution is validated over a numerical solution. Finally, the HPM- and OHAM-based analytical solutions are validated by comparing them with the numerical solution.

4.1 Selection of expressions and parameters

For the assessment of solutions, it can be seen that the dispersion coefficient $D\left(x,t\right)$ and seepage velocity $v_0\left(x,t\right)$ are taken as given by Eq. 2. Using these equations, we have:

$\frac{\partial D}{\partial x}=2{\alpha }_2{D}_0\left({\alpha }_1+{\alpha }_{2}x\right),\frac{\partial v}{\partial x}={\alpha }_2{v}_0(66)$

For the validation of solutions, space and time domains are considered as $0\le x\left(\mathrm{c}\mathrm{m}\right)\le 3000$ and $0\le t\left(\mathrm{h}\mathrm{r}\right)\le 40$. The required parameters are chosen as (Simpson and Ellery 2014; Chaudhary and Singh 2020): $a_1=2$, $a_2=0.0025$, $k_1=0.05\left(/hr\right)$, $k_2=0.01\left(/\mathrm{h}\mathrm{r}\right)$, $v_0=0.2\left(\mathrm{c}\mathrm{m}/\mathrm{h}\mathrm{r}\right)$, and $D_0=0.5\left({\mathrm{c}\mathrm{m}}^2/\mathrm{h}\mathrm{r}\right)$. Also, the other parameters are considered as $r_1=2$, $r_2=2.5$, ${\alpha }_1=1$, and ${\alpha }_2=0.05\left(/\mathrm{c}\mathrm{m}\right).$

4.2 Numerical convergence and validation of the HAM solution

There are two ways to handle the convergence of HAM: one is based on the $\hslash -$ curves, and the other is by finding the squared residual error (Liao 2012). It can be noted that for the system of equations considered here, we have chosen ${\hslash }_1={\hslash }_2=\hslash$. One can indeed choose two different auxiliary parameters; however, it is better to first consider a common $\hslash$ to see the accuracy of solutions, as this assumption simplifies the problem. Here, we calculate the so-called $\hslash -\mathrm{c}\mathrm{u}\mathrm{r}\mathrm{v}\mathrm{e}\mathrm{s}$ to find a suitable choice for the auxiliary parameter $\hslash$ , which determines the convergence of the HAM-based series solution. In this regard, for a particular order of approximation, we plot the approximate solution $C_i\left(x,t\right)$ (or its derivatives) at some point within the domain. The flatness of the $\mathrm{\hslash }-$ curve determines a suitable choice for the auxiliary parameter $\mathrm{\hslash }$. In Figure 4, we plot the $\mathrm{\hslash }-$ curves for 5^th and 10th order HAM-based approximations for the value $C_i\left(\mathrm{1,1}\right)$. From the exact solution, those quantities can be calculated, and it is observed from the figures that the curves exhibit the flat nature for a specific range of $\mathrm{\hslash }$. Any choice of $\mathrm{\hslash }$ within this range determines an optimal value for which the series solutions converge (Abbasbandy et al., 2011).

FIGURE 4

FIGURE 4. $\mathrm{\hslash }-$ curves for $C_i\left(\mathrm{1,1}\right)$ of the HAM solution: (A) 5^th order approximation, and (B) 10th order approximation.

Here, we test the performance of the HAM-based analytical solutions for the present problem by comparing them with a numerical solution. The numerical solution is obtained using an efficient MATLAB tool, namely ‘pdepe’. This MATLAB routine uses the method of lines by discretizing a parabolic PDE (single or system) in one space direction (Skeel and Berzins 1990). For the parameters described in the previous section, the pollutant concentrations are computed using pdepde for $t=$ 0, 20, and 40 h. On the other hand, HAM solutions are computed for each of the cases. The auxiliary parameter $\hslash$ is chosen as −0.87. Figure 5 depicts the pollutant concentration values for the selected cases. It can be seen from the figure that the 3rd order HAM solution agrees very well with the numerical solution.

FIGURE 5

FIGURE 5. HAM-based pollutant concentrations versus distance for both species at: (A) $t=0$ h, (B) $t=20$ h, and (C) $t=40$ h.

4.3 Validation of HPM-based solution

The HPM-based analytical solutions for the selected test case are validated over the numerical solution obtained using pdepe of MATLAB. It is seen that the five-term HPM series produces accurate results over the selected domain. Figure 6 compares the numerical solutions and the HPM-based approximations.

FIGURE 6

FIGURE 6. HPM-based pollutant concentrations versus distance for both species at: (A) $t=0$ h, (B) $t=20$ h, and (C) $t=40$ h.

4.4 Validation of OHAM-based solution

For the assessment of the OHAM-based analytical solution, one needs to calculate the constants $D_{i,j}$’s. For that purpose, we calculate the residual as follows:

$R_i\left(x,t,D_{i,j}\right)={\mathcal{L}}_i\left[C_{i,OHAM}\left(x,t,D_{i,j}\right)\right]+{\mathcal{N}}_i\left[C_{i,OHAM}\left(x,t,D_{i,j}\right)\right]+h_i\left(x,t\right),j=\mathrm{1,2},\dots ,s(67)$

where $C_{i,OHAM}\left(x,t,D_{i,j}\right)$ is the approximate solution. When $R_i\left(x,t,D_{i,j}\right)=0$, $C_{i,OHAM}\left(x,t,D_{i,j}\right)$ becomes the exact solution to the problem. One of the ways to obtain the optimal $D_{i,j}$’s, for which the solution converges, is the minimization of squared residual error, i.e.,

$J\left(D_{i,j}\right)={\int }_{\left(x,t\right)ϵD}{\displaystyle \sum }_{i=1}^2{R_i}^2\left(x,t,D_{i,j}\right)dxdt,j=\mathrm{1,2},...,s(68)$

where $D$ is the domain of the problem. The minimization of Eq. 68 leads to a system of algebraic equations as follows:

$\frac{\partial J}{\partial D_{\mathrm{1,1}}}=\frac{\partial J}{\partial D_{\mathrm{1,2}}}=...=\frac{\partial J}{\partial D_{1,s}}=0(69)$

Solving this system, one can obtain the optimal values of the parameters. We obtain the optimal values using the MATLAB routine fminsearch, which minimizes an unconstrained multivariate function. Three-term OHAM solution is computed and compared in Figure 7 with the numerical solution obtained using pdepe to see the effectiveness of the proposed approach. It can be seen that the OHAM-based approximation agrees well with the numerical solution for all cases.

FIGURE 7

FIGURE 7. OHAM-based pollutant concentrations versus distance for both species at: (A) $t=0$ h, (B) $t=20$ h, and (C) $t=40$ h.

4.5 Comparison between different solutions

This work develops three kinds of analytical solutions for two-species pollutant transport equation using the homotopy-based methods. In this section, we validated the performances of different solutions. It is shown that all of them agree well with the corresponding numerical solutions to the non-linear problem. However, the methods have their own advantages or disadvantages. Specifically, the HAM provides a great freedom in choosing linear, non-linear operators, and auxiliary functions subject to the choice of base functions. Further, the convergence-control parameter present in HAM greatly enhances the radius of convergence of the series and monitors the rate of convergence. On the other hand, the OHAM solution is an improved version of homotopy-based method with the aim of obtaining an accurate approximation with just two-three terms of the series. For a quantitative assessment, the performances of different solutions are compared numerically in Tables 1, 2, 3.

TABLE 1

TABLE 1. Comparison between HAM-, HPM-, and OHAM-based approximations and numerical solution for the selected case for $t=0$ h: (a) $C_1\left(x,t\right)$ and (b) $C_2\left(x,t\right)$. (a) Species, $C_1\left(x,t\right)$ (b) Species, $C_2\left(x,t\right)$

TABLE 2

TABLE 2. Comparison between HAM-, HPM-, and OHAM-based approximations and numerical solution for the selected case for $t=20$ hr: (a) $C_1\left(x,t\right)$ and (b) $C_2\left(x,t\right)$. (a) Species, $C_1\left(x,t\right)$ (b) Species, $C_2\left(x,t\right)$

TABLE 3

TABLE 3. Comparison between HAM-, HPM-, and OHAM-based approximations and numerical solution for the selected case for $t=40$ h: (a) $C_1\left(x,t\right)$ and (b) $C_2\left(x,t\right)$. (a) Species, $C_1\left(x,t\right)$ (b) Species, $C_2\left(x,t\right)$

5 Concluding remarks

Here, we derive the HAM-, HPM-, and OHAM-based analytical solutions for two-species transport equations with spatially varying dispersion coefficient and seepage velocity. The same equations were studied analytically using HAM by Chaudhary and Singh (2020). A numerical solution is also developed using the MATLAB routine pdepe. The proposed methods produce accurate solutions for all the cases. This work shows the potential of HAM, HPM, and OHAM in the context of obtaining a series solution for a system of parabolic PDEs. The theoretical as well as numerical convergence of the series solutions are provided.

Data availability statement

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.

Author contributions

MK: conceptualization, methodology, writing—original draft. VS: conceptualization, supervision, writing—review and editing, project administration.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Abbasbandy, S., Shivanian, E., and Vajravelu, K. (2011). Mathematical properties of h-curve in the framework of the homotopy analysis method. Commun. Nonlinear Sci. Numer. Simul. 16 (11), 4268–4275. doi:10.1016/j.cnsns.2011.03.031

CrossRef Full Text | Google Scholar

Arnold, J., Duddu, R., Brown, K., and Kosson, D. S. (2017). Influence of multi-species solute transport on modeling of hydrated Portland cement leaching in strong nitrate solutions. Cem. Concr. Res. 100, 227–244. doi:10.1016/j.cemconres.2017.06.002

CrossRef Full Text | Google Scholar

Bai, B., Jiang, S., Liu, L., Li, X., and Wu, H. (2021b). The transport of silica powders and lead ions under unsteady flow and variable injection concentrations. Powder Technol. 387, 22–30. doi:10.1016/j.powtec.2021.04.014

CrossRef Full Text | Google Scholar

Bai, B., Nie, Q., Zhang, Y., Wang, X., and Hu, W. (2021a). Cotransport of heavy metals and SiO2 particles at different temperatures by seepage. J. Hydrology 597, 125771. doi:10.1016/j.jhydrol.2020.125771

CrossRef Full Text | Google Scholar

Batu, V. (2005). Applied flow and solute transport modeling in aquifers: Fundamental principles and analytical and numerical methods. Florida, United States: CRC Press.

Google Scholar

Bear, J. (1988). Dynamics of fluids in porous media. Massachusetts: Courier Corporation.

Google Scholar

Chamkha, A. J. (2005). Modeling of multi-species contaminant transport with spatially-dependent dispersion and coupled linear/non-linear reactions. Int. J. Fluid Mech. Res. 32 (1), 1–20. doi:10.1615/interjfluidmechres.v32.i1.10

CrossRef Full Text | Google Scholar

Chaudhary, M., and Singh, M. K. (2020). Study of multispecies convection-dispersion transport equation with variable parameters. J. Hydrology 591, 125562. doi:10.1016/j.jhydrol.2020.125562

CrossRef Full Text | Google Scholar

Chen, J. S., Lai, K. H., Liu, C. W., and Ni, C. F. (2012). A novel method for analytically solving multi-species advective–dispersive transport equations sequentially coupled with first-order decay reactions. J. Hydrology 420, 191–204. doi:10.1016/j.jhydrol.2011.12.001

CrossRef Full Text | Google Scholar

Clement, T. P., Sun, Y., Hooker, B. S., and Petersen, J. N. (1998). Modeling multispecies reactive transport in ground water. Groundw. Monit. Remediat. 18 (2), 79–92. doi:10.1111/j.1745-6592.1998.tb00618.x

CrossRef Full Text | Google Scholar

Domenico, P. A. (1987). An analytical model for multidimensional transport of a decaying contaminant species. J. Hydrology 91 (1-2), 49–58. doi:10.1016/0022-1694(87)90127-2

CrossRef Full Text | Google Scholar

Fujikawa, Y., and Fukui, M. (1990). Adsorptive solute transport in fractured rock: Analytical solutions for delta-type source conditions. J. Contam. Hydrology 6 (1), 85–102. doi:10.1016/0169-7722(90)90013-7

CrossRef Full Text | Google Scholar

He, J. H. (1999). Homotopy perturbation technique. Comput. Methods Appl. Mech. Eng. 178 (3-4), 257–262. doi:10.1016/s0045-7825(99)00018-3

CrossRef Full Text | Google Scholar

Liao, S. (2003). Beyond perturbation: Introduction to the homotopy analysis method. Boca Raton: Chapman and Hall/CRC.

Google Scholar

Liao, S. (2012). Homotopy analysis method in nonlinear differential equations. Beijing: Higher education press, 153–165.

Google Scholar

Liao, S. J. (1992). The proposed homotopy analysis technique for the solution of nonlinear problems. Doctoral dissertation, PhD thesis. Shanghai: Shanghai Jiao Tong University.

Google Scholar

Lunn, M., Lunn, R. J., and Mackayb, R. (1996). Determining analytic solutions of multiple species contaminant transport, with sorption and decay. J. Hydrology 180 (1-4), 195–210. doi:10.1016/0022-1694(95)02891-9

CrossRef Full Text | Google Scholar

Marinca, V., and Herişanu, N. (2008). Application of optimal homotopy asymptotic method for solving nonlinear equations arising in heat transfer. Int. Commun. Heat Mass Transf. 35 (6), 710–715. doi:10.1016/j.icheatmasstransfer.2008.02.010

CrossRef Full Text | Google Scholar

Natarajan, N., and Kumar, G. S. (2010). Finite difference approach for modeling multispecies transport in porous media. Int. J. Eng. Sci. Technol. 2, 3344–3350.

Google Scholar

Natarajan, N., and Kumar, G. S. (2018). Spatial moment analysis of multispecies contaminant transport in porous media. Environ. Eng. Res. 23 (1), 76–83. doi:10.4491/eer.2016.147

CrossRef Full Text | Google Scholar

Simpson, M. J., and Ellery, A. J. (2014). Exact series solutions of reactive transport models with general initial conditions. J. Hydrology 513, 7–12. doi:10.1016/j.jhydrol.2014.03.035

CrossRef Full Text | Google Scholar

Skeel, R. D., and Berzins, M. (1990). A method for the spatial discretization of parabolic equations in one space variable. SIAM J. Sci. Stat. Comput. 11 (1), 1–32. doi:10.1137/0911001

CrossRef Full Text | Google Scholar

Slodička, M., and Balážová, A. (2010). Decomposition method for solving multi-species reactive transport problems coupled with first-order kinetics applicable to a chain with identical reaction rates. J. Comput. Appl. Math. 234 (4), 1069–1077. doi:10.1016/j.cam.2009.04.021

CrossRef Full Text | Google Scholar

Slodička, M., and Balážová, A. (2008). Singular value decomposition method for multi-species first-order reactive transport with identical decay rates. Transp. Porous Media 73 (2), 161–172. doi:10.1007/s11242-007-9175-7

CrossRef Full Text | Google Scholar

Sposito, G., Gupta, V. K., and Bhattacharya, R. N. (1979). Foundation theories of solute transport in porous media: A critical review. Adv. Water Resour. 2, 59–68. doi:10.1016/0309-1708(79)90012-5

CrossRef Full Text | Google Scholar

Sun, Y., and Clement, T. P. (1999). A decomposition method for solving coupled multi–species reactive transport problems. Transp. Porous Media 37 (3), 327–346. doi:10.1023/a:1006507514019

CrossRef Full Text | Google Scholar

Sun, Y., Petersen, J. N., and Clement, T. P. (1999). Analytical solutions for multiple species reactive transport in multiple dimensions. J. Contam. Hydrology 35 (4), 429–440. doi:10.1016/s0169-7722(98)00105-3

CrossRef Full Text | Google Scholar

Vajravelu, K., and Van Gorder, R. (2013). Nonlinear flow phenomena and homotopy analysis. Berlin: Springer.

Google Scholar

Van Genuchten, M. T. (1985). Convective-dispersive transport of solutes involved in sequential first-order decay reactions. Comput. Geosciences 11 (2), 129–147. doi:10.1016/0098-3004(85)90003-2

CrossRef Full Text | Google Scholar