Bedrock mapping of buried valley networks using seismic reflection and airborne electromagnetic data

Journal of Applied Geophysics 128 (2016) 191–201

Contents lists available at ScienceDirect

Journal of Applied Geophysics journal homepage: www.elsevier.com/locate/jappgeo

Bedrock mapping of buried valley networks using seismic reflection and airborne electromagnetic data G.A. Oldenborger a, C.E. Logan a, M.J. Hinton a, A.J.-M. Pugin a, V. Sapia b, D.R. Sharpe a, H.A.J. Russell a a b

Geological Survey of Canada, Natural Resources Canada, 601 Booth St., Ottawa, Ontario, Canada K1A 0E8 Istituto Nazionale di Geofisica e Vulcanologia, Via di Vigna Murata 605, 00143, Rome, Italy

a r t i c l e

i n f o

Article history: Received 23 December 2015 Received in revised form 22 February 2016 Accepted 16 March 2016 Available online 21 March 2016 Keywords: Airborne electromagnetics seismic reflection hydrogeophysics geological modelling

a b s t r a c t In glaciated terrain, buried valleys often host aquifers that are significant groundwater resources. However, given the range of scales, spatial complexity and depth of burial, buried valleys often remain undetected or insufficiently mapped. Accurate and thorough mapping of bedrock topography is a crucial step in detecting and delineating buried valleys and understanding formative valley processes. We develop a bedrock mapping procedure supported by the combination of seismic reflection data and helicopter time-domain electromagnetic data with water well records for the Spiritwood buried valley aquifer system in Manitoba, Canada. The limited spatial density of water well bedrock observations precludes complete depiction of the buried valley bedrock topography and renders the water well records alone inadequate for accurate hydrogeological model building. Instead, we leverage the complementary strengths of seismic reflection and airborne electromagnetic data for accurate local detection of the sediment-bedrock interface and for spatially extensive coverage, respectively. Seismic reflection data are used to define buried valley morphology in cross-section beneath survey lines distributed over a regional area. A 3D model of electrical conductivity is derived from inversion of the airborne electromagnetic data and used to extrapolate buried valley morphology over the entire survey area. A spatially variable assignment of the electrical conductivity at the bedrock surface is applied to different features of the buried valley morphology identified in the seismic cross-sections. Electrical conductivity is then used to guide construction of buried valley shapes between seismic sections. The 3D locus of points defining each morphological valley feature is constructed using a path optimization routine that utilizes deviation from the assigned electrical conductivities as the cost function. Our resulting map represents a bedrock surface of unprecedented detail with more complexity than has been suggested by previous investigations. Our procedure is largely data-driven with an adaptable degree of expert user input that provides a clear protocol for incorporating different types of geophysical data into the bedrock mapping procedure. Crown Copyright © 2016 Published by Elsevier B.V. All rights reserved.

1. Introduction Buried valleys are common across formerly glaciated parts of Canada, the northern United States, and northern Europe, among other places (Kehew and Boettger, 1986; Parks and Bentley, 1996; Piotrowski, 1997; Huuse and Lykke-Andersen, 2000; Cummings et al., 2012; Kehew et al., 2012). In the North American prairies, they are often considered high yield sources of groundwater where low yield bedrock and mud-rich diamicton otherwise dominate the hydrogeological setting (Winter et al., 1984; Shaver and Pusc, 1992). Buried valleys may be regional depressions associated with pre-glacial or glacial drainage patterns, that are filled and buried by preglacial, glacial and postglacial sediments, thereby lacking surface expression. Furthermore, buried valley systems exhibit complicated network geometries, multiple scales, and longitudinal and cross sectional variability

E-mail address: [email protected] (G.A. Oldenborger).

http://dx.doi.org/10.1016/j.jappgeo.2016.03.006 0926-9851/Crown Copyright © 2016 Published by Elsevier B.V. All rights reserved.

such that they may be incompletely mapped or characterized by traditional hydrogeological techniques (Jørgensen and Sandersen, 2009; Abraham et al., 2010; Stewart and Lonergan, 2011). Despite productive buried valley aquifers being common in Canada, knowledge of their distribution and groundwater resource potential is often inadequate (Russell et al., 2004; Betcher et al., 2005; Ahmad et al., 2009; van der Kamp and Maathius, 2012). Detailed mapping of the bedrock surface is a critical first step in the hydrogeological investigation of buried valley networks. Bedrock topography reveals buried valley locations, valley morphology and erosional events, such as cross-cutting valleys that provide information on process formation of aquifer systems (Kehew and Boettger, 1986) and buried valley aquifer discontinuities (Shaver and Pusc, 1992). Bedrock typically serves as the major basal aquitard in the prairie Quaternary aquifer sequence (Cummings et al., 2012). In prairie regions where bedrock is siliceous shale, the uppermost portion of bedrock may be heavily fractured. In these regions, the bedrock surface defines the top of a large-scale low-yield aquifer itself (the contact aquifer) sitting

192

G.A. Oldenborger et al. / Journal of Applied Geophysics 128 (2016) 191–201

above the regional aquitard. This contact aquifer often marks a transition to increasing salinity with depth and formation-specific potability (Grasby et al., 2014). Methods of subsurface investigation utilizing water wells, borehole coring and logging, or ground geophysics provide localized information on bedrock depth, but are generally not of sufficient density for detailed regional buried valley detection, delineation and mapping (Jørgensen and Sandersen, 2009). To achieve both regional coverage and high data density, techniques involving airborne electromagnetics (AEM) have been developed to map Quaternary aquifer systems including buried valleys (Gabriel et al., 2003; Sørensen and Auken, 2004; Steuer et al., 2009; Siemon et al., 2009). However, electromagnetic methods have finite penetration depth that depends on the geological conditions, in addition to suffering from depth-, conductivity- and thicknessdependent resolution (Goldman et al., 1994; Auken et al., 2008) and system-specific limitations that may affect interpretations (Christiansen et al., 2011). Several studies have investigated the integration of additional high-resolution data types, such as boreholes or seismic data, to reduce ambiguity associated with electrical conductivity models and geological interpretations derived from AEM surveys (Jørgensen et al., 2003; Høyer et al., 2011; Jørgensen et al., 2012; Foged et al., 2013; Pugin et al., 2014; Reninger et al., 2014; Sapia et al., 2014a). In some cases, the additional ground-based geophysical data are used as local “ground truth” for the AEM results, and in other cases, the supporting data are formally incorporated into inversion of the AEM data. Regardless, the area of influence of the supporting data is necessarily limited, such that the regional result is a cognitive interpretation of all available data (Jørgensen et al., 2013; Høyer et al., 2015; Sapia et al., 2015; Steinmetz et al., 2015). Such cognitive or knowledge-driven approaches are useful for dealing with multiple disparate data types, data inconsistencies, and problems associated with AEM resolution and non-unique correlation of conductivity to lithology in complex sedimentary environments. Unfortunately, cognitive interpretations can be subjective, laborious, difficult to reproduce and require a high level of expertise, all of which cause problems for technology transfer. Alternatively, we develop a largely data-driven bedrock mapping scheme which leverages highresolution, but localized seismic reflection data that are extrapolated using regionally extensive AEM data to determine the bedrock surface. When applied to data from the Spiritwood buried valley aquifer system in Manitoba, the resulting bedrock topography map reflects a complex buried valley network that cannot be resolved with water well records, boreholes, or any single geophysical data set.

