Abstract: Exploring the interaction between hydraulic fractures and complex geological conditions is critical for multilayered commingling production in the laminated continental shale oil reservoirs. In this paper, a 2D hydro-mechanical coupling numerical model is developed to investigate the fracture propagation behavior affected by the joint interaction of multi-lithologic stack and shale anisotropy. The model adopts a smear approach to reproduce the mechanical anisotropy of shale observed from laboratory experiments with powerful finite-discrete element method (FDEM) to precisely capture the transition from elastic deformation to failure during fluid injection in layered heterogeneous media. The results indicate that shale anisotropy affects hydraulic fracture initiation and propagation behavior. The estimated breakdown pressure is 15% higher than that in horizontal homogeneous shale oil reservoirs. The elastic anisotropy alters the stress trajectory and magnitudes, while the strength anisotropy affects the failure mode and fracture morphology. Under the combined two factors, the established fracture network reveals potent cross-layer abilities with less activation of weak planes. Additionally, the sedimentary structure of thin interlayers hinders fracture height extension, resulting in a limited stimulated reservoir volume (SRV). Optimization of engineering and geological parameters could mitigate this limitation and efficiently co-develop the multiple sweet-spot pay zones. For field application, it is proposed to select a modest stress difference formation (Kv around 0.75–1.00) for stimulation. Then, an alternated high/low injection rate can be employed to improve the cross-layer ability and activate the underlying weak planes, finally realizing an ideal SRV. The key findings are expected to provide new insights into the fracture propagation mechanism and guide reservoir stimulation in continental shale oil.