Spatial variation of groundwater response to multiple drivers in a depleting alluvial aquifer system,northwestern India

Abstract

Unsustainable exploitation of groundwater in northwestern India has led to extreme but spatially variable depletion of the alluvial aquifer system in the region. Mitigation and management of groundwater resources require an understanding of the drivers behind the pattern and magnitude of groundwater depletion, but a regional perspective on these drivers has been lacking. The objectives of this study are to (1) understand the extent to which the observed pattern of groundwater level change can be explained by the drivers of precipitation, potential evapotranspiration, abstraction, and canal irrigation, and (2) understand how the impacts of these drivers may vary depending on the underlying geological heterogeneity of the system. We used a transfer function-noise (TFN) time series approach to quantify the effect of the various driver components in the period 1974–2010, based on predefined impulse response functions (θ). The dynamic response to abstraction, summarized by the zeroth moment of the response M₀, is spatially variable but is generally large across the proximal and middle parts of the study area, particularly where abstraction is high but alluvial aquifer bodies are less abundant. In contrast, the precipitation response is rapid and fairly uniform across the study area. At larger distances from the Himalayan front, observed groundwater level rise can be explained predominantly by canal irrigation. We conclude that the geological heterogeneity of the aquifer system, which is imposed by the geomorphic setting, affects the response of the aquifer system to the imposed drivers. This heterogeneity thus provides a useful framework that can guide mitigation efforts; for example, efforts to decrease abstraction rates should be focused on areas with thinner and less abundant aquifer bodies.

Keywords

Transfer function-noise model impulse response function groundwater level fall northwestern India groundwater abstraction canal irrigation

I Introduction

In regions with large aquifer systems that undergo frequent water stress, groundwater is often used as an additional water source. If groundwater abstraction exceeds the natural groundwater recharge over an area for long periods of time, over-exploitation or persistent groundwater depletion occurs (Wada et al., 2010). The unsustainable exploitation of groundwater resources is now a very significant problem globally and requires urgent attention (Aeschbach-Hertig and Gleeson, 2012; Gleeson et al., 2010). As Famiglietti (2014) describes, whilst groundwater is a critically important global water resource, it is given insufficient attention within management systems compared to visible surface water resources. This is particularly true in countries where water governance is weak or absent and monitoring of aquifers is inadequate. Consequently, for many heavily exploited aquifers there is a lack of knowledge about (1) how groundwater storage responds to various drivers, such as precipitation, evapotranspiration, abstraction and irrigation, (2) how storage variations relate to aquifer heterogeneity, and (3) how future changes in groundwater levels might be anticipated and mitigated on the basis of these various drivers. Such knowledge is critical for management, especially if the heterogeneity of the aquifer system leads to substantially different responses to future stresses in different parts of a region.

The Indo-Gangetic foreland basin is one of the world’s most prominent hotspots of groundwater depletion, especially in northwestern India (Chen et al., 2014; Kumar et al., 2006; MacDonald et al., 2016; Richey et al., 2015; Rodell et al., 2009; Shah, 2009; Tiwari et al., 2009). This has been caused by the increased use of groundwater for irrigation since the mid-1960s as part of the “Green Revolution,” the popular name for the agricultural strategy that aimed to make India self-sufficient in food grain production. Groundwater abstraction for irrigation has become a particularly severe issue in the states of Punjab and Haryana, whose contribution to national food grain production increased from 3% before the Green Revolution to approximately 20% at the end of the 20th century (Singh, 2000). Rapid groundwater level decline associated with groundwater pumping has been documented across these states by MacDonald et al. (2016). However, in parts of the region, significant groundwater level rise has also occurred due to infiltration from canal irrigation return flow and canal seepage (Joshi et al., 2018). Recent work by Asoka et al. (2017) showed that groundwater storage changes in northwestern India between 2002 and 2013 can be explained by both abstraction and precipitation, although the former appears to have been somewhat more important. They considered only a single time series of storage change, however, and did not assess spatial variations in either groundwater level changes or in the potential drivers of those changes.

Here we address the urgent societal issue of groundwater depletion in northwestern India by applying a time series analysis to understand the spatial variation in groundwater level change, focusing in particular on the area between the Sutlej and Yamuna Rivers where historical decline has been greatest. The objectives of this study are to: (1) assess whether a time series model using parsimonious, physically-related impulse-response functions can reproduce changes in the spatial pattern of groundwater levels since 1974; (2) understand the extent to which the observed pattern of groundwater level change can be explained by the drivers of precipitation, potential evapotranspiration (PET), abstraction, and canal irrigation; and (3) investigate the extent to which the impacts of these drivers depend on the geological heterogeneity of the aquifer system. First, we present the available climatic and groundwater level data for the region, and describe the time series model and its implementation. We examine the influence of individual driver components on the aquifer system over time, and then use a single measure of the system response to quantify spatial variations in response across the study area. Finally, we relate our results to the geological framework of the study area, and explore the potential implications of the results for groundwater management and their application to other depleting aquifer systems.

II Study area

This study focuses on the northwestern region of India, which is part of the Indo-Gangetic foreland basin and is fed by the Sutlej River in the west and the Yamuna River in the east of the study area (Figure 1(a)). Recent studies have identified sediment deposits in the study area that were sourced from the Yamuna and Sutlej catchments (Clift et al., 2012; Singh et al., 2017), and geophysical profiles have verified the existence of large paleochannels within the subsurface (Khan and Sinha, 2019; Sinha et al., 2013). The alluvial aquifers in this study area are formed by sediments eroded from the Himalaya and redistributed by the Sutlej and Yamuna rivers, forming two major sedimentary fan systems (Geddes, 1960; Singh et al., 2016; Van Dijk et al., 2016a). The distribution of the channel sand deposits that form the primary aquifer bodies within these fan systems is controlled by river avulsion, sedimentation rate, and the stacking pattern of fluvial channel-belt units (Allen, 1978, 1984; Bridge, 1993; Bridge and Leeder, 1979; Heller and Paola, 1996; Holbrook, 2001; Leeder, 1978; Sheets et al., 2002; Straub et al., 2009; Van Dijk et al., 2016b). This alluvial stratigraphy, in turn, determines the characteristics and productivity of the aquifer (Anderson, 1989; Fogg, 1986; Weissmann et al., 1999), in terms of (1) the percentage of sand-rich aquifer bodies in the subsurface; (2) the geometry and dimensions of the aquifer bodies; (3) their hydraulic conductivity and specific yield; and (4) their vertical and horizontal connectivity (Flood and Hampson, 2015; Larue and Hovadik, 2006; Renard and Allard, 2013).

Figure 1.

(a) Geomorphological map and the subsurface aquifer percentage of the study area in northwestern India, modified after Van Dijk et al. (2016a). The study area covers the Sutlej and Yamuna fans (light green) and the interfan area between them (pink). Dark green areas show the incised valleys of the modern Sutlej and Yamuna rivers and the Ghaggar–Hakra paleochannel (Singh et al., 2017). Light blue lines show major canals. Dots show bulk aquifer percentage in the upper 200 m of the subsurface, based on CGWB aquifer-thickness logs (see Van Dijk et al., 2016a for description). Long-dashed line shows the medial transect in (b), while the short-dashed line shows the distal transect in (c). The study area is outlined by the thick black line and divided into 10 x 10 km grid cells. (b) Medial transect illustrating aquifer (yellow) and non-aquifer (green) units in the subsurface. Note abundant stacked aquifer bodies. (c) Distal transect showing the lower abundance of aquifer bodies compared to the medial part of the fan. Transects in panels (b) and (c) modified from Van Dijk et al. (2016a).

Van Dijk et al. (2016a) established that aquifer body thickness and percentage vary in systematic and predictable ways between proximal and distal parts of the fan systems, and between the fans and the interfan or marginal areas between them (Figure 1(a)). The elongated channel deposits that form aquifer bodies are highly connected in the down-fan direction but are less connected in the lateral direction. The bulk percentage of aquifer bodies decreases in the down-fan direction (Figure 1(b) and (c)), although the distribution of aquifer-body thickness remains the same (Van Dijk et al., 2016a). The geomorphic distinction between fan and interfan settings within the Indo-Gangetic basin also introduces an important large scale lateral heterogeneity. Aquifer bodies are generally thinner and less abundant for the interfan area, whereas the fans consist of abundant stacked channel sand bodies as result of the formation and subsequent filling of incised valleys across the fan surface (Van Dijk et al., 2016a). An example of subsequent filling of an incised valley is the Ghaggar–Hakra paleochannel (Singh et al., 2017), which is at the border of the Sutlej and Yamuna fan system. Thus, the observations provide evidence that the alluvial aquifers in northwestern India are highly spatially heterogeneous, and that the surface geomorphology provides a clear guide to the subsurface architecture of the alluvial aquifer system. Van Dijk et al. (2016b) further argued that the geological framework imposed by these fan systems and the interfan areas should be used as a basis for understanding and relating aquifer properties, groundwater level changes, and potential groundwater management approaches.

