Regional analysis of groundwater droughts using hydrograph classification

Groundwater drought is a spatially and temporally variable phenomenon. Here we describe the development of a method to regionally analyse and quantify groundwater drought. The method uses a cluster analysis technique (nonhierarchical k-means) to classify standardised groundwater level hydrographs (the standardised groundwater level index, SGI) prior to analysis of their groundwater drought characteristics, and has been tested using 74 groundwater level time series from Lincolnshire, UK. Using the test data set, six clusters of hydrographs have been identified. For each cluster a correlation can be established between the mean SGI and a mean standardised precipitation index (SPI), where each cluster is associated with a different SPI accumulation period. Based on a comparison of SPI time series for each cluster and for the study area as a whole, it is inferred that the clusters are independent of the driving meteorology and are primarily a function of catchment and hydrogeological factors. This inference is supported by the observation that the majority of sites in each cluster are associated with one of the principal aquifers in the study region. The groundwater drought characteristics of the three largest clusters, which constitute ∼ 80 % of the sites, have been analysed. There are differences in the distributions of drought duration, magnitude and intensity of groundwater drought events between the three clusters as a function of autocorrelation of the mean SGI time series for each cluster. In addition, there are differences between the clusters in their response to three major multi-annual droughts that occurred during the analysis period. For example, sites in the cluster with the longest SGI autocorrelation experience the greatest-magnitude droughts and are the slowest to recover from major droughts, with groundwater drought conditions typically persisting at least 6 months longer than at sites in the other clusters. Membership of the clusters is shown to be related to unsaturated zone thickness at individual boreholes. This last observation emphasises the importance of catchment and aquifer characteristics as (non-trivial) controls on groundwater drought hydrographs. The method of analysis is flexible and can be adapted to a wide range of hydrogeological settings while enabling a consistent approach to the quantification of regional differences in response of groundwater to meteorological drought.


Introduction
Groundwater drought is a type of hydrological drought characterised by sustained low groundwater levels, reduced base flow and reduced flows to springs and groundwater-fed rivers and wetlands (Van Lanen and Peters, 2000;Tallaksen and Van Lanen, 2004;Mishra and Singh, 2010;Van Loon, 2015).Like other hydrological aspects of drought, groundwater droughts are not a simple function of meteorological drivers.The impact of droughts on regional groundwater resources can vary in space and time.This is because the response of groundwater systems to meteorological droughts, through changes in groundwater levels and baseflow to groundwatersupported rivers, is influenced by spatial variations in intrinsic catchment and aquifer characteristics and processes.These include highly non-linear unsaturated zone processes, recharge, and saturated groundwater storage, flow and discharge over a range of spatial and temporal scales (Tallaksen Published by Copernicus Publications on behalf of the European Geosciences Union.et al., 2009;Bloomfield and Marchant, 2013;Van Lanen et al., 2013;Van Loon and Laaha, 2015).
In order to improve the design and operation of groundwater drought monitoring networks, the analysis and interpretation of data from such networks, and, more generally, water resource management at the onset, during and after episodes of groundwater drought, there is a need for a much better understanding of the heterogeneous spatio-temporal response of aquifers to major meteorological droughts (Bloomfield and Marchant, 2013).This includes the need for robust methods to systematically characterise and quantify the heterogeneous response of groundwater to meteorological droughts at a regional scale prior to investigation and attribution of the causes of any heterogeneous response.Despite extensive work on the regional analysis of meteorological and other hydrological droughts, to date there has been no systematic investigation of heterogeneities in groundwater droughts at the regional scale.This paper describes the application of one such suite of methods to regionally analyse groundwater level hydrographs and to assess variations in the spatial response of groundwater to meteorological droughts using a case study from the UK.

Controls on spatial heterogeneity in groundwater drought
A few previous studies have presented evidence for the spatially heterogeneous response of groundwater to meteorological droughts.To help develop an optimal monitoring network for groundwater resources under drought conditions, Chang and Teoh (1995) described the heterogeneous response of groundwater levels at 13 observation boreholes to meteorological droughts across a basin in Ohio, USA, although they did not investigate the hydrogeological causes of the heterogeneity.Van Lanen (2005) and Van Lanen and Tallaksen (2007) observed that drought characteristics derived from groundwater levels have "spatial effects" and noted that these spatial effects on groundwater drought are an important consideration when monitoring droughts using groundwater levels.Van Lanen and Tallaksen (2007) compared modelled groundwater recharge and discharge for a humid continental climate (Missouri, USA) and a tropical savannah climate (Guinea) for quick-and slow-responding catchments and showed that both climatology and the responsiveness of the catchment as defined by the aquifer characteristics have an influence on drought generation.Peters et al. (2006) investigated the propagation and spatial distribution of aspects of modelled groundwater drought, including recharge, groundwater level and groundwater discharge in the Pang catchment in the UK.They found that short droughts in groundwater levels were most severe near streams and were attenuated with distance from the streams, that longer periods of belowaverage recharge had more effect on suppressing groundwater levels on interfluves near groundwater divides, and that droughts in groundwater discharge are more attenuated up-stream and less so downstream in the catchment.Tallaksen et al. (2009) also modelled the spatio-temporal response of the Pang catchment to drought events and found large differences between the spatio-temporal response of groundwater recharge, level and discharge and the driving meteorological droughts, where droughts in groundwater recharge and levels were found to cover relatively small areas but last longer than the meteorological droughts.Mendicino et al. (2008) developed a groundwater resource index for drought monitoring and forecasting based on a simple distributed run-off/water balance model, and they evaluated the use of the index in three catchments in southern Italy.They found that the groundwater resource index was highly spatially variable and related it to variations in hydraulic conductivity across the catchments.Using a newly developed groundwater drought index, the standardised groundwater level index (SGI), Bloomfield and Marchant (2013) also investigated hydrogeological controls on groundwater drought.Based on 14 observation boreholes in different catchments across England, UK, they showed that groundwater drought duration depended on the autocorrelation structure of SGI time series.This was in turn inferred to be a function both of spatially varying recharge processes and of saturated flow processes within the local aquifer systems.
Although not previously applied to groundwater drought, CA and/or PCA techniques have been used to classify groundwater level hydrographs for a range of purposes.Winter et al. (2000) classified groundwater hydrographs from three small lake-dominated catchments to investigate groundwater recharge and differences in the hydrographs as a function of the geology of the catchments.Similarly, Moon et al. (2004) applied PCA to 66 groundwater level hydrographs from South Korea to characterise the spatial variability in groundwater recharge.Upton and Jackson (2011) used CA and PCA (following a methodology developed by Hannah et al., 2000) with 52 groundwater level hydrographs from the Pang and Lambourn catchments in the UK to produce regional or "master" hydrographs for modelling the spatial distribution of groundwater flooding.Here we present the first systematic regional analysis of groundwater droughts using a case study from Lincolnshire, UK.The case study consists of 74 groundwater hydrographs from an area of approximately 8000 km 2 that includes three regionally important aquifers, the Lincolnshire Limestone, the Chalk and the Spilsby Sandstone aquifers, each with contrasting aquifer characteristics (Sect.2).The groundwater hydrographs have been normalised using the SGI technique of Bloomfield and Marchant (2013), and groups or clusters of similar groundwater hydrographs have been identified using CA, where hydrogeologically meaningful clusters are identified by explicitly searching for groups of hydrographs that can be explained by a posteriori knowledge of the groundwater system (Sect.4.2).The drought characteristics of the clusters have been quantified in terms of drought event duration, magnitude and intensity, and the impact of three major, multi-annual droughts on the SGI time series has been investigated (Sect. 4.4).Controls on the groundwater drought response in each of the clusters have been explored and the results briefly discussed in terms of the implications for monitoring and managing groundwater droughts (Sect.5).

