Bedrock surface topography of Latvia

ABSTRACT A new map and digital bedrock surface elevation model of Latvia is presented with a horizontal resolution of 250 m. The local bedrock comprises largely undisturbed layers of Palaeozoic and Mesozoic sedimentary rocks covered by up to 200 m thick Quaternary strata composed mostly of glacigenic and marine sediments. The bedrock surface elevation model and the corresponding map were created by tessellation interpolation of the bedrock surface based on more than 20,000 boreholes and constrained by the present land surface. The known locations of paleo-incisions buried under Quaternary cover are indicated. The new digital map has applications in hydrogeological, environmental, and civil engineering, and fundamental and applied research.


Introduction
The bedrock surface, as a major sedimentary discordancy, has important practical and scientific applications ranging from prospecting of mineral resources (Paulen et al., 2006), mapping karst-sensitive areas (Chen et al., 2017), spatial planning of development, paleological reconstructions (Kalm & Gorlach, 2014), and hydrogeological investigations (Virbulis et al., 2013). The bedrock surface represents a radical shift in physical conditions that have shaped the composition and texture of sediments and sedimentary rocks.
In Latvia, bedrock topography maps have been produced routinely as a part of a set of maps from geological mapping campaigns at a scale of 1:50,000 (Aleksans, 1991;Aleksans & Ginters, 1988;Drikis, 1980;Gavrilova & Ginters, 1975Ginters et al., 1985Ginters et al., , 1986Jansons et al., 1965Jansons et al., , 1967Jansons et al., , 1971Juskevics & Vihots, 1978;Murnieks & Murniece, 1979;Podgurskis et al., 1985;Tracevskis, 1989Tracevskis, , 1993Tracevskis et al., 1984;Ulgis, 1981Ulgis, , 1983 and 1:200,000 (Freimanis et al., 1966;Gavrilova et al., 1962Gavrilova et al., , 1966Gavrilova et al., , 1967Gavrilova & Straume, 1963;Jankins et al., 1969Jankins et al., , 1973Jankins et al., , 1975Kajak, 1976;Mironovs et al., 1962;Tracevskis et al., 1964Tracevskis et al., , 1965Tracevskis et al., , 1967Tracevskis et al., , 1969Vetrennikovs & Bendrupe, 1976) during the second part of the twentieth century. The latest edition was elaborated in the 2000s by reviewing the historical materials (Brangulis et al., 2000;Juskevics et al., 1997Juskevics et al., , 1998Juskevics et al., , 1999Juskevics et al., , 2003Misans et al., 2001Misans et al., , 2002Murnieks et al., 2002Murnieks et al., , 2004. These maps have been published only on paper and split as separate sheets. Nominally the scale of these maps corresponds to the scale of mapping; however, the actual resolution based on input data was probably poorer. In addition, almost two decades have passed since the compilation of the latest bedrock maps, and a wealth of new geological well data has accumulated. Since the end of the 2000s, approximately 22,000 borehole descriptions have been supplemented with more than 5000 new records. We have produced a new map of the bedrock surface topography of Latvia. We use the reported bedrock surface depth compiled in the national geological well database (Takčidi, 1999) and tessellation based on a Voronoi diagram as an interpolation grid. The map is complemented with reconstructed locations of buried valleys. It is available as a map completed with reference grids and a digital elevation model of the bedrock topography. It can be used as an input for a range of studies, for example, a recent investigation of groundwater pollution with agrochemicals (Kalvāns et al., 2021).

Geological setting
Latvia is located in the northwestern corner of the East European plain within the last Scandinavian ice sheet (Zelčs et al., 2011). The modern topography was mostly shaped by the multiple advances and retreats of the Scandinavian Ice sheet during the Pleistocene and modified by the more recent action of the Baltic Sea, fluvial and other terrestrial processes. The landscape can be described as a combination of plains corresponding to ice streams or ice flows and hilly uplands (Kalm, 2012). The present surface elevation ranges from a few metres below to 312 m above sea level ( Figure 1).
The Quaternary sediments are mostly glacial, glaciofluvial and glaciolimnic sediments (Dreimanis & Kārklin, 1997). The Quaternary cover is thinnest in lowland areas, often just a few metres but rarely absent, and thickest in upland regions, particularly in the central and eastern part of the country, reaching 200 m in the uplands (Zelčs et al., 2011).
At the bedrock surface, most Devonian sedimentary rocks are exposed except for the southwestern corner where more recent sedimentary rocks of Paleozoic and Mesozoic are present (Lukševičs et al., 2012). Generally, the sedimentary sequence has close to horizontal bedding with a few dislocations (Brangulis & Kaņevs, 2002). The Paleozoic and Mesozoic sedimentary strata are slightly inclined towards the south and southwest. As a result, in the northern part of the country, middle Devonian sandstones, siltstones, and clays, as well as dolomites, are exposed at the bedrock surface. In the central and southern regions, upper Devonian carbonate and mixed carbonate and terrigenous sedimentary formations and Mesozoic, mostly terrigenous sediments are exposed.
The largest (up to hundreds of kilometres) features of the bedrock surface topography largely reflect surface topography featuresglacial uplands and lowlands (Gorlach et al., 2015). These features reflect denudation by major ice flows of numerous advances of the Scandinavian Ice Sheet during the Pleistocene. The same ice flows redeposited the material eroded upstream within the glacial uplands, masking the underlying bedrock from further glacial erosion (Kalm & Gorlach, 2014). On a kilometre scale, bedrock surface topography is determined by the bedrock composition. Cuesta-like escarpments or cliffs can be identified in places where hard Devonian dolomites are poorly consolidated terrigenous sediments.
There is widespread evidence of an existing network of paleo-incisions that are often referred to as buried valleys in the study region (Meirons et al., 1974;Rattas, 2007). These are deep and narrow valley-like incisions in the bedrock surface often filled with glacial, glaciofluvial and glaciolimnic deposits, sometimes distinguishable in the modern topography. According to the geological borehole database, the paleo-incision depths can exceed 300 m (average 70-80 m). The width of these features is from 10 to 20 m to a few kilometres, and they can be traced from a few hundred metres to up to 30 km. Initially, it was proposed that the buried valley systems are relict interglacial or pre-Quaternary river incisions (Skuodis, 1959). Later it was suggested that some of the incisions formed as tunnel valleys below the Scandinavian Ice Sheet by high-pressure subglacial streams (Vegt et al., 2012). In addition, it has been suggested that the pre-existing valley system was enhanced and modified by subglacial erosion (Vaher et al., 2010).

Borehole data
A geological borehole data set evaluated by Popovs et al. (2015), originally obtained from the Latvian Environment, Geology and Meteorology Centre and updated with the latest information (as of September 2020), was used. The dataset includes stratigraphic records from a borehole database covering the territory of Latvia (Takčidi, 1999). As supporting information, borehole data from similar databases for Estonia (Estonian Land board, 2015) and Lithuania (Belickas, 1999) were used. Lithuania and Estonian datasets were limited to a 50 km buffer zone along the Latvian border. Bedrock surface elevation from 17,481 wells in the Latvian mainland and 3071 wells from Lithuania and Estonia were used. No data from bordering areas of Belarus and Russia were available.

Data evaluation
Comprehensive quality control of the good dataset was carried out considering errors of recorded location and unrealistic stratigraphic sequences. The viability of borehole vertical reference was checked by a simple topological rulethe borehole cannot deviate more than the maximum possible vertical error against the Earth's surface for the national digital elevation model with a 20 m resolution (denoted maximal error in areas with poor ground surface visibility) (LĢIA, 2017 ). Such a large error margin accommodated decreased DEM accuracy along steep slopes of river valleys (Zhu et al., 2005). Unrealistic stratigraphic sequences are recorded in boreholes where either the order of stratigraphic units is mixed or where there are contradictions of entered thickness or depth values of stratigraphic units, resulting in a corrupted stratigraphic column and questionable legitimacy of the full borehole log. Such problems were screened algorithmically by checking the ascending sequence of the stratigraphic descriptors (Popovs et al., 2015). A total of 604 boreholes (2.94% from the initial dataset) were excluded from the dataset.

Buried paleo-incisions
Buried paleo-incisions pose a significant challenge for constructing a map of bedrock surface elevation. Knowledge about their distribution and morphology is limited and spatially uneven (Gorlach et al., 2015). Wells recovering buried incisions can be identified by comparing the Quaternary sequence with other nearby wells, but the incision morphology and extent remain unknown. Regardless of bedrock surface interpolation methods, since an insufficient number of wells within valley structures is available, holetype effects will be observed around such wells. Therefore, the dataset was screened through cross-correlation and comparison of borehole descriptions of nearby wells, where 1223 wells (5.95% from the initial dataset) identified as recovering buried incisions were excluded from the interpolation of bedrock surface topography.

Surface construction
A multi-step euclidean distance tessellation method was used to construct a sub-Quaternary surface, forcing the bedrock surface to be equal to or less than the modern land surface elevation (LĢIA, 2017). This methodology assumes that the bedrock and land surface topography is uncorrelated except at the locations where Quaternary sediments are absent.
Standard geostatistical methods or simple linear interpolation were not suitable for the construction of a digital elevation model of the bedrock surface due to the spatial clustering of boreholes and the heterogeneity of the pre-Quaternary surface itself.
The tessellation was based on standard Thiessen polygons (Thiessen, 1911), widely used in many fields of modelling (Bodin et al., 2012;Ju et al., 2011). For each point, a polygon is constructed with a numeric value of the input where the input point serves as the centroid (the arithmetic centre of mass) (Figure 2(A)). In such a way, each of the two nearest polygons shares one edge section. However, we expand it further by calculating the mean value between each two nearest polygons and generating ). This approach bypasses the geometric effects of linear interpolation when applied to irregularly scattered data, gives a visually realistic surface, and produces a good fit against input data. Five-step tessellation was performed where the average area of the resulting polygon approximately coincides with the chosen resolution of the target surface ( Figure 3).

Target grid resolution
The chosen nominal grid cell size of the target model was 250 m. The target grid resolution was determined following the approach described by Hengl (2006). Accordingly, the sampling density of the input dataset corresponds to 205 m, while the finest legible grid resolution (a finest meaningful resolution that represents the most (95%) of surface and spatial objects) was 246 m. The density of available data varies radically, with the highest density along densely populated areas (Figure 4(A)). The median value between the nearest neighbouring wells is 492 m. In certain areas, the resolution could be increased; however, in large areas with sparse borehole coverage, the map resolution should be coarser than this. Therefore, the chosen nominal resolution corresponds to the rounded finest legible grid resolution proposed by Hengl (2006) for the regions with the largest density of boreholes. We maintain the overall finer resolution to preserve the spatial information where it is available.

Results
The presented map (Main map, Figure 5) and bedrock surface digital elevation model provide a new digital dataset for future geological investigation and wide application in geological and hydrogeological modelling. Implementing the dataset in models will help advance geological research and integrate geological and ecological fields, providing the knowledge base for investigating groundwater-dependent ecosystems and interconnected networks of water flows at the surface and subsurface.

Uncertainty assessment
The uncertainty of the map was assessed in locations where more than one observation well was present within one 250 m grid cell. On average, the observed difference in bedrock surface elevation within a single grid cell mean error was 1.1 m ( Figure 6). This is the average uncertainty of this map at grid cells with observation points. The average difference between interpolated and actual (unknown) bedrock surface elevation in grid cells with no observation points is likely to be higher.

Data
The presented map is prepared as a formatted printable map (Main map) and georeferenced dataset comprising four raster bands: elevation of the bedrock surface above sea level, locations of boreholes used for the construction of the elevation of the bedrock surface, locations, where buried valleys in geological boreholes were uncovered, and thickness of Quaternary sediments.

Discussion
We have developed a new digital bedrock surface topography map of Latvia. In this discussion, we draw the attention of the reader to potential applications of the map and highlight some of the known issues. We use the case study of Kazu Leja (Kalvāns    (Juškevičs, 2000;Krievāns, 2015). It comprises a few interesting bedrock surface topography features: a cuesta-like escarpment of Devonian dolomites, partly buried paleo-incision and nearby river valleys as well as springs emerging from a karst aquifer. A 3D geological model of the Kazu Leja is available at the web page https://www. gprm.lu.lv/kazugrava-3d/

The modern river valleys
Many modern river valleys are incised in the bedrock surface, and in most cases they are well represented in our bedrock surface model (Figure 7). However, the methodology of representing these valleys in the bedrock surface elevation model has inherent flawsthese were cut into the interpolated surface of the bedrock by the intersection with the land surface model. Therefore, this approach disregards any alluvial sediments present in the river valley. Partly, the significance of the problem was curtailed by the relatively coarse resolution of the modelmost of the river valleys were barely represented at 250 m resolution. Nevertheless, it remains an issue for the largest rivers like Gauja and Venta.

The cuesta-like features and steep escarpments
The cuesta-like features and steep escarpments are among the geological features best represented in the presented bedrock surface map as long as these are expressed at the land surface. The nature of their  formationthe erosion of soft basement rock forming a pair of a steep slope and a gently inclined plateaurenders them suitable for representation on the map. At the first stage of map production, the smooth bedrock surface was interpolated between existing observation points. This was followed by the second stage, where the surface was forced to be at or below the land surface model (DEM), forming the escarpment front. The case study site, Kazu Leja, is located at the distribution margin of the dolomites of the Upper Devonian Pļaviņas formation (D 3 pl), whih forms an escarpment or cuesta-like feature projected in the modern topography (Figure 8).

Buried valleys (paleo-incisions)
Detailed mapping of paleo-incisions remains a major problem. On the main map, we show the paleoincisions mapped before at a scale of 1:200,000 (Brangulis et al., 2000;Juskevics et al., 1997Juskevics et al., , 1998Juskevics et al., , 1999Juskevics et al., , 2003Misans et al., 2001Misans et al., , 2002Murnieks et al., 2002Murnieks et al., , 2004 as well as locations where these incisions are recovered in boreholes (Main map). However, examination of borehole log data shows numerous locations where the geological record suggests paleo-incisions, but their configuration is difficult to interpret at the presented mapping scale. We did not attempt to re-interpret previously mapped buried incisions as the resolution of available data was insufficient.
In the case of Kazu Leja (Figure 8), a buried incision was uncovered in a borehole at -30 m a.s.l. filled with till and sand and gravel deposits (Bendrupe & Arharova, 1981). The trajectory of the incision is marked by modern topography and proved by geological outcrops. However, apart from a single borehole, there were no data on its depth profile. A mixed fluvial and subglacial origin of the paleo-incision was proposed (Vaher et al., 2010). A fluvial valley is expected to have a smooth longitudinal profile, making possible interpolation between observation points. However, if the dominant form of their formation was pressurized subglacial tunnel valleys (Vegt et al., 2012), then their longitudinal profile is likely to be variable, complicated with sills and troughs; thus, extrapolation of a single observation along an imagined channel axis is likely to be speculative. Thus, the issue of paleo-incisions remains unresolved.

Applications
In the study region, the bedrock surface is the major discordance in the geological record, representing a sedimentary hiatus of up to about 380 million years from Middle Devonian to Quaternary. Any geological model needs to account for the striking differences in sedimentary environment and processes during the Paleozoic and Quaternary, delineated by this sedimentary hiatus. Thus, the presented map will be a useful input for future work in geological modelling in the region.

Maps and models
The bedrock surface elevation map can be useful in geological construction models and deriving new high-resolution bedrock geology maps from these models. Pre-Quaternary strata formed in mostly shallow epicontinental seas (Lukševičs et al., 2012), with the sub-horizontal arrangement of the geological strata disrupted by a relatively small number of faults and low-amplitude folds (Brangulis & Kaņevs, 2002). These properties of regional bedrock geology were exploited by Virbulis et al. (2013) and Popovs et al. (2015) to create regional geological and hydrogeological models using interpolating marking layers among geological boreholes. Combining the presented map of bedrock surface topography and the geological modelling approach of Popovs et al. (2015) opens opportunities to produce high-resolution bedrock surface geological maps by intersecting the two data products. Furthermore, a bedrock surface lithology map can be produced by integrating the lithological descriptions from well databases into geological models and intersecting the result with bedrock surface topography. Such a map is currently is unavailable for the study region but would have many applications, for example, for geological prospecting and hydrogeological and ecological studies. Such a map would mimic the currently available Quaternary geology maps that integrate stratigraphical and lithological features.

Digital elevation analysis
The analysis of digital elevation maps is an established field in Earth sciences; however, it is primarily restricted to land surface topography. However, such features as coastline escarpments into bedrock, cuesta-like features along the margins of bedrock, carbonate formations, and elevated bedrock foundations underlying some of the uplands are exhibited in the bedrock topography map (see cross-sections in the main map). The resolution of the presented map does not match that of current land surface models, but it enables the rudimentary analysis of the bedrock surface elevation for geomorphology research.

Hydrogeology
Hydrogeological research can benefit from the availability of the presented bedrock elevation map. The groundwater flow is constrained by the underlying geology setting. For example, karst systems developed in carbonate aquifers can quickly funnel surface contaminants to the streams (Bakalowicz, 2005;Stevanović, 2018). To manage these hydrological systems properly, ensuring good groundwater and surface water status, detailed knowledge of the overburden and outcropping areas of these aquifers need to be known. Furthermore, the data must be presented in a consistent digital format, enabling seamless integration into modelling systems. It has been shown that detailed mapping of groundwater denitrification potential in cultivated lands can be a cost-effective measure in designating relevant, site-specific fertilizer application restrictions (Jacobsen & Hansen, 2016). Location of the interface between Quaternary and bedrock sediments can be an important parameter for assessing denitrification potential due to their contrasting hydrological properties or guiding the estimation of the depth of the redox interface (Hansen et al., 2014), i.e. the depth where anaerobic conditions facilitate denitrification.
Bedrock aquifer intersection by the land surface in the study region often supports groundwater-dependent ecosystems (Kalvāns et al., 2021;Koit et al., 2021) that are resilient due to consistent water supply but which might be vulnerable if groundwater suffers overexploitation or pollution. The presented map can be a valuable tool for predicting the locations of groundwater-dependent ecosystems and elaborating the conceptual hydrogeological models for such sites. Thus, the bedrock elevation map can be a useful input for water resources research.

Comparison of historic and new maps
A comparison was made between the newly generated map with newest historical map in scale 1:200,000 (Brangulis et al., 2000;Juskevics et al., 1997Juskevics et al., , 1998Juskevics et al., , 1999Juskevics et al., , 2003Misans et al., 2001Misans et al., , 2002Murnieks et al., 2002Murnieks et al., , 2004. The elevation isolines of the historical paper bedrock surface elevation map were digitized and a DEM created by the spline interpolation method. The historical DEM then was compared to geological well data available at the time when the historical map was finalized, i.e. before 2000 and well data established after 2000 and thus not used to create this historical map, and to the new map presented in this article ( Figure 9).
The comparison demonstrates that in all three cases, the error statistics are rather similar ( Figure  9). The historical map is relatively robust, capturing all larger-scale features of the bedrock surface and accuracy of the map at the locations of new observations, which although relatively poor, have approximately the same accuracy as the accuracy at locations of observations available when the map was created. The error statistic between the new map and the old map is comparable to the old map and geological wells. The algorithm used for the creation of the new map that ensures almost a one-to-one match between the observation point and mapped elevation.
Overall, the new map is superior to the historical one as it better uses borehole data and better represents more local scale features allowing the new DEM model to be used not only at the regional scale but also in more local studies. While both maps give a similar insight to general characteristics of the bedrock surface, a large deviation of the historic map from borehole data for more than ±100 m limits the practical use of this map.

Conclusions
A bedrock surface elevation map for Latvia was elaborated using a geological borehole database complemented by known locations of buried paleoincisions. The map is a useful data source for geology and hydrogeology-related research. The map is based on a bedrock surface digital elevation model that can be used as an input to modelling exercises to solve problems related to the exploration and management of geological resources and ecosystems.

Software
The Main Map was produced using ArcGIS Pro software from ESRI. To prepare the regional cross-sections, QGIS has been used. Final editing and PDF creation were carried out with Adobe Illustrator. Various map images of the manuscript were made using R-project (R Core Team, 2020).

Disclosure statement
No potential conflict of interest was reported by the author(s).

Data availability statement
A developed digital relief model of the Bedrock surface of Latvia used as a basis for the Bedrock surface map is available within the article's supplementary materials. Restrictions apply to the availability of the primary data which were used under the license for this study. The primary data are property of relevant state authorities and can be requested from the Latvian Environmental, Geology and Meteorology Centre (www.meteo.lv), Geological Survey of Estonia (www.egt.ee) and Geological Survey of Lithuania (www.lgt.lt).