Integrating storage and spatial variability into shallow groundwater balances: moving towards water security in hard rock coastal areas

Identifying sustainably extractable volumes of water from a groundwater system is vital from a regional planning viewpoint, particularly in regions which experience water scarcity. In Europe, the Water Framework Directive (EU WFD 2000) specifies that abstraction in excess of the overall recharge ‘not needed by the ecology’ be avoided. The Swedish Parliament established an environmental goal in 1999 ensuring public access to groundwater of good quality. However, extremely limited storage due to lack of soil cover and the inherently low-storage characteristics of the hard rock (Engqvist & Fogdestam 1984; Knutsson & Morfeldt 2002) implies that the timing of recharge throughout the year will play an important role in how the local reservoirs are influenced by spatially distributed, private water-supply wells. Thus, tools which account for the limited storage capacity which is prevalent in many of these areas are needed in order to manage groundwater resources appropriately and identify areas which are vulnerable.

Quantifying overall recharge, let alone the portion needed by the ecological system, is problematic and prone to error (Gaye & Edmunds 1996; Gleeson et al. 2009; Chesnaux 2013). However, in terrains with hard crystalline bedrock and limited soil cover, such as much of Scandinavia, the issue becomes much more problematic where field and regional estimates of recharge (Jie et al. 2011) may not be representative of local settings in terrains with elevated levels of heterogeneity (Earon & Olofsson 2018) and deterministic models may be too uncertain and resource-heavy to calibrate locally (Dripps & Bradbury 2007; Earon 2019).

Often lacking traditional productive aquifers (e.g. glaciofluvial deposits), coastal terrains in this region are often characterized by frequent hard rock outcrops, with some glacial till and clay or peat in topographical depressions. Reservoirs in these areas are often limited to the kinematic porosities (namely the porosity of the in situ volume which is connected to the flow network and excluding stagnant pores) of bedrock, which often are in the range of 0.01–0.08% (Singhal & Gupta 1999; Olofsson et al. 2001; Chesnaux et al. 2009) and the thin soil cover, if present. As this soil cover is often isolated to valleys and the surrounding hydraulic conductivity of the bedrock is very low (Engqvist & Fogdestam 1984; Chesnaux et al. 2009; Earon et al. 2015), these soil-filled valleys provide an extremely important function for local residents with privately drilled wells. Due to the low-population density in Swedish, peri-urban areas and the difficulty in connecting residences in these areas to municipal water systems due to limited soil cover, local drilled wells are often the only viable source of sustainable drinking water for many residents. It is thus exceedingly important to account for these small reservoirs both quantitatively and spatially in water balance calculations, as kinematic (or effective) porosity values are often several orders of magnitude higher in common soils (e.g. till) than in crystalline rock (Freeze & Cherry 1979). Thus, as a reservoir, a 1 m thick layer of glacial till could serve equally as well as 100 m of granitic bedrock and will likely act as a recharging reservoir which balances volumes extracted from local drilled wells. Hydraulic heterogeneity and anisotropy in crystalline bedrock make regional characterization based on point data, such as well tests, extremely difficult (Shapiro & Hsieh 1998; Wladis & Gustafson 1999; Earon et al. 2015; Earon 2019). Numerical modelling in these types of terrain is therefore data- and resource-intensive, thus necessitating the use of surrogate models and heuristic solutions for groundwater management, even if only to identify areas where more robust hydrogeological modelling and testing is needed.

Further complicating the matter is the elevated risk for increased chloride concentrations with depth in much of the most populated regions of Sweden (Boman & Hanson 2004), due to these areas being submerged in brackish or saltwater following the last glacial episode and the subsequent presence of fossil seawater. Quantifying impacts to storage via groundwater abstraction should therefore allow for identification of areas with increased risks for groundwater contamination, where greater proportions of maximum freshwater storages are used throughout the year. Many coastal regions in Sweden are experiencing increasing permanent residency (e.g. Tyresö 2017), which increases the risk for salinization from brackish water bodies, in this case, the Baltic Sea. As many of these regions are not connected to Municipal Water and Waste Networks, it is essential to quantify volumes of water which may be extracted in order to avoid increased risk for environmental stress on existing groundwater aquifers. Often, water budgets and sustainably extractable volumes are not estimated or calculated after the fact, increasing the risk that regions are mining groundwater (Voudouris 2006).