The case study
The case study area of Lincolnshire is situated in the east of England, UK.It is bounded by the North Sea to the east, the Wash estuary to the south and the Humber estuary to the north (Fig. 1).The area is predominantly rural with highly productive agricultural and horticultural land, fens and estuarine wetlands.Lincoln, Boston and Scunthorpe are the principal small conurbations in the study area.The land is generally flat and low-lying, typically less than 30 m a.s.l.(above sea level), apart from the Chalk of the Lincolnshire Wolds and the Lincolnshire Limestone outcrop, which form northwest-southeast-trending escarpments that reach elevations of approximately 150 and 70 m a.s.l.respectively.

Hydrometeorology and drought history
As a first-order approximation, it is assumed that the broad meteorological drought history of the study area is spatially homogeneous.This assumption means that any relative differences in drought histories between sites or clusters need to be explained in terms of catchment or hydrogeological factors, rather than differences in the drought climatology.This assumption is tested as part of the analysis of correlations between precipitation and regional groundwater levels (see Sect. 4.2).It is also supported by the observations that the whole study area is governed by the same broad climatic patterns, i.e. rain-bearing low-pressure systems from the Atlantic and high-pressure systems leading to a lack of rainfall, with only small variation in annual precipitation across the region (Marsh and Hannaford, 2008).The assumption is also consistent with the previously documented spatial coherence of major hydrological (surface water) droughts in the UK (Hannaford et al., 2011;Fleig et al., 2011;Folland et al., 2015) where the current study area falls within a homogeneous drought region ("region 4" of Hannaford et al. (2011), "region GB4" of Fleig et al. (2012) and Kingston et al. (2013), and the "English Lowlands" of Folland et al. (2015)), although it is noted that the effects of landscape processes can cause heterogeneous meteorological signals to become attenuated (Van Loon, 2015).
Mean annual rainfall varies across the study area from about 600 to 700 mm (Marsh and Hannaford, 2008).The groundwater hydrographs used in the study have been analysed from 1983 to 2012.During this period, three multiannual episodes of drought have previously been documented by Marsh et al. (2007Marsh et al. ( , 2013)), Kendon et al. (2013), Parry and Marsh (2013) and Folland et al. (2015) as follows: 1988 to 1992, 1995 to 1997 and 2010 to 2012 respectively.All are known to have been major drought events causing reduced surface flows and suppressed groundwater levels throughout large areas of central, eastern and southern UK as well as over parts of northwestern Europe (Lloyd-Hughes and Saunders, 2002;Lloyd-Hughes et al., 2010;Hannaford et al., 2011;Fleig et al., 2012;Kingston et al., 2013).

Geology and hydrogeology
The study area consists of a sequence of Jurassic and Cretaceous aquifers separated by low-permeability clay and shale units.The whole sequence generally dips gently eastwards, and where each of the aquifer units passes under an overlying low-permeability formation they typically become confined.The whole sequence is unconformably overlain by Quaternary superficial deposits.Figure 1 shows the distribution of the three main aquifers in the region -the Jurassic Lincolnshire Limestone; the Lower Cretaceous-Upper Jurassic Spilsby Sandstone; and the Upper Cretaceous Chalkand includes a schematic cross section of the hydrostratigraphy of the study area.These aquifers are hydrogeologically distinct from each other, and two of them, the Lincolnshire Limestone and the Chalk, have previously documented spatial variability.Below we summarise these features as they inform the heuristic rules used in Sect.3.2.2 to guide the selection of clusters as part of the CA.
The Lincolnshire Limestone Formation is an oolitic limestone with fine-grained, micritic and peloidal units (Allen et al., 1997), and it is up to 40 m thick at outcrop in the west.It dips and thins to the east, where it becomes confined and eventually pinches out down-dip.Maximum unsaturated zone thickness is up to about 45 m towards the southwest of the outcrop.Groundwater movement is almost entirely by fracture flow along well-developed bedding plane fractures and joints.Abstraction takes place mainly from the region immediately to the east of the outcrop.It has highly variable transmissivities and storage coefficients typical of a fractured limestone.Allen et al. (1997) have reported a wide range of transmissivity values for the Lincolnshire Limestone with an interquartile range of 260 to 2260 m 2 day −1 and a geometric mean of 660 m 2 day −1 , with slightly higher transmissivities being reported from the south of the region, and a very wide range of storage coefficients from 2 × 10 −7 to 0.58.
The Spilsby Sandstone aquifer is up to about 30 m thick, consisting of a variably, but often poorly cemented, pebbly quartz sandstone with alternating thin clays and marls (Whitehead and Lawrence, 2006).It outcrops along the foot of the Wolds escarpment (Fig. 1), where it is associated with springs and maximum unsaturated zone thickness is about 30m.It dips to the east and away from outcrop and is generally confined by clays above and below (Fig. 1).Jones et al. (2000) reported transmissivity values in the range 130 to 170 m 2 day −1 and a geometric mean of 140 m 2 day −1 , with storage coefficients ranging from 1 × 10 −4 to 1 × 10 −3 and with a geometric mean of 4 × 10 −4 .
The Chalk is a microporous fractured limestone (Bloomfield et al., 1995).Storage and transmissivity are controlled by local sub-karstic development of the fracture network (Bloomfield, 1996;Maurice et al., 2006).The Chalk group reaches a thickness of over 250 m.Groundwater flows from the recharge areas in the west, eastward down-dip towards and into the confined Chalk to the east.The Chalk bedrock surface was significantly altered during the Ipswichian interglacial of the Quaternary.As a result of glacial activity a cliff line and wave-cut platform were eroded into the Chalk (Fig. 1).The Chalk to the east of the palaeo-cliff line is now buried beneath a covering of till, sand and gravel superficial deposits (Whitehead and Lawrence, 2006).Maximum unsaturated zone thickness occurs towards the northwest of the Chalk outcrop and is about 60 m, contrasting with the relatively thin unsaturated zone to the east of the palaeo-cliff line.Allen et al. (1997) and Whitehead and Lawrence (2006) have reported that transmissivity values differ between the northern and southern Chalk in Lincolnshire.In the northern part of the region, transmissivity has an interquartile range of 1020 to 6070 m 2 day −1 with a geomet-ric mean of 2350 m 2 day −1 , whereas in the southern area, in the region of the eroded Chalk, transmissivity is slightly reduced and has an interquartile range of 850 to 3010 m 2 day −1 with a geometric mean of 1380 m 2 day −1 .Similarly, Allen et al. (1997) report storage coefficients with an interquartile range of 3.5 × 10 −5 to 1.5 × 10 −3 and with a geometric mean of 2 × 10 −4 for the northern Chalk and 6.1 × 10 −5 to 2.7 × 10 −3 and with a geometric mean of 1.5 × 10 −3 for the southern Chalk.
The Quaternary superficial deposits in the study area comprise glaciofluvial sand and gravels and tills; peat; tidal flat deposits; river terrace sands and gravels; and overlying alluvium.The Lincolnshire Limestone Formation and the western part of the Chalk outcrop are largely absent of superficial cover.

