Thaw-induced large deformation in granular media has attracted an increasing attention in geotechnical engineering, such as thaw-induced landslides in permafrost due to global warming and submarine landslides due to the dissociation of gas hydrates, which involves complex thermal-mechanical (TM) coupling effects. It remains challenging for conventional continuum-based numerical schemes to capture the TM effect that may cause sudden crash of a granular skeleton with interparticle cementation vanishing. This study presents a hierarchical multiscale modeling framework to model the thaw-induced large deformation in granular media. The proposed framework features a hierarchical conjugate of a continuum theory with a micromechanics-based homogenization technique. Rather than using a contunuum-based constitutive model to capture the thaw-induced response, a micromechanics-based method, i.e., discrete element method (DEM), is employed to model the response of a representative volume element (RVE) that is then affiliated with a material point of the entire continuum domain. Specifically, the mesh-free material point method (MPM) is applied to solving the thaw-induce large-deformation problem at the macroscopic scale, while the DEM-simulated RVE response takes the place of the consititutive model in MPM. To offer speedup for engineering-scale simulations, a GPU-based thread-block-wise parallelization scheme is proposed by leveraging the hierarchical computing structure of the framework. Numerical simulations are showcased to demonstrate the effectiveness and efficiency of the proposed framework.