Spatial approaches to groundwater balance calculations have been successfully applied to heterogeneous areas as a tool which can improve sustainable groundwater resources use (Batelaan & De Smedt 2007). These approaches, particularly using geographic information systems (GIS), allow for the incorporation of existing geological and meteorological databases with ease. However, traditional tools for groundwater balance and flow calculations often assume an unlimited storage capacity (El-Hames 2004; Portoghese et al. 2005; Voudouris 2006; Batelaan & De Smedt 2007; Tilahun & Merkel 2008; Stoeglehner et al. 2011; Antonakos et al. 2014), with few exceptions (Olofsson 2002). Statistical tools for estimating water quality (Cameron et al. 2018) may fail due to extremely high levels of heterogeneity found in these types of geological environments (Earon 2014; Earon & Olofsson 2018; Earon 2019), although promising advances for using surrogate models are being seen when applying deep-learning algorithms where extensive time series measurements are present (Müller et al. 2020) and geostatistics (Wang et al. 2020). Due to the limited storage capabilities of recently glaciated terrains, precipitation during spring and autumn months could potentially fill such aquifers very quickly, causing any excess precipitation to be lost as either baseflow runoff or overland runoff. While runoff volumes in hard rock terrains are often very high (Spence & Woo 2002), these runoff volumes will usually accumulate in topographically lower pockets of soil, where they likely infiltrate causing these soil pockets to act as a recharging reservoir to the surrounding area, which could mean that simple models for recharge such as precipitation-evapotranspiration balances might be an adequate approximation.

The purpose of this paper is to present a straightforward methodology for integrating a storage parameter with high levels of heterogeneity into a spatial, shallow (∼0–200 m) water balance with the intended application being as a heuristic surrogate for numerical groundwater modelling to aid in groundwater resources management and increase water security in hard rock coastal terrains. Specific aims include identifying how spatial heterogeneity in storage and extraction parameters will influence hydrogeological response in coastal areas. Additionally, the effects of varying the storage parameter, residency and groundwater flow will be investigated in order to assess the sensitivity of these parameters.

The study area, Östhammar Municipality, lies roughly 100 km north of Stockholm, Sweden (Figure 1). The soil cover in the area is limited to glacial and post-glacial clay or till-filled valleys punctuated by frequent bedrock outcrops (Figure 2). Bedrock in the area mainly consists of gneissic granite (1.87–1.96 GA or Giga-annum) or younger granites (1.75–1.83 GA). The area also has some outcropping regions of sedimentary gneisses, metaryolites and metadiorites. Sedimentary gneissic and gneiss-granitic bedrock hydraulic conductivity values in the region have been previously estimated to range from 2 × 10⁻⁶ to 5 × 10⁻¹¹ m/s (Earon et al. 2015).

Data used in this study were obtained from the Geological Survey of Sweden (SGU) (well archive, bedrock type, soil data) (SGU 2017) and the National Land Survey of Sweden (elevation) (Lantmäteriet 2016). The well archive data include data from private wells including location, estimated capacity and usage-type as a vector (point) shapefile. Location of individual wells had a spatial uncertainty of either 100 m, 250 m or unknown. Required input data for the model (Figure 3) were geological soils, bedrock type, well locations and surface elevation. Geological maps (soil and bedrock) were acquired from the Geological Survey of Sweden (SGU) in vector (polygon) format. Soil and bedrock data had a spatial resolution of 1:50,000 and include information such as rock type and estimated age of the bedrock and soil maps showing superficial soil cover. Elevation data were available from the Swedish Land Survey (Lantmäteriet) as continuous data (raster) with a cell resolution of 2 m and an average in-plane accuracy of 0.25 m (Lantmäteriet 2016). Elevation data were resampled to a resolution of 10 m for cell size consistency within the model. A flow chart of data sources, their use and limitations is presented in Figure 3.

