A new pore network extraction scheme is presented that uses geometrical, topological, and mineralogical features of natural porous media to construct digital, faithful representations of the pore space. This newly-developed methodology allows users to construct disordered three-dimensional pore network systems that can be used to model physically-based displacement mechanisms at the pore scale in natural porous media. This work has three main components: (1) digital image analysis, (2) equivalent pore network extraction, and (3) pore network modeling of single- and two-phase flow processes. The pore network extraction approach reads as input high-resolution images of pore space in a given porous medium. It uses advanced image processing techniques to identify and quantify important features such as global heterogeneity and local anisotropy to construct representative digital replicates of the medium in three dimensions. To this end, an automated computational method is introduced to replace subjective manual image analyses thereby mitigating uncertainty. Modern computational algorithms are incorporated to optimize data management. This allows us to scale-up the implemented algorithms to construct large, core-scale pore networks. These are digital porous media with dimensions identical to those of the rock samples used in physical flow experiments. The extracted pore networks are subsequently used to model immiscible two-phase flow in different sandstone and carbonate rock samples and to predict relative permeabilities and residual saturation. To handle the numerical analysis with this level of computational complexity effectively, we employ: (1) customized data structures for optimal data management, (2) network analysis for optimal multidimensional graph analysis, (3) high-performance computing techniques for maximizing the cache efficiency, and (4) data analytics and Big Data analogies to address simulation bottlenecks by strategically collecting, sorting, and mapping information.