Hydrate saturation ( $S_{\mathrm{h}}$ ) is a major factor governing the creep behavior of hydrate-bearing sediments (HBSs). During long-term natural gas hydrate (NGH) extraction, creep deformation in HBSs may induce significant geological and environmental hazards, making accurate constitutive modeling essential. In this paper, a creep model is established by integrating a variable-order fractal derivative with elastic–viscoplastic theory, with the objective of establishing linear quantitative relationships between model parameters and hydrate saturation. The proposed model is validated against available experimental data and demonstrates strong capability in reproducing the nonlinear time-dependent creep behavior of HBSs. The model is theoretically capable of predicting rheological responses under arbitrary loading paths and varying saturation conditions. A systematic range analysis is further conducted to evaluate the sensitivity of individual parameters and their interactions on the creep response. The results provide clear mechanistic insights into how hydrate saturation controls the time-dependent deformation characteristics of hydrate-bearing sediments. Overall, this investigation offers a robust theoretical framework for improved assessment of reservoir deformation, wellbore stability, and geo-environmental risks associated with NGH exploitation.
