Over the last decade, Airborne Electromagnetics (AEM) has become a widespread tool for groundwater applications. Besides the demand of acquiring good quality AEM data, there are two other fundamental steps to obtain a robust geological and hydrogeological model: accurate processing/inversions and advanced interpretation.

We show the results achieved adopting this approach, on a subset of a large AEM survey within the Peace Project, in British Columbia, Canada (www.geosciencebc.com/s/PeaceProject.asp). This project is a collaborative effort involving Geoscience BC, the Ministry of Forests, Lands and Natural Resource Operations, the Ministry of Environment, the BC Oil and Gas Commission, the Ministry of Natural Gas Development, BC Oil and Gas Research and Innovation Society (BC OGRIS), Progress Energy Canada Ltd. (now Petronas Canada), and ConocoPhillips Canada, with additional support from the Peace River Regional District and the Canadian Association of Petroleum Producers.

This paper (largely based on Viezzoli et al., 2018) is the natural continuation of the Best et al. (2017) December 17 RECORDER issue, in which they presented the results of the integration between AEM inversion results and gamma logs.

AEM data were acquired by SkyTEM Aps on behalf of Geoscience BC. Approximately half of the survey was later processed and inverted by Aarhus Geophysics Aps. Finally, GEUS (the Geological Survey of Denmark and Greenland) carried out a hydrogeological interpretation on a yet smaller subset. The whole prospect covers approximately 21,000 line km and we present the hydrogeological modeling for a subset of the survey area, located in the north-west (Figure 1) that corresponds to about 2000 km2 (i.e. approximately 4500 line km).

From a hydrogeological point of view, the main target is the detection of buried valleys, since they hold the most promising potential for sustainable groundwater production in the area.

Fig. 01
Figure 1. Right. Canada. Adapted from www.google.it/maps. Centre. The green polygon is the Peace Project Main Area. Adapted from www.geosciencebc.com. Left. The red polygon outlines the area of the dataset located in the N-W part of the Peace River Project. The black polygon (Area 1), within the red one, outlines the area involved for the detailed 3D interpretation and modelling. Adapted from Petrel Robertson, 2015.

Geological setting

Restricting our attention to the geological formations present within the expected penetration depths of the SkyTEM method (approximately 2 to 400 m in this geology), we identify, from older to younger: the Nordegg Formation (calcareous mudstones), the Buckinghorse Formation (shales, silty mudstones and siltstones), the Sikanni Formation (sandstones and siltstones), the Sully Formation (shales and siltstones) and the Dunvegan Formation (sandstone with subordinate conglomerates). The potential bedrock aquifers are represented by the sandstone units within the Sikanni and Dunvegan, together with the conglomerate lenses within the latter.

The Quaternary sequence is composed of glacio-fluvial deposits (sands, pebbles and gravels), glacio-lacustrine deposits (silty clays or clayey silts), tills, colluvial (sandy and clayey silts) and fluvial units (mainly sands and gravels). The sequence is generally quite thin outside the paleovalleys. The paleovalleys show thick sequences of often coarse-grained infill sediments that might constitute promising aquifers. However, we must pay great attention to the relationships between these shallower aquifers and any potential contaminant coming from the surface, as they are more vulnerable from an environmental point of view. At the same time, many of the paleovalleys are situated directly beneath the modern riverbed. Groundwater extraction from these areas might therefore influence the rivers markedly. Hence, correct hydrogeological modeling is mandatory to responsibly manage the groundwater resources.

The table below reports the expected resistivity ranges (from downhole resistivity logs) of some of the main formations encountered in the area.

Formation Resistivity
Glaciofluvial deposits ~ 100 Ωm
Glaciolacustrine deposits ~ 10 Ωm
Dunvegan Formation ~ 50 Ωm
Sully Formation ~ 15 Ωm
Sikanni Formation ~ 30 Ωm