Data
Groundwater level data for the 74 observation boreholes (Fig. 1) has been provided by the Environment Agency from their groundwater level monitoring network database (Environment Agency, 2014).Prior to the study none of the sites were believed to be significantly impacted by abstraction, although all three regional aquifers are used for public water supply, abstractions for agricultural irrigation and industrial use (Allen et al., 1997;Whitehead and Lawrence, 2006).Where observation boreholes penetrate both the Chalk and underlying Spilsby Sandstone aquifer, the boreholes are completed with screens so that they monitor water levels in only one of the two aquifers.Groundwater levels have been recorded over a range of frequencies, but typically at weekly to monthly time steps.Based on the raw groundwater level data, mean monthly groundwater levels have been estimated.If no observations were available for a given month, then a linear interpolation was used to estimate the monthly groundwater levels following the method described by Bloomfield and Marchant (2013).
Precipitation data have been taken from the Centre for Ecology and Hydrology's Continuous Estimation of River Flows (CERF) 1 km gridded precipitation data set (Keller et al., 2005;Dore et al., 2012;Bloomfield and Marchant, 2013).CERF daily gridded precipitation data are generated from rain gauge data held in the UK Met Office national precipitation monitoring network.A triangular planes methodology is used to produce a daily 1 km 2 grid based on a weighted average (inverse distance) of the three nearest rain gauges.Daily rainfall is then summed to give total monthly gridded rainfall.The precipitation data that are used with each groundwater level observation site are the monthly total for the CERF 1 km 2 grid square that contains the given groundwater observation borehole.

Hydrograph normalisation using the SGI method
The groundwater level hydrographs have been normalised to the SGI of Bloomfield and Marchant (2013).This is a non-parametric normalisation of data that assigns a value to the monthly groundwater levels based on their rank within groundwater levels for a given month from a given hydrograph.The normal scores transform is undertaken by applying the inverse normal cumulative distribution function to n equally spaced p i values ranging from 1/(2n) to 1 − 1/(2n).The values that result are the SGI values.They are then reordered such that the largest SGI value is assigned to the i for which p i is largest, the second-largest SGI value is assigned to the i for which p i is second largest and so on.In summary, for each of the 74 study sites, normalised indices are estimated from the groundwater level data for each calendar month using the normal scores transform.These normalised indices are then merged to form a continuous SGI.Precipitation records for each site have also been normalised.At each site a version of the standardised precipitation index (SPI) after McKee et al. (1993) has been estimated for precipitation accumulation periods of 1, 2, . . ., 36 months.For consistency between groundwater and precipitation indices, SPIs are estimated using the normal scores transform applied to accumulated precipitation data for each calendar month.

Cluster analysis
Cluster analysis attempts to identify clusters of similar individuals amongst a multivariate data set.In the context of this paper CA is used to form clusters of groundwater level hydrographs which exhibit similar fluctuations in their SGI time series.A wide range of CA algorithms exist.They are most coarsely distinguished according to whether or not they assume that the resultant clusters are hierarchical.Given the wide variety of algorithms, it is difficult to decide upon the best approach to cluster a particular data set.Webster and Oliver (1990) stress that this decision is rather subjective, although previous studies that have used CA to cluster hydrographs have typically justified their choice of algorithm by claiming that some produce more physically interpretable groupings.For example, Hannah et al. (2000) used the agglomerative hierarchical average linkage algorithm as they thought it was more interpretable than alternatives such as the centroid and Ward's clustering procedures.Webster and Oliver (1990) recommend that multiple clustering algorithms should be applied and expert knowledge of the system being investigated used to decide which set of clusters is most relevant.In this paper we adapt this approach by applying one hierarchical and one non-hierarchical method.
Hierarchical classifiers require a measure of the similarity (or dissimilarity) between each pair of individuals.Common examples include the Euclidean distance or the correla-tion between the measurements of the individuals.The pairwise similarities between s individuals are expressed in a s × s matrix B. A mathematical criterion is then used to allocate the individuals to different clusters in a manner that maximises the similarity between the individuals within the groups whilst minimising the similarity between individuals in different clusters.For our hierarchical clusters we measure the similarity between groundwater level hydrographs by the correlation matrix of their SGI time series and then apply the agglomerative hierarchical complete-linkage strategy (Webster and Oliver, 1990) to merge the boreholes into clusters.
We also apply the commonly used non-hierarchical kmeans clustering algorithm.It is widely used in spatial analysis studies; for example, Santos et al. (2010), Raziei et al. (2012) and Sadri et al. (2014) have all used the k-means clustering algorithm to investigate the regional characteristics of droughts.The approach partitions the individuals into a specified number of clusters.A numerical optimisation routine is used to select the partitioning which maximises the similarity between each individual and the centroid of the cluster in which it is contained.Again there is flexibility in the choice of similarity measure and the manner in which the centroid of a cluster is calculated.We use the squared Euclidean distance between the vectors of time series observations from each site to assess similarity and define the centroid of a cluster as the multi-dimensional mean of the time series within the cluster.
Clustering methods do not produce a unique partitioning of a given data set on their own, and for both the hierarchical and non-hierarchical approaches there remains the issue of deciding upon the optimal number of clusters.This can be achieved by asking an expert on the system in question to compare the attributes of clusterings consisting of a different number of groups.Here we use a rule-based approach to help identify the number of clusters based on knowledge of the general hydrogeology of the study area.Bloomfield and Marchant (2013) have previously shown that groundwater drought characteristics are a function of unsaturated zone thickness in fractured aquifers such as the Lincolnshire Limestone and Chalk aquifers, and that when a broader range of aquifer types are considered groundwater drought characteristics are also a function of the hydraulic diffusivity of aquifers.Here we use these observations and knowledge of the spatial variation in these features across the three aquifers in the study area (Sect.2.2) to design rules to aid in the selection of clusters.The rules adopted for the current study are to identify the smallest number of clusters that (i) broadly resolve the spatial distribution of the three aquifers across the study region; (ii) distinguish more than one region of the Lincolnshire Limestone, given the previously documented N-S variation in aquifer properties and unsaturated zone thickness across the Lincolnshire Limestone aquifer (Allen et al., 1997); and (iii) distinguish more than one region of the Chalk, given variations in aquifer properties and unsaturated zone thickness across the Chalk aquifer both N-S and across the buried cliff line (Allen et al., 1997).Note that this set of rules is specific to the current study; however, for any given study area the target number of classes and hence the rules used can be adapted to reflect the regional hydrogeology and in particular any knowledge of heterogeneity in the aquifer systems under investigation.However, mathematical criteria can also be used as a guide to clustering.We also calculate the RMSSD, the square root of sum of the squared Euclidean distance between each individual and the centroid of the group to which it is allocated.In combination with expert judgement related to the system under consideration, it is common practice to inform the choice of the number of clusters using plots of RMSSD versus cluster number.Since RMSSD decreases non-linearly as the number of clusters increases, a cluster number is selected associated with a decrease in the rate of RMSSD decline.

