Pore-network modeling is a widely used predictive tool for pore-scale studies in various applications that deal with multiphase flow in heterogeneous natural rocks. Despite recent improvements to enable pore-network modeling on simplified pore geometry extracted from rock core images and and its computational efficiency compared to direct numerical simulation methods, there are still limitations to modeling a large representative pore-network for heterogeneous cores. These are due to the technical limits on sample size to discern void space during X-ray scanning and computational limits on pore-network extraction algorithms. Thus, there is a need for pore-scale modeling approaches that have the natural advantages of pore-network modeling and can overcome these limitations, thereby enabling better representation of heterogeneity of the 3D complex pore structure and enhancing the accuracy of prediction of macroscopic properties. This paper addresses these issues with a workflow that includes a novel pore-network stitching method to provide large-enough representative pore-network for a core. This workflow uses micro-CT images of heterogeneous reservoir rock cores at different resolutions to characterize the pore structure in order to select few signature parts of the core and extract their equivalent pore-network models. The space between these signature pore-networks is filled by using their statistics to generate realizations of pore-networks which are then connected together using a deterministic layered stitching method. The output of this workflow is a large pore-network that can be used in any flow and transport solver. We validate all steps of this method on different types of natural rocks based on single-phase and two-phase flow properties such as drainage relative permeability curves of carbon dioxide and brine flow. Then, we apply the stochastic workflow on two large domain problems, connecting distant pore-networks and modeling a heterogeneous core. We generate multiple realizations and compare the average results with properties from a defined reference pore-network for each problem. We demonstrate that signature parts of a heterogeneous core, which are a small portion of its entire volume, are sufficient inputs for the developed pore-network stitching method to construct a representative pore-network and predict flow properties.