Figure 2 shows the schematic of the overall workflow adopted in this project. Notice how it is not unidirectional, but iterative, with preliminary results from later steps feeding back into previous steps, for fine tuning. Especially important is the close cooperation between geophysicists and geologists. The goal is to obtain a geological model consistent, at once, with all sources of information (and their respective uncertainties): geophysical data, geophysical models, prior geological information (e.g., drilling), and conceptual geological model.

Fig. 02
Figure 2. Schematic workflow applied. Notice the feedback loops disseminated throughout.

Data acquisition, processing and inversions

SkyTEM is a time-domain helicopter electromagnetic system designed for hydrogeophysical, environmental and mineral investigations. The SkyTEM system is composed of a transmitter (12 turns, 337 m2 eight-sided loop) and a z-component receiver (placed approximately 2 m above the frame). The system makes use of both a Low dipole Moment (LM) and High dipole Moment (HM) that yield both near surface and deep resolution.

Data are subject to QA/QC procedures at the end of each flight. This survey (21000-line km) was finalized in approximately 6 weeks, with 2 SkyTEM systems. After acquisition, the data is processed. The processing largely follows Auken at el. (2009).

Soundings were averaged over 1.4 s of flying, corresponding to approximately 30 m separation between soundings. The resulting data density greatly helps the geological-hydrogeological interpretation. Figure 3 shows an example of a TEM sounding, shown as apparent resistivity vs time, with different time gates displaying different noise levels for the low and high moment.

Fig. 03
Figure 3. Example of averaged sounding (transformed into late time apparent resistivity), with different time gates displaying different noise levels, for low (left curve) and high (right curve) moment.

The AEM and navigation data are first processed to eliminate cultural artefacts and improve signal to noise ratio. The fundamental task, carried out utilizing a combination of automated and manual procedures, prevents artefacts in the resistivity and hydrogeological models, as shown by Viezzoli et al. (2013).