Autocorrelation structure of the SGI time series
Bloomfield and Marchant (2013) demonstrated the importance of the autocorrelation structure of SGI time series for groundwater drought studies by establishing a relationship between the range of significant autocorrelation in the SGI series, m max , and corresponding SPI.They showed that m max scales linearly with q max , where q max is the SPI accumulation period which leads to the strongest correlation between SGI and SPI.Both m max and q max are also used here to characterise and quantify groundwater droughts within each of the clusters of groundwater hydrographs and have been estimated as follows.
If the mean SGI for a borehole is denoted by SGI, then the kth sample autocovariance coefficient is defined to be and the kth sample autocorrelation coefficient is where g 0 reduces to the population variance function (see Eq. ( 1) when k = 0).The correlogram is a plot of r k against k.If there is no correlation between the SGI(i) observed k months apart and if the SGI values are normally distributed, then r k is approximately normally distributed, with mean 0 and variance 1/n.Therefore values of r k with magnitude greater than 2/ √ n indicate significant correlation at approximately the 5 % level.We define the range of significant temporal correlation of a SGI time series to be the largest m, m max , for which r k > 2/ √ n for all k ≤ m.Since all of our groundwater records are of n = 355 months, the threshold on r k is equal to 0.11.To estimate q max , Pearson correlation coefficients are calculated between SGI and SPI with accumulation periods of q = 1, 2, . . ., 36 months, and the accumulation period associated with the maximum correlation gives q max .

Identification of regional droughts from average SPI and SGI time series
Before undertaking the regional drought analysis, the correlation between mean SPI and SGI for the entire region, based on all 74 sites, was investigated and the large-scale drought history of the study area were defined.Figure 2a is a heat map showing the correlation coefficient between SPI for precipitation accumulation periods q = 1 to 36 months and SGI for lags between SPI and SGI of 0 to 5 months based on average values of SPI and SGI for all 74 sites.Dark blue denotes zero correlation and dark red a perfect correlation.Figure 2a shows that there is a good correlation between SPI and SGI.The strongest correlation (0.84, denoted by the closed black circle in Fig. 2a) is for a precipitation accumulation period (q max ) of 12 months (SPI 12 ) with no lag between the SGI and SPI time series.This is consistent with the observations of Bloomfield and Marchant (2013), who previously reported q max for a variety of groundwater hydrographs from the UK with an average of 13 months, and Folland et al. (2015), who reported a q max of 12 months for aggregated time series representing the English Lowlands.Figure 2b and c, the average SPI 12 and SGI time series respectively, have similar features.For example, episodes of high groundwater levels in 1983, 1994, 2002, and 2008 correspond with high values of SPI 12 .Three episodes of regionally significant groundwater drought associated with prolonged low groundwater levels from October 1988 to November 1993, May 1995 to February 1998, and from August 2010 to August 2012 correspond closely with episodes of meteorological drought in the SPI 12 time series and are consistent with those identified by previous studies (Lloyd-Hughes and Saunders, 2002;Marsh et al., 2007Marsh et al., , 2013;;Kendon et al., 2013;Hannaford et al., 2011;Parry and Marsh, 2013;Folland et al., 2015).It is inferred from these observations that the large-scale drought history of the study area is represented well by the average SPI 12 and SGI time series.

