Journal of Geophysical Research: Earth Surface Linking the morphology of ﬂuvial fan systems to aquifer stratigraphy in the Sutlej-Yamuna plain of northwest India

The Indo-Gangetic foreland basin has some of the highest rates of groundwater extraction in the world, focused in the states of Punjab and Haryana in northwest India. Any assessment of the eﬀects of extraction on groundwater variation requires understanding of the geometry and sedimentary architecture of the alluvial aquifers, which in turn are set by their geomorphic and depositional setting. To assess the overall architecture of the aquifer system, we used satellite imagery and digital elevation models to map the geomorphology of the Sutlej and Yamuna fan systems, while aquifer geometry was assessed using 243 wells that extend to ∼ 200 m depth. Aquifers formed by sandy channel bodies in the subsurface of the Sutlej and Yamuna fans have a median thickness of 7 and 6 m, respectively, and follow heavy-tailed thickness distributions. These distributions, along with evidence of persistence in aquifer fractions as determined from compensation analysis, indicate persistent reoccupation of channel positions and suggest that the major aquifers consist of stacked, multistoried channel bodies. The percentage of aquifer material in individual boreholes decreases down fan, although the exponent on the aquifer body thickness distribution remains similar, indicating that the total number of aquifer bodies decreases down fan but that individual bodies do not thin appreciably, particularly on the Yamuna fan. The interfan area and the fan marginal zone have thinner aquifers and a lower proportion of aquifer material, even in proximal locations. We conclude that geomorphic setting provides a ﬁrst-order control on the thickness, geometry, and stacking pattern of aquifer bodies across this critical region.


Introduction
Rivers entering sedimentary basins distribute their sediment and water in major sediment fans which have been recognized in stratigraphic records around the world [DeCelles and Cavazza, 1999;Leier et al., 2005;Hartley et al., 2010;Weissmann et al., 2010;Fontana et al., 2014]. The development of alluvial stratigraphy is controlled by river avulsion, sedimentation rate, and the stacking pattern of fluvial channel-belt sand bodies [Leeder, 1978;Allen, 1978;Bridge and Leeder, 1979]. This alluvial stratigraphy, in turn, determines the characteristics and productivity of groundwater aquifers 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 connectivity [Larue and Hovadik, 2006;Renard and Allard, 2013;Flood and Hampson, 2015]. Understanding the stratigraphic architecture of large alluvial aquifer systems is particularly critical because these systems are major repositories for groundwater and are a primary source of fresh water in large parts of the world. Depletion of groundwater resources in alluvial aquifers is now a very significant international problem [Wada et al., 2010], and unsustainable exploitation of groundwater resources requires urgent attention [Gleeson et al., 2010]. We must first understand the spatial pattern and organization of aquifer bodies in order to predict aquifer performance, evolution, and sustainability. It is, however, difficult to do this for most sedimentary basins, due to the very limited subsurface data available in most parts of the world.
A promising way to obtain insights into subsurface stratigraphy and heterogeneity is through an understanding of the geomorphic setting of the aquifer system and the physical constraints that this setting, and the processes that were active during aquifer deposition, place on aquifer body geometry and distribution.
Experimental studies provide insights into the link between sediment transport processes, fan dynamics, and the resulting depositional stratigraphy and large-scale geomorphology of such sediment routing systems [Sheets et al., 2002;Paola et al., 2009;Straub et al., 2009;Van Dijk et al., 2009;Straub et al., 2012]. For example, most laboratory-scale experimental fan deposits fall somewhere between random basin filling that is uninfluenced by topography and purely compensational filling in which deposition always fills topographic lows . The compensation index ( ) is a measure of the relative importance of these different filling patterns . It can be related to the process and frequency of channel avulsion [Sheets et al., 2002]. These experiments provide a framework for understanding likely spatial relationships between channel bodies, but it is not always clear how to link this understanding to field-scale settings. This is because (1) processes and behaviors that are important at experimental scales may not be relevant at the field scale and (2) it is virtually impossible to obtain detailed information on the spatial variations in bed thickness, deposition rates, or avulsion frequency over field length scales of tens or hundreds of kilometers. There is thus a pressing need for analysis of channel body geometry and stacking patterns at these scales, but with a few exceptions [e.g., Rittersbacher et al., 2014;Flood and Hampson, 2015;Owen et al., 2015] this has not been done.
Such an analysis would be particularly useful in northwest India, because it is one of the world's most prominent hot spots of groundwater depletion [Kumar et al., 2006;Rodell et al., 2009;Shah, 2009;Chen et al., 2014]. Groundwater forms the largest supply of irrigation in the states of Punjab and Haryana, which have a combined population of more than 50 million people. The alluvial aquifers in northwest India were deposited by sediment routing systems, dominated by the Sutlej and Yamuna Rivers, which have deposited fluvial sediments in the Indo-Gangetic foreland basin [Geddes, 1960]. Understanding the geometry and evolution of the Sutlej and Yamuna fan systems should therefore give some insight into the spatial distribution of aquifer bodies in the region. Despite this recognition, there have been almost no regional or integrated stratigraphic studies of the aquifer systems in northwest India (see United Nations Development Programme (UNDP) [1985] for an exception), and studies of groundwater dynamics or age have been limited to small spatial scales [e.g., Kumar et al., 2007;Kumar and Gupta, 2010]. Rapid water level decline at the regional scale has been documented by analysis of data from the Gravity Recovery and Climate Experiment [Rodell et al., 2009;Chen et al., 2014]. These studies, however, have very low spatial resolution (1 ∘ by 1 ∘ ) and so cannot be directly related to spatial variability in the aquifer system or used to map detailed patterns of depletion. Thus, it remains unclear (1) how groundwater loss varies in detail across the region, (2) how this variation may relate to geological and geomorphogical heterogeneity in the alluvial aquifer system, and (3) how future changes in groundwater levels might be anticipated and mitigated on the basis of this heterogeneity.
Here we begin to address this urgent societal issue by using a geomorphic framework and available stratigraphic data to understand the large-scale architecture of the aquifer system in northwest India, focusing in particular on the area of the Sutlej and Yamuna Rivers. The objectives of this study are to (i) establish the geomorphic setting of the study area; (ii) explore the degree to which geomorphic setting correlates with, and controls, spatial variability in aquifer properties; and (iii) derive a conceptual model for aquifer body dimensions and how they vary across the region. We first give a detailed description of the study area and describe the methods and data that were used for geomorphological mapping and quantification of aquifer dimensions. Then, we present the geomorphological setting of the region and use that as a framework for analysis of aquifer thickness variations in space and depth. Finally, we develop a conceptual model of aquifer body thickness distribution and fan development in the study region and explore its potential implications for groundwater resources and management. VAN [Yashpal et al., 1980]. Faults are modified from Barnes et al. [2011]: HFT, Himalayan Frontal Thrust; BT, Bilaspur Thrust (BT); MBT, Main Boundary Thrust (MBT). (b) Locations of Central Groundwater Board (CGWB) aquifer thickness logs used in this study. Background is regional topography from SRTM data, with 3 arcsec resolution. Blue shading indicates total depth of the log below ground level. Boreholes for which both aquifer thickness and lithological logs were available are circled. Two representative logs (Pb 100, near the Sutlej River; Hr 579, near the Yamuna River) are labeled and shown in Figure 7a. VAN