Final inversions are carried out using the quasi 3-D Spatially Constrained Inversion (SCI) (Viezzoli et al., 2008). Using constraints across model parameters, which are tied together with a spatially dependent covariance, the SCI returns the spatial variability expected in the local geological setting. The flight altitude is included as an inversion parameter with a prior value calculated from the tilt corrected laser altitudes. The Depth of Investigation (DOI, Christiansen and Auken, 2012, based on the sensitivity to the model parameters, was calculated as the last step of the inversion. It is a crucial ancillary QC and is most useful when inspecting and interpreting results at depth.

Geophysical Results

The inversion results are compared with the existing boreholes. The resistivity cross-sections of Figure 4 show a couple of examples, with the stratigraphic data represented essentially by the thickness of Quaternary sediments and depth to and type of bedrock.

Fig. 04
Figure 4. Vertical resistivity cross sections compared to borehole data (depth to and – color coded- type of bedrock), overlain with preliminary qualitative interpretation of main structures. Distance (m) from drillhole location is shown above the boreholes. Models are shaded with DOI. The associated data misfits are shown at the bottom of the profiles (black line, to be read against right axis).

Even though the resistivity values for the different lithological formations estimated from the logs only show small differences (Table 1), the lithological formations seem to be easily distinguished in the SkyTEM’s derived resistivities. Thus, the Sully and Buckinghorse formations show low resistivities due to the presence of shale, whereas Sikanni and Dunvegan show high resistivities. The capability to distinguish different bedrock formations allowed delineating a complex structural framework (Figure 4).

Depth to bedrock is generally well resolved both shallow and deep. See for example the good agreement between the boreholes around 9000 m for the top cross-section and at 7000 m for the bottom one. The Quaternary units can display a sharp conductive feature where glaciolacustrine deposits prevail, or a rather resistive response, linked to buried valleys filled with gravels and sands.

Figure 5 displays the correlation between outcropping geology and shallow resistivity (0-5 m depth), confirming the electrical response of the geological formation: note the elongated resistive deposits associated with the fluvio-glacial or the alluvial deposits within the Halfway River and its tributary Chowade River. On the contrary, the finer glacial till is characterized by lower resistivity. In the northwest corner of the map it is possible to recognize the typical resistive response of the Nordegg Formation that is the only bedrock unit outcropping in this area and that is marked as “bedrock” in this geological map.

Fig. 05
Figure 5. Resistivity slice map at depth of 0-5 m (left) shown next to the surface geology map (right).

Geological interpretation

Interpretation of the AEM data reveals a highly complex structural setting, where the bedrock formations show large-scale folds and thrust structures (Figure 4). Additional structural information can be derived from the results (see Jørgensen et al., 2017, 2018), but this paper will mainly focus on the mapping and modelling of paleovalleys, due to their groundwater extraction potential.

A number of paleovalleys have been mapped in the study area (Figure 6a). The paleovalleys are recognized in the AEM data as elongate resistivity structures that broaden upwards. Most of the paleovalleys are situated along the modern valleys stretches, but others are detected in hilly areas. The buried valleys occur in at least two generations. The first generation has overdeepened sections along their thalwegs and may therefore be interpreted as subglacially formed tunnel valleys that later have been infilled with younger sediments. This type of buried paleovalleys are frequently found within all formerly glaciated areas (Jørgensen and Sandersen, 2006; Kehew et al., 2012). The formation of the second generation of buried valleys is connected to the recent fluvial erosion and deposition of pebbles and gravels in the riverbed. The second generation of valleys are limited in depth and are found beneath many stretches of modern rivers, such as the Halfway River and its tributaries (Figure 5).

A profile along the southern paleovalley is shown in Figure 6b. The valley shows an undulating bottom profile that is characteristic of subglacially formed valleys. The valley is up to 100 m deep at its deepest depressions and shows high-medium resistivities at depth. In the northern part of the valley (profile distance 0 – 20 km), very high resistivities are seen in the upper part. Here, the older valley has been cut by the younger, fluvial erosion and deposition of coarse-grained gravels and pebbles. Low resistivity deposits are recognized on top of the paleovalley structure in many places. These low resistivity deposits are glaciolacustrine clay that covers large parts of the modern river valley (Figure 5b, Petrel Robertson, 2015).

Fig. 06
Figure 6. Paleovalleys interpreted in the study area (note: due to marked topography in the area, the resistivity signatures are apparent at different horizontal slices), a) horizontal slice through the 3D resistivity grid at elevation 725 m, b) profile through the 3D resistivity grid following the thalweg of the southern paleovalley, shown from north (left) to south (right). The DOI is shown as a thin grey line. Note the undulating bottom profile.

The second long paleovalley is found beneath the northernmost stretches of the Halfway River Valley (Figure 7). It can be followed until it crosses Pink Mountain in the north. This valley is up to about 75 m deep in its northern part (Figure 7b, profile distance c. 10 – 15 km). These valleys have subsequently been cut by modern river streams along large sections of their length and in these sections it is difficult to distinguish between the modern river sediments and the older paleovalley fill. The uncertainty of the interpretations is thus relatively high here.

Fig. 07
Figure 7. Paleovalleys interpreted in the study area (note: due to marked topography in the area, the resistivity signatures are apparent at different horizontal slices), a) horizontal slice through the 3D resistivity grid at elevation 850 m, b) profile through the 3D resistivity grid following the thalweg of the northern paleovalley, shown from north (left) to south (right).

Geological modelling

A geological model was constructed for a subset of the area. The model is constructed based on a cognitive approach, in which all available data (surface geological map, stratigraphic boreholes, AEM data, terrain, etc.) and the overall conceptual geological knowledge are combined and translated into a manually interpreted 3D voxel model (see also Høyer et al., 2015; Jørgensen et al., 2013).

The 3D geological model was constructed using the geological software modelling tool ‘Geoscene3D’ from I-GIS (2018). The model area (Area 1 in Figure 1) covers a section of the long paleovalley that follows the southern part of the Halfway River Valley (profile distance 10 km to 18 km in Figure 6b). In the following, this valley is named ‘Buried Valley 1’. Figure 8 shows a close-up of the surface geology map in the model area together with the outline of the interpreted buried valleys and the position of the profiles in Figures 10.

