One of the most important challenges facing academic, industrial and policy-making sectors is meeting with the increasing energy demand while preserving the affordability of the energy and maintaining the quality of the planet earth (including the environment, specially by reducing its greenhouse gas footprint). This global challenge demands for exploiting the subsurface formations as giant storage space for industrial by prod- ucts, e.g., CO2, while a future utilisation is found for them. Exploitation of subsurface formations for CO2 storage depends on our capacity to “accurately” and “efficiently” simulate the multiphase nonlinear flow in heterogeneous large-scale natural formations. An accurate and efficient simulation would allow for proper predictive understanding for several operational aspects including: (1) capacity of the storage; (2) life-time long behaviour of the injected CO2 in the subsurface formation, and (3) safety and integrity of the capturing procedure. The challenge of extreme scale dissimilarity between the heterogeneous coefficients (rock con- ductivity and fluid physics, e.g., mixing) and the reservoir has been known in the scientific literature as one of the main simulation challenges. Classically, to resolve this mismatch, reservoir models have been exces- sively upscaled to a resolution which can be solved affordably with the state-of-the-art commercial-grade simulators. However, as the rapid extension of the computational capacity, as well as the newly developed multiscale scalable nonlinear solvers, one needs to revisit our simulation frameworks to provide a scalable (to field scale) and accurate CO2 simulation framework. This thesis work focuses on extension of the recently developed Algebraic Dynamic Multilevel Method (ADM) to highly-nonlinear simulations of CO2-Brine multi- phase flow problems. ADM maintains simulation scalability by imposing fine-scale grids only where needed; and imposing coarser grids paired with accurate prolongation operators far from the sub-domains with sharp gradients. In this work, to provide a systematic study, the developed ADM framework is formulated such that it allows for both multiscale-based basis functions and also upscaling-based procedure. This is crucially im- portant as it reveals how the "solution-based” interpolations based on fine-scale heterogeneity can impact the simulation results; compared with the “effective-coefficient-based” approach. In brief, the novelty and contribution of this research is two fold: (1) development of a multiscale-based ADM method for CO2-Brine simulator and (2) systematic study to find the best model order reduction strat- egy when it comes to complex multiphase characteristics and rock heterogeneity. Two different classes of simulations are studied: (1) injection and evolution of the injected CO2 phase into the brine with capillary effects; and (2) miscible mixing process in which resolving front fingers into the residing less-mobile fluid imposes computational challenges. The study, as such, intends to reveal the details in “upscaling-based” vs. “multiscal-based” approach; and tends to provide a generic framework in which the scalability of the simu- lations are preserved by employing fine-scale grids only where and when needed. Several challenging sensitivity analysis are performed in which fine-scale solutions are compared with the multi-scale ADM and upscaling ADM solutions for immiscible and miscible fluid flow displacement. From these numerical experiments can be concluded that Multiscale-ADM outperforms the Upscaling-based ADM approach. More precisely with similar active grid cells, it provides a more accurate simulation method. As such, Multiscale-ADM casts a promising approach for next-generation simulators for CO2-Brine systems.