Accurate characterization of the threshold pressure gradient (TPG) in shale oil reservoirs is essential for oil production. However, oil−water two-phase flow through shale oil reservoirs with narrow pores involves thermal−hydrological−mechanical coupling, and the relationship between TPG and pore structures under multiple mechanisms has not been adequately modeled mathematically. In this paper, an analytical TPG model of shale oil reservoirs is derived first by considering the combined effects of effective stress, temperature, capillary pressure, pore structures, and boundary layer. After the established model is adequately verified by the obtainable experimental data, it is embedded into a multistage fractured horizontal well productivity model to investigate the influence of relevant parameters on shale oil production. The simulation results suggest that the TPG for the oil phase increases with the increase in effective stress, leading to the decrease in shale oil production. The cumulative oil production with dynamic TPG is much greater than that considering fixed TPG. Moreover, the wettability of the oil−water two phase (e.g., oil−water interfacial tension and equilibrium contact angle) has crucial influences on well performance; thus, changing the wettability of the oil−water two phase by chemical injection is an effective method to enhance oil recovery. The research results can provide technical support for the efficient development of shale oil reservoirs from the perspective of pore-scale two-phase nonlinear flow characteristics and macroscopic productivity evaluation.