(Nicolas et al., 2010; Bluemle, 1983). Buried valley fill consists of a series of mud-rich till units interstratified with gravel, sand, silt and clay of variable thickness and extent (Winter et al., 1984). The Spiritwood aquifer system is used for agricultural, domestic and municipal purposes and is among the highest yielding aquifers in North Dakota (Paulson, 1983). The current study area extends approximately 50 km north from the Canada-USA border within a till plain of little topographic relief (Fig. 1). Seismic reflection and AEM data collected for the Spiritwood buried valley in this area have revealed a complicated bedrock surface and aquifer geometries with multiple erosional surfaces and valley morphologies (Oldenborger et al., 2013). Borehole and water well records in the region do not reflect the spatial heterogeneity observed in the geophysical data and the buried valley network was heretofore mapped as a broad 10–20 km wide bedrock valley with variable distribution of productive and low-permeability units and poorly understood recharge (Sie and Little, 1976; Randich and Kuzniar, 1984). Conceptual geological models based on the borehole information alone would not accurately support analysis of groundwater flow, yield and resource potential. Construction of an effective 3D hydrogeological model requires a clear protocol to incorporate high-resolution and spatially-extensive geophysical data into mapping bedrock topography. 2. Methodology 2.1. Data sources The data utilized for bedrock mapping consist of water well records from Manitoba Conservation and Water Stewardship and the North Dakota State Water Commission, approximately 3000 km of AeroTEM III data over 1060 km2, and 63 km of high-resolution landstreamer seismic reflection data (Oldenborger et al., 2013). Borehole logs, interpreted cross-sections and surficial geological maps were also used where available (Randich and Kuzniar, 1984; Crow et al., 2012a, 2012b). The water well records are normalized to a common lithostratigraphic legend based on textural descriptions and location in the stratigraphic sequence (Logan et al., 2015). In addition to bedrock, the lithologic

1.1. Study area The Canadian Prairie landscape has been shaped by two continentalscale events during the past 50 million years. Tertiary uplift of the Canadian Rockies was associated with large drainage systems that flowed north-eastward, carved valleys in the strata and transported material from the sedimentary basin. Over the Quaternary, continental glaciers traversed the Prairies and eroded and deposited a succession of sedimentary deposits up to 300 m thick of predominantly mud-rich till with glacial-lacustrine sediment and inter-till sands and gravels (Cummings et al., 2012). Buried valleys are a product of these two continental events and have been identified across the Canadian Prairies from the foothills to the eastern edge of the Western Canadian Sedimentary Basin with outcrops, water well records, borehole logs, seismic shot holes, and petroleum wells comprising most of the data for existing maps (Grasby et al., 2014). Poorly consolidated Cretaceous shale and minor sandstone form the main bedrock substrate of Prairie buried valleys (Maathuis and Thorleison, 2000). The Spiritwood is a trans-border buried valley aquifer system extending approximately 500 km from Manitoba, across North Dakota and into South Dakota. The Spiritwood buried valley is eroded into fractured siliceous shale from the Odanah Member of the Pierre Formation

Fig. 1. Location and DEM of the Spiritwood buried valley survey area in southern Manitoba, Canada including locations of water wells, seismic reflection survey lines and boreholes with lithological information. Landsat 7 orthoimage © Natural Resources Canada.

G.A. Oldenborger et al. / Journal of Applied Geophysics 128 (2016) 191–201

units in order of abundance in the water well records are till, inter-till sand and gravel, surficial glaciofluvial, alluvial, sand and gravel, and glaciolacustrine. The wash-boring technique used to economically drill most residential water wells is subject to a large degree of uncertainty in terms of the textural observations and depth of contacts which are often based on cuttings, drilling rates and resistance. Identification of bedrock is likely the most accurate water well observation, but is still subject to uncertainty that is difficult to quantify because of the high variability in skill, experience and rigor of individual drillers. Water well locations are recorded based on quarter-sections which results in a practical spatial precision of approximately 0.6 km. The AeroTEM data were filtered for cultural noise and inverted using a psuedo-3D spatially constrained inversion to yield electrical conductivity (Viezzoli et al., 2008; Sapia et al., 2014b). The inversion results consist of a collection of irregularly-spaced 1D models with 29 logspaced layers to 200 m depth. To contextualize and visualize the AEM results in terms of regional patterns, the 1D models are interpolated (natural neighbour) to a regular horizontal grid (100 × 100 m2) for each layer depth to yield a 3D electrical conductivity model (e.g., Pryet et al., 2011). The interpolated 3D AEM electrical conductivity model is shown in Fig. 2. The shale bedrock represents a conductive target relative to the overlying sediment package, and variations in conductivity are interpretable in terms of complex buried valley geometry. To compare electrical conductivity to other observations, we use 3D linear interpolation to project the electrical conductivity model to observation locations. Seismic reflection data were acquired with a Minivib vibratory seismic source and a landstreamer receiver array with 48 geophones spaced at 1.5 m (Oldenborger et al., 2013). The data are processed to produce compressional-wave and/or shear-wave reflection sections with time-depth conversion performed using interval velocity obtained from semblance analysis of the common-shot gathers (Pugin et al., 2009a, 2009b). In many areas, particularly where the near-surface sediments are water-saturated or consist of dense till, good quality P-wave reflections are visible from the glacial sediments and the bedrock surface. In areas of poor P-wave transmission, such as for a deep water table, often only S-wave reflection data are recovered which may not exhibit strong bedrock reflections. Bedrock reflections may also be subtle when dense till overlies shale with only a small contrast in velocity, or when internal multiples within the sedimentary sequence interfere with bedrock reflections. For these data, the bedrock surface is interpreted based on a transition in reflection facies from high

