Improved numerical evaluation of the radial groundwater flow equation

The closed-form solution for the hydraulic head, derived for the radial groundwater flow equation subject to the constant-head boundary condition at the wellbore, is in an integral form that covers an integration range from zero to infinity. The integral is difficult to evaluate due to the integrand not only consisting of the product and the square of the Bessel functions but also having a singularity at the origin. A unified numerical method is proposed to evaluate the solution with accuracy to five decimal places and for very wide ranges of dimensionless distances and times. This approach includes a singularity removal scheme, Newton's method, the Gaussian quadrature, and Shanks' method. It gives the dimensionless heads in tabular forms with better accuracy while comparing to those given by other approaches. A formula describing the flow rate across the wellbore, derived by the authors based on Darcy's law, is proved to equal those presented by Jaeger [Proc Royal Soc Edinburgh 61 (1942) 223] and Jacob and Lohman [Am Geo Union 33 (4) (1952) 559]. The same singularity removal scheme and the Gaussian quadrature are also employed to evaluate the wellbore flow rate. Computed values of dimensionless flow rate versus dimensionless time expressed in tabular forms are also correct to five decimal places. Both the tabular results of dimensionless hydraulic head and dimensionless flow rate may be useful in engineering applications.