Anisotropy of hydraulic properties of the subsurface geologic formations is an essential feature that has been established as a consequence of the different geologic processes that undergo during the longer geologic time scale. With respect to subsurface reservoirs, in many cases, anisotropy plays significant role in dictating the direction of flow that becomes no longer dependent only on driving forces like the pressure gradient and gravity but also on the principal directions of anisotropy. Therefore, there has been a great deal of motivation to consider anisotropy into the subsurface flow and transport models. In this dissertation, we present subsurface flow modeling in single and dual continuum anisotropic porous media, which include the single-phase groundwater flow coupled with the solute transport in anisotropic porous media, the two-phase flow with gravity effect in anisotropic porous media, and the natural gas flow in anisotropic shale reservoirs. We have employed the multipoint flux approximation (MPFA) method to handle anisotropy in the flow model. The MPFA method is designed to provide correct discretization of the flow equations for general orientation of the principal directions of the permeability tensor. The implementation of MPFA method is combined with the experimenting pressure field approach, a newly developed technique that enables the solution of the global problem breaks down into the solution of multitude of local problems. The numerical results of the study demonstrate the significant effects of anisotropy of the subsurface formations. For the single-phase groundwater flow coupled with the solute transport modeling in anisotropic porous media, the results shows the strong impact of anisotropy on the pressure field and the migration of the solute concentration. For the two-phase flow modeling with gravity effect in anisotropic porous media, it is observed that the buoyancy-driven flow, which emerges due to the density differences between the phases, migrates upwards and the anisotropy aligns the flow directions closer to the principal direction of anisotropy. Lastly, for the gas flow modeling in anisotropic shale reservoirs, we observe that anisotropy affects the pressure fields and the velocity fields of the matrix and fracture systems as well as the production rate and cumulative production. It is observed from the results that all of the anisotropic cases produce higher amount of gas compared to isotropic case during the same production time. Furthermore, we have also examined the performance of MPFA with respect to mixed finite element (MFE) method over the lowest-order Raviart-Thomas (RT0) space and the first-order Brezzi-Douglas-Marini (BDM1) space. From the comparison of the numerical results we observe that MPFA method show very good agreement with the BDM1 than RT0. In terms of numerical implementation, however, MPFA method is easier than BDM1 and it also offers explicit discrete fluxes that are advantageous. Combining MPFA with the experimenting pressure field approach will certainly adds another advantage of implementing MPFA method as compared with RT0 and BDM1. Moreover, the computational cost (CPU cost) of the three different methods are also discussed.