This study applied a mass balance of freshwater resources using a 10 m cell size in ArcGIS 10.4.1. This cell size was selected to maintain a model resolution over a fairly large area but still be capable of accounting for important geological features, such as small pockets of till. Geological data were converted to raster sets also with a 10 m cell size, after which outcrops and elevation were used to generate a soil-depth model as presented in Karlsson et al. (2014). Additionally, a separate soil depth model based on interpolated depth measurements in separate soil layers developed by SGU (Daniels & Thunholm 2014) using well archive and other geological data (SGU 2017) was used in order to assess parametric sensitivity within the storage parameter.

For each cell in the model, soil layer properties were inferred (Table 1) based on typical geological strata found in the region (Fredén 1994) and modelled soil depth. For example, if the soil was estimated to be deeper than 1 m and was classified as glacial or post-glacial clay in the geological soils map, it was assumed that the clay would be underlain by 1 m of till. Another example of assumed stratification was that glacially washed sand was underlain by glacial clay and till, thus often separating the storage in the sand layer via a low- to semi-permeable layer as it would likely act as a perched aquifer. Kinematic porosity values were set based on those commonly found in literature (Freeze & Cherry 1979; Singhal & Gupta 1999). Measured surficial fractures in a separate area with similar bedrock (Earon 2014; Earon et al. 2015) were used as an estimate to assign kinematic porosity values to bedrock. Based on the distance to tectonic features identified in the bedrock maps and data presented in Earon et al. (2015), kinematic porosity values were decreased by 25% at distances of 0–50 m from a fracture zone and increased by 50, 25 and 10% for distances of 50–75, 75–100 and 100–150 m, respectively. This is due to the possibility of the fractures very close to the centre of the fracture zone having lower transmissivities (Earon et al. 2015), due to, for example, mineral intrusion. The total maximum storage for each cell was assumed to be equal to the elevation of the cell (m.a.s.l.) multiplied by a factor of 1.5 (Figure 2), which was assumed to be a reasonable reservoir thickness in near-coast regions for avoiding saltwater intrusion or relic salt water intrusion, where likelihood of chloride contamination increases significantly with depth (Lång et al. 2006). While this results in more shallow depths to salt/freshwater interfaces than traditional estimates might dictate (Fetter 2001), this estimate of the approximate height of freshwater column is set in order to account for the risk for brackish water infiltration in addition to relic salt water infiltration (Lång et al. 2006; SGU 2006). Based on these assumptions, available maximum storage volumes were calculated for each cell based on soil and rock type combined with elevation and soil depth (Figure 4(a)).

An assumed extraction for each well was assigned based on use. For a household well, a value of 300–495 l/day was applied (SCB 2010), with summer cottages at the lower end of this estimate. For recreational houses, extraction was limited to summer months but information was not available as to the water-use standards of individual households. The code presented here allows for modelling of permanent and seasonal residency through randomly assigning each extraction point to one of these categories (either seasonal or permanent) based on a uniform probability density function and assigned residency rate in order to model scenarios of increasing degree of permanent residency. During summer months (May through September), it was assumed that all wells would be in use, but during winter months only those wells with permanent residents would be. Due to the large numbers of wells in these areas, a well extraction radius was assumed, where an extraction rate was evenly distributed over an area corresponding to the defined radius. These radii for all wells were projected, and discrete polygons created from the overlapping sections. The centroids of these polygons were assigned an extraction rate per square metre based on the total number of the contributing overlapping polygon centroids multiplied by the rate per square metre. These values were then interpolated in order to create an extraction map (Figure 4(b)).

