Contribution of bedrock structures to the bedrock surface topography and groundwater flow systems within deep glaciofluvial aquifers in Kurikka, Western Finland

One of the key challenges related to glaciofluvial aquifers is understanding how the underlying bedrock structures and the associated bedrock surface topography contributed to the deposition of the glaciofluvial sediments and the generation of groundwater flow pathways. In this study from Western Finland, we present a new digital elevation model of the Precambrian crystalline bedrock surface buried up to 100 m beneath the glaciofluvial sediments along the Kurikka depression. We link bedrock topography to structural anisotropy using additional data from rock outcrops, boreholes and gravity data. Brittle deformation zones are seen in the bedrock‐digital elevation model (DEM) as sharp breaks in the rock surface. These vertical fracture zones contributed to the development of interlinked elongate bedrock depressions and delimit rock blocks with different elevations. A narrow, zigzagging trench, with a stepped floor, follows WNW‐ESE and NE‐SW oriented fracture zones and forms a major hydraulic connection between two major parts of the aquifer. Furthermore, the conductive fracture zones are conduits which connect the shallow glacifluvial aquifer system to deep groundwater in the bedrock. Understanding bedrock structures and buried topography is critical for successful groundwater modelling in crystalline, fractured bedrock.


| INTRODUCTION
Groundwater flow within glacial drift aquifers is known to depend not only upon the thickness, physical properties and internal structure of the Quaternary sediments but also on the topography of the underlying bedrock (e.g. Godbout et al., 2011) and the hydraulic connectivity between the bedrock and the overlying glaciofluvial deposits (Cassidy et al., 2014;Kempton et al., 1991;Lipponen, 2006). Despite the recognised role of the bedrock, systematic analysis of the bedrock structures as a base layer for the glacial deposits has been frequently overlooked in regional aquifer modelling, with exceptions, for example, on the Canadian Shield (DesRoches et al., 2014;Gao, 2011;Gleeson & Novakowski, 2009;Richard & Chesnaux, 2014;Steelman et al., 2018). For this reason, digital elevation models representing the bedrock surface (hereafter bedrock-DEM) and acting as the depositional basement for the glaciofluvial aquifers do not always adequately honour the bedrock structure. Moreover, most existing geological and hydrological models lack sufficient understanding of the linkage between the crystalline basement and the sediment cover.
Bedrock structures contribute to the bedrock-DEM morphology, but brittle deformation zones within the crystalline bedrock also provide groundwater flow pathways which largely define the hydrogeological conductivity of the bedrock. Distributed flow along individual planar discontinuities may be localised into channelised pathways (Brown et al., 1998;Sanderson & Peacock, 2019;Stephens et al., 2015) generated by interlinked linear 'tubes' of pronounced conductivity within fault bends, intersections and step-overs (Cox, 1999;Dimmen et al., 2017;Micklethwaite et al., 2015;Peacock et al., 2018).
In this paper, we show that the bedrock structure and surface morphology play a much more fundamental role in governing the groundwater flow in large-scale glaciofluvial aquifers than previously considered. Specifically, zones of brittle deformation act as high permeability conduits, as well as contribute to the generation of bedrock depressions, which further localise groundwater flow into channelised networks. Further, we show that even limited datasets of the bedrock structure available from poorly outcropping areas are useful in addressing the morphology of the bedrock surface topography.
We conducted the investigation within the Kurikka buried valley aquifer complex, Western Finland ( Figure 1). The aquifer complex is characterised by laterally extensive coarse-grained sediment assemblages comprising alternating layers of till, with significant groundwater resources comparable with the buried valley aquifers in Canada (Ross et al., 2005;Russell et al., 2004). At large scales, deposition of the subglacial sediments is controlled by the ice sheet thickness, glacier bed and ice sheet topography (Benn & Evans, 2010;Boulton, 1986;Boulton et al., 2001;Lovell et al., 2011;Shreve, 1972Shreve, , 1985 which govern both the localisation and the stratigraphy of the deposits. Furthermore, spatial association between individual glaciofluvial deposits and bedrock structures such as brittle fractures and lithological contacts and particularly the morphologic features they induced upon the bedrock surface have been recognised (Barr & Lovell, 2014;Burbank & Fort, 1985;Okko, 1962;Palmu, 1999;Skyttä et al., 2015;Sugden & John, 1976). Brittle structures have a direct linkage to the hydrogeology by the groundwater filling the pore spaces of the sediment and rock. Those voids play the major role in this type of aquifer system. Recent (Putkinen et al., 2020;Putkinen et al., 2012;Rashid et al., 2020) and ongoing studies by the Geologi- and hydrogeological properties of the Quaternary deposits alone cannot explain the large amount of groundwater recharge and that some of the bedrock depressions with fractures likely provide hydrogeological linkages between different parts of the aquifer.
However, the bedrock structure has not been accounted for when interpolating the bedrock-DEM from dominantly gravity data, a typical dataset used in bedrock surface modelling (Møller et al., 2007;Steelman et al., 2018), and no further investigations on the bedrock structure have been conducted.
In this study, we addressed the relationship between the bedrock structure and bedrock surface topography through conducting structural analysis of the bedrock within the main recharge area of the Kurikka aquifer system (Figure 1c). We used the results to reinterpret the morphology of the existing bedrock-DEM which was previously generated without utilizing the structural information of the bedrock.
We placed specific attention upon recognition of steeply dipping zones within the bedrock-DEM which likely reflect the presence of discrete vertical to sub-vertical deformation zones but where the surface geometries are not realistically represented due to the characteristic smooth, curvilinear shape of the bedrock surface interpretation derived from gravity data. Subsequently, we re-interpreted the gravity data from selected survey lines and used the revised interpretation together with the field observation data in modifying the bedrock-DEM to better honour the character of the structure of the bedrock.
We also mapped bedrock lineaments from ground surface DEM for the regional scale comparison. Moreover, we generated a 3Dphotogrammetric model from a bedrock outcrop in the study area to illustrate the effect of a discrete shear zone and associated fracturing on the topography of the bedrock surface and compare the results with the improved bedrock-DEM to reduce the conceptual uncertainty of the models.

| GEOLOGY OF THE STUDY AREA
The western and southern parts of the study area (Figure 1a,b) comprise post-kinematic, largely undeformed 1.88-1.87 Ga granodioriticgranitic-dioritic intrusive rocks of the Central Finland Granitoid Complex (CFGC) (Nikkilä et al., 2016;Nironen et al., 2000). The northern and eastern parts of the study area comprise supracrustal 1.90-1.89 Ga rocks of the Ostrobothnia Schist Belt (OSB). OSB is lithologically dominated by biotite paragneisses and structurally characterised by the ENE-WSW striking, dextral Vittinki high-strain zone (Pääkkönen, 1966) occurring approximately 30 km north of the study area and NE-SW striking ductile shear zones abutting the Vittinki zone at the central parts from the Southwest (Figure 1b; Vaarma & Kähkönen, 1994;Lehtonen et al., 2005). Structural signatures within the study area and its surroundings have been attributed to the coupled Bothnian oroclines (Lahtinen et al. (2017) bounded by the migmatitic Vaasa Complex (VC; Chopin et al., 2020) in the west and the CFGC in the east (Figure 1a).
The bedrock of the study area is characterised by negligible variability in overall topography and ultra-slow erosion of <2.5 m/Ma as revealed by tiered unconformities and impact structures (Hall et al., 2021). However, the presence of Mesoproterozoic sedimentary rocks comparable with the sandstones in the Bothnian Sea and Satakunta half-graben (Pokki et al., 2013) is found within the 100-m-deep valley systems in Ostrobothnia and point towards the presence of long-lived valleys since 1.5 Ga (Hall et al., 2021). In Kurikka area, the bedrock surface morphology is characterised by one major elongate NNE-trending depression located within the Kyrönjoki river valley (KY in Figure 1c) and minor but distinct depressions beneath the approximately NNW-SSE-trending Paloluoma (PA) and Häjyluoma (HÄ) Valleys. Since their formation, a long burial history persisted until the valleys were finally exhumed at about 40 Ma ago (Gibbard & Lewin, 2016). Since then, the evolution is controlled by the river system operation across the Late-Paleocene and Neogene periods while the Pleistocene glaciers made the final touch-up for the first-and second-order valley morphologies. Besides the larger depressions, a connective E-W-trending depression at Kuusistonloukko (KU) occurs in-between the HÄ and PA depressions and has been considered to provide hydrogeological linkage of aquifers associated with the HÄ and PA depressions (Putkinen et al., 2012). However, the character of the depression is so far unknown. Even smaller valleys are present in the Lehtivuori-Karhuvuori areas (third-order valleys).
However, Pitkäranta (2009)  Infiltrated groundwater recharge routes from the topographically higher areas (in the west) to the aquifers are controlled by the atmospheric pressure and gravity, and the ground elevation difference of over 100 m from the main recharge area to the destination of aquifer areas constitutes up to 100-m difference to the potentiometric surface for the groundwater that has artesian character in sediments and bedrock. In the KU depression, Putkinen et al. (2012) attributed the hydrogeological connectivity of the different aquifer parts to the presence of the conductive middle sand horizon. Putkinen et al. (2012) also recognised the role of the bedrock fracture zones as potential carriers of groundwater, and recent finite element groundwater flow modelling (Rashid et al., 2020) has further confirmed that flow systems within the Kurikka aquifer system are complex, and the bedrock and its fracture zones apparently have a major role in controlling the inflow and maintaining the groundwater system ( Figure 2d).

| DATA AND METHODS
Understanding the brittle bedrock structures in areas of limited outcrop must be primarily developed using conceptual models over the regional fracture patterns, with particular focus on the zones of localised brittle deformation (faults, fracture zones) comprising intensely fractured cores and surrounding damage zones (DZs) (Childs et al., 2009;Kim et al., 2004;Peacock et al., 2017). Further, lithological composition of the crust, structural inheritance from older (ductile) structures and reactivation of older deformation zones will likely contribute to the generation of brittle structures, including the architecture of the fault zones (Korme et al., 2004;Mancktelow & Pennacchioni, 2005;Pizzi & Galadini, 2009;Skyttä et al., 2021;Skyttä & Torvela, 2018), and should hence be studied to decrease the uncertainties related to the build-up of the brittle structural model.
Bearing in mind the above, we conducted systematic analysis of the ductile and brittle structures of the bedrock within the study area using regional scale geological and geophysical datasets (aeromagnetic map, 1:200 000 Bedrock map form the Geological Survey of Finland; LiDAR elevation data and topographic map from National Land Survey) aiming at defining the structural framework of the area comprising (i) the major discontinuities, representing the deformation zones, and (ii) ductile trends inferred to represent foliation and lithological contacts. We verified and updated the regional interpretation with field observation data over the lithology and bedrock structure from 274 outcrops. The structural field data comprise observations over the variations in apparent foliation intensities, measured orientations of the ductile (foliation, fold, mineral lineation) and brittle (fractures and outcrop walls) structural elements. Furthermore, we observed the topological relationship of the dominant fracture sets, as well as how they control the topography of the exposed bedrock surface on highly elevated, well-outcropping areas. We defined fracture densities on outcrops and categorised the areas into seven classes (0-6) based on the number of observed fractures proportional to a reference outcrop size of 5 Â 5 m to highlight the relative changes of the fracture intensity within the study area. We studied the relationship between the brittle and ductile structures to recognise any potential structural inheritances and used the resulting understanding to provide guidelines for probable fracture orientations within the areas where no outcrops are present but which still need to be covered by the bedrock-DEM modelling.
To reduce the conceptual uncertainty over the relationship For the final part of the work, development of a structurally constrained bedrock-DEM, we used the dataset described in Putkinen et al. (2012), updated with new data (gravity, refraction seismic studies, groundwater monitor data) and models (bedrock-DEM) until December 2020 as a result from ongoing investigations by GTK. The gravity data comprise 71 survey lines, with a total length of approximately 68.5 km (20-m measurement point spacing), 32 groundwater well drilling sites, 2 intersecting seismic lines trending NE-SW and E-W and a bedrock surface DEM generated from these data using Topo to Raster and IDW interpolation methods in ArcGIS and MOVE, respectively.
Gravity data have proven successful to assess the topography of bedrock covered by sediment deposits especially for buried valley type formations (Ahokangas et al., 2020;Valjus et al., 2004), given that the density of the bedrock and the soil cover and the depth to bedrock at the end of the gravity survey lines are known. For reinterpretation of gravity data, we selected five crossing survey lines located within the area of intersection of the NNW-SSE-and E-Wtrending PA and KU bedrock depressions (Figures 1 and 10). Reinterpretation involved promoting the presence of abrupt vertical breaks in the gravity data as initially recognised from over 10-m elevation contrasts between neighbouring interpretation points. For the selected lines, we conducted a 2.5-D interpretation separately for each line, using densities of 1.9 and 2.67 g/cm 3 for the saturated sediment layer and the underlying crystalline bedrock, respectively. Reinterpretation used the same parameters (densities and the drilling depth to bedrock) as the original model but additionally used the sharp sub-vertical breaks as input constraints.
We integrated the new gravity interpretation lines into an improved detail-scale bedrock-DEM using the orientation of the dom-   (Figures 7 and 11). These re-interpreted profiles provide the baseline for improving the bedrock-DEMs (Section 5.3).

| The 3D-photogrammetry model
The outcrop used for the detail-scale 3D-photgrammetric modelling is characterised by a narrow E-W-trending, north-dipping shear zone

| INTERPRETATION AND DISCUSSION
Based on the presented results (Section 4), it is evident that the initial bedrock-DEM with smooth elevation transitions from the elevated areas to the elongate depressions ( Figure 1d) does not realistically represent the true fracture-controlled character of the bedrock surface topography (Figure 8). For this reason, in the following sections, we place the observations of the brittle structures in a regional context and further use as orientation constraints for re-interpreting the gravity data and improving the bedrock-DEM, so that the new DEM better honours the structural heterogeneity of the bedrock. The improved detail-scale bedrock-DEM targets a limited area within the critical intersection of E-W and NNW-SSE bedrock depressions, but we discuss the applicability of the gained results in the scale covering the available bedrock depth data.
5.1 | The relationship between the valleys, bedrock surface depressions, regional structures and observed fracture patterns The bedrock surface topography within the study area is uneven and characterised by several elongate bedrock depressions. Considering The structural characterisation of the bedrock allows us to make a regional interpretation of the major bedrock deformation zones in the area and shed further light into the origin of the studied valleys.
Considering the bedrock surface morphology, the NNE-trending major depression that spatially coincides with the KY river valley also coincides with the SZ1 deformation zone and the distinct NNE ductile structural trends of the bedrock (Figure 9). The smaller NNW-trending PA and HÄ depressions terminate against the KY depression in SSE and are recognisable 4-6 km towards NNW. The PA bedrock surface depression is a well-defined narrow (500-700 m) linear feature, which  Figure 9). On the above basis, we categorised bedrock depression into major (first-order), minor (second-order) and connective (thirdorder) groups by their size, rate of matureness and prominent orientation trend and topological relationships (Figure 9), to highlight morphological differences of the bedrock depressions as successfully conducted by, for example, Hall and Gillespie (2017).  (Figure 6), but they typically abut against the longer N-Strending ones and hence play a smaller role in governing the bedrock surface morphology. The areas with the highest relative fracture intensities occur to the east of the KY zone and at the intersection area of F I G U R E 7 Initial bedrock-digital elevation model (DEM) (Topo to Raster & IDW interpolation). Gravity data lines are black solid lines (bolded: re-interpreted gravity lines in Figure 10). White points represent locations where are 10-m elevation contrast between the neighbouring interpretation points. The rectangular area represents the scale of the improved bedrock-DEM in Figure 11. [Color figure can be viewed at wileyonlinelibrary.com] the KU, HÄ and PA zones (Figure 9b). Such spatial distribution supports that the bedrock depressions are expressions of tectonic faults or zones of brittle deformation, with surrounding DZs characterised by increased fracturing (Choi et al., 2016;Torabi et al., 2018) and heterogeneity of fracture densities (as caused by the presence of secondary slip surfaces; Ostermeijer et al., 2020). Moreover, the apparent lack of fractures may be indicative of heterogeneous DZs and the mechanical anisotropy of fault-stage structures and its controlling effect on the subsequent regional fractures (Skyttä et al., 2021).
Direct observations about fractured rocks representing zones of brittle deformation were found from two boreholes (located at the ends of the gravity line L57; Figure 11b,c), but as no drill core or drill hole video was acquired, we have no further control on the character of fracturing. However, using the recognised intimate relationship between the bedrock surface topography and the fracture system (Figure 8), we argue that the fracture density within the unexposed cores of the deformation zones underlying the bedrock depressions is much higher than within the surrounding bedrock.
Consequently, the localised character of fracturing contributed not only to the bedrock surface topography but also to the hydraulic transmissivity of the bedrock, which is pronounced within the brittle deformation zones.

| Re-interpretation of selected gravity survey lines
We re-interpreted gravity survey lines L57, L59, L61 and L63, which are located within the intersection of the WNW-and NNE-trending KU and PA bedrock depressions (Figures 1d and 9). Reinterpretation of the gravity data involved replacement of steeply dipping bedrock-DEM segments in the initial bedrock surface, resulting from the >10 m elevation contrast between the neighboring gravity interpretation points, with distinct sub-vertical segments.
These segments correspond to the fracture-controlled abrupt subvertical breaks of the bedrock surface (Figure 8), and dip constraints (80-90 dips) for these breaks were defined by the dominant fracture sets, represented in Figure 6, occurring in the vicinity (Figure 9).
The dominant sub-vertical fracture attitudes within the area covered by the re-interpreted gravity lines and the bedrock-DEM (Figures 10   and 11) are supported by the change from inclined to subvertical fractures across the Iso Karhuvuori-Pikku Karhuvuori transect ( Figure 6b). We infer that the ground surface depression is underlain by a bedrock surface depression which is defined by a zone of brittle bedrock deformation, and this deformation zone also acts as a boundary separating bedrock domains with contrasting fracture orientations.
Prior to the re-interpretation of the gravity data, just five such >10-m elevation contrast were recognised from lines L57, L59, L61 and L63 (all from L59; Figure 7), which were interpreted with the conventional approach placing no focus on potential structural discontinuities. Implementing the new approach, the re-interpretation introduced nine additional sub-vertical breaks that define either isolated thresholds but more commonly graben-like structures where bedrock surface locates 10-25 m lower than within the surrounding areas ( Figure 10). Introducing the sub-vertical breaks but still honoring the original gravity data points, the intervening segments between the breaks became more horizontal and the overall topography more fragmented and discontinuous with respect to the initial, smoothly varying geometry. Moreover, the revised interpretation generally improved the quality of the interpretation on the selected gravimetric lines or at least resulted in equal goodness-of-fit as the conventional interpretation, as shown by the match of the interpretation curve with the Bouguer anomaly curve (Figure 9). This indicates that the new bedrock models with the distinct fracture-controlled vertical breaks, comparable with the geometry within the 3D-photogrammetry model in Figure 8, are at least as plausible as the smooth ones. Models based on gravity data, like the potential field models in general (e.g. Møller et al., 2007), quickly turn smooth in appearance if a critical depth is exceeded, but this was not the case in the present study. F I G U R E 1 0 Re-interpretation of gravity data along survey lines L57, L59, L61 and L63. The fracture-controlled breaks within the bedrock-digital elevation model (DEM) marked are as purple solid lines, and the numbers indicate their approximated dips as inferred from re-interpreted bedrock surface (for location, see Figure 7). Purple dotted lines on L59 are breaks recognised before re-interpretation; equivalent with white dots in Figure 7). We compiled a structural interpretation developing the isolated bedrock breaks (Figures 7,10) into a network of bedrock discontinuities (fracture systems, faults) by correlating the breaks with fracture data (Section 4.2) and distinct linear slopes on the original bedrock-DEM (Figures 9 and 11b). The most pronounced slopes within the original bedrock-DEM within the whole study area trend approximately N-S (Figures 9 and 11a), whereas the detail-area additionally comprises distinct NE and WNW-trending slopes, which we consequently used to define the trends of the block boundaries. The resulting network of brittle deformation zones hence suggests that the KU zone linking the HÄ and PA depressions is not a linear feature but rather composed of approximately orthogonal, NW-and WNWtrending segments of connected bedrock structures (Figure 11b).
Using the above sub-division of bedrock blocks, we modelled each block separately using the elevation point data available for the The bedrock depression defined by the elongate Block 1 is asymmetric, with the lowest elevations along the southern margin of the block and progressive increase in elevation towards the northern block margin (Figure 11c,d). This is comparable with the asymmetric fracturecontrolled bedrock depression imaged by the photogrammetric methods ( Figure 8) and supports the used approach in improving the bedrock-DEM. Our focus in this first-pass conceptual bedrock-DEM reinterpretation was placed on recognising and modelling the domains with the most pronounced elevation contrasts, that is, the 'highs and lows'. However, elevation contrasts are present also within the lowelevation blocks, which indicates that significant but smaller-scale fragmentation occurred within the bedrock depressions (Figure 11c).

