Holistic 3D Model of an Urban Area in Norway: An Integration of Geophysical, Geotechnical, Remote Sensing,and Geological Methods

Abstract: Understanding the distribution and characterization of natural and non-natural materials on the surface and near-subsurface is important for the development of infrastructure projects. This may be a challenge in highly urbanized areas, where outcrops are scarce, and anthropogenic activities have altered the morphological expression of the landscape. This study tests the integration of ground-penetrating radar (GPR), borehole drilling, aerial imagery, geological mapping, and aerial laser scanning as complementary mapping tools for determining the stratigraphy of glacial and post-glacial Quaternary sediments, the depth to the bedrock, and the distribution of anthropogenic material in Mosvatnet, a lake in Stavanger, Norway. The integration proved to be efficient and enabled the generation of a 3D holistic model, which provided a broad understanding of the subsurface geology and the induced anthropogenic changes in the area through time. Bedrock, till, fluvioglacial, and lacustrine geological units were modeled. Accumulations of post-glacial organic matter were mapped, and the distribution of non-natural infill material was determined. The interpreted dataset suggests that the eastern shoreline of Mosvatnet has artificially prograded about one hundred meters westward since the 1930s and the elevation of the corresponding area has increased by about ten meters relative to the lake level.


Introduction
More than 50% of Norwegian land areas consist of bare bedrock or bare bedrock with a thin cover of Quaternary sediments [1]. In the remaining regions, Quaternary sediments are thicker, more continuous, and mostly consist of unconsolidated glacial tills deposited during the transition between glaciation and deglaciation periods about 15 ka years BP. Understanding the distribution of Quaternary sediments and their positions over the underlying bedrock in mainland Norway has been fundamental for developing anthropogenic activities, such as agriculture, construction of cities, and mining and quarrying. These activities have greatly transformed the natural landscape through time, particularly within and around urban agglomerations, where it is sometimes difficult to discern between natural and anthropogenic materials that constitute the structures on the surface and the near subsurface.
The characterization of geological structures within an urban landscape may be challenging due to the scarcity of natural rock exposures and the absence of the original geomorphological landforms. In Norway, geotechnical drilling is traditionally used in urban infrastructure projects to map and characterize unconsolidated Quaternary sediments in the subsurface. Although geotechnical boreholes provide hard data (i.e., tangible evidence that provides the ground truth) from the subsurface, the information is derived from sparse observation points. Therefore, given the variability and heterogeneity of Quaternary sediments, using only geotechnical drilling may introduce a high degree of uncertainty in the subsurface analysis. Furthermore, geotechnical drilling is an invasive method limited by logistical and environmental restrictions. Geophysical methods, such as ground-penetrating radar (GPR), increase confidence by filling the information gaps between geotechnical boreholes. GPR is a non-invasive method that enables the collection of geological information by performing near-surface geophysical investigations. The method has been widely used in different subjects in Norway (archaeology, e.g., [2,3]; glacial and Quaternary studies, e.g., [4,5]; geotechnics and geohazards, e.g., [6][7][8]; and construction and urban infrastructure mapping, e.g., [9]). Compared to geotechnical drilling, GPR is a cheaper and more urban-friendly alternative that does not imply big logistical efforts. However, GPR is an indirect method, which delivers soft data (i.e., non-ground truth information that requires interpretation), and its imaging contrast and depth of ground penetration capabilities are limited by the varying physical properties of sediments [10] and local environmental conditions [11]. Although geophysical methods may help to reduce costs and environmental impacts by reducing the number of drilled boreholes, a minimum amount of boreholes is required as ground truth to correlate geophysical responses with natural and non-natural materials in the subsurface.
Combining geotechnical drilling and GPR methods with other techniques traditionally applied in geosciences as complementary tools have proven to be efficient when characterizing Quaternary deposits on the near surface. Pellicer and Gibson, 2011 [12], combined GPR, geotechnical boreholes, electrical resistivity tomography (ERT), and geological fieldwork to map poorly exposed glacial and post-glacial Quaternary unconsolidated sediments in the Irish Midlands. Similarly, Sauvin et al., 2014 [13], integrated in situ geotechnical measurements and GPR along with other geophysical methods, such as ERT, seismic refraction tomography, and seismic reflection profiling, to map the distribution of quick clay in Hvittingfoss, Norway. A lithofacies characterization and the description of the depositional geometry of a Holocene prograding delta were conducted by Eilertsen et al., 2011 [14], in Norway, using GPR, auger drilling, sedimentary logs, and aerial imagery. Regli et al., 2002 [15], used GPR, drill cores, and outcrop methods to generate a probability estimation approach for interpreting the stratigraphy of the Rhine/Wiese aquifer, Switzerland. Burki et al., 2008 [16], and Lerche et al., 2014 [17], combined geophysical surveying with remote sensors and field measurements and observations to study the architecture of moraines. To explain the formations of sawtooth moraines at Bødalen, Norway, the former employed GPR, geomorphological mapping on aerial imagery, and sedimentological studies. The latter interpreted parallel ridge topographic features as ribbed moraines in Himmerland, Denmark, based on morphological observations made on an aerial laser scanning-based digital elevation model (ALS-based DTM), geological field mapping, and the interpretation of GPR transects.
Workflows that integrate diverse techniques, such as geotechnical drilling, remote sensing, geophysical surveying, and geological field mapping, have not been widely explored for studying Quaternary glacial sediments in urban areas, where there is poor availability of outcrops, and the morphological expression of the landscape has been largely altered by anthropogenic activities. This study aims to test the integration of the abovementioned techniques as complementary mapping tools for determining the stratigraphy of glacial and post-glacial sediments along with the distribution of anthropogenic material associated with the development of Stavanger, Norway.
This study combines datasets acquired by applying diverse techniques, at different epochs, on different media (land, water, and ice), and in different seasons. The study used datasets collected by different parties: (1) onshore and offshore GPR and geological field observations collected by the authors; (2) publicly available datasets collected by the Norwegian government: historical aerial imagery and ALS-based DTM; and (3) historical geotechnical drilling information provided by the city government offices of Stavanger. Except for geotechnical drilling, the techniques employed for data collection are non-invasive and non-destructive, making them highly suitable for urban area surveys. The study area focuses on Mosvatnet, the third-largest lake in Stavanger, located 1.5 km SW of downtown Stavanger. Mosvatnet is considered a national bird sanctuary. An islet that acts as a bird nesting ground lies on the eastern side of the lake. With the expansion and urbanization of Stavanger, the construction of transport infrastructure has increased. Such urban development has altered the lake footprint, the surrounding landscape, and the geological record, mostly on the eastern side. Nonetheless, most of the lake has not undergone noticeable human-induced changes, enabling the geological mapping of sediments lying in the lakebed using remote sensing and geophysical techniques.
The main objective of this study is to develop a workflow that enables the combination of diverse techniques used in geosciences to generate a holistic model that facilitates the unraveling of the stratigraphy of the subsurface, and at the same time, assists in the detection of morphological changes induced by humans in an urbanized landscape. In order to accomplish the objective, this study focused on answering the following research questions:

1.
What methodological approach is required to combine remote sensing, geotechnical data, geological survey, and geophysical tools to create a holistic 3D model in a highly urbanized area? 2.
How can a holistic model be created that incorporates the different remote sensing, geological, geotechnical, and geophysical data from onshore to offshore and from subsurface to surface? 3.
What are the determining factors to identify anthropogenic landscape changes in the Mosvatnet area? 4.
What are the identifiable geologic units in the area?
The datasets of this study are presented together with the methods employed for data collection. Steps for data processing, interpretation, and integration are described. Subsequently, the 3D holistic model of Mosvatnet is presented as a result derived from the data integration. Finally, the model is discussed along with its geological uncertainties associated with the limitations of the employed techniques and how these techniques complement each other to fill gaps in information.

Geological Setting and Location
The present-day Norwegian landscape results from cumulative geomorphological processes associated with the repeated advancing and retreatment of glaciers during the last three million years [1]. However, most of the Quaternary glacial landforms and deposits that remain onshore are derived from Weichselian glaciation [18]. After deglaciation (about 15 ka BP), more than 90% of the post-glacial sediments have been reworked, transported, and deposited offshore by Holocene fluvial, gravitational, and coastal erosional processes [19][20][21]. The remaining Quaternary onshore sediments are mainly dominated by unconsolidated glacial tills, of which, less than 10% are older than 115 ka [19].
Quaternary geological maps of mainland Norway produced by the Norwegian Geological Service (Norges Geologiske Undersøkelse (NGU); http://www.ngu.no/kart/ losmasse/, accessed on 27 December 2020) [22] separate areas covered by continuous till deposits thicker than 0.5 m from areas with exposed bedrock. The latter includes areas where the bedrock is overlain by discontinuous deposits of till sediments thinner than 0.5 m (inset map of Figure 1a).   [19], calculated a 2-m average thickness for glacially derived unconsolidated deposits everywhere in mainland Norway; 25% of the area covered by these sediments has a thickness of 6 m, on average, whereas the sediment thickness in the remaining 75% of the area is 1 m.
The Jaeren region in southwestern Norway is one of the areas characterized by continuous and thick sediment cover (inset map of Figure 1a). This region marks the transition zone between the North Sea shelf and mainland Norway [23]. Sediments record the sea level rise associated with ice melting during glacial retreatment. In the Stavanger area, north of Jaeren, the post-glacial marine limit is approximately 30 m above the present-day sea level (asl) [24][25][26]. Above this limit, glacial sediments have remained unaffected by marine processes. Areas shaded in warm colors in Figure 1a represent regions above the post-glacial marine limit in the Stavanger area.  Stavanger has been settled on Quaternary till sediments that overlie the Cambro-Ordovician Ryfylke phyllite [28]. The phyllite constitutes the lower unit of the Scandinavia Caledonides Allochthonous complex that overthrusts the Baltic Crystalline Shield (BCS) [28][29][30]. In Stavanger exposures, the phyllite is gray to greenish-gray, sometimes black with quartz veins and intense microfolding. The foliation strikes approximately northeast-southwest (NE-SW) and predominantly dips towards the northwest [28]. Quaternary till sediments in the Stavanger area have been altered by anthropogenic activities and mixed with non-natural infill materials associated with the city's development. However, to the north of Hålandsvatnet and Stokkavannet, glacial landforms (as indicated by black arrows shown in Figure 1a) remain mostly "unaltered" by anthropogenic processes. These smooth, oval-shaped hills, known as drumlins [31,32] are oriented SWW-NEE and indicate a main glacial flow direction from the northeast [26,33].
The study area is located on the eastern side of Mosvatnet, a lake situated within the urban area of Stavanger (Figure 1b). In the south and east, the lake is bounded by Ullandhaug and Vålandshaugen highs, respectively (Figure 1a). It has an average depth of 2 m (Figure 1b), and its water surface is 35.7 m.a.s.l, i.e., ca. 6 m above the late post-glacial marine limit. According to NGU maps, the current lake footprint overlies the contact between thick (>0.5 m) till deposits and peats located at the northern shoreline (Figure 1b). A forest and a walking trail of ca. 3 km long surround the water body with an area of ca. 0.45 km².

Materials and Methods
The datasets employed in this study were acquired in five different surveys at different times and by different parties. Aerial imagery, ALS-based DTMs, local and regional geology, and groundwater and geotechnical boreholes have been acquired and produced by the Norwegian government and are in the public domain. Geophysical and additional local geologic data were acquired by the authors in October 2020 and February 2021. The locations and spatial coverage of the datasets are shown in Figure 2. Since the datasets were derived from different sources operating in different coordinate reference systems, the whole project was georeferenced and projected in WGS84 UTM zone 32 North. The local Norwegian vertical datum NN2000 was set as the Z coordinate system.   Groundwater boreholes -NGU ! .

Aerial Imagery
Recent and historical digital orthophotos taken from aircraft were downloaded from the Norwegian aerial imagery repository (https://norgeibilder.no/, accessed on 24 February 2021) [34]. The orthophotos are organized by year in the historical record. The newest photos were taken in the spring of 2020. They showcase the present-day lake footprint embedded in an urban landscape (Figure 2a). The oldest photos were taken in the fall of 1937. They depict a more "natural" lake footprint surrounded by a rural landscape (Figure 2b). Until the 1970s, the record contained single-band imagery with pixel-size resolutions ranging from 0.5 m (e.g., Figure 2d) to 0.20 m. On the other hand, the most recent imagery (from the 1980s) consists of RGB orthophotos with resolutions ranging from 0.125 to 0.04 m (e.g., Figure 2c) pixel sizes.
RGB orthophotos taken in the spring of 2007 show similar landscape characteristics to those observed in the 2020 imagery, particularly the geometry of the lake shoreline. Moreover, photos collected in 2007 with 0.08 m pixel resolution, allow better visibility of the underwater sediments lying in the lakebed. To enhance the visibility of lakebed sediments, image analysis tools were employed. These tools enable the exploitation of meaningful information and display capabilities of orthophotos. The following are the four main steps used to obtain the image that best highlights the distribution of the sediments in the lakebed:

1.
Mosaicking: Consists of merging orthophotos with the same characteristics into a raster dataset.

2.
Image classification: Water body and terrain pixels were separated based on the spectral characteristics of each area. These characteristics are defined by the attributes inherent to each color band. The image-supervised classification workflow was applied in the ArcGIS Desktop (version 10.7.1) [35]. The output was a classified raster with two classes: water body and ground. The water body was extracted from the classified raster by using the polygon that outlines the water body.

3.
Image Enhancing: Application of image analysis tools in ArcGIS [36] to enhance the appearance of the extracted raster data representing the water body. Three parameters were equally adjusted in the red, green, and blue bands: the gamma correction, the contrast stretch, and the display resampling.

4.
Image denoising: Noise filtering was required to remove random brightness variations. The image from step 3 was imported into Matlab (version R2021a) [37] and separated into RGB bands. Then, a pre-trained denoising neural network, DnCNN, was applied individually to each band. After applying the filter, the bands were recombined to deliver the final image, which was used for further interpretation (Section 4.5).

Airborne Laser Scanning (ALS) Derived DTM
ALS uses light detection and ranging (LiDAR) technology to produce highly accurate 3D point clouds that contain both ground and non-ground (e.g., vegetation, buildings) spatial information. Ground points are the only ones required to generate DTMs.
This study used two DTMs derived from point clouds collected by the Norwegian government in the 2014 and 2019 ALS surveys. The data are public and can be accessed online (https://hoydedata.no/LaserInnsyn/, accessed on 10 December 2020) [38]. The 2014 ALS survey was acquired with a Leica ALS70-HA scanner in April 2014. The point density of the clouds was approximately 5 points/m². The accuracy was 0.04 m, which corresponds to the standard deviation obtained after comparing measured geodetic ground control points and the point clouds.

Groundwater and Geotechnical Borehole Data
Groundwater and geotechnical borehole data provide insight into either natural or non-natural materials in the subsurface. Drilling reports and GIS shapefiles of nine groundwater boreholes were downloaded from the Norwegian National Groundwater Database (GRANADA; https://geo.ngu.no/kart/granada_mobil/, accessed on 27 December 2020) [39]. The groundwater boreholes are located south of the study area ( Figure 2) and only provide the depth to the bedrock.
The information derived from geotechnical boreholes was extracted from a geotechnical data report provided [40] by the Stavanger municipality. The boreholes were drilled between 1999 and 2013 as part of the Eiganes tunnel project [40]. The drilling information details the vertical distribution and composition of the subsurface layers and the depth to the bedrock. A table containing borehole coordinates and bedrock depths is included in the report, along with geotechnical plots and laboratory essays. Table 1 summarizes the drilling information available for this study. A total of forty-eight geotechnical boreholes located on the eastern side of Mosvatnet ( Figure 2) were used in this study. Most boreholes were drilled over an area occupied by the water body until the 1970s (Figure 2a,b). Single-and total-sounding techniques were applied to drill the geotechnical boreholes. Single-sounding and total-sounding techniques are Norway's most popular drilling techniques for geotechnical investigations. The standards and parameters defined under the Norwegian Geotechnical Society (Norsk Geoteknisk Forening (NGF); https://ngf.no/, accessed on 2 March 2021). Guidelines [41][42][43][44] are mainly tailored for mapping poorly consolidated sediments and their contacts with the underlying bedrock. The single-sounding technique consists of only percussion hammering for determining solid layers, rock boulders, or the assumed bedrock [40,42,43]. The total-sounding technique combines rotary pressure sounding and bedrock control-sounding [40][41][42][43]45]. This method is intended to map the soil layers and determine the depth of the bedrock. The employed drilling system delivers geotechnical plots (Figure 4), which were interpreted under the NGF guidelines [41][42][43][44]. The penetration resistance force (FDT, kN) is the most important parameter since it helps to estimate the firmness of the soil layers [45]. High-resistance materials, such as coarse silt, sand, and gravel, exert significant friction when in contact with the rod. They are characterized by a serrated curve character alternating between low-and high-resistance values. Conversely, clay and fine silt sediments exert low friction on the rod, and in general, they are characterized by smooth curves with resistance values lower than 5 kN at a depth shallower than 10 m. A hard rock is represented by resistance values ranging between 5 (sometimes lower) and 10 kN, and a serrated curve character ( Figure 4).
In addition to the geotechnical drilling information, eight laboratory essays consisting of particle size distributions and organic matter analyses were provided for this study. The analyses were carried out on the samples collected in the geotechnical boreholes under NGF standards [44]. The results of these analyses were used to define lithofacies under the Folk (1954) scheme [46].

Geologic Observation Points (GOP)
The geological survey aims to identify lithological units that crop out within the study area in order to establish a correlation with those units interpreted in the boreholes, the geophysical survey, and remote sensors. Lithological description and geopositioning using GNSS with RTK correction were carried out on each GOP. Structural data of geological features were collected when possible. GOPs in the main water body were accessed by walking on top of the frozen lake. In addition, a hand auger was used to collect two samples of organic matter from the lake. In total, nineteen GOPs were identified and recorded ( Figure 2).
More than eight kilometers of GPR data were acquired in the study area at different seasons (fall and winter), on different mediums (water, ice, and ground), and with different levels of positioning accuracy (sub-meter and sub-centimeter). The temporality and length of the survey, as well as the setup of the equipment used for data acquisition, are summarized in Table 2. The GPR survey was performed with a single channel common offset Malå GX HDR system equipped with an 80 MHz shielded antenna mounted on a skid plate. An encoder wheel attached to the system enabled distance measuring when the data were collected over the ground and the ice. The system was switched to time-trigger mode when the data collection took place on the water. The built-in DGPS provided geographical positioning with sub-metric accuracy. The GPR system uses the European Geostationary Navigation Overlay Service (EGNOS) to enhance the DGPS accuracy to less than one meter. sub-centimeter accuracy was achieved by attaching a high-accuracy external GNSS rover to the system. The first acquisition occurred in October 2020, when the environmental permit for data collection was granted. After detaching the skid plate, the GPR was mounted on the rear part of a rowing boat. The survey was intended to map the area around the islet and the area between the islet and the closest points in the shoreline ( Figure 2). Although the boat was rowed across the lake following a random pattern, the surveying lines were tracked in real-time while navigating. The tracking was performed using a tablet device synchronized with a GNSS onboard the boat. The built-in DGPS was used to obtain the geopositioning of the survey. Rock boulders and vegetation prevented access to shallow areas near the islet and shoreline.
Due to the low temperatures in February 2021, the lake froze to a safe thickness (ca. 0.2 m) for walking, enabling the second GPR acquisition. This acquisition includes an on-land survey. Both the on-ice and on-land surveys took place on the same day and were performed in the traditional manner: the skid plate was attached to the GPR system, which was subsequently dragged along the ice and ground. A GNSS rover was externally connected to the GPR via the COM port. The rover delivered centimeter and sub-centimeter accuracy when synchronized via GSM mobile with the Norwegian RTK network, CPOS. The second survey was intended to ( Figure 2): (i) deliver a better definition of the transition between the lakebed and the islet, and the lakebed and the shoreline; (ii) to help correlate GPR responses with ground truth rock exposure on GOPs; and (iii) to help correlate GPR responses with geological information provided by the geotechnical boreholes.
Due to the amount of collected information and the diversity of the data acquired in different mediums, the processing workflow had to be automated. Some steps were tailored to the data acquired in each medium accordingly. The processing workflow was implemented in Matlab (version R2021a) using in-house scripts that process multiple GPR lines simultaneously. GPR lines were systematically processed as follows: 1.
Extraction of data from .rd3 and .cor files; 2.
Binning to 12.5 cm distance intervals using the median of traces and coordinates; 3.
Application of the amplitude gain function to preserve anomalies; 7.
Kirchhoff migration taking care of time-zero correction and topographic values relative to a specific datum. Two dielectric permittivity constants, 80 (33.52 m/µs) for the lake (water and ice) and 25 (60 m/µs) for the land, were used to migrate the data with a summation width of 40 traces; and finally 8.
Writing SEG-Y files to be imported into the geological interpretation software.
Binning and interpolation are very important steps for migration. They ensure that the data are well distributed or spaced in relation to the data coverage, regardless of sampling discrepancies during the data acquisition. The velocity model used for migrating on-water data was also assumed to migrate the data acquired on ice/water. This assumption is based on the thickness of the ice layer, which is relatively thin (between 20 and 25 cm) compared to the average total water column thickness (about 2 m). The wavelengths were calculated using the nominal frequency of the antenna (80 MHz). Considering the applied velocity models, electromagnetic wavelengths are 0.42 m and 0.75 m for the surveys performed in the lake and on the land, respectively. Therefore, the surveys can resolve layers with a minimum thickness of ca. 0.2 m in the lake and ca. 0.4 m on land. It is important to note that the 0.42 m wavelengths and 0.2 m resolution are under the assumption that microwaves travel through water, and it will not be valid as soon as the microwaves pass the lakebed. Below the lakebed, the wavelengths are comparable to those on land.
SEG-Y files were imported into the Petrel software package (version 2021) for the geological interpretation of radar reflection profiles or radargrams. Vertical datum and sample interval values were manually set for each SEG-Y file. Sample interval values of 0.014 m, 0.04 m, and 0.023 m were set for profiles acquired on water, ice/water, and land, respectively. All radargrams were displayed together in 2D and 3D windows, where a quality control check of the data was performed. Discrepancies of the order of 10 cm in the vertical position of the lake's radargrams were observed. The discrepancies are associated with the time delay generated by the increment of velocity when the electromagnetic waves travel through the ice layer. A static correction (depth shift) in Petrel was performed in order to remove the depth discrepancies.
Two attribute workflows from Petrel were applied separately on the original radargrams to reveal concealed geological features and ultimately assist with the interpretation. The automatic gain control (AGC) attribute is a gain recovery method in seismic analysis that enhances weak signals [62,63]. The attribute equalizes amplitude values by considering the normalized RMS amplitude over time [64]. It was applied to boost weak amplitudes in the GPR radargrams. The cosine of the instantaneous phase attribute is also known as the Hilbert transform. It is calculated from a complex signal known as the analytic signal, where the real part consists of the time-varying radar traces and the complex part is built by applying a 90 • phase shift filter to the real part [65]. It is a non-dependent amplitude attribute that helps to better delineate structural features [64]. In this study, the cosine of the instantaneous phase attribute was instrumental in determining the contact between the bedrock and the overlying sediments.

Data Integration and Digital Mapping
The objective of the data integration is to gather into the same environment (either 2D or 3D), and under the same coordinate system, geospatial datasets, to identify plausible correlations supporting the further reconstruction of 3D surfaces, which ultimately constitute the holistic model. The model surfaces represent geologic units, the lake bathymetry, and the surrounding topography. Geological surfaces resulted from interpolating interpreted geospatial data into a 3D environment.
The datasets were interpreted by performing digital mapping, which is the process of generating vector data (points, polylines, and polygons) that delineate the structure of each geological unit and define the contact between them. Polylines representing geologic features in the lakebed were manually digitized in 2D onto the aerial orthophotos in ArcMap. The units were mapped based on the observed image color and tonalities, and pixel textures. After digitization, elevation values were assigned to the resulting polylines. For example, polylines representing the interception between the water surface and a geological exposure were assigned 35.75 m.a.s.l, the current elevation of the water surface. The same elevation value of the water surface was assumed when interpreting the 1937 image because the depth of the lake is not historically or regularly monitored. Elevation values of polylines representing geological contacts in the lakebed were assigned when integrated with the geophysical data. Not all vector data were generated by manual digitization; this is the case with the interpretation of the bedrock in the land. GOPs reporting the bedrock were filtered out and displayed on top of the 2014 DTM in ArcMap. Contour lines were subsequently generated every 0.5 m from the DTM. Those contours honored by bedrock GOPs were chosen for further modeling.
The 2D digital mapping was also performed on the GPR profiles. The digitization of radar surfaces that ultimately represent the geological contacts between the radar packages was carried out in Petrel. Likewise, reflections depicting the internal geometries or radar facies that constitute the radar package were digitized. It is important to note that only lines representing radar surfaces were used for modeling.
Vector data interpreted in ArcMap were imported into Petrel software, where they were further interpolated with the geophysical and borehole vector data to generate the geological surfaces. The convergent interpolation method was used to interpolate the data. First, the grids were computed at a 1 × 1 m cell resolution within the study area. Where there were no data, grids were extrapolated using the normal extrapolation method, which tended to flatten the surfaces. Subsequently, a low-pass filter was applied to smooth grid nodes. Only one iteration performed over a 4 × 4 cell neighborhood was enough to remove spikes and random noise from the surfaces. Surface-to-surface operations were applied based on superposition and lateral continuity laws to establish a realistic stratigraphic relationship between the units. Finally, to exclude areas of high uncertainty, a polygon that enveloped the concentration of the data was used to clip the surfaces. The

Geotechnical Boreholes
Five lithological units ( Figure 6) were interpreted by combining the geotechnical plots and the sedimentological information of the available samples: (i) bedrock; (ii) till; (iii) fluvioglacial; (iv) gyttja; and (v) infill material. The samples provided hints about the grain sizes and compositions of the sediments, from which three facies were established (C: muddy sandy gravel; F: gravelly muddy sand; and N: fine sandy mud). However, samples were not definitive when establishing clear criteria to differentiate between till and fluvioglacial units. Figure 6 (top) shows how drilling information was interpreted using geotechnical plots of boreholes 3116 and 3402; the boreholes' location is shown in the inset map at the bottom. In both boreholes, the bedrock is represented by the serrate-blocky shape of the FDT curve that records resistance values between 5 and 10 kN. In addition, the top of the bedrock is marked by the symbology used in the total sounding method (see Section 3.3). At the top of the bedrock, the FDT curve dramatically increases, recording the maximum resistance values with a serrated pattern representing till sediments in both boreholes. Fluvioglacial sediments are not present in borehole 3116. Instead, the alternation of troughs and peaks ranging from low to medium resistance values suggests infill material at the top of the glacial till. On the other hand, FDT serrated curves ranging between 5 and 15 kN depict fluvioglacial sediments on top of the tills in borehole 3402. Moreover, very low resistance values (close to zero) suggest poorly consolidated, muddy material overlying fluvioglacial sediments. This material is classified as gyttja, an organic-rich thick mass composed mainly of dead plant remains, clay, silt, and fine sand. Fluvioglacial sediments and gyttja in borehole 3402 suggest subaqueous conditions under fluvial and lacustrine domains. Conversely, the absence of both units in borehole 3116 may be associated with the proximity to the former lake shoreline where subaerial processes were more dominant (Figure 2b). The 12 boreholes at the bottom of Figure 6 are representative results obtained using the same interpretation criteria employed for boreholes 3116 and 3402. All of them were drilled using the total-sounding technique and have geotechnical plots. In addition, 5 of them (3145, 3418, 3430, 3426, and 3424) have particle size distributions and organic matter analyses.

Interpretation of Geologic Observation Points
UTM coordinates of each GOP, together with their lithological classification and geologic measurements, are provided in Table 3. Plotting GOPs on both the 1937 and the 2020 aerial photos (Figure 2) enables the analysis of rock and sediment exposures in a historical context. Geological features observed at each GOP can be visually correlated by comparing their shape and color tonalities on the imagery record. Some rock exposures (e.g., GOPs 2,4,9,and 15) classified in the field as either in situ bedrock or non-in situ erratic block (Table 3) are easily identified in the 2020 RGB photo (Figure 2a,c). Using the GOP coordinates, the same rock exposures can be visually identified in the 1937 black-and-white photo. Such exposures exhibit rounded to sub-rounded geometries with whitish and light grayish tonalities (Figure 2b,d). Another example is GOP 17 (Table 3), where collected clay-rich organic material was classified as gyttja with marshy grass on top. Marshy grass is seen in brownish colors in the 2020 RGB image (Figure 2a), whereas it can be associated with dark gray irregular shapes in the 1937 image (Figure 2b). Location and geometric and color descriptors can be used to cross-correlate geological features through the historic imagery record. In this way, the correlation criteria can be extrapolated to visually classify those geological features that were once subaerially exposed (observed in the 1937 image) but are not visible in the lake today.
Examples of rock exposures (in situ and non-in situ) are shown in Figure 7.  Figure 7a shows evidence of glacial striation, probably caused when erratic blocks were transported by the glacier ground the underlying bedrock. The striation trend (NE/SW) coincides with the regional trend of glacial transportation (Figure 1a). Figure 7b-d show the foliation and internal deformation of the phyllite. The average dip direction (301 • ) of foliation measured in the field coincides with that reported by Birkeland, 1983 [28].
Although layers of till sediments were identified in the boreholes, there is no sufficient evidence of these sediments outcropping within the study area. However, the non-in situ blocks emerging from the lake at GOPs 1, 5-9 ( Figure 2) are interpreted as glacial erratic blocks embedded within glacial till sediments in the lakebed. These are allochthonous blocks of gneiss (e.g., Figures 7f,g and Table 3) transported by glaciers from northeastern Norway [26,33] and deposited in the study area, which is dominated by the Ryfylke phyllite. The source of these blocks may be the middle to upper Scandinavian Caledonides Allochthonous complex or the BCS. The layers of till sediments were further confirmed in the GPR profiles within the lake. Table 3. Summary description of rocks and sediments observed in GOPs. See the location of GOPs in Figure 2 by correlating with consecutive numbers in the first column.

Interpretation of GPR Data
About 8 km of GPR data (available as Supplementary Material) [70] split into 210 radargrams were processed and subsequently interpreted using the proposed radar stratigraphy coding scheme ( Figure 5).
The position and orientation of three radargrams representative of the GPR survey acquired within the lake are shown in Figure 8. The uninterpreted and interpreted radargrams are displayed in Figures 9-11. The GPR depth of penetration within the lake is greater than 2.7 m from the lakebed, which is at 2.5 m depth on average from the water surface. In general, the lakebed is almost planar and horizontal, and the underlying layers of sediments are mostly continuous. The exception is observed near the islet where the lakebed is rough due to rock boulders, and the underlying sediment layers are truncated against the islet. Four GPR packages separated by three GPR surfaces were interpreted in the radargrams within the lake. The GPR packages represent the following geological units: (i) bedrock; (ii) till; (iii) fluvioglacial; and (iv) lacustrine. The lacustrine unit was not identified in the boreholes, whereas the organic matter (i.e., the gyttja unit observed in the boreholes) was not detected by the GPR microwaves.
The performance of the GPR survey on the onshore side was not as expected due to the characteristic of the medium over which the survey was collected. Only the infill material unit was interpreted in a single onshore radargram. The geotechnical boreholes and remote sensors confirm the lateral extension of this unit. A representative radargram of the onshore GPR survey is shown in Figure 12.

Bedrock Unit
The bedrock exhibits different GPR facies depending on the thickness of the overlying sediments. The reflections depict a moderately continuous pattern at shallower areas near the shoreline and the islet, where thin layers of sediments cover the bedrock (Figure 9b). On the other hand, reflections are chaotic in deeper areas toward the center of the lake, where a thicker column of sediments overlies the bedrock (Figure 9b). The bare bedrock exposed at the shoreline exhibits continuous, horizontal, and slightly wavy radar facies (Figure 9b) with very low amplitudes (Figure 9a). However, once the medium changes (i.e., the air-to-ice transition across the shoreline), the reflections become less continuous, and the amplitudes sharply increase towards the lake. Downlap and onlap radar reflection geometries suggest an irregular GPR surface (magenta line in Figure 9b), which represents the top of the bedrock and the geological contact with the overlying unit composed of till sediments.

Till Unit
The radar package representing the till unit exhibits continuous to moderately discontinuous dipping and horizontal reflections, occasionally wavy and sometimes planar. This is interpreted as incipient bedding (Figures 9 and 11). The hyperbolic shape derived from wave diffractions suggests the presence of boulders within the tills (Figure 10b). Tills are interpreted as variable thickness layers of poorly sorted sediments, which were reworked and redistributed by fluvial currents during deglaciation periods. In some places, till sediments are absent due to fluvial erosion (Figure 10c). Erosional truncation on the tills (Figure 9b,c), as well as onlap terminations observed in the overlying sediments (Figures 10b,c and 11b,c), mark the top of the till unit and the bottom of the fluvioglacial unit.

Fluvioglacial Unit
Generally, the fluvioglacial unit is characterized by continuous and horizontal radar reflections (Figures 9-11). The gradual vertical increase of amplitudes from the bottom to the top is interpreted as a coarsening upward stacking pattern (Figures 9a, 10a and 11a). Layers of finer sediments at the bottom are mostly planar, whereas the layers of coarser sediments at the top are observed as wavy radar facies. In the south of the lake, a progradational pattern is observed (Figure 10c). When prograding, layers of fluvioglacial sediments exhibit reflections that gently dip toward the lake associated with cross-stratification; notice the wedge facies (Figures 10b,c). Post-glacial fluvial processes begin by eroding the till unit, followed by the deposition of fluvioglacial sediments along the paleotopography until they fill up the area. Finally, these sediments prograde towards the lake (Figure 10c). Discontinuous and high-amplitude reflections within the bottom of the fluvioglacial package are interpreted as occasional remaining till boulders embedded in fine fluvioglacial sediments (Figure 10a).

Lacustrine Unit
The uppermost radar package is interpreted as layers of lacustrine sediments. The unit is characterized by mostly horizontal, planar, parallel, and continuous radar reflections (Figures 9-11). Lacustrine sediments drape the underlying units and constitute most of the lakebed. The contact with the fluvioglacial unit is mostly characterized by concordant facies, although onlap ( Figure 10c) and downlap (Figure 11c) terminations are also observed. Towards shallower areas, the lacustrine unit is in lateral contact with the till unit (Figures 9c and 11c). Amplitude values progressively increase towards the top, suggesting a coarsening upward stacking pattern. Moreover, amplitude values tend to increase laterally, particularly towards shallower areas. This is interpreted as an increase in the grain size adjacent to the contact with the till unit (e.g., x-x' dots, Figure 9a) (see also Section 4.5).

Infill Material Unit
Although the boreholes confirm the Quaternary sediments and bedrock on the presentday onshore side, the reflections observed in the radargram (Figure 12a,b) do not depict them clearly. These reflections strongly differ from those observed in the lake radargrams (Figures 9-11).   One GPR surface separating two packages of infill material is defined by interpreting together the facies and anomalies observed in the radargram and the position of the DTMs and the well tops (Figure 12a,b). The former topography represented by the 2014 DTM matches with the GPR surface. The reflections observed between the 2014 and 2019 DTMs represent the upper package, which exhibits GPR facies that vary from continuous on top, with higher amplitudes, to moderately continuous at the bottom with low amplitudes (Figure 12a,b). The top of the lower infill package is clearly defined by the 2014 DTM, from which the 3406 and 3404 boreholes were drilled. However, the definition of the lower package's base and lateral extension is unclear. From the right end (south) to the center of the radargram, a ca. three-meter-thick package, represented by planar, parallel, continuous, dipping, and very low amplitude GPR facies, is observed (Figure 12a,b). Below this package, facies are chaotic, presumably corresponding to gyttja, according to borehole 3406. Vertical changes of GPR facies may correlate with the position of the top of the bedrock and the lower package of infill material observed in borehole 3404 (Figure 12a,b). GPR facies are moderately continuous above the bedrock borehole top, whereas facies are wavy and continuous below the well top. This observation is confirmed in Figure 12d, where the 1937 aerial image shows the position of borehole 3404 adjacent to the exposure of the bedrock before being buried by the infill material.
The strong reflection objects observed in the center of the radargrams are interpreted as pipelines buried into the upper and lower infill packages (Figure 12a,b). The low-amplitude smiling artifact that represents the pipeline embedded in the upper package exhibits a lowfrequency signal signature. This anomaly may also be interpreted as a metallic pipeline or electric wires. Similarly, the pipeline within the lower package exhibits a concave-upward semicircular artifact, though its amplitudes and frequencies are higher.

Aerial Imagery Enhancing
By visually comparing the digital orthophotos year by year, major geomorphological terrain changes can be observed from the 1930s to the present (Figure 2). They depict the transition from a rural to an urban landscape, where the construction of civil infrastructure has greatly altered the lake footprint and its vicinity along with the geological record. A historical aerial imagery record was instrumental in the following aspects: (1) the detection of noticeable human-induced morphological changes along the shoreline through time; (2) the identification and analysis of geological features that were once subaerially exposed and subsequently buried by the artificial shoreline progradation; and (3) characterization of the underwater sediments laying in the lakebed. Steps to achieve results on the latter aspect are explained in detail in Section 3.1. The results of the sequential processing of the 2007 RGB aerial imagery are shown in Figure 13. The geological interpretation of the resulting imagery is presented in Section 4.5.  On top of the imagery, colored polylines outlining the interpreted units within the lake were drawn. Since 1937, the shoreline has artificially prograded about one hundred meters to the west. Since then, exposed geological features, such as the bigger islet close to the eastern shoreline, have been buried.

Integrated Data Interpretation within the Lake
The till unit (outlined by blue lines, Figure 14) is represented by orange to reddish colors in the 2007 imagery. It exhibits a coarse texture interpreted as poorly sorted sediments with the presence of boulders, which sometimes correspond to erratic blocks as confirmed by GOPs 6-9, for instance. The contact (magenta line) between the bedrock and the glacial till is observed along the shoreline, especially on the northern side. This contact is also observed around the islet and in the underwater mound about 70 m west of the islet.
Lacustrine sediments (outlined by green lines, Figure 14) vary from aquamarine to olivine-green colors and have a smoother texture that may suggest fine grains. Color changes may suggest grain size variations across lacustrine layers. This observation is supported by the DAT_0009 radargram (Figure 9a). Color changes marked in the image along the DAT_0009 survey line (x, x', and x" black dots, Figure 14) were plotted in the radargram (x, x', and x" black dots, Figure 9a). x and x" indicate the contact between the till and lacustrine sediments. The x-x' segment, representing olivine-green sediments, matches stronger GPR amplitudes ( Figure 9a) and is, therefore, interpreted as coarser lacustrine sediments. Aquamarine sediments in the segment x'-x" are interpreted as finer sediments and are represented by low amplitudes in the radargram (Figure 9a). Other internal variations across the lacustrine unit were observed in the 2007 imagery. At the eastern shoreline, note the dark-colored apron-shape feature, interpreted as lacustrine fine sediments combined with organic matter. The east apron end seems to have been buried by the progression of the shoreline. The organic matter consists mostly of vegetation remains and depicts different fashions. Organic matter concentrations are outlined by dark red polylines northeast of the lake (Figure 14). In the 1937 photo, dark grayish features that emerged from the lake with irregular shapes were interpreted as benches of gyttja with marshy grass on top. This correlates with what is observed in GOPs 17 and 18. Black speckles and amorphous-shaped, dark bodies randomly distributed across the lakebed are also interpreted as accumulations of organic matter in the colored 2007 photo. Finally, along the eastern shoreline, yellowish tonalities outlined by gray polylines were classified as the infill material that humans have historically added to artificially prograde the shoreline to the west (Figure 14).

Onshore Integrated Data Interpretation
Based on the rock exposures observed at the GOPs 11-16 and 19 (Figure 15), bedrock is the only unit that crops out on the onshore side, particularly on the east where the Våland high stands out from the surrounding topography. The solid black contours extracted from the 2014 DTM represent the relief expression of the bedrock, whose minimum elevation on the surface is 39 m and its highest point is above 56 m a.s.l. Information from the geotechnical boreholes at the western flank of the Våland High suggests that the bedrock deepens about 6 m under the 2014 surface and then flattens towards the lake. In this area, borehole tops of bedrock average 30 m a.s.l, though occasionally small bedrock highs, such as the islet stick out from the water level according to GOPs 2 and 4 and the digital mapping carried out on to the aerial imagery (magenta polylines, Figures 14 and 15). The bedrock also rises above the water level towards the north and south of the lake, as confirmed by the rock exposures ( GOPs 15,16,and 19) and bedrock tops recorded in the groundwater boreholes, respectively. Groundwater boreholes -NGU ! .

Holistic Model of Mosvatnet
The stratigraphic interpretation of Mosvatnet is ultimately illustrated by the 3D holistic model (Figure 16). The NW-SE perspective view of the 3D model provides a first overview of the position of the bedrock, the lateral distribution of the sediments in the lake, and the spatial relationship between the natural and non-natural infill materials (Figure 16a). The four main geologic units (bedrock, till, fluvioglacial, and lacustrine) are displayed along with the RGB-textured 2019 DTM, which represents the study area's present-day topography. A more detailed representation of the lateral and vertical distribution of the geologic layers and the morphology on the onshore side are shown in the NW-SE cross-section of Figure 16b. Differences in the geometry of the top of the bedrock and the overlying sediments are observed on both sides of the islet. Towards the northwest of the islet, where the bedrock is flatter and shallower, the layers are dominantly horizontal and exhibit a thinner and more constant thickness. Conversely, towards the southeast, the bedrock top is deeper and exhibits a depression between the islet and Våland High. The depression was filled by glacial, fluvioglacial, and lacustrine sediments. In this area, the layers are less horizontal and thicker, and the thickness is more irregular relative to the other side of the islet.  (Figure 16a,b).
Structural maps depict the morphology of each geologic unit across the study area ( Figure 17). The top of the bedrock (Figure 17a) is interpreted as the paleotopography over which the glacial sediments were deposited. The bedrock depicts the most significant morphological changes. It records the highest and lowest elevations in the entire study area. The elongated depression striking NE-SW corresponds to that basin observed in the cross-section (Figure 16b). It is interpreted as a glacial incision generated by the abrasion of the glacier upon the bedrock. The morphological expression of the depression is inherited by the glacial till sediments, which partially fill the basin (Figure 17b). In general, the till sediments are well distributed over the bedrock. Blank patches (holes) in Figure 17b indicate areas where till sediments were either not deposited or eroded by subsequent fluvial processes. The distribution of fluvioglacial sediments over tills is less homogeneous (Figure 17c). Their deposition is more or less restricted to the elongated basin, where fluvial currents incised, reworked, and transported the till sediments, ultimately redeposited as fluvioglacial sediments. Finally, lacustrine sediments homogeneously drape the underlying units, except for highs, such as the islet (Figure 17d). Lacustrine sediments mostly make up the lakebed, which depicts a very gentle, almost flat bathymetry.

Interpretation of Mosvatnet Model
Mosvatnet was initially a confined paleo-depression overdeepened by sub-glacial erosion exerted by glaciers moving from the northeast on the underlying bedrock. The depression was filled with subglacial melt-out till sediments during the period spanning from the last glacial maximum to the deglaciation around 12.9 ky-11.5 ky BP [18]. After deglaciation, glacial meltwater channels reworked, eroded, and redistributed till sediments within the lake. Since then, fluvial processes have contributed mainly to filling Mosvatnet with sediments. Besides erosional and depositional processes acting within the lake, external fluvial currents have transported and dumped fluvioglacial sediments on top of the till sediments within the lake. For example, the small creek in the southern shoreline, between the groundwater boreholes and GOP 10 ( Figure 2b) is responsible for the progradational structure observed in radargram DAT_0002_4 (Figure 10). Similarly, input sediments that consist of silt and clay have been transported in suspension everywhere within the lake and gradually settled to form layers of post-glacial lacustrine sediments that drape the tills and fluvioglacial sediments. Lacustrine sediments primarily constitute the lakebed away from the lake shoreline. Coarser lacustrine sediments are derived from the constant erosion exerted by waves and internal lake currents on pre-existing sediments close to the shoreline. For example, the olivine-green coarser lacustrine sediments resulted when waves reworked till deposits in the northern shoreline ( Figure 14). Finally, post-glacial organic matter (gyttja) accumulated as isolated patches on top of the lacustrine sediments, mostly in the northeast of the study area.
The description of the bedrock unit observed in the GOPs on the onshore side and the islet correlates well with the Ryfylke phyllite reported by previous authors [26,[28][29][30]. There are discrepancies between the geological units interpreted in this study and those reported by the surficial deposit map generated by the NGU (Figure 1b). The discrepancies are mainly observed in the eastern onshore side, where anthropogenic activities have greatly modified the landscape through time by adding infill material for construction over the glacial Quaternary deposits. In this area, NGU reports thick and thin moraine material. The Våland High is interpreted as thin moraine material or bedrock covered by thin moraine material. However, based on field observations from outcrops and the morphology exhibited by this structure in the ALS-derived DTM (Figure 15), the Våland High is interpreted as bare bedrock in this study. Finally, GOPs collected in this study confirmed that the bare bedrock constitutes the islet, conversely as proposed in the NGU map, where the islet is interpreted as part of an accumulation of organic matter or peat (gyttja in this study).
With the urbanization of Stavanger, natural materials have been removed or covered by non-natural infill materials, especially in the eastern shoreline, where the correlation between GPR and multitemporal remote sensing enables the identification of anthropogenic events that occurred in different time windows. The events suggest the development of the E39 motorway in two epochs, represented by the lower and upper infill packages mapped in the onshore GPR (Figure 12b): Epoch 1: artificial progradation of the eastern shoreline in about one hundred meters westward between 1937 and 2014 (Figures 14 and 16b); and Epoch 2: the prograded area increased its elevation by about eight meters relative to the lake level between 2014 and 2018 due to the construction of the Mosvatnet culvert. The two pipelines buried within each package confirm the two anthropogenic epochs.
The pipelines were confirmed by comparing the raw and processed GPR data ( Figure 12). The employed velocity model properly migrated the primary reflections of the bedrock, sediments, and infill material. Nonetheless, the model produces an over-migration effect on the pipelines, resulting in smile artifacts from the multiples of the hyperbolas.

Advantages and Limitations of the Employed Techniques and Model Uncertainty
None of this study's techniques can be used as a single solution to optimally map natural and non-natural features in both the surface and subsurface. However, these techniques are complementary when combined into a workflow to generate a holistic model that integrates anthropogenic and geological features in an urban landscape. Gaps in knowledge due to a technique's limitation are filled by the capabilities of other technique(s). Borehole drilling and GPR surveying can collect geospatial data in the subsurface where surficial methods, such as geological field mapping and remote sensing cannot, and viceversa. For example, by combing the enhanced orthophoto and GPR radargrams, it was possible to map in 2D the geological units within the lakebed and in the subsurface, respectively, and ultimately generate a 3D characterization of the lateral and vertical distribution of such units.
GPR surveying and remote sensing are non-invasive techniques during the acquisition of data. Both techniques are well adapted to urban environments and deliver centimeterscale resolution and accuracy datasets. As in non-urban areas, non-invasive geological field mapping may also provide valuable data in urban environments. On the other hand, geotechnical drilling is an invasive technique. It may destroy the sediment records and hampers the integrity of the samples.
Human-induced changes in the landscape, particularly within or adjacent to urban areas, are remarkable and happen too quickly relative to morphological changes derived from natural processes. Assisting in detecting and mapping these changes through time is another advantage of remote sensing. The analysis of ALS-derived DTMs and systematically collected aerial imagery over time has enabled the identification and mapping of not only the areas where Mosvatnet has undergone significant anthropogenic changes but also the natural areas that have remained stable. Furthermore, characterizing morphological changes through time helped to understand better why the capabilities of methods such as GPR are limited when applied to map certain materials under certain medium conditions. For example, the strong signal attenuation observed in the GPR survey collected over the infill material at the eastern onshore side of Mosvatnet ( Figure 12). The GPR signal attenuation can be associated with the presence of silt, clay, and moist saline conditions in the studied terrain [10]. Conversely, the GPR survey conducted offshore was able to image the bedrock top and resolve the centimeter-scale structures of the overlaying sedimentary layers. Even better results were obtained for radargrams collected over the frozen lake.
Limitations of the employed methods may increase the uncertainty in the geological interpretations and hinder the degree of confidence in the Mosvatnet model. Other aspects intrinsically related to the employed methods, such as type, quantity, coverage, and distribution of the collected datasets, may also affect the geological uncertainty. Based on the degree of confidence, the data can be divided into hard and soft data. This study considers hard data to be direct observations and measurements made on the outcrops and samples collected from the geotechnical boreholes. In contrast, indirect measurements, such as ALS-derived DTMs, aerial imagery, and GPR, are considered soft data.
Combining hard and soft datasets may help reduce the model uncertainty, especially when geographically overlapped. For example, GOPs helped validate that the Våland high can be easily mapped from the DTM corresponds to the bedrock unit ( Figure 15). Besides being sparse, another limitation of GOPs is that not all the units can be observed in the exposures for detailed descriptions. In this study, bedrock is the only unit fully exposed in outcrops. Exposures of allochthonous glacial erratic blocks emerging from the lake are just part of the available tangible evidence suggesting a till unit in the study area. Additional evidence derived from radargrams was required to determine more unit constituents and confirm the existence of the glacial till in the subsurface within the lake.
Hard and soft datasets derived from drilling and GPR surveys provide information about the subsurface. However, the overlapping between these datasets is minimum ( Figure 2) since borehole data are restricted to the present-day onshore side, whereas highquality GPR data are restricted to the lake. In this scenario, the data combination does not help reduce the model uncertainty on both the onshore and offshore sides. It is important to note that the temporal aerial imagery has shown that geotechnical boreholes were drilled in an area originally occupied by the water body. This is why the stratigraphy observed (below the E39 motorway) in the geotechnical plots and sediment samples correlate well with that interpreted in the radargrams within the present-day water body lake. The onshore hard data from geotechnical boreholes were instrumental in the understanding of the stratigraphy of Mosvatnet. However, discerning the vertical position and composition of lacustrine sediments was not possible in any of the boreholes. The employed drilling technique may not be optimal for preserving the integrity of thin layers of lacustrine sediments, which may be too soft, poorly consolidated, and mixed with either gyttja, infill material, or both. The model surface representing the lacustrine unit underneath the infill material ( Figure 16b) is a mere extrapolation of the interpreted surface in the lake where the unit is well resolved by GPR and mapped on the radargrams. Hence, the lacustrine unit on the onshore side is the modeled surface with the highest geological uncertainty.

Conclusions
This study employed a qualitative methodological approach to combine remote sensing, geotechnical data, geological field mapping, and geophysical tools to create a holistic 3D model in a highly urbanized area. After unifying the coordinate reference system, the datasets derived from the aforementioned methods were interpreted separately using geological and geophysical criteria. Firstly, the interpretation enabled discernment between natural and anthropogenic materials. Secondly, the interpretation criteria allowed us to differentiate and characterize the geological units within the natural materials. All datasets and their respective interpreted vector entities were integrated into a common 3D space to find patterns that allow the correlation between materials from the onshore, offshore, surface, and subsurface.
None of the methods employed offers a complete solution when mapping natural and anthropogenic materials. However, the methods complemented each other well. Their integration proved to be more efficient and enabled the generation of the 3D holistic model, which provided a broad understanding of the subsurface geology and the induced anthropogenic changes in the area through time.
Both groundwater and geotechnical boreholes provided an initial understanding of the stratigraphy on the present-day onshore side. GPR data interpreted within the lake using the seismic stratigraphy approach supported estimating radar facies' lateral and vertical distributions, which correlated well with the stratigraphy observed in the boreholes. Moreover, GPR interpretation tied well with the lithological contacts in the lakebed, which were mapped onto the aerial imagery. Due to the high energy absorption produced by the dielectric and conductivity properties of the non-natural infill material, very poor ground penetration and sediment imaging were obtained in the GPR survey conducted on the eastern onshore side.
The holistic model allowed the reconstruction of the depositional environment of glacial and post-glacial units. Firstly, the glacial ice sheets moving from the northeast abraded and polished a paleo-depression on the Ryfylke phyllite bedrock over which subglacial melt-out till sediments were deposited. Subsequently, glacial meltwater channels eroded the glacial tills within the depression during the glacial retreatment and deposited sandy-rich fluvioglacial sediments on top. Simultaneously, fluvial sediments transported by external currents entered the lake and contributed to building up the fluvioglacial unit. External input consisting of silt and clay sediments moved in suspension over the lake and settled to form thin layers of the lacustrine unit draping the underlying units and forming the lakebed. Finally, accumulations of post-glacial organic matter (gyttja) have developed and were distributed as patches, primarily in the northeast of the study area.
Onshore GPR and multitemporal remote sensing enabled the detection of prominent human-induced topographic changes in Mosvatnet. These changes are split into two epochs of anthropogenic activities associated with the development of the E39 motorway. The eastern shoreline artificially prograded about one hundred meters toward the west, and the corresponding area increased its elevation by about ten meters relative to the lake level in 82 years.
This study demonstrated that natural city lakes, such as Mosvatnet, with relatively undisturbed conditions on and under its lakebed, may act as a window to the past, enabling deciphering the areal and vertical distribution of subsurface structures. Investigating the Mosvatnet lakebed provides an excellent opportunity to access geological information that would otherwise be impossible to obtain. Future work will include filling the gaps in geophysical information: (i) collecting additional GPR surveying in the offshore site, preferably conducted in winter when the lake is frozen; and (ii) conducting other geophysical investigation methods, such as seismic refraction tomography and ERT, might help to improve the subsurface knowledge in the onshore side and contribute to improving the geological confidence of the model in this area.