Water balance calculations from Östhammar show that the vast majority of the land area is unaffected by extraction from private wells (Figure 5). However, in more localized areas, volumes remaining within the modelled cells reached lower than 25% of the maximum storage of the cell, particularly in close proximity to the coast where available storage of freshwater would be very low due to low topographical values and shallow depth to saline water. Typical problem regions, or regions where remaining volumes were estimated to reach below 50% in 8-month and 50% permanent residency scenarios, occurred within roughly 500 m of the coastline (Figure 5). Often, these areas showed values which could be considered to be critical, having less than 25% of their maximum storage volumes. However, impacts of extraction were seen more than a kilometre inland in some areas, usually in areas with low elevation as they would have a more limited storage according to the model. Model results showed that pockets of soil in close proximity to areas of withdrawal were often largely unaffected, often having remaining reservoir ratios of more than 90% at the end of August. In the study areas, this is usually the end of the summer season immediately before the autumn groundwater recharge events begin. The extent of impacts of water extraction progress rapidly inland with increasing water use and residency, although even at low percentages of permanent residency and low rates of extraction reservoirs will dip below 50% in particularly vulnerable areas.

Increasing permanent residency rates within the model scenario had fairly evenly distributed effects. An increase to 50% from 15% residency leads to extensive decreases in reservoir volumes remaining after 8 months by 15–20% (Figure 6). More extreme impacts were found to be localized to areas which were already at risk, increasing the extent of these risk zones. Increasing permanent residency rates amplified the extent and intensity of problematic areas identified at lower rates. While the model results may appear quite similar, local comparisons of the model showed that differences between the models could be greater than 75% (Figure 6). As residency changes may often occur as conversions of seasonal domiciles (cottages) to permanent residences, additional environmental impacts due to increased water use will likely also be constrained to these localized areas.

The effect of subsurface flows within the model was also examined by direct comparison of two scenarios with identical parameters with the exception of medium (K = 1 × 10⁻⁸ m/s) or low (K = 5 × 10⁻¹¹ m/s) regional hydraulic conductivities, with a third baseline scenario with no subsurface flows, in addition to evaluating model stability (Figure 7). Average differences between the models were zero and −2%, with a modelled remaining storage standard deviations of 1 and 10%, for low and medium K scenarios, respectively. Model stability criteria were not fulfilled for higher estimates of hydraulic conductivity (2 × 10⁻⁶ m/s). Results indicate that subsurface recharge to depleted reservoirs has limited influence at great horizontal distances from the soils. The most substantial effects were seen in close proximity to pockets of soil or till, where groundwater reservoir levels were largely unaffected by the extraction whereas the surrounding bedrock were markedly affected. In these areas, reservoirs could be replenished by 20–50% and in some cases much more. However, the spatial extent of these changes was seldom greater than 50–100 m, where influenced areas typically showed differences around the perimeter of the original influenced zone. However, the area consisted mostly of hard rock outcrops and mean value of the Euclidean distance to soil cover in the study area was 502 m, although the distribution of distances appeared to be log-normal and most cells in very close proximity to the coast had soil pockets within roughly 100 m.

Extraction and storage maps themselves reveal the spatial heterogeneity of the hydrogeological parameters which were used in the groundwater balance calculations (Figure 4(a) and 4(b)). Modelled extraction volumes were largely concentrated in areas with the most residents. Extraction rates were seen to double or potentially triple over the span of a few hundred metres and ranged from 0 to 4 × 10⁻⁶ m³ per day per m². Similarly, when considering a 100 m² model resolution, estimated storage values can vary several orders of magnitude within a few hundred metres. Isolated pockets of till and even clay areas covering till may increase the storage capacity within a cell by more than 0.200 m³ from perhaps as little as 0.05 m³. Local pockets of soil likely act as a recharging reservoir when volumes are extracted from hard rock and have the additional function of increasing hydraulic head locally, resulting in a decreased likelihood of saltwater groundwater upconing or intrusion (Figure 8). The storage volumes within these soil pockets are likely largely unaffected by the extraction (Figure 8) as the storage volumes within the model cell far exceeds the quantity of the modelled extraction, due to the large differences in kinematic porosity between the local soils and bedrock.

