Analytical Approximations of Well Function by Solving the Governing Differential Equation Representing Unsteady Groundwater Flow in a Confined Aquifer

Abstract

A solution of the governing equation representing the drawdown in a horizontal confined aquifer, where groundwater flow is unsteady, was first provided by Theis and is famously known as the Theis solution. This solution was given in terms of an exponential integral, also called the well function, for which simple and reliable approximations are preferred due to their practical applications. To that end, several approximations are available in the literature, of which some are based on series approximations for the integral, and others are numerical approximations. This study employs three kinds of homotopy-based methods, namely the homotopy analysis method (HAM), homotopy perturbation method (HPM), and optimal homotopy asymptotic method (OHAM), for analytically solving the governing partial differential equation (PDE). For convenience, the PDE is first converted to a boundary value problem (BVP) using a similarity transformation. Comparing the series approximations obtained using these methods with the Theis solution, it is found that the 10th-order HAM, and just three terms of OHAM-based solutions, provide accurate approximations. On the other hand, the HPM-based solution is found to be accurate only within a small domain. Further, the proposed approximations are compared with several series and numerical approximations available in the literature using the percentage error. The proposed methodology provides accurate approximations of the well function by directly solving the governing differential equation in a general framework and thus can be adapted to other practical situations arising in groundwater flow.

1. Introduction

Pumping test analysis in groundwater flow is one of the key methods for determining aquifer parameters, such as transmissivity and storativity. Unsteady groundwater flow in confined aquifers is described by coupling the continuity equation and the Darcy flux law. The resulting equation is a partial differential equation (PDE) of diffusion type. Under simplified assumptions, approximate analytical solutions of the PDE have been derived. One common assumption is that the confined aquifer is of uniform thickness and is homogeneous and isotropic. Using Boltzmann transformation, Theis [1] transformed the partial differential equation to an ordinary differential equation and then derived a solution in terms of an exponential integral, which came to be known as the well function, and gave tabular values of this function for different values of the Boltzmann variable or the argument of the well function. Theis found the unsteady flow of groundwater to be analogous to the unsteady flow of heat in a homogeneous solid and derived the desired solution.

In general, the exponential integral can occur in many applications of transient groundwater flow, hydrological problems, mathematical physics, and applied mathematics ([1,2,3]). The direct explicit evaluation of this integral may not be analytically tractable. Thus, most of the evaluations are based on approximations using either series expansion or numerical behavior. In groundwater studies, this integral is commonly known as the well function. Expressing this function as a power series and considering only two terms of the series, an approximate solution has also been proposed [4]. However, this series approximation is valid only within a small domain. The asymptotic (divergent) series can also be proposed for large values of the argument [5]. Many studies that are valid within a restricted domain are available for finding the approximation, based on polynomial or rational approximation or series expansions ([6,7,8,9,10,11]). Swamee and Ojha [12] combined several approximations valid for a specific region of the argument to provide an approximation to the well function. Barry et al. [13] constructed an approximation using the interpolation between large and small asymptotes. In recent work, Vatankhah [14] proposed a simple and very accurate approximation for the well function by combining the approximations of small and large values of the argument. However, these studies are based on the computational aspects of the exponential integral, not the solution methodology of the governing equations.

Differential equations play one of the most important roles in modelling in science and engineering. Therefore, the methods for solving these equations analytically are key topics of research. Most of the mathematical tools available are valid only for a specific class of nonlinear equations. On the other hand, Liao [15] proposed a methodology, namely the homotopy analysis method (HAM), using the basic concept of homotopy from topology to solve nonlinear equations analytically in series form. It was shown that the method provided great freedom for solving nonlinear problems, and its applicability was not confined to a small class of problems [16,17]. Following the usefulness of the method, two other variants, namely, the homotopy perturbation method (HPM) and optimal homotopy asymptotic method (OHAM), have also been proposed [18,19]. One of the key advantages of these methods is that they do not depend on the presence of a small/large parameter in the governing equation or boundary conditions. While these methods have been applied in several disciplines of science and engineering, their application to water engineering problems is limited. Further, apart from the well-known series expressions, the approximation of the well function using the closed-form formula mentioned in the previous paragraph is based on empirical observation or the characteristics of the function. Thus, the change in flow configuration or incorporation of other flow factors, which changes the solution to the problem, cannot be explained by the proposed approaches. Following this, the objective of this work is to apply the homotopy-based methods to unsteady groundwater flow in a confined aquifer and obtain an approximation for the well function by solving the governing differential equation itself. For that purpose, first, the methods are briefly described in a general framework and then applied to the considered problem. The convergence of the series solution is also discussed.

2. Brief Overview of Homotopy-Based Methods

For the convenience of the reader, here, we describe a general framework for three variants of homotopy-based methods. It may be noted that all of these methods are based on the fundamental concept of homotopy from topology, which describes the continuous deformations between two mathematical objects. Specifically, two objects can be called homotopic if one can be continuously deformed into the other. The circle and square in 2D and the doughnut and coffee cup in 3D are two examples of homotopy. In the context of differential equations, this idea was extended by [15], as they represent curves (or functions) in a mathematical sense. Liao [15] developed the homotopy analysis method (HAM), which was further extended by deriving two variants of this method, namely the homotopy perturbation method (HPM) and the optimal homotopy asymptotic method (OHAM). These three methods are employed in this paper.

2.1. Homotopy Analysis Method

Let us write a general differential equation in the following form:

$$\mathcal{N}\left(u\left(x,t\right)\right)=0$$

(1)

where $\mathcal{N}$ is the nonlinear operator or the original operator of the equation; $u$ is the dependent (unknown) variable; and $x$ and $t$ are the independent variables, e.g., space and time. Now, the zeroth-order deformation equation is constructed as follows [16]:

$$\left(1-q\right)ℒ\left[\Phi \left(x,t;q\right)-u_0\left(x,t\right)\right]=qћH\left(x,t\right)\mathcal{N}\left[\Phi \left(x,t;q\right)\right]$$

(2)

subject to the boundary conditions:

$$ℬ\left(\Phi ,\frac{\partial \Phi }{\partial x},\frac{\partial \Phi }{\partial t}\right)=0$$

(3)

Here, $q$ is the embedding parameter; $\Phi \left(x,t;q\right)$ is the representation of solution across $q$; $u_0\left(x,t\right)$ is the initial approximation; ħ is the auxiliary parameter; $H\left(x,t\right)$ is the auxiliary function; and $ℒ$ and $\mathcal{N}$ are the linear and nonlinear operators, respectively. Following Equation (2), the core idea of HAM is that as $q$ varies from 0 to 1, $\Phi \left(x,t;q\right)$ transforms from the initial approximation to the final solution. Mathematically, at $q=0$, $\Phi \left(x,t;0\right)=u_0\left(x,t\right)$, and at $q=1$, $\Phi \left(x,t;1\right)=\mathrm{u}\left(x,t\right)$. Now, the higher-order terms need to be determined. For that purpose, the following higher-order deformation equation is used [16]:

$$ℒ\left[u_m\left(x,t\right)-{\chi }_m{u}_{m-1}\left(x,t\right)\right]=ћH\left(x,t\right)R_m\left({\overrightarrow{u}}_{m-1}\right), m=1,2,3,\dots$$

(4)

where

$${\chi }_m=\{\begin{array}{c} 0 \mathrm{when} m=1,\\ 1 \mathrm{otherwise}\end{array}$$

(5)

and

$$R_m\left({\overrightarrow{u}}_{m-1}\right)={\frac{1}{\left(m-1\right)!}\frac{{\partial }^{m-1}\mathcal{N}\left[\Phi \left(x,t;q\right)\right]}{\partial q^{m-1}}|}_{q=0}$$

(6)

where $u_m$ for $m\ge 1$ are the higher-order terms. The derivation of Equation (4) requires the successive differentiation of the zeroth-order deformation Equation (2). The final solution can now be obtained as follows:

$$u\left(x,t\right)=u_0\left(x,t\right)+{\displaystyle \sum }_{m=1}^{\infty }{u}_m\left(x,t\right)$$

(7)

For the assessment of the solution, Equation (7) is truncated in order to have an approximate solution. In the framework of HAM, several operators and functions need to be chosen to obtain the solution. Liao [16] proposed some fundamental rules, namely the rule of solution expression, the rule of coefficient ergodicity, and the rule of solution existence, which will be discussed in the next section.

2.2. Homotopy Perturbation Method

Let us rewrite the differential equation as follows:

$$\mathcal{N}\left(u\left(x,t\right)\right)=f\left(x,t\right)$$

(8)

Now, the homotopy that satisfies [18] is constructed:

$$\left(1-q\right)\left[ℒ\left(\Phi \left(x,t;q\right)\right)-ℒ\left(u_0\left(x,t\right)\right)\right]+q\left[\mathcal{N}\left[\Phi \left(x,t;q\right)\right]-f\left(x,t\right)\right]=0$$

(9)

where the symbols carry the same meaning as mentioned in the previous section. Additionally, similar to HAM, Equation (9) shows that at $q=0$, $\Phi \left(x,t;0\right)=u_0\left(x,t\right)$, and at $q=1$, $\Phi \left(x,t;1\right)=\mathrm{u}\left(x,t\right)$. Let us now express $\Phi \left(x,t;q\right)$ as a series in terms of $q$, as follows:

$$\Phi \left(x,t;q\right)={\Phi }_0+q{\Phi }_1+q^2{\Phi }_2+q^3{\Phi }_3+q^4{\Phi }_4+\dots$$

(10)

where ${\Phi }_m$ for $m\ge 1$ are the higher-order terms. As $q\to 1$, Equation (10) produces the final solution as:

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

(11)

First, the series Equation (10) is substituted into Equation (9), and then the like powers of $q$ are equated to obtain the HPM-based solution, Equation (11). In comparison with HAM, the HPM does not contain any additional auxiliary function and auxiliary parameters, which restricts its applicability and also the rate and region of convergence of the series [16,17]. Indeed, the HPM is a special case of HAM, subject to the same set of linear and nonlinear operators and unit auxiliary function when the auxiliary parameter $ћ=-1$.

2.3. Optimal Homotopy Asymptotic Method

In some cases, HPM and HAM require several terms of the series solution in order to obtain a good approximation. Therefore, the optimal homotopy asymptotic method (OHAM) was developed with the aim of obtaining an accurate solution with just two to three terms of the series. To that end, Marinca and Herisanu [19] proposed OHAM by using asymptotic expansions of the functions and operators, which is described below. We consider the differential equation in the following form:

$$ℒ\left[u\left(x,t\right)\right]+\mathcal{N}\left[u\left(x,t\right)\right]+h\left(x,t\right)=0$$

(12)

subject to the boundary conditions:

$$ℬ\left(u,\frac{\partial u}{\partial x},\frac{\partial u}{\partial t}\right)=0$$

(13)

where symbols denote the same variables as discussed in the previous section. Following HAM, one can construct the homotopy as:

$$\left(1-q\right)\left[ℒ\left(\Phi \left(x,t,C_j;q\right)\right)+h\left(x,t\right)\right]=H\left(x,t,C_j;q\right)\left[ℒ\left(\Phi \left(x,t,C_j;q\right)\right)+h\left(x,t\right)+\mathcal{N}\left(\Phi \left(x,t,C_j;q\right)\right)\right]$$

(14)

where symbols have their usual meaning, and $C_j$ are the unknown parameters that need to be determined. The auxiliary function is defined as:

$$\begin{gathered} H\left(x,t,C_j;q\right)=0 \mathrm{for} q=0 \\ \ne 0 \mathrm{for} q\in \left(0,1\right] \end{gathered}$$

(15)

It can be verified from Equation (14) that at $q=0$, $\Phi \left(x,t,C_j;q\right)=u_0\left(x,t\right)$, and at $q=1$, $\Phi \left(x,t,C_j;q\right)=u\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. The initial approximation $u_0\left(x,t\right)$ should be determined by solving the following equation:

$$ℒ\left(u_0\left(x,t\right)\right)+h\left(x,t\right)=0$$

(16)

subject to the boundary conditions:

$$ℬ\left(u_0,\frac{\partial u_0}{\partial x},\frac{\partial u_0}{\partial t}\right)=0$$

(17)

Equation (16) can be constructed after setting $q=0$ in Equation (14). In the framework of OHAM, one of the key steps is to expand the auxiliary function in terms of $q$, as follows:

$$H\left(x,t,C_j;q\right)=q{H}_1\left(x,t,C_j\right)+q^2{H}_2\left(x,t,C_j\right)+q^3{H}_3\left(x,t,C_j\right)+\dots$$

(18)

where $H_i\left(x,t,C_j\right)$ are the auxiliary functions that depend on parameters $C_j$. Now, the final solution is expressed in the following form:

$$\Phi \left(x,t,C_j;q\right)=u_0\left(x,t\right)+{\displaystyle \sum }_{i=1}^{\infty }{u}_i\left(x,t,C_j\right)q^i$$

(19)

Substituting Equation (19) into Equation (14), and equating the like powers of $q$, the following equations are obtained [where $q^0$corresponds to Equations (16) and (17)]:

$$ℒ\left(u_1\left(x,t,C_j\right)\right)=H_1\left(x,t,C_j\right){\mathcal{N}}_0\left(u_0\left(x,t\right)\right)$$

(20)

$\mathrm{subject}\mathrm{to}$the boundary condition:

$$ℬ\left(u_1,\frac{\partial u_1}{\partial x},\frac{\partial u_1}{\partial t}\right)=0$$

(21)

For $k=2,3,4,\dots$,

$$\begin{array}{ll}ℒ[u_k\left(x,t,C_j\right)& -u_{k-1}\left(x,t,C_j\right)]\\ & =H_k\left(x,t,C_j\right){\mathcal{N}}_0\left(u_0\left(x,t\right)\right)\\ & +{\displaystyle \sum _{j=1}^{k-1}}{H}_j\left(x,t,C_j\right)[ℒ\left[u_{k-j}\left(x,t,C_j\right)\right]\\ & +{\mathcal{N}}_{k-j}\left[u_0\left(x,t\right),u_1\left(x,t,C_j\right),\dots ,u_{k-j}\left(x,t,C_j\right)\right]{\displaystyle ]} \end{array}$$

(22)

subject to the boundary conditions:

$$ℬ\left(u_k,\frac{\partial u_k}{\partial x},\frac{\partial u_k}{\partial t}\right)=0$$

(23)

where the term ${\mathcal{N}}_{k-j}\left[u_0\left(x,t\right),u_1\left(x,t,C_j\right),\dots ,u_{k-j}\left(x,t,C_j\right)\right]$ is the coefficient of $q^m$, which is obtained by expanding $\mathcal{N}\left(\Phi \left(x,t,C_j;q\right)\right)$ as follows:

$$\begin{array}{ll}\mathcal{N}\left(\Phi \left(x,t,C_j;q\right)\right)& \\ & ={\mathcal{N}}_0\left(u_0\left(x,t\right)\right)+q{\mathcal{N}}_1\left(u_0\left(x,t\right),u_1\left(x,t,C_j\right)\right)\\ & +q^2{\mathcal{N}}_2\left(u_0\left(x,t\right),u_1\left(x,t,C_j\right),u_2\left(x,t,C_j\right)\right)+\dots \end{array}$$

(24)

It may be noted that the choice of auxiliary functions $H_j\left(x,t,C_j\right)$ strongly influences the convergence of the series Equation (19). According to [20], $H_j\left(x,t,C_j\right)$ should be chosen in such a way that the product $H_j\left(x,t,C_j\right)\left[ℒ\left[u_{k-j}\left(x,t,C_j\right)\right]+{\mathcal{N}}_{k-j}\left[u_0\left(x,t\right),u_1\left(x,t,C_j\right),\dots ,u_{k-j}\left(x,t,C_j\right)\right]\right]$ and $H_j\left(x,t,C_j\right)$ are of the same form. Now, if the series Equation (19) converges at $q=1$, then we have:

$$u\left(x,t,C_j\right)=u_0\left(x,t\right)+{\displaystyle \sum }_{j=1}^{\infty }{u}_j\left(x,t,C_j\right)$$

(25)

Finally, the approximate solution can be obtained by considering a finite number of terms of the series Equation (25). The choices for parameters $C_j$ and auxiliary function will be discussed in the next section.

3. Governing Equation and Solution Methodologies

Let us consider a horizontal confined aquifer having a constant thickness, which is infinitely extended horizontally and is homogeneous and isotropic. Some of the further assumptions are that the aquifer has a single pumping well with a constant rate with respect to time and a negligibly small diameter. Further, it is assumed that the well penetrates the entire aquifer, and the hydraulic head in the aquifer, before pumping, is uniform throughout the aquifer. The continuity equation and Darcy’s law are combined together to obtain the following governing partial differential equation, representing saturated flow in a horizontal confined aquifer [21]:

$$\frac{{\partial }^{2}s}{\partial x^2}+\frac{{\partial }^{2}s}{\partial y^2}=\frac{S}{T}\frac{\partial s}{\partial t}$$

(26)

where $x$ and $y$ are the spatial variables; $t$ is the temporal variable; $s$ is the hydraulic head; and $T$ and $S$ are the transmittivity and storativity, respectively. It is convenient to convert Equation (26) into radial coordinates because of the radial symmetry of the hydraulic-head drawdown around a well. Therefore, using $r=\sqrt{x^2+y^2}$, where $r$ is the radial coordinate, the governing equation becomes:

$$\frac{{\partial }^{2}s}{\partial r^2}+\frac{1}{r}\frac{\partial s}{\partial r}=\frac{S}{T}\frac{\partial s}{\partial t}$$

(27)

The radial coordinate allows us to convert two space dimensions into a single coordinate, specifically a one-dimensional line starting from the well center at $r=0$ to the infinite extremity at $r=\infty$. The flow configuration is schematically provided in Figure 1.

The initial condition is given as:

$$s\left(r,0\right)=s_0\mathrm{for}\mathrm{all}r$$

(28)

where $s_0$ is the constant initial hydraulic head. For the boundary conditions, no hydraulic-head drawdown is assumed at the infinite extremity, i.e.,

$$s\left(\infty ,t\right)=s_0\mathrm{for}\mathrm{all}t$$

(29)

Again, a constant pumping rate is assumed at the well:

$$\underset{r\to 0}{\lim }\left(r\frac{\partial s}{\partial r}\right)=\frac{Q}{2\pi T}\mathrm{for}\mathrm{all}t$$

(30)

The second boundary condition Equation (30) is invoked by the application of Darcy’s law at the well face. The solution of Equation (27), together with the initial and boundary conditions Equations (28)–(30), provides the hydraulic head at any radial distance $r$ and time $t$ after the start of pumping. In practice, the solutions are often presented in terms of the drawdown expressed in the head $s_0-s\left(r,t\right)$.

3.1. Theis Solution

An analytical solution for the governing Equation (27), together with the given conditions Equations (28)–(30), was first provided by Theis [1], who followed an analogy with heat conduction in solids to arrive at the desired solution. The derived solution can be expressed as:

$$s_0-s\left(r,t\right)=\frac{Q}{4\pi T}\underset{v}{\overset{\infty }{{\displaystyle \int }}}\frac{\exp \left(-v\right)}{v}dv$$

(31)

where $v=\frac{r^{2}S}{4Tt}$. The solution Equation (31) is popularly known as the Theis solution. It may be noted that the governing PDE can be solved using different techniques, such as the Laplace transform, Fourier transform, similarity transformation, etc. In the mathematics literature, the integral on the right side of Equation (31) is known as the exponential integral. However, in relation to the problem considered here, it is generally called the well function and is denoted by $W\left(v\right)$. Accordingly, Equation (31) becomes:

$$s_0-s\left(r,t\right)=\frac{Q}{4\pi T}W\left(v\right)$$

(32)

The function $W\left(u\right)$ follows several interesting properties, e.g.,

$$W\left(-\infty \right)=-\infty ,W\left(0\right)=+\infty ,W\left(+\infty \right)=0,W\left(v\right)=\mathsf{\Gamma }\left(0,v\right)$$

(33)

where $\mathsf{\Gamma }\left(x,v\right)$ is the upper incomplete gamma function, defined as:

$$\mathsf{\Gamma }\left(x,v\right)=\underset{v}{\overset{\infty }{{\displaystyle \int }}}{t}^{x-1}\exp \left(-t\right)dt$$

(34)

There are two convergent series for the integral $W\left(v\right)$. One of them can be expressed as follows ([4,8]):

$$W\left(v\right)=-\gamma -\ln v+{\displaystyle \sum }_{i=1}^{\infty }\frac{{\left(-1\right)}^{i+1}{v}^i}{i i!} \left|Arg\left(v\right)\right|<\pi$$

(35)

where $\gamma$ is the Euler–Mascheroni constant, and its value (up to four decimal places) is 0.5772. For calculating the well function using Equation (35), there are some drawbacks. For example, this series produces inaccurate results for $v>2.5$ if one considers a few terms. This occurs due to the cancellation. Further, the convergence rate of the series is slow, hence, the approximation may not be desirable in practical situations. A more rapidly convergent series, attributed to S. Ramanujan, is given by [22]:

$$Ei\left(v\right)=\gamma +\ln v+\exp \left(v/2\right){\displaystyle \sum }_{k=1}^{\infty }\frac{{\left(-1\right)}^{k-1}{v}^k}{k!2^{k-1}}{\displaystyle \sum }_{n=0}^{\lfloor \left(k-1\right)/2\rfloor }\frac{1}{2n+1}$$

(36)

where $\lfloor \left(k-1\right)/2\rfloor$ denotes the floor function. The integral $W\left(v\right)$ is related to Equation (36) as follows:

$$W\left(v\right)=-Ei\left(-v\right)$$

(37)

Equation (37) provides an accurate estimate for the integral for small values of $v$. However, both Equations (35) and (36) are not useful for larger values of $v$. For large values of $v$, there is an asymptotic series approximation as follows [5]:

$$W\left(v\right)=\frac{\exp \left(-v\right)}{v}\left({\displaystyle \sum }_{k=0}^{n-1}\frac{k!}{{\left(-v\right)}^k}+\mathcal{O}\left({\left|v\right|}^{-n}\right)\right)$$

(38)

where $\mathcal{O}$ denotes the ‘Big-O’ notation. Equation (38) can be obtained by expanding the integral of $W\left(v\right)$ by parts.

The previous discussions are based on the series approximation for the integral and found to be accurate for either small or large values of $v$. However, for practical applications, one needs to have an accurate expression for $W\left(v\right)$ for a wide range of $v$. To that end, several approximations have been proposed in the literature. Here, we mention some of them. Swamee and Ojha [12] proposed the following approximation:

$$W\left(v\right)={\left[f_1{\left(v\right)}^{-7.7}+f_2\left(v\right)\right]}^{-0.13}$$

(39)

where

$$f_1\left(v\right)=\ln \left[\left(0.65+\frac{0.56146}{v}\right)\left(1+v\right)\right]$$

(40)

$$f_2\left(v\right)=v^4\exp \left(7.7v\right){\left(2+v\right)}^{3.7}$$

(41)

Barry et al. [13] offered the following full-range solution:

$$W\left(v\right)=\frac{\exp \left(-v\right)}{a_1+\left(1-a_1\right)\exp \left(-\frac{v}{1-a_1}\right)}\ln \left[1+\frac{a_1}{v}-\frac{1-a_1}{{\left(f_3\left(v\right)+a_{2}v\right)}^2}\right]$$

(42)

where

$$a_1=\exp \left(-\gamma \right), a_2=\sqrt{\frac{2\left(1-a_1\right)}{a_1\left(2-a_1\right)}}, \gamma \mathrm{is}\mathrm{the}\mathrm{Euler}–\mathrm{Mascheroni}\mathrm{constant}$$

(43)

$$f_3\left(v\right)=\frac{1}{1+v\sqrt{v}}+\frac{\widehat{a_2}\widetilde{a_2}}{1+\widetilde{a_2}}, \widehat{a_2}=\frac{\left(1-a_1\right)\left(a_1{}^2-6a_1+12\right)}{3a_1{\left(2-a_1\right)}^2{a}_2}, \widetilde{a_2}=\frac{20}{47}{v}^{\sqrt{\frac{31}{26}}}$$

(44)

In a recent study, Vatankhah [14] proposed the following approximation:

$$W\left(v\right)={\left\{{\left[\left(1+b_1{v}^{b_2}\right)\ln \left(\frac{b_3}{v}+b_4\right)\right]}^{-p}+{\left[\frac{1}{v\exp \left(v\right)}\left(\frac{v+b_5}{v+b_6}\right)\right]}^{-p}\right\}}^{-\frac{1}{p}}$$

(45)

where $p=2$, $b_1=-0.19$, $b_2=0.7$, $b_3=0.565$, $b_4=4$, $b_5=0.444$, and $b_6=1.384$.

3.2. Homotopy-Based Solutions

It may be noted that the homotopy-based methods are easier to apply to ODEs than the PDEs. This is the case because the PDE involves more than one independent variable, which makes it difficult to choose the operators and initial approximations for these methods. Therefore, first, we convert the governing PDE Equation (27) into an ODE, using the following similarity transformation:

$$v=\frac{r^{2}S}{4tT}$$

(46)

Then

$$\frac{\partial s}{\partial r}=\frac{ds}{dv}\frac{\partial v}{\partial r}=\frac{2rS}{4tT}\frac{ds}{dv} \Rightarrow \frac{1}{r}\frac{\partial s}{\partial r}=\frac{2S}{4tT}\frac{ds}{dv}$$

(47)

$$\frac{{\partial }^{2}s}{\partial r^2}=\frac{\partial }{\partial r}\left[\frac{2rS}{4tT}\frac{ds}{dv}\right]=\frac{2S}{4tT}\frac{ds}{dv}+\frac{4r^2{S}^2}{16t^2{T}^2}\frac{d^{2}s}{d{v}^2}$$

(48)

$$\frac{\partial s}{\partial t}=\frac{ds}{dv}\frac{\partial v}{\partial t}=-\frac{r^{2}S}{4t^{2}T}\frac{ds}{dv}$$

(49)

Using Equations (29)–(47), the governing Equation (27) becomes:

$$v\frac{d^{2}s}{d{v}^2}+\left(1+v\right)\frac{ds}{dv}=0$$

(50)

Following Equation (46), the conditions Equations (28) and (29) become:

$$s\left(v\to \infty \right)=s_0$$

(51)

and Equation (30) changes to:

$$\underset{r\to 0}{\lim }\left(r\frac{\partial s}{\partial r}\right)=\frac{Q}{2\pi T}\Rightarrow \underset{v\to 0}{\lim }\left(v\frac{ds}{dv}\right)=\frac{Q}{4\pi T}$$

(52)

To be in line with Theis’s method, we intend to obtain the solution as $s_0-s\left(v\right)$. For that, we use the transformation:

$$s\left(v\right)=s_0+\overline{s}\left(v\right)$$

(53)

Accordingly, Equation (50) becomes:

$$v\frac{d^2\overline{s}}{d{v}^2}+\left(1+v\right)\frac{d\overline{s}}{dv}=0$$

(54)

The boundary conditions (51) and (52) become:

$$\overline{s}\left(v\to \infty \right)=0$$

(55)

$$\underset{v\to 0}{\lim }\left(v\frac{d\overline{s}}{dv}\right)=-\frac{Q}{4\pi T}$$

(56)

3.2.1. HAM-Based Solution

Here, we apply HAM to Equation (54) together with the conditions Equations (55) and (56). Following the discussion in Section 2.1, we consider the nonlinear operator for the problem as follows:

$$\mathcal{N}\left[\Phi \left(v;q\right)\right]=v\frac{{\partial }^2\Phi \left(v;q\right)}{\partial v^2}+\left(1+v\right)\frac{\partial \Phi \left(v;q\right)}{\partial v}$$

(57)

It may be noted that Equation (57) is the original governing equation. Using Equation (57), terms $R_m$ are obtained as follows:

$$R_m\left({\overrightarrow{\overline{s}}}_{m-1}\right)=v\frac{d^2{\overline{s}}_{m-1}}{d{v}^2}+\left(1+v\right)\frac{d{\overline{s}}_{m-1}}{dv}$$

(58)

Now, we define the base functions to represent the solution of the governing equation. To that end, the following set of base functions is chosen for the considered problem:

$$\left\{\exp \left(-v\right)\ln \left(v\right)+\left[W\left(v\right)+\exp \left(-v\right)\right]v^m | m=0,1,2,\dots \right\}$$

(59)

so that

$$\overline{s}\left(v\right)=a_0\exp \left(-v\right)\ln \left(v\right)+\left[W\left(v\right)+\exp \left(-v\right)\right]{\displaystyle \sum }_{n=0}^{\infty }{b}_n{\left(v\right)}^n$$

(60)

where $a_0$ and $b_m$ are the coefficients of the series. Equation (60) provides the so-called rule of solution expression. Following the rule of solution expression, the linear operator and the initial approximation are chosen, respectively, as follows:

$$ℒ\left[\Phi \left(v;q\right)\right]=\frac{{\partial }^2\Phi \left(v;q\right)}{\partial v^2} \mathrm{with}\mathrm{the}\mathrm{property} ℒ\left[C_{2}v+C_3\right]=0$$

(61)

$${\overline{s}}_0\left(v\right)=-\frac{Q}{4\pi T}\exp \left(-v\right)\ln \left(v\right)$$

(62)

where $C_2$ and $C_3$ are integral constants. It may be noted that the initial approximation Equation (62) is chosen in accordance with the given boundary conditions Equations (55) and (56). Using Equation (61), the higher-order terms can be obtained from Equation (62) as follows:

$${\overline{s}}_m\left(v\right)={\chi }_m{\overline{s}}_{m-1}\left(v\right)+ћ{{\displaystyle \int }}_0^v{{\displaystyle \int }}_0^{v}H\left(x\right)R_m\left({\overrightarrow{\overline{s}}}_{m-1}\right)dxdy+C_{2}v+C_3, m=1,2,3,\dots$$

(63)

where $R_m$ is given by Equation (68), and constants $C_2$ and $C_3$ can be determined from the boundary conditions for the higher-order deformation equations given by Equation (63).

Now, the rule of coefficient ergodicity determines the auxiliary function $H\left(v\right)$. It may be noted that the function $H\left(v\right)$ can be of many forms, based on the rule of solution expression Equation (60). However, for simplicity, one can select $H\left(v\right)=1$ in order to avoid unnecessary difficulty in computation. Further, this is consistent with the theory given in [23]. Finally, the approximate solution can be obtained as follows:

$$\overline{s}\left(v\right)=s_0-s\left(v\right)\approx {\overline{s}}_0\left(v\right)+{\displaystyle \sum }_{n=1}^M{\overline{s}}_n\left(v\right)$$

(64)

3.2.2. HPM-Based Solution

Following the discussion in Section 2.2, we select the linear and nonlinear operators as follows:

$$ℒ\left(\Phi \left(v;q\right)\right)=\frac{{\partial }^2\Phi \left(v;q\right)}{\partial v^2}$$

(65)

$$\mathcal{N}\left(\Phi \left(v;q\right)\right)=v\frac{{\partial }^2\Phi \left(v;q\right)}{\partial v^2}+\left(1+v\right)\frac{\partial \Phi \left(v;q\right)}{\partial v}$$

(66)

Using these operators, Equation (9) takes on the form:

$$\left(1-q\right)\left[\frac{{\partial }^2\Phi \left(v;q\right)}{\partial v^2}-\frac{d^2{\overline{s}}_0}{d{v}^2}\right]+q\left[v\frac{{\partial }^2\Phi \left(v;q\right)}{\partial v^2}+\left(1+v\right)\frac{\partial \Phi \left(v;q\right)}{\partial v} \right]=0$$

(67)

Substituting the series expression Equation (10) into Equation (67), we have:

$$\begin{array}{ll}\left(1-q\right)\lceil \frac{{\partial }^2\Phi \left(v;q\right)}{\partial v^2}-\frac{d^2{\overline{s}}_0}{d{v}^2}\rceil & \\ & =\left(\frac{d^2{\Phi }_0}{d{v}^2}-\frac{d^2{\overline{s}}_0}{d{v}^2}\right)+q\left(\frac{d^2{\Phi }_1}{d{v}^2}-\left(\frac{d^2{\Phi }_0}{d{v}^2}-\frac{d^2{\overline{s}}_0}{d{v}^2}\right)\right)+q^2\left(\frac{d^2{\Phi }_2}{d{v}^2}-\frac{d^2{\Phi }_1}{d{v}^2}\right)+q^3\left(\frac{d^2{\Phi }_3}{d{v}^2}-\frac{d^2{\Phi }_2}{d{v}^2}\right)\\ & +q^4\left(\frac{d^2{\Phi }_4}{d{v}^2}-\frac{d^2{\Phi }_3}{d{v}^2}\right)+\dots \end{array}$$

(68)

$$\begin{array}{ll}v\frac{{\partial }^2\Phi \left(v;q\right)}{\partial v^2}+\left(1+v\right)\frac{\partial \Phi \left(v;q\right)}{\partial v}& \\ & =\left[v\frac{d^2{\Phi }_0}{d{v}^2}+\left(1+v\right)\frac{d{\Phi }_0}{dv}\right]+q\left[v\frac{d^2{\Phi }_1}{d{v}^2}+\left(1+v\right)\frac{d{\Phi }_1}{dv}\right]\\ & +q^2\left[v\frac{d^2{\Phi }_2}{d{v}^2}+\left(1+v\right)\frac{d{\Phi }_2}{dv}\right]+q^3\left[v\frac{d^2{\Phi }_3}{d{v}^2}+\left(1+v\right)\frac{d{\Phi }_3}{dv}\right]+\dots \end{array}$$

(69)

Using the boundary conditions $\overline{s}\left(v\to \infty \right)=0$ and $\underset{v\to 0}{\lim }\left(v\frac{d\overline{s}}{dv}\right)=-\frac{Q}{4\pi T}$, we have:

$${\Phi }_0\left(v\to \infty \right)=0,{\Phi }_1\left(v\to \infty \right)=0,{\Phi }_2\left(v\to \infty \right)=0,\dots .$$

(70)

$$\underset{v\to 0}{\lim }\left(v\frac{d{\Phi }_0}{dv}\right)=-\frac{Q}{4\pi T}, \underset{v\to 0}{\lim }\left(v\frac{d{\Phi }_1}{dv}\right)=0, \underset{v\to 0}{\lim }\left(v\frac{d{\Phi }_2}{dv}\right)=0, \dots$$

(71)

Using Equations (70) and (71) and equating the like powers of $q$ in Equation (67), the following system of differential equations is obtained:

$$\frac{d^2{\Phi }_0}{d{v}^2}-\frac{d^2{\overline{s}}_0}{d{v}^2}=0 \mathrm{subject}\mathrm{to} {\Phi }_0\left(v\to \infty \right)=0, \underset{v\to 0}{\lim }\left(v\frac{d{\Phi }_0}{dv}\right)=-\frac{Q}{4\pi T}$$

(72)

$$\frac{d^2{\Phi }_1}{d{v}^2}-\left(\frac{d^2{\Phi }_0}{d{v}^2}-\frac{d^2{\overline{s}}_0}{d{v}^2}\right)+v\frac{d^2{\Phi }_0}{d{v}^2}+\left(1+v\right)\frac{d{\Phi }_0}{dv}=0 \mathrm{subject}\mathrm{to} {\Phi }_1\left(v\to \infty \right)=0, \underset{v\to 0}{\lim }\left(v\frac{d{\Phi }_1}{dv}\right)=0$$

(73)

$$\frac{d^2{\Phi }_2}{d{v}^2}-\frac{d^2{\Phi }_1}{d{v}^2}+v\frac{d^2{\Phi }_1}{d{v}^2}+\left(1+v\right)\frac{d{\Phi }_1}{dv}=0 \mathrm{subject}\mathrm{to} {\Phi }_2\left(v\to \infty \right)=0, \underset{v\to 0}{\lim }\left(v\frac{d{\Phi }_2}{dv}\right)=0$$

(74)

$$\frac{d^2{\Phi }_3}{d{v}^2}-\frac{d^2{\Phi }_2}{d{v}^2}+v\frac{d^2{\Phi }_2}{d{v}^2}+\left(1+v\right)\frac{d{\Phi }_2}{dv}=0 \mathrm{subject}\mathrm{to} {\Phi }_3\left(v\to \infty \right)=0, \underset{v\to 0}{\lim }\left(v\frac{d{\Phi }_3}{dv}\right)=0$$

(75)

$$\frac{d^2{\Phi }_4}{d{v}^2}-\frac{d^2{\Phi }_3}{d{v}^2}+v\frac{d^2{\Phi }_3}{d{v}^2}+\left(1+v\right)\frac{d{\Phi }_3}{dv}=0 \mathrm{subject}\mathrm{to} {\Phi }_4\left(v\to \infty \right)=0, \underset{v\to 0}{\lim }\left(v\frac{d{\Phi }_4}{dv}\right)=0$$

(76)

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

$$\frac{d^2{\Phi }_m}{d{v}^2}=\left(1-v\right)\frac{d^2{\Phi }_{m-1}}{d{v}^2}-\left(1+v\right)\frac{d{\Phi }_{m-1}}{dv} \mathrm{subject}\mathrm{to} {\Phi }_m\left(v\to \infty \right)=0, \underset{v\to 0}{\lim }\left(v\frac{d{\Phi }_m}{dv}\right)=0 \mathrm{for} m\ge 2$$

(77)

The initial approximation can be chosen as ${\Phi }_0=-\frac{Q}{4\pi T}\exp \left(-v\right)\ln v$. Using this initial approximation, we can solve the equations iteratively using symbolic software. Finally, the HPM-based solution can be approximated as follows:

$$\overline{s}\left(v\right)\approx {\displaystyle \sum }_{i=0}^M{\Phi }_i$$

(78)

3.2.3. OHAM-Based Solution

When dealing with the problem using OHAM, it is observed that due to boundary conditions, it may not be possible to solve the governing Equation (54) using OHAM. To that end, first, we convert the boundary value problem Equation (54) into an ODE, as follows:

$$v\frac{d^2\overline{s}}{d{v}^2}+\left(1+v\right)\frac{d\overline{s}}{dv}=0\Rightarrow \frac{d}{dv}\left(v\frac{d\overline{s}}{dv}\right)+v\frac{d\overline{s}}{dv}=0\Rightarrow \frac{du}{dv}+u=0$$

(79)

subject to

$$\underset{v\to 0}{\lim }\left(u\right)=-\frac{Q}{4\pi T}$$

(80)

Solving Equation (79) together with Equation (80), one obtains:

$$v\frac{d\overline{s}}{dv}+\frac{Q}{4\pi T}\exp \left(-v\right)=0 \mathrm{subject}\mathrm{to} \overline{s}(\infty )=0$$

(81)

For applying OHAM, Equation (81) can be written in the following form:

$$ℒ\left[\overline{s}\left(v\right)\right]+\mathcal{N}\left[\overline{s}\left(v\right)\right]+h\left(v\right)=0$$

(82)

We select $ℒ\left[\overline{s}\left(v\right)\right]=\frac{d\overline{s}}{dv}$, $\mathcal{N}\left[\overline{s}\left(v\right)\right]=\left(v-1\right)\frac{d\overline{s}}{dv}$, and $h\left(v\right)=\frac{Q}{4\pi T}\exp \left(-v\right)$. Using these expressions, solution of the zeroth-order representation Equation (16) becomes:

$${\overline{s}}_0\left(v\right)=\frac{Q}{4\pi T}\exp \left(-v\right)$$

(83)

The following relation is obtained for the expressions of ${\mathcal{N}}_0$, ${\mathcal{N}}_1$, ${\mathcal{N}}_2$, etc.:

$$\left(v-1\right)\frac{d\overline{s}}{dv}=\left[\left(v-1\right)\frac{d{\overline{s}}_0}{dv}\right]+q\left[\left(v-1\right)\frac{d{\overline{s}}_1}{dv}\right]+q^2\left[\left(v-1\right)\frac{d{\overline{s}}_2}{dv}\right]+q^3\left[\left(v-1\right)\frac{d{\overline{s}}_3}{dv}\right]+\dots$$

(84)

Using Equation (84), the first-order representation Equation (20) reduces to:

$$\frac{d{\overline{s}}_1}{dv}=H_1\left(v,C_i\right)\left\{\left(v-1\right)\frac{d{\overline{s}}_0}{dv}\right\} \mathrm{subject}\mathrm{to} {\overline{s}}_1\left(v\to \infty \right)=0$$

(85)

The auxiliary function is chosen as $H_1\left(v,C_i\right)=C_1+C_2\exp \left(-v\right)$. Putting $k=2$ and $H_2\left(v,C_i\right)=C_3+C_4\exp \left(-v\right)$, the second-order representation Equation (22) becomes:

$$\frac{d{\overline{s}}_2}{dv}=\frac{d{\overline{s}}_1}{dv}+\left(C_3+C_4\exp \left(-v\right)\right)\left[\left(v-1\right)\frac{d{\overline{s}}_0}{dv}\right]+\left(C_1+C_2\exp \left(-v\right)\right)\left[\frac{d{\overline{s}}_1}{dv}+\left(v-1\right)\frac{d{\overline{s}}_1}{dv}\right] \mathrm{subject}\mathrm{to} {\overline{s}}_2\left(v\to \infty \right)=0$$

(86)

The first three terms of the OHAM-based solution can be obtained by solving the above equations using symbolic computation software, such as MATLAB. Further, we restrict our analysis up to $k=2$, as the objective is to obtain an accurate solution with just a few terms of the OHAM-based series. Finally, the approximate solution can be found as:

$$\overline{s}\left(v\right)\approx {\overline{s}}_0\left(v\right)+{\overline{s}}_1\left(v,C_1,C_2\right)+{\overline{s}}_2\left(v,C_1,C_2,C_3,C_4\right)$$

(87)

where the terms are given by Equations (83), (85), and (86).

4. Results and Discussion

Here, first, we discuss the validity of the approximations for the well function. Then, the HAM-, HPM-, and OHAM-based approximations are verified with the numerical solution to the problem. Finally, the proposed solutions are compared with the existing approximations in order to have a numerical assessment.

4.1. Validation of the Well Function’s Approximations

In Figure 2, we validate the series approximation given by Equation (35) by comparing it with the numerical solution of the corresponding exponential integral. The numerical values of $W\left(v\right)$ were calculated using the MATLAB script ‘integral’, which uses the global adaptive quadrature rule [24]. We compared 10, 20, and 40 terms of the series Equation (35) and observe that one needs more than 40 terms of the series to obtain an accurate approximation for $v$ up to 6.5. Further, the convergence rate is too slow as we increase the order of approximation.

Figure 3 compares the numerical solution and the series approximation (38) of $W\left(v\right)$, considering some terms of the series. The series is an asymptotic expansion that is valid for large values of $v$. As expected, it is observed from the figure that the series performs well only for large values of $v$. However, it produces inaccurate results for small values of $v$. For a comparative assessment, we plot together the numerical solution and both the series Equations (35) and (38) in Figure 4. It is observed from the figure that one of them performs well for smaller values, and the other is accurate only for large values of $v$.

4.2. Numerical Convergence and Validation of the HAM-Based Solution

In the framework of HAM, the auxiliary parameter $ћ$ plays the key role in determining the convergence of the series solution. An optimal value for the parameter can be obtained by minimizing the following squared residual error of Equation (54), which can be calculated as follows:

$${\Delta }_m=\underset{v\in \Omega }{{\displaystyle \int }}{\left(\mathcal{N}\left[\overline{s}\left(v\right)\right]\right)}^{2}dv$$

(88)

where $\Omega$ is the domain of the equation. The HAM-based series solution converges when the corresponding residual error Equation (88) becomes zero. A test case is considered here, where the parameters were chosen as $Q=4\times {10}^{-3}$ m³s⁻¹, $T=0.0023$ m²s⁻¹, and $S=7.5\times {10}^{-4}$. Using these parameter values, we assessed the HAM-based solution by calculating the squared residual errors for different orders of approximations. The squared residual errors for different orders of approximation are plotted in Figure 5, where it is seen that the error decreases with the increasing order of approximation. Thus, the numerical convergence was established, and the choice of operators and parameters was validated. For a quantitative assessment, the numerical results are also reported in

Table 1. Squared residual error (${\Delta }_m$) and computational time versus different orders of approximation ($\mathrm{m}$) for the selected case.

$\mathbf{Order}\mathbf{of}\mathbf{Approximation}\left(\mathbf{m}\right)$	$\mathbf{Squared}\mathbf{Residual}\mathbf{Error}\left({\mathbf{\Delta }}_{\mathit{m}}\right)$	Computational Time (s)
2	7.15 $\times$ 10−4	0.212
4	5.23 $\times$ 10−5	1.135
6	1.98 $\times$ 10−5	2.481
8	1.01 $\times$ 10−5	5.074
10	8.86 $\times$ 10−6	6.843
12	6.23 $\times$ 10−6	10.149

, along with the computational time taken by the computer to produce the corresponding order of approximations. It can be seen from the table that even though HAM is an analytical series approximation technique, it still does not involve time complexity. On the other hand, for the selected case, we compared the Theis solution (with the integral computed using the MATLAB script ‘integral’, which uses the global adaptive quadrature rule [24]) and the 10th-order HAM-based approximate solution in Figure 6. It may be noted that the same method was used for the numerical solution related to HPM and OHAM. An excellent agreement is found between the computed and observed values. Additionally, for a comparative idea, 4th-, 7th-, and 10th-order approximations were considered and compared with the Theis solution, as shown in

Table 2. Comparison between HAM-based approximation and numerical solution for the selected case.

$\mathit{u}$	Numerical Solution	HAM-Based Approximation
		4th Order	7th Order	10th Order
0.1	2.523 $\times$ 10−1	3.009 $\times$ 10−1	2.821 $\times$ 10−1	2.726 $\times$ 10−1
1	3.036 $\times$ 10−2	3.995 $\times$ 10−2	3.548 $\times$ 10−2	3.305 $\times$ 10−2
2	6.768 $\times$ 10−3	9.302 $\times$ 10−3	8.011 $\times$ 10−3	7.228 $\times$ 10−3
3	1.806 $\times$ 10−3	2.560 $\times$ 10−3	2.146 $\times$ 10−3	1.886 $\times$ 10−3
4	5.230 $\times$ 10−4	7.603 $\times$ 10−4	6.213 $\times$ 10−4	5.377 $\times$ 10−4
5	1.589 $\times$ 10−4	2.359 $\times$ 10−4	1.881 $\times$ 10−4	1.627 $\times$ 10−4
6	4.983 $\times$ 10−5	7.530 $\times$ 10−5	5.870 $\times$ 10−5	5.150 $\times$ 10−5
7	1.598 $\times$ 10−5	2.453 $\times$ 10−5	1.871 $\times$ 10−5	1.689 $\times$ 10−5
8	5.213 $\times$ 10−6	8.111 $\times$ 10−6	6.066 $\times$ 10−6	5.688 $\times$ 10−6
9	1.723 $\times$ 10−6	2.713 $\times$ 10−6	1.993 $\times$ 10−6	1.955 $\times$ 10−6
10	5.753 $\times$ 10−7	9.160 $\times$ 10−7	6.617 $\times$ 10−7	6.814 $\times$ 10−7

. It can be observed from the table that the higher the order of approximation, the better the accuracy. All computations were performed using the BVPh 2.0 package developed by [25]. A flowchart containing the steps of HAM for the present problem is provided in Figure 7. The theoretical convergence analysis is provided in Appendix A.

4.3. Validation of HPM-Based Solution

For the selected case, the HPM-based analytical solution was validated over the solution given by [1], where the integral is performed numerically. It may be noted that, unlike HAM, HPM solutions are often valid only within a small domain [16]. Therefore, for our case here, we considered the domain $v\in \left[0.1,2\right]$. After assessing the solution within this domain, we observe that the series with four terms is more accurate than the other lower-order terms. The comparison between the Theis solution and the four terms of HPM approximation is given in Figure 8, where it can be seen that as the domain increases, the accuracy of the solution decreases. This is a default problem with HPM, as the methodology does not contain any convergence-control parameter like HAM. Indeed, one may try with different combinations of initial approximations and operators, which might then work in producing a more accurate solution to the problem.

Table 3. Comparison between four terms of the HPM-based approximation and Theis solution for the selected case.

$\mathit{u}$	Numerical Solution	Four Terms of the HPM-Based Approximation
0.1	2.523 $\times$ 10−1	3.135 $\times$ 10−1
0.3	1.253 $\times$ 10−1	1.533 $\times$ 10−1
0.5	7.747 $\times$ 10−2	9.362 $\times$ 10−2
0.7	5.173 $\times$ 10−2	6.246 $\times$ 10−2
0.9	3.601 $\times$ 10−2	4.346 $\times$ 10−2
1.1	2.574 $\times$ 10−2	3.097 $\times$ 10−2
1.3	1.875 $\times$ 10−2	2.237 $\times$ 10−2
1.5	1.384 $\times$ 10−2	1.563 $\times$ 10−2
1.7	1.033 $\times$ 10−2	8.139 $\times$ 10−3
2.0	6.768 $\times$ 10−3	1.298 $\times$ 10−2

shows a numerical comparison between the HPM-based values and the Theis solution. For the convenience of the readers, a flowchart containing the steps of HPM is provided in Figure 9.

4.4. Validation of OHAM-Based Solution

It can be seen from Equation (87) that the OHAM-based solution contains parameters $C_i$, which need to be calculated. For that purpose, one can construct the residual as follows:

$$R\left(v,C_i\right)=ℒ\left[{\overline{s}}_{OHAM}\left(v,C_i\right)\right]+\mathcal{N}\left[{\overline{s}}_{OHAM}\left(v,C_i\right)\right]+h\left(v\right), i=1,2,\dots ,s$$

(89)

where ${\overline{s}}_{OHAM}\left(v,C_i\right)$ is the approximate solution. When $R\left(v,C_i\right)=0$, ${\overline{s}}_{OHAM}\left(v,C_i\right)$ becomes the exact solution to the equation. One of the ways to obtain the optimal $C_i$ is the minimization of squared residual error, i.e.,

$$J\left(C_i\right)={{\displaystyle \int }}_{vϵD}{R}^2\left(v,C_i\right)du, i=1,2,\dots ,s$$

(90)

where $D=\left[0.1,10\right]$ is the domain of the problem. The minimization of Equation (90) leads to a system of algebraic equations as follows:

$$\frac{\partial J}{\partial C_1}=\frac{\partial J}{\partial C_2}=\dots =\frac{\partial J}{\partial C_s}=0$$

(91)

One can obtain the optimal values of parameters by solving this system of equations numerically. Here, we used the MATLAB script fminsearch, which minimizes an unconstrained multivariable function. It was found that only a three-term solution produces accurate values. Therefore, a three-term OHAM solution was computed and compared with the Theis solution, as shown in Figure 10. It can be seen that just three terms of the OHAM-based series agree well with the corresponding analytical solution to the problem. For a quantitative assessment, we also compare the numerical values of the solutions in

Table 4. Comparison between three terms of the OHAM-based approximation and Theis solution for the selected case.

$\mathit{u}$	Numerical Solution	Three Terms of OHAM-Based Approximation
0.1	2.523 $\times$ 10−1	2.466 $\times$ 10−1
1	3.036 $\times$ 10−2	3.143 $\times$ 10−2
2	6.768 $\times$ 10−3	6.539 $\times$ 10−3
3	1.806 $\times$ 10−3	2.112 $\times$ 10−3
4	5.230 $\times$ 10−4	8.502 $\times$ 10−4
5	1.589 $\times$ 10−4	3.066 $\times$ 10−4
6	4.983 $\times$ 10−5	1.000 $\times$ 10−4
7	1.598 $\times$ 10−5	3.013 $\times$ 10−5
8	5.213 $\times$ 10−6	8.461 $\times$ 10−6
9	1.723 $\times$ 10−6	2.283 $\times$ 10−6
10	5.753 $\times$ 10−7	5.978 $\times$ 10−7

. A flowchart containing the steps of OHAM is provided in Figure 11.

4.5. Comparison between Different Approximations

In this section, we compare the approximations given by series approximations Equations (35), (38), (39), (42), (45), (64), (78) and (87). For that purpose, we considered the same test parameters, i.e., $Q=4\times {10}^{-3}$ m³s⁻¹, $T=0.0023$ m²s⁻¹, and $S=7.5\times {10}^{-4}$. Moreover, for each of the cases, we computed the well function numerically using ‘integral’ of MATLAB to obtain the main solution. The series Equation (35) with 40 terms, Equation (38) with 10 terms, 10th-order HAM-based solution, four-term HPM solution, and three-term OHAM solution were considered. Importantly, it may be noted that the numerical values of solutions are very small, which can make the computations ill-posed or produce numerical instabilities. To that end, logarithmic form for the error was considered. Specially, we checked the performances of the approximations by calculating the percentage error as PE (%) = 100 $\times \frac{\left(\ln W_{num}-\ln W_{apprx}\right)}{\ln W_{num}}$, where $W_{num}$ and $W_{apprx}$ are the values of $W\left(v\right)$ obtained from the Theis solution and the corresponding approximation, respectively. The percentage errors were calculated for the approximations and compared, as shown in Figure 12. It can be seen that among the series approximations, HAM- and OHAM-based approximations provide accurate approximations for the problem. On the other hand, the HPM-based solution is shown only within a small domain, as the solution provides accurate values there. Further, series approximation given by [12,13,14] are reasonably accurate within the domain. The different homotopy-based methods provide solutions that are valid within a certain range of the domain. The HAM- and OHAM-based approximations are more accurate and valid for larger domains, as they contain convergence-control parameters, which monitor the rate and radius of convergence of the series solutions. Further, the OHAM-based solution is more preferable due to its ability to provide an accurate approximation with just two–three terms of the series. Finally, it is concluded that while the homotopy-based methods do not produce as accurate solutions as do the closed-form formulae available in the literature, they are better than the series expansions and also may be improved further using different sets of base functions, linear operators, and initial approximations. It may also be noted that the proposed study differs from the existing empirical formula-based work from the viewpoint of its derivation, which starts from the governing differential equation. Therefore, the approach is flexible to use when the flow configuration is different, and the model parameters vary.

5. Concluding Remarks

Theis first derived the analytical solution for the governing partial differential equation representing saturated unsteady flow in a horizontal confined aquifer. The solution was presented in the form of an exponential integral and is popularly known as the Theis solution. This exponential integral, also called the well function, requires an approximation that is simple to adapt and reliable for practical applications. Researchers have provided several approximations for the well function using series expansion, numerical approximation, etc. However, most of these approximations are valid either within a restricted domain or for a large number of terms of the series. This work directly solves the governing PDE analytically after converting it to a BVP using a similarity transformation. The HAM, HPM, and OHAM are used to obtain the solution analytically in the form of a series. It is found that both HAM and OHAM provide accurate solutions to the problem for a sufficiently large domain. Specifically, ten terms of the HAM series and three terms of the OHAM series provide the approximation well when compared with the original Theis solution. On the other hand, HPM also produces an accurate solution within a restricted domain—this is desirable, as the methodology does not contain a convergence-control parameter like HAM. Several series and numerical approximations are validated using a test case. Further, the proposed approximations are compared with the existing series and numerical approximations by calculating the percentage error. It is seen that the proposed approximations are reliable for predicting the drawdown. Because this study derives the solutions starting directly from the governing equation, the methodology adopted here can be extended to other kinds of flow configurations to find efficient solutions.