With the advancement of computational technology, many BHE models were conducted by numerical methods, particularly using the computational fluid dynamics (CFD) approach. Complex tube design [19] and the influence of groundwater flow [20] could be considered in CFD models. In addition, the thermal and hydraulic performance of BHEs in different diameter ratios could be evaluated numerically, which in turn facilitated the optimal design of BHEs [21]. [PERSON] et al. [22] also employed a parameter optimization approach using a full three-dimensional CFD model to find the optimum heat exchanger for a given system. Apart from the CFD approach, the finite-element method (FEM) had also been used on GSHP systems. A well-known computationally efficient FEM was formulated in conjunction with the Petrov-Galerkin and finite difference methods for both steady-state [23] and transient heat flow problems [24]. Short-term and long-term thermal performances in BHEs were recently evaluated by a 1-D FEM framework with the internal capacity of the borehole and the thermal resistance between the two legs of U-pipe [25]. To minimize the computational costs in 3-D numerical models, [PERSON] et al. [26] tackled such heat flow problems by firstly solving a 1-D domain and then extending it to a 3-D configuration via FEM. The second methodology is the hybrid or semi-analytical approach, which removes some restrictions in the analytical method and reduces the computational power and time in the numerical model. \(g\)-Functions, introduced by [PERSON] [27], offered a non-dimensional relation between the heat extracted from the ground at the borehole wall and the borehole wall temperature, thus making it fast-to-compute and practical in GSHP analyses. A few new contributions to the original \(g\)-function were made to adapt short-time responses [28], various boundaries [29; 30], and the effect of natural convection in the groundwater [31]. Another type of hybrid model couples analytical solution with numerical algorithm. Recently, [PERSON] and [PERSON] [32] developed a relatively cost-effective hybrid method for vertical BHEs along with groundwater flow, where the numerical results from a co-located FEM and finite volume method were marched with thermal resistance models. The hybrid approach can also march two spatial dimensions (so-called \"1+1D algorithm\") to greatly increase the computational efficiency compared with solving a full-scale model. This method has been employed in solar-assisted geothermal borehole systems for evaluating conjugate heat transfer [33]. While the numerical schemes in the numerical and hybrid models often require high computational power and time, the analytical solution is generally easy-to-compute and the developed closed-form final expression could be calculated instantaneously, thus gaining a significant amount of practical interest. Most analytical models assumed an infinite/finite line source or infinite cylindrical source as the borehole, and then exact solutions could be found in a transient cylindrical heat conduction problem [34]. Another type of analytical model took the assumptions of coaxial borehole heat exchangers, where the computational domain was considered as an annular space. Several scenarios subjected to heat flux boundaries were resolved analytically by [PERSON] and [PERSON] [35; 36]. As for multiple boreholes, these analytical solutions for a single borehole could be added into a multiple boreholes configuration [37; 38]. Recently, [PERSON] et al. [39] solved a one-dimensional transient heat equation in the cylindrical coordinate subjected to a constant heat flux at the borehole and extended the solution from a single borehole to N-by-N boreholes by the method of superposition (i.e., the use of symmetry). This paper fills the research gaps by extending [PERSON] et al.