This paper presented a coupled fluid flow-DEM approach to investigate the nucleation, propagation and coalescence of hydraulic fractures in glutenite reservoirs, where embedded gravel causes strong stress concentration at the interfaces between gravel and matrix and significantly affects local mechanical response. The geometry heterogeneity induced by embedded gravel and multi-sized grains was accurately characterized in DEM, while a flow algorithm was used to model fluid flow. Gravel with different sizes and shapes were randomly embedded in specimens, the interaction mechanisms between propagating fractures and encountering gravel were investigated, and the effects of geo-stress difference, gravel strength and grain size distribution were studied. As can be expected, the macroscale fracture morphology is dominated by far-field geo-stress state, but the microscale mechanical responses of glutenite specimens under hydraulic loading are significantly affected by stress heterogeneity induced by the geometry heterogeneity. Due to the stress heterogeneity in the vicinity of a wellbore, multiple initiation points usually form, causing multi-branched fractures, but some of them are soon arrested by encountering gravel due to lack of sufficient energy for further propagation. When propagating fractures encounter gravel, four interaction patterns (attraction, deflection, penetration and termination) are observed, and then high energy is required for hydraulic fractures to propagate forward, this causes local stress reorientation and new initiation points, and finally complicated fracture morphology forms with branched fractures. (C) 2019 Elsevier B.V. All rights reserved.