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 semi-implicit 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.