193

amplitude to low amplitude with poor penetration (Pugin et al., 2014). For both P- and S-wave sections, the phase and interval velocity are also used for interpretation of the bedrock surface. Example P-wave and S-wave reflection sections with the interpreted bedrock surface are shown in Fig. 3. The theoretical quarter-wavelength resolution for the seismic data is on the order of 3 m for the P-wave sections and 1.5 m for the S-wave sections. Where higher quality monitoring well or borehole information is available, we observe the interpreted bedrock surface to be in good agreement with lithological information. Although we are focussing on the bedrock surface, the seismic data also reveal significant heterogeneity within the glacial sediments. Full characterization of the reflection section aids in interpretation of the bedrock surface and is important for understanding the Quaternary aquifer system (e.g., Pugin et al., 2014; Steinmetz et al., 2015). Additional ground-based electrical and electromagnetic data supporting the interpretations are presented in Oldenborger et al. (2013) and Oldenborger and Brewer (2014). 2.2. Conductivity classification Despite visually apparent discrimination of different material types, the electrical conductivity model is not easily separable in terms of observed lithostratigraphy. The distribution of electrical conductivity over the model volume, and the relationship between the electrical conductivity and the water well lithostratigraphy are illustrated in Fig. 4. The electrical conductivity model is subject to the limitations of inversion including artefacts, non-uniqueness and regularization (Day-Lewis et al., 2005). Any conductivity-lithology correlation (including that of the bedrock surface) will be resolution-dependent as exemplified by the smooth nature of the conductivity distribution (Fig. 4a). Furthermore, electrical conductivity is a function of many variables including material composition, texture and pore water content such that the conductivity-to-lithology mapping is not be 1:1. A single conductivity may represent multiple lithologies, and a distinct lithology does not have a unique conductivity (Fig. 4b). Based on depth to bedrock reported in water well observations, the electrical conductivity of bedrock appears largely separable from that of the overlying glacial sediments due to the higher conductivity of the shale. To determine the electrical conductivity associated with the bedrock surface, we examine the distribution of electrical conductivity at the bedrock surface as derived from water well records and all seismic sections (Fig. 5). We see that the bedrock/sediment interface exhibits

Fig. 2. Constant-elevation slices of the Spiritwood buried valley AEM electrical conductivity model at a) 445 masl, b) 415 masl and c) 385 masl corresponding to approximate depths of a) 30 m, b) 60 m and c) 90 m. Locations of 807 water wells within the AEM survey area and seismic survey lines are shown along with interpreted features of the buried valley network.

194

G.A. Oldenborger et al. / Journal of Applied Geophysics 128 (2016) 191–201

Fig. 3. a) P-wave seismic reflection section S7E. Monitoring well Kilcart 8 encounters bedrock at 55 m depth. b) S-wave seismic reflection section S3; coherent deep reflections are most apparent near 4250 m and 5000 m line position. Monitoring well OA011 encounters bedrock at 56 m depth. c) P-wave seismic reflection section S2. c) P-wave seismic reflection section S1. Borehole GSCBHSW07 does not encounter bedrock at 97 m depth. Transparent line is the bedrock surface interpreted from the seismic data. Insets show section locations on the 415 masl elevation slice of the AEM electrical conductivity model (Fig. 2b).

a significant range in conductivity across the survey area, possibly due to physical variations in pore water salinity and fracture density, but also due to vertical and lateral smoothing of the conductivity model. From these independent water well and seismic observations, we select a target electrical conductivity of approximately 66 mS/m that best represents the bedrock surface. A section of the electrical conductivity model along S7E in Fig. 6 shows that the target conductivity at the bedrock interface is intermediate between the conductivity of the overlying sediment and that of the bedrock itself, which is likely greater than 100 mS/m. Below the resistive fill of the buried valley, the target conductivity occurs below the depth of model reliability as measured by inversion sensitivity (Christiansen and Auken, 2012). Using the depth of the target electrical conductivity, we derive the bedrock surface illustrated in Fig. 7. As expected, the bedrock surface determined from a single target conductivity results in under- and over-prediction of bedrock depth in comparison to observations at water wells and along seismic sections. However, we also have a bias in our prediction of target depth. Of the 807 water wells within the AEM survey area, only 144 encounter bedrock and only a few of those encounter bedrock within the deeper regions of the buried valley network. This is a natural consequence of a drilling campaign that stops

at the first useable groundwater resource, but it results in limited representation of buried valleys and a bias of the target conductivity since shallow bedrock often appears more conductive in the AEM model. We also have additional model-dependent limitations in using a target conductivity in that the depths of incised valley floors are consistently overestimated. This is apparent along section S7E where the target conductivity occurs very near the bedrock surface on the valley flanks, but several tens of meters below the bedrock surface at the deepest part of the valley floor (Fig. 6). Furthermore, along much of the thalweg of the deep incised valley, the target conductivity occurs below the reliable depth of investigation, resulting in an unreliable or undefined bedrock surface in critical areas (Fig. 7). In some regions, the target conductivity may represent a lithology other than bedrock (Fig. 4). In other regions where the conductivity model exhibits conductive surface artefacts, the model may never fall below the target conductivity and the bedrock depth is undetermined or estimated as zero. Christensen et al. (2015) illustrate a methodology to interpolate bedrock topography from AEM conductivity for a bedrock interface with (apparent) spatially variable electrical conductivity. However, this method relies on a network of dense, evenly distributed

G.A. Oldenborger et al. / Journal of Applied Geophysics 128 (2016) 191–201

195

in key areas, correct for consistent over-prediction of incised valley depth, or compensate for finite depth of investigation. 2.3. Bedrock mapping

Fig. 4. a) Distribution of electrical conductivity for all 1D models on a log-based scale (no compensation for variable model layer thickness). The bimodal distribution is largely attributable to shale bedrock and glacial sediments. b) Boxplot showing the relationship between water well lithology and electrical conductivity projected to well locations.

geotechnical boreholes, as opposed to sparse water wells with a shallow bedrock bias. Furthermore, while potentially ameliorating some underand over-prediction, it would not address our paucity of borehole data

Fig. 5. Distribution of electrical conductivity projected to the bedrock surface from a) water well records and b) seismic reflection sections. For seismic sections, observations are counted by common-shot gather points.