3 Methods

1 Methodology

Here, we focus on trying to understand the relative role of different drivers, imposed on the heterogeneous regional framework mentioned above, in explaining the rates and patterns of historical groundwater level change. This requires a groundwater modeling approach. A diverse range of models have been applied to simulate groundwater dynamics, and the appropriate degree of model complexity depends on the goal of the modeling, the amount of data with which to constrain the model, and practical project constraints (Guthke, 2017). Distributed groundwater models, which solve physical laws governing groundwater flow by discretizing the aquifer domain, remain the most widely used. This approach allows for complex heterogeneous and anisotropic fields of aquifer properties, but typically results in models with large numbers of parameters, for which a careful assessment of uncertainty is required (Hill and Tideman, 2007; Refsgaard et al., 2012; Rojas et al., 2010). Given the poor spatial resolution across northwestern India in both measurements of aquifer properties (UNDP, 1985; Van Dijk et al., 2016a) and water levels (MacDonald et al., 2016), it is difficult to justify the use of such a complex approach in this region.

Conceptual, lumped-parameter modeling is an alternative approach that simplifies the representation of processes incorporated in physically-based models, but maintains some fundamental physical principles from our conceptual understanding of groundwater systems (e.g. Kazumba et al., 2008; Mackay et al., 2014; Park and Parker, 2008). These models incorporate parameters that can be associated with measurements of properties made in the field, such as hydraulic conductivity or specific yield. For example, Mackay et al. (2014) applied a lumped-parameter groundwater model, AquiMod, to simulate groundwater levels in observation boreholes in different aquifers across the UK. The model was driven by rainfall and PET time-series, and contained three conceptual stores representing soil, an unsaturated zone, and a saturated aquifer, which were used to simulate infiltration to the groundwater table. As with other lumped-parameter groundwater models (Barrett and Charbeneau, 1997; Hong et al., 2017; Long, 2015; Pozdniakov and Shestakov, 1998), however, AquiMod typically required the specification of a large number of parameters, although many of these could be fixed based on prior information and expert judgment. Given the limited prior information for this region, there are no valid ways to constrain all of those parameters, and therefore we did not adopt this approach.

Another means of simulating groundwater level fluctuations is by the use of a statistical time series approach, such as transfer function-noise (TFN) models (Berendrecht et al., 2003; Von Asmuth et al., 2002), which are especially useful for modeling systems whose behavior cannot easily be described by physical processes or quantified by physically-observable parameters. TFN models have been widely used in hydrology and can be divided into three types: (a) models that start from a geo-statistical methodology, applying space-time kriging or co-kriging (Van Geer and Zuur 1997); (b) models that are based on multivariate time series analysis, where multiple time series are correlated in space (Von Asmuth et al. 2002); and (c) models that combine elements of methods (a) and (b) (Bierkens et al., 2001; Yuan et al., 2008). For example, Von Asmuth et al. (2008) used a time series analysis method to predict fluctuations in groundwater level from multiple drivers. They described a class of parsimonious TFN models that implement predefined impulse response functions in continuous time (PIRFICT), which circumvents a number of limitations of discrete TFN models linked to time discretization and model identification. The PIRFICT methodology has since been used in a number of studies to determine the effect of multiple drivers on groundwater levels in different regions of unconfined aquifers (Obergfell et al., 2013; Shapoori et al., 2015a). Application of the Von Asmuth et al. (2008) model is potentially useful for regions like northwestern India where (1) we lack a detailed physical understanding of the aquifer system and high-resolution subsurface data with which to constrain the key parameters (Van Dijk et al., 2016b), and (2) groundwater level is likely to be controlled by a small number of major driver components. This approach is ideal for understanding how groundwater responses vary spatially and how the response is linked to the underlying aquifer characteristics and geology. The model can approximate regional heterogeneity, but is not set up to deal with finer scale heterogeneity that is observed within individual channel bodies (Donselaar and Overeem, 2008; Holbrook, 2001; Miall, 1985; Van de Lageweg et al., 2016a, 2016b; Willis and Tang, 2010).

2 Data acquisition

For the model inputs we collected district-wise data on precipitation and PET from the Indian Meteorological Department (IMD) for the period 1951–2010. We obtained district-averaged abstraction values from the Central Groundwater Board (CGWB) for the years 2004, 2009, and 2011, and reconstructed the monthly groundwater abstraction for irrigation by the deficit between crop water requirements and effective precipitation. Groundwater level data were collected for the period 1974-2010 from borehole databases maintained by the state groundwater departments of Haryana and Punjab, and by the CGWB. Measurements were made twice yearly (pre- and post-monsoon) by the state groundwater departments, and four times yearly by the CGWB. For more detail on the data acquisition see the Supplemental Material.

Initial analysis of the climate and groundwater level data indicated a potential relation between groundwater level decline and total abstraction over the period of observations. This is not surprising (Asoka et al., 2017), although declines in groundwater level will be controlled by the amount by which abstraction exceeds groundwater recharge, rather than by the magnitude of abstraction itself. The analysis, and that by Asoka et al. (2017), does not illuminate the reasons for the spatial patterns in decline (Supplementary figure 3), however. In addition, the degree of scatter that is visible in Supplementary figure 4 suggests that abstraction alone is not an adequate predictor of groundwater level change. To investigate the response of groundwater levels to the combined effects of precipitation, PET, abstraction, and irrigation, we implemented a transfer function-noise time-series model.

3 Model setup

We adopted the formulation for multiple drivers presented by Von Asmuth et al. (2008) in the TFN method. This formulation relies on predefined impulse response functions for each driver in continuous time (Von Asmuth et al. 2002). The estimated groundwater level at time t in response to driver i, $h_{i} (t)$ , is given by

h_{i} (t) = \int_{- \infty}^{t} R_{i} (τ) θ_{i} (t - τ) d τ

where $R_{i}$ is the value of driver i at time t, $θ_{i}$ is the impulse response function of driver i, and τ is time after the impulse is applied. The groundwater level in response to the four drivers considered in this study, $h (t)$ , can then be written as

h (t) = h_{p} + h_{e} + h_{w} + h_{c} + d + m (t)

where $h (t)$ is the estimated level at time t, and is computed by the summing the contributions from precipitation (p), PET (e), groundwater abstraction (w), and canal irrigation (c), the local drainage level relative to a reference level d, and the residual series m, which is the difference between observed and modeled.

Equation (1) shows that impulse response functions must be defined for each of the drivers that we impose. Our impulse response function for precipitation, $θ_{p}$ , takes the form of a Pearson type III distribution function, as used by several previous studies (Von Asmuth et al., 2002, 2008) in a variety of hydrogeological settings:

θ_{p} (t) = A \frac{b^{n} t^{n - 1} e x p (- b t)}{Γ (n)}

where A, b, and n are parameters that define the shape of $θ_{p}$ , and

Γ (n) = \int_{0}^{\infty} e x p (- t) t^{n - 1} d t .

The physical interpretation of equation (3) is a series of coupled linear reservoirs where n is the number of reservoirs, b is the inverse of the reservoir constant normally used, and A adjusts the area of response (Von Asmuth et al., 2002). Von Asmuth et al. (2008) argued that PET should have a similar effect as precipitation on h, although it will have an opposite sign and a different magnitude. Consequently, water level in response to PET was modeled as

h_{e} (t) = \int_{- \infty}^{t} - e (τ) f θ_{p} (t - τ) d τ

where e is the PET time series and $θ_{p}$ is the response of the system to precipitation given by equation (3). The PET factor f accounts for a reduced dependence of h on e compared to p, and should depend on soil and land cover conditions that vary through time; for simplicity we assumed that it was constant.

Earlier PIRFICT-based TFN modeling studies that incorporated groundwater abstraction (Obergfell et al., 2013; Shapoori et al., 2015a, 2015b; Von Asmuth et al., 2008) implemented a three-parameter impulse response function based on the Hantush (1956) solution to the drawdown in a leaky confined aquifer. Following Shapoori et al. (2015a, 2015b), we used a two-parameter impulse response function based on the Theis solution for the drawdown in a confined aquifer of the form

θ_{w} (t) = - \frac{γ}{t} e x p (- \frac{α}{t})

where α and γ are calibration parameters. The parameter γ is equivalent to $\frac{1}{4 π T}$ with T (L² T⁻¹) being the transmissivity of the aquifer. The parameter α is defined by $\frac{r^{2} S}{4 T}$ , where S (dimensionless) is the aquifer storage coefficient and r represents the distance between a pumped borehole and the observation point. However, given that we were simulating the cumulative impact of pumping from numerous tube wells on the gridded, spatially-averaged groundwater level over a multi-decadal period, and did not seek to explicitly identify values for T and S, we did not define these variables (T, S and r) and only considered the integrated parameters α and γ.

The final driver to be incorporated into the TFN model represents the combined effect of canal irrigation return flow and canal seepage. We modeled this response function with a simple exponential function of the form

