A comprehensive spectral element formulation for nonhomogeneous heat flow in a shallow geothermal system consisting of a borehole heat exchanger embedded in a multilayer soil mass is introduced. The spectral element method is utilized to solve the governing heat equations in the borehole heat exchanger and the soil mass simultaneously using the fast Fourier transform, the eigenfunction expansion, the Fourier Bessel series and the complex Fourier series, together with the finite element method. Only one spectral element is necessary to describe heat flow in a homogeneous domain. For a nonhomogeneous multilayer system, the number of spectral elements is equal to the number of layers. The proposed spectral element model combines the exactness of the analytical methods with an important extent of generality in describing the geometry and boundary conditions of the numerical methods. Verification examples illustrating the model accuracy, and numerical examples illustrating its capability to simulate multilayer systems are given. Despite the apparent rigor of the proposed model, it is robust, computationally efficient and easy to implement in computer codes. (C) 2016 Elsevier Ltd. All rights reserved.