To capitalize on the high-resolution, but localized nature of the bedrock contact interpreted from seismic data, we develop an integrated approach whereby we construct a set of 3D control lines describing the bedrock surface and buried valley morphology using both seismic reflection sections and the AEM electrical conductivity model. For each seismic survey line, we examine the reflection cross-section and assign four representative control points at inflection points along the buried valley profile. These control points define the lateral extents of the valley flank and floor and approximate the buried valley shape as shown in Fig. 8 for section S3. To form the complete buried valley, corresponding control points for each valley are connected between different seismic sections by control line segments that are determined using path optimization rules where adherence to the AEM electrical conductivity model provides the basis for the cost function. When building each control line, a constant-depth slice is extracted from the AEM model that best represents the average depth of a collection of control points for a particular morphological feature (e.g., the western flank of the eastern incised valley between sections S1 and S2). This 2D AEM depth slice forms the basis for the cost raster input into the “Cost Distance” and “Cost Path” tools in ArcGIS® Spatial Analyst®. The cost raster is defined as the absolute value of the difference between the 2D electrical conductivity and the average of the electrical conductivities at the start and destination control points. Starting from one assigned valley control point, the path optimization routine constructs a path to the next corresponding valley control point using an iterative node/link neighbour representation, a variant of the accumulated cost surface method (Douglas, 1994). The minimum cost path is defined as the shortest path which deviates least from the average of the corresponding electrical conductivities at the starting and destination control points as shown in Fig. 9 for sections S1 and S2. This method allows for arbitrary variation of the recovered electrical conductivity at the bedrock surface across a valley profile (Figs. 6 and 8), but assumes that electrical conductivity should be relatively constant for corresponding control points over short distances along the valley axis. The procedure is automated using Python scripting and repeated for all corresponding source and destination control points for each seismic section. In some areas, additional manual control points are required to guide path selection between widely-spaced seismic sections to prevent control lines from crossing valleys (Fig. 9a) or to prevent control lines from diverging at bedrock highs (Fig. 9b). In other cases, manual control points were added where incised buried valleys were interpreted to terminate or rise above the bedrock surface and be floored in the glacial stratigraphy (Fig. 9c). Furthermore, in addition to the main incised valley features interpreted from the AEM electrical conductivity model, several other continuous valley features are observed (Fig. 2). In particular, the broad, shallow buried valley is interpreted to pre-date all other channels and the smaller scale cross-cutting valleys are interpreted as the youngest valleys in the system (Pugin et al., 2014). To represent features that are observed in the AEM survey, but not captured by the water well records or seismic sections, additional manual bedrock control points were positioned in 3D based on joint interpretation of the AEM electrical conductivity model, horizontal gradients of the electrical conductivity model, and nearby seismic sections and/or water well observations. Control lines are built from these manual control points in the same manner as described above. However, selection of these manual control points represents a cognitive or knowledgedriven component of our procedure that can be carried out to varying degrees. We limit ourselves to the main valley features that are observed and marked in Fig. 2. The final result is a collection of control points and lines where each line represents a 3D locus of points that

196

G.A. Oldenborger et al. / Journal of Applied Geophysics 128 (2016) 191–201

Fig. 6. Electrical conductivity along seismic section S7E. Solid black line is the bedrock surface interpreted from the seismic data. Dashed black line is the sensitivity-based AEM depth of investigation. Grey line is the 66 mS/m contour.

defines the bedrock topography for the buried valley geometry as shown in Fig. 10. Once all bedrock control lines are established, they are discretized to points at 50 m spacing where depth along the 3D control line is assigned based on linear interpolation of the depths of the starting and destination control points. The bedrock surface interpreted from the seismic data (Fig. 3) is also converted to control points at every common-shot gather point and all bedrock control data are merged with the water well records. The occurrences of bedrock identified in the water well records and defined by the 3D control points/lines are jointly interpolated to form the bedrock surface illustrated in Fig. 11. 3. Results We compare the bedrock surface interpolated from water well data alone to that interpolated from water wells and the collection of control points and lines (Fig. 11). We utilize natural neighbour interpolation to avoid complication and bias associated with variogram estimation. For control data and the AEM electrical conductivity model, choice of interpolation algorithm is largely inconsequential due to data spacing. For water well data, the interpolation algorithm cannot compensate for the data sparsity. The water well bedrock surface represents a poor reconstruction of the buried valley network because few water wells

Fig. 7. Estimated bedrock depth derived from the AEM electrical conductivity model and the target conductivity of 66 mS/m. Bedrock estimates below AEM depth of investigation are outlined by the dashed white line.

reach bedrock and even fewer reach bedrock within the buried valley system. On their own, the water well records capture the existence, but not the boundaries of the broad Spiritwood buried valley (Fig. 2b) which is traditionally considered to be the Spiritwood buried valley aquifer. There are too few water well bedrock observations to suggest the presence of a complex nested buried valley system. In contrast, our control line methodology depicts the existence of the deep incised valleys interpreted from the seismic data as well as the broad Spiritwood valley and the numerous smaller-scale valley features that are introduced by the AEM-defined control points and lines. Equal weight is placed on all the data for interpolation, but for the most part, the control lines define the bedrock surface over the buried valley system where water wells are absent. The high-resolution seismic sections define precise lateral and vertical locations of valley features while the AEM conductivity model effectively dictates the path of those features between seismic sections. In areas where multiple observations are present, consensus between data is a natural consequence of our method since seismic observations are referenced to borehole observations and control lines are anchored at seismic observations and/or boreholes. To further increase the influence of the AEM survey on the bedrock map, we consider an additional step whereby, instead of a single bedrock target conductivity, we use the combination of water well records, and the control points/lines to interpolate a 2D bedrock target conductivity (Fig. 12a). At each point in the AEM electrical conductivity model, we then find the depth of the variable target conductivity (Fig. 12b). We filter this target depth to exclude regions below the depth of investigation. We also filter regions with target depth less than 5 m to exclude artefacts associated with early-time AEM data, and regions for which the target conductivity falls below 40 mS/m, the approximate minimum conductivity observed at the bedrock surface for the water well records

Fig. 8. Electrical conductivity along a portion of seismic section S3. Solid black line is the bedrock surface interpreted from the seismic data. Dashed black line is the sensitivitybased AEM depth of investigation. White dots represent the assigned control points defining the valley morphology. Grey line is the 66 mS/m contour.

G.A. Oldenborger et al. / Journal of Applied Geophysics 128 (2016) 191–201