θ_{c} (t) = - \frac{1}{S_{y}} e x p (- \frac{t}{λ})

where $S_{y}$ is the specific yield of the aquifer and λ is a decay constant. The monthly irrigation driver c is represented by the irrigation rate, I, divided by $S_{y}$ . When λ is large the effect of the canal irrigation persists for a long period of time, so that the contribution of irrigation during any time-step to $h_{c}$ decays only slowly. Conversely, when a small value of λ is used, the contribution of irrigation to $h_{c}$ decays rapidly to zero, in other words ( $I / S_{y}$ ) $\cdot θ_{c}$ is virtually 0.

4 Model application

The TFN model was applied using a monthly time step to simulate the time series of cell-averaged groundwater levels independently for each of the 664 10 × 10 km cells across the study area. The drivers $R_{i}$ for precipitation, PET and groundwater abstraction were defined in units of mm/day from the measured or estimated data, as described in Section III.3 above, whereas the canal irrigation driver, $R_{i}^{c}$ , was unknown. Consequently, this was defined as an extra calibration parameter, with the assumption that it was constant over time. Rather than calibrating both $S_{y}$ and λ in equation (7), however, we set $S_{y}$ equal to one and therefore sought to calibrate the canal irrigation driver scaled by the specific yield, i.e. $R_{i}^{c} / S_{y}$ . This was considered reasonable given the resolution of our modeling and that specific yield is relatively uniform across our study area; median values have been estimated to be between 0.11 and 0.15 by CGWB (2011) and MacDonald et al. (2016).

Each model was run for the 60-year period 1951–2010, and calibrated against the 1974–2010 groundwater level time series. Therefore, the simulation period contains a 23-year spin up period, during which time the effect of pre-1951 memory in the impulse response functions is lost. The local drainage parameter, d, was fixed to a level defined by extrapolating the groundwater level data over the period 1974–1984 back to 1951 using linear regression. The model was calibrated using a Monte Carlo procedure, within which values for the eight parameters (A, b, n, f, γ, α, λ, $R_{i}^{c} / S_{y}$ ) were sampled from uniform distributions with pre-defined lower and upper limits (Table 1). Each model was run 150,000 times, and model error was calculated using the Nash–Sutcliffe efficiency ( $N S E$ ) (Nash and Sutcliffe, 1970), which is given as

N S E = 1 - \frac{\sum_{i = 1}^{n} {(h_{o} - h_{m})}^{2}}{\sum_{i = 1}^{n} {(h_{o} - {\bar{h}}_{m})}^{2}}

where $h_{o}$ and $h_{m}$ are the observed and modeled groundwater levels over n time steps. If $N S E$ = 1, the model is a perfect match to the observations. If $N S E$ = 0 the model is as accurate as the mean of the observed data, and if $N S E$ < 0, it is worse than the observed mean. A value of $N S E$ = 0.2 was set a priori as the threshold between behavioral and non-behavioral model simulations; in other words, those simulations with $N S E$ > 0.2 were deemed to have produced acceptable fits to the observed groundwater level time series. The $N S E$ does not measure how good a model is in absolute terms (Criss and Winston, 2008; McCuen et al., 2006; Schaefli and Gupta, 2007). Therefore for further analysis we looked at the parameter values for the outcome with the best $N S E$ value, but also to the median parameter values of all acceptable models.

Table 1.

Ranges of parameter values sampled during Monte Carlo procedure.

Impulse response function	$θ_{p}$			$θ_{e}$	$θ_{w}$		$θ_{c}$
Parameter	A	b	n	f	γ	α	λ	$R_{i}^{c} / S_{y}$
Upper limit	500	0.2	2.0	0.4	100	1.0	1000	20
Lower limit	50	0.05	0.1	0.05	5	0.05	10	0

Parameter values for the impulse response functions were chosen to encompass a wide range of shapes and scales (Figure 2). An additional condition was imposed on the ratio $\frac{n}{b}$ in equation (3), which determines the number of months it takes to reduce the precipitation contribution to the groundwater level $h_{p}$ by half after a precipitation event. If this was greater than 24 months, then n and b were resampled until the ratio was less than or equal to 24. Without this condition, the $θ_{p}$ response curve can persist over many years, enabling the model to generate an unjustified long-term upward trend in groundwater level even with a stationary precipitation time series.

Figure 2.

Example impulse response functions for (a) precipitation $θ_{p}$ (with A = 1, so no magnification), (b) abstraction $θ_{w}$ , and (c) canal irrigation $θ_{c}$ , illustrating the lower and upper limits of the sampled parameter value ranges (Table 1). Note that $θ_{e}$ is assumed to be identical in form to $θ_{p}$ .

5 Data analysis

A primary goal of this work is to assess not only the relative importance of different drivers in determining water level $h (t)$ , but also how that relative importance varies in space across the study area. In other words, we wish to understand whether groundwater levels in some areas are more sensitive to one or another of the different input drivers. To quantify the importance of driver i, we used its zeroth moment $M_{0, i}$ , also known as the stationary response, defined as the integral of the calibrated impulse response function $θ_{i}$ over time:

M_{0, i} = \int_{0}^{\infty} θ_{i} (t) d t

A large value of M₀ means that the driver has a large effect on the groundwater level, due to a response that is large in magnitude, persistent in time, or both. A small value of M₀ means that the effect on the groundwater level by the driver is minimal.

To characterize the spatial variation in the relative importance of precipitation, abstraction, and canal irrigation drivers, we plotted the zeroth moment for the best fit solution with the highest $N S E$ value. We also plotted the median zeroth moment from all of the acceptable solutions with $N S E > 0.2$ and the coefficient of variation, $C V$ , defined as the standard deviation (σ) divided by the mean (μ) for those solutions. The zeroth moment for PET stress was not determined because $θ_{e}$ was assumed to be a fraction of $θ_{p}$ , and so this will show the same spatial pattern. The $N S E$ values give the efficiency of the model outcomes, but we were also interested in the goodness of the model outcome. A disadvantage of the $N S E$ is that larger values in a time series are strongly overestimated whereas lower values are neglected, because of the use of squared difference values (Legates and McCabe, 1999). Because of limited seasonal and daily time series of the groundwater level data, the model fails to reproduce any smaller time scale fluctuations but can still report a good NSE value (Schaefli and Gupta, 2007). A more natural measure of average error of the model is given by the mean absolute error ( $M A E$ ), which is a more unambiguous measure of the difference between the modeled and observed groundwater level.

IV Results

1 Groundwater level series

The measured groundwater level changes show a general groundwater level decline in the northeastern part of the study area and water level rise in the southwestern part over the period 1974–2010 (Figure 3(a)). The calibrated TFN model reproduces this overall pattern (Figure 3(b)). The model is able to capture the spatial distribution, although not necessarily the observed values, in areas of large groundwater level change. In detail, the modeled declines in groundwater level are more patchy than the observed pattern, and the absolute declines are somewhat under-predicted, leading to relatively high $M A E$ values in these areas (Figure 3(d)). The model does poorly in areas with no groundwater level change and along the Ghaggar–Hakra paleochannel (Figure 3(a)), and these areas were excluded because the $N S E$ values were < 0.2 (indicated by the cross-hatch pattern in Figure 3).

Figure 3.

(a) Observed and (b) modeled groundwater level changes for the period 1974–2010. Black square boxes indicate grid cells that are elaborated in Figures 4 and 5. Cross-hatched cells yielded no model solutions with $N S E > 0.2$ . (c) The Nash–Sutcliffe error ( $N S E$ ) value for the best model outcomes in each grid cell. (d) The mean absolute error ( $M A E$ ) of the best model outcome in each grid cell. Overall the $M A E$ is less than 4 m. Black dots show the major cities in northwestern India given in Figure 1.

We illustrate the results of the time series model from three example locations with different patterns of temporal evolution (see Figure 3(b) for locations): groundwater level rise, decline in the area of the Ghaggar–Hakra paleochannel, where we expect thinner and less abundant aquifer bodies (Van Dijk et al., 2016a), and decline on the Sutlej fan, where we expect thick and more abundant aquifer bodies (Figure 4). At each location, the modeled groundwater level is decomposed into four partial series to show the response of the level to precipitation, PET, abstraction, and canal irrigation. For the location with groundwater level rise (Figure 4(a)), the effect of canal irrigation on the groundwater level $h_{c}$ is constantly increasing. The precipitation component $h_{p}$ shows strong seasonal variation but its absolute value is much less than $h_{c}$ . The PET component also varies seasonally and counteracts the precipitation component. The abstraction component $h_{w}$ is monotonically increasing due to increasing abstraction, but again its magnitude is significantly less than that of the canal irrigation component.

Figure 4.

