Compiling a national resistivity atlas of Denmark based on airborne and ground-based transient electromagnetic data

Article history: Received 19 January 2016 Received in revised form 8 September 2016 Accepted 12 September 2016 Available online 13 September 2016 We present a large-scale study of the petrophysical relationship of resistivities obtained from densely sampled ground-based and airborne transient electromagnetic surveys and lithological information from boreholes. The overriding aim of this study is to develop a framework for examining the resistivity-lithology relationship in a statistical manner and apply this framework to gain a better description of the large-scale resistivity structures of the subsurface. In Denmark very large and extensive datasets are available through the national geophysical and borehole databases, GERDA and JUPITER respectively. In a 10 by 10 km grid, these data are compiled into histograms of resistivity versus lithology. To do this, the geophysical data are interpolated to the position of the boreholes, which allows for a lithological categorization of the interpolated resistivity values, yielding different histograms for a set of desired lithological categories. By applying the proposed algorithm to all available boreholes and airborne and ground-based transient electromagnetic data we build nation-wide maps of the resistivity-lithology relationships in Denmark. The presented Resistivity Atlas reveals varying patterns in the large-scale resistivity-lithology relations, reflecting geological details such as available source material for tills. The resistivity maps also reveal a clear ambiguity in the resistivity values for different lithologies. The Resistivity Atlas is highly useful when geophysical data are to be used for geological or hydrological modeling. © 2016 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
Aquifer mapping and vulnerability assessment is important for securing clean drinking water globally. Geophysical methods, in particular electric and electromagnetic methods, with resistivity as the main petrophysical property, have proven to be important tools for mapping groundwater resources and their vulnerability, when combined with hydrological and lithological information from boreholes (e.g. Slater, 2007;Siemon et al., 2009). Boreholes generally have a low spatial density, but are complemented by spatially dense geophysical data, e.g. Sandersen et al. (2009) and Nobes (1996). Combining lithological logs from boreholes and resistivity data from ground-based and airborne electromagnetic soundings provides the area coverage and detailed lithological information required for creating 3D hydrogeological models of the subsurface with adequate accuracy for the purpose of the models and the related 3D groundwater models (e.g. Thomsen et al., 2004;Sandersen et al., 2009;Bosch et al., 2009;Jørgensen et al., 2012;Oldenborger et al., 2013). When building hydrogeological models, knowledge or assumptions regarding the relation between resistivity and lithology is required, but this relation is rarely straightforward. Such relations have been subject to many studies and have been researched in the scientific community since the mid 1900's, where e.g. Archie (1942) and Jones and Buford (1951) published on relationships related to resistivity. Since then, scientists worldwide have studied and described the relations between resistivity and different earth materials under a variety of conditions resulting in a high number of mathematical formulas (e.g. Glover, 2015).
In sedimentary environments the resistivity is related to the clay content, clay type, porosity, saturation, and pore water chemistry of the medium. These parameters are closely related to hydrological properties, and lithological models forming hydrogeological models can therefore be translated into 3D groundwater models. When creating hydrogeological models for groundwater models, the modeler therefore needs to have a thorough understanding of both the geophysical data and the geology, and not least how they are related. To our knowledge, the resistivity and lithological relations has not been studied before on a large national scale, except in one similar large-scale study carried out by Beamish (2013aBeamish ( , 2013b, where the conductivity-lithology relations are presented for the entirety of the United Kingdom. The study presented by Beamish specifically focuses on relating 3 kHz apparent conductivity data obtained by an airborne frequency domain electromagnetic system to the existing UK geological lexicons containing lithostratigraphical and lithological maps. This paper presents a general method for relating resistivity models to lithological samples from boreholes applied on a dataset of a national scale. The method is centered on a procedure, which categorizes the resistivity model data according to a representative lithology from the boreholes. From the categorized resistivity data normalized histograms are created for a specific lithology or a group of lithologies. The resistivity models originate from airborne and ground-based Transient Electromagnetic (TEM) soundings stored in the Danish national geophysical database, GERDA (Møller et al., 2009), while the lithological data originate from the Danish national borehole database, JUPITER (Møller et al., 2009;Hansen and Pjetursson, 2011). Combined, the data from these databases provide the statistical basis needed to describe the resistivity-lithology relationships. The large-scale resistivity-lithology relations are portrayed as a compilation of a national Resistivity Atlas of Denmark (Fig. 1). The Resistivity Atlas visualizes the histograms in a grid covering all of Denmark, with 10 × 10 km square cells. One histogram is calculated for each 10 × 10 km cells for each of the dominating lithologies. The final Resistivity Atlas portrays the resistivity-lithology relations in a number of maps, of which an excerpt is presented here.

Geological background
The major part of the Danish subsurface consists of sedimentary basins, which started forming during the Carboniferous and Permian. Primarily, the shallow subsurface consists of soft sediments spanning from Cretaceous chalk, Paleogene limestone, marls and clays, Neogene clay and sand, to Quaternary glacial and interglacial deposits (Liboriussen et al., 1987). The Danish Pre-Quaternary stratigraphic succession is tilted away from the Scandinavian shield. Thereby the youngest sediments from Miocene are found in the southwestern part of Jutland and the deposits become older towards North and East (Fig. 2a).
The thickness of the Quaternary sediments varies between ca. 0 and 300 m and they are deposited in a series of glaciation and deglaciation events (Houmark-Nielsen, 1987;Houmark-Nielsen and Kjaer, 2003). The sediments deposited in these glacial cycles are dominated by till and meltwater deposits and show a high degree of heterogeneity and complexity related to the formation of e.g. buried valleys or glaciotectonic complexes as described by Jørgensen and Sandersen (2006) and Høyer et al. (2013). The youngest Quaternary sediments in the western part of Denmark are dominated by coarse-grained glacial outwash deposits from the Scandinavian ice sheet at the glacial maximum in Late Weichselian (Fig. 2b). The ice sheet expanded from Norway and Sweden in south, southwest and westward directions towards Denmark and reached its maximum position at the Main Stationary Line in the central part of Jutland (Fig. 2b), forming large outwash plains in the low-lying areas of the Saalian landscape Sand-Jensen et al., 2006). To the north and east of the Main Stationary Line, the youngest sediments are dominated by clay till formed by the ice sheet at the latest glacial maximum or younger re-advances in the Late Weichselian.

Data
The data utilized in this study originates from the comprehensive national geophysical and borehole databases of Denmark. The GEophysical Relation DAtabase, GERDA (Møller et al., 2009) comprises a majority of the geophysical data collected onshore in Denmark in relation to groundwater and raw material mapping. In particular, the database includes all ground based and airborne TEM data collected since the early 1990's, which is the main source of data used in this study. As of 2015, GERDA contains about 75,000 ground based TEM soundings and 1.5 million SkyTEM soundings (Sørensen and Auken, 2004), which cover close to 40% of the Danish onshore area (Fig. 3a). The database contains both raw and processed data as well as the resulting resistivity models. There have been several generations of the SkyTEM system, the earliest system is described by Sørensen and Auken (2004). The SkyTEM system has undergone technological developments, which have improved the resolution and depth of investigation of the system over time. The ground based TEM systems used in Denmark include: the Geonics PROTEM 47 system (Geonics, 2013), the High Moment TEM system (HiTEM) and the Pulled Array TEM System (PATEM), as described by Danielsen et al. (2003).
The end result of the data processing and inversion is a 1D resistivity model for each sounding. The oldest data have been inverted either as single site soundings or with a laterally constrained inversion approach , while more recent data are inverted using a spatially constrained inversion approach (Viezzoli et al., 2008). The available inversion models are either, smooth multi-layer models (Constable et al., 1987), or few-layered models (Farquharson and Oldenburg, 1993). The database contains 135 SkyTEM surveys and 427 ground based TEM surveys, which have been acquired, processed, and inverted by different consultancy firms and research institutions primarily in relation to the National Groundwater Mapping Program (Thomsen et al., 2004). A set of standards and guidelines, as well as absolute calibration of the instrument systems at a national test site (Foged et al., 2013) ensure consistency in the data (Møller et al., 2009).
The borehole data used in this study comes from the Danish national borehole database, named JUPITER (Hansen and Pjetursson, 2011). The database currently contains 280,000 shallow wells, drilled in relation to groundwater, raw materials, geotechnical, and research purposes. The database contains information on geology, water levels, drilling process, water supply, abstraction information, groundwater/drinking-water chemistry etc.
The JUPITER and GERDA databases are both managed and hosted by GEUS. They are publicly available via the official homepage of GEUS (gerda.geus.dk; jupiter.geus.dk).

Methodology
The overall goal of the presented method is to produce an easy and statistically sound way of studying relationships between dense geophysical data and borehole lithological logs. Very large datasets are  (Håkansson and Pedersen, 1992) and b) the Quaternary sediments at 1 m below ground surface (Pedersen et al., 2011). The black line in b) displays the glacial maximum of the Scandinavian ice sheet during Late Weichselian. involved, which means that a large part of the procedure is related to importing and organizing data. Once the data have been organized, an algorithm is run. A conceptual diagram portraying the overall steps of the algorithm is seen in Fig. 4. For a given subset of the data the process can be divided into three key steps: 1) Based on the available boreholes and geophysical data the subsurface is divided into N horizontal planes, HP, for regular elevation intervals with a spacing of 0.2 m.
2) For each HP steps a) to c) are run: a) All layer resistivity values intersecting the given HP are retrieved and a semivariogram is created based on the logarithmic resistivity values. b) The semivariogram from a) is applied to interpolate the intersecting resistivity values to the intersecting borehole positions, i.e. one interpolated resistivity value for each borehole intersecting the HP. Only resistivity values within a radius of 400 m from the boreholes are used. This gives an estimate of the resistivity values at the boreholes for the given HP. The utilized interpolation method is Simple Kriging (Journel and Huijbregts, 1978;Pebesma and Wesseling, 1998). c) The interpolated resistivity values from b) are then assigned to the lithology found in the borehole lithology logs. Fig. 5 shows the typical data setup for steps a)-c) for one lithology from one borehole close to a SkyTEM flight line.
3) A summarizing histogram describing the resistivity to lithology relation is created for each lithological category or set of categories, e.g. as seen in Fig. 8.
To compile the national atlas, steps 1-3 are repeated over a 3D grid with 10 × 10 km cells laterally (Fig. 6a) and three layers vertically creating the summarizing histograms from step 3) for each cell in each layer. Vertically, the three layers cover the intervals 12.5-25 m, 25-55 m, and 55-100 m. The upper 12.5 m boundary is chosen to account for a poor near surface resolution of some of the TEM systems. The early TEM instruments have a poor resolution in the upper 15-20 m, while newer systems have poor resolution only in the upper of 5-10 m (Schamper et al., 2014a). Therefore, the 12.5 m value has been chosen as an intermediate value, taking into account that several generations of the ground-based TEM and SkyTEM systems are included in the dataset. The lateral and vertical size of the cells has been prescribed to fit with the overall data density of both the geophysical and borehole data. Following the grid, the spatially distributed histograms can now be visualized by displaying a pixelated map of e.g. the median resistivity or other key summary statistical values for desired lithologies.
The algorithm sketched above is straight forward and intuitive, but there are also a number of built-in issues to be dealt with. These issues are listed below and are handled by introducing weighting factors for the entries in the histograms: a) The resistivity models are interpolated to the borehole positions (item 2b in the list above). This is associated with an uncertainty depending on the geological heterogeneity and the distribution of the resistivity models relative to the boreholes. b) The resistivity data on which the interpolation is based are of varying quality. The models are non-unique, which maps into an uncertainty measure on the individual resistivities, which, in turn maps through the spatial interpolation to the borehole positions. c) The soil samples, and thereby the lithological descriptions of boreholes are of uneven quality, which can be expressed as an uncertainty for the lithological class description. d) Boreholes are commonly clustered in areas of heightened interest, e.g. drinking water abstraction areas, whereas less interesting areas are more sparsely sampled. This causes a very irregular spatial sampling of the lithologies of the model area. e) There are scale differences in the resolution capabilities of boreholes and geophysical data, which necessitate a scaling with depth in the sampling.
To address issue a), regarding the inconsistency in the location of boreholes and geophysical observations, two things are done. Firstly, only boreholes which have one or more TEM soundings within a radius of 400 m are selected (Fig. 6b). The 400 m threshold has been chosen to limit the size of the dataset while assuming that 400 m is larger than the overall correlation length in the typical geological architecture and more than the footprint of a SkyTEM sounding. In other words, boreholes without the required proximity of geophysical data points are discarded as a safety precaution due to the general high degree of geological heterogeneity within glacial deposits, e.g. as seen in Fig. 5. Secondly, the kriging interpolation ensures that uncertainties at the estimation points are estimated taking the interpolation distance into account.
Issue b) is also handled by the kriging algorithm. The resistivity models have a related model parameter uncertainty, which is a linearized error estimate obtained calculating the model covariance matrix (Auken et al., 2015). The model parameter uncertainty can be regarded as relative measures of uncertainty, since the inversion is carried out in logarithmic model space. Therefore, small values can be trusted and larger values are best viewed as guidelines and this is accounted for in the weighting scheme. First, the model parameter uncertainty is used as input to the Kriging routine together with the resistivity value itself. The total Kriging output variances accounting for both the interpolation distance uncertainties and the resistivity parameter uncertainties are converted to standard deviation factors, σ factor , which are then translated into weights, w σ , through a piecewise linear function, where σ max is the threshold standard deviation factor for which all standard deviation factors are assigned the minimum weight, w min , σ min is the minimum threshold value for the standard deviation below which all values are assigned a weight of 1, α ¼ ðw min −1Þ To account for challenge c) the uneven quality of the lithological data, a quality assessment of the boreholes has been computed to assign uncertainty to the borehole lithological samples. This approach has previously been presented by He et al. (2014) and Høyer et al. (2015). The quality assessment groups the lithological input into five quality categories based on a point system. The point system is based on a number of criteria: • Location method: The location method is determined; differential GPS, normal GPS, or manually marked on a map. The more precise the location method the more points assigned to the given borehole. • Sampling frequency: The number of samples taken per length unit (m) during the drilling process. The higher the frequency the more points are assigned. • Drilling contractor: Some drilling contractors are more reliable than others. The reliable contractors are assigned more points than the less reliable. • Drilling method: the sample quality is dependent on the drilling method. Drilling methods producing high quality samples are given more points. • Borehole purpose: The quality of the lithological descriptions often depends on the purpose of the borehole. E.g. seismic shot holes often have poor lithological descriptions, and are therefore assigned fewer points.
Each borehole is assigned a certain amount of points depending on these criteria. The boreholes are then grouped in five quality groups, Q1-Q5, based on their final score. The Q1 group boreholes are of the highest quality, while Q5 boreholes have little to no relevant information. Each quality group is then weighted accordingly. Full weight is generally assigned to Q1 boreholes while reduced or no weights are assigned to the other groups, Q2-Q5, accordingly. The weights, w QA , used in this paper are listed in Table 1.
To account for issue d) the clustering of boreholes, we propose a declustering routine based on Voronoi diagrams, accounting for the sparse and clustered nature of boreholes. The routine down-weights entries from closely spaced boreholes based on their relative Voronoi polygon areas (Fig. 6c). This ensures that the histograms are less dominated by a cluster of closely spaced boreholes, and that the variability of the resistivity field is represented more evenly. The Voronoi de-clustering routine has two steps. First, a Voronoi diagram is made based on the positions of the selected boreholes (Voronoï, 1908) by creating a polygon around each borehole (Fig. 6c). The areas, A voronoi , are calculated for each of these Voronoi Polygons and are then translated into weights, w voronoi , by the piecewise function where, A min is the threshold area below which all Voronoi polygons are assigned the minimum weight value, w min , and A max is the threshold area above which all Voronoi polygons are assigned a weight of 1, α ¼ ð 1−w min Amax−A min Þ and β ¼ ðw min −ð 1−w min Amax−A min ÞA min Þ. The above mentioned weights are combined into one total weight matrix, containing the cumulated effect of all the individual weights. Combining the weights in such a scheme makes it easy to weight the histograms, and also to implement new weights. The final weighting matrix, w tot , created by combining the individual weights is where K is the number of sampled data points to be weighted in the given histogram, and w k is a hypothetical kth weight, implying that any number of desired weights might be applied in this flexible weighting scheme. The total weight matrix, w tot , is then applied to the histogram, which is then re-normalized to 1 by dividing the sum of weighted counts by the applied weights. This ensures that the weighted histograms are comparable to each other, and that the normalized counts are representative fractions of the total area of the histograms.
Finally, addressing issue e), concerning resolution discrepancies, we are dealing with a common problem when relating geophysical information to geological information. Generally, boreholes resolve thin geologic layers better than most geophysical methods. For airborne and ground-based TEM systems the resolution capability towards an individual layer is dependent on many factors, such as the conductivity and depth of the layer. The more conductive the layer, the better the resolution Schamper et al., 2014b), and the deeper the layer, the poorer the resolution will be. A solution to this is to include only borehole layers, which are considered resolvable by the geophysical method. The TEM and SkyTEM resolution can roughly be described as exponentially decreasing (Schamper et al., 2014b). In this study, threshold values are chosen as 12.5 m at a depth of 1 m, exponentially increasing to 100 m at a depth of 150 m, with an acceptance threshold of 80%. The proposed borehole filter is implemented by considering the percentage of the target lithology for exponentially increasing intervals. The filter is implemented to run for an entire borehole, including the first 12.5 m which, as explained at the beginning of this section, are disregarded due to the resolution capabilities of the Fig. 6. a) Map of the grid used in the Resistivity Atlas of Denmark with boreholes and TEM soundings; the grid cells are 10 × 10 km in size. b) the data selection within a single grid node (node 294 marked with green in a)), where boreholes with a TEM sounding within 400 m are marked in green, and discarded boreholes are marked in red. c) Voronoi polygons for a selection of boreholes containing a specific lithology, based on the boreholes which are marked in green in b).

