CO2 applications for enhanced oil recovery and storage in continental shale reservoirs are promising, and there is a need to evaluate the impact of porous structure on oil–water two-phase flow under CO2 environment. In this study, first, digital cores of quartz-rich, carbonate-rich, and clay-rich shales are established using Focused Ion Beam Scanning Electron Microscopy scanning data processed through generative adversarial networks. Subsequently, the pore networks generated by digital cores are quantitatively analyzed using the generalized extreme value distribution. Finally, pore network modeling is carried out to elucidate the effect of porous structural differences on oil–water flow considering CO2 dissolution and capillary forces. The results show that quartz-rich shale, characterized by nanopore intergranular dominance and the highest pore network connectivity, demonstrates the highest relative permeability of the oil phase. Carbonate-rich shale exhibits intermediate relative permeability of oil phase, while clay-rich shale exhibits the worst. The dissolution of CO2 reduces oil–water interfacial tension and oil viscosity, enhances oil mobilization within nanopores, and notably increases the relative permeability of the oil phase. The permeability of the oil phase is governed by pore structure, displaying positive correlations with core heterogeneity, pore radius, coordination number, and throat length, and negative correlations with throat radius.

The Cambrian Qiongzhusi Formation shale of the Sichuan Basin, buried below 4000 m, is a key deep-shale-gas target whose hydraulic fracturing is jointly controlled by bedding weak planes and high in situ stress. This study separates the two controls quantitatively and establishes a reproducible criterion for selecting favourable bedding attitudes and scaling treatment energy. A fluid–solid coupled discrete-element model was built on PFC2D, combining parallel-bond matrix contacts, smooth-joint bedding planes installed on a discrete fracture network, servo-controlled stress boundaries and a cubic-law fluid network, with micromechanical parameters calibrated against a triaxial test on a core from a 4601.76 m depth. Seven bedding inclinations (0–90°) were simulated under a unified stress state, and nine stress combinations were arranged in a 3 × 3 orthogonal design of basal horizontal stress (88.1–125 MPa) and vertical–horizontal stress difference (0.1–6 MPa) at a fixed 60° bedding, then analysed by quadratic response surface and two-way ANOVA. The peak fracturing pressure follows a U-shaped relation with bedding inclination spanning 110–160 MPa. Low-angle bedding (0–30°) favours lateral extension (ratios: 0.63–0.68) but limits vertical connectivity, whereas high-angle bedding (75–90°) promotes vertical breakthrough with lateral reach falling to 0.08–0.09. Ranked by a composite stimulation score built from four normalised sub-scores—initiation efficiency, lateral reach, vertical reach and an SRV proxy—60° is the optimum (S = 0.76), ahead of 15° (0.72) and 45° (0.67), combining the lowest fracturing pressure of 110 MPa with progressive branching rather than sudden breakthrough. Raising basal σh from 88.1 to 125 MPa lengthens propagation approximately fourfold and shifts the failure mode from tensile splitting to shear–tensile ductile rupture, with the tensile-fracture proportion falling from 1.0 to 0.25–0.35. The two-way ANOVA attributes 82–99% of the variance of every response to σh (p ≤ 0.031), whereas the stress difference contributes at most 8.2% and reaches significance for none (p > 0.17) The bedding optimum is placed on a stated objective function rather than qualitative judgement, and the two stress controls are shown to combine additively rather than multiplicatively within the Qiongzhusi stress envelope. Practically, treatment energy should be scaled with basal σ_h rather than with the stress difference, and target intervals selected at 45–60° bedding.