A mathematical model is developed to describe the groundwater inflow into a tunnel in a multi-layer aquifer system. Based on the model, the closed-form solution is derived to estimate the groundwater flow rate entering the multi-layer tunnel during progressive drilling. The solution has an integrand not only consisting of the product and square of the Bessel functions but also having a singularity at the origin. A unified numerical approach is proposed to evaluate the solution with accuracy to five decimal places. This approach includes a singularity removal scheme, the Gaussian quadrature, and the Shanks method. For a multi-layer formation, the results obtained from the solution based on the equivalent hydraulic conductivity and the newly derived solution differ significantly. This solution is capable of estimating the maximum flow rate inside the horizontal tunnel, and thus can be used as a tool for designing the drainage tunnel system in a multi-layer formation.