A borehole heat exchangers (BHEs) combined with pumping-injection well is established in areas where groundwater is shallow but the seepage velocity is weak, which sets up pumping and injection wells on both sides of the BHEs. According to the three-dimensional unsteady state heat transfer model in aquifer, we derive the convection-dispersion analytical solution of excess temperature in aquifer that considers groundwater forced seepage and thermal dispersion effects in aquifer and the axial effect of the BHEs. Then, we use dimensional analysis method and similarity criteria to build a controllable forced seepage sandbox. The theoretical analysis is combined with the indoor experiment test to verify the correctness and accuracy of the numerical simulation software FEFLOW7.1. On this basis, we perform the numerical simulation calculation to explore the effects of different pumping-injection flow volume on the Darcy flow velocity of the aquifer where the BHEs are located, the average heat transfer efficiency and the heat transfer rates per unit borehole depth of the BHEs. The results show that when the pumping flow volume is increased from 200 m3∙d-1 to 1200 m3∙d-1, the Darcy velocity correspondingly increases to about 10 times the previous velocity. The average heat efficiency coefficient of the BHEs increases by 11.5% in cooling stage, and by 7.5% in heating stage. When the pumping-injection flow volume is 400~600 m3∙d-1, the increment of heat transfer rates per unit borehole depth of the BHEs reaches 12.8~17.9 W∙m-1 and 3.6~4.2 W∙m-1 during the cooling stage and heating stage respectively, and then decreases as the flow volume increases gradually.