AbstractA Lagrangian method to simulate the advection, dispersion, and reaction of a single chemical, biological, or physical constituent within drinking water pipe networks is presented. This Lagrangian approach removes the need for fixed computational grids typically required in Eulerian and Eulerian-Lagrangian methods and allows for nonuniform computational segments. This makes the method fully compatible with the advection-reaction water quality engine currently used in EPANET. An operator splitting approach is used, in which the advection-reaction process is modeled before the dispersion process for each water quality step. The dispersion equation is discretized using a segment-centered finite-difference scheme, and flux continuity boundary conditions are applied at network junctions. A staged approach is implemented to solve the dispersion equation for interconnected pipe networks. First, a linear relationship between the boundary and internal concentrations is established for every pipe. Second, a symmetric and positive definite linear system of equations is constructed to calculate the concentrations at network junctions. Last, pipe internal concentrations are updated based on the junction concentrations. The solution generates exact results when the analytical solutions are available and leads to more accurate water quality simulations than advection-reaction-only water quality models, especially in the areas where dispersion dominates advection.