In this study, we determine the fast polarization direction (FPD) and delay time (\(\delta\)t) of shear wave splitting from 313 Bhutan Himalayan crustal earthquakes recorded by 44 stations of the GANSSER (Geodynamics ANd Seismic Structure of the Eastern Himalaya Region) network deployed and operational in Bhutan from January 2013 to December 2014. We acquired 1686 high-quality measurements by using the Multiple Filter Automatic Splitting Technique (MFAST), which enables rapid processing of millions of events. The delay-time was found to vary between 0.125-0.6 s, indicating the presence of strong crustal anisotropy with an increase in \(\delta\)t with depth. We also observe that the crustal anisotropy varies from \(\sim\)3% at 0-20 km depth to \(\sim\)5% in the 20-35 km depth interval, and more than 8% below 35 km hypocentral depth. The fast directions across the region is generally aligning with the regional stress directions, which is significantly different from the predominant trend of SW-NE to WSW-ENE as observed from the upper mantle splitting. An increase in the number of anisotropic measurements in the western and central Bhutan may be linked to Indian plate subduction. Bhutan Himalaya displays a noticeable along-strike variation in crustal thickness and geometry of Main Himalayan Thrust across a short lateral distance (\(\sim\)90 km), which is significantly reflected in the anisotropic patterns, due to the uneven advancement of the Indian plate below western and eastern Bhutan. Integrated patterns of crustal and mantle splitting parameters provide new insights into the deformation processes beneath the Bhutan Himalaya and their implications in regional geodynamics.
