The physical model is defined by a coupled system of seepage displacement for simulating chemical oil recovery numerically, formulated by three nonlinear partial differential equations concerning the pressure of Darcy-Forchheimer flow, the concentration and the component saturations. The pressure appears within the concentration equation and saturation equations and the Darcy-Forchheimer velocity controls the concentration and saturations. The flow equation is solved by the conservative mixed finite element method. The order of the accuracy is improved by the velocity. The conservative mixed volume element with characteristics is applied to compute the concentration, i.e., the diffusion is discretized by the mixed volume element and convection is computed by the method of characteristics. The method of characteristics has strong computational stability at sharp fronts and confirms high computational accuracy. The mixed volume element simulates diffusion, concentration and the adjoint vector function simultaneously. The nature of conservation is an important physical feature in the numerical simulation. The saturations of different components are computed by the method of characteristic fractional step difference and the computational work is shortened significantly by decomposing a three-dimensional problem into three successive one-dimensional problems and using the method of speedup. By using the theory and technique of a priori estimates of differential equations, convergence of the optimal second order in $l^2$ norm is obtained. Numerical examples are provided to show the effectiveness and viability of this method. This composite method can solve the challenging benchmark problem well.