The geometry of hydraulic fractures in deep shale facies is significantly affected by the longitudinal inhomogeneity of rock physical properties and stresses. Numerous studies have been conducted on the influence of the longitudinal inhomogeneity of rocks on fracture morphology. However, there is still a lack of research that simultaneously considers the reservoir dip, bedding plane interface, and longitudinal inhomogeneity of the reservoir. To fill this gap, a three-dimensional (3D) numerical model of multireservoir hydraulic fracturing, which takes into account the bedding plane interface, was developed using the finite element method (FEM). The Drucker−Prager elastic-plasticity criterion was incorporated to accurately represent the plasticity of deep shale. The research revealed the influence of the formation dip angle on fracture morphology. Additionally, the perforation layer position and pump rate were optimized based on the actual geological parameters in North Jiangsu shale reservoir. The study findings indicate that reservoir fractures with a formation dip are easily detected by the interface. However, it is not necessarily true that the larger the formation dip, the easier it is for fluids to enter the interface. Fracturing from high-strength and stress reservoirs to lower reservoirs promotes the propagation of fracture height and the connectivity of multiple reservoirs. On the other hand, fractures initiated from low-strength and stress reservoirs tend to be confined to adjacent reservoirs more easily. The pump rate significantly affects the vertical propagation of fractures. At high interface strength, fractures with pump rate below 2.4 m 3 /min can only propagate at the perforation layer. The limited fracture height in shale reservoirs is likely due to substantial energy consumption by the fracturing fluid at the bedding plane interface. These studies offer theoretical guidance for understanding the vertical propagation of fractures in a deep multilayer reservoir.