In this article, we consider a geothermal energy storage system in which the spatio-temporal temperature distribution is modeled by a heat equation with a time-dependent convection term. Such storage systems are often embedded in residential heating systems. The control and management of such systems requires knowledge of aggregated characteristics of the temperature distribution in the storage. These describe the input–output behavior of the storage, the associated energy flows, and their response to charging and discharging processes. Our aim is to derive an efficient, approximate description of these characteristics by using low-dimensional systems of ordinary differential equations (ODEs). This leads to a model order reduction problem for a large-scale linear system of ODEs resulting from the semidiscretization of the heat equation combined with a linear algebraic output equation. In a first step, we approximated the nonautonomous system of ODEs by a linear time-invariant system. Then, we applied Lyapunov balanced truncation model order reduction to approximate the output by a reduced-order system that has only a few state equations but almost the same input–output behavior. The results of our extensive numerical experiments show the efficiency of the applied model order reduction method. We found that only a few suitably chosen ODEs are sufficient to achieve good approximations of the input–output behavior of the storage.