Study Area
This study focuses on the area of the Himalayan foreland basin that is fed by the Sutlej River in the west and the Yamuna River in the east (Figure 1a). These rivers have drainage areas of 10,616 km 2 and 10,542 km 2 upstream of the Himalayan mountain front, respectively, and flow into the Indus and Ganga river systems [Sinha et al., 2013]. Uplift and erosion of the Himalaya have resulted in transport and deposition of large volumes of sediment in the Indo-Gangetic basin, but temporal variations in sediment supply and transport capacity have determined the detailed patterns and timing of erosional and depositional events in the Sutlej-Yamuna plain [Goodbred, 2003;Sinha et al., 2005;Gibling et al., 2005Gibling et al., , 2008Roy et al., 2012].
The smaller Ghaggar River drains an area of the Himalayan foothills (485 km 2 ) between the Sutlej and Yamuna catchments (Figure 1). Yashpal et al. [1980] identified a large paleoriver channel that is partly coincident with the location of the modern Ghaggar River. They interpreted the paleochannel, also known as the Ghaggar-Hakra paleochannel, as a former course of the Sutlej, now partly occupied by the underfit modern Ghaggar River. Recent studies have identified sediment deposits in the Ghaggar-Hakra paleochannel that were sourced from the Yamuna and Sutlej catchments [Clift et al., 2012], and geophysical profiles have verified the existence of a large paleochannel within the subsurface [Sinha et al., 2013]. These observations, and the fact that the modern Sutlej and Yamuna Rivers are confined to narrow incised valleys, provide evidence of a complex late Quaternary history of channel avulsion and incision in the Indo-Gangetic plain Tandon et al., 2006;Sinha et al., 2005Sinha et al., , 2007Roy et al., 2012]. Apart from the Ghaggar-Hakra paleochannel, however, further subsurface evidence of former courses of the Sutlej or Yamuna Rivers, or information on the depositional history and subsurface stratigraphy of the Sutlej and Yamuna fans, has not previously been documented.

Methods
This paper evaluates the relationship between the sedimentary deposits of the Sutlej-Yamuna plain, particularly the characteristics of their underlying aquifer bodies, and the geomorphic setting of those deposits.
To establish the geomorphic setting, extents, and dimensions of the major depositional systems, digital elevation models and satellite imagery are used to separate the region into its major constituent geomorphic units, including the major alluvial fans and interfan areas. On the fans, further subdivision is made between inactive fan surfaces and active channel belts, including floodplains, bars, and river channels. Stratigraphic data are then used to relate these geomorphic units with the subsurface stratigraphy and distribution of aquifer bodies.

Remote Sensing Data
To identify geomorphic units, we use mosaics of Landsat 5 and Landsat 8 satellite imagery, along with Shuttle Radar Topography Mission (SRTM) digital elevation data ( Figure 1). A combination of both data sets is needed, as the SRTM lacks the resolution necessary to identify alluvial features, such as abandoned river channels, which are visible on the Landsat images, whereas the Landsat images do not allow discrimination of topographic boundaries between geomorphic units. Both true and false color Landsat images are used to determine drainage patterns and near-surface soil-moisture content. SRTM data are used to distinguish regional patterns of relative elevations associated with different sediment fan units, as well as interfan areas between the major fans.
We use Landsat 8 Operational Land Imager (OLI) data to map active and abandoned channels within our study area. Prior work has used Landsat 4 Multi-Spectral Scanner (MSS) satellite images for mapping a major paleochannel on the Sutlej-Yamuna plain [Yashpal et al., 1980]. We combine nine individual OLI scenes, acquired between November and December 2013, to produce a relatively seamless color composite mosaic and to map channel features at a higher spatial resolution. Timing of image acquisition is critical to mapping ability, as vegetation cover should ideally be kept to a minimum. Imagery acquired just after the monsoon season is particularly useful because inundation of flood waters is affected by soil composition and surface topography. The visible bands (two blue, three green, and four red) are badly affected by atmospheric scattering (haze) so that true color and standard false color composites lack visual clarity and are difficult to interpret; we have therefore mainly used bands 5, 6, 7, and 10 for our analyses. It is well known that moisture content depresses the overall reflectance of soils and rocks [e.g., Price, 1990;Lobell and Asner, 2002], especially in the near and short-wave infrared (bands 5 and 6) and to a lesser extent in band 7. The Tasselled Cap Transform [Crist and Cicone, 1984]  Landsat bands. This combination reveals that the sediments in the Ghaggar-Hakra paleochannel are less reflective (darker) and wetter than the surrounding sediments and that these effects are not caused by the presence of vegetation. We interpret this to indicate that sediments within the paleochannel have higher moisture content and are less well drained than those outside. Our work also reveals that a color composite of bands 5, 6, and 10 (RGB, referred to below as 5610) best exploits this effect on the relative brightness of alluvial materials in this area. The thermal infrared (band 10) has lower reflectance for the wetter soils [Price, 1980;Wang and Qu, 2009], resulting in a dark blue color in Figure 1a. Dry regions are shown as yellow because of the high reflectance in bands 5 and 6, while the Thar Desert appears almost white due to high reflectance in band 10 ( Figure 1a). The margins of the Ghaggar-Hakra paleochannel have higher elevations and appear brighter in bands 5 and 6, giving them a lighter pale blue color (high reflectance in both red and green) against the dark blue tones of the paleochannel.
In addition, we use true color (bands 2, 3, and 4 RGB, referred to below as 234) Landsat 8 imagery from both premonsoon and postmonsoon periods to map channel, bar, and floodplain features of the active Sutlej and Yamuna Rivers. The floodplain is the area of land between the active river banks and the base of the valley walls and experiences flooding only during periods of high discharge. The active channel is the position of the modern channel of the Sutlej and Yamuna Rivers identified from Landsat 8 imagery. The channel bars are mapped as areas of bare sediment along the active channel, likely due to yearly flood inundation. For both the Sutlej and Yamuna, we distinguish between the active channel and channel belt of the total fluvial corridor or incised valley and measure the widths of both features at multiple locations to get a range of channel system widths across the study area.
To identify different geomorphic units, we use a subset of NASA's global Shuttle Radar Topography Mission (SRTM) elevation data with a base resolution of 1 arc-seconds (about 30 m). To reduce the noise in the data in low relief foreland areas, we apply median filtering with a window size of 3 by 3 pixels to the data and determine flow paths automatically in Matlab using the Topo-Toolbox 2 [Schwanghart and Scherler, 2014]. The flow paths are used to identify the river channels such as the Ghaggar River that drain the Himalayas but are not identified from the Landsat 8 (234 RGB) imagery.
The Sutlej and Yamuna fans are identified from the SRTM data by extracting concentric elevation profiles that are centered on the points at which the rivers exit the Himalayan Mountains and enter the alluvial basin. These profiles show quasi-uniform elevations (Figures 2a and 2b), indicating near-conical fan shapes. The Sutlej fan shows a fairly uniform gradient with distance from the apex, whereas the Yamuna fan is slightly concave up ( Figure 2b). The conical fan shapes imply that the locus of active deposition has shifted over time due to repeated migration or avulsion of the channel system.

