This paper introduces a semi-analytical model based on the spectral analysis method for the simulation of transient conductive-convective heat flow in an axisymmetric shallow geothermal system consisting of a double U-tube borehole heat exchanger embedded in a soil mass. The proposed 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. It calculates the temperature distribution in all involved borehole heat exchanger components and the surrounding soil mass using the fast Fourier transform, for the time domain; and the complex Fourier and Fourier-Bessel series, for the spatial domain. Numerical examples illustrating the model capability to reconstruct thermal response test data together with parametric analysis are given. The CPU time for calculating temperature distributions in all involved components, pipe-in, pipe-out, grout, and soil, using 16,384 FFT samples, for the time domain, and 100 Fourier-Bessel series samples, for the spatial domain, was in the order of 3s in a normal PC. The model can be utilized for forward calculations of heat flow in a double U-tube geothermal heat pump system, and can be included in inverse calculations for parameter identification of shallow geothermal systems.