The effect on numerical solution of different thermal boundary conditions at the ground surface was analysed in modelling HGHEs. Boundary conditions of the 1st, 2nd and 3rd kind have been alternately tested by means of a finite element numerical code, solving the unsteady-state heat transfer problem in a 2D domain. An energy balance equation at the ground surface (3rd kind BC) has been developed and implemented in the numerical model. A preliminary simulation has been carried out in absence of the HGHE operating using real weather data. The solution has been validated with experimental data, and assumed as reference. The calibrated GSEB equation proved to properly predict the temperature in the soil. The resulting heat flux and temperature at the top of the domain have been considered respectively as the 2nd and 1st kind of equivalent boundary conditions for two new models. Finally, all three models have been solved with the supposed HGHE operating, to analyse how the different BCs affected the numerical solution. The results have been compared in terms of average temperature at the HGHE wall surface and in the ground. The use of a heat flux as BCs at the ground surface appeared as an extremely precautionary approach due to the resulting thermal drift in the soil. On the contrary, to assign an energy balance equation or a temperature as BCs on the ground surface seemed to have a limited effect in terms of temperature at the heat exchanger and in the soil.