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

Abstract

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.

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 between two functions ; , where is a dimension (space or time), is itself a continuous function, defined as and satisfies when and when , where and are the topological spaces. This shows that as goes from 0 to 1, varies from to . The functions and 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:where are the non-linear operators or the original operators of the system of equations; are the unknown variables for say ; and are the independent variables. Now, the zeroth-order deformation equation can be constructed as follows (Liao 2003):

subject to the initial and boundary conditions:where is the embedding parameter, are the representations of solutions across , are the initial approximations, are the auxiliary parameters, are the auxiliary functions, and and are the linear and non-linear operators, respectively. Eq. 2 shows that as , , and as , . This means as varies from 0 to 1, 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):whereandwhere for 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:

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:

Now, we can construct a homotopy that satisfies (He 1999):where the symbols denote the same variables as mentioned in the previous section. Also, similar to HAM, Eq. 9 shows that at , , and at , . Let us now express as a series in terms of as follows:where for are the higher-order terms. As , Eq. 10 produces the final solutions as:

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 to obtain the terms of the series . 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 .

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:

subject to the initial and boundary conditions:where symbols denote the same variables as discussed in the previous section. Following HAM, one can construct a family of equations as follows:where symbols have their usual meaning and are the unknown parameters. The auxiliary functions are defined as:

It is easy to verify that when , and when , , which is the same as HAM and HPM, i.e., as goes from 0 to 1, we have the continuous deformation from the initial approximation to the final solution. Now, the initial approximations can be found by solving the following set of equations, which is obtained after substituting in Eq. 14:

subject to the boundary conditions:

In OHAM, expanding the auxiliary function with respect to the embedding parameter is the key step. The auxiliary function can be expressed as follows:where ’s are the auxiliary functions that depend on the independent variables and and parameters . 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:

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

the initial and boundary conditions:

For ,

subject to the initial and boundary conditions:where the term is the coefficient of , which is obtained by expanding as follows:

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 , and there exist many ways to choose it. According to Marinca and Herişanu (2015), one can choose in such a way that the product of Eq. 22 is of the same form as that of . 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 , then

Finally, the approximate solution can be obtained by considering a finite number of terms of the series Eq. 25. The unknown parameters 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):Here, and represents the number of species; and denote the spatial and temporal coordinates, respectively; is the retardation factor for the -th species; is the concentration strength for the -th species; is the decay rate coefficient for the -th species; and and are the dispersion coefficient and seepage velocity, respectively.

The present work considers two-species system (i.e., ) and includes the effect of spatially and temporally dependent transport parameters on pollutant migration. To that end, the model becomes (Chaudhary and Singh 2020):where and represent the pollutant concentration level for parent and daughter species, respectively.

The initial and boundary conditions for the model can be given as follows:where and are considered as (Simpson and Ellery 2014):

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 and , 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:where and are the initial dispersion and velocity, respectively; and are the parameters, where is knows as the heterogeneity parameter. Eq. 32 shows that dispersion is directly proportional to the square of velocity (Batu 2005). For and , i.e., is very small, the system converts into a homogeneous one, where the effect of space on dispersion and seepage velocity is absent. As 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:

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:

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 can be calculated from Eq. 6 as follows:

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 of the present problem:

so thatwhere 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:where 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:where are given by Eqs 37, 38, and constant can be determined from the initial condition for the higher-order deformation equations, which simply yields for all .

Now, the auxiliary functions can be determined from the rule of coefficient ergodicity. Based on the rule of solution expression, the general form of should be:where , , and are the integers. However, for simplicity, we can take (Vajravelu and Van Gorder 2013). Finally, the approximate solutions can be obtained as follows:

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

FIGURE 1

3.2 HPM-based solution

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

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 :

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

The initial approximation can be chosen as and . Using this initial approximation, we can solve the equations iteratively using a symbolic software. Finally, the HPM-based solutions can be approximated as follows:

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

FIGURE 2

3.3 OHAM-based solution

In relation to the discussion given in Section 2.3, we choose the following operators:andWith these considerations, we can solve the zeroth-order Eq. 16 to have the following solution:

Now to have the higher-order equations, we need the expressions , , , 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:

The auxiliary functions can be chosen in many ways. Here, we select and Putting in Eq. 22 and then using the expansion of , we have the following second-order equation:

Here, we choose and . 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 . Finally, the approximate solution can be found as: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

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 and seepage velocity are taken as given by Eq. 2. Using these equations, we have:

For the validation of solutions, space and time domains are considered as and . The required parameters are chosen as (Simpson and Ellery 2014; Chaudhary and Singh 2020): , , , , , and . Also, the other parameters are considered as , , , and

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 . 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 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 (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 5^th and 10th order HAM-based approximations for the value . 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).

FIGURE 4

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 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.

FIGURE 5

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

4.4 Validation of OHAM-based solution

For the assessment of the OHAM-based analytical solution, one needs to calculate the constants ’s. For that purpose, we calculate the residual as follows:where is the approximate solution. When , becomes the exact solution to the problem. One of the ways to obtain the optimal ’s, for which the solution converges, is the minimization of squared residual error, i.e.,where is the domain of the problem. The minimization of Eq. 68 leads to a system of algebraic equations as follows:

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

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.

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.

