An immersed boundary method (IBM) for the advection, diffusion, and reaction (ADR) of chemicals or species in incompressible flow is proposed. Suitable for complex geometries, the IBM with multidirect forcing is used in this paper to impose Dirichlet boundary conditions to surfaces of systems modeled by the Navier-Stokes, continuity, and ADR equations. Embedded in a Fourier Pseudospectral scheme, the method lends itself to high accuracy and computational efficiency. The numerical algorithm is verified using analytical solutions given by 2-D counter-rotating Taylor-Green vorticies, implemented with and without immersed boundaries. Processes of diffusion-reaction are verified using the Gray-Scott model. High order convergence and machine precision is observed. Two validation cases are then run. The first simulation validates the method’s performance applied to the Boussinesq approximation to model natural convection in a horizontal annulus. The next case simulates reactive flow over a circular cylinder. All results achieved are in good agreement with available data in literature.