With the increasing engineering applications of geothermal borehole heat exchangers (BHEs), accurate and reliable mathematical models can help advance their thermal design and operations. In this study, an analytical model with a time-dependent heat flux boundary condition on the borehole wall is developed, capable of predicting the thermal performance of single, double, and multiple closed-loop BHEs, with an emphasis on solar- and waste-heat-assisted geothermal borehole systems (S-GBS and W-GBS) for energy storage. This analytical framework begins with a one-dimensional transient heat conduction problem subjected to a time-dependent heat flux for a single borehole. The single borehole scenario is then extended to multiple boreholes by exploiting lines of symmetry (or thermal superposition). A final expression of the temperature distribution along the center line is attained for single, double, and multiple boreholes, which is verified with a two-dimensional finite-element numerical model (less than 0.7% mean absolute deviation) and uses much lesser computational power and time. The analytical solution is also validated against a field-scale experiment from the literature regarding the borehole and ground temperatures at different time frames, with an absolute error below 6.3%. Further, the thermal performance of S-GBS and W-GBS is compared for a 3-by-3 borehole configuration using the analytical model to ensure its versatility in thermal energy storage. It is concluded that our proposed analytical framework can rapidly evaluate closed-loop geothermal BHEs, regardless of the numbers of boreholes and the type of the heat flux on the borehole wall.