Aquifer Thickness Data
In order to understand the bulk sedimentary architecture and aquifer geometry, we use aquifer thickness logs obtained from the Central Groundwater Board (CGWB). The data set consists of the thicknesses of aquifer and nonaquifer units interpreted by the CGWB from the electrical logs taken from each borehole. The depth of the logs varies between 50 m and 500 m (Figure 1b), but 90% are at least 200 m deep; here we restrict our statistical analysis to the top 200 m of each log and discard those records that did not reach that depth to maximize the data coverage, leaving us with 243 logs. We also obtained 12 CGWB boreholes for which both aquifer thickness logs and lithological logs are available (indicated in Figure 1b), allowing direct comparison of the two data sets and enabling us to understand the relationship between aquifer units and actual subsurface stratigraphy. The lithological logs contain a description of the drill cuttings returned by a rotary bit at regular intervals (around 3-4 m) or where there is a notable change in formation, and are classified into clay, silt, sand, and gravel. Aquifer units are inferred from the lithological logs by classifying fine-coarse sand and gravel as aquifer material, while silt and clay were classified as nonaquifer material.

Aquifer Distribution
Understanding the spatial distributions of both aquifer body thickness and bulk aquifer percentage is essential for determining the likelihood of finding aquifer bodies of a given thickness in the subsurface and for understanding how aquifer thicknesses vary across different geomorphic units. Also, because of the grain-size difference between aquifer and nonaquifer layers, the bulk percentage of aquifer bodies is related to the overall specific yield of the subsurface [e.g., Johnson, 1967;Robson, 1993]. Compaction and dewatering of VAN DIJK ET AL. LINKING FAN MORPHOLOGY TO STRATIGRAPHY 205 nonaquifer layers may also affect the spatial subsidence rate associated with pumping [Higgins et al., 2014]. We analyze the bulk percentage of aquifer material within the top 200 m of the CGWB aquifer thickness logs and look at the spatial variability in aquifer percentage both within and between different geomorphic units-that is, between the fan surfaces, the interfan area between the fan heads, and the marginal zone along the boundary between the two fans. A two-sample t test is used to determine whether the mean aquifer percentages between different geomorphic units are equivalent. VAN  To quantitatively compare aquifer thickness patterns across space and depth, we compile the exceedance probability of aquifer thickness data-that is, the probability of finding an aquifer unit of at least a given thickness-for the entire region, for different geomorphic units, for varying distances from the fan apices, and from different depth intervals. The exceedance probability, or complementary cumulative distribution function, is more robust than a probability density function against fluctuations due to finite sample size, particularly in the tail of the distribution (as suggested by Clauset et al. [2009]). We apply the maximum likelihood methods of Clauset et al. [2009] on aquifer thickness data, on 1 m bin intervals, to evaluate the likelihood that the aquifer body thicknesses follow a heavy-tailed distribution and where appropriate to fit a power law function to the tail of the distribution. We calculate a p value, which indicates if the power law hypothesis is a plausible one for the aquifer body thickness data and assume that power law behavior can be ruled out if p ≤ 0.1. The exceedance probability asymptotes to 1 as x approaches zero, so that the power law behavior cannot hold for all x ≥ 0 and there must be some lower bound, x min , to the regime of potential power law behavior. Here we focus our attention on the tail of the distribution, as it gives an indication of the likelihood of finding thick aquifers within the subsurface. The tail is described by a truncation value or lower bound, x min , and a slope or scaling parameter, . We also compile the exceedance probabilities of both the aquifer body thickness data and the bed thicknesses from the full depth extent of the 12 boreholes, in order to quantitatively understand the relationship between the two data sets.

Aquifer Persistence Analysis
The sediment filling or stacking pattern determines the spatial persistence of the channel system over time, or equivalently its propensity to occupy different parts of the basin. It is thus ideal for assessing the degree to which individual aquifer bodies are stacked vertically during deposition and thus the likelihood of vertical connectivity between those bodies. Straub et al. [2009] defined compensational stacking or filling by the time scale over which the sediment routing system occupies every spot in the basin to "compensate" the subsidence. This time scale can be identified by examining the standard deviation of sedimentation rate over subsidence rate ( ss ): where r(T; x, y) is the local sedimentation rate measured over a stratigraphic time difference T, x and y are horizontal coordinates, A is area measured parallel to stratal surfaces, andr is the long-term average sedimentation rate. The value of ss approaches zero for increasingly large time intervals, over which subsidence must eventually balance deposition. Straub et al. [2009] showed that this decay of ss with increasing time interval T is expected to follow a power law function of the time window T, with the compensation index ( ) defined as the power law exponent: where a and are empirical coefficients. A compensation index of 1.0 indicates that the deposits stack in a purely compensational manner, meaning that the depocenter shifts progressively to fill the lowest point in the basin and sedimentation rates rapidly approach the long-term subsidence rate over increasing time intervals Wang et al., 2011;Hajek and Wolinsky, 2012]. In contrast, a compensation index of 0.5 indicates random filling of the basin that is uncorrelated in time, and an index of 0 indicates perfect anticompensation-in other words, persistence of the channel along a single corridor through time. The compensation index is thus a measure of the tendency for channels to stack along one or several preferred channel pathways.
Here we adopt a modified version of equation (1), because we lack the stratigraphic and age data required to reconstruct true depositional thickness and sedimentation rates over time within the Sutlej-Yamuna region. Instead, we are interested in the persistence of channel deposits and thus their potential for vertical connectivity. We assume that aquifer units in the CGWB aquifer thickness logs are likely to represent either individual or amalgamated channel deposits and can therefore be treated in the same way as distinct beds-recognizing that a single aquifer unit may be composed of one or several different beds. By analogy with Straub et al.
[2009], we examine the standard deviation of the fraction of aquifer material (f ) over progressively larger stratigraphic thickness intervals (D). We expect channel persistence to be shown by values of f that are relatively uniform over particular thickness ranges ( ∼0). We define the standard deviation of the aquifer fraction ( f ) at a single point as wheref is the average fraction of aquifer material in a single aquifer thickness log and B is the stratigraphic thickness. Instead of calculating f along a transect for different stratigraphic intervals, as in Straub et al. [2009] (equation (1)), the aquifer fraction is calculated within individual logs for different thickness intervals (D) ranging from 1 m to 100 m (with logarithmic bin intervals). As before, we limit our analysis to the top 200 m of the aquifer thickness logs and divide the available logs by geomorphic unit into the Sutlej fan, Yamuna fan, and interfan areas. The value of f approaches zero for increasing stratigraphic thicknesses, as the aquifer fraction approaches the average aquifer fraction for that log. Again by analogy with Straub et al. [2009], we observe that f decays as a power law with increasing D, with a power law exponent that defines the aquifer persistence index f : where a f and f are empirical coefficients and f is analogous to the compensation index of Straub et al. [2009].
We plot median f values against D for the Sutlej fan, Yamuna fan, and interfan area. Random, uncorrelated thicknesses of aquifer units should result in f decreasing as the square root of stratigraphic thickness for increasing thickness intervals, i.e., f = 0.5 ]. If f is less than 0.5, then the f is relatively independent of stratigraphic interval, indicating persistence of the aquifer fraction (although note that local values of f can still be quite different from the overall borehole average). If f is greater than 0.5, then the standard deviation decreases rapidly with increasing stratigraphic interval, approaching the overall borehole average.

