Abstract: The glutenite reservoir is strongly heterogeneous due to the random distribution of gravels, making it challenging to perform hydraulic fracturing effectively. To solve this issue, it is essential to study interaction behavior between hydraulic fractures (HFs) and gravels. A coupled hydro-mechanical model is proposed for HF propagation in glutenite using a grain-based discrete element method. This paper first investigates the dynamic evolution of HFs in glutenite, then analyzes the influences of various factors such as horizontal stress difference (Δσ), minimum horizontal stress (σh), gravel content (Vg), gravel size (dg), and stiffness ratio of gravel to matrix (Rs) on HF propagation geometries. Results show that penetrating the gravel is the primary HF-gravel interaction behavior, which follows sequential and staggered initiation modes. Bypassing the gravel is the secondary behavior, which obeys the sequential initiation mode and occurs when the orientation of the gravel boundary is inclined to the maximum horizontal stress (σH). An offset along the gravel boundary is usually formed while penetrating gravels, and the offsets may cause fracture widths to decrease by 37.8%–84.4%. Even if stress dominates the direction of HF propagation, HFs still tend to deflect within gravels. The deviation angle from σH decreases with rising Δσ and increases with the increase of dg and Rs. Additionally, intra-gravel shear HFs (IGS-HFs) are prone to be generated in coarse-grained glutenite under high Δσ, while more gravel-bypassing shear HFs (GBS-HFs) tend to be created in argillaceous glutenite with high Rs than in sandy glutenite with low Rs. The findings above prompt the emergence of a novel HF propagation pattern in glutenite, which helps to understand the real HF geometries and to provide theoretical guidance for treatments in the field.