In large scale Ground Source Heat Pump (GSHP) systems, multiple Borehole Heat Exchangers (BHEs) are often connected with the pipe network array to extract shallow geothermal energy. In this study a comprehensive numerical model was developed. The heat transport within and around the BHEs and the pipe network is explicitly quantified in a coupled manner. The model allows a dynamic heat extraction calculation on the individual BHE that is determined by the hydro-thermal processes in the pipe network. The model is thus capable of capturing the long-term thermal interference among BHEs. The model was verified against analytical solution with respect to its hydraulic and thermal balances. Based on it, a series of numerical experiments have been performed to quantitatively investigate the amount of shifted thermal extraction rate in large BHEs array. It is found that, the heat extraction rate on the central BHEs was gradually shifted towards those located at the edge in the long-term operation. Over different seasons, the strongest shifting phenomenon was observed in the month with the lowest thermal load. The shift becomes significant with the increasing number of BHEs installed. The result of numerical study suggests that traditional super-positioned based infinite line source approach with a constant heat flux is not accurate enough for long-term prognosis since it does not fully consider the thermal recharge and the thermal interference effects.