A particulate model is developed to analyze the effects of steady state and transient seepage of water through a randomly-packed coarse-grained soil as an improvement to conventional seepage analysis based on continuum models. In the new model the soil skeleton and pore water are volumetrically coupled. In the first phase of the study, the concept of relative density has been used to define different compaction levels of the soil layers of a completely saturated pavement filter system and observe the seepage response to compaction. First, Monte-Carlo simulation is used to randomly pack discrete spherical particles from a specified Particle Size Distribution (PSD) to achieve a desired relative density based on the theoretical minimum and maximum void ratios. Then, a water pressure gradient is applied across one two-layer filter unit to trigger water seepage. The pore water motion is idealized using Navier Stokes (NS) equations which also incorporate drag forces acting between the water and soil particles. The NS equations are discretized using finite differences and applied to discrete elements in a staggered, structured grid. The model predicted hydraulic conductivities are validated using widely used equations. The critical water velocities, hydraulic gradients and flow within the xi saturated soil layers are identified under both steady state and transient conditions. Significantly critical transient conditions seem to develop. In the second phase of the study the model is extended to analyze the confined flow through a partly saturated pavement layer and unconfined flow from a retention pond into the surrounding saturated granular soil medium. In partly saturated soil, the water porosity changes resulting from water flow is updated using the Soil Water Characteristics Curve (SWCC) of the soil. The results show how complete saturation develops due to water flow following the water porosity Vs pressure trend defined by the SWCC. Finally, the model is used to predict the gradual reduction in the water level of a retention pond and the location of the free-surface. The free-surface is determined by differentiating the wet and dry zones based on the Heaviside step function modified NS equations.