Mathematical models have been widely used for the study of water quality. Regarding the pollution of water resources by organic matter, dissolved oxygen (DO) and biochemical oxygen demand (BOD) are important indicators for their monitoring. In this way, this research presents a two-dimensional mathematical model in space that evolves over time, for the study of the variations in concentrations of these variables, through a coupled system of non-linear partial differential equations, with Holling type III reaction kinetics between DO and BOD. Under some simplifications, the stability of the equilibrium points of the model is analyzed. An approximate solution is proposed using the centered finite difference (CFD) method for spatial variables and the Crank–Nicolson method for the temporal variable. An upwind scheme was employed in the advective-term discretization. Computer simulations were performed using Matlab software, in a rectangular domain. For the advective transport, the influences of the wind and the current given by a parabolic profile were considered. We performed a series of space-time numerical simulations and found that the model allows the analysis of regions with higher and/or lower DO and BOD concentrations, as well as the temporal variations of concentrations at specific points in the domain, and the influence of diffusive-advective transport in the mass transfer process.