A multi-region model for solute transport through saturated soils has been developed to describe preferential flow. The model consists of numerous discrete pore groups, which are characterized by a discrete dispersion coefficient, flow velocity, and porosity. The hydraulic properties for each pore group are derived from a soil's hydraulic conductivity and soil water characteristic functions. Flow in pore group is described by the classical advection-dispersion equation (ADE). An implicit finite difference scheme was applied to the governing equation that results in a block-tridiagonal system of equations that is very efficient and allows the soil to be divided into any number of pore groups. The numerical technique is derived from methods used to solve coupled equations in fluid dynamics problems and can also be applied to the transport of interacting solutes. The results of the model are compared to the experimental data from published papers. This paper contributes on the characteristics of the method when applied to the parallel porosity model to describe preferential flow of solutes in soil.