Terahertz (THz) communication is now being considered as one of possible technologies for the sixth generation (6G) communication systems. In this paper, a novel three-dimensional (3D) space-time-frequency non-stationary massive multiple-input multiple-output (MIMO) channel model for 6G THz indoor communication systems is proposed. In this geometry-based stochastic model (GBSM), the initialization and evolution of parameters in time, space, and frequency domains are developed to generate the complete channel transfer function (CTF). Based on the proposed model, the correlation functions including time auto-correlation function (ACF), spatial crosscorrelation function (CCF), and frequency correlation function (FCF) are investigated. The results show that the statistical properties of the simulation model match well with those of the theoretical model. The stationary intervals at different frequencies are simulated. The non-stationarity in time, space, and frequency domains is verified by theoretical derivations and simulations.