Cryogenic liquid nitrogen (LN2) fracturing or treatment in nearly impermeable reservoirs has attracted widespread attention due to its strong thermal shock and frost heave effects. We derive a fully coupled thermo-hydro-mechanical (THM) model considering the ice-water phase change to simulate the process of LN2 injection into shale gas reservoirs. The novel model accommodates the following dominant phenomena: (1) LN2 driving unfrozen water flow in porous media; (2) pore water freezing into ice below the freezing temperature; (3) fully coupled thermo-poroelastic effect incorporating frost heave pressure; (4) effective fluid permeability associated with unfrozen water content and phase saturation; (5) temperature and pressure-dependent fluid thermodynamic properties. The proposed model is verified against a classical thermo-poroelastic model and a coupled THM model for freezing rocks. Simulation results show that LN2 injection into water-bearing shale reservoirs will result in ultra-high injection pressure. A high injection rate is beneficial for the fracture initiation, but it is detrimental from phase solidification. Reducing the injection rate is not conducive to the subsequent penetration of LN2 into the reservoir. The evolution of effective permeability, as well as phase saturation, can be reasonably captured during LN2 injection. The comparative advantages between fluid infiltration and temperature propagation are different at high and low injection rates, which leads to significant differences in phase saturation and effective permeability changes. The multiple mechanisms, including pore pressure, frost heave pressure, and thermal stress that induce rock damage, can be efficiently revealed. At relatively low temperatures, the deformation of the reservoir at high and low injection rates is dominated by high-pressure expansion and cryogenic contraction, respectively. The research results can provide some insights for cryogenic non-aqueous fracturing.