Sediment Routing Systems and Geomorphology
Observations of Landsat imagery and the DEM enable us to distinguish the major sediment routing systems and their deposits ( Figure 3a). Broadly, the region comprises two major sediment fan systems associated with the Sutlej and Yamuna Rivers (as originally identified by Geddes [1960]), separated by an interfan area. These fans are bounded by the faults of the Himalayan Frontal Thrust (HFT) to the northwest and by the deposits of the Thar desert and crystalline bedrock of the Indian craton to the southwest and south, respectively ( Figure 3). Most of the current surface area of the fans is disconnected from the Sutlej and Yamuna Rivers, as both rivers flow within incised valleys that are cut into older fan deposits. At their distal margins, about 250 km from the Sutlej fan apex and 200 km from the Yamuna fan apex, the fan surfaces are covered by dune deposits of the Thar Desert. Perpendicular to the mountain front, the slope of the Sutlej fan decreases from 0.066% near the apex to 0.027% at 150 km from the mountain front, whereas the slope of the Yamuna fan decreases from 0.057% near the apex to 0.017% at 150 km from the mountain front.
The surfaces of both the Sutlej and Yamuna fans show elongated, discontinuous ridges oriented northeast to southwest, especially in proximal and medial areas of the fan (Figure 3a). The ridges are 10-100 km long and 650-2300 m wide (Table 1) and show local relief of up to 5 m. The ridges appear to radiate from the fan apices and are largely coincident with higher relative reflectance (i.e., low soil moisture content) zones visible on the Landsat 8 5610 (RGB) mosaic and on the Landsat 5 image of bands 5, 3 and 1 (RGB, Figures 3b-3d). The elevated topography, radial distribution about the fan apices, and low moisture content of these features lead us to interpret them as abandoned, sand-rich paleochannel deposits, preserved on the surfaces of both fans. Similar features have been noted in other alluvial channel belts and have been ascribed to older channel deposits that are picked out by differences in sediment grain size, leading to variable compaction and subsidence [e.g., Berendsen and Volleberg, 2007]. They are also observed on other fan surfaces of the Ganga sediment routing system, where they have been interpreted as paleoriver channels that are later infilled by eolian sediments after abandonment [Srivastava et al., 2000;Gibling et al., 2005]. They are potentially very useful as analogues for buried channel bodies within the Sutlej and Yamuna fans, whose dimensions are much harder to constrain. These inferred paleochannel locations should, however, be tested in the field with lithological data to determine if the deposits are fluvial or eolian.
Between the conical fan surfaces lies an interfan area of c. 4000 km 2 that occupies the region adjacent to the mountain front. It is characterized by smaller river channels compared to the Sutlej and Yamuna Rivers, and lacks the elongate ridges or other surficial evidence of paleochannel positions found on the fans. The boundary of the interfan area is determined by the Landsat 8 image as well as the DEM. On the Landsat 8 5610 (RGB) mosaic, the interfan is characterized by relatively high, and uniform, soil moisture. The interfan area is relatively high compared to the Sutlej and Yamuna fan surfaces, and is planar rather than conical, as shown by elevation contours that are parallel to the Himalayan mountain front (Figure 2a).
VAN DIJK ET AL. LINKING FAN MORPHOLOGY TO STRATIGRAPHY 208   Yashpal et al. [1980], is characterized by low reflectance and thus high soil moisture content. The paleochannel is about 5000-8000 m wide, while the present-day Ghaggar River is only 60-100 m wide (Table 1).
The dimensions of these channel features visible on the fan surfaces, including the incised valleys and ridges, are important, because they illustrate the typical widths of recent channel deposits in these sediment routing systems, and provide a first-order constraint on the dimensions of older channel bodies within the subsurface. The width of the paleochannel ridges may be more appropriate analogues to use than the incised valley dimensions, as the ridges were formed under net aggradational conditions on the fan rather than reflecting the dimensions of the sediment routing system during incision and excavation. On the other hand, the presence of incised active and inactive channels indicates that at least some of the buried channel bodies in the Sutlej and Yamuna fan systems are likely to consist of incised valley fills.

Subsurface Architecture
In this section, we quantify spatial variations in the dimensions and the persistence of the aquifer bodies across the fan surfaces and within the different geomorphic units (Sutlej fan, Yamuna fan, and interfan area). Because we lack detailed subsurface data around the boundary between the Sutlej and Yamuna fans, we assume for simplicity that the surface boundary between the two fan systems has persisted throughout deposition of the upper 200 m of sediment. It is certain that this boundary must have shifted over time, leading to interfingering between Sutlej and Yamuna fan deposits, but at the moment we are unable to quantify the extent of this variability.