Varying different parameters within the model in order to assess parametric sensitivity was also carried out. Comparisons were made between soil-depth models estimated from well archive and geological data (SGU 2017) and a soil-depth model estimated from the slope of outcropping bedrock (Karlsson et al. 2014). Differences between the storage models, using either the simplified regolith model (Karlsson et al. 2014) or SGU's soil depth model (SGU 2017), varied between null and more than an order of magnitude from the original potential storage (Figure 9(a)). Similarly, varying the kinematic porosity values for till (0.02 or 0.05) and clay (0.001 or 0.01) caused large variations from the original model (Figure 9(b)) as did altering the kinematic porosity of the bedrock (granites from 0.0005 to 0.001 and gneisses from 0.0003 to 0.0008) (Figure 9(c)). Finally, the thickness of the till layer was doubled (1–2 m) where relevant soils, e.g. clay or post-glacial sand, were present on the soil type map (Figure 9(d)). Doubling the thickness of the till layer had extremely localized effects since the area consisted mostly of exposed bedrock, whereas in the other scenarios, the storage characteristics were more evenly distributed throughout the study area. When changing the kinematic porosity and soil depth model, variations of 50–100% were common, whereas changing the kinematic porosity of the bedrock had similar effects which were of the magnitude 50–100% where no soil was present and 0–10% changes in areas which were soil covered. While no data are available to directly validate the storage maps, Spearman's correlation analyses with the different storativity maps showed significant correlations between storage and estimated specific capacity with more than 90% confidence, and between storage and groundwater level taken from the well archive at the time of drilling with more than 99% confidence. Distribution of the storage variable showed strong bimodal characteristics, regardless of how the storage variable was modified during the sensitivity analysis.

The results from the groundwater balance calculations indicate that spatial variance of storage in hard rock coastal terrains is often the determining factor governing sustainable water usage. Inclusion of a simple estimate of subsurface flows shows that flows within the fracture network recharged extracted volumes to a limited extent, but low hydraulic conductivity values of the hard rock and heterogeneous spatial distribution of soil cover imply that the influence of subsurface flows is limited. The emphasis of this groundwater balance model is an applied tool for groundwater resources management and should not be taken as a robust steady state or transient numerical flow model. Several important sources in the water budget have been neglected as their quality would be undesirable in terms of drinking water supplies, e.g. reintroduction of infiltrated water in septic systems or infiltration from saline water bodies. Municipalities in hard rock coastal regions, which often lack resources and/or hydrogeological knowledge with which to construct, calibrate and validate robust hydrogeological models will likely greatly benefit from surrogate tools such as this one which identify areas of likely environmental stress. Identification of these areas could allow for more resilient resource management strategies and/or a better allocation of limited resources to decrease uncertainty in vulnerable areas.

