Despite their potential influences on surface water and climate, groundwater processes are generally not represented in climate models. Here, a simple groundwater scheme including two-dimensional flow dynamics and accounting for groundwater–river exchanges is introduced into the global Total Runoff Integrated Pathways (TRIP) river routing model coupled to the Météo-France climate model. This original scheme is tested in offline mode over France at high () and low (0.5°) resolution against a dense network of river discharge and water table observations over the 1970–2010 period, and is compared to the fine-tuned Système d’Analyze Fournissant des Renseignements Atmosphériques à la Neige (SAFRAN)–Interactions between Soil, Biosphere, and Atmosphere (ISBA) coupled hydrometeorological model (MODCOU). In addition, the simulated terrestrial water storage (TWS) variations are compared to the TWS estimates from the Gravity Recovery and Climate Experiment (GRACE) satellite mission. The aquifer basins over France are defined using the World-wide Hydrogeological Mapping and Assessment Programme (WHYMAP) groundwater resources map, a simplified French lithological map, and the International Geological Map of Europe (IGME). TRIP is forced by daily runoff and drainage data derived from a preexisting simulation of the ISBA land surface scheme driven by the high-resolution SAFRAN meteorological analysis. Four simulations are carried out with or without groundwater at both resolutions. Results show that the groundwater scheme allows TRIP to better capture the spatiotemporal variability of the observed river discharges and piezometric heads. Summer base flows are particularly improved over the main rivers of France. Decreasing the horizontal resolution has a limited impact on the simulated discharges, while it slightly degrades the simulation of water table variations.
In climate models, the land surface hydrology has a major influence on the terrestrial water and energy budgets and, thereby, on the simulated weather and climate (Dirmeyer 2000, 2001; Douville 2003, 2004; Koster et al. 2000, 2002). It can affect the temperature and ocean salinity at the mouths of the largest rivers (Durand et al. 2011), the water and energy exchanges at the land surface, and the climate, at least at the regional scale (Gedney et al. 2000; Douville et al. 2000a,b; Molod et al. 2004; Lawrence and Slater 2007; Alkama et al. 2008). These land surface processes are parameterized in continental hydrological systems (CHSs) based on two components: 1) the land surface models (LSMs), which provide realistic lower boundary conditions of temperature and moisture in atmospheric general circulation models (AGCMs), and 2) the river routing models (RRMs), which convert the total runoff provided by LSMs into river discharges in order to evaluate the simulated water budget and/or to transfer continental freshwater to the oceans, thereby closing the global hydrological cycle.
However, many LSMs used in climate modeling still neglect the representation of the groundwater processes. Groundwater constitutes about 30% of the world’s total freshwater resources (Shiklomanov and Sokolov 1983)—much more than soil moisture (0.05%) and rivers (0.006%). It interacts with surface water and is therefore likely to influence the surface energy and water exchanges with the lower atmosphere. Its slow response to climate variations helps to maintain base flows in humid climates during dry periods, while it receives the river seepage in arid climates. Water table rise and fall can also interact with the soil moisture profile and thereby affect evapotranspiration and the land surface energy budget (e.g., Dingman 1994).
During the last decade, several studies have pointed out the importance of including groundwater processes in CHSs. Van den Hurk et al. (2005) analyzed seven regional climate models with respect to the land surface hydrology over the Rhine basin. They found that insufficient water storage led to overestimation of the seasonality of the simulated runoff compared to the observations. Through observations and model simulations carried out in Illinois, Yeh and Eltahir (2005) demonstrated that the free-drain or no-drain soil bottom conditions commonly used in LSMs could significantly affect the simulated soil water budget and river discharges. Alkama et al. (2010) compared global hydrological outputs from the Interactions between Soil, Biosphere, and Atmosphere–Total Runoff Integrating Pathways (ISBA–TRIP) CHS to in situ river discharges and terrestrial water storage (TWS) variations derived from the Gravity Recovery and Climate Experiment (GRACE). They found that an underestimation of continental evaporation and an overestimation of the annual discharges were likely due to the lack of a groundwater reservoir in ISBA–TRIP.
In this context, several attempts have been made to introduce groundwater processes into CHSs. To represent the groundwater flow contribution to the river, some studies proposed the addition of a simple pseudogroundwater reservoir into RRMs using a time delay factor to delay the flow to the river (Arora and Boer 1999; Decharme et al. 2010). Though useful for a better evaluation of CHSs against TWS estimates and/or discharge observations, such a method does not account for groundwater dynamics.
Other studies introduced a groundwater component in one-dimensional LSMs for global climate applications (Gedney and Cox 2003; Liang et al. 2003; Maxwell and Miller 2005; Yeh and Eltahir 2005; Niu et al. 2007; Ngo-Duc et al. 2007; Lo et al. 2010). Most of them represent the groundwater reservoir as a new deep layer under the soil column, which interacts with the unsaturated zone. Gedney and Cox (2003) added a very deep soil layer using a prescribed depth in order to represent shallow aquifers. Niu et al. (2007) used a single layer to resolve the water table depth as the lower boundary condition of the soil column. Liang et al. (2003) and Maxwell and Miller (2005) proposed a more realistic approach by explicitly coupling an LSM and a groundwater scheme in order to simulate the saturated zone and the overlying unsaturated zone of the soil as a continuum soil column and to explicitly compute the shallow water table position.
Nevertheless, most of these models are not coupled with RRMs and neglect the interactions with the river network, making it difficult to evaluate them against in situ river discharges.
Regional hydrometeorological studies use more detailed two-dimensional groundwater schemes. Gutowski et al. (2002) and York et al. (2002) coupled a single-column atmospheric model directly with a high-resolution land surface scheme itself coupled with a detailed two-dimensional groundwater model. Lateral flows occurred in each terrestrial grid cell and groundwater reservoirs exchanged water with the river network. Their findings confirm the potential feedback between groundwater and atmospheric forcing and suggest the feasibility of incorporating physically based representations of aquifers into LSMs. Recently, Fan et al. (2007) and Miguez-Macho et al. (2007) developed a more sophisticated two-dimensional hydrological model based on a groundwater diffusive scheme taking the interactions with the land surface into account. This model was applied at a very fine resolution (~1 km) over the United States. The results showed good agreement between observed and simulated river discharges, and revealed that groundwater enhanced the memory of soil moisture processes.
At Météo-France, the operational Système d’Analyze Fournissant des Renseignements Atmosphériques à la Neige (SAFRAN)–ISBA coupled model (MODCOU) (SIM) hydrometeorological system is used to monitor real-time water resources over France at an 8-km horizontal resolution (Habets et al. 2008). SIM is composed of three modules: the SAFRAN meteorological analysis provides high-quality atmospheric forcings (Durand et al. 1993), the ISBA land surface model computes the surface water and energy budgets (Noilhan and Planton 1989), and the MODCOU hydrogeological model computes the evolution of the aquifers, the river flows, and the exchanges between them (Ledoux et al. 1989). Two aquifer basins are represented: the Seine and the Rhone basins (Habets et al. 1999; Etchevers et al. 2001; Rousset et al. 2004). The representation of the groundwater processes appears especially relevant over the Seine basin to simulate a more realistic summer base flow (Rousset et al. 2004).
Nevertheless, these regional, fine-tuned models present an important limitation for large-scale applications: many parameters have to be calibrated using high-resolution geological and topographic data as well as daily streamflow observations. Such a method cannot be transposed to the global scale given the lack of data. This is the reason why global applications using simplified groundwater parameterizations in LSMs generally assume uniform parameters across different climatic and hydrological regions (Arora and Boer 1999; Niu et al. 2007; Decharme et al. 2010).
The main objectives of this paper are to describe an original simple groundwater scheme, to evaluate its influence on the daily river discharges simulated at a relatively high resolution by the TRIP RRM used at the Centre National de Recherches Météorologiques (CNRM), and to demonstrate its robustness and suitability for lower-resolution applications using a simple methodology based on available global datasets to estimate the aquifer geometry and parameters. The model is tested over France (550 000 km2) at high ( × ) and low resolution (0.5° × 0.5°) over the 1970–2010 period. High-resolution simulations, with and without the groundwater scheme, are first compared to observations and to the SIM benchmark, and then provide a reference for the lower-resolution integrations. All TRIP simulations are forced by the same daily surface runoff and deep drainage derived from a preexisting SIM experiment carried out by Vidal et al. (2010a,b). These inputs are not affected by the MODCOU hydrogeological model, which does not interact with the ISBA land surface model within SIM. Evaluation is made against two dense networks of river gauging stations and piezometric gauges. The simulated TWS variations are also compared to the estimates from GRACE.
2. The surface-groundwater representation
a. Overview of the TRIP RRM
The basic concepts and equations of TRIP are described in Decharme et al. (2012) and are simply summarized here. The original TRIP RRM was developed at Tokyo University by Oki and Sud (1998). It is used at Météo-France to convert the total runoff simulated by ISBA into river discharges using a river network at 0.5° or 1° resolution in global simulations. TRIP is based on a prognostic mass balance equation for the stream water mass, which is solved at a 30-min time step on each cell of the river network:
where S (kg) is the stream water mass, (kg s−1) is the sum of the surface runoff from ISBA within the grid cell with the water inflow from the upstream neighboring grid cells, (kg s−1) is the deep drainage from ISBA, and (kg s−1) is the river discharge into the downstream cell, which is computed as follows:
where υ (m s−1) is the streamflow velocity and L (m) the length of the river inside the cell. The streamflow velocity υ is computed via Manning’s formula as proposed by Arora and Boer (1999) and described in detail in Decharme et al. (2010). Note that unlike in Decharme et al. (2010), no pseudogroundwater reservoir using a time delay factor to delay the flow to the river is used in this study.
b. Description of the new groundwater scheme
The proposed groundwater scheme is based on the MODCOU hydrogeological model (Ledoux et al. 1989). The main difference is that the two-dimensional groundwater flow equation expressed in m s−1 is solved in spherical coordinates to simulate the dynamics of the water table H:
Only one layer representing an unconfined aquifer is simulated in the model. The specific yield, ω (m3 m−3), corresponds to the effective porosity; θ and ϕ are the longitude and latitude coordinates, respectively; r (m) is the mean radius of the earth; Tθ and Tϕ (m2 s−1) are the transmissivities along the longitude and latitude axes, respectively; qriv (m s−1) the groundwater–river flux; and qsb (m s−1) is the deep drainage from ISBA per unit area. In other words, qsb = Qsb/Acell where Acell (m2) is the gridcell area. As in the ISBA–TRIP CHS, spherical coordinates are used because of the spherical form of the earth, especially in the high southern and northern latitudes. Omitting this detail could lead to the simulated lateral flows being underestimated along the latitude axis and overestimated along the longitude axis. It could also impact the mass conservation in the high-latitude grid cells [see the discretized form of Eq. (3) in the appendix]. Equation (3) is then solved in m3 s−1 with a time step of 1 day using an implicit finite-difference numerical method. More details can be found in the appendix.
In TRIP, each grid cell is considered as a river cell. As a result, each of them can potentially exchange water with the river (gaining or losing streams). This conceptual approach appears relevant at the resolutions considered ( and 0.5°). These exchanges are represented through the concept of a river coefficient RC commonly used in a majority of regional groundwater models such as MODCOU (Ledoux et al. 1989) or the Modular Three-Dimensional Groundwater Flow Model (MODFLOW; McDonald and Harbaugh 1988). The fundamental assumption of this approach is to consider that the head losses between the stream and aquifer are limited to those across the streambed itself. The river coefficient for a river channel of width W (m) and length L (m) with a riverbed of thickness b (m) and hydraulic conductivity Kriv (m s−1) equals LW(Kriv/b). Because b and Kriv are generally poorly known, uncertainties in estimating the riverbed properties make it necessary to adjust this coefficient through model calibration. The b/Kriv (s) quantity represents the duration of water flow through the riverbed. In TRIP, this quantity is approximated to a coefficient, τ(s), representing the time transfer coefficient between river and groundwater. The groundwater–river flow is therefore parameterized as follows:
where Zbed (m) is the riverbed elevation, which is the elevation in the grid cell minus the river bankfull height hc (m), defined in Decharme et al. (2012), and Hriv (m) is the river stage elevation, calculated as the sum of Zbed (m) with the river water height hs (m) (Fig. 1).
Equation (4a) corresponds to the case where the water table is connected to the river. If the water table falls under the riverbed elevation, Eq. (4b) is applied and the river feeds the groundwater reservoir. If the river height hs falls under 10 cm and Hriv is lower than H, Qriv is set to zero in order to avoid a completely empty river and/or negative discharges. The drainage term Qsb from ISBA in Eq. (1) now feeds the groundwater scheme. Consequently, this term is replaced by the river–groundwater exchange flux Qriv and Eq. (1) is thus modified as follows:
3. Experimental design
a. Model parameters
1) TRIP parameters
The elevation at resolution is derived from the Global 30 Arc-Second Elevation Data Set (GTOPO30) digital elevation model (DEM; http://eros.usgs.gov/#/Find_Data/Products_and_Data_Available/gtopo30_info). Because each grid cell in TRIP is defined as a river cell, the elevation appears as a critical parameter and is derived from a specific interpolation. Each gridcell topography at resolution is computed as the average value of the first decile of the actual 30 arc-second-resolution topographic values within the grid cell, ranked in ascending order. The result is corrected using GIS processing in order to remove elevation anomalies that could interfere with a hydrological correct flow. As described previously, this elevation helps us to compute the riverbed elevation Zbed, and reflect the altitude of the river in the grid cell. This approach allows more realistic river slopes, flow directions, stream orders, and watershed areas to be computed (Fig. 2a).
Moving to the low-resolution parameter estimation, elevation and slope are derived from the data previously calculated by taking the average of the gridcell values within each 0.5° cell. Indeed, the data are expected to be representative of the river network characteristics in each 0.5° grid cell. This approach computes elevation and slope along the river at 0.5° resolution with a better accuracy than direct use of the GTOPO30 data does. Flow directions are based on the 0.5° × 0.5° TRIP network, including some manual adjustments to correct the river basin boundaries (Fig. 2b).
At both resolutions, the river width is also an important parameter because it is used in the river flow velocity computation as well as in the groundwater–river exchanges parameterization. It is calculated through the empirical formulation described in detail in Decharme et al. (2012):
where Qyr is the annual mean discharge in each grid cell computed at each resolution using the annual mean total runoff simulated by ISBA over the 1958–2010 period. In this study, the β coefficient is fixed to 15, according to the values usually used for temperate basins (Arora and Boer 1999; Moody and Troutman 2002; Decharme et al. 2012). This method leads to a river width of 548 and 557 m at the Loire outlet for the and 0.5° resolutions, respectively, to 404 and 394 m at the Seine outlet, to 732 and 726 m at the Rhone outlet, and to 527 and 537 m at the Garonne outlet. Details about the other TRIP specific parameters such as the river length, L, or the height of the riverbed, hc, can be found in Decharme et al. (2012).
2) Aquifer characteristics
Since the groundwater scheme presented in this study was developed for global climate applications, only major regional groundwater basins are taken into account (Figs. 2c,d). The question is how these domains should be delimited. Basically, the main features that characterize the geometry and the properties of an aquifer are the topography, the lithology, and the geology. Only widespread unconfined aquifers concerned by diffusive groundwater movements are simulated by the model. These processes are preferentially located in sedimentary basins with regional aquifer formations constituted by permeable porous and fractured rocks, and mainly in alluvial plains characterized by gravel and sand materials with high permeability. The French national database on water systems (http://www.sandre.eaufrance.fr) provides a map of such aquifers that lie in the Paris and the Aquitaine basins and the Rhone and the Rhine valleys (Fig. 3). However, as the Référentiel Hydrogéologique Français (BDRHF) database is only available over France, whereas the model is intended to be used at the global scale, it was necessary to build a methodology to derive the aquifer limits in relation with the databases available at global scale. The methodology developed was the following.
A global map of the groundwater resources of the world edited by the World-wide Hydrogeological Mapping and Assessment Programme (WHYMAP; http://www.whymap.org) was used as the primary information to delimit the main aquifer basins (Fig. 4a). The main aquifer units of the BDRHF correspond to the classes’ “major groundwater basins” and “complex hydrogeological structures,” thus confirming the good consistency between the finescale BDRHF and the large-scale WHYMAP databases. So the classes’ “local and shallow aquifer” and “karstic areas” could be neglected. However, the two selected classes of the WHYMAP dataset cover too large a domain because of the low accuracy of the WHYMAP data. This includes some karstic areas in the south of France, as well as the east of the Paris basin, mainly composed of siliciclastic rocks such as claystone or sandstone and characterized by mixed confined and unconfined aquifers not intended to be simulated by TRIP. Consequently, to refine the groundwater basin ;boundaries in TRIP, the International Geological Map of Europe (IGME) provided by the Federal Institute for Geosciences and Natural Resources in Hannover (BGR) (Asch 2006; Fig. 4b) and the simplified lithological map of the Bureau des Ressources Géologiques et Minières (BRGM; Fig. 4c) were used as the second criterion in the proposed methodology.
The comparison of the BDRHF map (Fig. 3) with the two selected classes of WHYMAP, the IGME geological map, and the BRGM lithological map show that the main aquifers may be limited to the youngest geological formations and that the older areas—for example, dating from the middle Jurassic and the Triassic, with a majority of siliciclastic rocks (Figs. 4b,c)—can be removed to refine the aquifer domains. This filter is also well correlated with the old platforms of Brittany and the Massif Central already removed in WHYMAP. Finally, the third criterion was based on the slope and was defined in order to eliminate the remaining mountainous areas. Using the slope from the GTOPO30 elevation dataset, the cells defined as mountainous at the TRIP resolution ( or 0.5°) were those having at least 70% of slopes greater than or equal to 10% (0.1 m m−1) in the 30 arc-second grid cell. This criterion, when applied over France, masked the remaining part of the Alps. The final aquifer domains are shown in Figs. 2c,d. It is important to note that the methodology used to delimit the main groundwater basins over France is also suitable for global applications since sufficient geological and lithological data are available (Dürr et al. 2005; Bouysse 2009; Gleeson et al. 2011).
The time transfer coefficient τ varies arbitrarily from 30 days in major river streams to 5 days in the upstream grid cells through a linear relationship with the river stream order, SO, given by the TRIP network in each grid cell of a given basin:
where τmax and τmin are the maximum (30 days) and the minimum values (5 days) of τ chosen to be consistent with the time delay factors used in previous studies (Arora et al. 1999; Decharme et al. 2010). Here, SOmax is equal to 10 and 5 at high and low resolution, respectively, and SOmin is the minimum stream order in each basin of the TRIP network, equal to 1.
To estimate the transmissivity and the effective porosity in each grid cell, the simplified lithological map of France was used (Fig. 4c). Five main units of lithology material were selected in the domains where groundwater flows were simulated: clay, chalk, limestone, sandstone, and sand. Mean values of transmissivity and effective porosity were chosen among the usual values so as to be physically realistic considering this lithology and are summarized in Table 1.
TRIP was integrated at high () and low (0.5°) resolution using a 30-min time step. For each resolution, an offline hydrological simulation with the groundwater scheme was compared to a control experiment without groundwater. The simulations were named as follows:
NOGW12: control simulation without groundwater at high resolution.
NOGW05: same as NOGW12 but at low resolution.
GW12: groundwater and surface water simulation at high resolution.
GW05: same as GW12 but at low resolution.
In addition to both GW12 and GW05, parallel experiments were also performed using τ ± 75% in order to explore the model sensitivity to this empirical parameter.
TRIP is forced by the surface runoff and the deep drainage coming from an independent high-resolution (8 km) ISBA simulation covering the 1958–2010 period. Further details can be found in Vidal et al. (2010a,b). These forcing fields were produced in a 50-yr study by Habets et al. (2008) aimed at evaluating the SIM system over France and were interpolated at the TRIP resolution using a nearest-neighbor interpolation method via GIS processing. The forcing fields at the 0.5° resolution were computed by aggregating the runoff and drainage obtained at . Every day, TRIP computes piezometric heads and river discharges. To start the model at equilibrium, the water table was initialized to the topography and a first 60-yr spinup simulation was done by repeating the 1958–70 period five times, until the water table reached equilibrium. The model was then evaluated over the 1970–2010 period.
c. Evaluation datasets
In situ river daily discharge measurements at 318 gauging stations were selected to evaluate the simulated discharges at high resolution. Only the time series with a minimum period of 10 years were used. Moreover, when several gauging stations were located in one cell, the largest observed drainage area was kept. The same approach was applied for evaluating the low-resolution simulation. Only 99 gauging stations were selected. When several stations were present in the same grid cell, only the station with the largest upstream drainage area was conserved.
Piezometric heads are well monitored over France and numerous data are available in the Accès aux Données sur les Eaux Souterraines (ADES) database (http://www.ades.eaufrance.fr/). More than 500 water table observations were selected, corresponding to an unconfined aquifer, with a continuous time series lasting at least 10 years, and not directly affected by pumping.
In addition, TRIP was compared to the 50-yr simulation of the SIM hydrological forecast system performed at 8-km resolution by Vidal et al. (2010b) in order to study the robustness of the proposed simple groundwater parameterization against a more detailed model.
Finally, the simulated TWS variations were compared to the GRACE estimates. GRACE provides monthly TWS variations in terms of anomalies (ΔTWS) based on highly accurate maps of the earth’s gravity fields at monthly intervals and at spatial length scales of about 300–400 km (Wahr et al. 2004; Swenson et al. 2003). The instrumentation and onboard instrument processing units are described in detail by Haines et al. (2003). GRACE data can be used in a hydrological application to estimate ΔTWS from basin (Crowley et al. 2006; Seo et al. 2006) to continental scale (Schmidt et al. 2006; Tapley et al. 2004), as well as groundwater storage variations (Yeh et al. 2006; Rodell et al. 2009). Moreover, recent studies have used GRACE in order to evaluate their groundwater schemes (Niu et al. 2007; Lo et al. 2010; Decharme et al. 2010). Here, we used 85 months (from July 2003 to July 2010) of the Release 04 data produced by the Center for Space Research (CSR at The University of Texas at Austin) averaged over the aquifers defined over France in TRIP in order to evaluate the model’s capability to simulate the TWS variability at high and low resolution (Swenson and Wahr 2006).
a. River discharges
The simulated river discharges with and without groundwater were first compared to gauge measurements. This evaluation was made with the help of popular skill scores used in hydrology: the annual discharge ratio criterion (Ratio = Qsim/Qobs), the RMSE, and the efficiency (Eff) criterion (Nash and Sutcliffe 1970), which measures the ability of the model to capture the daily discharge dynamics (Nash and Sutcliffe 1970). This last skill score is defined as follows:
where represents the observed temporal mean and Eff can be negative if the simulated discharge is very poor, and is above 0.5 for a reasonable simulation.
Figure 5a shows the spatial distribution of the daily discharge efficiencies in the control experiment (NOGW12) while Fig. 5c compares these scores with (GW12) and without (NOGW12) a groundwater scheme for each river gauge located in the TRIP aquifer domains (gray-colored zone in Fig. 5c). The NOGW12 efficiency scores underline some weaknesses in the unmodified TRIP model. The scores are indeed negative in most of the Seine basin, in the north of the Loire River basin, in the Adour basin, and in the western part of the Alps. Adding a groundwater scheme improves the efficiency scores at 62% of the gauging stations (Eff difference greater than 0.05 in Fig. 5c) located inside the aquifer mask, while only 15% are deteriorated (Eff difference lower than −0.05). The groundwater scheme improves the simulation over the Aquitaine basin, the Rhone valley, and especially over the Paris basin, in which the widespread chalk aquifer is well connected to the river. Nevertheless, some weaknesses appear in the east of this basin and concern the Meuse River and also some tributaries of the Seine River such as the Aisne or the Marne.
Figures 5b and 5d show the comparison between TRIP and SIM. SIM takes two groundwater domains into account: the multilayer aquifer system beneath the Seine basin and the alluvial plain of the Rhone River. Not surprisingly, SIM outperforms NOGW12 for all river gauges located in the Seine aquifer (Fig. 5b). However, the introduction of the groundwater processes in TRIP (GW12) tends to reduce the gap with SIM (Fig. 5d). The GW12 scores can even exceed those of SIM over the Loire, the Somme basin, or north of the Garonne, pointing out the importance of groundwater processes in these regions. Both models give similar results in Brittany and in the Massif Central where no groundwater processes are simulated. Finally, the comparison of cumulative distributions of the efficiency scores, in which the GW12 curve is superior to NOGW12 and quasi-similar to SIM (Fig. 5e), confirms these results.
Moving to the discharge analysis at low resolution (NOGW05 and GW05; Fig. 6), the same kinds of results are found. Among the 61 gauging stations shown in Fig. 6c, 44% are better simulated by introducing the groundwater scheme, and only 16% are deteriorated. Note that the scores are not modified over the Rhone basin. Improvements in GW05 are, however, less pronounced than in GW12 (Fig. 6e). Even though the GW05 curve is still above the NOGW05, SIM, with a finer resolution (8 km) and a fine-tuned calibration of parameters, logically gives better results than GW05. Nevertheless, the efficiency differences shown in Fig. 6d confirm that, even at low resolution, accounting for groundwater processes has a positive impact on TRIP scores compared to the fine-tuned SIM simulation over the Loire and the north of the Aquitaine basin.
Figure 7 compares observed and simulated daily river flows at the outlets of the four main watersheds of France. Locations of the four gauging stations are circled in Fig. 5a. Observed river flows are satisfactorily reproduced by NOGW12 at the Loire, the Garonne, and the Rhone outlets, while the summer base flow is underestimated at the Seine outlet. The groundwater scheme improves scores over the Seine by increasing the efficiency from 0.04 to 0.82 (Table 2), and over the Loire (0.74 to 0.93). The Rhone and the Garonne river discharge scores are not deteriorated.
For all simulations, the mean annual cycle of daily river flows is also presented in Fig. 7. The scores of the low-resolution simulations are close to those obtained at high resolution. Only the Garonne discharge is impacted by the resolution. A spurious peak for NOGW05 and GW05 is observed in spring, which is related to unrealistic runoff and drainage from a few cells in the southwest of the Garonne watershed because of the aggregation of the forcing fields at low resolution. These cells encompass a part of the Spanish side of the Pyrenees not belonging to the Garonne drainage area. Elsewhere, GW12 and GW05 compare well with SIM, except over the Loire outlet, where TRIP obtains better results because of the fact that it simulates the existing groundwater explicitly (although with coarse parameters) while there is no regional modeling of the Loire aquifer included in SIM yet.
b. Piezometric head
It is rather difficult to compare local observations of the piezometric head with a simulation at a resolution larger than some kilometers, since the model cannot reproduce the impact of the local topography or the finescale variations of the geology. However, the temporal evolution of the piezometric head can be expected to be captured by the model, allowing an assessment of the model consistency at least for GW12. Two statistical criteria were used to quantify the ability of the model to reproduce the observations: the correlation and the RMSE. Figures 8a,c show GW12 correlation and RMSE scores for each selected piezometer in terms of monthly evolution of the head. More than 70% of them have a correlation greater than 0.5, and 68% have a RMSE lower than 2 m. A comparison with the simulated piezometric heads from SIM is shown in Figs. 8b,d in terms of correlation and RMSE, respectively, and only for the piezometers located in the SIM aquifer domains. In terms of correlation, the time evolution of the piezometric head is simulated better by GW12 on 27% of the selected wells and with a similar correlation on 25% of the measurements. Regarding the RMSE scores (without mean bias), 55% are in favor of GW12 and the two models are comparable for 5% of the scores. No clear pattern appears in the spatial distribution of these scores, confirming the difficulty of evaluating the water table simulation against in situ piezometric data. Nevertheless, these results suggest that the proposed simple groundwater model represents the time evolution of the water table as well as a more sophisticated, calibrated hydrogeological model used for operational regional applications.
Figure 9 compares the observed and simulated monthly evolution of the piezometric heads for the six piezometric gauges circled in Fig. 8a. These gauges were selected in order to represent different geological domains. Four piezometers belong to the SIM aquifer domains. The shape of the observed piezometric heads is well captured by the model. Two of them are located in the Seine basin and correspond to two kinds of geology: chalk and limestone. In both cases, GW12 is able to reproduce the observed interannual variability with good accuracy. Nevertheless, GW12 overestimates the amplitude of the annual cycle, especially in the limestone where the fine-tuned SIM simulation is closer to the observations. This result will be discussed in section 5. Two other piezometers located in the alluvial aquifer of the Rhone River valley and in the Saône watershed are shown in Fig. 9. The amplitude of the observed annual cycle is well captured by GW12 for the first one (RMSE = 0.61) but underestimated for the second one (RMSE = 4.6). In this basin, SIM was calibrated using the river flow only, and the piezometric heads are poorly reproduced, showing the need to use both types of observations in constraining hydrogeological models.
As already mentioned for the comparison with GW12, it is difficult to evaluate the simulated piezometric heads against local water table observations here, considering the coarse resolution of GW05 (about 50 km). Nevertheless, the comparison between GW12 and GW05 (Fig. 10) averaged over each large aquifer defined in TRIP (Figs. 2c,d) reveals that our groundwater scheme is not very sensitive to horizontal resolution. The two simulations give similar results on the 40-yr simulation. More details can be seen in the monthly mean seasonal cycles. Both the high- and low-resolution cycles are very close even if the GW05 amplitude is larger than GW12 over the Paris and the Aquitaine basins. The main differences appear over the upper Rhine plain, where a shift between the two curves is observed, together with a weaker amplitude of variations for GW05. This is due to more intensive exchanges between river and groundwater, which tend to smooth the low-resolution water table variations. Other reasons could be the low-resolution geometry of the aquifer, which encompasses only four grid cells with some subgrid areas not included in the high-resolution domain. Despite this problem, the impact of the resolution is generally weak in terms of water table variability, at least when a domain average is considered. Even though the low resolution induces some additional uncertainties in parameter estimation (porosity and transmissivity), these results show that the supposed weaker precision of these parameters compared to the high resolution is nevertheless sufficient to capture the first-order variations of the water table.
Figure 11 presents a detailed analysis of the groundwater scheme behavior in terms of simulated water table depth and groundwater–river exchange budget [Eq. (4)]. Results are shown for the Seine and the Rhone aquifers simulated in SIM (Fig. 5b), which is considered as a benchmark here. The monthly piezometric head averaged over each domain and the groundwater–river exchange budget are plotted for each aquifer. A 12-month running average is also plotted to focus on interannual variations. Over the Seine aquifer, GW12 captures the same water table interannual variations as simulated by SIM, although the amplitude is generally less pronounced. This amplitude seems to be resolution dependent since it is less pronounced for GW05 than for GW12. The associated mean annual cycles reveal that the seasonal variations compare well with the one simulated by SIM, regardless of the resolution. In terms of groundwater–river exchange variations, the GW12 signals are well correlated with SIM, while the GW05 interannual variability is slightly more pronounced. At both resolutions, the TRIP seasonal variations differ from those simulated by SIM. The TRIP seasonal maximum occurs in March while it occurs in February in SIM. This fact can be related to the same difference as observed on the mean annual cycle of the simulated discharge shown in Fig. 7 at the Seine outlet. During summer and autumn, the TRIP groundwater–river exchanges are less pronounced than in SIM. This is due to the fact that TRIP allows the river to feed the groundwater reservoir [negative Qriv in Eq. (4)], while such a process is neglected by SIM in this basin.
The comparison with SIM over the Rhone is less obvious because of the deficiency of SIM over this domain characterized by larger amplitude of the water table variations compared to the observations (see Fig. 9). Both interannual and seasonal variations of the simulated water table over the Rhone aquifer are less dependent on the resolution. The RMSE of the interannual variations between GW05 and GW12 is equal to 0.91 for the Seine, whereas it is reduced to 0.34 for the Rhone aquifer. This is also the case in terms of the groundwater–river exchange budgets for both the interannual and the seasonal variations. Note that the y axis of the mean annual cycle over the Seine and the Rhone basins are quite different.
c. Comparison with GRACE
The evaluation of the monthly TWS variations (ΔTWS) simulated by ISBA–TRIP against GRACE data is shown in Fig. 12. The simulated ΔTWS is computed as the sum of all land surface components, from ISBA, and the river storage (ΔS) and water table (ΔH) variations, from TRIP. The land surface components in ISBA are soil moisture (ΔWG), vegetation interception (ΔWR), and snow water equivalent (ΔSWE). Since GRACE data are anomalies relative to a reference geoid, each component is also calculated in terms of anomaly over the GRACE period (2003–10). The comparison is made in water equivalent height (cm) and only over the TRIP groundwater domains (Figs. 2c,d).
The ΔTWS simulated with and without groundwater are plotted for each resolution. The mean annual cycle of each component is also shown. Table 3 provides the contribution of each land surface component to ΔTWS. The soil moisture component is the main contribution (~90%) when groundwater is not considered, but it decreases to around 68% when the groundwater reservoir is taken into account. The groundwater contribution to the ΔTWS signal is significant at both high (~26%) and low (~30%) resolution. Finally, the introduction of the groundwater processes allows TRIP to simulate ΔTWS with better accuracy.
d. Sensitivity tests
The parameter τ regulates the groundwater–river exchanges. As defined previously, τ varies from 5 to 30 days from upstream to downstream grid cells within the river basin. To assess the sensitivity of the groundwater scheme to this empirical coefficient, two additional experiments were performed at low and high resolution using τ ± 75%, thereby keeping the values within a reasonable range.
Figures 13 and 14 show the mean annual cycles of the simulated and observed daily discharges at different gauging stations at high and low resolution, respectively. In each figure, four simulations are presented: NOGW12 (or NOGW05) and GW12 (or GW05), as well as both sensitivity experiments with τ ± 75%. The spatial distribution of the efficiency difference between these experiments and GW12 at each gauging station is also shown in Fig. 13.
The time transfer coefficient τ shows a limited impact on the discharge scores simulated at the outlets of the four main rivers of France whatever the model resolution. However, some regional impacts can be observed on the spatial distribution of efficiency scores. To illustrate these findings, three other gauging stations are presented at the high resolution (Fig. 13) to evaluate the groundwater scheme behavior over small upstream subbasins that are more difficult to analyze at the low resolution. Two subbasins, the Avre and the Aisne Rivers, are located in the Paris basin while the third, the Isle River, is in the Aquitaine basin (Fig. 13).
The discharge scores are improved by using τ − 75% over the Aisne River and more generally in the eastern part of the Paris basins. A weaker τ favors groundwater–river exchanges during autumn and winter to the detriment of groundwater water storage and then decreases the water flux to the river during spring and summer compared to GW12. As already mentioned, this region is mainly composed of siliciclastic rocks, characterized by mixed confined and unconfined aquifers, so the river flow is less impacted by groundwater processes. This could explain why a weaker τ improves the simulated discharge scores over this region, and underlines some limitations of the proposed definition of the aquifer basin boundaries.
Using τ + 75% improves the discharge scores over the Avre River and, to lesser extent, over the Isle River. A larger τ favors groundwater storage during the rainy season to the detriment of groundwater–river exchanges and thus more water remains available for the summer base flow. So, the river discharge decreases in autumn and winter, and increases during the dry season compared to GW12. This process is obvious over the Avre subbasin where the amplitude of the seasonal streamflow variations is generally weak. This result underlines another limitation of our simple groundwater scheme regarding parameterization of the groundwater–river flow.
Nevertheless, the uncertainties related to the empirical coefficient τ remain limited compared to the clear improvement of the river discharge scores provided by the representation of groundwater processes in TRIP.
Figure 15 presents the same detailed analysis of the groundwater scheme behavior in terms of simulated water table depth and groundwater–river exchange budget as in Fig. 11. Only the running averages are drawn for GW12, GW05, SIM, and the four sensitivity experiments. Using τ − 75% favors the groundwater–river water exchanges and so decreases both interannual and seasonal variations of the water table, and conversely for τ + 75%. This result confirms the previous discharge analysis and also reveals that the sensitivity to τ is more pronounced for water table variations than for river discharges.
In this study, the hourly surface runoff and deep drainage produced by the ISBA land surface model driven by the SAFRAN meteorological analysis have been used to feed the TRIP river routing model in order to simulate the daily discharges of the main French river basins at high and low resolution. Besides control experiments based on Decharme et al.’s (2010) version of TRIP, parallel simulations have been performed in which a simple groundwater scheme has been coupled with TRIP in order to represent the groundwater dynamics and its interaction with the river flow. After spinup, this new model leads to an improved simulation of daily river discharges over the 1970–2010 period, whatever the model resolution used. Indeed, results obtained at low resolution confirm the robustness of the model and its suitability for regional or even global applications.
As expected, the main improvement is found over the Seine watershed. This is due to the presence of a chalk aquifer over a large area, which is connected to the rivers in the valleys and so delays the base flow contribution to the river discharges. Groundwater acts as a buffer reservoir over this region and consequently the base flows simulated during summer are more realistic (Rousset et al. 2004; van den Hurk et al. 2005; Fan et al. 2007; Alkama et al. 2010). The river discharges simulated over the Rhone are also slightly improved, but the base flow remains underestimated during summer. The main reason is that the model does not consider the numerous dams used for hydroelectric power in the Alps, which tend to sustain the summer flow (Habets et al. 2008).
As mentioned in Habets et al. (2008), aquifers are widespread in France and the Seine and the Rhone basins are not the only ones that benefit from the introduction of a groundwater reservoir. The new model also improves the discharge scores over regions where groundwater processes are not yet represented in the operational SIM model: the Somme basin, the northern part of the Aquitaine basin, and especially the Loire basin. It is interesting to note that the positive impact found over these regions compared to SIM is also clear at low resolution (Fig. 6d).
The comparison between simulated and observed piezometric heads indicates that the simple methodology proposed in this study is sufficient to capture the interannual and seasonal variations of the water table over France as a whole. In contrast to the simulated river discharges, however, water table variations seem to be more sensitive to horizontal resolution. The main reason for this is that the hydraulic gradients between adjacent cells [Eq. (3)] are lower at 0.5° resolution than at resolution. This lower diffusive system favors the seasonal variations of the water table to the detriment of the long-term variability (Figs. 10 and 11). In a recent study, Lam et al. (2011) found that the lateral flows could be neglected at the resolutions and time scales of climate models. However, in this study, the spatial average of the absolute value of the lateral groundwater flow balance is 0.23 m3 s−1 for GW12 and 0.83 m3 s−1 for GW05. Even though the groundwater diffusion is lower at 0.5° because of the largest size of the cells, the mass of water transferred between the grid cells appears nonnegligible. Therefore, neglecting lateral groundwater flow remains questionable. Nevertheless, these findings need to be confirmed in future global hydrological and/or climate applications since the French domain is small compared to the largest river basins of the world.
The comparison between simulated and observed piezometric head variations in Fig. 9 also suggests some deficiencies, especially over the chalk aquifer of the Somme valley and the limestone aquifer of Beauce. In situ observations show a weak annual cycle while TRIP simulates larger seasonal variations. For at least two reasons, this can be related to the fact that each grid cell in TRIP is considered as a river cell. First, the unconfined aquifers over these regions are connected to a river network that is sparser than in most of the other aquifer basins of France (Martin 2000). Secondly, the piezometric gauges may be located relatively far from the river network. It is well known that, the nearer the observed piezometers are to the river, the larger are their seasonal variations because the water table variations are constrained by the seasonal variations of the river water height.
In addition, several studies have pointed out the impact of the unsaturated zone, which can be relatively deep in local plateaus, on the groundwater dynamics through the delay it introduces in the transfer of water between surface and groundwater (Pinault et al. 2005; Price et al. 2000; Habets et al. 2010). Such a process, not represented in TRIP, could explain some difference between observed and simulated piezometric head as shown over the Somme basin by Habets et al. (2010). Finally, the simple estimation of specific yield could also play a nonnegligible role. Nevertheless, over the alluvial plains of the Rhone and the Rhine, the RMSE scores between the observed and simulated piezometric head evolution highlight the suitability of the chosen parameterization to represent water table variations as well as groundwater–river exchanges. The piezometric level is constrained by the river, leading to a more realistic amplitude of the water table variations (Fig. 9).
These groundwater–river exchanges are controlled by the river conductance. Its simple parameterization using the τ parameter is obviously very empirical. In addition, since each TRIP grid cell is considered as a river cell, the system behavior is more sensitive to this parameter than to the geological characteristics of the aquifer formations. Our results show, however, that the choice of the τ value has a limited impact on the model results relative to the general improvement because of the representation of the groundwater processes in TRIP. Globally speaking, τ could be tuned so as to improve the discharge scores over some regions where discrepancies between the simulations and the observations are large. For global hydrological applications, such tuning is certainly not advisable for at least two reasons. First, this strategy would not be efficient over many regions of the globe where in situ observations are too sparse to usefully constrain the model. Secondly, this tuning could also compensate for systematic biases in the land surface model and/or uncertainties in the prescribed atmospheric forcings (especially precipitation). Another way to limit the impact and the tuning of the τ parameter could be to explicitly represent the unsaturated zone in TRIP.
The methodology used to define the geometry of the aquifers also shows some limitations by embracing the entire basin without distinction between confined and unconfined aquifer zones. This is especially true over the Aquitaine basin where the BDRHF hydrogeological database indicates that shallow, unconfined aquifers are only defined along the rivers while confined aquifers are present elsewhere. This could explain the score deteriorations observed over some gauging stations in the south of the Aquitaine basin. The same limitation appears in the eastern part of the Paris basin (as mentioned in section 4d) where the standard TRIP model simulates more realistic river discharges than GW12 or GW05. This result is confirmed by the comparison with the fine-tuned SIM simulation (Figs. 5 and 6).
However, the proposed simple definition of the aquifer geometry appears acceptable, especially in the perspective of global applications. Thanks to the WHYMAP hydrogeological resources map, it allows the major unconfined aquifers to be distinguished from the shallow, local aquifers located over old platforms composed of rocks with low permeability, such as metamorphic or crystalline rocks. This observation is confirmed by the good results obtained with NOGW12 over mountains and old basement massifs such as the Massif Central and Brittany. Such a distinction was also used with success at fine resolution by Fan et al. (2007), who considered two kinds of geological domains in order to calibrate their model: bedrock and regolith.
Finally, the aquifer parameters have been estimated using a simple methodology potentially applicable at the global scale given the required inputs. For example, Dürr et al. (2005) compiled a global lithological map and, more recently, Gleeson et al. (2011) proposed a map of permeability covering the whole surface of the earth.
While the direct evaluation of low-resolution simulations against local water table observations remains questionable because of the heterogeneity of the geology and the topography inside the cell, the successful comparison with GRACE TWS estimates confirms the suitability of this satellite product for evaluating water table variations at the global scale, which is in line with the studies by Niu et al. (2007), Lo et al. (2008), and Decharme et al. (2010).
This study describes a simple groundwater scheme that has been coupled to the TRIP river routing model for global hydrological and climate applications. Here, it is evaluated over France using offline simulations over the 1970–2010 period at high () and low (0.5°) resolution. This scheme accounts for the groundwater–river exchanges and the water table dynamics explicitly. The simulated river discharges are evaluated against 318 in situ gauging measurements distributed all over France, while the simulated water table heads are compared to more than 500 piezometers. The simulated ΔTWS are also compared to the GRACE satellite-derived ΔTWS estimates.
Without a groundwater scheme, and whatever the horizontal resolution used, the standard TRIP model shows significant deficiencies in simulating daily discharges compared to observations, with overestimated river flows during winter and underestimated base flows during summer. These problems are mainly due to the absence of an explicit coupling between the groundwater and the river and are partly overcome by the new model, which is able to capture the spatiotemporal variability of both river discharges (especially the summer base flow) and water table variations over the main French river and aquifer basins. These positive results are obtained by using a simple methodology based on geological data available at the global scale.
Nevertheless, some deficiencies appear throughout the model evaluation, for which there are several possible explanations: uncertainties in the SAFRAN meteorological forcing, errors in the ISBA land surface model, possible anthropogenic influences on the observed discharges and piezometric heads, but also obvious deficiencies in the modified TRIP hydrogeological model. Besides uncertainties in the standard TRIP parameters (river geometry and slope), some aspects of the groundwater scheme are obviously questionable, such as the simple method to delimit the aquifer boundaries or estimate the values of certain parameters, such as the transmissivities or the porosities. The groundwater–river exchange coefficient is potentially tunable region by region. However, such tuning is not yet an option for at least two reasons. First, it might be sensitive to the experimental conditions, especially the quality of the prescribed atmospheric forcing or the land surface parameterization. Second, it could, to some extent, compensate for the absence of hydrological processes, such as floods or water capillary rises in the soil, in ISBA.
While the model skill could be improved through a calibration of these parameters, the reasonable scores obtained in the present study and their limited sensitivity to the empirical parameter, τ, suggest that the groundwater scheme is robust and is suitable for global hydrological and climate applications. The positive impacts found on the simulated river discharges over France will have also to be confirmed over the other large river basins of the world. As already suggested, TRIP will be explicitly coupled with the ISBA land surface model in order to simulate the interaction between the deep water table dynamics and the overlying unsaturated soil. The goal will be to represent the impact of water capillary rise on the land surface energy and water budgets. The ultimate objective will be to introduce this new land surface component into the CNRM global climate model (Voldoire et al. 2012) in order to assess the relevance of groundwater processes for the simulation of both recent and future climates.
The authors wish to thank the BGR and UNESCO for providing the WHYMAP and IGME data, and the BRGM for making the simplified lithological map of France available. Discharge observations were provided by the French Hydro database (Ministère de l’écologie; http://www.eaufrance.fr) and piezometric head observations by the ADES database (Ministère de l’écologie; http://www.ades.eaufrance.fr). GRACE land data were processed by Sean Swenson, supported by the NASA MEASURES Program, and are available at http://grace.jpl.nasa.gov. This work is supported by the Centre National d’Études Spatiales (CNES), the Centre National de Recherches Météorologiques (CNRM) of Météo-France, and the Centre National de la Recherche Scientifique (CNRS) of the French Research Ministry. The authors thank Emmanuel Ledoux (MINES Paris Tech) for providing precious advice on this work.
Discretized Groundwater Diffusive Equation
The implicit finite-difference numerical scheme and the Gauss–Seidel iterative method used in MODCOU (Ledoux et al. 1989) solve the discretized spherical form of the diffusive groundwater equation [Eq. (3)]. This discretized form can be obtained by following several algebraic steps.
Figure A1 shows five neighboring cells of the latitude–longitude domain. The term corresponding to the lateral flows along the latitude axis is discretized as follows:
where Qec and Qwc are the approximations of between the nodes e and c, and w and c, respectively (Fig. A1):
The coefficients Tec and Twc are equal to the geometric mean of the transmissivity values of two adjacent neighboring cells, while ϕec and ϕwc are the latitude between the nodes e and c, and w and c, respectively (Fig. A1). Similarly for the longitude axis, we have
Since , after several algebraic steps, we obtain the following discretized form of Eq. (3) expressed in m3 s−1:
where is the grid cell area and Qsb and Qriv are the deep drainage from ISBA and the groundwater–river exchanges, respectively.