The partial differential equations that describe species transport and flow in permeable media in subsurface are highly nonlinear. The thermodynamic functions that provide local equilibrium conditions are also strongly nonlinear. The transport and the local equilibrium with specified boundary and initial conditions that define production from hydrocarbon reservoirs are solved numerically. We have advanced numerical solution of the flow problem, based on finite element discretization. A combination of the mixed finite element and discontinuous Galerkin method are advanced to solve multiphase-multicomponent flow in unfractured and fractured media. Despite orders of magnitude increase in conjunctional speed from our work in the last few years, many challenges remain to be addressed.