In this paper, we present a new semi-analytical method based on the B-spline approximation for solving 2D time-dependent convection-diffusion-reaction equations to model transfer in the anisotropic inhomogeneous medium. The mathematical model is expressed as the initial-boundary value problem for quasi-linear parabolic equation with the second order elliptic spatial operator with mixed derivatives and variable coefficients. The time-stepping Crank-Nicolson scheme transforms the original equation into a sequence of quasi-linear elliptic partial differential equations. The approximate solution is sought as series over basis functions which are taken in the form of the tensor products of the  B-splines with centers distributed inside the solution domain. Due to the modification of spline basis, the final approximate solution satisfies the boundary conditions of the initial problem with any choice of the coefficients of the series. The numerical examples demonstrate the high accuracy of the proposed method in solving 2D convection-diffusion-reaction problems in single- and multi-connected domains.