Table 1
Quality assessment weights, w QA , are applied to all histograms throughout this paper, unless otherwise stated.
Quality group Q1 Q2 Q3 Q4 Q5 w QA 1 0.8 0.7 0.6 0 SkyTEM system. The upper 12.5 m are removed by the Resistivity Atlas grid layer constraints. The filter scans a given borehole from top to bottom with a particular lithology in focus. The filter starts by selecting the first HP intersecting the given borehole as the focus point, along with the borehole lithology. Afterwards all borehole lithologies found within the first filter interval, stretching from 0 m to 12.5 m including the focus point, are selected. If the selection consists of more than 80% of the target lithology, then the focus point is accepted; if not the focus point lithology is discarded. The filter then moves to the next focus point, being the next successive HP intersecting the borehole, and enforces the 80% threshold value on a slightly increased filter interval thickness. This procedure is continued for all intersecting HP's and as the filter moves deeper, the interval thickness increases significantly, e.g. at a depth of 75 m the interval thickness is 32.7 m, and at a depth of 125 m the thickness is 70.5 m. The end result, once the entire borehole has been scanned, is a borehole lithology log with approximated SkyTEM resolution, e.g. case 3 in Fig. 7.

Results
The resistivity-lithology relation within one 10-by-10 km cell (outlined in green on Fig. 6a) is presented as histograms for two lithological categories (Fig. 8), which are the two most frequently occurring lithologies from the borehole lithology logs from the JUPITER database and are not to be viewed as a binary sand/clay classification. The first category consists of meltwater sand, gravel, and pebbles; the second category contains clay tills. The definition of the lithological categories is based on the fact that the processes and elements defining the resistivity-lithology relations are quite complex. One important element is the physics behind the TEM method, which in the end is what defines the resistivity models. The TEM method is especially sensitive to conductive materials, such as clays, while not as sensitive towards resistive structures such as sands and gravels. In fact, TEM instruments have a difficult time resolving details in resistive structures above 150-200 Ω m Jørgensen et al., 2005). This means that it is difficult to differentiate two coarse-grained layers, such as dry sand and gravel, therefore meltwater sand, gravel, and pebbles have been lumped together into one category. Clay tills, on the other hand, generally have lower resistivity values and are abundant in the Danish geological setting and can therefore form a lithological category of their own.
Histogram summary statistics are stated in a table for each of the histograms in Fig. 8. The statistics include the number of interpolated resistivity data points, the histogram median, interquartile range, number of unique boreholes, and the core length. The interquartile range, IQR, is used to describe the width of the histogram. It is given by the difference between the 25th and 75th percentile, calculated for the log 10transform of the resistivity data. Fig. 8a shows the histograms for the two lithological groups without any weighting or filtering, whereas weighting and filtering have been applied in Fig. 8b. It is clear that the coarse meltwater sediments and the clay tills have overlapping resistivity values. By applying the filtering and weighting schemes the overlap between the two histograms is reduced and the overall shape is altered. For this specific cell, the median of the clay till is reduced from 46 Ω m before filtering to 30 Ω m after filtering, and the IQR of the log data is reduced from 0.31 to 0.23. For the coarser meltwater category, the median resistivity value is increased from 77 Ω m to 90 Ω m, and the IQR is slightly decreased from 0.32 to 0.31. Many of the high resistivity values (50-200 Ω m) in the clay till histogram are removed by the weighting and filtering procedure. The weighting and filtering does not have the same clear effect on the meltwater sand, gravel and pebbles category, where the median is increased, but the width of the histogram is only slightly decreased. This can be partially explained by the relatively poor sensitivity towards resistive materials when using EM methods. Removing layers by filtering the boreholes still means that the remaining layers are less determined compared to more conductive materials i.e. spans a wider range due to low sensitivity. Another explanation could be that the sands/gravels actually do span a wider range in the log(resistivity) scale, which would make sense as the EM method is actually sensitive to lin(conductivity).
Overall, the combined filtering and weighting procedures reduce the overlap of the resistivity-lithology histograms for the two lithological categories. Although not shown, we conclude from several experiments that the general effect of the filtering is to separate the two lithology categories and reduce the width. The weighting scheme, on the other hand, does not reduce the overlap of the histograms significantly, but by reflecting over the different types of uncertainty and spatial problems it cleans up the resistivity-lithology histograms. The effects of the different weighting procedures differ across the dataset. The borehole declustering has the largest effect in areas where boreholes are more clustered. If a model area contains evenly distributed boreholes the declustering weights will become equal and therefore have no effect. The borehole quality assessment has the largest effect in areas with old boreholes. Here, the old boreholes, which are often of poor quality, are down-weighted. Finally, the combined Kriging interpolation and TEM standard deviation weights depend on the model parameter uncertainty and the average distance between TEM soundings and boreholes, respectively. The average distance between soundings and boreholes is large in areas where only ground-based TEM exist and in urbanized areas, where soundings are processed and removed due to a general high degree of noise . In these areas the Kriging weights will have a high influence.
An excerpt from the full Resistivity Atlas is shown in Figs. 9 and 10. The figures portray the median value of the histograms of the first and second Resistivity Atlas layers, ranging from 12.5 to 25 m (Figs. 9a and Fig. 7. Borehole 108.154 displayed in Fig. 5 shown in three different cases: 1) the selection of the clay till layers in green, 2) the intersecting horizontal planes with a spacing of 1 m, again the selected clay tills are marked in green and 3) same as 2) but after filtering related to the clay till layers. 10a) and 25 to 55 m (Figs. 9b and 10b). The third and final layer of the Resistivity Atlas, spanning 55-100 m, contains very few data, both due to the general lack of boreholes with depth and the filtering procedure, e.g. as seen in Fig. 7. The excerpt shows the two aforementioned lithological categories, the first being the meltwater sand, gravel, and pebble category (Fig. 9) and the second being the clay till category (Fig. 10). For each colored 10-by-10 km pixel a histogram is computed, provided data is present in the given cell. In some cases very few data are present, and the lower limit for showing the given pixel has been chosen as at least four unique boreholes and more than 50 interpolated resistivity values must be present. Generally the IQR is low, with 90% of the nodes having an IQR below 0.45, distributed in a pseudo exponentially decreasing fashion, meaning the majority of the values lie around 0.2-0.3.
Figs. 9 and 10 reveal that the overall resistivity patterns vary geographically, for both of the lithological categories. The overall variations are closely related to large-scale geological processes and the general geology of the Danish area. However, depending on the sediment type different controlling factors dominate the bulk resistivity value. In  short, for coarse and porous sediments, such as the meltwater sand, gravel, and pebbles category, the pore water ion content is essential, since most of the electrical current is conducted in the pore water, e.g. Glover (2015). However, if the sediment matrix contains clay minerals the current flow additionally takes place along the clay mineral surface, and the bulk resistivity value also depends on the clay content and the clay mineral composition and therefore on the matrix materials (Waxman and Smits, 1968). Hence, for clay-rich sediments, such as the clay till category, water chemistry has a relatively small influence, while source material and depositional processes dominate the bulk resistivity. The bulk resistivity is controlled mainly by the source materials that were available to the depositional processes depositing the sediments. The geological processes, which in this case are glacial, rework the material found at the substrate-ice interface. In areas of coarse source materials, i.e. abundant amounts of sand, the clay tills inherently contain more sand and gravel, as opposed to areas with more clay rich source material, where more clay would be found in the clay till matrix.
By closer inspection of the Resistivity Atlas (Figs. 9 and 10), it is seen that there are generally more cells containing meltwater sand, gravel and pebbles than clay till; in particular to the west of the main stationary line at the last glacial maximum presented in Fig. 2b. The sand deposits found in the western and northern part of Jutland has a higher resistivity value on average, with a majority of the median values lying between~100 and 500 Ω m. On the islands, Funen and Zealand, the median resistivity of the meltwater sands, gravels, and pebbles are generally lower, between~35 and 125 Ω m (Fig. 9). It is also seen generally that the number of data decreases with depth, i.e. fewer colored pixels are found in the second layer (Figs. 9b and 10b). This decrease is partly due to the shallow surface nature of boreholes, where many shallow boreholes with depths b 25 m exist, but fewer boreholes with depths N 25 m exist and partly due to the filtering skipping layers assumed to be below the TEM resolution. The overall patterns seem to persist with depth, though the global average of the median resistivity values seems to decrease in the second layer (Figs. 9b and 10b) for both lithological categories.
The high resistivity trend seen in the western part of Jutland for the meltwater sands, gravels and pebbles is partly explained by the ice margin seen in Fig. 9. To the west and south of this line, outwash deposits along with meltwater sands and gravels dominate as cover layers.
Looking at Fig. 2b, this trend also seems to persist northeast of the main-stationary line where the cover layers generally consist of meltwater sands and gravels.
A different trend is seen for the clay till category in Fig. 10. Not a lot of clay till is found west and south of the ice margin at last glacial maximum. The clay tills found in this area, striking southwest-northeast, generally have relatively high median resistivity values ranging be-tween~35 and 250 Ω m. The clay tills located in the eastern part of Jutland, east of the glacial maximum, show lower mean resistivity values between~20 and 50 Ω m. This corresponds to the presence of hemipelagic Paleogene clay just below the Quaternary deposits (Fig. 2a) and along the southeastern and eastern coast of Jutland. The lower median resistivity values can therefore be explained by the proximity of Paleogene clays which has been worked into the clay till matrix during deposition. However, on Zealand higher resistivity median values are seen with mean resistivity values generally ranging betweeñ 40 and 125 Ω m. These trends can, again, be related to the overall geologic setting. If chalk with a particle size similar to clay is worked into the till matrix the bulk resistivity value is larger due to a smaller amount of clay minerals present for electrical conduction along the clay mineral surface. We recall from the geological map, Fig. 2a, that on Zealand the pre-Quaternary sediments are dominated by chalk and limestone, corresponding to the generally higher median resistivity values. Therefore, it is likely that limestone and/or chalk is worked into the clay till matrix, increasing the bulk resistivity values.

