Hydrate　distribution　heterogeneity　is　often　observed　in　natural　and　artificial　hydrate-bearing　sediments　(HBSs).　To　capture　hydrate　distribution　heterogeneity,　a　pore-scale　method　is　proposed　to　model　cemented　HBSs　using　the　discrete　element　method　(DEM).　This　method　takes　into　account　the　quantitative　effect　of　hydrate　saturation　in　a　sediment　pore　on　the　contact　bond　parameters　surrounding　the　pore.　A　series　of　DEM　specimens　with　different　macroscopically　and　microscopically　heterogeneous　hydrate　distributions　are　prepared.　The　mechanical　behavior　of　heterogeneous　HBSs　is　investigated　by　performing　biaxial　compression　tests　with　flexible　boundaries.　The　simulation　results　show　that　both　macroscopic　and　microscopic　hydrate　distribution　heterogeneity　can　influence　the　mechanical　properties　of　HBSs.　The　shear　strength　is　promoted　in　both　macroscopically　and　microscopically　heterogeneous　HBSs.　Longitudinally　heterogeneous　HBSs　have　a　higher　secant　modulus,　while　transversely　heterogeneous　HBSs　have　a　lower　secant　modulus　than　homogeneous　HBSs.　The　secant　modulus　of　microscopically　heterogeneous　HBSs　first　increases　and　then　decreases　with　increasing　pore　hydrate　saturation.　It　is　found　that　the　deformation　behavior　and　bond　breakage　evolution　of　HBSs　depend　on　hydrate　distribution　heterogeneity.　These　findings　can　provide　insights　into　understanding　the　mechanical　behavior　of　natural　HBSs　with　heterogeneous　hydrate　distributions.