Complexity in shale-gas reservoirs lies in the presence of multiscale networks of pores that vary from nanometer to micrometer scale. Scanning electron microscope (SEM) and atomic force microscope imaging are promising tools for a better understanding of such complex microstructures. Obtaining 3D shale images using focused ion beam-SEM for accurate reservoir forecasting and petrophysical assessment is not, however, currently economically feasible. On the other hand, high-quality 2D shale images are widely available. In this paper, a new method based on higher-order statistics of a porous medium (as opposed to the traditional two-point statistics) is proposed in which a single 2D image of a shale sample is used to reconstruct stochastically equiprobable 3D models of the sample. Because some pores may remain undetected in the SEM images, data from other sources, such as the pore-size distribution obtained from nitrogen adsorption data, are integrated with the overall pore network using an object-based technique. The method benefits from a recent algorithm, the cross- correlation-based simulation, by which high-quality, unconditional/conditional realizations of a given sample porous medium are produced. To improve the ultimate 3D model, a novel iterative algorithm is proposed that refines the quality of the realizations significantly. Furthermore, a new histogram matching, which deals with multimodal continuous properties in shale samples, is also proposed. Finally, quantitative comparison is made by computing various statistical and petrophysical properties for the original samples, as well as the reconstructed model.