Percentage of Aquifer Bodies
The mean percentage of aquifer bodies across all CGWB aquifer thickness logs is 39%, but values for individual logs range from 0% to 80% (Figure 4), with major variations between adjacent wells (Figure 3a). The percentage of aquifer bodies within each fan body does not noticeably vary laterally, although a general down fan decrease in aquifer percentage is observed in both fans (Figure 3a). In contrast, the interfan area and the fan marginal area at the boundary between the Sutlej and Yamuna fans (Figure 3a) both show lower percentages of aquifer bodies compared to the fans themselves (Figures 3 and 4 and Table 2), especially in the deeper parts of the section. Two-sample t tests show that the mean aquifer body percentages of the Sutlej and Yamuna fans are indistinguishable (p = 0.97) but that both are significantly larger than the mean aquifer body percentage in the interfan area and fan marginal area (p < 0.05). There is a small decrease of the mean percentage of aquifer bodies width depth within both fans (Table 3), except for the top 50 m.
To illustrate the fan-scale variability in aquifer body thickness and depth, we compile two representative transects of aquifer thickness logs at medial and distal positions down fan (Figures 5 and 7a). There is no clear relationship visible between aquifer body thickness and depth for adjacent logs and no evidence that aquifer bodies are laterally connected or correlatable at the length scale of the log spacing (median ∼7000 m). This result is perhaps not surprising, as this median log spacing is larger than the widths of the channel features identified on the Sutlej and Yamuna fan surfaces ( Figure 5)  (Table 3). Aquifer body thickness varies across both transects, but most aquifer bodies are less than 10 m thick. Because of this lack of spatial correlation, we focus our analysis on statistical descriptions of the spatial variability of aquifer body thickness.

Spatial Variability in Aquifer Thickness Distributions
The mean thickness of aquifer bodies across the study area is about 6 m, with individual values that range between 1 and 100 m. A two-sample t test shows that mean aquifer body thicknesses from the two fans are similar (p = 0.14). In contrast, the mean thickness of aquifer bodies in the interfan area and fan marginal area is less than that of the fans (p < 0.05). The median aquifer body thickness of the fan marginal area, however, is similar to that of both fans ( Table 2).
The aquifer thickness data from all geomorphic units (both fans and the interfan area) are well characterized by heavy-tailed exceedance probability distributions using the criteria of Clauset et al. [2009], with values of p > 0.1 indicating that heavy-tailed behavior cannot be ruled out [Clauset et al., 2009] (Table 2 and Figure 6a). The x min value is comparable for the two fans, but for the Sutlej aquifer units is larger, meaning that it is  somewhat less likely to find aquifer bodies thicker than 17 m (x min ) in the deposits of the Sutlej fan compared to the Yamuna fan. The interfan area has a comparable value as the Yamuna fan, but the x min value is lower, meaning that there are fewer thick aquifer bodies. Aquifer bodies from the fan marginal area do not follow a heavy-tailed distribution, and the data in Figure 6 show that there are fewer thick aquifer bodies in the fan marginal area compared to the interfan or the fans themselves.
The variations in aquifer body thickness distributions measured over different depth intervals and distance from the fan apices give an indication of potential changes in depositional characteristics of the Sutlej and Yamuna fan systems over time and space. In general, the distributions of aquifer body thickness for different depths and at different distances are comparable, as both the and x min values are relatively invariant with distance from the apex as well as depth below the surface for most intervals (Table 3 and Figure 6b). This is also observed in the quantiles of the aquifer body distribution (25th, 50th, and 75th) as these remain relatively invariant for different depth or distance intervals. Some intervals, however, have a lower value, meaning that thicker aquifer bodies should be more frequent, but these intervals typically also have a lower x min value which offsets this trend. Although the distribution of aquifer body thickness does not change appreciably with distance from the apex, the overall aquifer body percentage does decrease down fan for both the Sutlej and Yamuan fans (Table 3). These findings indicate that while aquifer bodies are less common in the distal parts of the fan systems, those bodies that are present follow similar thickness distributions as seen in more proximal locations. In other words, aquifer bodies are less common in distal settings, but if found are just as likely to be of at least a given thickness as in proximal parts of the system.

Accuracy of Aquifer Body Thickness Data
Because the CGWB aquifer thickness data are interpreted from geophysical (electrical) logs rather than from lithological information, it is important to establish the relationship between the aquifer body thicknesses and their constituent lithologies. Cross comparison of aquifer body thicknesses derived from electrical logs with the lithological logs for the 12 boreholes where both records are available shows that aquifer units VAN  Both the Sutlej and Yamuna fans are characterized by aquifer-rich and aquifer-poor zones. In both panels, the lack of correlation between adjacent wells in both transects, even where they are closely spaced, argues for limited lateral dimensions of channel bodies, as expected in a fan sediment routing system. generally correspond to material that is recorded as fine-grained sand or coarser, while nonaquifer units generally correspond to silt and clay (Figure 7a). This relationship is not always consistent; in particular, units within the top 20 m are often recorded as nonaquifer material by the CGWB. To assess the effects of the relationship on our aquifer body thickness distributions and thus on the potential uncertainty in our statistical descriptions of aquifer body thickness, we classified the 12 available lithological logs into aquifer (fine-grained sand and coarser) and nonaquifer (silt and clay) units. This yielded a total of 101 distinct lithology-based aquifer units in the 12 boreholes, compared to 146 geophysically based aquifer units in the corresponding aquifer thickness logs. Comparison of the exceedance probability distributions of these two different aquifer data sets shows that the geophysically based aquifer bodies are slightly thinner compared to those derived from lithological data (Figure 7b). Thus, the "true" aquifer bodies in the study area are likely to be slightly thicker, but less numerous, than indicated by the CGWB aquifer thickness data, and our analysis of aquifer body thickness distributions is thus slightly conservative. Encouragingly, the mean percentage of aquifer bodies in the two data sets is essentially identical (38% in the geophysically based aquifer thickness data, 39% in the lithology based data).
VAN DIJK ET AL. LINKING FAN MORPHOLOGY TO STRATIGRAPHY 213 Figure 6. Exceedance probability curves of aquifer body thickness for each geomorphological unit, separated into the (a) Sutlej and Yamuna fans, the interfan area, and the fan marginal area, and exceedance probabilities of aquifer body thickness for the (b) proximal and distal parts of the fans. Dashed lines show best fit heavy-tailed distributions as determined by maximum likelihood [Clauset et al., 2009], along with the corresponding value of the scaling exponent . Solid vertical lines show the median (50th percentile) thicknesses for each distribution. Line color is tied to symbol color for each unit. Note in Figure 6a that aquifer body thicknesses for the interfan and fan marginal area are consistently smaller than those in the two fans. Thicknesses in the fan marginal area deviate substantially from a heavy-tailed distribution, with a p value of 0.06 indicating that such a distribution is unlikely [Clauset et al., 2009]. Note in Figure 6b that aquifer body thicknesses for the distal part of the Sutlej fan are slightly thinner than for the proximal part, but that both parts of the Yamuna fan have similar probabilities.

