Fractures are potential pathways for subsurface fluids. In the context of geologic CO2 sequestration, fractures could threaten the security of reservoirs since they are pathways for leakage. However, the constitutive equations describing the pressure‐saturation and relative permeability‐saturation, key inputs for flow modeling and prediction, are poorly understood for two‐phase fracture flow at the field scale where each fracture might be an element of a complex fracture network or dual‐porosity domain. This knowledge is required for the safe storage of CO2 under potentially fractured caprocks. To address this, we numerically simulated two‐phase viscous flow through horizontal heterogeneous fractures across a broad range of roughness and aperture correlation length. The numerical experiments, which solved the modified local cubic law, showed a weak to strong phase interference during the drainage process; that is, supercritical CO2 displaces brine if the imposed pressure + local viscous force > local capillary force. The drainage process ended when no new brine cells can be further invaded by CO2. Afterward, we froze the saturation field for each individual phase and calculated the relative permeability by single‐phase flow modeling. Thousands of simulated pressure‐saturation and relative permeability curves enabled us to estimate and connect the parameters of the Van‐Genuchten model and of the generalized ν‐type model to fracture roughness and aperture correlation length. This empirical connection will be useful for assessing and predicting immiscible two‐phase flow in horizontal rough fractures at the large scale.