Fig. 9. Example bedrock control line paths between seismic sections S1 and S2. Seismic control points are shown on the seismic survey lines. Addition manual control points are required to a) prevent control lines from crossing valleys, b) prevent control lines from diverging at conductivity (bedrock) highs, and c) force control lines to close where buried valleys are interpreted to terminate or surface. Bottom panel shows a region around S2 with and without the control points/lines (dashed white rectangle depicts the region). AEM electrical conductivity model is shown at 385 masl (approximately 90 m depth).

197

(Fig. 5). The latter limit excludes some incorrect identification of other lithologies as bedrock (Fig. 4). Target depth is then used in conjunction with water wells and control points/lines to interpolate the bedrock surface (Fig. 12c). Due to the disparate data density between the grid-based target depth versus the point-based water wells and control points/lines, it is necessary to use a gridded version of the bedrock surface in Fig. 11b as representative of the water well and control point/line data. Without this step, the composite interpolant simply honours water well or control point depths at those points exactly and then quickly reverts to target depth. The result is the addition of detail to the bedrock depth map between control lines and seismic sections, but the additional heterogeneity is likely inconsequential at the scale of hydrogeological modelling (Figs. 11b vs. 12c). In the absence of control line data, a 2D target conductivity built from water wells and seismic results in a target depth that is an improvement over that obtained with a single target conductivity, but still suffers from limitations of model resolution and over-prediction of buried valley depths. Within Manitoba, the northward depth trajectory of the central incised valley suggested that buried valley fill should outcrop in a tributary of the Souris River (Fig. 13). Satellite images and subsequent site inspections confirmed the presence of sand and gravel outcrops and deposits within the tributary valley. Baseflow measurements during drought conditions revealed sustained discharge at this location, which is interpreted as evidence of groundwater supply from the incised buried valley aquifer system as mapped herein. Given these observations, additional off-survey control lines were constructed to extend the incised valley from the northern edge of the survey area to the Souris River and its tributary, thereby hydraulically connecting the buried valley aquifer to surface water (Fig. 13). These additional control lines were combined with water well and borehole data to construct the buried valley bedrock topography to the limits of the regional watershed (Fig. 14). The bedrock topography map was extended into North Dakota in accordance with water well records and interpreted cross-sections (Randich and Kuzniar, 1984). The North Dakota cross sections are derived from sub-parallel borehole transects across the Spiritwood buried valley at 3–6 mile N-S intervals with less than 1 mile E-W spacing between boreholes. The bedrock surface along these cross sections was defined as control lines and incorporated along with North Dakota water wells, Manitoba water wells, control points/lines, and Pelican

Fig. 10. Bedrock control points (grey) and control lines (white) over the AEM survey area that represent the 3D locations of valley flanks and floors shown in relation to the AEM electrical conductivity model at a) 415 masl and b) 385 masl.

198

G.A. Oldenborger et al. / Journal of Applied Geophysics 128 (2016) 191–201

Fig. 11. Estimated bedrock depth over the AEM survey area derived from natural-neighbour interpolation of a) water well records and b) water well records and control points/lines.

Lake bathymetry into the interpolation of the bedrock topography over the regional watershed (Fig. 14). Even with the high quality and regularly spaced North Dakota boreholes, without AEM data, we are unable to guide control lines from section-to-section. The detail of the bedrock surface drops dramatically as we transition out of the survey area and it becomes difficult to delineate and correlate continuous incised valley features between sections. Several North Dakota boreholes encounter bedrock at depths comparable to those observed for the central incised valleys in Manitoba (Fig. 11), but they cannot be convincingly associated with those features. The incised valleys leave the survey area to the south and appear to merge once in North Dakota (Fig. 14). It is unknown if this is an artefact of disparate resolution between the two data regions, or whether this is in fact a real geological transition into an apparent outwash plain style geometry.

4. Discussion Incorporation of the seismic reflection data and the AEM conductivity model for mapping the bedrock topography allows us to define critical features of the hydrogeological system. The high resolution of the seismic reflection data allows us to compensate for variations in the electrical conductivity of the bedrock surface resulting from both real changes in pore water or fracture density, and/or resolution limitations of the AEM data and inversion. Nevertheless, there are still some areas where valley bottoms are identified as occurring below the AEM depth of investigation (Figs. 6 vs. 8). We retain these points in our modelling with the rationale that the seismic data provide the depth control and that regardless of the value of electrical conductivity at the control point, or its veracity, the lateral variation in conductivity which guides the path optimization remains valid.

Fig. 12. a) Target conductivity derived from natural-neighbour interpolation of electrical conductivity projected to observations of bedrock from water well records and control points/ lines. b) Depth of target conductivity filtered for depth of investigation, minimum target depth of 5 m and minimum target conductivity of 40 mS/m. c) Estimated bedrock depth over the AEM survey area derived from natural-neighbour interpolation of target depth, water well records and control points/lines.

G.A. Oldenborger et al. / Journal of Applied Geophysics 128 (2016) 191–201

199

Fig. 13. Manual extension of buried valley geometry to a discharge region along a tributary of the Souris River. Insets show photographs of coarse-grained aquifer material observed on-site. Satellite image © 2015 DigitalGlobe, © 2015 Google.

Despite the high AEM data density, the seismic sections are still localized and the accuracy of our control lines and bedrock mapping likely decreases with distance from the seismic sections. If the conductivity structure of the valley cross section changes significantly between seismic sections, the AEM-based control line may not accurately represent a unique morphological feature (i.e., a valley flank) along its entire path. Furthermore, if a buried valley is not identified in a seismic section, it is either not represented as a control line and it does not inform our interpolation of the bedrock surface, or it must be added through

Fig. 14. Estimated bedrock surface for the regional watershed containing the AEM survey area derived from interpolation of a) water well and borehole records, and b) water well and borehole records and control points/lines. Hatched region indicates the presence of sandstone bedrock outside of the AEM survey area.

