3D geological modelling of a complex buried-valley network delineated from borehole and AEM data
Buried and open tunnel valleys are common features in former glaciated regions (Kehew et al., 2012, van der Vegt et al., 2012). The buried valleys are often a key factor for groundwater recharge and flow and in areas with buried valleys, delineation of the valleys and their infill is therefore important in relation to groundwater modelling (Andersen et al., 2013, Sandersen and Jørgensen, 2003). In places where the buried valleys are eroded into an impermeable substratum, the coarser infill of the valleys typically constitute important aquifers that tend to be deep seated and well protected. However, buried valleys can also be filled with fine-grained deposits with low hydraulic conductivity, in which case they can act as barriers for the groundwater flow. The buried valleys have typically been influenced by repeated erosional incisions and depositions, and the resulting valley infills are therefore often highly heterogeneous in nature. This heterogeneity affects the groundwater flow paths and may result in potentially scattered recharge areas with greater risk of downward transport of contaminants.
The buried valleys in Denmark are typically between 50–350 m deep and 0.5–4 km wide (Jørgensen and Sandersen, 2006). The individual valleys show characteristic shapes, such as undulating thalwegs and abrupt terminations. The valleys were formed by erosion of the substrata by highly pressurised meltwater flow underneath the ice sheet (Jørgensen and Sandersen, 2006). Generally, in areas with impermeable substrata, the subglacial meltwater cannot drain as groundwater and channels were consequently eroded deeply into the underlying sediments (Sandersen and Jørgensen, 2012). The same process has been observed in other places, as for instance, in Poland (Piotrowski et al., 2009). When tunnel valleys have been eroded and subsequently filled with coarser sediments, the hydraulic properties of the substrata were changed and the tunnel valleys eroded by the subglacial meltwater from one ice advance tend to be re-used by the meltwater of following ice advances (Jørgensen and Sandersen, 2006, Sandersen and Jørgensen, 2012). The result is complex networks of buried valleys with very heterogeneous infill caused by variations in both sedimentation history and history of erosional events. Especially the formation of cross-cutting valleys adds a high degree of complexity to the entire architecture of the valley systems.
In order to map the buried valleys, dense data coverage is essential. Typically, the borehole density is relatively small and especially deep boreholes are sparse. Often, borehole data are therefore not sufficient by themselves (Jørgensen and Sandersen, 2009a). A number of studies demonstrate that the spatially dense resistivity information from airborne electromagnetic (AEM) data can improve the mapping of buried valleys (Bosch et al., 2009, Jørgensen et al., 2012, Oldenborger et al., 2013) and that buried valleys can be successfully mapped when combined with stratigraphic and lithological information from boreholes. While 2D mapping of buried valleys have been performed and described from several studies (e.g. Gabriel et al., 2003, Jørgensen and Sandersen, 2009a, Jørgensen and Sandersen, 2009b), 3D modelling of the valleys has not yet been reported in literature to the same degree. The handling and interpretation of the very large data sets from areas with complex buried valley networks require formulation of an advanced modelling strategy. Thus, the cross-cutting nature of the valleys, their abrupt terminations, and generally heterogeneous infill facies make the valleys difficult to model with a traditional layered modelling approach, and more advanced methods therefore have to be used.
In the current study, we operate in an area that is characterised by a complex network of buried valleys. The geological modelling is based on densely spaced AEM data and borehole information available prior to this study. Here, we present a 3D modelling strategy with which we have produced a detailed 3D geological model where a large number of cross-cutting buried valleys are resolved and included as discrete elements.
The study area covers 45 km2 and is located in the central part of Denmark (Fig. 1). This part of Denmark has been covered by numerous ice sheets from different directions during the Pleistocene (Houmark-Nielsen, 2004). Within the study area, the terrain surface varies from sea level to almost 100 m above sea level.
Based on earlier mapping projects (e.g. Jørgensen and Sandersen, 2006), the geology in the uppermost 200 m can be divided into two major geological elements: i) thick deposits of hemipelagic Palaeogene clays and marls that generally occur at shallow depths and ii) deeply (up to 130 m) incised tunnel valleys arranged in a complex network. In most of the area, Paleogene clays and marls directly underlie the Quaternary deposits, but sandy Miocene sediments occur locally. Most of the buried valleys are incised into the Paleogene deposits, which accordingly show a considerable variation in surface relief. The Quaternary deposits are mainly of glacial origin and composed of clay till and fine to coarse-grained meltwater deposits. Meltwater sand and gravel within the buried valley structures constitute important aquifers used for regional drinking water supply.
The airborne transient electromagnetic data (Fig. 1) were acquired with the SkyTEM system (Sørensen and Auken, 2004). In this survey, the early-time data have been corrected for the primary field response, using a method comparable to the one in Schamper et al. (2014). This allowed unbiased data to be obtained from as early as 7 μs after beginning of ramp down (equivalent to 3 μs from end of ramp) for the super low moment. The last gate was centred at 9 ms. The very early times carry information about thin layers, in the order of a few meters, in the upper c. 20 m, similar to the performance achieved previously with other inversion techniques (Schamper et al. (2014)). At the same time, layers in more than 250 m depth can be mapped, although with both decreasing lateral and vertical resolution and highly dependent on the resistivity structure. The survey covers 333 line km with a line spacing of 100 m. Due to the dense infrastructure in the area, many of the soundings suffered from couplings to man-made conductors and were removed during the data processing, and prior to the inversion. Failure to eliminate these data would have inevitably resulted in artefacts in the resistivity models recovered (Viezzoli et al., 2013). We estimate that approximately 30% of the acquired data had to be culled. Consequently, about 9500 soundings were used for the inversion. The inversion was based on the spatially constrained (SCI) quasi-3D approach: a local, 1D exact forward response (Auken et al., 2014), with model parameters spatially constrained in 3D (Viezzoli et al., 2008). The inversion was performed with 25 layers that provided smooth transitions in the resulting models. The resistivity models were then further interpolated into horizontal grids and stacked, resulting in a 3D resistivity grid.
Borehole information from 435 boreholes (Fig. 1) was extracted from the Danish national borehole database, Jupiter (Møller et al., 2009). Most of the boreholes are shallow, with 63% of the boreholes being no more than 35 m deep and with only 3% of the boreholes deeper than 100 m. The majority of the boreholes have been drilled for water abstraction, and the deep boreholes are therefore typically clustered in areas around waterworks. The boreholes vary considerably in quality due to e.g. age and drilling method. The impact from these factors is evaluated during the geological modelling process.
Mapping and modelling strategy
When modelling in areas with complex geology, we typically operate with three main phases: performing a basic geological interpretation and setting up a general conceptual model, performing the layer modelling, and finally, performing the voxel modelling. When modelling buried valleys, it is necessary to follow a dedicated modelling strategy in order to address the highly complex geology. This strategy encompasses five main steps: i) constructing a conceptual geological model of the buried valleys, ii) determining the age relationships between the valleys, iii) layer modelling of the valley floors and sides, iv) trimming of the valley floor surfaces in accordance with the age relationships, and v) building a voxel model and populate it with lithofacies and stratigraphic units
Basic geological interpretation and conceptual model
The resistivity models recovered from the inversion of the SkyTEM data were thoroughly interpreted together with borehole data and the mapping was carried out by following the guidelines described in Jørgensen et al. (2003). Apart from being visualised as a 3D grid, the resistivity data were jointly displayed in vertical sections and horizontal slices (Fig. 2A). The buried tunnel valleys were detected by the significant resistivity contrasts between the infill and the incised surroundings. Particularly good results were achieved where the valleys were incised into the low-resistive substratum, because the valleys are filled with deposits characterised by higher resistivities. To identify the shapes (floor and walls) of these valleys, a surface map of the deepest conductive layer was extracted (Fig. 2B). Other valleys that are cut into sandy sediments were possible to identify and delineate if they are filled with low-resistive sediments. Some valleys with sandy fill and without any significant contrast to the substratum were somewhat imprecisely mapped and some valleys may remain undetected. The valleys were in general seen as elongate, often irregular, anomalies in horizontal slices through the resistivity data, either as high-resistive bodies in low-resistive environments or as low-resistive bodies in high-resistive environments. A basic, conceptual understanding of the area’s geological history was developed and used as a guideline during the subsequent construction of the 3D geological model (Fig. 3). All the mapped valleys are shown in map view in Fig. 2C.
The geological modelling was made in the software, Geoscene3D (I-GIS, 2014). The SkyTEM resistivity data were displayed both as 1D models and as an interpolated 3D resistivity grid. During the basic geological interpretation, the resistivity grid was sliced horizontally at different levels and a surface map of the deepest conductive layer was studied in order to find the exact extent and position of the valley structures. The modelling was mainly performed along vertical sections, where lithological information from boreholes were shown together with SkyTEM resistivity data within a predefined buffer zone, typically 100 m (Fig. 4A,B). The profiles were initially drawn perpendicular to the valleys and gradually moved 100 m along the valleys. They were continuously checked against profiles drawn along the thalwegs of the valleys as mapped during the basic interpretation (Fig. 2C).
The first step was to model the surface of the valley floors and walls. The appearance of the resistivity structures and the lithological information from the boreholes were used to place the interpretation points. These points were subsequently gridded using kriging to obtain the floor of each valley. Most of the buried valley structures reach the terrain surface, and the gridded valley surfaces were accordingly cut at this level. However, due to later erosional events, the buried valleys are not possible to trace in the present-day terrain. A thin sheet of glacial sediments (less than 2–4 m) seems to exist above most of the buried valleys that reach the surface, but this thin sheet does not clearly occur in the data due to limitations in resolution capability of the TEM method and degree of detail in the borehole descriptions (Schamper et al., 2014).
During modelling of the valley floors, the age relationships between the individual valleys were determined. This was done based on the structural and compositional characteristics of the infill, cross-cut relationships, orientation of the valleys, and on the background geological knowledge from the area. The age relationship was imperative during the subsequent trimming (cutting) of the valley surfaces. During the trimming, the surfaces of older valleys surfaces were cut by valley surfaces of younger valleys. This is for instance illustrated in Fig. 4B, where valley 5 cuts down into valley 2. The final, trimmed valley surfaces are shown in Fig. 5.
The final step of the layer modelling was to construct surfaces of the top of the Paleogene clay and the Miocene sand. As for the valley surfaces, this was done using interpretation points and subsequent kriging interpolation. In areas with incised buried valleys, the bottoms of the valleys also constitute the top of the Paleogene. The Paleogene clay shows very low resistivities and is easy to recognise in the models obtained from SkyTEM (Fig. 4B). The Miocene is on the other hand impossible to discern in the resistivity models, because these sandy sediments do not show resistivity contrasts against the Quaternary sandy sediments. Modelling of the Miocene unit is therefore solely based on borehole information.
The voxel model covers the entire model area from 180 m below sea level to the ground level at 0–100 m above sea level. The voxel grid discretisation is 50 m laterally and 5 m vertically. These dimensions were chosen to achieve a proper resolution of the mapped valley structures.
The property units of the voxel model were selected on the basis of the basic geological interpretation and on the encountered lithofacies in the boreholes. The SkyTEM resistivity information was interpreted into lithologies by considering the typical formation resistivities, knowing that the natrium chloride content of the groundwater in the area is low (Søndergaard et al., 2004). The Pre-Quaternary sediments were grouped into Paleogene clay and Miocene sand. The Quaternary was grouped into deposits occurring outside the valleys (sand till, clay till, and meltwater sand) and valley infill deposits (meltwater sand, meltwater clay, and clay till). The units were thus defined both on basis of lithology and valley stratigraphy, which resulted in a total of 42 geological units.
The voxel modelling was performed using a manual approach, where the geological units were interpreted without any automated routines. This was done by using a set of specific tools in GeoScene3D, by which it is possible to select and populate the voxels in different ways (Jørgensen et al., 2013). In the current work, we have primarily used a tool that enables the selection of voxels within volumes delineated by surfaces, and a tool where voxels were selected within digitised polygons on cross-sections. Using these tools, major volumetric bodies were populated by picking groups of voxels that were constrained by the initially modelled surfaces, whereas minor bodies were manually digitised.
The basic interpretation of the data resulted in a conceptual model for the area with several buried tunnel valleys cutting down into the Palaeogene clay and into one another, giving a “cut-and-fill” architecture in vertical section view (Fig. 3). Many of the valleys are deep, but also shallow valleys occur. The resistivity distribution and borehole information indicate in both map view (Fig. 2A) and section view (Fig. 4A,B) that the valleys are filled with both sandy (high resistivities) and clayey deposits (low resistivities). For instance, valley nos. 2, 11, and 14 (Fig. 4B) are filled with deposits with low to very low resistivities corresponding to clay, whereas nos. 1, 5, and 12 are filled with resistive, sandy deposits. The resistive layers (> 50–60 ohmm) are found to be coarse-grained glaciofluvial deposits, whereas the low-resistive layers correspond to glacial clay, either clay till or glaciolacustrine clay. Most valleys are filled with varying sediments both vertically (e.g. no 1 in Fig. 4B) and laterally. The very low resistivities (< 12 ohmm) outside the valleys correspond to Palaeogene clay.
A total of 20 buried valleys were integrated in the layer model. The dominant orientation, depth, width, and interpreted generation of the valleys are presented in Table 1. When shown in 3D (Fig. 5), it is seen how the valley surfaces cut into one another resulting in a complex buried-valley network. Most of the valleys show a preferred orientation from SW-NE (Table 1, Fig. 2C). The depths of the valleys range from a few tens of metres to 170 m, and the width generally falls within the range of 400–1100 m (Table 1). The lengths of the valleys are difficult to estimate, since most of the valleys are expected to continue outside the modelled area. The shapes of the valleys are clearly irregular with hollows and thresholds along their thalwegs (Fig. 2B and 5). As noted above, most valley surfaces have been trimmed at ground level indicating that this constitutes a major erosional boundary.
According to the analysis of the age relationships, the buried valleys fall within at least eight different generations (Fig. 6). Different valley generations may have been formed during the same glaciation by different ice advances (Jørgensen and Sandersen, 2006), but since we see a number of generations with very different orientations it is likely that they were formed during more than just one glaciation. This is in accordance with (Jørgensen and Sandersen, 2006, Jørgensen and Sandersen, 2009b), who argue that the valleys in the area are of Weichselian to Elsterian age, perhaps even pre-Elsterian. The number of valleys belonging to each generation varies. Some of the valleys (e.g. nos. 3 and 8, Fig. 2C and 6) show abrupt bends of up to 90 degrees over short distances, possibly because the valleys tend to form where the substratum is most easily eroded (Janszen et al., 2012, Sandersen and Jørgensen, 2012), which in this case is the older coarse-grained valley infill instead of the less erodible Paleogene clay deposits.
The 42 different geological units in the 3D voxel model can be thematised with regard to either valley stratigraphy (Fig. 4C, 7) or lithofacies (Fig. 4D, 8). An overview of the voxel model results is shown in Fig. 7, Fig. 8 as a fence diagram with a number of vertical N-S sections through the final voxel model. The base of the model is constituted by Paleogene clay in the entire model area. Above this, a thin layer of Miocene sand is present in the central part. This is followed by the heterogeneous Quaternary which covers the entire area. The glacial deposits are grouped into valley fill and deposits that occur outside the buried valleys (Fig. 8). Outside the valleys, the deposits are mainly composed of clay till, but with frequent occurrences of meltwater sand inclusions. Each valley is typically filled with a predominant lithology with minor occurrences of other lithofacies (Fig. 8). Slightly more than half of the valleys contain meltwater sand as the main lithofacies, whereas the rest primarily contain clay till. Only one valley (no. 2) has meltwater clay as the predominant lithofacies.
Mapping and modelling complex geological environments like buried-valley networks require dense data coverage, which can be obtained through airborne electromagnetic data (Jørgensen and Sandersen, 2009a, Sapia et al., 2014). Such data are typically very complex, provide lots of information, and are difficult and very time consuming to understand and interpret. Proper geological modelling of the buried-valley networks using large-scale AEM data therefore requires a clear strategy with a well-defined workflow. Buried valleys are difficult to model with a traditional layer-based approach due to a number of factors. Firstly, the cross-cutting nature of the valleys makes them complicated to interpret and model with respect to the age relationships of the valley generations. Secondly, the abrupt terminations and sudden sharp bends of the valley structures make them difficult to follow and map. Finally, as most valleys are filled with rapidly shifting material, several internal layers within the valley fill would be required to model the actual geological setting. To be able to include these strong internal variations the manual voxel modelling methodology proved to be advantageous for modelling the valley fill.
Recently, research has focused on developing automatic or semi-automatic routines for hydrostratigraphical modelling of AEM data with the objective to create a faster routine and avoid the subjective element inherent in the interpretations (e.g. Foged et al. (2014) and Gunnink and Siemon (2015)). Høyer et al. (2015) made a comparison of the cognitive, manual approach and two different automatic modelling methods. Here, the manual model was the only one that could provide stratigraphical information, but all the model approaches resulted in highly comparable results if translated into binary sand-clay models. It must be noted though that in Høyer et al. (2015), the geological setting was relatively homogeneous compared to the buried valley-network in the current study, where the stratigraphical information is vital for the model output. Since the stratigraphical information is unobtainable using automatic methods, we therefore expect more significant differences when operating in this type of environment. The geological knowledge required for making the geological model includes subjects like regional geology and the glacial history of the area, without which it would be impossible to infer the different valley generations. Another important input is knowledge regarding glacial processes and the formation of the buried valleys, which is essential for delineating realistic shapes and extents of the structures. This is especially crucial when modelling the shallow valleys that are incised into other glacial deposits and therefore only show low resistivity contrasts. Also the geophysical knowledge regarding resolution capabilities is vital when interpreting the AEM information. Generally, the vertical and the horizontal resolution diminish drastically with depth, such that deep geological structures have to be considerably larger in order to be resolved than structures located at shallower depths. Resistivity variations within resistive deposits such as meltwater sand or gravel are more difficult to identify using AEM compared to variations in conductive deposits, due to the inductive nature of the methodology (Fitterman and Stewart, 1986, Schamper et al., 2014). Another important aspect is the appearance of sloping structures in the geophysical sections. These tend to be underestimated in the inverted resistivity sections if the data are inverted using a 1D inversion approach. This effect is most pronounced for the steepest slopes, which also can cause 2D/3D effects (Danielsen et al., 2003). Most of these considerations are impossible to take into account in automatic modelling approaches, which is why manual interpretation is preferred when the product is a geological model with stratigraphic content.
In this study, the deepest valley structures were easily delineated by investigating the depth to the deepest conductive layer, which corresponds to the boundary between the Paleogene clay and the Quaternary deposits. Modelling of the deep valley structures is therefore associated with only little uncertainty and would also be outlined using a proper automatic modelling method. The shallow valleys, on the other hand, were considerably more difficult to interpret and model and could only be delineated by using the geological knowledge of the area. The lower resistivity contrasts between the shallow valleys and their surroundings will induce model uncertainties. In this study, we have been able to map and model several valleys, but even more valleys may be present in the area. First of all, buried valleys that are too small to be resolved in the data may be found, and secondly, buried valley structures can be unresolved in both resistivity data and borehole data, if the sediments lack resistivity and lithology contrasts. Modelling these valleys require other data types such as reflection seismic data, by which it is possible to resolve the erosional boundaries of the valleys (Jørgensen et al., 2003).
The lithological information delivered by the geological model is crucial, when the model is going to be used for groundwater modelling. While some of the sand-filled valleys constitute important aquifers, the clay-filled valleys might act as barriers for groundwater flow. An example of this is seen in the western part of the area, where a clay-filled buried valley (no. 11) cuts through a valley filled with sandy deposits (valley no. 1, Fig. 4). The narrowing of the aquifer is expected to have a significant influence on the resulting transmissivity, such that the water flow is diminished across the cross-cutting valley (no. 11). Also, the stratigraphic information can be vital with regard to groundwater modelling, and this information is only obtainable using the manual modelling approach. A certain sediment type deposited under different conditions will typically not show the same sorting degree or grain size distribution. Consequently, the hydraulic conductivity of a deposit like meltwater sand is expected to differ, when deposited under different valley infill events. Furthermore, the depositional history in the individual valleys can hold information regarding anisotropy of the valley infill, and thereby anisotropy for the resulting hydraulic conductivities. The detailed information contained in the geological model is consequently highly valuable, when transferring the model into a hydrogeological model for further groundwater modelling.
In this study, we have formulated and followed a strategy for 3D geological modelling of a geological environment characterised by a network of buried valleys. The approach is based on 5 steps: i) The AEM-derived resistivity models are examined to identify the buried valleys in the area, and the resistivity data and the geological background knowledge are utilised to build a conceptual geological model. ii) The internal age relationship between the valleys is determined from their direction and cross-cut relationships. iii) The valley floors and sides are modelled using a layer modelling approach. iv) The internal age relationship is used to trim the valley floor surfaces. v) Finally, a voxel model is built and populated with lithofacies and stratigraphical units. The result is a 3D geological voxel model that, for this case study, contains 42 different geological units. A high number of units have been chosen, because the Quaternary deposits are grouped into deposits occurring outside and inside the individual valleys, and most of the valleys comprise more than one lithology. This grouping was made, since the hydraulic conductivities of the different units are influenced by their depositional characteristics. The model contains 20 buried valleys that show a complex cross-cut relationship, which suggests at least eight different valley generations. Some of the valleys constitute important aquifers, but about half of the valleys contain clayey deposits that in some cases may act as barriers for the groundwater flow. Consequently, the modelling approach developed in this study resulted in a highly detailed model of both lithologies and stratigraphy. Apart from being applicable in buried valley settings, the same methodology may be applicable in other complex geological environments, where the conceptual geological understanding is crucial for the model output.
Source: 3D geological modelling of a complex buried-valley network delineated from borehole and AEM data
Authors: A.-S.Høyer, F.Jørgensen, P.B.E.Sandersen, A.Viezzoli, I.Møller