OR/17/068 Distributed models

Collins, S. 2017. Incorporating groundwater flow in land surface models: literature review and recommendations for further work. British Geological Survey. (OR/17/068).

Many LSMs ignore lateral groundwater flow on the basis that lateral fluxes between grid cells are very small. Krakauer et al. (2014) showed that significant groundwater flow (>10% of local recharge or 10 mm/year) occurs over 42% of the global land area at a resolution of 0.1° (~1 km), but that this drops to 1.5% at a resolution of 1˚ (~100 km). There are two principal advantages in using a distributed model: (1) a more accurate representation of groundwater–surface water interactions; and (2) a more accurate simulation of water table depth, and thus the effect of groundwater on evapotranspiration and climate (Anyah et al., 2008). There is, however, another reason why groundwater flow is often not represented in LSMs: namely a paucity of hydrogeological data. In this section, distributed models that have been incorporated into LSMs, as well as global distributed groundwater models, are compared. Table 2 summarises these models.

Models in land surface models

Gutowski et al. (2002) and York et al. (2002) were the first to demonstrate that a distributed groundwater model could be coupled with a single-column land surface–atmosphere model. Gutowski et al. (2002) developed their own simple groundwater model with 1D groundwater flow towards a central river running through the middle of a cell (cell boundaries were no flow boundaries). York et al. (2002) replaced the soil–vegetation and groundwater–surface water modules of the same land surface–atmosphere model with routines integrated into 3D MODFLOW (Harbaugh et al. 2000) for watershed-scale simulations (cell size 50–500 km). The soil–vegetation zone interacts with the aquifer through recharge to the aquifer and evapotranspiration directly from the water table; flow at the catchment outlet is the sum of streamflow and total leakage to or from the aquifer, calculated based on the head difference between the stream and the aquifer.

Fan et al. (2007) incorporated a simple 2D steady state groundwater flow model into an LSM at the continent scale (LEAF-Hydro). The groundwater model assumes that hydraulic conductivity decreases exponentially with depth beneath 1.5 m below ground. The decay factor in the relationship is a function of terrain slope, which the authors parameterised for regolith and bedrock based on concepts of weathering profiles (see Section 2.3). The same authors (Miguez-Macho et al., 2007) extended the model to a transient model with improved treatment of groundwater–surface water interactions. In their original work (Fan et al., 2007), groundwater was discharged to surface water when the water table reached the ground surface and surface water did not discharge to groundwater. This was built upon by Miguez-Macho et al. (2007), who used a statistical approach looking at mean geomorphological parameters across a grid cell (12.5 km resolution), the model resolution being too low for explicit treatment of individual channels. Owing to a lack of geomorphological data, they lumped the river bed hydraulic conductivity, river bed thickness, channel width and channel segment length into one ‘river conductance’ parameter (as in regional groundwater flow model codes such as MODFLOW), for which an equation was derived comprising equilibrium and dynamic parts. The dynamic part is a function of water table elevation and terrain slope and is calibrated based on river discharge observations. River elevation was fixed from the ‘naturally occurring’ rivers in the steady state run (1.25 km resolution) (Fan et al., 2007) or at the lowest ground surface elevation, when no river cells occurred within a 12.5 km cell.

Vergnes et al. (2012) added a 2D transient groundwater component to a hydrological model, which they later coupled with an LSM with a scheme for the unsaturated zone (Vergnes et al., 2014). The authors applied the model to France with a parameterisation technique that could be upscaled to the global scale (see Section 2.3). In contrast to most authors, Vergnes et al. (2012, 2014) applied the groundwater model only to the areas where major aquifers had been identified. The representation of the groundwater-surface water interaction is similar to that used in LEAF-Hydro (Miguez-Macho et al., 2007), in that all groundwater cells can exchange water with a ‘river’ and the rate is determined by a lumped parameter. Vergnes et al. (2012), however, lump only river bed conductivity and thickness into a ‘transfer time’ parameter and river width and length are calculated with empirical models (see Decharme et al. 2012; Vergnes et al. 2014). The transfer time parameter is dependent on its maximum and minimum values (taken from the literature) and on stream order. Vergnes et al. (2012) note that the system behaviour of the model is more sensitive to the transfer time parameter than to the hydrogeological properties of the aquifer, but also that changes in model performance achieved by varying this parameter are limited compared with the improvement in performance from including groundwater in the hydrological model. In their second paper (Vergnes et al., 2014), which considers capillary flux in the unsaturated zone, the authors reduced capillary flux based on the spatial variability of topography in a grid cell, such that more capillary rise occurs in flatter terrains.

