We describe a two­phase pore network simulator of drainage and imbibition, which integrates a realistic representation of pore connectivity and morphology, a quasi­static descrip­tion of uid displacement mechanisms, and a sound representation of the wetting properties of a sedimentary rock and of their alteration. The simulator works with 3D disordered net­ works of cylindrical ducts with triangular, square and circular cross­sections obtained directly from the analysis of micro­CT images of rock samples. All pore­level displacement mecha­nisms: piston­type, snap­off, and cooperative pore body lling are considered with arbitrary receding and advancing contact angles. Bond invasion percolation description is used in pri­mary drainage, while bond­site invasion percolation with ordinary percolation on a dual net­ work and compact cluster growth are used in secondary imbibition. In the paper, we resolve how to calculate the relative permeability of NAPL in the quasi­static approximation of im­bibition, and illustrate spatial distribution of the clusters of trapped NAPL using our gener­alization to disordered networks of the Hoshen­Kopelman cluster­labelling algorithm. To un­derstand the impact of wettability alteration on the capillary pressure and relative permeabil­ity functions, we perform a series of drainage and imbibition simulations by changing the range of advancing contact angles. Our study indicates that in imbibition, transport properties of a permeable solid are quite sensitive to the hysteresis between the receding and advancing con­tact angle. This sensitivity reects competition among the different displacement mechanisms, which shapes the relative permeabilities, capillary pressures, and the distribution of the trapped NAPL.