We propose a novel numerical methodology to compute the advective transport and diffusion-reaction of tracer quantities. The tracer advection occurs through flow map composition and is super-accurate, yielding numerical solutions almost devoid of compounding numerical errors, while allowing for direct parallelization in the temporal direction. It is computed by implicitly solving the characteristic evolution through a modified transport partial differential equation and domain decomposition in the temporal direction, followed by composition with the known initial condition. This advection scheme allows a rigorous computation of the spatial and temporal error bounds, yields an accuracy comparable to that of Lagrangian methods, and maintains the advantages of Eulerian schemes. We further show that there exists an optimal value of the composition timestep that yields the minimum total numerical error in the computations, and derive the expression for this value. We develop schemes for the addition of tracer diffusion, reaction, and source terms, and for the implementation of boundary conditions. Finally, the methodology is applied in three flow examples, namely an analytical reversible swirl flow, an idealized flow exiting a strait undergoing sudden expansion, and a realistic ocean flow in the Bismarck sea. New benchmark problems for advection-diffusion-reaction schemes are developed and used to compare and contrast results with those of classic schemes. The results highlight the theoretical properties of the methodology as well as its efficiency, super-accuracy with minimal numerical errors, and applicability in realistic simulations.