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 description of fluid displacement mechanisms, and a sound representation of the wetting properties of a sedimentary rock and of their alteration. The simulator works with 3-D disordered networks of cylindrical ducts with triangular, square, and circular cross sections obtained directly from the analysis of microfocused computed tomography (CT) images of rock samples. All pore-level displacement mechanisms (piston type, snap off, and cooperative pore body filling) are considered with arbitrary receding and advancing contact angles. Bond invasion percolation description is used in primary drainage, while bond site invasion percolation with ordinary percolation on a dual network and compact cluster growth is used in secondary imbibition. In this paper, we resolve how to calculate the relative permeability of nonaqueous phase liquid (NAPL) in the quasi-static approximation of imbibition and illustrate spatial distribution of the clusters of trapped NAPL using our generalization to disordered networks of the Hoshen-Kopelman cluster-labeling algorithm. To understand the impact of wettability alteration on the capillary pressure and relative permeability 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 contact angle. This sensitivity reflects competition among the different displacement mechanisms, which shapes the relative permeabilities, capillary pressures, and the distribution of the trapped NAPL.