Fractures are widely present in the subsurface, often representing primary channels for fluid flow in low permeability rocks. While fracture surfaces are composed by different minerals and are rough by nature, mathematical models to predict flow properties rarely take in account these heterogeneities. Therefore, the pore-scale mechanisms of flow through fractures are not well understood. Because characterizing multiphase flow phenomena in these geometries has received limited attention, this thesis aims to address this issue, by studying the effect of surface roughness and heterogeneous wettability in immiscible displacement through single fractures. Since analytical solutions are restricted to simple domains and obtaining data from laboratory experiments is unpractical, a 3D direct simulation approach via the lattice Boltzmann method was selected. This was chosen based on its rigorous kinetic derivation, its ability to simulate immiscible displacement, and its versatile boundary conditions. To study the effects of surface heterogeneities, synthetic domains exhibiting geometrical mineral arrangements, and self-affine fractures were created to carry out drainage and imbibition simulations with different input parameters. The relationships of different wetting/non-wetting patterns and surface roughness, with interfacial areas, capillary pressure, and residual fluid saturation were quantified. It has been shown that there is an effective heterogeneous feature size related to the fracture dimensions that modifies the capillary pressure behavior, and the shape of an invasive fluid front. We further found that for increasingly rough surfaces, there is a linear relation between the residual non-wetting saturation and capillary pressure with the aperture distribution. Thus, the shape, mineral size ratio, and surface roughness can have a significant effect on flow behavior. The results of this work can be used to better inform field simulations, by providing physically-accurate input parameters to characterize fracture network models, enhanced flow rate predictions for naturally fractured reservoirs can be obtained.