| Applicability of the new approach to regional bedrock-DEM signatures
The main problems with the initial bedrock-DEM from the study area are associated with the standard interpolation approaches, which do not adequately address (i) the structural anisotropy of the of bedrock and its spatial variation and (ii) the distinct sub-vertical breaks controlled by the brittle deformation zones or larger fractures. Instead, the interpolation results in smooth and rounded shapes within the bedrock-DEM in-between the input data points, not honoring the anisotropy of controlling bedrock structures. Examples of significance for the PA-KU area include two areas in Figure 12; the area of elevated topography within the KU depression fails to represent the suggested WNW-ESE-and NE-SW-trending sub-valley system (circle 1) and the smooth rounded shapes occurring on the eastern edge of the PA depression (circle 2; Figure 12a), where the edge does not honor the sharp geometry of the bedrock topography controlled by spatial anisotropy of fracturing.
We tackled the above issues within the initial bedrock-DEM by constructing the improved detail-scale bedrock-DEM which comprises both the sub-division of the crust into structurally controlled blocks, as well as honors the controlling structural anisotropy in each block (Figure 11c). For the regional dataset, we tested the applicability of using specific structural anisotropies within interpolation of the larger bedrock-DEMs. The first interpolation These examples also indicate that no single anisotropy value in surface interpolation algorithms will result geologically realistic bedrock-DEMs, but a new approach involving a dataset and an algorithm accounting also for the spatial variability in anisotropy will be needed. Moreover, the role of the major bedrock discontinuities and their contribution to the heterogeneity of the bedrock-DEM should be assessed in the studies of glaciofluvial aquifers as generating structurally controlled bedrock-DEMs may result in such local variations in the morphology of the bedrock-DEM that has major significance to the groundwater flow pathways.