The geological modelling is performed using a combination of layer- and voxel modelling tools using a similar approach as the geological modelling of buried valleys in Høyer et al. (2015) and of Sapia et al. (2015). The bottom of each geological unit is modelled using interpretation points that have been kriged into surface layers. The bottom of four overall stratigraphic boundaries have been modelled; Buckinghorse, Sikanni, Sully and Dunvegan Formations. In the valley, the bottom of six Quaternary boundaries have been modelled; Buried Valley 1, Buried Valley 2, Buried Valley 3, the Meltwater Plain, the Modern Valley Fill and the Glaciolacustrine Deposits.

Fig. 08
Figure 8. The model area (grey) shown on top of the surface geology map by Petrel Robertson, 2015. Outline of the interpreted buried valleys are shown with blue.

Subsequently, a voxel model is constructed by using the tools described in Jørgensen et al. (2013). The final voxel model is discretized with 100 m x 100 m laterally and 5 m vertically and covers the elevation interval from 300 m a.s.l. to the terrain surface (up to 1210 m a.s.l.), resulting in about 7.7 million voxels. The voxel model is constructed by selecting and populating voxels within volumes delineated by the modelled layer surfaces and defined regions. Overall, the units correspond to the ones defined by the layer surfaces, but two of the units (‘Buried Valley 1’ and the ‘Dunvegan Fm’) are further subdivided into lithological facies based on the resistivity values within the units (with 60 ohm-m as a cut-off value). The stratigraphic formations are known to be heterogeneous, but due to the greater burial depth of these formations, the SkyTEM data are not able to resolve the small lithological variations, and resistivity values are therefore not used for lithological characterization within these units.

Fig. 09
Figure 9. 3D view of the model results seen from southeast. 5 x vertical exaggeration, a) 3D view of the voxels, b) E-W and N-S slices through the model shown together with the layer surface outlining Buried Valley 1. The bedrock formations are shown in grey colors (legend embedded on figure a) and Quaternary deposits are shown with colors (legend at bottom).

The resulting voxel model consists of 13 categories, which are shown together with the 3D view of the modelling results in Figure 9. The Quaternary deposits have only been modelled within the modern valley structure (Figure 9), since the remaining part of the area generally shows a thin cover of Quaternary sediments (Figure 8), with thicknesses below the vertical discretization in the voxel grid (<5 m). In the modern valley structure, the Quaternary deposits have thicknesses between 20 m to approximately 200 m, where the buried valleys are deepest.

Figure 10 shows the resistivity data together with the model results along the northernmost flightlines within the study area (Figure 8). In the northern part of the model area (Figure 10a) a synclinal structure is recognized east of the valley, whereas a thrust fault structure is recognized below the Buckinghorse Formation in the southern part (Figure 10b). The Sully and Dunvegan Formations are only present in the southeastern part of the model.

Fig. 10
Figure 10. a) SkyTEM resistivities and the DOI along the northernmost flightline within the 3D model area (for location see Figure 7), b) geological model results along the same profile. 5x vertical exaggeration.

Three buried valleys are modelled within the model area. The largest is the one named ‘Buried Valley 1’. This corresponds to the southernmost long tunnel valley (Figure 5) that is located directly below the modern river in the southern part of the model (Figure 9a,b). In the north, it is located below an erosional remnant of glaciolacustrine deposits (Figure 10, profile distance 7700 m). The lower part of the valley is mainly filled with clayey deposits, whereas the upper part has coarse-grained infill. Close to the hill-side in the eastern side of the modern valley, Buried Valley 2 is located and filled with mainly coarse-grained deposits. Buried Valley 3 is located in the northwestern part of the model area (Figure 10) and is filled with clayey deposits.

