Computer simulation of borehole ground heat exchangers used in geothermal heat pump systems was conducted using three-dimensional implicit finite difference method with rectangular coordinate system. Each borehole was approximated by a square column circumscribed by the borehole radius. Borehole loading profile calculated numerically based on the prescribed borehole temperature profile under quasi-steady state conditions was used to determine the ground temperature and the borehole temperature profile. The two coupled solutions were solved iteratively at each time step. The simulated ground temperature was calibrated using a cylindrical source model by adjusting the grid spacing and adopting a load factor of 1.047 in the difference equation. With constant load applied to a single borehole, neither the borehole temperature nor the borehole loading was constant along the borehole. The ground temperature profiles were not similar at different distances from the borehole. This meant that a single finite difference scheme was not sufficient to estimate the performance of a borefield by superposition. The entire borefield should be discretized simultaneously. Comparison was made between the present method and the finite line source model with superposition. The discrepancies between the results from the two methods increased with the scale of borefield. The introduction of time schedule revealed a discrepancy between the load applied to the ground heat exchanger and that transferred from the borehole to the ground, which was usually assumed to be the same when using analytical models. Hence, in designing a large borefield, the present method should give more precise results in dynamic simulation. (C) 2007 Elsevier Ltd. All rights reserved.