Three examples of TFN model outcomes for grid cells with (a) groundwater level rise, (b) groundwater level decline in the area of the Ghaggar–Hakra paleochannel, and (c) decline in the area of the Sutlej fan. See Figure 3 for locations. In each column, the top panel shows observed (dots) and predicted (black line) water levels over time; the gray area indicates the model spin-up time from 1951 to 1974. In the lower panels, the red lines indicate the input time series of precipitation, PET, and abstraction. Black lines indicate the contribution of each component to the water level. Groundwater level rise in (a) can be explained by the monotonically-increasing canal irrigation component ( $h_{c}$ ), whereas the precipitation and abstraction components are much smaller and dominated by seasonal variations. Groundwater level decline in the paleochannel (b) is dominated by the abstraction component ( $h_{w}$ ). In comparison, decline on the Sutlej fan (c) shows a similar pattern, although the magnitude of the abstraction component is only about half of the value at the paleochannel site in (b). The canal irrigation component $h_{c}$ becomes stable after the spin-up time for the areas with groundwater level fall (b and c), indicating that it has no effect on the long-term groundwater level.

Groundwater level decline in the area of the Ghaggar–Hakra paleochannel (Figure 4(b)) is dominated by the abstraction component $h_{w}$ . The influence of precipitation is relatively constant after the spin-up time, and while there is some seasonal variation as well as yearly variation because of precipitation variation over the years, the system responds rapidly to these inputs. PET has an essentially constant impact after the spin-up time, as does the canal irrigation component $h_{c}$ , so these components effectively do not influence long-term groundwater levels. Similar behavior is observed for the area of groundwater level fall on the Sutlej fan (Figure 4(c)). There, abstraction remains the dominant influence on the long-term decline in water levels, although the influence of that component does not increase after 1990 as abstraction rates appear to have stabilized after that time.

2 Impulse response functions

The impulse response function describes the dynamic response of the groundwater level after a sudden input or change of a driver i. The calibration parameters determine the shape of this response in terms of its amplitude and duration. In our implementation of the TFN model, impulse response functions are generated for each grid cell independently based on the water-level history and stresses imposed in that cell. Here we illustrate the best-fit impulse response functions for the same three locations that were shown in Figure 4. The response to precipitation is rapid for the location experiencing groundwater level rise (Figure 5(a)), but is delayed by a few months for the locations experiencing decline (Figure 5(b) and (c)). These relatively short-term responses explain why $h_{p}$ fluctuates with both seasonal and inter-annual variations in precipitation (Figure 4). Recall that the model is run at a monthly time step, so individual storm events are not included in the calculation, and that the precipitation response is effectively truncated by limiting $\frac{n}{b}$ to 24 months. Response to abstraction is very rapid for all locations (Figure 5(d) to (f)). The response to canal irrigation is highly variable between the three locations (Figure 5(g) to (i)); there is a fairly rapid response in areas of groundwater level decline, but a protracted response for the location with groundwater level rise.

Figure 5.

Best-fit impulse response functions, shown as black lines, for precipitation (top row), abstraction (middle row), and canal irrigation (bottom row) drivers for the three locations shown in Figure 3. The area under each response function curve gives the zeroth moment for precipitation ( $M_{0, p}$ ), abstraction ( $M_{0, w}$ ), and canal irrigation ( $M_{0, c}$ ) as a quantitative scalar measure of the response. Grey bars show the distributions of best-fit zeroth moment values from all model outcomes at these locations with $N S E >$ 0.2. Left column (a, d, and g) shows results from the location with groundwater level rise shown in Figure 4(a). Middle column (b, e, and h) shows results from the location with groundwater level decline in the area of the Ghaggar–Hakra paleochannel shown in Figure 4(b). Right column (c, f, and i) shows results from the location with groundwater level decline on the Sutlej fan shown in Figure 4(c). Note changes in vertical and horizontal scales.

There is substantial variability in M₀ values for all three drivers, even when we consider only model outcomes with $N S E >$ 0.2, as visible as the histograms in Figure 5. The locations with groundwater level decline generally behave similarly, with well-defined modal values for $θ_{w}$ and $θ_{h}$ in particular, and are distinct from the location with groundwater level rise. The response function $θ_{p}$ for precipitation is notably variable for all sites, showing a wide range of permissible $M_{0, p}$ values that indicate a range of both amplitudes and time delays. This suggests that the model is not particularly sensitive to the details of the precipitation response, as a wide range of response functions can yield acceptable model behavior. While the form of the response function $θ_{w}$ for abstraction is similar across all three locations (Figure 5(d) to (f)), it varies substantially in amplitude. The amplitude is greater for the paleochannel location (Figure 5(e)) and the histogram is skewed toward high negative values of $M_{0, w}$ , indicating a strong negative response to abstraction at this location. The Sutlej fan location also shows a negative response with high negative values of $M_{0, w}$ , but the best-fit amplitude is somewhat smaller (Figure 5(f)). Finally, $M_{0, c}$ shows a clear difference between areas with groundwater level rise and fall: values are predominantly $<$ 50 in the areas with groundwater level decline but are much larger in the location with groundwater rise. This demonstrates that the canal irrigation driver has a dominant influence on groundwater level rise but is much less important in areas of decline. Overall, the best-fit values and distributions of M₀ fit our general expectations for model behavior: large values of $M_{0, w}$ for the abstraction driver in areas where groundwater level has significantly declined (Figure 5(e) and (f)), relatively similar values of $M_{0, p}$ for the precipitation driver irrespective of location, and larger values of $M_{0, c}$ for the canal irrigation driver in areas that have experienced groundwater level rise.

3 Spatial variation of the zeroth moment

There are substantial and systematic variations in M₀ for the precipitation, abstraction, and canal irrigation drivers across the study area (Figure 6). Recall that a high value of M₀ means that the driver has a large influence on the groundwater level, and a low value means there is little influence. The response to precipitation is fairly uniform across the study area, except for a hotspot on the central part of the Sutlej fan (Figure 6(a) and (b)). This area has the highest $M_{0, p}$ values, indicating high sensitivity to precipitation, along with a low coefficient of variation (Figure 6(c)). It is not clear why this area should be distinct from adjacent parts of the fan, although it is worth noting that groundwater level decline in this area is both high and not well predicted according to $M A E$ (Figure 3(d)). It is thus possible that some of the mismatch between observed and modeled groundwater levels will be due to a precipitation response that is either too large or too delayed.

Figure 6.

Spatial variation in the zeroth moment M₀ of the impulse response functions for precipitation (top row), abstraction (middle row), and canal irrigation (bottom row) drivers. In each row, the left panel shows M₀ for the best-fit model result in each grid cell, the middle panel shows the median M₀ for all model runs with $N S E >$ 0.2, and the right panel shows the coefficient of variation ( $C V$ ) for those runs. The zeroth moment for the precipitation driver $M_{0, p}$ is spatially fairly uniform apart from high values on the central Sutlej fan (a and b), corresponding to low to moderate $C V$ values. In contrast, the zeroth moment for the abstraction driver $M_{0, w}$ is strongly negative across the central Sutlej and Yamuna fans, and especially in the interfan area and along the Ghaggar–Hakra paleochannel (d and e). The $C V$ in those areas is less than one which indicates a significant adverse effect of abstraction on groundwater level. The zeroth moment for the canal irrigation driver $M_{0, c}$ is small in the northwestern part of the study area but large in the southwest where groundwater levels have risen (g and h). The $C V$ is less than one for locations with groundwater level rise, but much higher where $M_{0, c}$ is small and small differences between model runs lead to large $C V$ values.

The zeroth moment for the abstraction stress shows a strong negative response (Figure 6(d) and (e)) in the central part of the study area, and most notably in the interfan area between the Sutlej and Yamuna fans and along the Ghaggar–Hakra paleochannel (Figure 1). This high sensitivity to abstraction is visible in both best-fit and median model results, with a low coefficient of variation (Figure 6(f)), and is centered on the areas of greatest groundwater decline (Figure 3(a)). The distal parts of both fans show much lower sensitivity to abstraction, corresponding to lower overall groundwater depletion. Distal areas have some high values of the coefficient of variation; these areas have very small zeroth moments so the coefficient of variation is sensitive to small variations between model runs.

High values of the zeroth moment for canal irrigation $M_{0, c}$ are limited to distal parts of the study area (Figure 6(g) and (h)) and are strongly associated with areas of groundwater level rise. These areas also have low values of the coefficient of variation, indicating consistency between acceptable model outputs (Figure 6(i)). As with the abstraction driver, very low $M_{0, c}$ values in the northeastern part of the study area are associated with high coefficients of variation, indicating that groundwater levels in this region are not sensitive to the canal irrigation driver.

V Discussion

The TFN model yields insights into the response of groundwater levels to the four most common drivers that determine the groundwater depletion rate in northwestern India. Here we first discuss each individual parameter that sets the impulse response function and the zeroth moment of that response. Second, we link the spatial variation of the responses (as represented by the zeroth moment and the parameters) with the underlying geological heterogeneity in the aquifer system. Third, we discuss the implications of our model outcomes for understanding groundwater level changes, along with model limitations. Finally, we provide our recommendations for future sustainable groundwater management.

