This　paper　presents　a　meso-scale　numerical　modelling　method　to　investigate　the　stability　of　tunnels　in　sandy　cobble　stratum.　Using　Monte　Carlo　simulation,　Python　routine　is　implemented　to　randomly　generate　and　filter　the　location　information　of　blocks　including　center　coordinates,　sizes　and　orientation　angles,　then　extract　and　import　them　into　CAD　modelling　software　in　batches　to　mesh.　The　problem　of　overlapping　detection　for　elliptical　and　ellipsoidal　blocks　is　solved　based　on　affine　transformation,　achieving　higher　efficiency　and　applicability　than　the　existing　methods.　Then,　using　the　proposed　numerical　modelling　method,　convergence-confinement　analysis　of　circular　tunnels　with　different　cover　depths　is　performed　in　plane　strain　condition.　Meanwhile,　face　stability　analysis　of　shield　tunnels　is　carried　out　in　three-dimensional　condition.　Results　show　that　sandy　cobble　stratum　may　not　be　in　accordance　with　the　scale-independent　soil-rock　mixture.　Blocks　with　small　sizes　are　geotechnically　significant　and　contribute　much　to　the　overall　strength.　When　the　cover　depth　ratio　reaches　a　certain　level,　limit　support　pressure　remains　basically　unchanged　due　to　the　strong　arching　effect　above　tunnel　crown,　while　permissible　support　pressure　increases　approximately　linearly　to　control　deformation.　Compared　with　the　multi　scale　model,　hardening　soil　model　considering　the　overall　strength　is　suitable　for　the　stability　analysis　of　tunnels　in　sandy　cobble　stratum　at　the　macro-scale　with　high　precision,　while　the　macro　homogeneous　model　ignoring　the　presence　of　blocks　provides　a　conservative　(i.e.　uneconomical)　prediction　for　the　design　of　tunnel　stability.