Lateral subsurface flow plays an essential role in sustaining the terrestrial ecosystem, but it is not explicitly represented in most Earth System Models. In this study, we implemented a hybrid-3D hillslope hydrology model into the Energy Exascale Earth System Model land model (ELM). Our model explicitly describes lateral flow in the saturated zone by representing, for each model grid, an idealized hillslope consisting of 5 hydrologically connected columns. We conducted three model experiments driven by 0.125° atmospheric forcing data during 1980–2015 over California including: the default ELM, a modified version of ELM to improve the infiltration modeling, and our model with the treatment of lateral flow. The simulated runoff, evapotranspiration, and terrestrial water storage anomaly (TWSA) from the three simulations were evaluated against available observations and our model performs best. The new model produces greater gridcell-averaged evapotranspiration especially over the mountainous regions with large topographic slopes and wet climates. Most importantly, it improves the modeled seasonal variations, interannual variabilities, and the recent decadal decline of TWSA. Many of these improvements can be attributed to enhanced ecosystem resilience to droughts demonstrated by the transpiration increases in this newly developed model. Model sensitivity experiments suggest that subsurface runoff is most sensitive to the ratio between horizontal and vertical saturated hydraulic conductivity, followed by hillslope planforms (convergent, divergent, and uniform), number of columns, and lower boundary conditions. Future work should effectively characterize hillslopes in global models and explore the long-term influences of lateral water movement on modeled biogeochemical cycle.