The ice generation is one of the challenges facing the methane hydrate depressurization, which, however, has not been fully addressed by existing numerical models for hydrate-bearing sediments (HBS). In this study, we develop a high-fidelity, fully coupled thermo-hydro-mechanical-chemical numerical model that incorporates the effect of ice. The model, developed using COMSOL, takes into account water–ice phase change, thermally induced cryogenic suction and constitutive relation in HBS. It is verified well against the temperature, pressure and cumulative gas production of Masuda’s experiment. The model is then employed to investigate multiphysical responses and gas/water production when ice generation is induced by setting a low outlet pressure. The results reveal that ice forms near the outlet boundary of the specimen center, leading to a reduction in intrinsic permeability and fluid velocity and an increase in the bulk modulus of ice-HBS. This enhanced bulk modulus results in higher porosity under axial load. Although the exothermic effect of ice generation promotes the hydrate dissociation, the effect on cumulative gas production is negligible after the ice melts. A negative correlation between ice saturation and water production rate is observed, indicating that a higher gas–water ratio can be achieved by adjusting the ice duration during hydrate production. The developed coupled model proves to be crucial for understanding the effect of ice on hydrate exploitation.