Dynamic pore-network model (PNM) has been widely used to model pore-scale two-phase flow. Numerical algorithms commonly used for dynamic PNM including IMPES (implicit pressure explicit saturation) and IMP-SIMS (implicit pressure semi-implicit saturation) can be numerically unstable or inaccurate for challenging flow regimes such as low capillary number (Ca) flow and unfavorable displacements. We perform comprehensive analyses of IMPES and IMP-SIMS for a wide range of flow regimes under drainage conditions and develop a novel fully implicit (FI) algorithm to address their limitations. Our simulations show the following: (1) While IMPES was reported to be numerically unstable for low Ca flow, using a smoothed local pore-body capillary pressure curve appears to produce stable simulations. (2) Due to an approximation for the capillary driving force, IMP-SIMS can deviate from quasi-static solutions at equilibrium states especially in heterogeneous networks. (3) Both IMPES and IMP-SIMS introduce mass conservation errors. The errors are small for networks with cubic pore bodies (less than 1.4% for IMPES and 1.2% for IMP-SIMS). They become much greater for networks with square-tube pore bodies (up to 45% for IMPES and 46% for IMP-SIMS). Conversely, the new FI algorithm is numerically stable and mass conservative regardless of the flow regimes and pore geometries. It also precisely recovers the quasi-static solutions at equilibrium states. The FI framework has been extended to include compressible two-phase flow, multicomponent transport, and phase change dynamics. Example simulations of two-phase displacements accounting for phase change show that evaporation and condensation can suppress fingering patterns generated during invasion.