Front. Earth Sci. Frontiers in Earth Science Front. Earth Sci. 2296-6463 Frontiers Media S.A. 1064110 10.3389/feart.2022.1064110 Earth Science Original Research Approximate analytical solutions for multispecies convection-dispersion transport equation with variable parameters Kumbhakar and Singh 10.3389/feart.2022.1064110 Kumbhakar Manotosh * † Singh Vijay P. Department of Biological and Agricultural Engineering, Zachry Department of Civil and Environmental Engineering, Texas A and M University, College Station, TX, United States

Edited by: Yingfang Zhou, University of Aberdeen, United Kingdom

Reviewed by: Qingzhi Hou, Tianjin University, China

Bing Bai, Beijing Jiaotong University, China

 *Correspondence: Manotosh Kumbhakar, manotosh.kumbhakar@gmail.com

Present address: Department of Civil Engineering, National Taiwan University, Taipei, Taiwan

This article was submitted to Hydrosphere, a section of the journal Frontiers in Earth Science

 01 02 2023 2022 10 1064110 07 10 2022 28 12 2022 Copyright © 2023 Kumbhakar and Singh. 2023 Kumbhakar and Singh

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

Multispecies pollutant migration often occurs in polluted groundwater systems. Most of the multispecies problems that have been dealt in the literature assume constant transport parameters, primarily because analytical solutions for varying parameters become a challenge. The present study analytically solves a two-species convection-dispersion transport equation, considering spatially varying dispersion coefficient and seepage velocity, which corresponds to the steady migration in a steady flow domain. Indeed, the methodology can be adopted for other cases, such as the transient migration in a steady flow domain and transient migration in an unsteady flow domain, without any difficulty. Three kinds of homotopy-based methods, namely the homotopy perturbation method (HPM), homotopy analysis method (HAM), and optimal homotopy asymptotic method (OHAM), are used to derive approximate analytical solutions in the form of truncated series. In homotopy analysis method, the convergence-control parameter ℏ plays a key role in the convergence of the series solution. It is observed that for a specific case of this parameter, namely ℏ=−1, the HAM-based solution recovers the HPM-based solution. For the verification of homotopy-based solutions, we utilize the MATLAB routine pdepe, which efficiently solves a class of parabolic PDEs (single/system). An excellent agreement is found between the derived analytical solutions and the numerical solutions for all three methods. Further, a quantitative assessment is carried out for the derived solutions. Also, convergence theorems are proposed for the series solutions obtained using HAM and OHAM.

 pollutant transport multispecies homotopy analytical solution series solution 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 H t ; q between two functions f t ; g t , where t is a dimension (space or time), is itself a continuous function, defined as H : T × 0,1 → U and satisfies H t ; q = f t when q = 0 and H t ; q = g t when q = 1 , where T and U are the topological spaces. This shows that as q goes from 0 to 1, H t ; q varies from f t to g t . The functions f t and g t 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: N i y i x , t = 0   where 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 = 1,2 , … , n ; x and t are the independent variables. Now, the zeroth-order deformation equation can be constructed as follows (Liao 2003): 1 − q L i Φ i x , t ; q − y i , 0 x , t = q ħ i H i x , t N i Φ i x , t ; q ,   i = 1,2 , … , n  

subject to the initial and boundary conditions: B Φ i , ∂ Φ i ∂ x , ∂ Φ i ∂ t = 0 ,   Ι Φ i , ∂ Φ i ∂ x , ∂ Φ i ∂ t = 0   where q is the embedding parameter, Φ i x , t ; q are the representations of solutions across q , y i , 0 x , t are the initial approximations, ħ i are the auxiliary parameters, H i x , t are the auxiliary functions, and L i and N i are the linear and non-linear operators, respectively. Eq. 2 shows that as q = 0 , Φ i x , t ; 0 = y i , 0 x , t , and as q = 1 , Φ i x , t ; 1 = y i x , t . This means as q varies from 0 to 1, Φ i x , t ; q 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): L i y i , m x , t − χ m y i , m − 1 x , t = ħ i H i x , t R i , m y → i , m − 1 ,   m = 1,2,3 , . . .   where χ m =   0   w h e n   m = 1 , 1   o t h e r w i s e and R i , m y → i , m − 1 = 1 m − 1 ! ∂ m − 1 N i Φ i x , t ; q ∂ q m − 1 q = 0 where y i , m for m ≥ 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 x , t = y i , 0 x , t + ∑ m = 1 ∞ y i , m x , t  

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: N i y i x , t = f x , t  

Now, we can construct a homotopy that satisfies (He 1999): 1 − q L i Φ i x , t ; q − L i y i , 0 x , t + q N i Φ i x , t ; q − f x , t = 0 ,   i = 1,2 , … , n where the symbols denote the same variables as mentioned in the previous section. Also, similar to HAM, Eq. 9 shows that at q = 0 , Φ i x , t ; 0 = y i , 0 x , t , and at q = 1 , Φ i x , t ; 1 = y i x , t . Let us now express Φ i x , t ; q as a series in terms of q as follows: Φ i x , t ; q = Φ i , 0 + q Φ i , 1 + q 2 Φ i , 2 + q 3 Φ i , 3 + q 4 Φ i , 4 + . . . where Φ i , m for m ≥ 1 are the higher-order terms. As q → 1 , Eq. 10 produces the final solutions as: y i x , t = lim q → 1 Φ i x , t ; q = ∑ k = 0 ∞ Φ i , k

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 Φ 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: L i y i x , t + N i y i x , t + h i x , t = 0 ,   i = 1,2 , … , n  