The model results also clearly illustrate the importance of a spatial approach in water resource management questions where soil cover is limited, seen in the spatial variability of the remaining reservoir ratio in 8-month scenarios (Figure 5), largely due to heterogeneity in the hydraulic parameters, in agreement with Batelaan & De Smedt (2007). As reservoir volumes on a regional scale are often largely unaffected by the extraction made in a few local housing areas and a limited number of private wells, if a spatial approach or a too-coarse model resolution is adopted the erroneous conclusion could be made that the water resources in a given region are capable of supporting far greater numbers of residents than would be possible in reality. Water resources in terrains where soil cover is limited must be considered at quite a local scale, as even after a few hundred metres resources which are present may not be obtainable due to low bedrock transmissivities. Exacerbating the problem is the observation of non-uniform spatial distribution of extraction. Contrary to watershed- and subwatershed-based mass-balance approaches where usage is assumed to be dispersed equally over the study area, houses will often be clustered within a specific and limited area. Thus, the extraction from the reservoirs in the centres of these areas will be greater than if residents were uniformly distributed over the area. Increasing permanent residency rates may not influence the regional water balance particularly heavily, but instead the effects may be constrained to local areas and may amplify stresses in areas which are already the cause for some concern. If water use is more evenly distributed, these effects might be avoided. In the model, reintroduction of untreated water (through, for example septic system leakage) and surface water flows are neglected. This former assumption is based solely on quality issues, where shorter residence times due to low-storage volumes indicate that this infiltrated water would likely be undesirable as a water source in problem areas. The latter assumption of neglecting surface water interactions in the model is based on two factors: (1) that the smaller streams in Sweden usually exhibit gaining behaviour where groundwater typically outflows to the stream rather than the opposite and (2) that streamflow in these types of terrain is often very limited, especially during summer months. Similar to methods such as the SWB-model (Dripps & Bradbury 2007), the spatial groundwater balance strives to balance sufficient conceptual descriptive power of the model with applicability, which is particularly useful in areas where relevant time series or calibration data are present, such as in Müller et al. (2020). This leads to what is perhaps the greatest benefit of this approach, that the model might be best applied as a tool for resource management. Using this approach, local and municipal environmental engineers can identify problem areas and allocate resources to mitigate these problems. Applying more simple methods (subwatershed-level balances) might misrepresent the actual storage and extraction behaviours in the region, and more complex numerical methods (e.g. 3D Finite Element Modelling) might require extensive groundwater level measurements, experience and skill which many coastal municipalities do not have or lack resources to acquire. This underscores the importance of developing applied, heuristic solutions such as the method presented here to increase water security in areas which lack sufficient resources and/or hydraulic data coverage for more traditional hydrogeological solutions.

Local till- and clay-filled valleys, which would likely often be overlooked in regional groundwater balances, have vital roles in mitigating environmental stresses which may occur as a result of overextraction from private wells. The importance of such geological features in groundwater resources modelling has been established in earlier research such as Dripps & Bradbury (2007) and Batelaan & De Smedt (2007). While perhaps not contributing much to the water balance in the form of groundwater flow in the model, in reality should one of these valleys be located between an extraction zone and a saline water body or contamination source they could act as a groundwater divide, meaning that saline water would be less likely to flow into the extraction zone. This mitigating role would be emphasized even more during dry years, when recharging volumes are even more limited, and would essentially serve as a protective barrier (Figure 8). However, it is still possible that wells in these extraction areas would be hydraulically linked through the fracture network to saline water sources and thus be at risk. The extent of influence occurring from these local reservoirs appears to be limited to roughly 50–100 m over the course of a summer season. However, while the average Euclidean distance from a single cell in Östhammar to a soil-filled valley was more than 500 m, most of the cells in the study area had soil-filled cells within 100 m or less. Additionally, as the extraction was simplified to an area of influence around each well, groundwater flows very near the wells themselves is under- or misrepresented. Modelled flows assuming an equivalent porous medium with a uniform hydraulic conductivity over the study area do not entirely represent the true dynamic of the flow regime within the fracture network, where the majority of flows occur within a very small percentage of the fracture network and the hydraulic conductivity field would behave in a far more stochastic way as model cells are nowhere near large enough to approximate representative equivalent volume. However, as topographical depressions in the study area may often coincide with increases in fracture frequency, it is not unreasonable to assume a significant groundwater exchange occurs between the fracture network and soil-filled valleys. As the intended use of the model is to be a decision support tool, the simplified flow model was intended not to be a robust approximation of the groundwater transport network, but instead to account for some of the exchange between the soil and fracture reservoirs.

