Abstract T1-weighted imaging holds wide applications in clinical and research settings; however, the challenge of inter-scanner variability arises when combining data across scanners, which impedes multi-site research. To address this, post-acquisition harmonization methods such as statistical or deep learning approaches have been proposed to unify cross-scanner images. Nevertheless, how inter-scanner variability manifests in images and derived measures, and how to harmonize it in an interpretable manner, remains underexplored. To broaden our knowledge of inter-scanner variability and leverage it to develop a new harmonization strategy, we devised a pipeline to assess the interpretable inter-scanner variability in matched T1-weighted images across four 3T MRI scanners. The pipeline incorporates ComBat modeling with 3D superpixel parcellation algorithm (namely SP-ComBat), which estimates location and scale effects to quantify the shift and spread in relative signal distributions, respectively, concerning brain tissues in the image domain. The estimated parametric maps revealed significant contrast deviations compared to the joint signal distribution across scanners (p < 0.001), and the identified deviations in signal intensities may relate to differences in the inversion time acquisition parameter. To reduce the inter-scanner variability, we implemented a harmonization strategy involving proper image preprocessing and site effect removal by ComBat-derived parameters, achieving substantial improvement in image quality and significant reduction in variation of volumetric measures of brain tissues (p < 0.001). We also applied SP-ComBat to evaluate and characterize the performance of various image harmonization techniques, demonstrating a new way to assess image harmonization. In addition, we reported various metrics of T1-weighted images to quantify the impact of inter-scanner variation, including signal-to-noise ratio, contrast-to-noise ratio, signal inhomogeneity index, and structural similarity index. This study demonstrates a pipeline that extends the implementation of statistical ComBat method to the image domain in a practical manner for characterizing and harmonizing the inter-scanner variability in T1-weighted images, providing further insight for the studies focusing on the development of image harmonization methodologies and their applications.