Operator splitting is often used for solving advection-dispersion-reaction (ADR) equations. Each operator can be solved separately using an algorithm appropriate to its mathematical behavior. Although a lot of research has been done in operator splitting for solving ADR equations, numerical approaches for the reaction operator are computationally expensive. To meet the convergence criteria of ODE (ordinary differential equation) or DAE (differential algebraic equation) solvers, a transport time step has to be subdivided into a large number of reaction time steps. Additional computation effort is also required to reduce the splitting error. In this paper, we develop exact solutions of various first-order reaction networks for the reaction operator and couple those solutions with numerical solutions of the transport operator. The reactions are treated as local phenomena and simulated using exact solutions that we develop, while advection and dispersion are treated as global processes and simulated numerically. The proposed method avoids the numerical error from the reaction operator and requires a single-step calculation to solve the reaction operator. Compared to conventional operator-splitting methods, the proposed method offers both computational efficiency and simulation accuracy.