The dynamic solution of multiphase flow through porous media is of special interest to several fields of science and engineering, such as petroleum, geology and geophysics, bio-medical, civil and environmental, chemical engineering and many other disciplines. A natural application is the modeling of the flow of two immiscible fluids (phases) in a reservoir. Others, that are broadly based and considered in this work include the hydrodynamic dispersion (as in reactive transport) of a solute or tracer chemical through a fluid phase. Reservoir properties like permeability and porosity greatly influence the flow of these phases. Often, these vary across several orders of magnitude and can be discontinuous functions. Furthermore, they are generally not known to a desired level of accuracy or detail and special inverse problems need to be solved in order to obtain their estimates. Based on the physics dominating a given sub-region of the porous medium, numerical solutions to such flow problems may require different discretization schemes or different governing equations in adjacent regions. The need to couple solutions to such schemes gives rise to challenging domain decomposition problems. Finally, on an application level, present day environment concerns have resulted in a widespread increase in CO2 capture and storage experiments across the globe. This presents a huge modeling challenge for the future. This research work is divided into sections that aim to study various inter-connected problems that are of significance in sub-surface porous media applications. The first section studies an application of mortar (as well as nonmortar, i.e., enhanced velocity) mixed finite element methods (MMFEM and EV-MFEM) to problems in porous media flow. The mortar spaces are first used to develop a multiscale approach for parabolic problems in porous media applications. The implementation of the mortar mixed method is presented for two-phase immiscible flow and some a priori error estimates are then derived for the case of slightly compressible single-phase Darcy flow. Following this, the problem of modeling flow coupled to reactive transport is studied. Applications of such problems include modeling bio-remediation of oil spills and other subsurface hazardous wastes, angiogenesis in the transition of tumors from a dormant to a malignant state, contaminant transport in groundwater flow and acid injection around well bores to increase the permeability of the surrounding rock. Several numerical results are presented that demonstrate the efficiency of the method when compared to traditional approaches. The section following this examines (non-mortar) enhanced velocity finite element methods for solving multiphase flow coupled to species transport on non-matching multiblock grids. The results from this section indicate that this is the recommended method of choice for such problems. Next, a mortar finite element method is formulated and implemented that extends the scope of the classical mortar mixed finite element method developed by Arbogast et al (12) for elliptic problems and Girault et al (62) for coupling different numerical discretization schemes. Some significant areas of application include the coupling of pore-scale network models with the classical continuum models for steady single-phase Darcy flow as well as the coupling of different numerical methods such as discontinuous Galerkin and mixed finite element methods in different sub-domains for the case of single phase flow (21, 109). These hold promise for applications where a high level of detail and accuracy is desired in one part of the domain (often associated with very small length scales as in pore-scale network models) and a much lower level of detail at other parts of the domain (at much larger length scales). Examples include modeling of the flow around well bores or through faulted reservoirs. The next section presents a parallel stochastic approximation method (68, 76) applied to inverse modeling and gives several promising results that address the problem of uncertainty associated with the parameters governing multiphase flow partial differential equations. For example, medium properties such as absolute permeability and porosity greatly influence the flow behavior, but are rarely known to even a reasonable level of accuracy and are very often upscaled to large areas or volumes based on seismic measurements at discrete points. The results in this section show that by using a few measurements of the primary unknowns in multiphase flow such as fluid pressures and concentrations as well as well-log data, one can define an objective function of the medium properties to be determined, which is then minimized to determine the properties using (as in this case) a stochastic analog of Newton's method. The last section is devoted to a significant and current application area. It presents a parallel and efficient iteratively coupled implicit pressure, explicit concentration formulation (IMPEC) (52-54) for non-isothermal compositional flow problems. The goal is to perform predictive modeling simulations for CO2 sequestration experiments. While the sections presented in this work cover a broad range of topics they are actually tied to each other and serve to achieve the unifying, ultimate goal of developing a complete and robust reservoir simulator. The major results of this work, particularly in the application of MMFEM and EV-MFEM to multiphysics couplings of multiphase flow and transport as well as in the modeling of EOS non-isothermal compositional flow applied to CO2 sequestration, suggest that multiblock/multimodel methods applied in a robust parallel computational framework is invaluable when attempting to solve problems as described in Chapter 7. As an example, one may consider a closed loop control system for managing oil production or CO2 sequestration experiments in huge formations (the "instrumented oil field"). Most of the computationally costly activity occurs around a few wells. Thus one has to be able to seamlessly connect the above components while running many forward simulations on parallel clusters in a multiblock and multimodel setting where most domains employ an isothermal single-phase flow model except a few around well bores that employ, say, a non-isothermal compositional model. Simultaneously, cheap and efficient stochastic methods as in Chapter 8, may be used to generate history matches of well and/or sensor-measured solution data, to arrive at better estimates of the medium properties on the fly. This is obviously beyond the scope of the current work but represents the over-arching goal of this research.