References

• 1

  AbbasbandyS.ShivanianE.VajraveluK. (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

  • CrossRef
  • Google Scholar

• 2

  ArnoldJ.DudduR.BrownK.KossonD. 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

  • CrossRef
  • Google Scholar

• 3

  BaiB.JiangS.LiuL.LiX.WuH. (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

  • CrossRef
  • Google Scholar

• 4

  BaiB.NieQ.ZhangY.WangX.HuW. (2021a). Cotransport of heavy metals and SiO2 particles at different temperatures by seepage. J. Hydrology597, 125771. 10.1016/j.jhydrol.2020.125771

  • CrossRef
  • Google Scholar

• 5

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

  • Google Scholar

• 6

  BearJ. (1988). Dynamics of fluids in porous media. Massachusetts: Courier Corporation.

  • Google Scholar

• 7

  ChamkhaA. 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

  • CrossRef
  • Google Scholar

• 8

  ChaudharyM.SinghM. K. (2020). Study of multispecies convection-dispersion transport equation with variable parameters. J. Hydrology591, 125562. 10.1016/j.jhydrol.2020.125562

  • CrossRef
  • Google Scholar

• 9

  ChenJ. S.LaiK. H.LiuC. W.NiC. F. (2012). A novel method for analytically solving multi-species advective–dispersive transport equations sequentially coupled with first-order decay reactions. J. Hydrology420, 191–204. 10.1016/j.jhydrol.2011.12.001

  • CrossRef
  • Google Scholar

• 10

  ClementT. P.SunY.HookerB. S.PetersenJ. N. (1998). Modeling multispecies reactive transport in ground water. Groundw. Monit. Remediat.18 (2), 79–92. 10.1111/j.1745-6592.1998.tb00618.x

  • CrossRef
  • Google Scholar

• 11

  DomenicoP. A. (1987). An analytical model for multidimensional transport of a decaying contaminant species. J. Hydrology91 (1-2), 49–58. 10.1016/0022-1694(87)90127-2

  • CrossRef
  • Google Scholar

• 12

  FujikawaY.FukuiM. (1990). Adsorptive solute transport in fractured rock: Analytical solutions for delta-type source conditions. J. Contam. Hydrology6 (1), 85–102. 10.1016/0169-7722(90)90013-7

  • CrossRef
  • Google Scholar

• 13

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

  • CrossRef
  • Google Scholar

• 14

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

  • Google Scholar

• 15

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

  • Google Scholar

• 16

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

  • Google Scholar

• 17

  LunnM.LunnR. J.MackaybR. (1996). Determining analytic solutions of multiple species contaminant transport, with sorption and decay. J. Hydrology180 (1-4), 195–210. 10.1016/0022-1694(95)02891-9

  • CrossRef
  • Google Scholar

• 18

  MarincaV.HerişanuN. (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

  • CrossRef
  • Google Scholar

• 19

  NatarajanN.KumarG. S. (2010). Finite difference approach for modeling multispecies transport in porous media. Int. J. Eng. Sci. Technol.2, 3344–3350.

  • Google Scholar

• 20

  NatarajanN.KumarG. S. (2018). Spatial moment analysis of multispecies contaminant transport in porous media. Environ. Eng. Res.23 (1), 76–83. 10.4491/eer.2016.147

  • CrossRef
  • Google Scholar

• 21

  SimpsonM. J.ElleryA. J. (2014). Exact series solutions of reactive transport models with general initial conditions. J. Hydrology513, 7–12. 10.1016/j.jhydrol.2014.03.035

  • CrossRef
  • Google Scholar

• 22

  SkeelR. D.BerzinsM. (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

  • CrossRef
  • Google Scholar

• 23

  SlodičkaM.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

  • CrossRef
  • Google Scholar

• 24

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

  • CrossRef
  • Google Scholar

• 25

  SpositoG.GuptaV. K.BhattacharyaR. 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

  • CrossRef
  • Google Scholar

• 26

  SunY.ClementT. P. (1999). A decomposition method for solving coupled multi–species reactive transport problems. Transp. Porous Media37 (3), 327–346. 10.1023/a:1006507514019

  • CrossRef
  • Google Scholar

• 27

  SunY.PetersenJ. N.ClementT. P. (1999). Analytical solutions for multiple species reactive transport in multiple dimensions. J. Contam. Hydrology35 (4), 429–440. 10.1016/s0169-7722(98)00105-3

  • CrossRef
  • Google Scholar

• 28

  VajraveluK.Van GorderR. (2013). Nonlinear flow phenomena and homotopy analysis. Berlin: Springer.

  • Google Scholar

• 29

  Van GenuchtenM. T. (1985). Convective-dispersive transport of solutes involved in sequential first-order decay reactions. Comput. Geosciences11 (2), 129–147. 10.1016/0098-3004(85)90003-2

  • CrossRef
  • Google Scholar