Terahertz (THz) communication is now being considered as one of possible technologies for the sixth generation (6G) wireless communication systems. In this paper, a novel three-dimensional (3D) space-time-frequency non-stationary theoretical channel model is first proposed for 6G THz wireless communication systems employing ultra-massive multiple-input multiple-output (MIMO) technologies with long traveling paths. Considering frequency-dependent diffuse scattering, which is a special property of THz channels different from millimeter wave (mmWave) channels, the relative angles and delays of rays within one cluster will evolve in the frequency domain. Then, a corresponding simulation model is proposed with discrete angles calculated using the method of equal area (MEA). The statistical properties of the proposed theoretical and simulation models are derived and compared, showing good agreements. The accuracy and flexibility of the proposed simulation model are demonstrated by comparing the simulation results of the relative angle spread and root mean square (RMS) delay spread with corresponding measurements.