Regional analysis of the SGI hydrographs
CA has been used to analyse the heterogeneous response of groundwater to droughts across the study region.Clustering has been undertaken using both an agglomerative hierarchical complete-linkage algorithm and a non-hierarchical k-means clustering algorithm, and the resulting clusters searched for those that are hydrogeologically meaningful and that can be explained by known features of the catchment and groundwater systems.Figure 3a is a dendrogram that fully illustrates the level of similarity between individuals within the clusters formed by the hierarchical clustering.The number of clusters is controlled through the threshold on the distance between groups.For example, a threshold of 0.62 leads Figure 3b and c show that the spatial distribution of sites as a function of the clusters formed by the hierarchical and nonhierarchical approaches are broadly similar, so the choice of clustering algorithm is based on a plot of RMSSD against number of clusters.Figure 4 shows that the RMSSD for the k-means clustering is systematically lower than that for the hierarchical clustering algorithm where there are three clusters or more, so we have chosen to use the non-hierarchical k-means clustering approach.Note also that both clustering algorithms are better than a clustering scheme based solely on the three classes of aquifer (e.g.Lincolnshire Limestone, Chalk and Spilsby Sandstone).However, an optimal number of k-mean clusters is not clearly evident in Fig. 4.After careful inspection of the clusters formed by a range of k-means clustering classes and a consideration of the study-specific clustering rules described in Sect.3.2.2,k = 6 was selected.Based on k-means clustering where k = 6, Fig. 3c shows the distribution of sites between the six clusters (cluster 1 to cluster 6, or CL1-CL6).
It can be seen from Fig. 3c that the resulting k-means clusters have a degree of spatial coherency.We have previously assumed that such spatial correlations in the SGI time series are primarily a function of catchment and hydrogeological factors and not a consequence of heterogeneity in the driving meteorology.Here we test if this is the case, prior to further exploration of the features of each cluster, by investigating if precipitation associated with each cluster is substantially different from regional average precipitation.To do this, we first need to identify a representative accumulation period, q max , for precipitation for each cluster.
Figure 5 is a set of heat maps, similar to Fig. 2a, showing the correlation between SPI for precipitation accumula-  tion periods, q, 1 to 36 months, and SGI for lags between SPI and SGI time series of 0 to 5 months for each of the six clusters.Dark blue denotes zero correlation, and dark red a perfect correlation, with the strongest correlation for each cluster marked by the closed black circle.Table 1 gives q max for each cluster and also gives the maximum associated correlation coefficient.In all cases except CL2, the maximum correlation between SPI and SGI is found where there is no lag between the two time series.For CL2 it is found at a lag of 1 month.The highest correlations are for CL2, CL4 and CL1 at 0.86, 0.82 and 0.74 respectively.The correlations for CL3 and CL5 are moderate (0.36 and 0.53), and for CL6 there is effectively no correlation (0.09).This is consistent with the observations made in Sect.4.3 below that linear trends in CL3 and CL5 appear to affect the SGI time series and that the SGI hydrograph for CL6 appears to be anomalous, departing from the mean regional SGI and SPI signals.Values of q max for CL1 to CL5 from Fig. 5 are 4, 16, 15, 9, and 17 months respectively.Based on these, Fig. 6 shows SPI time series for each cluster, where black lines are the mean SPI for the cluster and the red lines are average SPI across the study area based on the same cluster-specific q max .Since Fig. 6  we infer that heterogeneity in the driving meteorology across the study region, or at least between the clusters as defined here, does not play an important role in the clustering process and that membership of clusters is dominated by catchment or hydrogeological factors.

Characteristic features of the SGI hydrograph clusters
Figure 7 shows the mean SGI time series for each cluster.Two main qualitative observations can be made regarding the SGI hydrographs.Five of the six clusters have a similar overall form to the mean SGI hydrograph for the region (Fig. 2c) showing common patterns of low (and high) groundwater level stand.However, CL6 appears to be an exception with a different overall form to the SGI hydrograph -it also exhibits an anomalous step change in SGI from drought to high groundwater level stand over an 8-month period from May 1990 to December 1990.Secondly, two of the clusters, CL3 and CL5, appear to show declining linear trends in SGI, making direct comparison of drought histories between these and other clusters problematic.Bloomfield and Marchant ( 2013) have previously shown that m max , a measure of the significant autocorrelation length of SGI time series, relates to features of groundwater drought.A similar analysis of autocorrelation structure of SGI time series for each cluster is presented here.Figure 8 shows autocorrelation plots for SGI hydrographs for each of the six clusters.In each figure the pale grey lines are autocorrelation plots for individual sites and the solid black line is the autocorrelation plot for the mean SGI time series for the cluster, with the horizontal dashed line indicating the significant level of autocorrelation based on the record length.Based on these plots, values of m max for the mean SGI time series for each cluster are given in Table 1.Values of m max for CL3, CL5 and CL6 are anomalously large, consistent with the anomalous features of these SGI hydrographs described above.For the remaining clusters, Fig. 8 and Table 1 show that CL1 has the shortest autocorrelation of 15 months.In comparison, CL2 has an autocorrelation of 23 months and CL4 is intermediate at 18 months.
These contrasting characteristics between the clusters can be seen clearly in Fig. 9a, which illustrates SGI time series for all sites within each cluster, grouped in their respective clusters and presented in the form of a heat map where low values of SGI (associated with drought conditions) are in shades of green to red (increasing drought intensity) and episodes of high groundwater level stand are in shades of green to blue (increasingly high groundwater levels).The  -CL2 is dominated by sites from the northern part of the Chalk.The cluster has the longest mean SGI autocorrelation (m max of 23 months), and hydrographs within CL2 are highly correlated, indicating a high degree of coherency in groundwater levels across the northern part of the Chalk in the study area.
-CL3 is a relatively small cluster of six sites, four of which are from the confined Spilsby Sandstone and two from the Lincolnshire Limestone.The main feature of the cluster is a trend in decreasing SGI across the observational record.This trend is consistent with a previous water balance assessment for the Spilsby Sandstone (Whitehead and Lawrence, 2006), where annual groundwater deficits have been reported.The sites in this cluster are inferred to be possibly variably impacted by long-term abstraction.Given this inference and the small size of the cluster of sites, CL3 is not included in the subsequent analysis of groundwater droughts.
-CL4 is dominated by sites from the southern Lincolnshire Limestone and also includes five unconfined sites on the southern Chalk and one site located in the northern Lincolnshire Limestone.It has a moderate autocorrelation, m max , of 18 months.Individual SGI hydrographs within the cluster show a moderate degree of coherency.
-CL5 is a small cluster of five sites all from the southeastern Chalk to the east of the palaeo-wave-cut platform, and they are the five sites closest to the coast.It has a moderately long autocorrelation, m max , of 28 months that may be affected by an apparent weak trend in declining SGI -there is only a weak correlation between SPI and SGI.Given the small size of the cluster and the apparent trend in mean SGI, CL5 is not included in the subsequent analysis of groundwater droughts.
-CL6 consists of three SGI hydrographs from the confined Spilsby Sandstone aquifer.The hydrographs are characterised by an anomalous step change in SGI from drought to high groundwater level stand over an 8month period from May 1990 to December 1990.The mean SGI hydrograph shows no correlation with the other five clusters, and there is no correlation between SPI and SGI within the cluster.All three sites are within a radius of about 3 km of a public water supply borehole, and it is inferred that groundwater levels may be influenced by abstraction.So, as with CL3 and CL5, this very small cluster is not included in the subsequent analysis of groundwater droughts.