Tian et al. (2012) coupled the AquiferFlow (Wang, 2007) groundwater model to an LSM, using this code to simulate both the saturated and unsaturated zones. They applied the model at the regional scale (grid resolution 3 km, total area ~13 000km²), incorporating an unconfined aquifer, an aquitard and a confined aquifer. The Heihe river, north-eastern China, and its major tributaries were represented as fixed head boundaries, and the model boundaries as well as the hydraulic conductivity were parameterised through calibration against groundwater level data.

The fully integrated groundwater–surface water platform ParFlow (Maxwell et al., 2017) has been coupled to the Terrestrial Systems Modeling Platform (TerrSysMP) atmospheric and land surface model (Shrestha et al., 2014), the Weather Research and Forecasting (WRF) atmospheric model (Maxwell et al., 2011), and the regional-scale meteorological model Advanced Regional Prediction System (ARPS) (http://www.caps.ou.edu/ARPS/). A key advantage of ParFlow is its explicit treatment of the groundwater–surface water interaction, which is either treated as a one-way drainage or parameterised with simple relationships (i.e. a functional relationship between river head and water table depth) in other models. The physically based approach requires little parameterisation, but is computationally extremely expensive, as shown by Maxwell et al. (2015), who applied the model at the continent scale.

Global groundwater models

Fan et al. (2013) presented the first global groundwater model. The model is highly simplified and has a number of limitations. The connection between surface water and groundwater is modelled only implicitly, in that when the groundwater table is above the land surface the excess water is removed. There is no recharge from surface water to groundwater. Moreover, as in their previous work (Fan et al., 2007; Miguez-Macho et al., 2007), no hydrogeological information, such as aquifer thickness or hydraulic conductivity, was used for parameterisation. Instead, a soil database (FAO, 1974) was used to obtain hydraulic conductivity close to the surface and it was then assumed to decrease exponentially. The model is also only steady state, requires calibration to head observations and does not take account of human influences, i.e. pumping, irrigation or drainage.

de Graaf et al. (2015, 2017) presented two global groundwater models of increased sophistication, which built on the group’s regional model (Sutanudjaja et al., 2011). They used global datasets on permeability and lithology and a digital terrain model to parameterise hydraulic conductivity and aquifer thickness as well as to delineate confining layers (see Section 2.3) (Gleeson et al., 2011; Hartmann and Moosdorf, 2012). As in LEAF-Hydro (Miguez-Macho et al., 2007), permeability was assumed to decrease exponentially with depth and the rate of decrease is controlled by the terrain slope. The models have a more sophisticated representation of groundwater–surface water interaction with three variations. (1) For larger rivers (width >10 m in first paper, >20 m in second), the MODFLOW river (RIV) package was used (same method as used by Miguez-Macho et al., 2007; Vergnes et al., 2012, 2014): that is, the rate of recharge/discharge is based on the head difference between the river and groundwater. The bed resistance parameter (combining bed thickness and conductivity) was set as a constant throughout the model (for which no justification was provided), river width was calculated with an empirical model and river length was assumed equal to the diagonal cell length. (2) Small rivers (width £10 m in first paper, £20 m in second) were simulated by a head dependent leakage function, similar to the MODFLOW drain (DRN) package, with water leaving the groundwater system only when the groundwater level exceeds the surface elevation. (3) An extra term, based on the digital elevation model and estimated storage, was added to account for rivers and springs in mountainous areas for which the model is too coarse to capture.

The first model (de Graaf et al., 2015) was steady state only, comprised a single-layer unconfined aquifer and failed to account for human impacts, i.e. abstraction and irrigation return flow. The model achieved good accuracy (coefficient of determination R² = 0.95, regression coefficient α = 0.84) against observed water levels in sedimentary basins, but tended to overestimate groundwater levels in mountainous areas, which the authors attribute to the exclusion of perched aquifers from the model. The second model (de Graaf et al., 2017) was transient, comprised a confined layer as well as the unconfined layer and included abstraction. de Graaf et al. (2017) achieved only a slight improvement in performance by simulating a confined layer. However, they claimed the fact that their estimate of global depletion is closer to that calculated by Konikow (2011) using a volume-based approach demonstrates that the groundwater–surface water interaction is better represented when the confined aquifer is included. However, it should be noted that there is considerable variation between different estimates of global groundwater depletion (from 113 km³/year [for period 2000–2009] to 330 km³/year [for year 2000]; Wada et al. 2010; Konikow, 2011; Döll et al., 2014; Pokhrel et al., 2015; de Graaf et al., 2017). Besides the simple approximation of the 3D hydrogeology, necessitated by a lack of data on the global scale, a major limitation of the model is that it is only one-way coupled: that is, the hydrological model is run for the entire simulation period and then time series of surface water levels, net recharge and groundwater abstractions from the hydrological model are passed to the groundwater model. Thus, there is no capillary rise from groundwater and the effects of pumping on surface water levels cannot be included.

Table 2    Distributed groundwater models on a global scale or coupled to land surface models.

Author	Year	Groundwater model	Parameterisation		Steady/ transient	GW–SW interaction	Resolution	Description
			Hydraulic conductivity and porosity	Aquifer extent and thickness
Fan et al.	2007	LEAF2-Hydro (LSM)	Vertical K from soil database (FAO, 1974). Anisotropy factor based on soil class to find lateral K. Exponential decrease with depth	Thickness represented by e-folding depth, for which an equation was derived based on slope	Steady	Implicit. When water table is above land surface, water removed as river discharge	1.25 km	North America model
	2013	LEAF2-Hydro	Vertical K from soil database (FAO, 1974). Anisotropy factor based on soil class to find lateral K. Exponential decrease with depth	Thickness represented by e-folding depth, which is function of slope	Steady	Implicit. When water table is above land surface, water removed as river discharge	30 arc-sec (roughly 1 km)	Global model
Miguez-Macho et al.	2007	LEAF2-Hydro (LSM)	Vertical K from soil database (FAO, 1974). Anisotropy factor based on soil class to find lateral K. Exponential decrease with depth	Thickness represented by e-folding depth, which is function of slope	Transient	Statistical approach. River conductance requires calibration against observed river flows	12.5 km	Applied to the USA
de Graaf	2015	MODFLOW	Gleeson et al. (2011) permeability map	Statistical method for thickness based on topography	Steady	Three different types depending on river size. MODFLOW river (RIV) package used for large rivers	5’	Single, unconfined aquifer. No capillary rise, groundwater pumping or recharge through irrigation return flows. Global model
	2017	MODFLOW	Gleeson et al. (2011) permeability map	Statistical method for thickness based on topography	Transient	Three different types depending on river size. MODFLOW river (RIV) package used for large rivers	5’	Unconfined and confined aquifers. No capillary rise from groundwater. Includes abstraction. Global model
Vergnes	2012	Based on MODCOU (Ledoux et al., 1989)	Transmissivity and effective porosity chosen based on typical values for given lithology	WHYMAPa, IGMEb and lithological maps used to delineate main aquifer basins. Slope also used	Transient	All cells are river cells. RC as in MODFLOW. Based on head in river and aquifer. River width calculated with empirical formula	0.5° and 1/12°	Single layer, unconfined aquifer. GW flow equation solved in spherical coordinates
	2014	Based on MODCOU (Ledoux et al., 1989)	Transmissivity and effective porosity chosen based on typical values for given lithology	WHYMAPa, IGMEb and lithological maps used to delineate main aquifer basins. Slope also used	Transient	All cells are river cells. RC as in MODFLOW. Based on head in river and aquifer. River width calculated with empirical formula	0.5° and 1/12°	Same as Vergnes et al. (2012), but includes capillary rise and unsaturated zone
Maxwell	2015	ParFlow	Gleeson et al. (2011) permeability maps	Assumed 100 m aquifer thickness everywhere	Steady	Modelled explicitly	1 km (over total area ~6.3 M km2)	Computationally expensive. No transient dynamics, human activities (e.g. pumping)
York et al.	2002	MODFLOW	From literature (assumed constant across watershed)	Assumed aquifer base	Transient	Proportional to aquifer-river head difference and river conductance (standard MODFLOW)	Variable	Single column LSM-atmosphere model. Watershed scale, single layer, groundwater flow within a cell
Gutowski et al.	2002	1D model	From literature (assumed constant)	Assumed from geological knowledge of the area	Transient	River at centre of model cell; river channel extended to aquifer base, no river conductance (flow controlled by aquifer K)	Variable	Single column LSM-atmosphere model. Groundwater flow within cell
Tian et al.	2012	AquiferFlow (unsaturated flow)	Treated as calibration parameters	Borehole logging data	Transient	Fixed head boundaries	3 km	Regional model. Unconfined aquifer, aquitard and confined aquifer. Fully two-way coupled (~ 13 000 km2)

^a WHYMAP is available at: http://www.whymap.org.
^b IGME (International Geological Map of Europe and Adjacent Areas) is available at www.bgr.bund.de.

Parameterisation

Global datasets that have been used for global groundwater modelling can be found in Table 3.

1.3.1 Hydraulic conductivity and porosity

Gleeson et al. (2011) developed global maps of permeability and porosity for consolidated and unconsolidated geological units up to a depth of 100 m (http://spatial.cuahsi.org/gleesont01/). The maps have been used in global-scale (de Graaf et al., 2015, 2017), continent-scale (Maxwell et al., 2015) and regional-scale (Shrestha et al., 2014) models. de Graaf et al. (2015, 2017) used the maps in combination with the high-resolution global lithology map of Hartmann and Moosdorf (2012). Fan et al. (2013), however, mention technical difficulties in using the maps of Gleeson et al. (2011) and instead used the FAO global soil map (FAO, 1974). In their North America model (Fan et al., 2007; Miguez-Macho et al., 2007), the authors used a US soil database that gives conductivity in the vertical direction, so they had to assume an anisotropy factor to determine lateral conductivity and then assumed an exponential decrease in conductivity with depth. Vergnes et al. (2012, 2014) used WHYMAP (http://www.whymap.org), the International Geological Map of Europe (https://www.bgr.bund.de) and a simple lithological map to delineate different geological formations in France and then assumed typical values of hydraulic conductivity and porosity based on the lithology.

Although the assumption of exponentially decreasing hydraulic conductivity with depth – often referred to as the TOPMODEL approach (Beven and Kirkby, 1976) – appears unsatisfactory, it as an assumption that can be found widely in the literature (Niu et al., 2007, 2011; Fan et al., 2007; Miguez-Macho et al., 2007; Fan et al., 2013; Koirala et al., 2014; de Graaf et al., 2015, 2017; Maxwell et al., 2015). Many authors (e.g. de Graaf et al., 2015, 2017; Maxwell et al., 2015) use the e-folding depth1 (α) first proposed by Fan et al. (2007):

𝐾(𝑧) = 𝐾0𝑒−𝑧/𝛼

where K is hydraulic conductivity, K0 is known hydraulic conductivity at the bottom of the soil layer and z is depth. The authors use the principle that erosion, weathering and deposition determine the decrease in conductivity with depth and that these processes are controlled by slope, climate and bedrock lithology. In order to simplify the approach they base their equation for e- folding on slope alone, such that the steeper the slope, the thinner the regolith. Fan et al. (2007) derived two relationships for e-folding against terrain slope: one for regolith and one for bedrock. The authors mention that the relationships were determined by trial and error, but it is not clear what data were used in deriving them.

1.3.2 Aquifer extent and depth

Shangguan et al. (2017) recently presented the first global map of depth to bedrock (250 m resolution) (https://www.soilgrids.org). Using machine learning algorithms, they derived the map from a global compilation of soil profile data (at 130,000 locations), borehole data (1.6 million locations) and pseudo-observations consisting of remote sensing data, terrain slope and geological maps. The authors warn of low accuracy in extrapolated areas (borehole data are from only eight countries) as well as in areas where depth to bedrock is > 100 m. Cross validation of the data set suggested moderate performance for absolute depth to bedrock.

Pelletier et al. (2016) created a map of average thickness of soil, intact regolith and sedimentary deposits (30 arcsec or ~ 1 km resolution). They used geomorphological models (both process

1 Here e-folding is used to denote the time interval in which an exponentially growing quantity increases by a factor of e; it is the base-e analog of doubling time. based and empirical) to calculate these thicknesses based on topographic, climate and geological data. The models were calibrated with data sets from the USA.

de Graaf et al. (2015, 2017) used statistical methods to determine aquifer thickness based on the assumption that mountain ranges have negligible sediment thickness and sediment basins below river valleys contain thicker, more productive aquifers. They used the difference between surface elevation and floodplain elevation within a cell to distinguish between mountain ranges and sediment basins.

Vergnes et al. (2012, 2014) attempted to delineate the main aquifer basins in France using only the WHYMAP global groundwater map, but found it too coarse. They also used the International Geological Map of Europe and a simplified lithological map, removing mountainous areas based on slope.

1.3.3 Categorising confined and unconfined aquifers

de Graaf et al. (2017) delineated confining layers and categorised aquifers into confined and unconfined using information on grain size and sediment properties from the global lithological map GLiM (Hartmann and Moosdorf, 2012). The authors introduced their own method for coastal zones, which are not fully represented in GLiM, classifying coastal zones around large rivers as confined (in total ~11% of global coastline). de Graaf et al. (2017) also made the simplifying assumption that the thickness of the confining layer is always 10% of the estimated aquifer thickness. OR/17/068

Table 3 Global data sets useful in global groundwater modelling

Author

Year

Data

Availability

Gleeson et al.

2011, 2014

Maps of permeability and porosity available at different resolutions

Available at http://crustalpermeability.weebly.com/glhymps.html

Shangguan et al.

2017

Map of depth to bedrock

https://www.soilgrids.org

Pelletier et al.

2016

Maps of average thickness of soil, intact regolith and sedimentary deposits

Available from author on request

FAO

1974

Global soil map

Used by Fan et al. (2007) and Miguez-Macho et al. (2007) to estimate hydraulic conductivity. A digital soil map is available at: https://worldmap.harvard.edu/data/geonode:DSMW_RdY

Bundesanstalt für Geowissenschaften und Rohstoffe

International Geological Map of Europe

Available at: https://www.bgr.bund.de

WHYMAP

Groundwater basins of the world

Available at: https://www.whymap.org/whymap/EN/Home/whymap_node.html

Fan et al.

2013

Map of simulated depth to water table. 1,603,781 groundwater head observations from across the world but predominately in the USA

Online database not found, but data has since been used in other studies (e.g. de Graaf et al., 2015), so can presumably be obtained from the authors

Wada et al.

2010

Global mapping of average recharge and abstraction (in the year 2000, 0.5 x 0.5 degree resolution) with the hydrological model PCR- GLOBWB. Used demand modelling and abstraction data from IGRAC to determine abstraction

Döll and Fiedler

2008

Global long-term average recharge mapping

International Groundwater Resources Assessment Centre (IGRAC)

–

Data portal including groundwater abstraction, water level data, etc.

https://ggis.un-igrac.org/ggis- viewer/viewer/exploreall/public/default