All three buried valleys are interpreted to be first generation tunnel valleys that are older than the meltwater plain. The coarse-grained deposits from this unit therefore cut through the buried valley structures (Figures 9 and 10). Glaciolacustrine clays are modelled on top of the Meltwater Plain, and in the north, a marked erosional remnant of these glaciolacustrine deposits forms a small hill structure in the middle of the modern valley (Figure 9a). The youngest deposits in the geological model are the sand and gravel deposits from the modern riverbed, which is cut into the underlying deposits (e.g. Figure 10, profile distance 8500-10,000 m).


A coherent hydrogeological model of a subset of the Peace River SkyTEM survey was obtained through advanced processing, modelling and geological interpretation of the AEM data. Resistivity models provided a new, detailed insight into the local geology, improving the structural and lithological framework. Moreover, thanks to the sharp resistivity contrast, it was possible to map buried valleys, filled with permeable units that could represent important local aquifers.

The highest resistivities are found within the modern riverbed (>100 ohm-m), thus indicating coarse-grained sediments without much clay. The most interesting sites for groundwater extraction could therefore be within the sand-filled parts of the buried valley structures, or in the meltwater plain, where those units are not in direct hydraulic contact with the modern riverbed.

Future work

As Morgan et al. (2017) summarized in their final report, recommendations for future work may include: extending the geological model to the whole area; improving the understanding of the groundwater resource potential of bedrock aquifers; additional geophysical surveys (e.g. seismic, ground-based EM, a SkyTEM survey flown south of the Peace River); additional drilling, including pumping tests, to further geological understanding.



The authors acknowledge the in-kind and funding support that the partners of the Peace Project have made to this study.


About the Author(s)

Melvyn Best has more than 40 years of geophysical experience in industry and government. He worked for Shell and Teknica Resource Development Limited in Calgary, Shell Oil in Houston, and Royal Dutch Shell in the Netherlands. He then joined the Geological Survey of Canada in Dartmouth, NS and later in Sidney, BC as director of the Pacific Geoscience Centre. In 1997, he took early retirement from the GSC and started a geophysical consulting practice in oil and gas exploration and production as well as ground water and environmental studies. He is a past adjunct professor at the University of Victoria and at the University of Calgary. Mel has published over 50 journal and GSC papers.

Anne-Sophie Høyer is a Senior Researcher at the Geological Survey of Denmark and Greenland. She is mainly working with geological interpretation of geophysical data and 3D geological modelling.

Flemming Jørgensen is Ph.D. and Chief Consultant in the Central Denmark Region. He has a 20-years of research-based experience with integrated geological modelling of hydrogeological and geophysical data, in particular, electromagnetic data.

Antonio Menghini is working as a Senior Geophysicist for Aarhus Geofisica s.r.l. for the processing, inversion and interpretation of airborne electromagnetic (AEM) data. He is a member of SEG, EAGE, ASEG and EEGS, and Associate Editor of the Journal of Environmental and Engineering Geophysics.

Carlos Salas is a professional geologist with over thirty years of wide-ranging energy sector experience including senior positions and directorships in public and private companies. Initially schooled as a geologist specializing in ancient rivers, he began his professional career with a major oil and gas company, then as an explorationist in several ‘start-ups’ where he honed his entrepreneurial and business skills. He then went on to work at a medium-sized consultancy where he managed a large group of professionals. In his current role of Executive Vice President and Chief Scientific Officer at Geoscience BC, Carlos provides, technical leadership to initiatives promoting informed resource management decision-making through earth science.

Andrea Viezzoli earned his Ph.D. in Geophysics from Monash University, Australia. Author of dozens of publications, he is interested in all aspects of airborne electromagnetics (AEM), with special focus on hydrogeological applications and on induced polarization (IP) effects in AEM data.


Auken, E., A.V. Christiansen, J.H. Westergaard, C. Kirkegaard, N. Foged, and A. Viezzoli. 2009. An integrated processing scheme for high-resolution airborne electromagnetic surveys, the SkyTEM system. Exploration Geophysics 40: 184-192.