Analysis of droughts using the hydrographs from CL1, 2 and 4
Clusters CL1, CL2 and CL4 consist of 61 of the 74 hydrographs analysed.Here the characteristics of groundwater droughts in these clusters are quantified, and the response of the clusters to three major drought episodes is investigated.The duration, magnitude and mean intensity of groundwater drought events have been investigated based on an analysis of the SGI hydrographs where, following the convention of McKee et al. (1993), negative values of SGI denote drought conditions (note, however, that the current convention of the World Meteorological Organization for SPI refers to drought conditions where SPI is continuously negative and reaches and intensity of −1.0 or less and that negative values between 0 and −1 are classified as near normal and simply indicate less than a median precipitation; World Meteorological Organization, 2012).Groundwater drought duration, D, is taken to be the total number of consecutive months where SGI is negative.Groundwater drought magnitude, M, is taken to be the total cumulative value of monthly SGI for a given drought event, and mean drought intensity, I, is given by M / D. Summary drought statistics for CL1, CL2 and CL4 are given in Table 2. Table 2 shows that there are differences in the character of the groundwater drought events in the SGI hydrographs for clusters CL1, CL2 and CL3.For example, CL1 has more than twice the number of drought episodes (39 episodes) as CL2 (15 episodes), and the average and maximum duration of droughts in CL1 (4.6 and 27 months respectively) are less than half those of CL2 (11.3 and 61 months).The mean drought event magnitude in CL1 (−2.9) is less than half that in CL2 (−7.9), and the mean drought event intensity in CL1 (−0.43) is almost twice that of CL2 (−0.28).In all cases, the drought event statistics for CL4 fall between those for CL1 and CL2.In summary, CL1 exhibits shorter but generally more intense drought episodes compared with CL2, with CL4 drought events being of intermediate character.These relative drought phenomena are a consequence of the degree of autocorrelation in the respective SGI time series, where CL1 has a relatively short autocorrelation compared with relatively long autocorrelation for CL2.This observation is consistent with previous site-specific and modelling studies that noted a similar relationship between the "flashiness" or responsiveness of the groundwater system to meteorological divers and the number of droughts, where quickly responding groundwater systems typically experience more droughts than more slowly responding catchments (Peters et al., 2003;Van Loon and Van Lanen, 2012;Van Lanen et al. 2013).
There is a strong relationship between drought duration and magnitude for all three clusters (Fig. 10), where longer episodes of groundwater drought are associated with droughts of greater magnitude.However, there is no such regular or simple relationship between drought duration and intensity.Maximum drought intensity is similar for all three clusters -for CL1, CL2 and CL4 it is −1.10, −1.05 and −1.13 respectively (Table 2 and Fig. 11) -and is associated with two of the major drought events, i.e. with the latter part of the 1988-1993 drought for CL2 and the 2010-2012 drought for CL1 and CL4. Figure 11 shows the empirical distribution of D, M and I for clusters CL1, CL2 and CL4.Drought duration (Fig. 11) in all three clusters is highly positively skewed with many short drought events and relatively few long drought events.As previously noted, the longest duration droughts are associated with CL2, the cluster with the longest autocorrelation in the SGI time series.These observations are consistent with those of Hisdal and Tallaksen (2003), Tallaksen et al. (2009) and Fleig et al. (2011), who have also described strongly skewed distributions of hydrological drought durations.
Three major, multi-annual droughts have already been described from the regional (Fig. 2) and the cluster-specific (Figs. 7 and 9a) SGI time series.Table 3 summarises differences in the relationships between the driving meteorology and the drought characteristics of each cluster for the three major droughts.Each of the major drought episodes has been quantified using drought characteristics as applied to SPI 12 and SGI for each of the clusters.
The 1988-1993 event was the longest of the three major droughts and consequently had the greatest drought magnitude.The groundwater and meteorological droughts started approximately contemporaneously in the winter of 1988.In CL2 the drought was continuous with negative SGI from November 1988 to November 1993, whereas in CL4 there were two short breaks in the drought and numerous breaks in the drought in CL1.In CL2 there was a gradual intensification in the drought magnitude across the event, peaking in June 1992 at an SGI of −1.85 (4 months after the peak SPI 12 meteorological drought).In contrast, not only were there short breaks in the drought in CL1 and CL4 but there were approximately annual cycles of drought intensification and decline over the 4-year period -these were particularly pronounced in CL4.This is seen in Fig. 9a, where between 1988 and 1993 the drought status of CL4 is designated by the red tones in the heat map, but these tones show a series of approximately annual variations giving the appearance of vertical stripes during that period and within that cluster.However, the most pronounced differences in response to major droughts between clusters CL1, CL2 and CL4 is in the timing of the end of drought.Groundwater drought conditions ended in CL1 and CL4 in May 1993, 7 months after the end of the meteorological drought, but this was still 6 months before the groundwater drought ended in CL2 (Fig. 9a).
The 1995-1997 drought, although shorter than the 1988-1993 drought, followed a similar pattern, with groundwater drought starting approximately contemporaneously with the meteorological drought.Although it was a continuous event for all three clusters (there were no breaks in the drought for CL1 and CL4), CL1 and CL4 again show approximately annual intensifications and declines in drought status during the episode.Such approximately annual changes in drought status are not seen in CL2.The 1995-1997 drought had the greatest magnitude in CL2 due to the prolonged end to the drought in this cluster, with groundwater drought in CL1 and CL4 finishing approximately contemporaneously with the meteorological drought but 6 months later in CL2.The 2011-2012 drought was much shorter than the other two multi-annual droughts, lasting just over a year starting relatively abruptly in early 2012 and finished abruptly in CL1 and CL4 in May 2012 in response to an unusual episode of spring recharge (Parry et al., 2013).The groundwater drought in CL2 again finished relatively late, this time about 3 months later, in August 2012.The relatively short delay in the breaking of the groundwater drought in CL2 compared with CL1 and CL4 probably reflects the relatively smaller groundwater drought deficit accumulated due to the shorter duration and lower magnitude of the drought compared with the 1988-1993 and 1995-1998 drought episodes.

