Gas production from hydrates induced by depressurization is a complex thermal-hydrodynamic-mechanical–chemical (THMC) coupled process. In this paper, we present a THMC coupled model to simulate the fluid flow in hydrate-bearing sediments (HBS) and the geomechanical behavior of HBS. The model is made of two subsystems, which are the fluid part of non-isothermal multi-phase flow with hydrate kinetic and solid part of geomechanical deformation. It accounts for two-way coupling effects between these two subsystems, i.e. the effect of pore pressure and hydrate dissociation on the solid mechanical behavior and the effect of stress on the hydraulic behavior. A new numerical method based on the hybrid control volume finite element method (CVFEM)-finite element method (FEM) is developed to solve the mathematical models. The local conservative CVFEM is used for the fluid part, and the standard FEM for the solid part. In the framework of hybrid CVFEM-FEM, the local conservation is reserved and the primary variables for the two subsystem are co-located. A multi-point flux approximation (MPFA) is adopted without orthogonal meshes so that it is very flexible to build complex geometrical models. The accuracy and reliability of the newly developed simulator QIMGHyd-THMC are tested by comparing with two experimental examples and a large-scale benchmark problem of other popular simulators.