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

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 M0, 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.


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(Allen, , 1984Bridge, 1993;Bridge and Leeder, 1979;Heller and Paola, 1996;Holbrook, 2001;Leeder, 1978;Sheets et al., 2002;Straub et al., 2009;. 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). 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 H i m a l a y a n f r o n t 4 4.5 5 5.5 6 6.5 7 7.5 ). 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.  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.

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, Aqui-Mod, 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 functionnoise (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 geostatistical 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 , 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., 2016aWillis and Tang, 2010).

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 districtaveraged 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.

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 where R i is the value of driver i at time t, y i is the impulse response function of driver i, and t 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 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, y p , takes the form of a Pearson type III distribution function, as used by several previous studies (Von Asmuth et al., 2002Asmuth et al., , 2008 in a variety of hydrogeological settings: where A, b, and n are parameters that define the shape of y p , and 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 where e is the PET time series and y 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., 2015aShapoori et al., , 2015bVon 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. (2015aShapoori et al. ( , 2015b, we used a two-parameter impulse response function based on the Theis solution for the drawdown in a confined aquifer of the form where a and g are calibration parameters. The parameter g is equivalent to 1 4pT with T (L 2 T À1 ) being the transmissivity of the aquifer. The parameter a is defined by r 2 S 4T , 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 a and g.
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 where S y is the specific yield of the aquifer and l is a decay constant. The monthly irrigation driver c is represented by the irrigation rate, I, divided by S y . When l 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 l is used, the contribution of irrigation to h c decays rapidly to zero, in other words (I=S y ) Áy c is virtually 0.

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 c i , 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 l 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 c i =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 , g, a, l, R c i =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 (NSE) (Nash and Sutcliffe, 1970), which is given as where h o and h m are the observed and modeled groundwater levels over n time steps. If NSE ¼ 1, the model is a perfect match to the observations. If NSE ¼ 0 the model is as accurate as the mean of the observed data, and if NSE < 0, it is worse than the observed mean. A value of NSE ¼ 0.2 was set a priori as the threshold between behavioral and non-behavioral model simulations; in other words, those simulations with NSE > 0.2 were deemed to have produced acceptable fits to the observed groundwater level time series. The NSE 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 NSE value, but also to the median parameter values of all acceptable models. 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 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 y 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.

Data analysis
A primary goal of this work is to assess not only the relative importance of different drivers in Example impulse response functions for (a) precipitation y p (with A ¼ 1, so no magnification), (b) abstraction y w , and (c) canal irrigation y c , illustrating the lower and upper limits of the sampled parameter value ranges (Table 1). Note that y e is assumed to be identical in form to y p . 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 y i over time: A large value of M 0 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 0 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 NSE value. We also plotted the median zeroth moment from all of the acceptable solutions with NSE > 0:2 and the coefficient of variation, CV , defined as the standard deviation (s) divided by the mean (m) for those solutions. The zeroth moment for PET stress was not determined because y e was assumed to be a fraction of y p , and so this will show the same spatial pattern. The NSE values give the efficiency of the model outcomes, but we were also interested in the goodness of the model outcome. A disadvantage of the NSE 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 (MAE), which is a more unambiguous measure of the difference between the modeled and observed groundwater level.

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 MAE 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 NSE values were < 0.2 (indicated by the cross-hatch pattern in Figure 3).
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.
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 the areas with groundwater level fall (b and c), indicating that it has no effect on the long-term groundwater level.
component does not increase after 1990 as abstraction rates appear to have stabilized after that time.

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 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.
There is substantial variability in M 0 values for all three drivers, even when we consider only model outcomes with NSE > 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 y w and y h in particular, and are distinct from the location with groundwater level rise. The response function y 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 y 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 0 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.

Spatial variation of the zeroth moment
There are substantial and systematic variations in M 0 for the precipitation, abstraction, and canal irrigation drivers across the study area ( Figure 6). Recall that a high value of M 0 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 MAE (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.
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 The CV 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 CV 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 CV values.
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 0 are explained by the different parameters that determine y, 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 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).
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 a (Figure 7(d)), which is a multiple of the reciprocal of aquifer diffusivity (S=T ), shows no clear spatial pattern. However, the variable g (Figure 7(e)) shows the same spatial pattern as M o;w . High g 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 l, 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.

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 g 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 finegrained 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 grainsize 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.

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 y 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 0 which is a stationary measure of the response to a driver. This, in turn, suggests that scaling the total abstraction or precipitation by M 0 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)).
In contrast, total abstraction (w tot ) 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 Á w tot 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 (g).
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 Monthly canal irrigation, I=S y , for groundwater level rise areas where l is significant and canal irrigation is needed to fit the observations (equation (7)), shows that irrigation goes up to 6 mm/day. 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 firstorder 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 NSE > 0.2, and locations with larger disparities between observed and modeled groundwater levels as indicated by the MAE (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.

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 a. 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 for this article is available online.