1 Link to specific drivers

Spatial variations in M₀ are explained by the different parameters that determine θ, and so it is useful to examine the variations in those parameters for the median model solutions and their link to specific drivers. High values for $M_{0, p}$ in the center of the Sutlej fan, and to some extent in the Yamuna fan (Figure 6(b)), are mainly due to high values of A in equation (3) (Figure 7(a)). This parameter adjusts the area of response to precipitation, and may take high values in the central parts of the fans because of the abundant thick aquifer bodies in these areas compared to the distal or interfan areas (Van Dijk et al., 2016a). The ratio $\frac{n}{b}$ determines the number of months it will take to reduce the groundwater perturbation by half after a precipitation event; while we have limited it to a maximum of 24 months. Figure 7(b) shows no clear spatial pattern. Similarly, there is little evidence of spatial variation in the evapotranspiration factor f (Figure 7(c)), which sets the relation between p and e in equation (5).

Figure 7.

Model parameters associated with the median solution for all grid cells. (a) Parameter A is larger for the areas with groundwater level rise and especially for the central part of the Sutlej fan. (b) The ratio ( $\frac{n}{b}$ ) sets the response time to a precipitation event, and shows no clear spatial pattern. (c) The evapotranspiration fraction f shows no clear spatial pattern. (d) There is no spatial pattern observed for α, which indicates the distance from a well over which abstraction should be affected (note that each grid cells contains multiple wells). (e) The parameter γ is highest where $M_{0, w}$ is high, likely reflecting low hydraulic conductivity values in the interfan area and Ghaggar–Hakra paleochannel where groundwater levels have declined. (f) The parameter λ varies in conjunction with $M_{0, c}$ , and is highest in areas that have experienced groundwater level rise.

More interesting are the parameters that determine the variation of $M_{0, w}$ (equation (6)) because of the relationship between groundwater level change and w (Asoka et al., 2017). The parameter α (Figure 7(d)), which is a multiple of the reciprocal of aquifer diffusivity ( $S / T$ ), shows no clear spatial pattern. However, the variable γ (Figure 7(e)) shows the same spatial pattern as $M_{o, w}$ . High γ values correspond to low aquifer transmissivity T, which in turn is related to the hydraulic conductivity (where T is the integral of the hydraulic conductivity over the saturated aquifer thickness). This result is consistent with substantially lower hydraulic conductivity values in the interfan area and the Ghaggar–Hakra paleochannel compared to other parts of the study area. Finally, the response to canal irrigation depends solely on the parameter λ, which varies spatially in much the same way as $M_{0, c}$ (Figure 7(f)). The values are high only where $M_{0, c}$ is also high, corresponding to areas of observed groundwater level rise; in contrast, low values elsewhere yield a low stationary response ( $M_{0, c}$ ) irrespective of the value of the canal irrigation driver.

2 Spatial relations between abstraction and geomorphology

The zeroth moment of the response to the abstraction driver $M_{0, w}$ varies with distance from the Himalayas as well as along the strike of the foreland from northwest to southeast (Figure 6). This pattern is distinct from the annual abstraction pattern, which increases towards the Sutlej River at the mountain front. Intriguingly, the spatial variation in $M_{0, w}$ bears a striking resemblance to the regional geomorphic and geological heterogeneity of the study region, as documented by Van Dijk et al. (2016a). The largest negative values of $M_{0, w}$ correspond with the interfan area between the Sutlej and Yamuna fans, which is characterized by (1) lower bulk aquifer content and (2) thinner individual aquifer bodies compared to the sedimentary fans on either side (Van Dijk et al., 2016a). Our earlier work (Singh et al., 2017; Van Dijk et al., 2016a) also showed that aquifer bodies are not continuous in the subsurface across the interfan area, and that there are important lateral disconnections between the Sutlej and Yamuna fan (each of which is deposited by a distinct hinterland sediment source) and the interfan area (which is sourced only from the Himalayan foothills). These lateral disconnections are likely to place strong limitations on lateral groundwater flow and recharge, which in the TFN model is represented by high γ values (Figure 7(e)).

High levels of groundwater decline and high negative $M_{0, w}$ values are also observed along the Ghaggar–Hakra paleochannel, which runs down the boundary between the Sutlej and Yamuna fans. Singh et al. (2017) demonstrated that the paleochannel is underlain by an incised valley fill, consisting of a 30 m thick succession of coarse- to fine-grained sands capped by thin silts and clays. Weissmann et al. (2004) showed that groundwater moves faster through an incised valley fill than through an open alluvial fan system, which in turn affects the recharge rates. The direction of groundwater flow is mainly from the incised valley fill to the open-fan deposits (Weissmann et al., 2004), which makes the valley fill a good location for enhanced or artificial recharge. This also makes the incised valley fill more sensitive to abstraction, however, resulting in rapid groundwater level decline as observed here. This issue is likely exacerbated by the grain-size difference between the coarse valley fill along the Ghaggar–Hakra paleochannel and the finer-grained sediment on its flanks (Singh et al., 2017; Van Dijk et al., 2016a), leading to poor lateral hydraulic conductivity and limiting the amount of lateral recharge into the paleochannel from the surrounding fans. Similarly high negative $M_{0, w}$ values are seen locally along the incised present-day channels of the Sutlej and Yamuna rivers (Figure 6(d) and (e)) and may indicate similar relations, although we do not have direct evidence of the sedimentary architecture in those locations. More distal parts of the paleochannel, characterized by low abstraction rates but high values of groundwater level decline, are not well simulated by the TFN model. Elucidation of the drivers behind groundwater level decline in such a complex three-dimensional stratigraphic setting may require a more sophisticated model than our spatially-averaged one-dimensional approach.

3 Model implications and limitations

The time series approach was applied to study the spatial variation of groundwater level in response to multiple drivers for a regional hotspot of groundwater depletion in northwestern India. The model incorporated impulse response functions θ to four imposed drivers: precipitation, PET, abstraction, and canal irrigation. These response functions were calibrated against the observed groundwater levels to produce a cell-by-cell prediction of modeled groundwater levels through time. Supplementary figure 4 showed that there was only a weak relation between total abstraction or precipitation and the observed groundwater level changes on a cell-by-cell basis, irrespective of whether the level had risen or declined. The time series analysis demonstrates that there are strong spatial differences in the response to the four modeled drivers, as quantified by the zeroth moment M₀ which is a stationary measure of the response to a driver. This, in turn, suggests that scaling the total abstraction or precipitation by M₀ may improve the correlation between the total volumes and observed groundwater level change. However, when applying this to precipitation, a negative relation with groundwater level change is derived (Figure 8(a)).

Figure 8.

Observed groundwater level change on a cell-by-cell basis plotted against (a) the median zeroth moment for precipitation $M_{0, p}$ times total precipitation (p), and (b) the median zeroth moment for abstraction $M_{0, w}$ times total abstraction (w). The lines indicate best-fit relations derived from linear least-squares regression. Open symbols show all cells, while closed symbols show only those cells that experienced groundwater level decline over the period 1974–2010. Total precipitation scaled by $M_{0, p}$ indicates an unexpected relation to groundwater level decline. In contrast, there is a clear relation between abstraction and the groundwater level change, especially for those areas that experienced groundwater level decline. (c) Monthly canal irrigation, $I / S_{y}$ , for groundwater level rise areas where λ is significant and canal irrigation is needed to fit the observations (equation (7)), shows that irrigation goes up to 6 mm/day.

In contrast, total abstraction ( $w_{t o t}$ ) scaled by the stationary response $M_{0, w}$ shows a clear positive correlation with groundwater level change, such that a larger negative value for $M_{0, w} \cdot w_{t o t}$ corresponds to more groundwater decline (Figure 8(b)). This is not surprising, as we would expect a fall in groundwater level if either total abstraction increases or there is a stronger stationary response to abstraction. Overall, these results confirm that changes in groundwater levels are predominantly due to increases in abstraction compared to relatively small declines in annual rainfall, as was argued independently by Asoka et al. (2017). We further infer that observed groundwater level decline relates to abstraction volume combined with the stationary response, which in turn is controlled by the hydraulic conductivity of the aquifer (γ).

Initially, we focused on the abstraction driver to explain the groundwater depletion observed in GRACE (Chen et al., 2014; Kumar et al., 2006; MacDonald et al., 2016; Richey et al., 2015; Rodell et al., 2009; Shah, 2009), but return flow from canal irrigation and canal leakage make up an important component for artificial recharge of the aquifer in the southwestern part of the study area. The canal network in northwestern India was constructed during the nineteenth and twentieth centuries, and leakage from canals has historically been a significant source of recharge (Cheema et al., 2014; MacDonald et al., 2016; Raza et al., 2013). The southwestern part of the region has been more dominantly irrigated by canals and tube wells, whereas the districts closer to the Himalayas were mainly irrigated by tube wells (Jeevandas et al., 2008). Groundwater level rise can be explained well by the stationary response to canal irrigation. This behavior is seen in the parts of the study area that were identified as predominantly recharged by canal irrigation by Joshi et al. (2018) on the basis of groundwater chemistry observations. The calibrated TFN model provides an estimate of modeled monthly irrigation driver scaled by the specific yield, $R_{c} / S_{y}$ . Model values indicate that for the locations with groundwater level rise, recharge by canal irrigation is up to 6 mm/day (Figure 8(c)). The value I is scaled by $S_{y}$ , which means that assuming $S_{y}$ is between 0.03 and 0.32 for silt and medium sand, respectively (Johnson, 1967), I ranges between 0.18 to 1.92 mm/day. This range is in agreement with the findings of Cheema et al. (2014).

