The authors develop an Eulerian-Lagrangian localized adjoint method (ELLAM)-mixed finite element method (MFEM) solution technique for accurate numerical simulation of coupled systems of partial differential equations (PDEs) which describe complex fluid flow processes in porous media. The ELLAM, which was shown previously to outperform many widely used methods in the context of linear convection-diffusion PDEs, is employed to solve the transport equation for concentration. Since accurate fluid velocities are crucial in numerical simulations, the MFEM is used to solve the equations for pressure and Darcy velocity. This minimizes the numerical difficulties occurring in standard methods for approximating velocities caused by differentiation of pressure and then by multiplication by rough coefficients. The ELLAM-MFEM solution technique significantly reduces temporal errors, symmetrizes the governing transport equation, eliminates nonphysical oscillation and/or excessive numerical dispersion, conserves mass, and treats boundary conditions accurately. Numerical experiments show that the ELLAM-MFEM solution technique simulates miscible displacements of incompressible fluid flows in porous media accurately with fairly coarse spatial grids and very large time steps, which are one or two orders of magnitude larger than the time steps used in other methods. Moreover, the ELLAM-MFEM solution technique can treat large mobility ratios, discontinuous permeabilities and porosities, anisotropic dispersion in tensor form, and point sources and sinks.