Discussion
The goal of the Resistivity Atlas is to study the resistivity-lithology relations on a large scale, contrary to the many local and regional studies carried out over decades leading to an overall understanding of the resistivity-lithology relations among geoscientists. An understanding based on various types of geological interpretation and modeling based on geoelectrical and electromagnetic data. It is well-known that these relationships vary depending on a number of factors, such as pore water resistivity, matrix material, clay content, clay type, grain size, post depositional processes, etc. Therefore, the resistivitylithology relations must be regarded to some extent as site-specific. In this sense the resistivity-lithology relations are ad-hoc, since a new relation must be established every time a new area is investigated. The Resistivity Atlas addresses these large-scale spatial variations in the resistivity-lithology relations and attempts to describe the overall patterns on a national scale. As can be seen from the presented excerpt of the Resistivity Atlas of Denmark (Figs. 9 and 10), there are some significant spatial variations in the relationship between resistivity and lithology.
In this paper it was attempted to remove as much bias as possible from the resistivity-lithology histograms. There are, though, inherent problems that cannot be fixed in the weighting scheme. One is the manner in which geologic descriptions are carried out, which influences the final resistivity-lithology histograms. Generally, lithologies are classified according to the depositional process, which has deposited them, e.g. tills are simply defined as unsorted glacial sediment, hence any sediment, which has been deposited by a glacier. Clay tills are therefore clays, which are unsorted, but which contain enough clay to be classified as clay. Therefore, clay tills can be very sandy and can contain as much as 88% sand, or quite clayey with more than 20% clay with some of the most clay rich clay tills containing up to 40% clay. This is defined in the Danish scheme for describing borehole geology (Ditlefsen et al., 2008). Obviously the bulk resistivity signature will be very different for the clays in either end of this spectrum.
The lithologies chosen for this study were selected based on the fact that they are the most frequent in the Danish borehole lithology logs. Alternative Resistivity Atlas maps could also be created for any other lithologies found in the borehole lithology logs; however most of these maps are very sparse. This is partly due to the filtering scheme proposed in Section 4. Methodology, and partly due to the low occurrence rate in the lithology logs.
The high degree of heterogeneity present in Figs. 9 and 10 does not seem geologic in nature and in reality they are probably smoother transitions. A number of solutions could make these sharp boundaries less sharp, e.g., refining the grid or using a moving window to compute the histograms. The problem is that the data density is not high enough for either of these approaches and keeping the 10 by 10 km grid cells makes it easier to communicate exactly what data values are included in the histograms as the cells do not overlap.
Information regarding the coarseness of the given sediment is valuable information in hydrological applications, and the Resistivity Atlas might provide some insight into the coarseness of clay tills in a given area. If a given cell portrays a relatively high mean resistivity value, this might indicate that the clay till is coarser in this area, i.e. has a lower clay content, and vice versa. Although the overall trends in the geologic environment are also required, as for example the high resistivity clay tills on Zealand most likely being caused by chalk or limestone in the sediment matrix.
The statistical resistivity-lithology relation is not simple and significant overlaps are present, as seen in Figs. 8,9,and 10. Consequently, it is not trivial to interpret and translate resistivity models into lithological models. Lately, automatic and semi-automatic approaches have been suggested to build geological and/or hydrogeological models (e.g., Bosch et al., 2009;He et al., 2014;Foged et al., 2014;Gunnink and Siemon, 2015). The Resistivity Atlas may prove an efficient tool in building these models as it provides spatial input on the link between lithological and geophysical information. The Resistivity Atlas also reveals the large-scale resistivity-lithology patterns across the entirety of Denmark, making it easier to relate these patterns to the large-scale geologic structures, e.g., pore water ion content and source material.

Conclusion
A method for studying the resistivity-lithology relations on a large scale has been presented, along with an excerpt of a national Resistivity Atlas of Denmark for the two layers, in the interval 12.5-25.0 m and 25-55 m and two lithological classes. The presented method integrates lithological information from boreholes with geophysical resistivity models to create spatially distributed resistivity histograms for selected lithological groups. By applying weight functions the presented approach accounts for the clustering of boreholes, borehole uncertainty, resistivity model uncertainty and interpolation uncertainties. Resolution differences between boreholes and the geophysical methods are handled through a thickness filtering algorithm.
The presented Resistivity Atlas reveals varying patterns in the largescale resistivity-lithology relations, reflecting numerous geological details such as available source material for tills. The resistivity maps also reveal a certain amount of ambiguity in the resistivity values for different lithologies. In some areas this ambiguity has been revealed to be quite large, whereas it is much less in other areas. Such information is crucial when geophysical data are to be used for lithological/geological modeling or other derived products such as hydrological models.