Characterizing and analysing parametric uncertainty in groundwater storage is of vital importance in reservoirs where kinematic porosity values are limited and/or they are in close proximity to saline water sources. Assessing the sensitivity of the storage parameter by varying the soil depth, till thickness or kinematic porosity of the different subsurface materials (Figure 9) showed effects which varied both spatially and in magnitude. As the study area is mostly covered with till and clay where soil is present, it is logical that even a small change in these parameters would have significant impacts. However, kinematic porosity of the bedrock was shown to have a limited influence in areas which are covered by soil, implying that there is a hierarchical order of importance of reservoirs. Where pockets of soil are more frequent, fewer resources need be allocated to minimizing uncertainty in the fracture network, as the role of the bedrock in storage is limited. However, in areas with more sparse soil cover, uncertainty in this parameter will increase model uncertainty proportionally. Thickness of the till layer had a very local influence on the modelled groundwater storage sensitivity, primarily being centralized in areas with relatively thick soil cover. This may be due to the limited thickness of the soil depth in the area, as much of the soil cover may be less than the modelled change in thickness (from 1 to 2 m). Implications of the sensitivity analyses are that results may prove vital in identifying areas which are of interest for exploitation. If year-round occupancy of a particular region is forecasted to increase, applying the model will allow for a more targeted use of limited resources in decreasing uncertainty in strategic areas, such as more detailed local flow models or groundwater balances. For example, electrical resistivity tomography, depth sounding or drilling could be applied to a few strategic locations (for example, as identified in Figure 9(d)) to assess actual till thickness and drastically lower uncertainty in the model in order to better manage local water resources. In other areas which show lower model sensitivity, field investigations may not be necessary for providing an adequate conceptual basis for managing the local groundwater resources. Applying the explicit solution of the finite difference method to estimate flow (Equations (2) and (3)) in the model demonstrated that the model is sensitive to groundwater flows locally, but the extent of influence of the flows is rather limited, particularly at lower hydraulic conductivities, implying the application of the tool to identify areas of environmental stress due to excessive groundwater abstraction is viable.

Limitations of the model and sources of uncertainty should be addressed. Actual recharge volumes, water usage, hydraulic conductivity and kinematic porosity are all uncertain parameters used in this model. Recharging volumes of water were estimated as the difference between effective precipitation and evapotranspiration, neglecting slope, type of soil and temperature. This approach to recharge was utilized in order to maintain model simplicity and the assumption that while runoff would be a significant factor in areas with frequent hard rock outcrops, much of this runoff would infiltrate through the till at the peripheral of the outcrops of soil-filled valleys regardless. Limiting the storage of each model cell, while vital in this type of terrain, is also a source of uncertainty. In the model, volumes of water in excess of the storage capacity of the model cell were lost. In reality, these volumes would flow overland and perhaps be lost as runoff to the Baltic Sea but potentially could contribute to recharge in another area. Treatment of bedrock as a homogeneous medium also contributes to model uncertainty in both extraction maps and subsurface flow calculations. Extraction would not be distributed evenly over a circular area but would be dependent on the dominant fracture orientations and distribution, as well as having a strongly in situ radius of influence. Thus, radius of influence would be spatially anisotropic and its extent largely heteroscedastic. Sensitivity analyses of the impact of extraction radius have been shown to influence the extent of the areal impact of groundwater use but not identification of impact areas (Earon 2019). Similarly, kinematic porosity of the bedrock would also be dependent on the nature fracture network. Water reserves within the pockets of soil in the model were less affected by extraction than in the rock mass, and so the most important processes to capture in these settings would be the flow from soil reservoirs to the hard rock. The intention of this model is to provide a more holistic understanding of the local effects of water use in hard rock coastal terrains, and thus a rough approximation of subsurface flows was deemed appropriate.

