Summary: This paper extends the hierarchical multiscale approach developed earlier by the authors to model the coupled hydro-mechanical behaviour for saturated granular soils. Based on a hierarchical coupling of the finite element method (FEM) and the discrete element method (DEM), the approach employs the FEM to solve a boundary value problem while using the DEM to derive the required nonlinear material responses at each FEM Gauss integration point. It helps to bypass the phenomenological constitutive model required in the conventional FEM simulations and offers a natural pathway for scale bridging. The current extension features further key consideration of the coupled hydro-mechanical behaviour in a saturated granular soil due to the presence of pore fluid and its flow. By invoking Terzaghi’s effective stress principle, a \(\mathbf{u}-p\) formulation is proposed for the multiscale framework to derive the effective stress from DEM solution of the representative volume element (RVE) embedded at each Gauss point which is then superimposed with the pore fluid pressure to obtain the total stress in solving the coupled governing equations for a fluid-solid mixture. This extension greatly enriches the predictive capacity of the multiscale approach for simulating saturated granular soils relevant to a wide variety of important civil/mining applications. The new approach is first benchmarked by closed-form solutions to classical 1D and 2D consolidation problems. It is then applied to simulate globally undrained biaxial compression tests on both dense and loose soils subjected to large deformation. We further discuss interesting fluid flow patterns and the failure modes (both localisation and diffuse liquefaction) observed from the simulations and provide detailed cross-scale analyses.