In this paper a 2D pollution water model with an improved numerical method is considered. In order to reduce the approximation errors of the numerical scheme, a new approximation method is introduced to calculate the concentration flux between two cells (j-cell and l-cell) in the direction of the normal vector \(\vec {n}\) orthogonal to their common side. The advantage of this approximation is that the concentration flux \({\partial C}/{\partial \vec {n}}\) from j-cell to l-cell and the other one from l-cell to j-cell are only different by their signs but not by their absolute values. Therefore, the errors of concentration simulated by this scheme are reduced and less than those obtained by a normal differential implicit discretization. This improvement of the scheme will be illustrated by two test cases. In the numerical tests we will display the difference between the exact solution, the classical scheme and the proposed scheme. Numerical results demonstrate the improvement of this approach.