Mixed finite element discritizations for problems arising in flow in porous medium applications are considered. We first study second order elliptic equations which model single phase flow. We consider the recently introduced expanded mixed method. Combined with global mapping techniques, the method is suitable for full conductivity tensors and general geometry domains. In the case of the lowest order Raviart-Thomas space, quadrature rules reduce the method to cell-centered finite differences, making it very efficient computationally. We consider problems with discontinuous coefficients on multiblock domains. To obtain accurate approximations, we enhance the scheme by introducing Lagrange multiplier pressures along subdomain boundaries and coefficient discontinuities. This modification comes at no extra computational cost, if the method is implemented in parallel, using non-overlapping domain decomposition algorithms. Moreover, for regular solutions, it provides optimal convergence and discrete superconvergence for both pressure and velocity. We next consider the standard mixed finite element method on non-matching grids. We introduce mortar pressures along the non-matching interfaces. The mortar space is chosen to have higher approximability than the normal trace of the velocity spaces. The method is shown to be optimally convergent for all variables. Superconvergence for the subdomain pressures and, if the tensor coefficient is diagonal, for the velocities and the mortar pressures is also proven. We also consider the expanded mixed method on general geometry multiblock domains with non-matching grids. We analyze the resulting finite difference scheme and show superconvergence for all variables. Efficiency is not sacrificed by adding the mortar pressures. The computational complexity is shown to be comparable to the one on matching grids. Numerical results are presented, that verify the theory. We finally consider the mixed finite element discretizations for the nonlinear multi-phase flow system. The system is reformulated as a pressure and a saturation equation. The methods described above are directly applied to the elliptic or parabolic pressure equation. We present an analysis of a mixed method on non-matching grids for the saturation equation of degenerate parabolic type.