Aquifer Persistence Analysis
For aquifer bodies underlying all geomorphic units, the standard deviation of aquifer fraction f is approximately independent of the stratigraphic thickness D for small thickness intervals, and decays with increasing D. For small D, D is either dominated by aquifer or nonaquifer bodies and deviates the most with the mean aquifer fractionf . f increases monotonically from 0 to 1.0 with increasing D, which means that the aquifer fraction distribution changes from persistent stacking of aquifer units to a more random stacking pattern  (Figure 8). These values indicate that the threshold for f = 0.5 is around twice the median aquifer body thickness for the fans but around 4 times the median aquifer body thickness for the interfan area. The threshold for f = 1.0 is around 5 times the median aquifer body thickness for the fans and as much as 8 times the median thickness for the interfan area. Alternatively, the threshold for f = 0.5 is approximately equal to the 75th percentile of aquifer body thickness for the fans, and that for f = 1.0 is around 3 times the 75th percentile. These results indicate that the interfan area consistently shows more persistent, less compensated behavior, and that aquifer fraction must be averaged over greater stratigraphic thicknesses in the interfan area in order to observe the onset of compensational behavior.

Discussion
This study provides the first regional view of the spatial distribution and statistics of aquifer bodies in the subsurface of the Indo-Gangetic basin in northwest India. Importantly, our results show a generic link between aquifer body dimensions and distribution and geomorphic setting across the Sutlej-Yamuna plain. This means that separation of the surface geomorphology into sedimentary fans and interfan areas provides a first-order framework for understanding and therefore predicting aquifer body geometry and thickness variations. Below, we discuss how our observations fit within this framework of fan construction and alluvial aquifer VAN DIJK ET AL. LINKING FAN MORPHOLOGY TO STRATIGRAPHY 214  Figure 1b. For each borehole, the left-hand panel shows the lithological log as determined from drill cuttings, while the right-hand panel shows aquifer and nonaquifer units inferred from the geophysical log by CGWB. Kankar refers to carbonate nodules formed by pedogenetic processes or groundwater precipitation [Sinha et al., 2007]. For well Haryana 579, aquifer units generally correspond to fine-coarse sand or gravel beds, while nonaquifer units correspond to silt and clay layers; the main exceptions to this occur in the upper 20 m of the well, which has been interpreted as nonaquifer material by CGWB regardless of grain size. For well Punjab 100, most fine-medium sand layers correspond to aquifer units, but there are several exceptions to this rule. Note that the thickness of individual aquifer units in Punjab 100 is often less than the thickness of contiguous sand beds in the lithological log. (b) Comparison of the exceedance probability curves of aquifer body thickness from the aquifer thickness logs (black symbols) and thickness inferred from the lithological logs (grey symbols) for the 12 logs. Dashed vertical lines show the quartile thicknesses of each data set; line color is tied to symbol color. Aquifer bodies extracted from the CGWB aquifer thickness logs are consistently slightly thinner than those inferred from the lithological logs, meaning that the distributions and scaling relationships in Figure 6 are slightly conservative in terms of true aquifer body thickness.
VAN DIJK ET AL. LINKING FAN MORPHOLOGY TO STRATIGRAPHY 215 stratigraphy. We also compare our results to those of other studies that have characterized the statistics of fluvial channel bodies, discuss the hydrogeological implications of our key observations, and consider the major remaining gaps in our understanding of the northwest Indian aquifer system.

