Terahertz (THz) extremely large-scale MIMO (XL-MIMO) is considered a key enabling technology for 6G and beyond due to its advantages such as wide bandwidth and high beam gain. As the frequency and array size increase, users are more likely to fall within the near-field (NF) region, where the far-field plane-wave assumption no longer holds. This also introduces spatial non-stationarity (SnS), as different antenna elements observe distinct multipath characteristics. Therefore, this paper proposes a THz XL-MIMO channel model that accounts for both NF propagation and SnS, validated using channel measurement data. In this work, we first conduct THz XL-MIMO channel measurements at 100 GHz and 132 GHz using 301- and 531-element ULAs in indoor environments, revealing pronounced NF effects characterized by nonlinear inter-element phase variations, as well as element-dependent delay and angle shifts. Moreover, the SnS phenomenon is observed, arising not only from blockage but also from inconsistent reflection or scattering. Based on these observations, a hybrid NF channel modeling approach combining the scatterer-excited point-source model and the specular reflection model is proposed to capture nonlinear phase variation. For SnS modeling, amplitude attenuation factors (AAFs) are introduced to characterize the continuous variation of path power across the array. By analyzing the statistical distribution and spatial autocorrelation properties of AAFs, a statistical rank-matching-based method is proposed for their generation. Finally, the model is validated using measured data. Evaluation across metrics such as entropy capacity, condition number, spatial correlation, channel gain, Rician K-factor, and RMS delay spread confirms that the proposed model closely aligns with measurements and effectively characterizes the essential features of THz XL-MIMO channels.