Exploring groundwater drought responsiveness in lowland post-glacial environments

Groundwater drought response to meteorological forcing depends on initial hydrological conditions. This makes it difficult to characterise groundwater droughts and identify the drought vulnerability of aquifers. The objective is to increase the understanding of groundwater memory and response to meteorological forcing in lowland post-glacial environments. Eighty-one groundwater hydrographs are analysed, using the standardised groundwater level index (SGI) and the precipitation index. Memory and response times are assessed using auto- and cross-correlation functions. Response time is estimated by comparing two approaches: (1) the traditional use of the maximum cross-correlation, and (2) the alternative use of the cross-correlation slope. Results are interpreted for different hydrogeological settings. The analysis showed that sand aquifers have the longest memory and response times, particularly in confined settings where the memory could be over 4 years. Silts and tills have relatively short memories and response times, at less than 1 year, though median values for silt are higher than for unconfined sand aquifers. In this study, estimating response time using the correlation slope is superior at capturing the initial response time of groundwater to precipitation. However, the results showed that groundwater anomalies in lowland post-glacial environments are sometimes more influenced by climate teleconnections than concurrent forcing. This emphasises the need for a holistic approach for the characterisation and projection of groundwater drought, as it develops in simultaneous response to meteorological forcing at different timescales.


Introduction
A groundwater drought is a prolonged period of abnormally low groundwater levels, which occurs naturally in any climate (Van Lanen and Peters 2000). Hydrological drought develops in response to a meteorological drought, primarily caused by a period of lower precipitation rates than normal (Peters et al. 2006). Groundwater drought development and duration are functions of hydraulic properties which affect recharge processes and water storage (Tallaksen and Van Lanen 2004;Van Lanen et al. 2013;. Due to the high spatial variation of hydraulic properties, drought characteristics (i.e., propagation, duration, severity) can vary significantly, even for nearby locations (Bloomfield and Marchant 2013;Lloyd-Hughes 2013;Haslinger et al. 2014;Kumar et al. 2016;Hellwig et al. 2020).
Spatial differences in drought duration and development, regardless of the spatial extent of meteorological drought, also exist in other parts of the hydrological cycle due to groundwater influence-for example, variable baseflow influences spatial differences in drought development in streams Brooks et al. 2021), and drought development in soil moisture is affected by variable thickness of the capillary fringe and water-table depth (Han et al. 2019; Martínez-de la Torre and Miguez-Macho 2019). However, groundwater can also mitigate drought impacts in other parts of the hydrological cycle. Groundwater can greatly increase the drought resilience of streams and soil moisture (Sharma and Panu 2012;Martínez-de la Torre and Miguez-Macho 2019), subsequently impacting, for example, vegetation stress, crop yields and hydropower production (Mishra and Singh 2010;Kleine et al. 2021).
The ability of groundwater to mitigate drought impacts to society, streams and soil moisture, is reduced or lost when the drought has propagated to groundwater (Kleine et al. 2021); therefore, it is important to quantify the spatial variation of groundwater drought response. The variability in response can be described by changes in dynamic (i.e. active) groundwater storage as a function of initial hydraulic condition (Thober et al. 2015;Arnal et al. 2017;Sutanto et al. 2020;de Lavenne et al. 2022). Some terms used to describe this phenomenon, in terms of drought but also recharge response, are memory (Mangin 1984;Imagawa et al. 2012), persistence (Barco et al. 2010;Duvert et al. 2015) and responsiveness (Van Lanen et al. 2013). The term memory is used in this paper, and it is here defined as the duration of influence from antecedent groundwater level anomalies.
The autocorrelation function of a hydrograph from a groundwater observation well allows for the characterisation of memory of the groundwater system. Autocorrelation is a popular approach applied to karst spring discharge (e.g. Mangin 1984;Massei et al. 2006). More rarely, it has been applied to characterise bedrock (Lee and Lee 2000), coastal aquifers (Tularam and Keeler 2006), silt loam (Eltahir and Yeh 1999), and alluvial aquifers (Imagawa et al. 2012;Duvert et al. 2015;Le Duy et al. 2021). Moreover, the propagation of the precipitation signal to groundwater can be characterised using the cross-correlation function, comparing the precipitation time series to the groundwater hydrograph. Commonly, the lag or accumulation of precipitation with the maximum cross-correlation coefficient is used as the estimated propagation time of precipitation to groundwater, hereafter called response time. However, the cross-correlation function can fluctuate (as described in some studies of karst springs, e.g. Massei et al. 2006;Zhang et al. 2013;Dufoyer et al. 2019) or include a plateauing effect (e.g. shown in a drought study, Bloomfield and Marchant 2013; and in a karst study, Delbart et al. 2016). Further research is needed to explore the application of statistical methods for characterising the response time of groundwater systems with different hydraulic properties.
Unfortunately, data on hydraulic properties related to groundwater recharge and storage (e.g. hydraulic conductivity, storativity, transmissivity), which are correlated to groundwater drought response time and memory (Tallaksen and Van Lanen 2004;Hellwig and Stahl 2018;Hellwig et al. 2020), are seldom available. Similarly, few hydrographs are available because observation wells are often sparsely distributed. Therefore, it can be valuable to connect groundwater memory to typical hydrogeological settings that are expected to have different hydrogeological properties, as done in studies on till and alluvial and outwash plain aquifers of sand and gravel by Ferris et al. (2020) and Boutt (2017).
The purpose of this study is to provide a better understanding of the temporal dynamics that influence the development and propagation of drought, and to increase the understanding of the variability of groundwater drought development in lowland post-glacial environments. Groundwater resources in such cold and humid climates may be more vulnerable to climate-change-induced water deficits than in arid environments . These climate zones are predominant at high latitudes in the Northern Hemisphere (Beck et al. 2018). To the authors' knowledge, there are few studies (Bloomfield et al. 2015;Brakkee et al. 2022) on groundwater drought propagation in temperate and cold climate regions that focus on aquifer variability.
The objective is to (1) analyse the temporal extent and variability of drought memory and response time, and (2) increase the understanding of the relationship between groundwater drought development and the environment. The hydrogeology in large parts of these regions is characterised by unconsolidated Pleistocene post-glacial deposits (Ehlers et al. 2010;Richts et al. 2011). Accordingly, different study sites from Sweden (25 wells from 7 sites) and Finland (45 wells from 9 sites), and the lowland coastal area of the Lower Fraser Valley (11 wells from 1 site) in British Columbia (BC), Canada, are chosen for this study. The wells are classified according to hydrogeological setting, based on information considered as proxy data for dynamic groundwater storage (Tallaksen and Van Lanen 2004). These include sediment type/grain size, since this can be used to approximate hydraulic conductivity (Fetter 2002) and confinement, since glaciated environments commonly contain clay-rich layers that form aquicludes (Ferris et al. 2020). These settings are typical of lowland post-glacial environments.
Groundwater drought memory and drought response time are first characterised for individual wells, and then compared and connected to the different hydrogeological settings. In this study, autocorrelation functions are calculated from the standardised groundwater (level) index (SGI) time series and used to assess the aquifer memory of groundwater level anomalies. The groundwater drought response time to precipitation is evaluated using the cross-correlation function, which is calculated between the standardised precipitation index (SPI) and SGI. Hence, it represents the groundwater response time to precipitation anomalies. In this study, the traditional approach to estimating response time, by using the maximum crosscorrelation, is compared to an alternative approach based on the cross-correlation slope.

