<p>Moisture content critically influences the matric suction and stress distribution within subgrade soils, directly affecting their load-bearing capacity. Understanding the internal mechanisms by which moisture variations degrade subgrade performance is essential. This study proposes a novel coupled simulation method based on the Volume of Fluid (VOF) theory to fully resolve the relationship between matric suction and stress conditions among the air–water–solid phases in unsaturated subgrades. By coupling Computational Fluid Dynamics (CFD) with the Discrete Element Method (DEM), the proposed approach incorporates phase fraction modeling to capture particle interactions and quantify the drag forces between phases. The results demonstrate that using phase fractions enables effective modeling of matric suction variations at the pore scale. The proposed method accurately captures the intermediate suction range (10⁴–10⁶&#xa0;kPa), with errors below 2% in low suction scenarios. It provides detailed insights into particle-scale stress redistribution under varying moisture conditions. Furthermore, the numerical simulations of resilient modulus of subgrade soils matches well laboratory measurements, with deviations within 3–5%, confirming the reliability and predictive capability of the proposed numerical approach.</p> Graphical Abstract <p></p>