Link Between the Morphology and Stratigraphy of the Fan Aquifer System
The Sutlej and Yamuna sediment routing systems form a pair of laterally interacting fans within the Himalayan foreland basin [Geddes, 1960]. This leads to a conceptual model of fan morphology and stratigraphy that has some useful implications for interpreting their stratigraphic architecture and thus for understanding aquifer geometry. Here we link the results of our statistical analysis on aquifer distribution with the overall construction and architecture of the fan systems, illustrated in Figure 9.
Fluvial fans are deposited by channel systems that radiate downslope from the fan apex, such that water and sediment are distributed over a conical space but follow different transport pathways over time (Figure 9a). This means that individual channel deposits are likely to form elongate sand bodies that are highly longitudinally connected (in the down fan direction) but are less connected in the lateral direction. The aquifer thickness logs from our study area show that, consistent with this expectation, individual aquifer bodies cannot be correlated laterally between adjacent wells with a median spacing of ∼7 km ( Figure 5) and must therefore be narrower than this, on average. It is not possible, with our available data, to determine the widths of the aquifer bodies more precisely, but we can place some approximate constraints on likely aquifer body widths using (1) detailed characterization of the Ghaggar-Hakra paleochannel in a few locations, (2) observations of active and relict channel-belt widths from these and other fan surfaces, and (3) channel body thickness-width scaling relationships [e.g., Gibling, 2006]. Sinha et al. [2013] used coring and resistivity soundings to infer the presence of a composite sand body below the Ghaggar-Hakra paleochannel, with a width of >12 km. They interpreted this body as the amalgamation of multiple individual fluvial channel bodies deposited by a large river flowing along the paleochannel axis. Channel-belt widths of modern Sutlej and Yamuna Rivers show typical widths of up to 5 km (Table 1), while the ridges associated with aggradational paleochannel deposits on the fan surfaces are up to 2.3 km wide. Abandoned paleochannels on the Tista megafan in the eastern Ganga Basin show widths of up to 3.3 km [Chakraborty and Ghosh, 2010]. Finally, empirical relationships between channel body thickness and width [Gibling, 2006] show a common width-to-depth range of 30-1000, which means that the median aquifer body thickness of 6 m should correspond to a width of up to 6 km. Together, these disparate observations all suggest that maximum across-strike channel body widths in this setting are likely to be no more than ∼5-10 km, consistent with the lack of lateral correlation between our aquifer thickness logs along the medial and distal transects ( Figure 5). This upper limit imposes an inherent lateral length scale into the system which may influence hydrogeological connectivity and flow paths.
VAN DIJK ET AL. LINKING FAN MORPHOLOGY TO STRATIGRAPHY 216  (Table 3), and aquifer bodies may meander out of the plane of section. (e) Hypothetical EP curves for aquifer body thickness showing potential variations in the cross fan direction. (f ) For the same overall proportion of aquifer material, an exponential or thin-tailed distribution (green triangle) would yield a very low probability of finding thick aquifer units, implying discrete or single-storied aquifer bodies-perhaps due to frequent avulsions and compensational stacking. (g) In contrast, a power law or heavy-tailed distribution (orange triangle) would suggest a greater probability of finding very thick aquifer bodies, perhaps due to stacking of multistoried channel deposits or filling of incised valleys. Data from the Sutlej and Yamuna fans are consistent with the latter model, implying locally high vertical connectivity.
Down fan trends in aquifer percentage and aquifer body thickness distribution can also be understood in relation to the construction of these fan depositional systems. We observe that the scaling exponent on the thickness distribution is essentially uniform with distance from the fan apex but that the percentage of aquifer material decreases down fan. These results indicate little or no down fan decrease in aquifer body thickness; instead, the dominant variation in the down fan direction is a decrease in aquifer body volume as a proportion of overall fan sediment volume, which can be understood as a simple volumetric consequence of the conical fan shape. Rivers on fans are typically characterized by a distributive drainage system and thus lose or maintain rather than gain water and sediment discharge down fan [e.g., Nichols and Fisher, 2007;Weissmann et al., 2010;Hartley et al., 2010;Weissmann et al., 2015]. The near-uniform value on the thickness distributions is consistent with little down fan variation in water and sediment discharge during channel body deposition (Table 3 and Figures 9b-9d)-not surprising, given the relatively short length scales of the fan systems compared to total catchment sizes. We see no evidence in our aquifer body thickness distributions for regional down fan thinning or "feathering" of the aquifer bodies [e.g., UNDP, 1985].
The geomorphic distinction between fan and interfan settings also introduces an important large-scale lateral heterogeneity. Aquifer thickness data from the interfan area show that the aquifer bodies are consistently VAN DIJK ET AL.
LINKING FAN MORPHOLOGY TO STRATIGRAPHY 217 thinner than those in the fans and make up a smaller proportion of the upper 200 m, even close to the mountain front. This is because the interfan area is not fed by a major Himalayan sediment routing system. Because of this lateral heterogeneity in aquifer body dimensions, it is not possible to simply use proximity to the mountain front as a proxy for key aquifer properties, such as grain size or channel body thickness; knowledge of the geomorphic setting and proximity to major sediment entry points is required as well. We note that the variation in aquifer body percentage between the fan areas and interfan area documented in Figure 3a provides a close match to spatial variability in specific yield values tabulated by UNDP [1985], although that study did not provide an explanation for the observed patterns. It remains unclear, however, whether the lower specific yield values in the interfan area are the result of finer overall grain sizes or more poorly sorted material.
Our results also shed some light on channel body stacking patterns across the Sutlej and Yamuna fans. Aquifer body thickness and vertical connectivity will be strongly controlled by the channel-stacking pattern, which in turn results from the competition between avulsion rate and sedimentation rate [Bryant et al., 1995;Mackey and Bridge, 1995] and channel reoccupation [Stouthamer, 2005]. Our analysis shows that a transition to approximately random aquifer body stacking ( = 0.5) occurs over stratigraphic thicknesses that are approximately equal to the 75th percentile of aquifer body thickness and that the aquifer fraction approaches the borehole average value-indicating compensational behavior-beyond about 3 times the 75th percentile ( Figure 8 and Table 2). We interpret these results as indicating relative persistence of aquifer bodies over thickness intervals that are less than about ∼35 m on the Sutlej and Yamuna fans and impersistence over larger intervals. For example, if the upper 35 m of a borehole log is dominated by aquifer units, then the lower portion of the log is likely to be dominated by nonaquifer units in order to maintain a typical mean aquifer fraction f of ∼0.4. This break in aquifer thickness scaling behavior is reminiscent of that documented by Wang et al. [2011], who showed that full compensation in a section of clustered channel deposits occurred only over a stratigraphic interval of at least 4 times greater than the maximum channel body thickness.
While these results are necessarily tentative because of the limitations of our aquifer thickness data, we interpret them as indicating that over short time scales, locally persistent occupation of a single channel corridor can allow the deposition of thick aquifer units, leading to the heavy-tailed aquifer thickness distributions that we observe across the study area. These thick units are likely to represent stacked, multistoried channel bodies with thicknesses that are a multiple of the median aquifer body thickness (Figure 9g). In contrast, if the study area was dominated by simple or single-story channel deposits (Figure 9f ), then we would expect less evidence of local persistence and a thinner-tailed aquifer thickness distribution (Figure 9e). Chamberlin and Hajek [2015] showed that multistoried sand bodies are more likely to occur under conditions of persistent or random filling, rather than pure compensational stacking. Importantly, however, even these persistent aquifer bodies are limited in their total thickness, as we do not observe individual aquifer bodies that are >100 m thick. We infer that, on short time scales, the fan systems may have been dominated by local avulsions that allowed the construction of thick aquifer units composed of stacked channel deposits. Over longer time scales, however, larger-scale or regional avulsions have shifted the channel away into different depositional corridors. One way of creating these corridors is through the formation and subsequent filling of incised valleys across the fan surface [Weissmann et al., 2002;Fontana et al., 2008]. The Ghaggar-Hakra paleochannel represents a filled, abandoned incised valley, whereas the modern Sutlej and Yamuna valleys have incised but are not yet filled.
Overall, this conceptual model provides a plausible explanation for the occurrence of widespread, relatively thick aquifer units, as indicated by the heavy-tailed aquifer thickness distributions (Figure 6), without recourse to channels and thus channel deposits that are much larger than those that are active at the present day.

