Ground-source heat pump systems utilizing vertical borehole ground heat exchangers are energy efficient. The most accurate method for sizing the boreholes is based on simulations which make use of superposition of a unit-step thermal response g-function. Due to long calculation times, a library of g-functions was constructed in the 1990's and is still utilized in ground heat exchanger design tools. Literature has been silent concerning the expansion of the library or further improvements to the closed source g-function calculation program.

g-Functions can now be computed by utilizing a recently developed open-source program, pyfunction. Version 1.1 of pygfunction had excessive memory consumption for large borehole fields, motivating development of an alternative program written in C++, named cpgfunction. The program, when coupled with an adaptive discretization scheme, provided the ability to compute a new expanded library of g-functions. A pure C++ fork of cpgfunction, named cpgfunctionEP, has been integrated into the whole building energy simulation program, EnergyPlus, as a third-party g-function calculation tool. Early versions of cpgfunction were faster than pygfunction, though recent developments in pygfunction have broken the computation barrier for g-function calculation and reduced memory consumption. G-functions are computed live time in a newly developed Ground Heat Exchanger Design Toolbox (GHEDT), written in Python. As part of the work described in this thesis, many additional developments in pygfunction have been made.

A fast monthly hybrid time step simulation is developed. This simulation is enhanced and automated compared to an existing program capable of ground heat exchanger sizing; Ground Loop Heat Exchanger Professional (GLHEPRO). The GHEDT sizing results are validated against GLHEPRO. GHEDT provides a novel development of automated selection of a borehole field by placing boreholes within user-defined geometric constraints. The advanced methodology searches borehole fields constructed in domains. When simulated at a constant height over the domain, a discontinuous unimodal excess temperature function is formed. An integer bisection routine enables fast selection of a desired field. The most capable design routine automatically places boreholes within a constrained polygonal available land description excluding a polygonal unavailable drilling zone, where a building or utilities may be located.