This paper introduces a new methodology for the generation of thermal response factors of geothermal bore fields using the concept of g-functions introduced by Eskilson. Boreholes are divided into segments to consider the variation of the heat extraction rates along the length of the boreholes and the analytical finite line source (FLS) solution is used to calculate the temperature variations at the wall of each borehole segment along the axial direction. The proposed methodology accounts for the time variation of the heat extraction rates among boreholes and along the length of individual boreholes to obtain a uniform borehole wall temperature equal for all boreholes in accordance with the original boundary conditions proposed by Eskilson. In addition, the methodology is generalized to account for boreholes of different lengths and buried depths. g-Functions calculated with the proposed methodology are compared to the numerical technique used by Eskilson to derive the g-functions for fields of 1 to 12 x 12 boreholes. The difference between the two models is within 5% for all studied bore fields, except for fields of boreholes located on a single row. The variation of the heat extraction rates of individual boreholes along their length as well as in time also showed good agreement with the numerical model. It is shown that using 12 borehole segments is adequate to calculate the g-functions in most practical cases. For instance, the error on the g-function of a 10 x 10 bore field calculated using 12 borehole segments is 2.2% after 20 years and 4.7% at steady-state. (C) 2013 Elsevier Ltd. All rights reserved.