| Implications of the new approach to the hydrogeology of the HÄ and PA buried valleys
Considering the flow pathways of groundwater, the revised bedrock-DEM presented in Figure 11c suggests that flow is more localised into structurally controlled channels of lower elevation along the bedrock surface than previously considered. As such, the mechanism is comparable with channelised fluid flow controlled by the geometry of individual structural surfaces (Brown et al., 1998). The connected pathway, 'channel', from HÄ to PA depressions is particularly important as it defines a suitable site of deposition for the well-conductive glaciofluvial sediments (middle sand) considered responsible for the hydrogeological connectivity between these two depressions (Putkinen et al., 2012). Moreover, distinct bedrock depressions comparable with, for example, the deepest parts of the asymmetric depression in Figure 8 may host even coarser sediments with higher permeability and further improve the channelizing effect.
The recognised fracture signatures of the bedrock, associated with the elongate bedrock depressions and abrupt breaks along the bedrock-DEM (Figures 8, 9 and 11), suggest that not only the bedrock depressions but also the fracture zones within the underlying bedrock act as fluid pathways (Faulkner et al., 2010). The linkage of the Häjynluoma, PA and KU deformation zones and the resulting network of fracture zones likely contributed to the groundwater flow within the bedrock, with potential development of channelised flow pathways also within the bedrock (Cox, 1999;Dimmen et al., 2017;Micklethwaite et al., 2015;Peacock et al., 2017;Peacock et al., 2018;Seebeck et al., 2014). Previous studies have shown that the hydrogeology system of the HÄ-PA buried valleys is strongly controlled by the watershed-scale bedrock surface topography and the local bedrock topography and the resulting groundwater channelisation (Figure 11c) and also by the hydraulic properties of the bedrock and the sediments (Putkinen et al., 2020;Putkinen et al., 2012;Rashid et al., 2020). We infer that the groundwater flow route from the recharge region downward to the buried valleys selectively utilises the material of the highest conductivity-be it the sediment cover or the fractured bedrock within the contact zone of the bedrock-sediment contact.

| CONCLUSIONS
Weathering and erosion over the last 1.5 Ga have exploited bedrock deformation zones to form valleys at Kurikka that today hold glaciofluvial sediments up to 100 m thick, hosting major glaciofluvial aquifers. Re-interpretation of the gravity data used for imaging and modelling bedrock surface topography shows that the margins of the elongated bedrock depressions are sharp in geometry, reflecting their fracture-controlled character. Mapping structural anisotropy in bedrock exposures and gravity data allowed the development of improved DEMs of the bedrock surface. Vertical fracture zones connect the surface aquifer to deep groundwater. Furthermore, previously unrecognised pathways for water movement were recognised along trench floors. Modelling groundwater movement in crystalline