There are a growing number of large-scale, complex hydrologic models that are capable of simulating integrated surface and subsurface flow. Many are coupled to land-surface energy balance models, biogeochemical and ecological process models, and atmospheric models. Although they are being increasingly applied for hydrologic prediction and environmental understanding, very little formal verification and/or benchmarking of these models has been performed. Here we present the results of an intercomparison study of seven coupled surface-subsurface models based on a series of benchmark problems. All the models simultaneously solve adapted forms of the Richards and shallow water equations, based on fully 3-D or mixed (1-D vadose zone and 2-D groundwater) formulations for subsurface flow and 1-D (rill flow) or 2-D (sheet flow) conceptualizations for surface routing. A range of approaches is used for the solution of the coupled equations, including global implicit, sequential iterative, and asynchronous linking, and various strategies are used to enforce flux and pressure continuity at the surface-subsurface interface. The simulation results show good agreement for the simpler test cases, while the more complicated test cases bring out some of the differences in physical process representations and numerical solution approaches between the models. Benchmarks with more traditional runoff generating mechanisms, such as excess infiltration and saturation, demonstrate more agreement between models, while benchmarks with heterogeneity and complex water table dynamics highlight differences in model formulation. In general, all the models demonstrate the same qualitative behavior, thus building confidence in their use for hydrologic applications.

Groundwater is by far the largest unfrozen freshwater resource on the planet. It plays a critical role as the bottom of the hydrologic cycle, redistributing water in the subsurface and supporting plants and surface water bodies. However, groundwater has historically been excluded or greatly simplified in global models. In recent years, there has been an international push to develop global scale groundwater modeling and analysis. This progress has provided some critical first steps. Still, much additional work will be needed to achieve a consistent global groundwater framework that interacts seamlessly with observational datasets and other earth system and global circulation models. Here we outline a vision for a global groundwater platform for groundwater monitoring and prediction and identify the key technological and data challenges that are currently limiting progress. Any global platform of this type must be interdisciplinary and cannot be achieved by the groundwater modeling community in isolation. Therefore, we also provide a high‐level overview of the groundwater system, approaches to groundwater modeling and the current state of global groundwater representations, such that readers of all backgrounds can engage in this challenge.

A highly modular and scale-consistent Terrestrial Systems Modeling Platform (TerrSysMP) is presented. The modeling platform consists of an atmospheric model (Consortium for Small-Scale Modeling; COSMO), a land surface model (the NCAR Community Land Model, version 3.5; CLM3.5), and a 3D variably saturated groundwater flow model (ParFlow). An external coupler (Ocean Atmosphere Sea Ice Soil, version 3.0; OASIS3) with multiple executable approaches is employed to couple the three independently developed component models, which intrinsically allows for a separation of temporal–spatial modeling scales and the coupling frequencies between the component models.Idealized TerrSysMP simulations are presented, which focus on the interaction of key hydrologic processes, like runoff production (excess rainfall and saturation) at different hydrological modeling scales and the drawdown of the water table through groundwater pumping, with processes in the atmospheric boundary layer. The results show a strong linkage between integrated surface–groundwater dynamics, biogeophysical processes, and boundary layer evolution. The use of the mosaic approach for the hydrological component model (to resolve subgrid-scale topography) impacts simulated runoff production, soil moisture redistribution, and boundary layer evolution, which demonstrates the importance of hydrological modeling scales and thus the advantages of the coupling approach used in this study.Real data simulations were carried out with TerrSysMP over the Rur catchment in Germany. The inclusion of the integrated surface–groundwater flow model results in systematic patterns in the root zone soil moisture, which influence exchange flux distributions and the ensuing atmospheric boundary layer development. In a first comparison to observations, the 3D model compared to the 1D model shows slightly improved predictions of surface fluxes and a strong sensitivity to the initial soil moisture content.

Emphasizing the physical intricacies of integrated hydrology and feedbacks in simulating connected, variably saturated groundwater‐surface water systems, the Integrated Hydrologic Model Intercomparison Project initiated a second phase (IH‐MIP2), increasing the complexity of the benchmarks of the first phase. The models that took part in the intercomparison were ATS, Cast3M, CATHY, GEOtop, HydroGeoSphere, MIKE‐SHE, and ParFlow. IH‐MIP2 benchmarks included a tilted v‐catchment with 3‐D subsurface; a superslab case expanding the slab case of the first phase with an additional horizontal subsurface heterogeneity; and the Borden field rainfall‐runoff experiment. The analyses encompassed time series of saturated, unsaturated, and ponded storages, as well as discharge. Vertical cross sections and profiles were also inspected in the superslab and Borden benchmarks. An analysis of agreement was performed including systematic and unsystematic deviations between the different models. Results show generally good agreement between the different models, which lends confidence in the fundamental physical and numerical implementation of the governing equations in the different models. Differences can be attributed to the varying level of detail in the mathematical and numerical representation or in the parameterization of physical processes, in particular with regard to ponded storage and friction slope in the calculation of overland flow. These differences may become important for specific applications such as detailed inundation modeling or when strong inhomogeneities are present in the simulation domain.