Surface subsidence and seepage damage in surrounding soils induced by leakage from municipal water supply pipelines pose significant risks to urban infrastructure. To clarify how leakage water diffuses in unsaturated soils and to assess seepage damage potential, this study established a numerical model based on the Richards equation combined with the van Genuchten (VG) model. The model was validated against physical model tests using remolded Q3 loess, ensuring consistency in soil parameters and leakage conditions. Simulation results reveal that soil saturation evolution follows three stages—initial, rising, and stable—with preferential flow paths forming above the leakage point before gradually evolving into radial diffusion controlled by both pressure and gravity. The extent of the saturated zone increases with pipeline pressure, but the enhancement effect diminishes as pressure rises, reflecting the nonlinear water-retention characteristics of loess. Seepage damage risk was evaluated using the Terzaghi critical hydraulic gradient criterion. The results show that higher pressures enlarge the critical zone more rapidly, yet its ultimate radius stabilizes within approximately 2.3 m around the leakage point. Moreover, this study proposes that potential seepage damage may occur once effective saturation reaches about 85%, corresponding to the air-entry value of loess, thus providing a more conservative criterion for engineering risk assessment. Overall, the validated Richards-based numerical model reproduces the key features of leakage-induced unsaturated diffusion and offers practical guidance for identifying seepage-prone zones and mitigating subsidence hazards in municipal water supply systems.