Hydrogeological setting
The study sites are in lowland post-glacial environments, which persist in large parts of northern Europe and northern North America that were affected by Pleistocene glaciations (Ehlers et al. 2010). Sweden and Finland are located in northern Europe, and the Lower Fraser Valley in southwestern BC, Canada (Fig. 1a).
The landscape in Sweden and Finland was strongly influenced by millions of years of peneplanation and glacial erosion, resulting in the flat or rolling topography (Amantov 1995). Elevations are below 500 m above sea level (masl) in most of the region (Fig. 1b), except for the Caledonide mountain range along the Sweden-Norway border. The most recent glaciation in Sweden and Finland was the Last Glacial Maximum, which covered both countries completely (Stroeven et al. 2016). The stratigraphy is heavily influenced by the different stages of the Baltic Sea, and the hydrogeological setting can be roughly divided into two groups: supra-aquatic or subaquatic sediments, or above and below the ancient shoreline (Stejmar Eklund 2002). Generally, the crystalline bedrock is overlain by till and/or glaciofluvial deposits, which are covered by Holocene and varyingly active deposition of bog, Aeolian or alluvial sediments. In regions affected by glacio-lacustrine and glacio-marine activity, ice-marginal and glaciofluvial deposits are additionally interlayered with glaciolacustrine, glaciomarine and marine sediments of fine sand, silt and clay (Stejmar Eklund 2002). These glacial deposits are relatively thin, with an average thickness of around 10 m, consisting mainly of till (though in Scania in Sweden, till deposits can be 50 m thick). However, glaciofluvial and ice-marginal deposits and fracture valley saprolites can have thicknesses of more than 100 m (Fredén 2009;Stroeven et al. 2016;Ymparisto 2019 Ehlers et al. (2010), and scaled according to the maximum elevation (c) Due to the low storage capacity of the crystalline bedrock and the presence of thin and discontinuous supra-aquatic eskers, till and fine sediment, socioeconomically important aquifers in Sweden and Finland are found mainly in deltas, extensive subaquatic eskers, outwash plains and in endmoraines (Knutsson 2008;Kløve et al. 2017).
The Lower Fraser Valley forms the southwestern corner of the BC mainland and the adjoining northwestern corner of Washington State in the United States. It consists of flattopped and gently rolling hills which are separated by wide valleys. Elevation ranges from 15 to 122 masl. The Lower Fraser Valley is surrounded by mountain ranges to the north, east and south, and the Salish Sea to the west which connects to the Pacific Ocean (Fig. 1c). Hence, the surrounding hydroclimate is much more complex than in Sweden and Finland (Moore et al. 2010).
This lowland region was repeatedly invaded by glaciers from the adjacent high mountains during the Pleistocene, the most recent of which was the Fraser Glaciation during the Late Wisconsinan (Armstrong 1981;Clague and Ward 2011). Earlier Pleistocene glaciations are poorly preserved in the Quaternary stratigraphic record as they were eroded away during the Fraser Glaciation (Clague and Ward 2011). Approximately 300 m of alternating glacial and nonglacial sediments accumulated within the Fraser Lowlands (Clague 1994). These Quaternary deposits form a complex system of unconfined and confined aquifers comprised of glaciofluvial sands and gravels, lodgement and flow tills, glaciolacustrine clayey silt, silty clay and fine sand, and interbedded marine and glaciomarine sediments. Present-day (Holocene) geomorphological processes continue to shape the landforms observed around the Fraser Lowlands, including deltaic, fluvial, floodplain, bog and mass wasting events.

Climate
The climate is humid cold in northern Sweden and Finland and humid temperate in southern Sweden, according to the Köppen-Geiger climate classification (Beck et al. 2018). The climate zones are generally reflected in the seasonality of groundwater levels (Nygren et al. 2020). In the northern cold climate region, groundwater levels decline in winter during snowfall and frost conditions, and the main annual groundwater recharge event coincides with the spring snowmelt period. In the southern temperate climate region of Sweden, groundwater levels decline during the active growing season. There, groundwater recharge occurs mainly in winter, in response to rainfall and low or nil evapotranspiration rates (Nygren et al. 2020). From the north (cold climate) in a southward direction, the influence of the (short) summer growing season on groundwater level seasonality generally increases. From the south (temperate climate) in a northward direction, the influence of winter frost, snowfall and snowmelt on groundwater level seasonality generally increases (Nygren et al. 2021).
In the Lower Fraser Valley, the climate is humid temperate according to the Köppen-Geiger climate classification (Beck et al. 2018), and groundwater recharge is dominated by winter rainfall (Allen et al. 2014). Similar to the temperate climate region in southern Sweden, the annual groundwater level minima occur after the summer recession, in September and November. The following period of low temperatures lead to increasing groundwater recharge and rising groundwater levels, and the annual groundwater level maxima occur in February or March (Allen et al. 2010).

Data
Precipitation data were provided as the monthly sum [mm] on a 0.5° grid by the University of East Anglia Climatic Research Unit, UK (version TS 4.04; Harris et al. 2020). The global dataset underestimates the median annual sum of precipitation in Sweden and Finland, compared to regional models (Nygren et al. 2021); however, the variability in precipitation rates is similar (Nygren et al. 2021), which is sufficient for the purpose of this paper. A comparison between the global dataset and the main climate station in the Lower Fraser Valley, Abbotsford International Airport, showed that the global dataset is also appropriate to capture meteorological drought in the Lower Fraser Valley (Appendix 1).
The Finnish Environment Agency (SYKE), the Geological Survey of Sweden (SGU) and the Ministry of Environment and Climate Change Strategy of British Columbia (BC ECCS), provide groundwater level hydrographs for the 81 observation wells used in the study. These observation wells were selected primarily based on a period of record and completeness of data. The time series contains a minimum of 15 years of biweekly observations, with no more than 4 months missing. Hydrographs classified with poor data quality by the national agencies, or with irregular time series, were excluded from analysis. In Sweden, there are seven study areas with altogether 25 wells; nine study areas in Finland with altogether 45 wells; and a single study area in the Lower Fraser Valley with 11 wells.
Many observation wells in Sweden were abandoned in 2017, and the Finnish groundwater dataset lacks data after 2017. Therefore, the data period 2000-2015 is used in the analysis of the 70 wells there. The surface elevation at all wells is lower than 500 masl, and the median groundwater depth is between 0 and 18 m below ground surface (mbgs). The wells are screened in unconsolidated Quaternary sediment of primarily glacial origin, but some wells were drilled in Holocene alluvium or Aeolian sand. In BC, only observation well data after 2005 were used since the hydrographs contain features related to the data collection method: prior to 2005 groundwater levels were measured manually at monthly intervals, while post-2005, wells were equipped with dataloggers that record hourly. Thus, the period of analysis is 2005-2020 for the 11 BC wells.
Groundwater level time series from all observation wells were classified according to the hydrogeological setting in the post-glacial environment. The hydrogeological setting was defined by sediment type (silt, till or sand), and wells in sand were further separated by confinement (confined sand or unconfined sand). Since the hydraulic conductivity is largely determined by the minimum grain size, sandy gravel and gravelly sand deposits were classified together as sand. Further, where information on confinement status did not exist, it was inferred from proxy data such as surface sediment maps, or from metadata from wells in the same aquifer. Semi-confined aquifers were classified as confined, as it is difficult to deduce the extent of the confinement from proxy data.
For the time series analysis, the results from each different group are discussed in detail using eight examples from wells in the different settings. The wells were chosen to represent within-group variability in seasonality and groundwater level depth.

Anomalies in groundwater level and precipitation
For drought analysis, groundwater level and precipitation data were standardised to enable the comparison of anomalies across regions with different climates. Median monthly groundwater levels from 81 observation wells were normalised to the SGI following the method developed by Bloomfield and Marchant (2013). The SGI is a measure to describe temporal variability in groundwater levels, as related to the mean. It indicates a continuous change between wet, normal and dry conditions. The method can be applied to spatially distributed data and then used to evaluate the spatiotemporal distribution of groundwater drought development and severity. In this study, the method is used to study the characteristics of drought memory between sites.
The distribution of groundwater levels is not normal, but commonly has a gamma distribution, or can be skewed by extreme values (Bloomfield and Marchant 2013). This variability is largely site-dependent but can also vary by month (Bloomfield and Marchant 2013). Therefore, groundwater levels were transformed to a normal distribution using the nonparametric inverse cumulative distribution function. To exclude normal seasonal fluctuations in the standardised index, SGI values were assigned based on the rank within a calendar month for each well hydrograph given. This was done using the qqnorm function from the statistics package in R version 4.1.1.
For the analysis of meteorological drought, precipitation was normalised to the SPI according to the method developed by McKee et al. (1993), using the SPI function of the SPEI package (Vicente-Serrano et al. 2010) in R version 4.1.1. The gamma distribution was used in this study, as recommended for calculating the SPI in Europe (Stagge et al. 2015). Since the precipitation signal is attenuated, lengthened and pooled before reaching groundwater (Changnon Jr 1987), SPI is calculated from cumulative precipitation data from 1 to 60 months before standardisation. Each accumulation period then produces a SPI time series representative of the cumulative months of precipitation. The different accumulation periods of SPI are referred to as SPI-k. SPI-1 is SPI calculated from the cumulative precipitation of one month, SPI-2 for 2 months and so on.
The SGI and SPI values represent anomalies in groundwater levels and precipitation rates, respectively. Indices between -2 and -1 represent drought conditions, -1 to -0.5 represent dry conditions, -0.5 to 0.5 represent normal conditions, 0.5 to 1 represent wet conditions, and 1-2 represent 'flood' or very wet conditions. Hence, 'normal' precipitation amounts or groundwater levels are compared against wet and dry conditions.

Time series analysis
The autocorrelation function applied to a time series defines the degree to which the time series is dependent on itself for different lag times. It is quantified using a correlogram (Diggle 1990), which is the autocorrelation function coefficient for r k with lag k [months] defined as: plotted against lag k. The autocovariance C k is defined as: where n is the number of time steps in the record, k is the lag, t is the time step [month] in the hydrograph record, and SGI a is the average SGI. Hence, the function quantifies the effect of antecedent conditions, i.e. the length of the 'memory effect' (Massei et al. 2006). In this study, autocorrelation of the SGI is calculated for lag times from null to 60 months. (1) For a variable without a memory effect, the function decreases rapidly to or near zero for short lag times. For a variable with a long memory effect, the function remains high for long lag times. For example, for SGI correlated against itself for lag k = 0 the correlation is always perfect (r k = 1), and the further back in time, the lower the correlation. Commonly, the memory effect for a time series is defined as the lag time for which the autocorrelation function is above a threshold value. Correlation coefficients above the threshold are considered to be significant, while values below are interpreted to be random noise. Bloomfield and Marchant (2013) applied a threshold of 0.11 to autocorrelation functions of groundwater anomalies, based on the estimated confidence interval of the autocorrelation. The confidence interval for significant autocorrelation at the 5 % level is here defined as 1.9599 / √ n , where n is the record length. Hence, if r k > 1.9599 / √ n , the correlation is significant at approximately the p < 0.05 level. In this study, the resulting threshold for significant autocorrelation is 0.17, which was used to estimate groundwater drought memory. However, the decrease in autocorrelation is seldom linear, and alternatives to the threshold application have been suggested (Massei et al. 2006;Delbart et al. 2016). To assess what information may be extracted from a groundwater anomaly standpoint, the entire function (for k null to 60 months) is considered together with the memory estimates.
The cross-correlation function quantifies the similarity of one time series to another. Here, this is the similarity between groundwater level and precipitation anomalies for different accumulation periods. Similar to the autocorrelation, the cross-correlation function is quantified with a correlogram showing the cross-correlation function coefficient against the different SPI accumulation periods. There is a relationship between the time series if the cross-correlation function is asymmetric and has an absolute maximum or minimum value for lags above 0 (Box et al. 1994). In this study, Spearman's rho is used to calculate the cross-correlation. The entire function is considered also for the cross-correlation function, but only the significant output is evaluated.
The SPI accumulation period with the maximum crosscorrelation coefficient is commonly taken as an estimate of the groundwater response time to precipitation. However, in this paper, the traditional approach is compared to an alternative approach to estimate response time. The aim of the alternative approach is to identify the first peak in the correlation, and the start of the plateau, where there is a plateau-feature. The first step is to calculate the first-order derivative (i.e. the slope) of the cross-correlation function. The response time is then estimated as the previous accumulation period for when the derivative is equal to or less than zero (zero vs. declining slope), for positive correlations and for accumulation periods higher than 2 months (for cases where correlations are initially negative), written as: where rt CCF � is the response time calculated by the first derivative of the cross-correlation function, k is the accumulation period in SPI-k and CCF is the cross-correlation coefficient.
To account for cases with short response times, the traditional approach is combined with the alternative approach. The response time is thus defined as the minimum lag of the maximum cross-correlation coefficient compared to rt CCF � , written as: where rt is the response time. To separate the response time calculated from the two approaches, the approach using only the maximum cross-correlation coefficient is hereafter referred to as the response time by max. correlation. The alternative approach using the combination in Eq. (4) is referred to as the response time by correlation slope.
The results from the correlation analysis, including memory and response times, are studied in detail for eight wells (see information in Table 1, and corresponding SGI values are given in Appendix 3). There are two examples from each hydrogeological setting, which are composed of wells in (1) confined sand aquifers, (2) unconfined sand aquifers, (3) silt and (4) till.

Groundwater level seasonality
The average groundwater level seasonality from different settings is shown in Fig. 2. The eight example wells are highlighted, and each represents an aquifer with a higher vs. lower amplitude in seasonal variation. There is a clear distinction in amplitude between the settings, whereby the confined sand aquifers have a small seasonal variation (0.5 m) compared to unconfined sand aquifers (1-2 m). One of the confined aquifers, monitored by well 357 in the Lower Fraser Valley, has a much stronger seasonality, with an annual groundwater level fluctuation of around 1.5 m (Fig. 2a). Table 3 in Appendix 2 contains no information to set it apart from other wells in confined sand aquifers; however, borehole data from the Ministry of Environment in BC describe the confining clay layer as thin, overlaid by gravel (Ministry of Environment and Climate Change Strategy of British Columbia 2005); thus, the aquitard may be locally leaky, which could explain the measured difference in seasonality. There is minor seasonal variation between different silt hydrographs, but the amplitude is low (<1 m), similar to what is measured in confined sand aquifers. The till hydrographs have one group with low seasonal variability (0.5 m) and one with a strong seasonal dependency (1-2 m).

Autocorrelation
Most aquifers have a memory shorter than one year (Fig. 3). However, sand aquifers have a large variability in memory, especially confined sand aquifers, ranging between 5 to 53 months with a median of 11 months ( Fig. 3c; Table 2). The example wells, however, represent aquifers with similar memory estimates, of 9 ( Fig. 3a) and 10 ( Fig. 3b) months.
Unconfined sand aquifers have a smaller range in memory, between 1 and 21 months, with a median of 5 months ( Fig. 3f; Table 2). The example wells represent this range clearly in this case, where Fig. 3d represents an aquifer with a 16-month memory, while Fig. 3e represents a 7-month memory.
There are few wells in silt; hence, the range in the autocorrelation result is comparatively small. The estimated memory for silts is between 6 and 12 months, with a median that is slightly higher than for unconfined sand aquifers at 8 months (Fig. 3i, Table 2). The two example wells in silt are from the same location, Kälviä, in Finland, giving silt memory estimates of 8 months (Fig. 3g,h). The decline in the autocorrelation function is slightly different between the examples, which indicates that the memory estimate is sensitive to the threshold.
Tills have a longer range in memory compared to silt, from 1 to 13 months, but the median is the smallest at 3.5 months ( Table 2). The example wells give short memory estimates for the different tills, 6 months ( Fig. 3j) and 2 months (Fig. 3k).
Some of the examples have a secondary peak in the autocorrelation function, near or above the threshold value of 0.17 around lags between 35 and 45 months (Fig. 3g,h,k), which also appears in the median autocorrelation function for some settings at 40-45-months lag (Fig. 3c,i).
Further, in some cases, the autocorrelation has a variable slope, primarily when the autocorrelation has a gentle decline and produces a longer memory estimate. Initially, the decline is linear and the slope is high (indicated by the diagonal blue line in Fig. 3), followed by a decrease in slope (e.g. examples Fig. 3b,d,g,h,j, and the median in Fig. 3c,l).

Cross-correlation
The cross-correlation between SGI and SPI is high for a wide range of accumulation periods (Fig. 4). For most examples, there is generally not a single peak which would suggest one absolute response time that trumps other accumulation periods; the exception is Fig. 4j.
The cross-correlation functions for confined sand aquifers have similar correlation coefficients for many different SPI-n values, and local peaks (Fig. 4c). Figure 4a represents an aquifer with a low maximum correlation (compared to the other examples) between SGI and SPI (~0.4). Due to the plateau feature in the correlation, the response time by correlation slope gives a lower estimate (16 vs 6 months). The other example in a confined setting has similarly low cross-correlations for short response times, but a higher correlation of ~0.7 for the response time by max. correlation (60 months, Fig. 4b).
For the unconfined sand aquifers, the range is again too large to establish a general indication for the setting (Fig. 4f). Both examples have a cross-correlation of around 0.5 for many SPI-n. The difference is that the maximum correlation is higher (~0.7) in Fig. 4e.
The six silts have similar cross-correlations with two separate peaks (Fig. 4i). As a result, the example silts have very different response time estimates by max. correlation; for Fig. 4g, the response time estimate by max. correlation occurs within the first peak, and for Fig. 4h, within the second.
Tills have similar initial cross-correlations, then the range increases (though not as much as for sand aquifers; Fig. 4l). Figure 4j represents the cross-correlation for a till with identical response times by max. correlation and correlation slope. Although the response time by max. correlation and correlation slope are identical also for the other till example (Fig. 4k), the cross-correlation is more similar to the silts.
Overall, there are two periods which impact the SGI response time: 5-20 and 55+ months for confined sand aquifers (Fig. 4c), 5-15 and 50+ months for unconfined sand aquifers (Fig. 4f), and 5-10 and 50+ months for silts (Fig. 4i). Tills are the exception, for which the cross-correlation indicates a general response time around 1-10 months (Fig. 4l), though the result for one of the example tills (Fig. 4k) indicates it is correlated to two accumulation periods (4-20 and 55+ months). Figure 5 shows the relationship between estimated response time using the two different approaches and memory. In general, estimated response times are different from the memory of confined sand aquifers: in comparison to the memory, the max. correlation approach gives longer estimates, and the correlation slope gives shorter estimates (Fig. 5a). However, the density plots show that the approach using the slope is, in general, a better match to the distribution of the memory (Fig. 5b).

Comparison between memory and response time
For unconfined sand aquifers, the response time by max. correlation is nearly always longer than the memory (Fig. 5c). The density plot further emphasises that the max. correlation estimate has little relation to the memory (Fig. 5d), whereas the alternative approach using the correlation slope has a near linear relationship to the memory, and a similar distribution.
For silt and till, memory and response time are similar and near the 1:1 ratio for both approaches of estimating response time (Fig. 5e,g). However, the response time by max. correlation produces some outliers where the response time estimate by max. correlation is, as for sand aquifers, much longer than the memory. The silt outlier is represented by observation well K5, which is one of the example wells ( Figs. 3h and 4h). The density plots show a tendency for the response time by correlation slope in silts to be shorter than the memory, by 1-3 months (Fig. 5f).

Discussion
The purpose of this study was to increase the understanding of groundwater drought memory and response to precipitation, with a special focus on lowland post-glacial environments. For this purpose, the autocorrelation of groundwater level anomalies was evaluated together with the cross-correlation of groundwater to precipitation anomalies. This was done from the perspective of the four different hydrogeological settings defined as (1) confined sand aquifers, (2) unconfined sand aquifers, (3) silt and (4) till.

Autocorrelation
In drought research, the autocorrelation reflects the memory (self-correlation) of the aquifer, which may also be referred to as the responsiveness (e.g. Van Lanen et al. 2013) or the persistence (e.g. Sharma and Panu 2012). An aquifer with a longer memory has a higher capacity to store water (Hellwig et al. 2020), is subject to slow processes, has a poor connection to the atmosphere, or a combination of all (Duvert et al. 2015). In theory, this means that the aquifer is less sensitive to local and short-term changes in precipitation and evapotranspiration. Further, a developed drought persists for longer in an aquifer with long memory, compared to an aquifer with short memory (Bloomfield and Marchant 2013;Van Lanen et al. 2013;. In this study, only confined sand aquifers have memories longer than 21 months, with a range between 5 months and over 4 years. This shows that confined sand aquifers can respond very differently to drought, and that they may have very long groundwater drought memory. Thus, they can take years to respond to and recover from drought. Due to the range, it is difficult to attribute the large difference in memory for confined sand aquifers. Additional factors than those listed in Table 1, such as hydraulic properties, vadose zone thickness, connectivity to stream, land use and extraction, are influencing the memory (Eltahir and Yeh 1999; Kløve et al. 2014;Haaf et al. 2020;Hellwig et al. 2020;Riedel and Weber 2020;Le Duy et al. 2021). This kind of data is lacking on the regional scale. The examples from confined aquifers give similar estimates of 9 and 10 months, respectively (Fig. 3a,b), so discussing differences between these is not useful.
Groundwater drought memory in unconfined sand aquifers varies from 1 month to over 2 years; hence, in this study, they have a shorter memory and recover faster from drought compared to the confined counterpart. The autocorrelation of the examples gives memory estimates of 16 ( Fig. 3d) and 7 months (Fig. 3e), respectively. This may be explained by the difference in seasonal amplitude of 2.5 m (Fig. 2b), which gives an estimate of the atmospheregroundwater connection. Since groundwater level depth is related to seasonality and also gives an indication of the atmosphere-groundwater connection in unconfined sand aquifers, the difference in depth to groundwater level, of 5 m (17 vs 12 mbgs; Table 1), may explain the difference in seasonal amplitude and in memory.
Silt and till have a low hydraulic conductivity (Fetter 2002;Ferris et al. 2020), which is a probable cause for the lower range and relatively short drought memories (maximum one year) of these settings. Nevertheless, the median autocorrelations for the settings suggest that silts generally take longer to recover from drought than the unconfined sand aquifers, and less time to recover than confined sand aquifers ( Fig. 3i; Table 2). The silt examples have identical memory estimates, so differences between silts cannot be discussed.
Comparatively, tills recover from drought faster than other units in post-glacial environments ( Fig. 3l; Table 2). The memory estimates of the till examples are slightly different. One example till has a memory of 7 months (Fig. 3j), while the other has a memory of 2 months (Fig. 3k). While both examples have shallow water tables, the one with a longer memory is deeper, at 3.5 mbgs, in comparison to the other well with a water table at 0.6 mbgs (Table 1). Paradoxically, the shallower example with a shorter memory has a lower amplitude in seasonality (Fig. 2d), which suggests other factors than water-table depth, such as diffusivity/dynamic groundwater storage, are influencing the autocorrelation.
Some of the example silts and tills (Fig. 3g,h,k), and the median autocorrelation for confined sand aquifers and silt (Fig. 3c,i), have secondary 'peaks' in the autocorrelation at lags around 3-3.5 years. This shows that, commonly, groundwater in post-glacial environments is affected by groundwater drought conditions years earlier, which can be an indication of influence from climate teleconnection patterns. In Sweden, the North Atlantic Oscillation (NAO), Arctic Oscillation and Scandinavian Oscillation influence streamflow with a periodicity of 2-4 years (Uvo et al. 2021). In BC, the El Niño Southern Oscillation influences streamflow with a periodicity of 2-7 years (Perez-Valdivia et al. 2012). Connecting groundwater drought to streamflow drought was outside the scope of this study, yet the two are inherently connected. Thus, it is reasonable that factors influencing streamflow propagate to influence the groundwater level, or vice versa where groundwater discharges to streams. Indeed, Rust et al. (2019) analysed drought frequency in the UK and found that teleconnections strongly influence aquifer drought occurrence at periodicities of 7 and 16-32 years. They contributed up to 40% of inter-annual groundwater storage variability to the 7-year NAO cycles (Rust et al. 2019). Therefore, it is likely that teleconnections in Sweden, Finland and the Fraser Valley influence the autocorrelation of groundwater anomalies, particularly in confined sand and silt.
Another feature shown in Fig. 3 is the variable slope in the decline of autocorrelations with a longer memory. Autocorrelations hence decrease with variable speed, which suggests the autocorrelation function for SGI time series reflects multiple processes, as for karst aquifers (e.g. Delbart et al. 2016). The decline is here interpreted to be a result of different signals superimposed on each other. The initial, linear and often steep decline in SGI autocorrelation is probably a product of direct, local and short-term precipitation  anomalies, i.e. recharge (Eltahir and Yeh 1999), depending on the atmosphere-groundwater connection (Eltahir and Yeh 1999;Duvert et al. 2015). The gentler slope which sometimes follows the linear decline, could be a result of slow, diffuse propagation of precipitation anomalies, depending on the dynamic groundwater storage (Imagawa et al. 2012;Bloomfield and Marchant 2013;Duvert et al. 2015).

Cross-correlation
The cross-correlation between groundwater level and precipitation is used to estimate a response time. Response time calculated by max. correlation is a good estimate for a textbook example of the cross-correlation that has a gradual increase, followed by a peak and a decrease. However, this is uncommon for the time series used in this study, due to irregularities, such as plateauing and multiple peaks, in the cross-correlation. The range in cross-correlation for confined sand aquifers is so large that the median cannot give a good representation of the response time in general (Fig. 4c), as mentioned for the autocorrelation. The examples can instead provide some perspective. One of the examples has a plateau feature (Fig. 4a), and the response time by correlation slope seems a more suitable response time estimate, as that is when near-maximum correlations are first attained and reflect the initial influence from precipitation. However, due to the difference in correlation for the other example (Fig. 4b), the response time by max. correlation seems more representative of drought propagation to groundwater. The variability in cross-correlation and response time is connected to the memory, or inertia, of the aquifer. As for memory, this is influenced by hydraulic properties, dynamic groundwater storage and the groundwater-atmosphere connection (e.g. Duvert et al. 2015). The long response time seen in Fig. 4b is in accordance with the assumed (comparatively) high transmissivity of the aquifer monitored by well V10.
Unconfined sand aquifers have a large variety in crosscorrelation values (Fig. 4f), like in the confined counterpart (Fig. 4c), emphasising that other factors are influencing variability in response time. Similar to the other well in the Fraser Valley (Fig. 4a), Fig. 4e has a plateau-feature, which may be due to the high specific yield of the aquifer (Allen et al. 2010). Due to the similar correlation of the two estimates, the response time by correlation slope is a better approach for Fig. 4e (31 vs 18 months). Figure 4d is more reminiscent of the classic example, though the response time estimates are different. Both are long, which reflects the median groundwater level depth at 17 mbgs. In this case, the maximum correlation is higher (~0.7) than for the response time (~0.5). Hence, the response time by max. correlation gives a better estimate for the aquifer represented in Fig. 4d (31 vs 18 months).
For silts, the cross-correlation shows strong influence by SPI from two different time intervals (Fig. 4g,h,i). Thus, the response time estimated by the correlation slope approach is superior in capturing the first peak and reflects the response time to contemporary precipitation anomalies.
The median cross-correlation for tills is the only setting which has the shape of a textbook example (Fig. 4f,j); hence, the SGI responds to one accumulated precipitation signal, identical for the two approaches. Figure 4k is similar to the examples in silt, and the response time by correlation slope captures the initial response to precipitation.
However, some of the examples (Fig. 4g,h,k) are clearly influenced by a periodicity, which is reflected in the median crosscorrelation for sand aquifers and silts (Fig. 4c,f,i). The first peak in the cross-correlations, at around 1 year, is likely due to forcing from contemporary precipitation anomalies (recharge), while the second peak could be due to influence from teleconnection patterns (Rust et al. 2019), as suggested for features in the autocorrelation function of similar periodicity. This is reasonable, since the influence of climate teleconnections on groundwater is driven by the input, i.e. precipitation. Note that in some cases, as for wells K1, K5 (in silt) and AT8 (in till), the strength of the cross-correlation is similar in the initial plateau as in the second. This suggests that the influence of climate teleconnection patterns on SGI in post-glacial environments could be nearly as high, or higher, as concurrent precipitation anomalies (recharge). Note that this varies on an individual basis.

Comparison of memory and response time
Based on a study of chalk, limestone and sandstone aquifers, Bloomfield and Marchant (2013) concluded that response time and groundwater memory have a linear relationship. Comparing their results to those in this study shows that the response time by max. correlation is not linearly correlated to memory for wells in post-glacial environments (Fig. 6a). The approach produces much longer estimates for response time compared to memory, particularly for wells in sand aquifers (confined and unconfined, see Fig. 5).
The approach using correlation slope, instead of the maximum correlation, produces estimates of response time more similar to memory (Fig. 6b). In general, this approach is far superior at capturing the distribution of memory estimates in this study (Fig. 5b,d,f,h). Outliers, which are fewer with this method, instead have shorter response times compared to the memory. This result is similar to findings by Bloomfield and Marchant (2013) and Bloomfield et al. (2015), where their comparisons also tended toward response times being similar or shorter than the memory. However, the autocorrelation results (Fig. 3) clearly show that memory estimates are very sensitive to the threshold used to determine the significant autocorrelation range. The cross-correlation results show that the appropriate method to estimate response time varies on an individual basis. Additionally, SGI anomalies can be significantly and highly correlated to show the distribution of the estimated response time and memory for all wells in b confined sand aquifers, d unconfined sand aquifers, f silt, and h till many SPI accumulation periods, and not necessarily driven by only one (Fig. 4). The benefits of estimating response time using the correlation slope are largest where the correlation has plateau features (minor differences in the correlation for multiple accumulation periods) and where the setting is affected by inputs containing strong periodicities. In these cases, the response time by correlation slope captures the initial groundwater drought response. Nonetheless, the auto-and cross-correlation functions should be reviewed, to ensure the appropriate interpretation of these estimates.

Approach limitations
One key limitation of the approach in this study is that the data collection process varies between the three countries. Since 2005 for the wells in the Fraser Valley, groundwater level data collection is mainly automatic with an hourly resolution. In Sweden and Finland, the start date of the installation of automatic data is unknown, and daily resolution is available only for some wells since around 2015. In Finland, the highest data resolution is biweekly. Therefore, monthly average values may be skewed disproportionately by extreme values, particularly in Sweden and Finland as compared to the Fraser Valley.
The Lower Fraser Valley has a complex hydroclimate due to its proximity to the sea and surrounding mountains (Moore et al. 2010). This probably has implications for the estimation of SPI, the cross-correlation function and the estimation of the drought response time.
While care was taken to avoid data that were influenced by human activity, it was difficult to certainly rule out such effects. Particularly seasonal pumping for irrigation was hard to identify, which potentially influences drought recovery. The approach of standardising groundwater level and precipitation data is also sensitive to the data record length. The time period used in this study is short (15 years) and this impacts the ranking of groundwater levels and precipitation, which results in the SGI and SPI values. This may also influence the correlation functions. However, since the groundwater level autocorrelation, and therefore also the cross-correlation to precipitation, varies over time (Eltahir Fig. 6 Plot of memory against response time, by a max. correlation and b correlation slope for all wells used in the study (Fig. 5), compared to two previous studies. Note that the data points for Bloomfield et al. (2015) are representative of clusters of wells, and not indi-vidual wells, and that one cluster lacked representative values and is therefore not included here. The diagonal dashed line shows a 1:1 relationship for reference and Yeh 1999; Lee and Lee 2000), the effect on the autocorrelation and cross-correlation may be unavoidable. Lastly, sediment type is a limited proxy for hydraulic conductivity, as this can vary within a sediment type, particularly with depth (e.g. Lachassagne et al. 2011;Ferris et al. 2020).

Conclusions
The purpose of this study was to increase the understanding of precipitation drought propagation to groundwater in lowland post-glacial environments. Drought conditions were estimated by calculating the SGI and SPI, following the methods developed by Bloomfield and Marchant (2013) and McKee et al. (1993). Auto-and cross-correlation analysis was used to quantify groundwater drought memory and groundwater response time to meteorological drought for wells in four common post-glacial environments. The estimates of memory were compared against two different approaches to estimates of response time.
1. The results of the autocorrelation function analysis, with a variable slope, suggest that hydrographs from different wells are variably affected by different signals superimposed on each other: • The initial fast decline is interpreted to be the groundwater response to concurrent meteorological forcing (recharge). • The gentler decline at longer lag times is interpreted to reflect the influence from large dynamic (active) groundwater storage.
2. Wells in confined sand aquifers can have a groundwater drought memory and response time over 4 years, and nearly 2 years in unconfined sand, but the range is large and factors such as vadose zone thickness and diffusivity are influencing the variability. Wells in silt and till can have a groundwater drought memory and response time of up to 1 year, and this is attributed to slow processes and the comparably low-dynamic groundwater storage. 3. Many wells, particularly in sand aquifers and silt, showed periodicity in the autocorrelation and crosscorrelation at around 3-4 years. This indicates that groundwater drought in the investigated post-glacial environments is influenced by teleconnection patterns, in some cases nearly as much or more than to concurrent climate forcing. 4. In this study, the approach of calculating response time by correlation slope was generally superior at identi-fying the initial groundwater response to precipitation, compared to the approach using the maximum correlation. This is because the correlation slope approach captures the initial, first, groundwater drought response to the input. 5. The correlation analyses show that memory estimates are very sensitive to the choice of threshold, and that groundwater drought is influenced by precipitation anomalies for many different accumulation periods. Therefore, memory and response time estimates, and the assumption that they are linearly correlated, should be reviewed on a case-by-case basis.
To understand the influence of hydrogeology on groundwater drought memory and propagation, further studies are needed to investigate the change in dominating signals for different study sites and types of drought, but also for nondrought and flood conditions. Particularly, the interaction between streamflow and groundwater drought could provide insight into drought propagation effects. In addition, more research at a smaller scale is needed to understand the variability in drought propagation in these settings, particularly focussing on confined aquifers and including additional features such as aquifer depth and aquifer thickness.

Appendix 1
For this study, meteorological drought was estimated using precipitation data from a global dataset, provided by the University of East Anglia Climatic Research Unit (CRU) (version TS 4.04; Harris et al. 2020). This was for consistency in the methods and also because the Abbotsford International Airport climate station, the main climate station in the Lower Fraser Valley, has data gaps, with the record ending in 2011.
To assess the feasibility of using the CRU dataset to analyse local drought, a comparison with the climate station was performed (Figs. 7 and 8). The CRU data generally underestimate annual precipitation rates, compared to the climate station. These results are in line with findings from a study comparing the CRU dataset to two regional models in Sweden and Finland (Nygren et al. 2021).
However, the variability in precipitation is more important for the feasibility analysis, since the focus of this paper is on drought. Figures 7 and 8 clearly show that precipitation rates for the two datasets have a similar variability; hence, the CRU dataset can be used to characterise meteorological drought patterns in the Lower Fraser Valley.