Advection-dispersion equations (ADEs) are an important class of partial differential equations (PDEs) that are used to describe transport phenomena in the fields of hydrology (Ingham & Ma, 2005) and hydrogeology (Patil & Chore, 2014).We consider both forward and backward ADEs. In the former case, the initial conditions at time t 0 as well as boundary conditions are specified, and the solutions of forward ADEs are found for later times. This is a well-posed problem with well-established numerical methods. Nevertheless, there are some challenges in numerically solving forward ADEs, mainly associated with the advection-dominated problems. In backward ADEs, the concentration is known at later (terminal) times and solutions are sought for earlier times. Backward ADEs arise in the source identification problems (Neupauer & Wilson, 1999;Wilson & Liu, 1994) and could lead to numerically unstable grid-based solutions that require some form of regularization (Xiong et al., 2006) or should be solved as an inverse problem that is computationally more expensive because it requires solving forward problems multiple times (Atmadja & Bagtzoglou, 2001).Numerical discretization-based methods, including the finite element (FE) and finite difference (FD) methods, are commonly used for solving the Darcy flow equation and ADE. Discretization-based methods approximate the PDE solution with its values at a set of grid points distributed over the spatiotemporal domain. The discrete solution is obtained by discretizing the time and spatial derivatives of state variables. It is worth noting that the space-dependent parameters such as hydraulic conductivity are usually given not as a continuous field but as a set of values at the grid points.The combination of an advection (first-order) term and a dispersion (second-order) term in ADEs present several challenges for numerical methods. For example, for advection-dominated transport (i.e., Péclet number Pe ≫ 1), the numerical solutions can develop oscillations (over-or undershoot) or numerical dispersion (Huyakorn, 2012;Pinder & Gray, 2013). These two numerical issues are closely related, and a numerical scheme developed to reduce numerical dispersion generally causes oscillation, whereas the suppression of oscillation comes at the cost of increased numerical dispersion (Wang & Lacroix, 1997).