We present a novel multiscale framework that integrates the single -point multiphase material point method (MPM) and the discrete element method (DEM) to model the complex freeze -thaw behavior of ice -bonded granular media. The proposed numerical framework is featured by (a) employing the continuum -based MPM to solve the macroscopic governing equations for granular systems involving thermo-hydro-mechanical (THM) coupling and phase transitions, and (b) using the grain -scale discontinuum-based DEM to capture the thermodynamically sensitive mechanical behaviors of ice -bonded granular media. The multiscale framework is constructed by attaching a DEM-based representative volume element (RVE) at each material point in MPM. This RVE serves as a live sample of each material point to track the state -dependent effective stress with respect to the local deformation and thermodynamic conditions like ice saturation, bridging the macroscopic phenomena and the underlying microstructural evolution. In particular, we implement a semiimplicit staggered integration scheme for the macroscale THM-coupled MPM to boost computational efficiency and enhance numerical stability. We also propose an innovative ice saturation -dependent bond contact to effectively reproduce the thermodynamically sensitive mechanical behaviors. The new multiscale framework is first benchmarked against analytical solutions for 1D non -isothermal consolidation problems. We then demonstrate its exceptional capability in simulating intricate freeze -thaw behavior of granular media through a boundary value problem involving cyclic freeze -thaw actions. Further cross -scale analyses reveal its potential in capturing key loading- and state -dependent THM responses with explainable microstructural mechanisms during complex freezing and thawing loading conditions.