manual addition of control points. Subjectivity of our bedrock mapping increases with manual addition of control points to represent features observed in the AEM results, but not captured by seismic sections. Extension of bedrock topography to the entire watershed outside of the AEM survey area is necessary for hydrogeological modelling, but is also subjective and of lower resolution. As such, we have not mapped every buried valley. Furthermore, we are necessarily limited to detecting buried valleys that erode bedrock. Detection of inter-till buried valleys is more complicated due to limited separability of electrical conductivity for the glacial sediments (Oldenborger et al., 2014; Logan et al., 2015) but it may be possible based on seismic reflection data (Pugin et al., 2014; Steinmetz et al., 2015) and electrical resistivity (Oldenborger et al., 2013). With the water wells and seismic data providing information in different regions (and depths) of the survey area, it is difficult to perform a robust test of the accuracy of our results. The limited overlap of the data sets precludes their utilization for cross validation, and omitting seismic lines from the construction of control lines necessarily results in poorer maps. Instead, we perform an analysis of the prediction error distribution associated with each of our estimated bedrock surfaces using a leave-one-out approach applied to the water well data set. Given all water wells with bedrock observations, we sequentially set one aside and perform the bedrock mapping procedure described above. We rebuild the estimated bedrock surfaces (Figs. 11a, 11b and 12c) and compare to the single water well to form a prediction error distribution with the number of samples equal to the number of bedrock water wells. We do not rebuild the control lines since they are largely independent of water wells, but we do rebuild the target conductivity surface (Fig. 12a) and the target depth (Fig. 12b). We recognize that characterisation of the absolute prediction error is not possible with this analysis since bedrock water wells occur predominantly for shallow bedrock and very rarely within the buried valley system, such that our prediction error is specific to these regions which are the least influenced by the control line construct. Nevertheless, our analysis provides some comparison between results. Incorporating control points and lines results in a reduced mean and standard deviation of the prediction error (Fig. 15). The composite bedrock surface exhibits a prediction error distribution essentially equivalent to that of the control data since the addition of target depth is largely responsible for small-scale structure between control

200

G.A. Oldenborger et al. / Journal of Applied Geophysics 128 (2016) 191–201

Fig. 15. Prediction error at bedrock water wells for estimated bedrock surfaces using a) water well records, b) water well records and control points/lines, c) target depth, water well records and control points/lines, and d) target depth. Insets show distribution means, standard deviations and coefficients of skewness.

data. On its own, target depth provides the worst estimation of bedrock depth at water well locations. The large magnitudes of the prediction error reflect both the importance of a single datum in influencing the bedrock surface and the high degree of spatial variability in the water well data set (some part of which is unknown error). In general, we see an over-prediction of bedrock depth for shallow water wells and an under-prediction of bedrock depth for deeper water wells (not shown). Normalized prediction error exhibits under-prediction to intermediate depths and an exponential increase for shallow depths, perhaps indicating a near-surface limitation of the AEM electrical conductivity model. For the Spiritwood buried valley, valley floors are observed to be relatively flat bottomed, but our method would work equally well for V-shaped valleys. Valleys with more complicated shapes may require more control lines and the presence of axially-continuous morphological features on which to base those lines. Similarly, we use constantdepth slices of the AEM conductivity model to guide control lines because we observe little axial elevation change of morphological features over the survey area, and we use an average conductivity between control points to define path cost because changes between corresponding control points (between seismic sections) for the same feature are much less than the changes between control points for different features in cross-section. Buried valleys with significant elevation changes of the thalweg, or significant axial changes in conductivity will require a more sophisticated technique of building the cost raster for path optimization. The spatial resolution of buried valley morphology captured by our geophysical-based bedrock mapping provides information that could not have realistically been obtained from a drilling program. Even the extensive drilling program in North Dakota could not provide the spatial resolution necessary to define the lateral extents, maximum depths and boundaries of the deep incised buried valleys. Given their dimensions, drill hole spacing of approximately 400 m or less would be necessary to define the buried valley depths and thalwegs. It is clear that the geophysical data along with bedrock mapping of the buried valley network has generated a more refined and complete picture and in-depth understanding of buried valley geometry and history (Pugin et al., 2014).

aquifer resources are warranted. Detailed mapping of the bedrock surface is crucial for understanding the hydrogeology of buried valley aquifer systems. For the Spiritwood buried valley aquifer, the limited spatial density of water well bedrock observations precludes detection of the buried valley bedrock topography and renders the water well records alone inadequate for accurate hydrogeological modelling (Fig. 11a). Our integrated use of seismic reflection and airborne electromagnetic data results in superior depiction of the bedrock surface at a level of detail commensurate with requirements for hydrogeological modelling (Fig. 11b). Seismic data provide high resolution information on the bedrock surface, but are limited in number and distributed over a regional area. Our bedrock mapping procedure combines the seismic reflection data, AEM survey results and water well records for a data-driven construction of buried valley morphology. Our method is unique in that the relatively robust interpretation of the bedrock surface is extrapolated from the localized seismic sections via control lines that are constructed from a rule-based analysis of the spatially extensive 3D AEM electrical conductivity model. Our control line dataset is augmented though cognitive interpretation of the AEM model in conjunction with supporting data and field observations for an integrated interpolation of watershed bedrock topography that hydraulically connects the buried valley system to surface waters. The manual augmentation of the control line dataset results in an adaptable balance between the data- and knowledge-driven components of the final bedrock map. Our procedure is semi-automated and reproducible while honouring the available data; new data are readily accommodated. The bedrock map produced over the survey area represents a critical component of geological and hydrogeological modelling and provides increased understanding of the aquifer system through identification of buried valley morphology and its connection to surface hydrology. The bedrock surface also defines the distribution of a regional low yield contact aquifer (the fractured bedrock surface) as well as the underlying regional aquitard. The buried valley network identified in the bedrock map frequently dictates that shape and potential connectivity of aquifers with hydrogeological implications for the groundwater flow dynamics of the Spiritwood buried valley aquifer system.

5. Conclusion

Acknowledgments

In light of the increased concern regarding competing water uses in the arid ecosystem of the Canadian Prairies, efforts to fully characterize

Groundwater Geoscience Program, Natural Resources Canada, contribution number 20150383. Use of ArcGIS® does not constitute

G.A. Oldenborger et al. / Journal of Applied Geophysics 128 (2016) 191–201

