Soil hydraulic properties relating saturation, water pressure, and hydraulic conductivity are known to exhibit hysteresis. In this paper, we focus on the determination of the water retention curve for a porous medium through a novel pore-scale simulation technique that is based on mathematical morphology. We develop an algorithm that allows for the representation of three-dimensional randomly packed porous media of any geometry (i.e. not restricted to idealized geometries such as spherical or ellipsoidal particles/pore space) so that the connectivity-, tortuosity-, and hysteresis-causing mechanisms are represented in both drainage and wetting processes, and their role in determining macroscopic fluid behavior is made explicit. Using this method, we present simulation results that demonstrate hysteretic behavior of wetting and non-wetting phases during both drainage and wetting cycles. A new method for computing interfacial surface areas is developed. The pore-morphology-based method is critically evaluated for accuracy, sample size effects, and resolution effects. It is found that the method computes interfacial areas more accurately than existing methods and allows for (i) examination of relationships between water pressure, saturation and interfacial area for hysteretic soils, and (ii) comparisons with previously developed theoretical models of soil hydraulic properties. The pore-morphology-based method shows promise for applications in vadose zone hydrology.