Best, M., Levson, V, and Salas, C., 2017, Integrating airborne electromagnetic and gamma log data for regional aquifer mapping: Montney play area, NE British Columbia, CSEG RECORDER, December, 2017, pp. 28-33.

Christiansen, A. V., and E. Auken, 2012, A global measure for depth of investigation, Geophysics, 77, 4, WB171-WB177.

Hayes, BJR., Levson, V, Carey, J, and Mykola, Y., 2016, Interpretation of Quateranry sediments and depth to bedrock, Peace River area, northeastern British Columbia: project update; in Geoscience BC Summary of Activities 2015, Geoscience BC, Report 2016-1, pp. 61-68.

Høyer, A.-S., Jørgensen, F., Sandersen, P. B. E., Viezzoli, A., and Møller, I., 2015, 3D geological modelling of a complex buried-valley network delineated from borehole and AEM data: Journal of Applied Geophysics, v. 122, pp. 94-102.

I-GIS, 2018, Geoscene3D: Aarhus, Denmark.

Jørgensen, F., Menghini, A., Vignoli, G., Viezzoli, A., Salas, C., Best, M. E., and Pedersen, S. A. S., 2017, Structural geology interpreted from AEM data: Folded terrain at the foothills of Rotial recky Mountains, British Columbia. Near Surface Geoscience, 3-7 September, 2017, Malmö.

Jørgensen, F., and Sandersen, P. B. E., 2006, Buried and open tunnel valleys in Denmark-erosion beneath multiple ice sheets: Quaternary Science Reviews, v. 25, no. 11-12, pp. 1339-1363.

Jørgensen, F., Møller, R.R., Nebel, L., Jensen, N.-P., Christiansen A.V. and Sandersen, P.B.E. 2013: A method for cognitive 3D geological voxel modelling of AEM data. Bulletin of Engineering Geology and the Environment. Vol. 72, 3, 421- 432. DOI: 10.1007/s10064-013-0487-2.

Jørgensen, F., Menghini, A., Kallesøe, A.J., Vignoli, G., Viezzoli, A. & Pedersen, S.A.S. 2016: Structural geology of folded terrain in the Rocky Mountains' foothills interpreted from SkyTEM – A preliminary study of SkyTEM data collected in the Peace region, NE British Columbia, Canada. Danmarks og Grønlands Geologiske Undersøgelse Rapport 2016/34. 21 p.

Kehew, A. E., Piotrowski, J. A., and Jørgensen, F., 2012, Tunnel valleys: Concepts and controversies – A review: Earth-science Reviews, v. 113, p. 33-58.

Morgan, S.E., Allen, D.M., Kirste. D., and Salas, C.J., 2017, Investigating the role of buried valley aquifer systems in the regional hydrogeology of the Peace River Region, Northeastern British Columbia. Geoscience BC Report, 2017-1, p. 63-68.

Sapia, V., Oldenborger, G. A., Jørgensen, F., Pugin, A. J. M., Marchetti, M., and Viezzoli, A., 2015, 3D modelling of buried valley geology using airborne electromagnetic data: Interpretation, v. 3, no. 4, p. SAC9-SAC22.

Viezzoli, A., Christiansen, A. V., Auken, E., and Sørensen, K., 2008, Quasi-3D modeling of airborne TEM data by spatially constrained inversion: Geophysics, 73(3), F105-F113.

Viezzoli, A., Jørgensen, F., Sorensen, C., 2013, Flawed processing of airborne EM data affecting hydrogeological interpretation. Ground Water. 51(2)191-202.

Viezzoli A., Menghini A., Jørgensen F., Høyer A.S. 2018. Processing and inversion of SkyTEM data leading to a hydrogeological interpretation of the Peace River North-Western Area. Report 2018-06, 2-27. Freely downloaded from Geoscience BC website: http://www.geosciencebc.com/s/PeaceProject.asp


Join the Conversation

Interested in starting, or contributing to a conversation about an article or issue of the RECORDER? Join our CSEG LinkedIn Group.

Share This Article