In general, the results show that abstraction and canal irrigation drivers can explain the first-order pattern of groundwater level change for a large part of the study area for the period 1974–2010 (Figure 3), but there are some limitations in the TFN model. As result of the limitations areas of weaker model performance can fall into two categories: locations where there were no acceptable model runs with $N S E >$ 0.2, and locations with larger disparities between observed and modeled groundwater levels as indicated by the $M A E$ (Figure 3). The latter is particularly noticeable in the center of the Sutlej fan, where modeled groundwater levels are higher than observations. Three reasons could explain this disparity.

First, areas with more stable groundwater levels are not explained well by the combination of the four drivers. This is specifically true for the transitional area between the regions of most pronounced groundwater level rise and fall. It may also be possible that important drivers have not been adequately included within the model (Von Asmuth et al., 2008).

Second, data availability, resolution, and accuracy are all highly variable between the different drivers, and affect the outcome of the TFN model. The TFN model includes monthly values for the four drivers, but groundwater level is only available twice per year. Although this does not pose a problem for the model, the impulse response functions are continuous in time and the predictions are therefore not fixed to a certain time interval (Von Asmuth et al., 2008).

Third, observed groundwater levels in many areas continued to decline after 2006 despite the apparent stabilization of abstraction rates (Figure 4(b) and (c)). The observed groundwater level time-series shows an increasing rate of decline which may be due to limited recharge, either vertically because of water loss before reaching the water table (Hoque et al., 2007), or horizontally as surrounding aquifer bodies are depleted as well. This non-linear behavior of the groundwater level is difficult to predict with our one-dimensional implementation of the statistical TFN model approach, in which each grid cell is modeled independently, and has important implications for the sustainability of the aquifer system. Therefore, future studies that investigate the sustainability of the groundwater resource should take into account lateral groundwater flow over distances of greater than 10 km.

4 Management recommendations and wider model application

While the time series approach outlined here is highly simplified, the congruence between our key results and independent estimates of regional-scale variations in aquifer properties (Singh et al., 2017; Van Dijk et al., 2016a) and recharge mechanisms (Joshi et al. 2018) suggests that it has some predictive skill to map out the areas of maximum abstraction and their geomorphic/stratigraphic controls. This will then determine the most appropriate management strategies (Sinha and Densmore, 2016), such as where to plan artificial recharge and where to advise crop management. So what lessons can be inferred for better management strategy of this regional-scale aquifer system?

Management interventions planned by the Government of India are limited to a decrease in groundwater abstraction (via piped water supply and crop management) and an increase in recharge (via artificial recharge pits and rainwater harvesting) (CGWB, 2013). Artificial recharge schemes are likely to be most effective where the response to abstraction is rapid, but not necessarily where the stationary response $M_{0, w}$ is large – because those places may have a long, drawn-out response, which is determined by α. Rapid response to abstraction is seen everywhere (Figure 4) but the zeroth moment is particularly high in the interfan area and along the Ghaggar–Hakra paleochannel. These areas are thus likely to be less effective locations for artificial recharge schemes, due to the thinner and less-abundant aquifer bodies in the subsurface (Figure 1). Conversely, the middle portions of the Sutlej and Yamuna fans appear better suited to artificial recharge as they combine a rapid response with more moderate values of the zeroth moment.

Groundwater level rise in the southwestern part of the study area is largely insensitive to temporal variations in precipitation or abstraction, and appears to be driven primarily by canal irrigation. While canal irrigation is estimated by our model rather than used as an input, the results (expressed in terms of irrigation stress normalized by specific yield) are spatially variable compared to the uniform estimates from Cheema et al. (2014). It thus appears likely that canal irrigation in this part of the basin has led to substantial aquifer recharge since 1974, which can provide insights for future management of the depleted aquifers in the study area. The sensitivity of distal areas to the canal irrigation driver also suggests that improved management of return flow, and efforts to decrease canal seepage, would help to limit water-level rises and consequent waterlogging.

It is widely recognized that there is a groundwater crisis in many of the Earth’s major aquifer systems (Richey et al. 2015; Richts et al. 2011) where unmanaged pumping of critically important groundwater resources has led to rapid rates of groundwater depletion (Famiglietti, 2014). Whilst groundwater use for irrigation has significantly increased crop yields and food security in many parts of the world (Khan and Hanjra, 2009; Siebert et al., 2010), this depletion threatens the sustainability of food production, and water and food security (Dalin et al., 2017). It is important to project changes in large-scale groundwater resources into the future to explore how best to adapt environmental policy and management (Green et al., 2011). However, most large-scale or global hydrological models, which could be used for this purpose, have a limited representation of groundwater. This has been due to the lack of groundwater level and subsurface property data, and inadequate representations of the interactions and feedbacks between human water use and management, and natural systems (Nazemi and Wheater, 2015). Addressing these issues is a current area of active research and reasonable representations of groundwater in large-scale models are beginning to appear (De Graaf et al., 2017; Pokhrel et al., 2015; Zeng et al., 2018). Parameterizing their subsurface properties, however, remains a challenge. As a first order approximation, aquifer hydraulic properties can be derived from national and global geological maps (Bhanja et al., 2016; Gleeson et al., 2014), but we consider that these will need to be refined if spatial variability in changes in groundwater levels are to be simulated adequately. Consequently, model results deviate strongly from groundwater level observations (Scanlon et al., 2018). To support the exploration of variability in aquifer properties, and investigate the response of groundwater levels to multiple-drivers, alternative parsimonious modelling approaches, such as that outlined here, will be valuable. Our study supports the conclusion of Shapoori et al. (2015c) that TFN models can provide important insights into hydrogeology, hydrological processes, and the response of the system to multiple drivers. Because they require minimal prior assumptions they are easily and rapidly transferable to other groundwater systems.

VI Conclusions

We have used a TFN time series approach to show that groundwater decline in the key regional-scale alluvial aquifer system of northwestern India over the period 1974–2010 is strongly influenced by abstraction, but that the spatial pattern of groundwater level decline is not simply based on abstraction rate alone. Time series analysis of 664 grid cells shows water levels can be predicted to first order by considering the time-varying response to four input drivers: precipitation, PET, abstraction, and canal irrigation. The results show that groundwater level decline across the northeastern part of the study area relates to the total volume of abstraction scaled by the modeled zeroth moment of the aquifer system to abstraction ( $M_{o, w}$ ). The zeroth moment is the inverse of the effective porosity of the aquifer for each individual driver at that location, and varies systematically across the study area. Much of that spatial variation in the response to abstraction can be explained by the underlying geological heterogeneity of the alluvial aquifer system, in which the storage capacity is controlled by the aquifer percentage. Large declines in groundwater level are observed in the interfan area between the Sutlej and Yamuna rivers, and along the Ghaggar–Hakra paleochannel, where aquifer percentage is low and disconnected from the large fan system. These areas show exceptionally large values of $M_{0, w}$ and low modeled hydraulic conductivity, corresponding to independent estimates of low aquifer abundance and thin aquifer bodies (Van Dijk et al., 2016a).

Our time series analysis provides a preliminary first-order approach for understanding the spatial controls of groundwater level changes in this critical region. The approach is effective in determining the relative importance of different stresses in driving the evolution of groundwater levels, but cannot reproduce fine-scale impacts from individual events or incorporate the complex three-dimensional stratigraphy of the alluvial aquifer system. We close by suggesting that an interdisciplinary approach that combines hydrology and geological heterogeneity should be considered in any future approaches to regional-scale aquifer management.

Supplemental material

Supplemental Material, VanDijketal_SI_SpatialVariationInGroundwaterStresses - Spatial variation of groundwater response to multiple drivers in a depleting alluvial aquifer system, northwestern India

Supplemental Material, VanDijketal_SI_SpatialVariationInGroundwaterStresses for Spatial variation of groundwater response to multiple drivers in a depleting alluvial aquifer system, northwestern India by Wout M van Dijk, Alexander L Densmore, Christopher R Jackson, Jonathan D Mackay, Suneel K Joshi, Rajiv Sinha, Shashank Shekhar and Sanjeev Gupta in Progress in Physical Geography: Earth and Environment

Footnotes

Acknowledgements

