Nowadays, utilizing shallow geothermal energy for heating and cooling buildings has received increased interest in the building sector. Among different technologies, large borehole heat exchanger arrays are widely employed to supply heat to various types of buildings. Recently, a 16-borehole array was constructed to extract shallow geothermal energy to provide heat to a newly-developed public building in Berlin. To guarantee the quality of the numerical model and reveal its sensitivity to different subsurface conditions, model simulations were conducted for 25 years with two finite element simulators, namely the open-source code OpenGeoSys and the widely applied commercial software FEFLOW. Given proper numerical settings, the simulation results from OpenGeoSys and FEFLOW are in good agreement. However, further analysis reveals differences with respect to borehole inflow temperature calculation implemented in the two software. It is found that FEFLOW intrinsically uses the outflow temperature from the previous time step to determine the current inflow temperature, which makes it capable of much faster simulation by avoiding iterations within a single time step. In comparison, OpenGeoSys always updates the inflow and outflow temperature based on their current time step values. Because the updates are performed after each iteration, it delivers more accurate results with the expense of longer simulation time. Based on this case study, OpenGeoSys is a valid alternative to FEFLOW for modeling ground source heat pump systems. For modellers working in this field, it is thus recommended to adopt small enough time step size, so that potential numerical error can be avoided.