In most ground-coupled heat pump systems, Borehole Heat Exchangers (BHE) represent the typical engineering solution for utilizing renewable energy from the ground. The design of a complex BHE field is a challenging task, due the inherent transient nature of the thermal interaction between the heat exchangers and the surrounding soil. A computation effective method for solving the 3D transient conduction equation describing the ground response to a variable heat load profile is the temporal superposition of pre-calculated temperature response factors or g-functions. In this study Comsol heat conduction models have been developed to calculate g-function values for a borehole field with 64 boreholes. The aim of the investigation is to get an insight on the numerical generation of temperature transfer functions and to some extent provide new information on the Finite Line Source method for analytically generated g-functions as well as on those existing behind existing design software such as EED. The results generally showed a good agreement in lower time ranges. Further in time, the Comsol model revealed to be influenced either by the domain dimensions or the simulation end time.