Finally, there are some variables which would influence the model results but which have not been included. While incorporating the simple flow model into the groundwater balance might address some of the natural variations in groundwater levels, it is unlikely that the simplicity of this estimate can accurately represent the complex behaviour of the flow systems in the fracture network. A result of this might be that the actual volumes of water which are sustainably extractable are overrepresented in the model. However, given the low-average hydraulic conductivity in the bedrock, large deviations from the model due to local anomalies in the hydraulic conductivity field may be few. It might be possible to incorporate a general statistics-based groundwater decline for each cell based on further analyses of climate and groundwater level time series combined with local vegetation, soil or rock type. Septic or infiltration beds for private waste water would influence the groundwater levels at a local scale and mitigate the actual volumes of water which were extracted, as these volumes are largely not transported from the site of extraction. Additionally, infiltration of water from saline water bodies was also not included and would likely occur due to the constant hydraulic head of the surrounding sea. However, it should be mentioned that this model addresses only freshwater reserves. While the overall groundwater levels will not be as low as those which were modelled, the water quality situation will likely be significantly worsened in the areas with greatest impact. Volumes of water which were shown to be inflowing from cells surrounding the area of influence of extraction would be less (due to a lower hydraulic gradient), meaning that freshwater is not replacing the water which was being extracted, but rather grey water or worse is being introduced. Due to the low porosity of the bedrock and limited soil depths, this in turn leads to local reservoirs which are highly vulnerable.

In heterogeneous terrains with low inherent storage and non-uniformly distributed water use, a spatial approach to groundwater resources inventory and management is needed to avoid excessive overextraction and/or potential water quality problems which may arise (e.g. from seawater intrusion). Spatial heterogeneity in extraction and subsurface storage characteristics result in local variations from regional averages, which may have negative environmental impacts but also may be missed if a model resolution is used which is too coarse. Variability in groundwater reservoir levels at a 100 m² resolution implies that vital information will be missed should regional groundwater balance estimates be applied rather than incorporating a higher spatial resolution, although more research is needed regarding the effects of spatial resolution on model results in different geological settings.

Incorporating limited storage within groundwater balance models and heuristic solutions for groundwater management has vital implications for resilient groundwater resources management in terrains with limited soil cover. Storage in terrains with such limited hydraulic conductivity values has an important role in this setting in governing the volumes of water which will be sustainably extractable by residents who are not reasonably able to connect to municipal water supplies. Given the heterogeneous spatial distribution of extraction, rock type and soil cover, it is therefore essential that groundwater balance methods be adapted with a limited storage component to compensate for these conditions where possible.

The approach presented in this study shows promise as a tool for successful management of groundwater resources in hard rock coastal terrains with limited soil cover, particularly in regions where municipal planners may have limited resources to apply. Results of this study imply that identification of extraction rates and storage volumes on a local scale can be accomplished using available geological, topographical and land-use GIS databases, which could lower the demand for expensive hydraulic tests and computational models which require extensive amounts of input data. This will lead to lower time and resource costs needed to identify and mitigate areas which may be susceptible to environmental stress and may have sparse hydraulic measurement data. The use of existing datasets which are available for most of Sweden and obtainable in many other countries in digital form is an enticing aspect of this approach to groundwater modelling. This allows the method presented in this study to be readily applied to all regions in Sweden and with little modification and to be applied to any region with general physical and geological information where storage in the subsurface may be the limiting factor in groundwater balance questions, such as in large parts of Canada or Scandinavia.

The authors would like to acknowledge the contributions of the anonymous reviewers for comments that have improved the quality of this study. Financial aid for this project was provided by the Geological Survey of Sweden (Dnr. 60-1640/2007) and the Swedish Research Council for the Environment, Agricultural Sciences and Spatial Planning (Dnr. 2017-01504). Digital geological data were provided by the Geological Survey of Sweden and digital topographical data were provided by the Swedish Land Survey.

Data cannot be made publicly available; readers should contact the corresponding author for details.