The project was initiated by W.M. van Dijk and A.L. Densmore. W.M van Dijk implemented all the data and performed the model runs and analysis together with C.R. Jackson. A.L. Densmore, J.D. Mackay, and C.R. Jackson contributed by discussing the results from the model and its implications. J.D. Mackay modeled the abstraction data. S.K. Joshi, S. Sheshank, R. Sinha and S. Gupta collected and processed the groundwater level data. The authors would like to thank Antoine Lafare and Edwin Sutanudjaja for earlier discussions. Constructive and positive reviews by two anonymous reviewers helped to clarify and strengthen the manuscript. The time-series analysis model AquiMod is free available online through the British Geological Survey (https://www.bgs.ac.uk/research/environmentalModelling/aquimod.html). The meteorological data of the Indian Meteorological Department is publicly available from a web portal (). To obtain the input data and model used in this paper, please contact the corresponding author.

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This project is supported by the UK Natural Environment Research Council and the Indian Ministry of Earth Sciences under the Changing Water Cycle - South Asia program (Grant NE/I022434/1). The abstraction modeling was supported by the UK Natural Environment Research Council (Grant NE/I022590/1).

ORCID iD

Wout M van Dijk

Supplemental material

Supplemental material for this article is available online.

References

Aeschbach-Hertig

Gleeson

(2012) Regional strategies for the accelerating global problem of groundwater depletion. Nature Geoscience 5: 853–861.

Allen

JRL

(1978) Studies in fluviatile sedimentation: an exploratory quantitative model for the architecture of avulsion-controlled alluvial studies. Sedimentary Geology 21(2): 129–147.

Allen

JRL

(1984) Sedimentary Structures, Their Character and Physical Basis. Amsterdam: Elsevier Scientific Publishing.

Anderson

(1989) Hydrogeologic facies models to delineate large-scale spatial trends in glacial and glaciofluvial sediments. Geological Society of America Bulletin 101(4): 501–511.

Asoka

Gleeson

Wada

, et al. (2017) Relative contribution of monsoon precipitation and pumping to changes in groundwater storage in India. Nature Geoscience 10: 109–119.

Barrett

Charbeneau

(1997) A parsimonious model for simulating flow in a karst aquifer. Journal of Hydrology 196(1–4): 47–65.

Berendrecht

Heemink

Van Geer

, et al. (2003) Decoupling of modeling and measuring interval in groundwater time series analysis based on response characteristics. Journal of Hydrology 278(1): 1–16.

Bhanja

Mukherjee

Saha

, et al. (2016) Validation of GRACE based groundwater storage anomaly using in-situ groundwater level measurements in India. Journal of Hydrology 543(Part B): 729–738.

Bierkens

MFP

Knotters

Hoogland

(2001) Space-time modeling of water table depth using a regionalized time series model and the Kalman filter. Water Resources Research 37(5): 1277–1290.

10.

Bridge

(1993) Description and interpretation of fluvial deposits: A critical perspective. Sedimentology 40(4): 801–810.

11.

Bridge

Leeder

(1979) A simulation model of alluvial stratigraphy. Sedimentology 26(5): 617–644.

12.

CGWB (2011) Dynamic ground water resources of India (as on 31 March 2009). Technical report, Central Groundwater Board, Ministry of Water Resources, Government of India, Faridabad.

13.

CGWB (2013) Master plan for artificial recharge to ground water in India. Technical report, Central Groundwater Board, Ministry of Water Resources, Government of India, Faridabad.

14.

Cheema

MJM

Immerzeel

Bastiaanssen

WGM

(2014) Spatial quantification of groundwater abstraction in the irrigated Indus basin. Ground Water 52(1): 25–36.

15.

Chen

Zhang

, et al. (2014) Long-term groundwater variations in Northwest India from satellite gravity measurements. Global and Planetary Change 116: 130–138.

16.

Clift

Carter

Giosan

, et al. (2012) U-Pb zircon dating evidence for a Pleistocene Sarasvati River and capture of the Yamuna River. Geology 40(3): 211–214.

17.

Criss

Winston

(2008) Do Nash values have value? Discussion and alternate proposals. Hydrological Processes 22(14): 2723–2725.

18.

Dalin

Wada

Kastner

, (2017) Groundwater depletion embedded in international food trade. Nature 543: 700–704.

19.

De Graaf

IEM

Van Beek

LPH

Gleeson

, et al. (2017) A global-scale two-layer transient groundwater model: Development and application to groundwater depletion. Advances in Water Resources 102: 53–67.

20.

Donselaar

Overeem

(2008) Connectivity of fluvial point-bar deposits: An example from the Miocene Huesca fluvial fan, Ebro Basin, Spain. AAPG Bulletin 92(9): 1109–1129.

21.

Famiglietti

(2014) The global groundwater crisis. Nature Climate Change 4: 945–948.

22.

Flood

Hampson

(2015) Quantitative analysis of the dimensions and distribution of channelized fluvial sandbodies within a large outcrop datdata: Upper Cretaceous Blackhawk formation, Wasatch Plateau, Central Utah, USA. Journal of Sedimentary Research 85(4): 315–336.

23.

Fogg

(1986) Groundwater flow and sand body interconnectedness in a thick multiple-aquifer system. Water Resources Research 22(5): 679–694.

24.

Geddes

(1960) Alluvial morphology of the Indo-Gangetic Plain: Its mapping and geographical significance. Transactions of Institute of British Geographers 28: 253–276.

25.

Gleeson

Moosdorf

Hartmann

, et al. (2014) A glimpse beneath earth’s surface: GLobal HYdrogeological MaPS (GLHYMPS) of permeability and porosity. Geophysical Research Letters 41(11): 3891–3898.

26.

Gleeson

VanderSteen

Sophocleous

, et al. (2010) Groundwater sustainability strategies. Nature Geoscience 3: 378–379.

27.

Green

Taniguchi

Kooi

, et al. (2011) Beneath the surface of global change: Impacts of climate change on groundwater. Journal of Hydrology 405(3–4): 532–560.

28.

Guthke

(2017) Degensible model complexity: A call for data-based and goal-oriented model choice. Groundwater 55(5): 646–650.

29.

Hantush

(1956) Analysis of data from pumping test in leaky aquifers. Transactions, American Geophysical Union 37(6): 702–714.

30.

Heller

Paola

(1996) Downstream changes in alluvial architecture: An exploration of controls on channel-stacking patterns. Journal of Sedimentary Research 66(2): 297–306.

31.

Hill

Tideman

(2007) Effective Groundwater Model Calibration: With Analysis of Data, Sensitivities, Predictions, and Uncertainty. Hoboken, NJ: John Wiley & Sons.

32.

Holbrook

(2001) Origin, genetic interrelationships, and stratigraphy over the continuum of fluvial channel-form bounding surfaces: An illustration from middle Cretaceous strata, southeastern Colorado. Sedimentary Geology 144(3–4): 179–222.

33.

Hong

Hama

Suenaga

, et al. (2017) Simplified lumped groundwater model to simulate Nitrate concentration dynamics. Journal of Hydrologic Engineering 22(10): 04017043.

34.

Hoque

Ahmed

(2007) Declining groundwater level and aquifer dewatering in Dhaka metropolitan area, Bangladesh: Causes and quantification. Hyrdogeology Journal 15(8): 1523–1534.

35.

Jeevandas

Singh

Kumar

(2008) Concerns of groundwater depletion and irrigation efficiency in Punjab agriculture: A micro-level study. Agricultural Economics Research Review 21(2): 191–199.

36.

Johnson

(1967) Specific yield – compilation of specific yields for various materials. Technical report, U.S. Geological Survey.

37.

Joshi

Rai

Sinha

, et al. (2018) Tracing groundwater recharge sources in the northwestern Indian alluvial aquifer using water isotopes δ¹⁸O, δ²H and ³H. Journal of Hydrology 559: 835–847.

38.

Kazumba

Oron

Honjo

, et al. (2008) Lumped model for regional groundwater flow analysis. Journal of Hydrology 359(1–2): 131–140.

39.

Khan

Sinha

(2019) Discovering ‘buried’ channels of the Palaeo-Yamuna river in NW India using geophysical evidence: Implications for major drainage reorganization and linkage to the Harappan Civilization. Journal of Applied Geophysics 167: 128–139.

40.

Khan

Hanjra

(2009) Footprints of water and energy input in food production: Global perspectives. Food Policy 34(2): 130–140.

41.

Kumar

Sahai

Kumar

, et al. (2006) High-resolution climate change scenario for India for the 21st century. Current Science 90(3): 334–345.

42.

Larue

Hovadik

(2006) Connectivity of channelized reservoirs: a modelling approach. Petroleum Geoscience 12(4): 291–308.

43.

Leeder

(1978) A quantitative stratigraphic model for alluvium, with spatial reference to channel deposit density and interconnectedness. In: Miall

(ed) Fluvial Sedimentology, Volume 5. Canadian Society of Petroleum Geologists Memoir, 587–596.

44.

Legates

McCabe

Jr (1999) Evaluating the use of “goodness-of-fit” measures in hydrologic and hydroclimatic model validation. Water Resources Research 35(1): 233–241.

45.

Long

(2015) RRAWFLOW: Rainfall-response aquifer and watershed flow model (v1.15). Geoscientific Model Development 8: 865–880.

46.

MacDonald

Bonsor

Ahmed

, et al. (2016) Groundwater quality and depletion in the Indo-Gangetic Basin mapped from in situ observations. Nature Geoscience 9: 762–766.

47.

Mackay

Jackson

Wang

(2014) A lumped conceptual model to simulate groundwater level time-series. Environmental Modelling and Software 61: 229–245.

48.

McCuen

Knight

Cutter

(2006) Evaluation of the Nash–Sutcliffe efficiency index. Journal of Hydrologic Engineering 11(6): 597–602.

49.

Miall

(1985) Architectural-element analysis: A new method of facies analysis applied to fluvial deposits. Earth-Science Reviews 22(4): 261–308.

50.

Nash

Sutcliffe

(1970) River flow forecasting through conceptual models. Part 1. A discussion of principles. Journal of Hydrology 10(3): 282–290.

51.

Nazemi

Wheater

(2015) On inclusion of water resource management in Earth system models. Part 2. Representation of water supply and allocation and opportunities for improved modeling. Hydrology and Earth System Science 19(1): 63–90.

52.

Obergfell

Bakker

Zaadnoordijk

, et al. (2013) Deriving hydrogeological parameters through time series analysis of groundwater head fluctuations around well fields. Hydrogeology Journal 21(5): 987–999.

53.

Park

Parker

(2008) A simple model for water table fluctuations in response to precipitation. Journal of Hydrology 356(3–4): 344–349.

54.

Pokhrel

Koirala

Yeh

PJF

, et al. (2015) Incorporation of groundwater pumping in a global land surface model with the representation of human impacts. Water Resources Research 51(1): 78–96.

55.

Pozdniakov

Shestakov

(1998) Analysis of groundwater discharge with a lumped-parameter model, using a case study from Tajikistan. Hyrdogeology Journal 6(2): 226–232.

56.

Raza

Latif

Shakir

(2013) Long-term effectiveness of lining tertiary canals in the Indus Basin of Pakistan. Irrigation and Drainage 62(1): 16–24.

57.

Refsgaard

Christensen

Sonnenborg

, et al. (2012) Review of strategies for handling geologcial uncertainty in groundwater flow and transport modeling. Advances in Water Resources 36: 36–50.

58.

Renard

Allard

(2013) Connectivity metric for subsurface flow and transport. Advances in Water Resources 51: 168–196.

59.

Richey

Thomas

, et al. (2015) Uncertainty in global groundwater storage estimates in a total groundwater stress framework. Water Resources Research 51(7): 5198–5216.

60.

Richts

Struckmeier

Zaepke

(2011) WHYMAP and the groundwater resources of the world 1:25,000,000. In: Jones

JAJ

(ed) Sustaining Groundwater Resources: A Critical Element in the Global Water Crisis. Berlin: Springer.

61.

Rodell

Velicogna

Famiglietti

(2009) Satellite-based estimates of groundwater depletion in India. Nature 460: 999–1002.

62.

Rojas

Kahunde

Peeters

, et al. (2010) Application of a multimodel approach to account for conceptual model and scenario uncertainties in groundwater modelling. Journal of Hydrology 394(3–4): 416–435.

63.

Scanlon

Zhang

Save

, et al. (2018) Global models underestimate large decadal declining and rising water storage trends relative to GRACE satellite data. PNAS 115(6): E1080–E1089.

64.

Schaefli

Gupta

(2007) Do Nash values have value? Hydrological Processes 21(15): 2075–2080.

65.

Shah

(2009) Climate change and groundwater: India’s opportunities for mitigation and adaptation. Environmental Research Letters 4(3): 035005.

66.

Shapoori

Peterson

Western

, et al. (2015a) Decomposing groundwater head variations into metereological and pumping components: A synthetic study. Hydrogeology Journal 23(7): 1431–1448.

67.

Shapoori

Peterson

Western

, et al. (2015b) Estimating aquifer properties using groundwater hydrograph modelling. Hydrological Processes 29(26): 5424–5437.

68.

Shapoori

Peterson

Western

, et al. (2015c) Top-down groundwater hydrograph time-series modeling for climate-pumping decomposition. Hydrogeology Journal 23(4): 819–836.

69.

Sheets

Hickson

Paola

(2002) Assembling the stratigraphic record: Depositional patterns and time-scales in an experimental alluvial basin. Basin Research 14(3): 287–301.

70.

Siebert

Burke

Faures

, et al. (2010) Groundwater use for irrigation – a global inventory. Hydrology and Earth System Science 14: 1863–1880.

71.

Singh

Paul

Sinha

, et al. (2016) Geochemistry of buried river sediments from Ghaggar Plains, NW India: Multi-proxy records of variations in provenance, paleoclimate, and paleovegetation patterns in the late quaternary. Palaeogeography, Palaeoclimatology, Palaeoecology 449: 85–100.

72.

Singh

Thomson

Sinha

, et al. (2017) Counter-intuitive influence of Himalayan river morphodynamics on Indus Civilisation urban settlements. Nature Communications 8: 1617.

73.

Singh

(2000) Environmental consequences of agricultural development: A case study from the Green Revolution state of Haryana, India. Agriculture, Ecosystems and Environment 82(1–3): 97–103.

74.

Sinha

Densmore

(2016) Focus on sustainable groundwater management. Economic and Political Weekly 51(52): 53–55.

75.

Sinha

Yadav

Gupta

, et al. (2013) Geo-electric resistivity evidence for subsurface palaeochannel systems adjacent to Harappan sites in northwest India. Quaternary International 308–309: 66–75.

76.

Straub

Paola

Mohrig

, et al. (2009) Compensational stacking of channelized sedimentary deposits. Journal of Sedimentary Research 79(9): 673–688.

77.

Tiwari

Wahr

Swenson

(2009) Dwindling groundwater resources in northern India, from satellite gravity observations. Geophysical Research Letters 36(18): L18401.

78.

UNDP (1985) Ground-water studies in the Ghaggar River Basin in Punjab, Haryana and Rajasthan. Final technical report, volume 1 (of 3 volumes), Surface and Ground Water Resources. New York: United Nations Development Programme.

79.

Van de Lageweg

Schuurman

Cohen

, et al. (2016a) Preservation of meandering river channels in aggrading channel belts. Sedimentology 63(3): 586–608.

80.

Van de Lageweg

Van Dijk

Box

, et al. (2016b) Archimetrics: A quantitative tool to predict three-dimensional meander belt sandbody heterogeneity. The Depositional Record 2(1): 22–46.

81.

Van Dijk

Densmore

Singh

, et al. (2016a) Linking the morphology of fluvial fan systems to aquifer stratigraphy in the Sutlej-Yamuna plain of northwest India. Journal of Geophysical Research – Earth Surface 121(2): 201–222.

82.

Van Dijk

Densmore

Sinha

, et al. (2016b) Reduced-complexity probabilistic reconstruction of alluvial aquifer stratigraphy, and application to sedimentary fans in northwestern India. Journal of Hydrology 541(Part B): 1241–1257.

83.

Van Geer

Zuur

(1997) An extension of Box-Jenkins transfer/noise models for spatial interpolation of groundwater head series. Journal of Hydrology 192(5): 339–382.

84.

Von Asmuth

Bierkens

MFP

Maas

(2002) Transfer function noise modeling in continuous time using predefined impulse response functions. Water Resources Research 38(12): 23-1–23-12.

85.

Von Asmuth

Maas

Bakker

, et al. (2008) Modeling time series of ground water head fluctuations subjected to multiple stresses. Ground Water 46(1): 30–40.

86.

Wada

van Beek

LPH

van Kempen

, et al. (2010) Global depletion of groundwater recources. Geophysical Research Letters 37(20): L20402.

87.

Weissmann

Carle

Fogg

(1999) Three-dimensional hydrofacies modeling based on soil surveys and transition probability geostatistics. Water Resources Research 35(6): 1761–1770.

88.

Weissmann

Zhang

Fogg

, et al. (2004) Influence of incised-valley-fill deposits on hydrogeology of a stream-dominated alluvial fan. In: Bridge

Hyndman

(eds) Aquifer Characterization. SEPM, 15–28.

89.

Willis

Tang

(2010) Three-dimensional connectivity of point-bar deposits. Journal of Sedimentary Research 80(5): 440–454.

90.

Yuan

Xie

Liang

(2008) Spatiotemporal prediction of shallow water table depths in continental China. Water Resources Research 44(4): W04414.

91.

Zeng

Xie

Liu

, et al. (2018) Global land surface modeling including lateral groundwater flow. Journal of Advances in Modeling Earth Systems 10(8): 1882–1900.

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.00 MB

6.28 MB