Abstract — Understanding the physics of two-phase flow of CO2 and brine in porous geological formations is essential to sequestration of carbon dioxide in deep saline reservoirs, as well as the older problem of enhanced oil recovery from hydrocarbon reservoirs by CO2 injection. A pilot CO2 injection in Decatur, Illinois, was undertaken, with the injection zone being the highly saline and heterogeneous Mt. Simon sandstone, in order to better understand the feasibility of full-scale sequestration process. This paper reports the results of an extensive study of the morphology of the sandstone and its heterogeneity, and simulation of single-phase and two-phase flow of CO2 and brine in the formation's three-dimensional images. As we demonstrate by extensive analysis, the formation is much more heterogeneous than the typical sandstone, such as Berea sandstone. In addition to characterizing the morphology of the sandstone and computing its important flow characteristics, an important goal of the study is to compare the accuracy and computational efficiency of three distinct simulation approaches, namely, the lattice-Boltzmann (LB) approach, direct numerical simulation (DNS) of the governing equations of fluid flow that uses the finite-volume method coupled with the OpenFOAM simulator, and pore-network (PN) simulation. After validating the simulators by comparing the computed relative permeabilities that they produce for Berea sandstone, we simulate displacement of brine by CO2 at low and relatively high capillary numbers, and compute the relative permeabilities and other quantities of interest. We demonstrate that all the three methods provide consistent relative permeability-saturation functions that are in close agreement with one another. However, although the LB and DNS both produce similar relative permeabilities, the DNS approach is computationally more efficient because it simulates drainage by only a single set of computations over the entire saturation range, whereas the LB simulation requires separate simulation for each set saturation. Thus, the question of what method to use for simulating such flow processes at the scale of core plugs should mainly be addressed based on the computational time that one can afford and the computational resources that one has access to. Another important question addressed is the effect of the resolution of the computational grids or lattices used, particularly when one uses the LB method with voxelized images of porous media. We show that, unlike many claims in the past, one may need many lattice units per voxel in order to obtain reliable, lattice-independent results.