We present an effective finite difference formulation for implementing and modeling multiple borehole heat exchangers (SHE) in the general 3-D coupled heat and flow transport code SHEMAT. The BHE with arbitrary length can be either coaxial or double U-shaped. It is particularly suitable for modeling deep BHEs which contain varying pipe diameters and materials. Usually, in numerical simulations, a fine discretization of the BHE assemblage is required, due to the large geometric aspect ratios involved. This yields large models and long simulation times. The approach avoids this problem by considering heat transport between fluid and the soil through pipes and grout via thermal resistances. Therefore, the simulation time can be significantly reduced. The coupling with SHEMAT is realized by introducing an effective heat generation. Due to this connection, it is possible to consider heterogeneous geological models, as well as the influence of groundwater flow. This is particularly interesting when studying the long term behavior of a single BHE or a SHE field. Heating and cooling loads can enter the model with an arbitrary interval, e.g. from hourly to monthly values. When dealing with large BHE fields, computing times can be further significantly reduced by focusing on the temperature field around the BHEs, without explicitly modeling inlet and outlet temperatures. This allows to determine the possible migration of cold and warm plumes due to groundwater flow, which is of particular importance in urban areas with a high BHE installation density. The model is validated against the existing BHE modeling codes EWS and EED. A comparison with monitoring data from a deep BHE in Switzerland shows a good agreement. Synthetic examples demonstrate the field of application of this model. (C) 2012 Elsevier Ltd. All rights reserved.