endorsement by the authors or Natural Resources Canada. The work presented herein was conducted with assistance from a number of people; we thank K. Brewer, T. Cartwright, H. Crow, A. Calderhead, M. Douma, L. Katz, P. Keating, R. Knight, B. Medioli, W. Miles, J. Oliver and M. White. We thank N. Benoit, C. Brandes and one anonymous reviewer for constructive critical review. Data for monitoring wells, water wells, and knowledge regarding the Spiritwood buried valley aquifer were provided by B. Betcher and G. Matile of Manitoba Conservation and Water Stewardship. Manitoba water well records from Manitoba Conservation and Water Stewardship are freely available through the Canadian Groundwater Information Network (http://gin.gw-info.net). North Dakota water wells were obtained from the North Dakota State Water Commission (http://www.swc.nd.gov/). Landsat images and AEM data are freely available from Natural Resources Canada GeoGratis (http:// geogratis.gc.ca/) and Geoscience Data Repository (http://gdr.agg. nrcan.gc.ca). References Abraham, J.D., Cannia, J.C., Peterson, S.M., Smith, B.D., Minsley, B.J., Bedrosian, P.A., 2010. Using airborne geophysical surveys to improve groundwater resource management models. Proc. Symp. Appl. Geophys. Environ. Eng. Probl. 309–314. Ahmad, J., Schmitt, D.R., Rokosh, C.D., Pawlowicz, J.G., 2009. High-resolution seismic and resistivity profiling of a buried quaternary subglacial valley: Northern Alberta, Canada. GSA Bull. 121, 1570–1583. Auken, E., Christiansen, A.V., Jacobsen, L.H., Sørensen, K.I., 2008. A resolution study of buried valleys using laterally constrained inversion of TEM data. J. Appl. Geophys. 65, 10–20. Betcher, R.N., Matille, G., Keller, G., 2005. Yes Virginia, there are buried valley aquifers in Manitoba. Proceedings of the 58th Canadian Geotechnical Conference, 6E-519. Bluemle, J.P., 1983. Geologic and topographic bedrock map of North Dakota. North Dakota Geological Survey, Miscellaneous Map 25. Christensen, C.W., Pfaffhuber, A.A., Anschütz, H., Smaavik, T.F., 2015. Combining airborne electromagnetic and geotechnical data for automated depth to bedrock tracking. J. Appl. Geophys. 119, 178–191. Christiansen, A.V., Auken, E., 2012. A global measure for depth of investigation. Geophysics 77, 171–177. Christiansen, A.V., Auken, E., Viezzoli, A., 2011. Quantification of modeling errors in airborne TEM caused by inaccurate system description. Geophysics 76, F43–F52. Crow, H.L., Brewer, K.D., Pugin, A.J.-M., Russell, H.A.J., 2012a. Downhole geophysical data from boreholes along the Spiritwood buried valley aquifer near Cartwright, Killarney, and southeast of Brandon, Manitoba. Geological Survey of Canada, Open File 7080. Crow, H.L., Knight, R.D., Medioli, B.E., Hinton, M.J., Plourde, A., Pugin, A.J.-M., Brewer, K.D., Russell, H.A.J., Sharpe, D.R., 2012b. Geological, hydrogeological, geophysical, and geochemistry data from a cored borehole in the Spiritwood buried valley, southwest Manitoba. Geological Survey of Canada, Open File 7079. Cummings, D.I., Russell, H.A.J., Sharpe, D.R., 2012. Buried valley aquifers in the Canadian prairies: geology, hydrogeology, and origin. Can. J. Earth Sci. 49, 987–1004. Day-Lewis, F.D., Singha, K., Binley, A.M., 2005. Applying petrophysical models to radar travel time and electrical resistivity tomograms: resolution-dependent limitations. J. Geophys. Res. 110, B08206. Douglas, D.H., 1994. Least-cost path in GIS using an accumulated cost surface and slopelines. Cartographica 31, 37–51. Foged, N., Auken, E., Christiansen, A.V., Sørensen, K.I., 2013. Test-site calibration and validation of airborne and ground-based TEM systems. Geophysics 78, E95–E106. Gabriel, G., Kirsch, R., Siemon, B., Wiederhold, H., 2003. Geophysical investigation of buried Pleistocene subglacial valleys in Northern Germany. J. Appl. Geophys. 53, 159–180. Goldman, M., Tabarovsky, L., Rabinovich, M., 1994. On the influence of 3-D structures in the interpretation of transient electromagnetic sounding data. Geophysics 59, 889–901. Grasby, S.E., Betcher, R.N., Maathius, H., Wozniak, P.R.J., 2014. Plains Region. In: Rivera, A. (Ed.), Canada's Groundwater Resources. Fitzhenry and Whiteside, pp. 358–413. Høyer, A.-S., Lykke-Andersen, H., Jørgensen, F., Auken, E., 2011. Combined interpretation of SkyTEM and high-resolution seismic data. Phys. Chem. Earth 36, 1386–1397. Høyer, A.-S., Jørgensen, F., Sandersen, P.B.E., Viezzoli, A., Møller, I., 2015. 3D geological modelling of a complex buried-valley network delineated from borehole and AEM data. J. Appl. Geophys. 122, 94–102. Huuse, M., Lykke-Andersen, H., 2000. Overdeepened quaternary valleys in the eastern Danish North Sea: morphology and origin. Quat. Sci. Rev. 19, 1233–1253. Jørgensen, F., Sandersen, P.B.E., 2009. Buried Valley mapping in Denmark: evaluating mapping method constraints and the importance of data density. Z. Dtsch. Ges. Geowiss. 160, 211–223. Jørgensen, F., Lykke-Andersen, H., Sandersen, P.B.E., Auken, E., Normark, E., 2003. Geophysical investigations of buried quaternary valleys in Denmark: an integrated application of transient electromagnetic soundings, reflection seismic surveys and exploratory drillings. J. Appl. Geophys. 53, 215–228. Jørgensen, F., Scheer, W., Thomsen, S., Sonnenborg, T.O., Hinsby, K., Wiederhold, H., Schamper, C., Burschil, T., Roth, B., Kirsch, R., Auken, E., 2012. Transboundary geophysical mapping of geological elements and salinity distribution critical for the assessment of future sea water intrusion in response to sea level rise. Hydrol. Earth Syst. Sci. 16, 1845–1862.

201

Jørgensen, F., Møller, R.R., Nebel, L., Jensen, N.-P., Christiansen, A.V., Sandersen, P.B.E., 2013. A method for cognitive 3D geological voxel modelling of AEM data. Bull. Eng. Geol. Environ. 72, 421–432. Kehew, A.E., Boettger, W.M., 1986. Depositional environments of Buried-Valley aquifers in North Dakota. Ground Water 24, 728–734. Kehew, A.E., Piotrowski, J.A., Jørgensen, F., 2012. Tunnel valleys: concepts and controversies — A review. Earth Sci. Rev. 113, 33–58. Logan, C.E., Hinton, M.J., Sharpe, D.R., Oldenborger, G.A., Russell, H.A.J., Pugin, A.J.M., 2015. Spiritwood buried valley 3D geological modelling – Part of a multidisciplinary aquifer characterization workflow. Geological Survey of Canada, Open File 7866. Maathuis, H., Thorleison, L.H., 2000. Potential impact of climate change on prairie groundwater supplies: Review of current knowledge. Saskatchewan Research Council, Publication No. 11304-2E00. Nicolas, M.P.B., Matile, G.L.D., Keller, G.R., Bamburak, J.D., 2010. Phanerozoic geology of southern Manitoba. Manitoba Geological Survey, Stratigraphic Map SM2010–1. Oldenborger, G.A., Brewer, K.B., 2014. Time-domain electromagnetic data for the spiritwood valley aquifer, Manitoba. Geological Survey of Canada, Open File 7593. Oldenborger, G.A., Pugin, A.J.-M., Pullan, S.E., 2013. Airborne time-domain electromagnetics, electrical resistivity and seismic reflection for regional three-dimensional mapping and characterization of the Spiritwood valley aquifer, Manitoba, Canada. Near Surf. Geophys. 11, 63–74. Oldenborger, G.A., Logan, C.E., Hinton, M.J., Sapia, V., Pugin, A.J.-M., Sharpe, D.R., Calderhead, A.I., Russell, H.A.J., 2014. 3D hydrogeological model building using airborne electromagnetic data. 20th European Meeting of Environmental and Engineering Geophysics, Tu-PA1–07. Parks, K.P., Bentley, L.R., 1996. Derivative-assisted evaluation of well yields in a heterogeneous aquifer. Can. Geotech. J. 33, 458–469. Paulson, Q.F., 1983. Guide to North Dakota's groundwater resources. United States Geological Survey, Water-Supply Paper 2236. Piotrowski, J.A., 1997. Subglacial hydrology in North-Western Germany during the last glaciation: groundwater flow, tunnel valleys, and hydrological cycles. Quat. Sci. Rev. 16, 169–185. Pryet, A., Ramm, J., Chilès, J.-P., Auken, E., Deffontaines, B., Violette, S., 2011. 3D resistivity gridding of large AEM datasets: A step toward enhanced geological interpretation. J. Appl. Geophys. 75, 277–283. Pugin, A.J.-M., Pullan, S.E., Hunter, J.A., Oldenborger, G.A., 2009a. Hydrogeological prospecting using P- and S-wave landstreamer seismic reflection methods. Near Surf. Geophys. 7, 315–327. Pugin, A.J.-M., Pullan, S.E., Hunter, J.A., 2009b. Multicomponent high resolution seismic reflection profiling. Lead. Edge 28, 936–945. Pugin, A.J.-M., Oldenborger, G.A., Cummings, D.I., Russell, H.A.J., Sharpe, D.R., 2014. Architecture of buried valleys in glaciated Canadian prairie regions based on high resolution geophysical data. Quat. Sci. Rev. 86, 13–23. Randich, P.G., Kuzniar, R.L., 1984. Geology of Towner County, North Dakota. North Dakota State Water Commission, County Groundwater Studies 36, Part III. Reninger, P.-A., Martelet, G., Lasseur, E., Beccaletto, L., Deparis, J., Perrin, J., Chen, Y., 2014. Geological environment of karst within chalk using airborne time domain electromagnetic data cross-interpreted with boreholes. J. Appl. Geophys. 106, 173–186. Russell, H.A.J., Hinton, M.J., van der Kamp, G., Sharpe, D., 2004. An overview of the architecture, sedimentology and hydrogeology of buried-valley aquifers in Canada. Proceedings of the 57th Canadian Geotechnical Conference, pp. 26–33. Sapia, V., Oldenborger, G.A., Viezzoli, A., Marchetti, M., 2014a. Incorporating ancillary data into the inversion of airborne time-domain electromagnetic data for hydrogeological applications. J. Appl. Geophys. 104, 35–43. Sapia, V., Viezzoli, A., Jørgensen, F., Oldenborger, G.A., Marchetti, M., 2014b. The impact on geological and hydrogeological mapping results of moving from ground to airborne TEM. J. Environ. Eng. Geophys. 19, 53–66. Sapia, V., Oldenborger, G.A., Jørgensen, F., Pugin, A.J.-M., Marchetti, M., Viezzoli, A., 2015. 3D modeling of buried valley geology using airborne electromagnetic data. Interpretation 3, SAC9–SAC22. Shaver, R.B., Pusc, S.W., 1992. Hydraulic barriers in Pleistocene buried-valley aquifers. Ground Water 30, 21–28. Sie, D., Little, J., 1976. Groundwater Availability Map Series, Brandon Area (62-G). Manitoba Natural Resources, Water Resources Branch. Siemon, B., Christiansen, A.V., Auken, E., 2009. A review of helicopter-borne electromagnetic methods for groundwater exploration. Near Surf. Geophys. 7, 629–646. Sørensen, K.I., Auken, E., 2004. SkyTEM – A new high-resolution helicopter transient electromagnetic system. Explor. Geophys. 35, 191–199. Steinmetz, D., Winsemann, J., Brandes, C., Siemon, B., Ullmann, A., Wiederhold, H., Meyer, U., 2015. Towards an improved geological interpretation of airborne electromagnetic data: A case study from the Cuxhaven tunnel valley and its Neogene host sediments (Northwest Germany). Neth. J. Geosci. 94, 201–227. Steuer, A., Siemon, B., Auken, E., 2009. A comparison of helicopter-borne electromagnetics in frequency- and time-domain at the Cuxhaven valley in Northern Germany. J. Appl. Geophys. 67, 194–205. Stewart, M.A., Lonergan, L., 2011. Seven glacial cycles in the middle-late Pleistocene of northwest Europe: geomorphic evidence from buried tunnel valleys. Geology 39, 283–286. van der Kamp, G., Maathius, H., 2012. The unusual and large drawdown response of buried-valley aquifers to pumping. Ground Water 50, 207–215. Viezzoli, A., Christiansen, A.V., Auken, E., Sorensen, K., 2008. Quasi-3D modeling of airborne TEM data by spatially constrained inversion. Geophysics 73, F105–F113. Winter, T.C., Benson, R.D., Engberg, R.A., Wiche, G.J., Emerson, D.G., Crosby, O.A., Miller, J.E., 1984. Synopsis of ground-water and surface-water resources of North Dakota. United States Geological Survey, Open File Report 84–732.