Hydrogeological Implications
The inferences about fluvial fan stratigraphy and fan architecture that we draw from our geomorphic and aquifer thickness observations are useful for understanding the hydrogeology of the Indo-Gangetic basin aquifer system in northwest India. Most critically, the aquifer bodies in the CGWB database appear to be dominated by sand-rich deposits that were deposited by the river systems that built the Sutlej and Yamuna fans, along with smaller distributive rivers across the fans and in the interfan area. By analogy with the modern Sutlej and Yamuna River systems, these deposits are continuous down fan but highly laterally discontinuous. We expect, therefore, that bulk hydraulic properties of the aquifer system should be strongly anisotropic [e.g., Anderson, 1989;Fogg et al., 2000] Yamuna fans (Figure 6b), although they make up a smaller proportion of the subsurface in distal settings ( Figure 5). These thick aquifer bodies are comprised of stacked, multistoried fluvial channel deposits, and we expect that vertical connectivity (and thus hydraulic conductivity) within such deposits should be locally high [e.g., Weissmann et al., 2004;Larue and Hovadik, 2006;Renard and Allard, 2013], especially in areas with low f values. Importantly, along-strike geomorphological variations between fan and interfan settings are closely correlated with differences in bulk aquifer percentage and in the statistical distribution of aquifer body thicknesses, as well as with independently compiled estimates of specific yield [UNDP, 1985]. Thus, simple proximity to the mountain front appears to be a poor predictor of aquifer properties. We suggest that assessment of across-strike aquifer variability is important and should account for position relative to major sediment entry points into the Himalayan foreland [Gupta, 1997].
The spatial variations in aquifer percentage and aquifer body thickness that we document here indicate that a laterally-uniform, "layer-cake" hydrogeological model is not applicable in fluvial fan systems like the Sutlej-Yamuna plain, as noted by previous workers [e.g., Fogg, 1986;Koltermann and Gorelick, 1996;Fontana et al., 2008Fontana et al., , 2014. The types of lateral and vertical heterogeneity that characterize fan systems, including variations in grain size, porosity, mineralogy, lithologic texture, and channel body structure, will cause variations in hydraulic conductivity, storage, and porosity, and thus control flow and transport through the subsurface [Fogg, 1986;Koltermann and Gorelick, 1996;Eaton, 2006]. Other studies of channel body aquifers have pointed out that ignoring the connectivity of permeable but spatially distinct channel deposits limits the ability to perform appropriate hydrogeological analysis [Anderson, 1989;Fogg et al., 2000;Burns et al., 2010;Van der Kamp and Maathuis, 2012]. Renard and Allard [2013] showed that connectivity is a key influence on a wide range of groundwater flow and transport processes but is most important in areas with moderate proportions of aquifer bodies. As our study area contains a bulk aquifer fraction of about 40%, the arrangement of aquifer bodies should be considered in future hydrogeological modelling of our study region.
Promisingly, however, we have shown that important characteristics of the aquifer system, including the percentage of aquifer bodies, the distribution of aquifer body thickness, and the stacking patterns of individual aquifer bodies, vary in systematic ways between fan and interfan geomorphic units. This raises the possibility that lateral variations in geomorphic setting within active fluvial fan systems-which can be easily assessed from surface characteristics-could serve as a useful proxy for subsurface hydrogeological heterogeneity at the basin scale, which is much more difficult to establish. The geomorphic model of the Sutlej-Yamuna aquifer system could, for example, be used as a broad framework for predicting likely bulk aquifer percentage or the probability of intercepting aquifer bodies of a given thickness in new boreholes, based only on the geomorphic setting of the borehole locality. The geomorphic setting could also be used to guide specific groundwater management approaches-for example, focusing artificial recharge schemes in proximal fan areas that are inferred to have abundant thick subsurface aquifer bodies, and thus a high specific yield [e.g., as applied in the Central Valley Aquifer of California ;Faunt, 2009]. Testing this approach will require more detailed information on channel body dimensions, depositional ages, and the extent of both vertical and lateral connectivity.
Finally, we note that a more refined and integrated depositional framework than hitherto achieved for the Indo-Gangetic plains is now possible with the combined use of satellite imagery and DEM data. When coupled with publicly available CGWB aquifer thickness logs, the aquifer geometry can now be linked to the surface-derived geomorphic framework. Thus, our approach of establishing a geomorphic framework to help understand, and potentially even predict, the subsurface distribution and thickness of aquifer bodies across the entire aquifer system could be applied to other alluvial aquifers in the Indo-Gangetic basin, or elsewhere. The framework could, of course, be refined by comparing predicted aquifer percentages or aquifer body thicknesses to new drilling results in poorly characterized parts of the system. It would also be highly instructive to compare the geomorphic framework to spatial variability in groundwater level change, abstraction, or recharge, to evaluate the large-scale effects of the aquifer body variations that we document here.

Key Unknowns
While the regional coverage of our borehole data is extensive, the results of this study are based nevertheless on relatively widely spaced data on aquifer body thickness. This raises an important issue, because the likely aquifer body widths that we infer on the basis of surface observations (5-10 km) are smaller than the median spacing between adjacent boreholes of ∼7 km. Thus, full characterization of aquifer body dimensions would require independent subsurface evidence of their widths or the ability to resolve individual channel bodies in VAN DIJK ET AL. LINKING FAN MORPHOLOGY TO STRATIGRAPHY 219 the stratigraphy. We are also limited to aquifer thickness data that have been classified from geophysical logs, yielding inferred aquifer body thicknesses that are somewhat different from true lithological units. Finally, we lack age control on the aquifer bodies, which would allow us to understand both the patterns and rates of fan construction and aquifer body deposition, and to correlate between different depositional units in the subsurface. The lack of depositional ages means that we have a very limited understanding of the vertical stacking pattern within the Sutlej and Yamuna fan systems and cannot constrain the avulsion frequency or avulsion magnitudes through time.

Conclusions
We have shown that the distribution of alluvial aquifer bodies in the Sutlej and Yamuna fans of northwest India depends at a broad scale on geomorphic setting and thus on the processes and patterns of deposition in the Himalayan foreland. Analysis of an extensive aquifer thickness data set shows that, across the Sutlej and Yamuna sediment fan systems, individual aquifer bodies have a median thickness of 6-7 m, and they are interpreted to be less than 5-10 km wide because of the lack of clear correlation between adjacent boreholes. The interfan area between the fan apices has both a lower overall percentage of aquifer bodies and thinner aquifer bodies, on average, than the Sutlej and Yamuna fans. The geomorphic setting-specifically, the distinction between fan, interfan, and fan marginal depositional units-thus provides a "framework" that defines clear differences in subsurface aquifer body dimensions and distributions.
The aquifer body thickness distribution remains approximately similar over different depth intervals, which suggests that the paleomorphology and depositional conditions of the sediment routing systems into the foreland have remained consistent over at least the time required to deposit the upper 200 m of stratigraphy. The percentage of aquifer material in individual aquifer thickness logs, however, decreases downstream, although the scaling exponent on the thickness distribution remains the same, indicating that aquifer bodies make up a smaller fraction of the basin fill in the down fan direction but do not thin appreciably. This indicates that rivers on the fan system likely maintained their water and sediment discharge over the lateral dimensions of the Sutlej and Yamuna fans (i.e., up to about 250 km from the mountain front). The aquifer body thickness distributions from the fans and the interfan area are heavy-tailed, and the aquifer body persistence index indicates that aquifer deposits in the fans show evidence for persistent channel positions over depth intervals of about 2-4 times the median aquifer body thickness or roughly the 75th percentile of thickness (that is, up to ∼14 m). Over larger stratigraphic thicknesses, the aquifer thickness logs show evidence of compensational behavior, perhaps related to large-scale avulsion and abandonment of channel corridors. We infer from these observations that the thickest aquifer units are likely to be stacked, multistoried sand bodies that were deposited during persistent reoccupation of particular corridors, possible associated with incised valleys. This inference is important because it implies high vertical connectivity within those stacked sand bodies but disconnection and low lateral (across fan) connectivity due to channel avulsion and abandonment of those corridors.
In conclusion, the geomorphic setting of the aquifer system provides a first-order control on the spatial distribution of aquifer bodies across the study area. The framework that we define here could be used to anticipate bulk aquifer characteristics, including volumetric percentage and likely thickness of aquifer bodies, even in regions without widespread borehole records. This geomorphic framework should be considered in any future approaches to regional-scale aquifer characterization and management. positive reviews by Gary Weissmann, one anonymous reviewer, Associate Editor Jason Kean, and Editor Giovanni Coco helped to clarify and strengthen the manuscript. To obtain the data used in this paper, please contact the authors.
VAN DIJK ET AL.
LINKING FAN MORPHOLOGY TO STRATIGRAPHY 220