subject to the initial and boundary conditions: B y i , ∂ y i ∂ x , ∂ y i ∂ t = 0 ,   Ι y i , ∂ y i ∂ x , ∂ y i ∂ t = 0 where symbols denote the same variables as discussed in the previous section. Following HAM, one can construct a family of equations as follows: 1 − q L i Φ i x , t , C i , j ; q + h i x , t = H i x , t , C i , j ; q L i Φ i x , t , C i , j ; q + h i x , t + N i Φ i x , t , C i , j ; q where symbols have their usual meaning and C i , j are the unknown parameters. The auxiliary functions are defined as: H i x , t , C i , j ; q = 0   f o r   q = 0 ≠ 0   f o r   q ∈ ( 0,1

It is easy to verify that when q = 0 , Φ i x , t , C i , j ; q = y i , 0 x , t and when q = 1 , Φ i x , t , C i , j ; q = y i x , t , 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 x , t can be found by solving the following set of equations, which is obtained after substituting q = 0 in Eq. 14: L i y i , 0 x , t + h i x , t = 0  

subject to the boundary conditions: B y i , 0 , ∂ y i , 0 ∂ x , ∂ y i , 0 ∂ t = 0 ,   Ι y i , 0 , ∂ y i , 0 ∂ x , ∂ y i , 0 ∂ t = 0

In OHAM, expanding the auxiliary function with respect to the embedding parameter is the key step. The auxiliary function H i x , t , C i , j ; q can be expressed as follows: H i x , t , C i , j ; q = q H i , 1 x , t , C i , j + q 2 H i , 2 x , t , C i , j + q 3 H i , 3 x , t , C i , j + . . . where H i x , t , C i , j ’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: Φ i x , t , C i , j ; q = y i , 0 x , t + ∑ j = 1 ∞ y i , j x , t , C i , j q j

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: L i y i , 1 x , t , C i , j   = H i , 1 x , t , C i , j N i , 0 y i , 0 x , t

s u b j e c t   t o the initial and boundary conditions: B y i , 1 , ∂ y i , 1 ∂ x , ∂ y i , 1 ∂ t = 0 ,   Ι y i , 1 , ∂ y i , 1 ∂ x , ∂ y i , 1 ∂ t = 0

For k = 2,3,4 , . . . , L i y i , k x , t , C i , j − y i , k − 1 x , t , C i , j = H i , k x , t , C i , j N i , 0 y i , 0 x , t + ∑ j = 1 k − 1 H i , j x , t , C i , j [ L i y i , k − j x , t , C i , j + N i , k − j [ y i , 0 x , t , y i , 1 x , t , C i , j , … , y i , k − j x , t , C i , j ] ]

subject to the initial and boundary conditions: B y i , k , ∂ y i , k ∂ x , ∂ y i , k ∂ t = 0 ,   Ι y i , k , ∂ y i , k ∂ x , ∂ y i , k ∂ t = 0   where the term N i , k − j y i , 0 x , t , y i , 1 x , t , C i , j , … , y i , k − j x , t , C i , j is the coefficient of q m , which is obtained by expanding N i Φ i x , t , C i , j ; q as follows: N i Φ i x , t , C i , j ; q = N i , 0 y i , 0 x , t + q N i , 1 y i , 0 x , t , y i , 1 x , t , C i , j + q 2 N i , 2 y i , 0 x , t , y i , 1 x , t , C i , j , y i , 2 x , t , C i , j + . . .

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 x , t , C i , j , and there exist many ways to choose it. According to Marinca and Herişanu (2015), one can choose H i , j x , t , C i , j in such a way that the product H i , j x , t , C i , j L i y i , k − j x , t , C i , j + N i , k − j y i , 0 x , t , y i , 1 x , t , C i , j , … , y i , k − j x , t , C i , j of Eq. 22 is of the same form as that of H i , j x , t , C i , j . 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 x , t , C i , j = y i , 0 x , t + ∑ j = 1 ∞ y i , j x , t , C i , j

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 ∂ C i ∂ t = ∂ ∂ x D x , t ∂ C i ∂ x − v x , t C i + k i − 1 C i − 1 − k i C i ,   i = 1 ,   2 , . . . , n 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 x , t and v x , t 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 ∂ C 1 ∂ t = ∂ ∂ x D x , t ∂ C 1 ∂ x − v x , t C 1 − k 1 C 1 r 2 ∂ C 2 ∂ t = ∂ ∂ x D x , t ∂ C 2 ∂ x − v x , t C 2 − k 2 C 2 + k 1 C 1 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 x , 0 = f 1 x ,   C 1 0 , t = f 1 0 ,   lim x → ∞ C 1 x , t = 0   C 2 x , 0 = f 2 x ,   C 2 0 , t = f 2 0 ,   lim x → ∞ C 2 x , t = 0 where f 1 x and f 2 x are considered as (Simpson and Ellery 2014): f 1 x = a 1 x e x p − a 2 x ,   f 2 x = 0

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 x , t and v x , t , 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 x , t = D 0 α 1 + α 2 x 2 ,   v x , t = v 0 α 1 + α 2 x where D 0 and v 0 are the initial dispersion and velocity, respectively; α 1 and α 2 are the parameters, where α 2 is knows as the heterogeneity parameter. Eq. 32 shows that dispersion is directly proportional to the square of velocity (Batu 2005). For α 1 = 1 and α 2 → 0 , i.e., α 2 is very small, the system converts into a homogeneous one, where the effect of space on dispersion and seepage velocity is absent. As α 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 ∂ C 1 ∂ t = D ∂ 2 C 1 ∂ x 2 + ∂ D ∂ x − v ∂ C 1 ∂ x − k 1 + ∂ v ∂ x C 1   r 2 ∂ C 2 ∂ t = D ∂ 2 C 2 ∂ x 2 + ∂ D ∂ x − v ∂ C 2 ∂ x − k 2 + ∂ v ∂ x C 2 + k 1 C 1

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: N 1 Φ i x , t ; q = r 1 ∂ Φ 1 x , t ; q ∂ t − D ∂ 2 Φ 1 x , t ; q ∂ x 2 − ∂ D ∂ x − v ∂ Φ 1 x , t ; q ∂ x + k 1 + ∂ v ∂ x Φ 1 x , t ; q N 2 Φ i x , t ; q = r 2 ∂ Φ 2 x , t ; q ∂ t − D ∂ 2 Φ 2 x , t ; q ∂ x 2 − ∂ D ∂ x − v ∂ Φ 2 x , t ; q ∂ x + k 2 + ∂ v ∂ x Φ 2 x , t ; q − k 1 Φ 1 x , t ; q

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 C → i , m − 1 = r 1 ∂ C 1 , m − 1 ∂ t − D ∂ 2 C 1 , m − 1 ∂ x 2 − ∂ D ∂ x − v ∂ C 1 , m − 1 ∂ x + k 1 + ∂ v ∂ x C 1 , m − 1 R 2 , m C → i , m − 1 = r 2 ∂ C 2 , m − 1 ∂ t − D ∂ 2 C 2 , m − 1 ∂ x 2 − ∂ D ∂ x − v ∂ C 2 , m − 1 ∂ x + k 2 + ∂ v ∂ x C 2 , m − 1 − k 1 C 1 , m − 1

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 x , t of the present problem: x m t p exp − n b x   |   m , n , p = 0,1,2 ,   3 , . . .  

so that C i x , t = ∑ p = 0 ∞ ∑ m = 0 ∞ ∑ n = 1 ∞ β m , n , p x m t p exp − n b x where β 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: L i Φ i x , t ; q = r i ∂ Φ i x , t ; q ∂ t   w i t h   t h e   p r o p e r t y   L i E 2 = 0 C 1,0 x , t = f 1 x = a 1 x ⁡ exp − a 2 x ,   C 2,0 x , t = f 2 x = 0 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 x , t = χ m C i , m − 1 x , t + ℏ i r i ∫ 0 t H i x , t R i , m C → i , m − 1 d t + E 2 ,   m = 1,2,3 , . . . 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 ≥ 1 .

Now, the auxiliary functions H i x , t can be determined from the rule of coefficient ergodicity. Based on the rule of solution expression, the general form of H i x , t should be: H i x , t = x n 1 t n 3 ⁡ exp − b n 2 x where n 1 , n 2 , and n 3 are the integers. However, for simplicity, we can take H i x , t = 1 (Vajravelu and Van Gorder 2013). Finally, the approximate solutions can be obtained as follows: C i x , t ≈ C i , 0 x , t + ∑ m = 1 M C i , m x , t

For the convenience of readers, a flowchart containing the steps of the method is provided in 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: L i Φ i x , t ; q = r i ∂ Φ i x , t ; q ∂ t   ,   i = 1,2

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 ∂ Φ 1,0 ∂ t − ∂ C 1,0 ∂ t = 0   s u b j e c t   t o   Φ 1,0 x , 0 = f 1 x ,   Φ 1,0 0 , t = f 1 0 ,   lim x → ∞ Φ 1,0 x , t = 0 r 2 ∂ Φ 2,0 ∂ t − ∂ C 2,0 ∂ t = 0   s u b j e c t   t o   Φ 2,0 x , 0 = f 2 x ,   Φ 2,0 0 , t = f 2 0 ,   lim x → ∞ Φ 2,0 x , t = 0 r 1 ∂ Φ 1,1 ∂ t − ∂ Φ 1,0 ∂ t − ∂ C 1,0 ∂ t + r 1 ∂ Φ 1,0 ∂ t − D ∂ 2 Φ 1,0 ∂ x 2 − ∂ D ∂ x − v ∂ Φ 1,0 ∂ x + k 1 + ∂ v ∂ x Φ 1,0 = 0   s u b j e c t   t o   Φ 1,1 x , 0 = 0 ,   Φ 1,1 0 , t = 0 ,   lim x → ∞ Φ 1,1 x , t = 0 r 2 ∂ Φ 2,1 ∂ t − ∂ Φ 2,0 ∂ t − ∂ C 2,0 ∂ t + r 2 ∂ Φ 2,0 ∂ t − D ∂ 2 Φ 2,0 ∂ x 2 − ∂ D ∂ x − v ∂ Φ 2,0 ∂ x + k 2 + ∂ v ∂ x Φ 2,0 − k 1 Φ 1,0 = 0   s u b j e c t   t o   Φ 2,1 x , 0 = 0 ,   Φ 2,1 0 , t = 0 ,   lim x → ∞ Φ 2,1 x , t = 0 r 1 ∂ Φ 1,2 ∂ t − ∂ Φ 1,1 ∂ t + r 1 ∂ Φ 1,1 ∂ t − D ∂ 2 Φ 1,1 ∂ x 2 − ∂ D ∂ x − v ∂ Φ 1,1 ∂ x + k 1 + ∂ v ∂ x Φ 1,1 = 0   s u b j e c t   t o   Φ 1,2 x , 0 = 0 ,   Φ 1,2 0 , t = 0 ,   lim x → ∞ Φ 1,2 x , t = 0 r 2 ∂ Φ 2,2 ∂ t − ∂ Φ 2,1 ∂ t + r 2 ∂ Φ 2,1 ∂ t − D ∂ 2 Φ 2,1 ∂ x 2 − ∂ D ∂ x − v ∂ Φ 2,1 ∂ x + k 2 + ∂ v ∂ x Φ 2,1 − k 1 Φ 1,1 = 0   s u b j e c t   t o   Φ 2,2 x , 0 = 0 ,   Φ 2,2 0 , t = 0 ,   lim x → ∞ Φ 2,2 x , t = 0

Proceeding in a like manner, one can arrive at the following recurrence relation: r 1 ∂ Φ 1 , m ∂ t − ∂ Φ 1 , m − 1 ∂ t + r 1 ∂ Φ 1 , m − 1 ∂ t − D ∂ 2 Φ 1 , m − 1 ∂ x 2 − ∂ D ∂ x − v ∂ Φ 1 , m − 1 ∂ x + k 1 + ∂ v ∂ x Φ 1 , m − 1 = 0   s u b j e c t   t o   Φ 1 , m x , 0 = 0 ,   Φ 1 , m 0 , t = 0 ,   lim x → ∞ Φ 1 , m x , t = 0   f o r   m ≥ 2 r 2 ∂ Φ 2 , m ∂ t − ∂ Φ 2 , m − 1 ∂ t + r 2 ∂ Φ 2 , m − 1 ∂ t − D ∂ 2 Φ 2 , m − 1 ∂ x 2 − ∂ D ∂ x − v ∂ Φ 2 , m − 1 ∂ x + k 2 + ∂ v ∂ x Φ 2 , m − 1 − k 1 Φ 1 , m − 1 = 0   s u b j e c t   t o   Φ 2 , m x , 0 = 0 ,   Φ 2 , m 0 , t = 0 ,   lim x → ∞ Φ 2 , m x , t = 0   f o r   m ≥ 2  

The initial approximation can be chosen as Φ 1,0 = f 1 x = a 1 x e x p − a 2 x and Φ 2,0 = f 2 x = 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 x , t ≈ ∑ k = 0 M Φ i , k

A flowchart containing the steps of the method is given in 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: L i C i x , t = r i ∂ C i x , t ∂ t   N 1 C i x , t = − D ∂ 2 C 1 x , t ∂ x 2 − ∂ D ∂ x − v ∂ C 1 x , t ∂ x + k 1 + ∂ v ∂ x C 1 x , t N 2 C i x , t = − D ∂ 2 C 2 x , t ∂ x 2 − ∂ D ∂ x − v ∂ C 2 x , t ∂ x + k 2 + ∂ v ∂ x C 2 x , t − k 1 C 1 x , t and h i x , t = 0   With these considerations, we can solve the zeroth-order Eq. 16 to have the following solution: C i , 0 x , t = f i x

Now to have the higher-order equations, we need the expressions N i , 0 , N i , 1 , 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 ∂ C 1,1 ∂ t = H 1,1 x , t , D 1 , i − D ∂ 2 C 1,0 ∂ x 2 − ∂ D ∂ x − v ∂ C 1,0 ∂ x + k 1 + ∂ v ∂ x C 1,0   s u b j e c t   t o   C 1,1 x , 0 = 0 ,   C 1,1 0 , t = 0 ,   lim x → ∞ C 1,1 x , t = 0 r 2 ∂ C 2,1 ∂ t = H 2,1 x , t , D 2 , i − D ∂ 2 C 2,0 ∂ x 2 − ∂ D ∂ x − v ∂ C 2,0 ∂ x + k 2 + ∂ v ∂ x C 2,0 − k 1 C 1,0   s u b j e c t   t o   C 2,1 x , 0 = 0 ,   C 2,1 0 , t = 0 ,   lim x → ∞ C 2,1 x , t = 0

The auxiliary functions can be chosen in many ways. Here, we select H 1,1 x , t , D 1 , i = D 1,1 and H 2,1 x , t , D 2 , i = D 2,1 . Putting k = 2 in Eq. 22 and then using the expansion of N , we have the following second-order equation: r 1 ∂ C 1,2 ∂ t = r 1 ∂ C 1,1 ∂ t + D 1,2 − D ∂ 2 C 1,0 ∂ x 2 − ∂ D ∂ x − v ∂ C 1,0 ∂ x + k 1 + ∂ v ∂ x C 1,0 + D 1,1 r 1 ∂ C 1,1 ∂ t − D ∂ 2 C 1,1 ∂ x 2 − ∂ D ∂ x − v ∂ C 1,1 ∂ x + k 1 + ∂ v ∂ x C 1,1   s u b j e c t   t o   C 1,2 x , 0 = 0 ,   C 1,2 0 , t = 0 ,   lim x → ∞ C 1,2 x , t = 0 r 2 ∂ C 2,2 ∂ t = r 2 ∂ C 2,1 ∂ t + D 2,2 − D ∂ 2 C 2,0 ∂ x 2 − ∂ D ∂ x − v ∂ C 2,0 ∂ x + k 2 + ∂ v ∂ x C 2,0 − k 1 C 1,0 + D 2,1 r 2 ∂ C 2,1 ∂ t − D ∂ 2 C 2,1 ∂ x 2 − ∂ D ∂ x − v ∂ C 2,1 ∂ x + k 2 + ∂ v ∂ x C 2,1 − k 1 C 1,1   s u b j e c t   t o   C 1,2 x , 0 = 0 ,   C 1,2 0 , t = 0 ,   lim x → ∞ C 1,2 x , t = 0  

Here, we choose H 1,2 x , t , D 2 , i = D 1,2 and H 2,2 x , t , D 2 , i = D 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 x , t ≈ C i , 0 x , t + C i , 1 x , t , D i , i + C 1,2 x , t , D i , i 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.

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 x , t and seepage velocity v 0 x , t are taken as given by Eq. 2. Using these equations, we have: ∂ D ∂ x = 2 α 2 D 0 α 1 + α 2 x ,   ∂ v ∂ x = α 2 v 0  

For the validation of solutions, space and time domains are considered as 0 ≤ x c m ≤ 3000 and 0 ≤ t h r ≤ 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   / h r , k 2 = 0.01   / h r , v 0 = 0.2   c m / h r , and D 0 = 0.5   c m 2 / h r . Also, the other parameters are considered as r 1 = 2 , r 2 = 2.5 , α 1 = 1 , and α 2 = 0.05   / c m .

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 ℏ − 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 ℏ 1 = ℏ 2 = ℏ . One can indeed choose two different auxiliary parameters; however, it is better to first consider a common ℏ to see the accuracy of solutions, as this assumption simplifies the problem. Here, we calculate the so-called ℏ − c u r v e s to find a suitable choice for the auxiliary parameter ℏ , 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 x , t (or its derivatives) at some point within the domain. The flatness of the ℏ − curve determines a suitable choice for the auxiliary parameter ℏ . In Figure 4, we plot the ℏ − curves for 5th and 10th order HAM-based approximations for the value C i 1,1 . 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 ℏ . Any choice of ℏ within this range determines an optimal value for which the series solutions converge (Abbasbandy et al., 2011).

ℏ − curves for C i 1,1 of the HAM solution: (A) 5th 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 ℏ 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.

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.

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 x , t , D i , j = L i C i , O H A M x , t , D i , j + N i C i , O H A M x , t , D i , j + h i x , t , j = 1,2 , … , s   where C i , O H A M x , t , D i , j is the approximate solution. When R i x , t , D i , j = 0 , C i , O H A M x , t , D i , j 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 D i , j = ∫ x , t ϵ D ∑ i = 1 2 R i 2 x , t , D i , j d x d t ,   j = 1,2 , . . . , s   where D is the domain of the problem. The minimization of Eq. 68 leads to a system of algebraic equations as follows: ∂ J ∂ D 1,1 = ∂ J ∂ D 1,2 = . . . = ∂ J ∂ D 1 , s = 0  

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.

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.

Comparison between HAM-, HPM-, and OHAM-based approximations and numerical solution for the selected case for t = 0  h: (a) C 1 x , t and (b) C 2 x , t . (a) Species, C 1 x , t (b) Species, C 2 x , t

x   c m t = 0  h
Numerical solution 3rd order HAM-based approximation 5th order HPM-based approximation Three-term OHAM-based approximation
0 0 0.0000 0.0000 0.0000
250 267.6307 267.6299 267.6299 267.6299
500 286.5048 286.5039 286.5039 286.5040
750 230.0325 230.0318 230.0318 230.0319
1000 164.1700 164.1696 164.1696 164.1696
1250 109.8423 109.8421 109.8421 109.8421
1500 70.5532 70.5531 70.5531 70.5531
1750 44.0585 44.0585 44.0585 44.0585
2000 26.9518 26.9518 26.9518 26.9518
2250 16.2295 16.2295 16.2295 16.2295
2500 9.6523 9.6523 9.6523 9.6523
2750 5.6831 5.6831 5.6831 5.6831
3000 3.3185 3.3185 3.3185 3.3185
x   c m t = 0  h
Numerical solution 3rd order HAM-based approximation 5th order HPM-based approximation Three-term OHAM-based approximation
0 0 0 0 0
250 0 0.0005 0.0005 0.0005
500 0 0.0006 0.0006 0.0005
750 0 0.0005 0.0005 0.0004
1000 0 0.0003 0.0003 0.0003
1250 0 0.0002 0.0002 0.0002
1500 0 0.0001 0.0001 0.0001
1750 0 0 0 0
2000 0 0 0 0
2250 0 0 0 0
2500 0 0 0 0
2750 0 0 0 0
3000 0 0 0 0

Comparison between HAM-, HPM-, and OHAM-based approximations and numerical solution for the selected case for t = 20 hr: (a) C 1 x , t and (b) C 2 x , t . (a) Species, C 1 x , t (b) Species, C 2 x , t

x   c m t = 20  h
Numerical solution 3rd order HAM-based approximation 5th order HPM-based approximation Three-term OHAM-based approximation
0 0 0 −1.6183 0
250 140.4908 140.6578 140.4567 141.3510
500 157.6763 157.7893 157.6617 158.9841
750 133.5315 133.5544 133.5323 134.9613
1000 101.2569 101.2147 101.2632 102.6167
1250 72.5128 72.4397 72.5181 73.7020
1500 50.2082 50.1308 50.2103 51.1863
1750 34.0330 33.9673 34.0324 34.7970
2000 22.7495 22.7023 22.7475 23.3180
2250 15.0663 15.0383 15.0643 15.4687
2500 9.9165 9.9046 9.9151 10.1857
2750 6.5007 6.5008 6.5003 6.6685
3000 4.2511 4.2586 4.2513 4.3454
x   c m t = 20  h
Numerical solution 3rd order HAM-based approximation 5th order HPM-based approximation Three-term OHAM-based approximation
0 0 0 −0.7254 0
250 70.9251 70.7487 70.9525 71.4407
500 79.2222 79.1049 79.2369 79.3371
750 66.7362 66.7144 66.7373 66.3508
1000 50.3063 50.3567 50.3017 49.5507
1250 35.7906 35.8784 35.7860 34.8341
1500 24.6059 24.7026 24.6035 23.5909
1750 16.5518 16.6392 16.5517 15.5772
2000 10.9746 11.0433 10.9759 10.0988
2250 7.2061 7.2536 7.2077 6.4560
2500 4.7005 4.7284 4.7018 4.0814
2750 3.0527 3.0646 3.0534 2.5566
3000 1.9770 1.9771 1.9771 1.5890

Comparison between HAM-, HPM-, and OHAM-based approximations and numerical solution for the selected case for t = 40  h: (a) C 1 x , t and (b) C 2 x , t . (a) Species, C 1 x , t (b) Species, C 2 x , t

x   c m t = 40  h
Numerical solution 3rd order HAM-based approximation 5th order HPM-based approximation Three-term OHAM-based approximation
0 0 0 −1.8523 0
250 73.5227 72.6647 72.2756 76.3031
500 86.1393 86.6023 85.6178 86.3587
750 76.4081 77.7250 76.5445 74.4875
1000 60.9565 62.3186 61.3086 58.3735
1250 46.1183 47.0389 46.3814 43.9499
1500 33.8653 34.2290 33.9435 32.5570
1750 24.4288 24.3471 24.3591 23.9755
2000 17.4318 17.0942 17.2960 17.6135
2250 12.3584 11.9382 12.2277 12.9082
2500 8.7296 8.3447 8.6433 9.4223
2750 6.1556 5.8654 6.1230 6.8376
3000 4.3387 4.1579 4.3498 4.9247
x   c m t = 40  h
Numerical solution 3rd order HAM-based approximation 5th order HPM-based approximation Three-term OHAM-based approximation
0 0 −1.7574 −1.8424 −2.4474
250 95.3208 96.0586 96.4530 92.9257
500 110.6372 110.1981 111.1681 110.9119
750 97.1614 95.9354 97.1046 99.6570
1000 76.6791 75.3736 76.3932 79.9126
1250 57.3454 56.4095 57.0978 60.1912
1500 41.5970 41.1632 41.4886 43.5274
1750 29.6238 29.6214 29.6432 30.5633
2000 20.8592 21.1365 20.9503 20.9755
2250 14.5865 14.9910 14.6929 14.1301
2500 10.1591 10.5758 10.2443 9.3708
2750 7.0611 7.4205 7.1103 6.1314
3000 4.9045 5.1760 4.9188 3.9647
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. 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. 10.1016/j.cnsns.2011.03.031 Arnold J. Duddu R. Brown K. 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. 10.1016/j.cemconres.2017.06.002 Bai B. Jiang S. Liu L. Li X. Wu H. (2021b). The transport of silica powders and lead ions under unsteady flow and variable injection concentrations. Powder Technol. 387, 22–30. 10.1016/j.powtec.2021.04.014 Bai B. Nie Q. Zhang Y. Wang X. Hu W. (2021a). Cotransport of heavy metals and SiO2 particles at different temperatures by seepage. J. Hydrology 597, 125771. 10.1016/j.jhydrol.2020.125771 Batu V. (2005). Applied flow and solute transport modeling in aquifers: Fundamental principles and analytical and numerical methods. Florida, United States: CRC Press. Bear J. (1988). Dynamics of fluids in porous media. Massachusetts: Courier Corporation. 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. 10.1615/interjfluidmechres.v32.i1.10 Chaudhary M. Singh M. K. (2020). Study of multispecies convection-dispersion transport equation with variable parameters. J. Hydrology 591, 125562. 10.1016/j.jhydrol.2020.125562 Chen J. S. Lai K. H. Liu C. W. 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. 10.1016/j.jhydrol.2011.12.001 Clement T. P. Sun Y. Hooker B. S. Petersen J. N. (1998). Modeling multispecies reactive transport in ground water. Groundw. Monit. Remediat. 18 (2), 79–92. 10.1111/j.1745-6592.1998.tb00618.x Domenico P. A. (1987). An analytical model for multidimensional transport of a decaying contaminant species. J. Hydrology 91 (1-2), 49–58. 10.1016/0022-1694(87)90127-2 Fujikawa Y. Fukui M. (1990). Adsorptive solute transport in fractured rock: Analytical solutions for delta-type source conditions. J. Contam. Hydrology 6 (1), 85–102. 10.1016/0169-7722(90)90013-7 He J. H. (1999). Homotopy perturbation technique. Comput. Methods Appl. Mech. Eng. 178 (3-4), 257–262. 10.1016/s0045-7825(99)00018-3 Liao S. (2003). Beyond perturbation: Introduction to the homotopy analysis method. Boca Raton: Chapman and Hall/CRC. Liao S. (2012). Homotopy analysis method in nonlinear differential equations. Beijing: Higher education press, 153–165. Liao S. J. (1992). The proposed homotopy analysis technique for the solution of nonlinear problems. Doctoral dissertation, PhD thesis. Shanghai: Shanghai Jiao Tong University. Lunn M. Lunn R. J. Mackayb R. (1996). Determining analytic solutions of multiple species contaminant transport, with sorption and decay. J. Hydrology 180 (1-4), 195–210. 10.1016/0022-1694(95)02891-9 Marinca V. 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. 10.1016/j.icheatmasstransfer.2008.02.010 Natarajan N. Kumar G. S. (2010). Finite difference approach for modeling multispecies transport in porous media. Int. J. Eng. Sci. Technol. 2, 3344–3350. Natarajan N. Kumar G. S. (2018). Spatial moment analysis of multispecies contaminant transport in porous media. Environ. Eng. Res. 23 (1), 76–83. 10.4491/eer.2016.147 Simpson M. J. Ellery A. J. (2014). Exact series solutions of reactive transport models with general initial conditions. J. Hydrology 513, 7–12. 10.1016/j.jhydrol.2014.03.035 Skeel R. D. 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. 10.1137/0911001 Slodička M. 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. 10.1016/j.cam.2009.04.021 Slodička M. 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. 10.1007/s11242-007-9175-7 Sposito G. Gupta V. K. Bhattacharya R. N. (1979). Foundation theories of solute transport in porous media: A critical review. Adv. Water Resour. 2, 59–68. 10.1016/0309-1708(79)90012-5 Sun Y. Clement T. P. (1999). A decomposition method for solving coupled multi–species reactive transport problems. Transp. Porous Media 37 (3), 327–346. 10.1023/a:1006507514019 Sun Y. Petersen J. N. Clement T. P. (1999). Analytical solutions for multiple species reactive transport in multiple dimensions. J. Contam. Hydrology 35 (4), 429–440. 10.1016/s0169-7722(98)00105-3 Vajravelu K. Van Gorder R. (2013). Nonlinear flow phenomena and homotopy analysis. Berlin: Springer. Van Genuchten M. T. (1985). Convective-dispersive transport of solutes involved in sequential first-order decay reactions. Comput. Geosciences 11 (2), 129–147. 10.1016/0098-3004(85)90003-2 Appendix A. Convergence theorems

The convergence of the HAM- and OHAM-based solutions are proved theoretically using the following theorems.

A.1 Convergene Theorem of HAM-Based Solution.

The convergence theorem for the HAM-based solutions given by Eq. 45 can be proved using the following theorems.

Theorem A.1: If the homotopy series ∑ m = 0 ∞ C i , m x , t ,   ∑ m = 0 ∞ ∂ C i , m x , t ∂ t ,   ∑ m = 0 ∞ ∂ C i , m x , t ∂ x , and ∑ m = 0 ∞ ∂ 2 C i , m x , t ∂ x 2 converge, then R i , m C → i , m − 1 given by Eqs. 37, 38 satisfies the relation ∑ m = 1 ∞ R i , m C → i , m − 1 = 0 .

Proof: The auxiliary linear operators was defined as follows: L i C i = r i ∂ C i ∂ t  

According to Eq. 4, we obtain: L i C i , 1 = ℏ R i , 1 C i , 0 →   L i C i , 2 − C i , 1 = ℏ R i , 2 C i , 1 → L C i , 3 − C i , 2 = ℏ R i , 3 C i , 2 → L i C i , m − C i , m − 1 = ℏ R i , m C i , m − 1 →

Adding all the above terms, we get: L i C i , m = ℏ ∑ k = 1 m R i , k C i , k − 1 → As the series ∑ m = 0 ∞ C i , m x , t ,   ∑ m = 0 ∞ ∂ C i , m x , t ∂ t ,   ∑ m = 0 ∞ ∂ C i , m x , t ∂ x , and ∑ m = 0 ∞ ∂ 2 C i , m x , t ∂ x 2 are convergent, lim m → ∞ C i , m x , t = 0 and lim m → ∞ C i , m ′ x , t = 0 . Now, recalling the above summand and taking the limit, the required result follows as ℏ ∑ k = 1 ∞ R i , k C i , k − 1 → = ℏ lim m → ∞ ∑ k = 1 m R i , k C i , k − 1 → = lim m → ∞ L C i , m = r i lim m → ∞ C i , m ′ x , t = 0

Theorem A.2: If ℏ is so properly chosen that the series ∑ m = 0 ∞ C i , m x , t ,   ∑ m = 0 ∞ ∂ C i , m x , t ∂ t ,   ∑ m = 0 ∞ ∂ C i , m x , t ∂ x , and ∑ m = 0 ∞ ∂ 2 C i , m x , t ∂ x 2 converge absolutely to C i x , t , ∂ C i x , t ∂ t , ∂ C i x , t ∂ x , and ∂ 2 C i x , t ∂ x 2 , respectively, then the homotopy series ∑ m = 0 ∞ C i , m x , t satisfies the original governing Eqs 33, 34.

Proof: Theorem A.1 shows that if ∑ m = 0 ∞ C i , m x , t ,   ∑ m = 0 ∞ ∂ C i , m x , t ∂ t ,   ∑ m = 0 ∞ ∂ C i , m x , t ∂ x , and ∑ m = 0 ∞ ∂ 2 C i , m x , t ∂ x 2 converge then ∑ m = 1 ∞ R i , m C → i , m − 1 = 0 .

Therefore, using the expressions Eqs 37, 38, and simplifying further lead to: r 1 ∑ m = 0 ∞ ∂ C 1 , m ∂ t − D ∑ m = 0 ∞ ∂ 2 C 1 , m − 1 ∂ x 2 − ∂ D ∂ x − v ∑ m = 0 ∞ ∂ C 1 , m − 1 ∂ x + k 1 + ∂ v ∂ x ∑ m = 0 ∞ C 1 , m − 1 = 0 r 2 ∑ m = 0 ∞ ∂ C 2 , m − 1 ∂ t − D ∑ m = 0 ∞ ∂ 2 C 2 , m − 1 ∂ x 2 − ∂ D ∂ x − v ∑ m = 0 ∞ ∂ C 2 , m − 1 ∂ x + k 2 + ∂ v ∂ x ∑ m = 0 ∞ C 2 , m − 1 − k 1 ∑ m = 0 ∞ C 1 , m − 1 = 0 which is basically the original governing Eqs 33, 34. Furthermore, subject to the initial and boundary conditions C i , 0 x , 0 = f i x ,   C i , 0 0 , t = f i 0 ,   lim x → ∞ C i , 0 x , t = 0 , and the conditions for the higher-order deformation equation C i , m x , 0 = 0 C i , m 0 , t = 0 lim x → ∞ C i , m x , t = 0 , for m ≥ 1 , we easily obtain ∑ m = 0 ∞ C i , m x , 0 = f i x ,   ∑ m = 0 ∞ C i , m 0 , t = f i 0 , and lim x → ∞ ∑ m = 0 ∞ C i , m x , t = 0 . Hence, the convergence result follows.

A.2 Convergence Theorem of OHAM-Based Solution.

Theorem A.3: If the series C i , 0 x , t + ∑ j = 1 ∞ C i , j x , t , D i , i ,   i = 1,2 , . . . , s converges, where C i , j x , t , D i , i are governed by Eqs 60, 61, 62, 63, 64, then Eq. 65 is a solution of the original Eqs 33, 34.

Proof: Based on the choice of the auxiliary function, suppose that the series Eq. 65 is convergent. Then, we get: lim j → ∞ C i , j x , t , D i , i = 0 ,   i = 1,2 , . . . , s

One can write: C i , j x , t , D i , i = C i , 0 x , t , D i , i + C i , 1 x , t , D i , i − C i , 0 x , t , D i , i + C i , 2 x , t , D i , i − C i , 1 x , t , D i , i + . . . + C i , j x , t , D i , i − C i , j − 1 x , t , D i , i = C i , 0 x , t , D i , i + ∑ k = 1 j C i , k x , t , D i , i − C i , k − 1 x , t , D i , i ,   i = 1,2 , . . . , s  

Using Eq. A11, one can obtain from Eq. A10: 0 = lim j → ∞ C i , j x , t , D i , i = C i , 0 x , t , D i , i + ∑ k = 1 ∞ C i , k x , t , D i , i − C i , k − 1 x , t , D i , i ,   i = 1,2 , . . . , s

Eq. A12 can be rearranged as: C i , 0 x , t , D i , i + h i x , t − h i x , t + C i , 1 x , t , D i , i − C i , 0 x , t , D i , i + ∑ k = 2 ∞ C i , k x , t , D i , i − C i , k − 1 x , t , D i , i ,   i = 1,2 , . . . , s

Using the property of the linear operator, i.e., L A 1 x , t + A 2 x , t = L A 1 x , t + L A 2 x , t and L 0 = 0 , we have: 0 = L i 0 = L i C i , 0 x , t , D i , i + h i x , t + L i C i , 1 x , t , D i , i − L i C i , 0 x , t , D i , i + h i x , t + ∑ k = 2 ∞ L i C i , k x , t , D i , i − L i C i , k − 1 x , t , D i , i = H i , 1 x , t , D i , i N i , 0 C i , 0 x , t , D i , i + ∑ k = 2 ∞ H i , k x , t , D i , i N i , 0 C i , 0 x , t , D i , i + ∑ l = 1 k − 1 H i , l x , t , D i , i L i C i , k − l x , t , D i , i + N i , k − l C i , 0 x , t , D i , i , C i , 1 x , t , D i , i , . . . , C i , k − l x , t , D i , i = ∑ k = 1 ∞ H i , k x , t , D i , i N i , 0 C i , 0 x , t , D i , i + ∑ k = 2 ∞ ∑ l = 1 k − 1 H i , l x , t , D i , i × L i C i , k − l x , t , D i , i + N i , k − l C i , 0 x , t , D i , i , C i , 1 x , t , D i , i , . . . , C i , k − l x , t , D i , i = H i x , t , D i , i N i , 0 C i , 0 x , t , D i , i + ∑ k = 2 ∞ ∑ l = 1 k − 1 H i , l x , t , D i , i × L i C i , k − l x , t , D i , i + N i , k − l C i , 0 x , t , D i , i , C i , 1 x , t , D i , i , . . . , C i , k − l x , t , D i , i = H i x , t , D i , i N i , 0 C i , 0 x , t , D i , i + ∑ k = 1 ∞ H i , k x , t , D i , i ∑ p = 1 ∞ L i C i , p x , t , D i , i + N i , p C i , 0 x , t , D i , i , C i , 1 x , t , D i , i , . . . , C i , p x , t , D i , i = H i x , t , D i , i N i , 0 C i , 0 x , t , D i , i + H i x , t , D i , i × L i ∑ p = 1 ∞ C i , p x , t , D i , i + ∑ p = 1 ∞ N i , p C i , 0 x , t , D i , i , C i , 1 x , t , D i , i , . . . , C i , p x , t , D i , i = H i x , t , D i , i N i , 0 C i , 0 x , t , D i , i + H i x , t , D i , i × L i C i x , t , D i , i − L i C i , 0 x , t , D i , i + N i C i x , t , D i , i − N i C i , 0 x , t , D i , i = H i x , t , D i , i N i , 0 C i , 0 x , t , D i , i + H i x , t , D i , i × L i C i x , t , D i , i − L i C i , 0 x , t , D i , i + h i x , t + h i x , t + N i C i x , t , D i , i − N i C i , 0 x , t , D i , i = H i x , t , D i , i L i C i x , t , D i , i + h i x , t + N i C i x , t , D i , i ,   i = 1,2 , . . . , s

Now, since H i x , t , D i , i ≠ 0 , from Eq. A14 we have L i C i x , t , D i , i + h i x , t + N i C i x , t , D i , i = 0 ,   i = 1,2 , . . . , s  

which shows that C i x , t , D i , i is the exact solution of Eqs 33, 34.