A numerical method based on the modified method of characteristics is developed for incompressible Darcy flow. Fluid elements modeled as grid cells are mapped back in time to their twisted forms and a strict equality of volumes is imposed between the two. These relations are then cast in terms of potentials using Darcy's law and a nonlinear algebraic problem is solved for potentials. Though a general technique for obtaining Darcy flow, this method is most useful when the solute advection problem also is solved with the modified method of characteristics. The combined technique (referred to as the characteristic-conservative method) using the same characteristics to obtain both velocities and concentrations is then a direct numerical approximation to the Reynolds transport theorem. The method is implemented in three dimensions and a few sample problems featuring nonuniform flow-fields are solved to demonstrate the exact mass conservation property. Inflow and outflow boundaries do not cause any problems in the implementation, in all cases, the characteristic-conservative method obtains velocities that preserve fluid volume and, concentrations that achieve exact local and global mass balance; a desirable property that usually eludes characteristics based methods for solute advection in multidimensional, nonuniform flowfields. (C) 1999 Elsevier Science Limited. Ail rights reserved.