This article presents yet another methodology for the calculation of dimensionless thermal response factors for vertical borehole ground heat exchanger (GHX) arrays, which is a concept introduced by Eskilson (1987). The presented method is based on a well-known solution to an analogous problem in the field of well hydraulics. This solution method, known mathematically as an incomplete Bessel function, and known in the field of well hydraulics as the 'leaky aquifer function', describes the hydraulic head distribution in an aquifer with predominantly radial flow to a well combined with vertical 'leakage' from geologic layers above and below the pumped aquifer. The solution is adapted to model heat transfer from an array of arbitrarily-placed vertical boreholes of finite depth. With proper expression of parameters in the incomplete Bessel function, we show that g-functions of previous researchers can be approximated. The proposed method has been implemented into Matlab and Excel/VBA for g-function generation and monthly GHX simulation.