Discussion
The results of the regional analysis of droughts based on cluster analysis are consistent with current conceptualisations of the dynamics of drought in hydrological systems.Propagation of drought through catchments and in particular through the groundwater compartment is well documented (Peters et al., 2003(Peters et al., , 2006;;Tallaksen et al., 2006), and four components of drought propagation are recognised, i.e. pooling, attenuation, lag and lengthening, three of which (attenuation, lag and lengthening) are associated with modifications of drought signals in groundwater (Van Loon, 2015).Attenuation results in smoothing of the maximum drought anomaly, lag describes the delay in the onset of the drought signal as it passes through the hydrological cycle (for example, see Figs. 3a and 4 of Van Loon, 2015), and lengthening extends the period of drought.Considering Table 3, which summarises the three multi-annual droughts, and comparing event magnitude for SPI 12 , CL1, CL2 and CL4, there is, as would be expected, evidence of a general attenuation of the SPI drought signal in the three clusters compared with SPI 12 .Lagging of the multi-annual groundwater droughts behind meteorological droughts is not so easy to quantify unambiguously.Clearly the nature and degree of the lag is sensitive to the rainfall accumulation period used to define the meteorological drought index most closely correlated with SGI.In the present case, accumulation periods of 4, 16, and 9 months are required for CL1, 2 and 4 respectively to achieve optimal correlation between the SPI and SGI time series.Finally, the results of the present study strongly support the concept of lengthening of groundwater drought relative to meteorological drought (Van Loon, 2015).The results demonstrate that lengthening is most pronounced following longer and deeper groundwater droughts.They serve to emphasise that there can be significant differences in the lengthening response between different clusters, even within with the same aquifer.It also appears that the degree of lengthening may also be related to SGI autocorrelation (the greatest degree of lengthening is observed in cluster CL2 associated with the largest SGI autocorrelation, m max ).
The results of the regional analysis add to our current understanding of the controls on groundwater droughts.Bloomfield and Marchant (2013) investigated how unsaturated zone thickness and the hydraulic diffusivity of aquifers may relate to m max .Using 14 SGI time series from four different aquifers around the UK (including one site from the Lincolnshire Limestone and nine sites on the Chalk, although none from the present study), they found that m max was broadly an inverse function of log hydraulic diffusivity, log D diff (where D diff is given by T /S and where T is aquifer transmissivity and S is specific storage of the aquifer).But they also noted that when fractured aquifers -such as the Lincolnshire Limestone and the Chalk, which have similarly high hydraulic diffusivities -are specifically considered there is no clear relationship between m max and log D diff .However, they did find a positive relationship between unsaturated zone thickness and m max for fractured such as the Chalk and Lincolnshire Limestone.Based on this observation, they proposed that unsaturated zone drainage and recharge processes were an important contributory factor in determining autocorrelation or "memory" in groundwater level hydrographs and by inference an influential factor on groundwater drought characteristics, particularly in fracture aquifer systems.Here we investigate if a similar relationship between m max and unsaturated zone thickness holds for CL1, CL2 and CL4, clusters dominated by fractured aquifers.
Figure 12 shows box plots of unsaturated zone thickness for CL1, CL2 and CL4 as a function of m max for each cluster (where unsaturated zone thickness is taken as the mean depth to groundwater recorded for sites in each cluster over the study period).In addition, corresponding observations for 10 boreholes in fractured aquifers from Bloomfield and Marchant (2013) are also shown for reference.The results of the present study are consistent with those of Bloomfield and Marchant (2013, Fig. 13a) and show increasing mean unsaturated zone thickness with increasing cluster m max ; increasing variability in unsaturated zone thickness with increasing cluster m max ; and increasing maximum unsaturated zone thickness with increasing cluster m max .Bloomfield and Marchant (2013) previously noted that such observations are consistent with the findings of Peters et al. (2005), since unsaturated zone thickness is a function of distance to streams.However, in the present study area (Fig. 1) surface drainage is virtually absent from the northern Lincolnshire Limestone that dominates CL1 and is limited over both the Chalk (CL2) and the southern Lincolnshire Limestone (CL4).Instead we postulate that unsaturated zone thickness, and hence m max , is affected by more general catchment characteristics such as extent of outcrop, topography, intrinsic aquifer characteristics and aquifer thickness, which all influence, through unsaturated zone drainage and saturated flow processes, the overall shape of the piezometric surface in the aquifers.For example, of the three aquifers in the study region the Chalk has the most extensive outcrop; it is the thickest aquifer, up to 5 times thicker than the Lincolnshire Limestone; it forms hills up to ∼ 150 m a.s.l., compared to hills about 70 m a.s.l.across the southern Lincolnshire Limestone; and it is associated (CL2) with the largest m max and the longest and largest magnitude droughts.As such, the relationships between unsaturated zone thickness, SGI autocorrelation and hence groundwater drought characteristics are not trivial and appear to reflect a number of fundamental catchment properties and processes that effect groundwater level dynamics and hence groundwater drought phenomena.
Although clustering of groundwater hydrographs is not novel in itself (Winter, 2000;Moon et al., 2004;Upton and Jackson, 2011), this is the first time these techniques have been systematically applied to investigate groundwater droughts.The approach described is generic and widely applicable, and here we briefly highlight some of the methodological considerations, and implications for monitoring and prediction of groundwater droughts.The k-means clustering has been performed on the complete SGI hydrographs, including periods of relatively high groundwater level stand, even though the aim of the hydrograph classification has been to investigate regional variations in groundwater droughts.Yet the resulting clusters have been shown to effectively identify distinct regional groundwater drought responses across the study area.For example, they reflect the major drought history across the study region (Figs. 2 and 7) and identify spatially coherent hydrographs that are consistent with know hydrogeological differences across the study area (Figs.3c and 9a).Eltahir and Yeh (1999) investigated the asymmetry of groundwater hydrographs to high and low groundwater level stands and noted that "droughts leave a significantly more persistent signature on groundwater hydrology than floods".They inferred that this phenomenon was because discharge of groundwater to streams is an efficient dissipation mechanism for wet anomalies and that this discharge is often strongly non-linear.This may explain, at least in part, why the hydrograph classification scheme based on full hydrographs provides such a good basis for analysis of the heterogeneous response of groundwater to drought at the regional scale.However, there is potential for future work to investigate if the hydrograph classification can be improved by focussing on, or giving more weight to, episodes of drought in the SGI time series.
In addition to identifying three clusters of SGI hydrographs -CL1, CL2 and CL4 -that exhibit different characteristic responses to meteorological drivers, the k-means clustering also identified three relatively small clusters of SGI hydrographs -CL3, CL5 and CL6 -where there were trends in the SGI time series; temporal anomalies expressed as anomalous phase relationships between cluster SGI and the regional SGI time series; or relatively poor coherency in SGI time series with a given cluster.In these three clusters it has been inferred that hydrographs may have been variably impacted by anthropogenic factors, such as groundwater abstraction.Although the CA was not specifically designed to identify anthropogenically impacted groundwater hydrographs, the classification scheme could be used to that end since it can differentiate between clusters showing trends superimposed on the regional signals (e.g.CL3 and CL5) and clusters with anomalous phase relationships with the regional signal (e.g.CL6).The presence of a trend in a cluster of hydrographs may be indicative of an anthropogenic impact, for example from unsustainable abstraction (declining trend) or from groundwater rebound (rising trend).
Where there is limited prior information regarding groundwater withdrawals across a region, a not uncommon situation in areas where abstraction is not highly regulated, cluster analysis could be used, either as it has been in the present study based on a set of heuristic rules to identify a suitable number of clusters or in an exploratory manner.If it is used in a more exploratory manner, either hierarchical or nonhierarchical clustering could be undertaken and then clusters searched to identify spatially coherent clusters that show significant downward trends in hydrographs (where significance of trends in a cluster could be tested and quantified using standard tests, such as Mann-Kendall and Sen's slope estimates).Any spatial coherence in clusters exhibiting downward trends may be taken as indicating the presence of potentially unsustainable abstraction.For the purposes of a study where the stationarity of the data is important, if trends in individual hydrographs are already known then either these hydrographs can be removed from an analysis or the trends could be identified and removed prior to standardisation and clustering of the hydrographs.
It has been shown that there can be pronounced differences in the characteristics of multi-annual drought episodes between aquifers within a region (Fig. 9a).During multi-annual droughts some clusters temporarily go out of drought conditions while others will continually show deepening drought conditions over 2 or more years, and some clusters stay in groundwater drought for many months after groundwater (and meteorological) drought has ceased in other clusters.If observations such as these or similar ones can be made for a region, they may have important implications for monitoring groundwater droughts and water resource management in multi-aquifer (cluster) systems.For example, at the end of a drought, sites in more quickly responding clusters may act as leading indicators of the end of groundwater drought at sites in more slowly responding clusters.In addition to the implications for groundwater monitoring particularly during long droughts, if there is sufficient understanding of regional variations in groundwater responses (i.e.relative differences in the timing and intensity of groundwater drought between different aquifers in a region or between sub-regions within an aquifer), then this understanding could be used to inform appropriate groundwater water resource management strategies and so may enable some of the worst impacts of the groundwater drought to be mitigated.
More generally we see a range of possible benefits to clustering groundwater hydrographs.For example, "sentinel" boreholes within each cluster, those that are closest to the mean behaviour of a group, could be identified and used as indicative of the groundwater response of a wider area.Missing data is a common issue with groundwater hydrographs, and clustering techniques could potentially be used to identify suitable boreholes from which groundwater levels could be infilled.However, more importantly, clustering could be used in combination with groundwater models to aid the prediction of groundwater droughts.A range of techniques can be used to model groundwater hydrographs at a site, i.e. non-distributed groundwater models, including statistical models (Ahn, 2000;Bloomfield et al., 2003), artificial neural network models (Sreekanth et al., 2009) and "blackbox" models (Mackay et al., 2014).The hydrograph cluster analysis could be used in combination with any of these techniques for groundwater drought forecasting.For example, forecasts of groundwater levels 1 to 3 months out are currently undertaken in the UK for selected sites using a black-box, lumped-parameter model (Jackson et al., 2013;Mackay et al., 2014;Hydrological Outlooks, 2015) driven by probabilistic estimates of future rainfall.Regional inferences of future groundwater levels are then based on qualitative interpretations of the individual sites.Applying similar modelling systems to mean cluster hydrographs that are representative of spatially coherent regions of groundwater drought response instead of individual site-specific hydrographs could enable more rigorous forecasts of the spatial distribution of groundwater drought.

Conclusions
Cluster analysis when applied to SGI time series of consistent length for multiple sites across a region has been shown to provide a robust approach to the regional analysis of groundwater droughts.In the present study an agglomerative hierarchical complete-linkage strategy and a k-means clustering strategy were tested.The k-means clustering was found to be most suitable.However, for any given case study a range of non-hierarchical algorithms and hierarchical classification schemes should be explored to see which is most appropriate.
A heuristic, rule-based approach was found useful in guiding the selection of the optimal number of clusters, where the rules applied prior knowledge of the hydrogeology of the study area, including information related to spatial variations in catchment and aquifer characteristics.For the present case study, both non-hierarchical algorithms and hierarchical classification schemes provide better clustering of SGI time series than a simple three-fold classification simply based on geology alone, with the k-means clustering providing the best clustering.Membership of the resulting k-means clusters is shown to be dominated by hydrogeological factors, and the effect of heterogeneity in precipitation over the study area on cluster composition is inferred to be negligible.
The clusters successfully discriminate different responses to groundwater drought, both in terms of drought metrics for the complete time series and with respect to the detailed response of sites in each cluster during episodes of major multiannual drought.Groundwater drought characteristics can be linked, through the autocorrelation structure of cluster hydrographs, to the distribution of unsaturated zone thickness.This reflects the role of a range of catchment and aquifer properties and processes that influence groundwater level dynamics, including topography, aquifer thickness and extent of outcrop, unsaturated zone drainage characteristics and saturated groundwater flow.
This approach to groundwater hydrograph clustering is flexible, can be applied in a wide range of hydrogeological settings where suitable hydrographs are available, and enables spatially variable responses of groundwater to drought to be quantified.

Figure 1 .
Figure 1.Case study area (left panel) and simplified geology map (right panel) showing locations of the observation boreholes.Cross section (bottom panel) illustrating the stratigraphy-depth relationships between the three major aquifers in the study region: the Lincolnshire Limestone, the Spilsby Sandstone and the Chalk.

Figure 2 .
Figure 2. (a) SPI-SGI correlation as a heat map, (b) mean SPI 12 time series and (c) mean SGI time series for all 74 hydrographs.

Figure 3 .
Figure 3. (a) Cluster dendrogram for hierarchical classification (k = 6) of SGI time series, (b) map showing the distribution of sites by clusters based on hierarchical classification (k = 6), and (c) map showing the distribution of sites by clusters formed by k-means clustering (k = 6).

Figure 4 .
Figure 4. RMSSD as a function of the number of clusters for the hierarchical and non-hierarchical k-means clustering algorithms and for a three-fold classification based on geology alone.

Figure 5 .
Figure 5. Heat maps of Pearson correlation between SGI and SPI for q = 1 to 36 months and for lags up to 5 months.Maximum correlation is denoted by the closed black circles.

Figure 6 .
Figure6.Mean SPI times series for each of the k-means clusters based on the accumulation period q max for each cluster.The black line is SPI based on gridded precipitation series for sites in a given cluster and the red line is SPI for the mean rainfall across the whole study area based on the different aggregation periods, q max , for each cluster.

Figure 7 .
Figure 7. Mean SGI time series for each of the six k-means clusters.

Figure 8 .
Figure 8. Correlograms for each of the mean SGI time series (bold) and individual site time series (grey) for each of the six k-means clusters showing variation in the autocorrelation function (ACF) for lags up to 60 months.

Figure 9 .
Figure 9. Heat maps showing (a) SGI varying with time for all 74 sites as a function of the six k-means clusters and (b) correlations between all pairs of sites sorted as a function of the six k-means clusters.

Figure 12 .
Figure 12.SGI autocorrelation (m max ) as a function of unsaturated zone thickness.

Table 1 .
Summary of features of the six k-means clusters.

Table 2 .
Summary of drought event statistics for clusters C1, C2 and C4.

Table 3 .
Summary of the 1988Summary of the  -93, 1995Summary of the  -98 and 2011-12 -12drought events for clusters CL1, CL2 and CL4 (where D event , M event and I event denote indices for drought event duration, magnitude and intensity respectively).