► A new, terrain following grid formulation is developed and implemented in the parallel, integrated hydrologic model, ParFlow. ► The full, analytical Jacobian for this approach is developed and is used as a preconditioner for the solution technique. ► The terrain following formulation is verified against the orthogonal formulation for simple and complex cases. ► This approach scales efficiently out to 16,384 processors for a problem of 8.1Billion unknowns. ► This numerical advancement enables large-scale, high resolution, integrated hydrologic modeling, e.g., continental scale. A terrain-following grid formulation (TFG) is presented for simulation of coupled variably-saturated subsurface and surface water flow. The TFG is introduced into the integrated hydrologic model, ParFlow, which uses an implicit, Newton Krylov solution technique. The analytical Jacobian is also formulated and presented and both the diagonal and non-symmetric terms are used to precondition the Krylov linear system. The new formulation is verified against an orthogonal stencil and is shown to provide increased accuracy at lower lateral spatial discretization for hillslope simulations. Using TFG, efficient scaling to a large number of processors (16,384) and a large domain size (8.1 Billion unknowns) is shown. This demonstrates the applicability of this formulation to high-resolution, large-spatial extent hydrology applications where topographic effects are important. Furthermore, cases where the analytical Jacobian is used for the Newton iteration and as a non-symmetric preconditioner for the linear system are shown to have faster computation times and better scaling. This demonstrates the importance of solver efficiency in parallel scaling through the use of an appropriate preconditioner.