Oil reservoirs have structural heterogeneities across multiple length scales and, particularly in carbonates, complexly distributed wettabilities. The interplay of structural and wettability heterogeneities is the fundamental control for sweep efficiency and oil recovery. This interplay must be captured in physically robust flow functions, such as relative permeability and capillary pressure functions. Such flow functions then allow us to choose the best improved-oil-recovery (IOR) or enhanced-oil-recovery (EOR) process and forecast oil recovery with adequate precision. Obtaining flow functions for reservoir rocks with varying wettability is a challenging task, especially when three fluid phases coexist. In this work, we use pore-network modeling, a reliable and physically based simulation tool, to predict three-phase flow functions. We have developed a new pore-scale network model for rocks with variable wettability. Unlike other models, this model combines three new and important features. (1) Our network model comprises a novel thermodynamic criterion for the formation and collapse of oil layers. This captures film/layer flow of oil adequately, which affects the oil relative permeability at low oil saturation. We can therefore predict residual oil more accurately. (2) We implemented multiple displacement chains, in which injection of one phase at the inlet triggers a chain of interface displacements throughout the network. This allows us to accurately model the mobilization of disconnected phase clusters that arise during higher-order [water-alternating-gas (WAG)] floods. Again, this feature is key to a better prediction of residual oil saturation (ROS). (3) Our model takes realistic 3D pore networks extracted from pore-space reconstruction methods and X-ray computerizedtomography (CT) images as input. This preserves both topology and pore shape of the rock, providing better estimates of phase conductivities and relative permeability. We have validated our model by use of available experimental data for a range of wettabilities and demonstrated the impact of single vs. multiple displacement on residual oil. We also used a proof-of concept study in which we use flow functions for different wettabilities that have been computed with our model in field-scale reservoir simulations to forecast oil recovery during tertiary gas injection. These results are compared with predictions that used empirical flow functions. Flow functions computed by our network model gave higher oil recovery than corresponding flow functions calculated by empirical models; oil recovery increases with decreasing water-wetness. This shows that the pore-scale physics encapsulated in our new network model leads to the right emergent behavior at the reservoir scale.