In the modeling of shallow geothermal systems, the accurate representation of the borehole field configuration is important for a proper estimation of the long term thermal behavior of borehole field systems. Modeling tools based on the so-called g-functions method, utilized for the design of borehole fields, assume that boreholes are vertical. This is a limitation since this condition might not apply in a real installation. This paper is focused on the calculation of g-functions of borehole fields featuring non-vertical boreholes. The strategy utilized consists in representing the boreholes as stacked finite line sources. The temperature along these finite lines, can be calculated by superposition of the effects of each linear heat source in the field. This modeling technique allows to approximate uneven heat distribution among the boreholes and along the axis of each individual borehole. This is a required feature for the calculation of g-functions according to Eskilson's boundary conditions. The test cases presented show that the method yields results that are compatible with the expected physical behavior of the system, and similar to previous results by Eskilson. The computational performance achieved indicates that the method proposed could be potentially utilized during the design phase of these systems. (c) 2015 Elsevier Ltd. All rights reserved.