Investigating peatland stratigraphy and development of the ˇ Sijec bog (Slovenia) using near-surface geophysical methods

Owing to their anoxic environment, peatlands play an important role in the preservation of records documenting past atmospheric depositions. To determine past records, data on peat stratigraphy and bog development are needed. Ground penetrating radar (GPR) was used to determine the peat thickness and morphology of the ˇ Sijec bog on the Pokljuka plateau in Slovenia, which will serve as a basis for further geochemical studies. Information on the stratigraphy below the peat/clay boundary was acquired by applying electrical resistivity tomography (ERT). The GPR results reveal four depressions within the peat bog, which are separated by elevated ridges. Within the depressions the peat reaches a depth of 6 – 9 m. The edges of the bog are flat, with peat thickness ranging from 2 to 4 m. The reach of the GPR was complemented with manual peat probing. A comparison of the depths obtained using GPR and the peat probe reveals that the results of both methods correspond well in most locations. The ERT indicated similar peat depths; peat responds with high electrical resistivity. In contrast, clayey sediments with low resistivity are found below the peat. The peat depressions are underlain with larger clayey depressions reaching more than 20 m in thickness and represent lake sediments. The complementary geophysical methods proved to be an efficient approach with which we can delineate the peat morphology and the under- lying stratigraphy. Both indicate bog formation from a lake with four deeper depressions, that are separated by glacial deposits. The results presented here show the potential for geophysical methods to infer formational processes in peatlands, showing the presence of a series of isolated basins that later coalesced into a single peat landform. This interpretation is consistent with previous conceptual models from studies in boreal regions.


Introduction
Most recent peatlands developed after the retreat of the glaciers and represent archives of past atmospheric depositions owing to the waterlogged anoxic and anaerobic conditions there (Charman, 2002). Ombrotrophic peatlands, commonly referred to as bogs, are domed peatlands that rise above the surrounding soils and above the influence of groundwater (Strzyszcz and Magiera, 2001;Cocozza et al., 2003). They receive water exclusively by precipitation (Clymo, 1987;Biester et al., 2014) and predominate in the northern hemisphere. In Europe, larger peatland areas are found in Northern Europe, Russia, Great Britain, and Ireland (Tanneberger et al., 2017). Peat thickness in ombrotrophic bogs has been successfully studied using GPR (e.g. Slater and Reeve, 2002;Lowry et al., 2009;Rosa et al., 2009;Parsekian et al., 2012), because the penetration depths are greater than those in minerotrophic peatlands (Malterer and Doolittle, 1984). This is due to the fact that ombrotrophic peatlands receive only precipitation inputs, which result in lower pH values and lower basic cation concentrations (Ca, Mg, Na, K) (Doolittle and Butnor, 2009).
Geophysical methods, particularly ground penetrating radar (GPR), hold particular potential to support the understanding of the peatland's stratigraphy, development and hydrology. Over the last 20 years, geophysical methods have been used increasingly in peatlands studies to determine the thickness of peat, estimate the volume, and to evaluate gas accumulation levels (Comas et al., 2005a;Sass et al., 2010;Kennedy et al., 2018). A major advantage of GPR compared to electrical resistivity imaging (ERT) is the ability to rapidly collect a large number of profiles and the good resolution of the data which is needed to detect the morphology of the peat layer. On the other hand, the results gathered using the ERT method provide information on the stratigraphy of the bog below the peat.
Peatlands are generally not ideal environments for GPR due to the higher water content, which can strongly attenuate the radar signal, but the low electrical conductivity and changes in water content and bulk density between the peat and mineral sediments can provide strong reflections at such boundaries (Slater and Reeve, 2002;Lowry et al., 2009). GPR is useful where changes in moisture content or bulk density are characteristic enough to generate measurable reflections (Doolittle and Butnor, 2009). In general, changes in the dielectric constant between the peat and the underlying sediments (gyttja, clayey sediments, sand) provide an interface that is clearly evident in the GPR data. Furthermore, changes in the moisture content of the peat and mineral soil results in large-amplitude reflections at the peat-mineral soil interface (Theimer et al., 1994). Some researchers (e.g. Meyer, 1989;Slater and Reeve 2002;Comas et al., 2004;2005b;2015;Kettridge et al., 2008;Walter et al., 2016;Legchenko et al., 2011) have applied, apart from the GPR method, induced polarization imaging, electromagnetic terrain conductivity and/or electrical resistivity imaging as alternative tools with which to determine peat thickness and stratigraphy. Even though geophysical methods are convenient as they are non-destructive, geophysical results combined with conventional methods (peat probing, coring) allow a more complete stratigraphic picture of the subsurface of peat deposits (Rosa et al., 2009). The bogs in Slovenia belong to the southernmost ombrotrophic peatlands in Europe (Kutnar, 2000). Peatlands on the Pokljuka plateau formed in a karst environment where surface water is rarely present. Nevertheless, in one area of the plateau, several smaller peatlands are present. One of the larger ones, Š ijec bog, shows a domed structure, indicating an ombrotrophic nature of the peatland. The bog surface is completely waterlogged in some areas, especially in the northwest area, which makes it almost impassable most of the year. Therefore, the use of geophysical methods in winter was used to test the condition needed to study Š ijec bog and to determine the formation and morphology of the bogs in the karst environment overlain by glacial sediments. The Š ijec bog has a rugged surface morphology and ridges covered with small bushes in its centre. Geophysical methods provide information on the stratigraphy, thus the relationship between surface and subsurface morphology and the reasons for the present morphology, as well as the development of the bog can be studied. Moreover, Š ijec bog has not yet been studied regarding morphology, development, hydrology (water inputs and outputs of the peatland, as well as hydrology of the surrounding area) and composition of peat organic and mineral matter. Based on the results of the presented study, we propose a conceptual model that follows on previous models for peatland initiation in domed Šijec bog with measured profiles. The GPR profiles mentioned in this study are marked in red and the ERT profiles in yellow. Arrows mark the direction of the profiles. Point data represent peat probing locations along the presented GPR profiles. Location of core data (white square) is from Budnar-Tregubov (1958). Basemap data acquired from LiDAR ARSO (http://gis.arso.gov.si/evode/profile.aspx?id=atlas_voda_Lidar@Arso). Photos of measurements of Profile 1 using the GPR and ERT method show the different terrain conditions needed for an optimal profiling of these methods. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) bogs (Comas et al., 2004). This study provides useful information for potential locations of further detailed geochemical and hydrological research, which can be used for a more accurate determination of the peatland type. In addition, the presents the first geophysical investigation of peat bogs undertaken in Slovenia. GPR is used to determine the thickness of the peat layer and simultaneously test the GPR's potential to analyse the variability of said thickness and to estimate the extent of the peat. The ERT is used to complement the GPR data on the peat/clay boundary and to determine the stratigraphy below the peat layer.

Study area
Pokljuka is a karst plateau in northeastern Slovenia. The Pokljuka plateau displays a flat topography at an altitude of 1200-1500 m and is surrounded by valleys that have been reshaped by glaciers in the south, southwest and east, while in the northwest it borders on higher mountains. Due to the karst environment, surface water can be found in many karst springs and many peatlands. Among them, the Š ijec bog, located 1194 m above sea level, is the largest and covers approximately 15 ha. Due to its altitude, the bog has a colder climate than the Slovenian average. Its average temperature in January and February 2020 was -2.5 and -0.8 • C respectively (measured at the Rudno polje meteorological station 5 km from the Š ijec bog), while the average temperature for January 2020 throughout Slovenia is 1.9 • C (ARSO, 2020). Precipitation in Pokljuka totalled 19.5 and 49.9 mm in January and February 2020 respectively (measured at the Rudno polje meteorological station). Š ijec bog is part of the Triglav National Park.
The Š ijec bog has a relatively flat topography with elevated edges, the highest of them along the southern border, where it reaches an elevation of 1198 m. The bog's lowest point lies in its north-western part, where a larger depression is apparent, at 1191 m above sea level (Figs. 1 and 2). This area has higher water content throughout the year due to many Sphagnum pools in the area, similar to the pools described in Comas et al. (2005b). The elevated areas are the driest and are covered with pine bushes. Flat areas are largely covered with Sphagnum moss and grasses. In the central part of the bog there are water pools and hollows, surrounded by small drier ridges that are covered with smaller pine bushes. The pools are relatively small, measuring approximately 4 m in length and 50 cm in depth. These shallower bodies of water are present throughout the year and freeze completely during the winter. The Š ijec bog is surrounded by dense spruce forest (Andrič et al., 2010).
Based on one borehole located roughly in the middle of the bog, the boundary between peat and lake sediment is found at a depth of 6 m, where the lake sediments are radiocarbon dated to 11,000 years (Maja Andrič, personal communication, 2019). Additional two cores are presented in Budnar-Tregubov (1958), one of which is described in detail. In the described core, peat is present to a depth of 4.7 m, and is followed by a mixture of peat and clay. At 5 m depth clay predominates, but still has a higher organic content. Sediments, similar to lake sediments, occur at 6 m depth. Finally, sand is reached at 8.7 m depth, attributed to moraine material. According to Jurkovšek (1987), glacial till sediments and Triassic carbonate rocks represent the bedrock of the Š ijec bog. Š ifrer (1952) constructed a map of Pokljuka with an emphasis on Quaternary landforms with a focus on the extent of the last ice age on the Pokljuka plateau. According to him, moraines from the last glaciation are located at the eastern and southern edges of the Š ijec bog, while moraines from older glaciation can be found in the western part, and partly in the peat bog area. The results of previous work suggest that after the last glacial maximum, when the Pokljuka plateau was presumably glaciated, a Lateglacial and early-Holocene lake formed at the study site (Andrič et al., 2010).

Methods
Since the Š ijec bog is part of the Triglav National Park, we used the non-invasive methods to determine the subsurface characteristics of the bog. In the presented study we used ground penetrating radar (GPR), electrical resistivity tomography (ERT) and manual peat probing. In order to cover the entire bog, 16 GPR, 2 ERT profiles were measured and 30 manual peat probe measurements along the GPR profiles were made (Fig. 2). As weather conditions limit the application of applied methods, we performed the geophysical procedures in winter, when the peat was frozen (maximum depth of 0.3 m). In contrast, peat probing is only feasible when the peat has thawed and was thus performed in late spring. The field conditions at the time of the GPR, ERT and peat probe measurements are summarized in Table 1. The temperature range for 10 days before the GPR and ERT survey is similar, but the conditions on the days of measurement indicate milder weather at the time of the ERT measurements, which can also be clearly seen from Fig. 2.

Ground penetrating radar
Ground penetrating radar (GPR) is a non-destructive, established method that emits short pulses of electromagnetic (EM) waves into the subsurface and detects reflected signals at the interfaces between materials with different electromagnetic properties (Blindow et al., 2007). In addition to various applications, GPR is also used for studies in peatlands. Compared to traditional surveying methods (peat probe, coring), GPR is faster and requires significantly less time and effort to obtain extensive and continuous information on the thickness, volume, and geometry of peatlands (Doolittle and Butnor, 2009). Even though EM wave velocities in peatland environments are low due to the high water content of the peat, changes in the dielectric constant at the boundary with underlying sediments and the low electrical conductivity of the water generally provide good results (Lowry et al., 2009;Sass et al., 2010). Significantly lower water content at the peat/mineral sediment interface results in high-amplitude reflections. Generally, the interface between the peat and the underlying material is clearly visible, allowing for continuous tracking of the basal of the peat, which provides detailed information on the depth and volume of it (Slater and Reeve, 2002;Sass et al., 2010;Kennedy et al., 2018;Ryazantsev and Mironov, 2018). GPR has also been used in hydrogeological studies of peatlands to determine the moisture content in the peat (Karušs and Bērziņš, 2015) together with water distribution (Legchenko et al., 2011), to study subsurface lithological and hydrological conditions, to delineate the water pathways in peatlands (Trappe and Kneisel, 2019) and to locate springs and ponds above the break and slopes to predict groundwater flow (Lowry et al., 2009). GPR has been used to detect oilfield pipelines in peatlands and identify subsurface piping (Smith and Jol, 1995;Holden et al., 2002), to study permafrost peatlands (Kaverin et al., 2018), estimate carbon stocks and peat soil characterization (Comas et al., 2015). Several studies have evaluated approaches used to estimate peat depth (Parry et al., 2014), study uncertainties in estimating peat volumes using GPR and probing (Parsekian et al., 2012), and to determine the number of manual measurements required to minimize EM wave velocity error (Rosa et al., 2009). Along the GPR, peat probing and coring have largely been used in the mentioned studies to calibrate and compare peat depth and stratigraphy with radar reflections. In general, lower-frequency (<200 MHz) antennas are typically used to profile peatlands (Doolittle and Butnor, 2009). Factors that have to be taken into consideration when choosing the type of GPR system are penetration depth and the resolution required to achieve a specific goal. In our study, the Mala ProEx GPR common-offset survey method was used with two different antennaean unshielded 50 MHz RTA (Rough Terrain Antenna) and a shielded 250 MHz antenna. At the study site, the 50 MHz antenna provided better results due to the better penetration depth required to detect the boundary between the peat and the underlying sediments. The penetration depth of the 250 MHz antenna proved insufficient to determine peat thickness at the deepest areas. Furthermore, in rough terrain, the unshielded 50 MHz RTA flexible (9.25 m-long tube) antenna generally proved to be more suitable than the rigid shielded 250 MHz antenna due to its easier maneuvering.

GPR survey, data processing and manual peat probing
The GPR measurements were performed during the winter, when the upper surface layers were frozen and equipment could be maneuvered (Fig. 2). During the other seasons, surface water pools and Sphagnum pools with a very high water content render GPR measurements more difficult. Profiles were acquired in longitudinal and cross lines in order to cover the entire surface of the bog. The GPR profiles were adjusted to the field conditions owing to the dense shrubbery inside the bog (Figs. 1 and 2).
We measured 13 profiles with the 50 MHz antenna (4 profiles, also in reverse direction) and 3 profiles using the 250 MHz antenna (Fig. 2). Due to vegetated and uneven terrain, all 13 profiles with a 50 MHz antenna were measured with a "hip-chain" encoder with a rotating wheel for data acquisition, which proved to be the only appropriate way on more challenging terrains such as karst environments (Čeru et al., 2018;Č eru and Gosar, 2019).
All GPR profiles were processed and analysed using Reflex-Win software by Sandmeier Geophysical Research. The following processing steps were employed to all profiles: editing, subtract-mean (dewow, time window: 20 ns), time zero correction (correct max. phase: 46-52 ns, move startime), background removal, gain correction (energy decay, scaling value: 0.2), bandpass filtering (low cut 25 MHz, low pass 50 MHz, high pass 150 MHz, high cut 300 MHz) and static correction (topography). The final process applied to all profiles was a 2D filter (xy average), which suppresses trace and time dependent noise and acts as a lowpass filter (Sandmeier, 2011). The processing steps applied to all profiles are shown in the case of Profile 1 (Fig. 3).
Time-depth conversion of the GPR data was also applied. Converting the two-way travel time into the depth to determine the velocity in the peat can be done in different ways: (1) using GPR with CPM (common midpoint surveys) or WARR (wide angle reflection and refraction) technique; (2) by measuring at locations where the depth of the peat was determined using the peat probe (manual probing); or (3) by hyperbolae fitting. All of these techniques and methods have their advantages and limitations and produce some error in estimating the peat depth due to the spatial variability of the water content, the porosity of the peat, bulk density and presence of woody layers (Parry et al., 2014). Based on the reported results of various studies, the velocity of peat as determined using various techniques ranges between 0.033 and 0.050 m/ns, with the most common velocity being around 0.035-0.038 m/ns (Legchenko et al., 2011;Plado et al., 2011;Sass et al., 2010;Parry et al., 2014;Comas et al., 2015).
Due to the fixed separation of the GPR system used in our study, the EM wave velocity could not be determined using a CMP survey. For converting time to depth, we used the average velocity v = 0.035 m/ns evaluated from the hyperbolae fitting, confirmed also by peat probe results. Similar values have been determined in studies of Sphagnum peat in other studies (e.g. Comas et al., 2005b;Parsekian et al., 2010;Plado et al., 2011). However, the application of a single velocity in GPR data results in some error in the estimation of peat depth. Therefore, we applied two different velocities (v = 0.035 and v = 0.045 m/ns) to show the difference in peat depth (Fig. 3d). The velocity of 0.045 m/ns was used as the maximum value to evaluate the difference in depth and the maximum error that can be made in depth conversion in peat material. The difference between the selected (v = 0.035) and the average lowest value of the dielectric constant for peat (v = 0.045) is less than 1 m (Fig. 3d).
Due to the winter conditions (frozen peat) during the geophysical investigations, the peat probing could not be applied at the same time, but was performed later in late spring, once the peat layer had thawed. Some published studies reported a correlation between GPR depth and manually probed depths. Parry et al. (2014) found that manually probed depth values are on average 35% shallower than those derived using GPR. The reasons for these variations may be that manual probes do not always reach the peat-mineral interface and the methods used to calibrate EM wave velocity may also introduce errors (Parry et al., 2014). Rosa et al. (2009) concluded that calibration can be more difficult when the peat-mineral contact is particularly irregular, and statistical analysis shows that at least 30 calibration points (data from the peat probe) are required in order to minimize EM wave velocity error when estimating peat thickness using GPR. Parsekian et al. (2012) revealed that GPR has great potential for estimating the depth and volume of the peat basin considering sitespecific peat properties with levels of uncertainty similar to those associated with invasive direct probe methods. Peat volumes estimated using GPR results were about 23% lower than the estimates made using the set of probe measurements (Parsekian et al., 2012). These results also reveal that a comparison of GPR and probe data corresponded well in regard to peat depth, although the deeper measurements saw greater variability, most probably due to the spatial variability of peat and its properties.
Peat probing in Š ijec bog was performed in May, when the ground had completely thawed. A peat probe was used to measure the depth at 30 locations (Fig. 2) where the probe could be easily pushed to the point of resistance contrast. Several measurements had to be taken at each location in order to reach the gradually stronger resistance at a depth corresponding to the GPR estimates. The peat probe often came to a sudden halt, probably when it encountered wood fragments. This indicates that although the peat probe can provide accurate results, several measurements are required to reach the peat/clay boundary, so a prior depth estimate from GPR data is particularly useful.

Electrical resistivity tomography
Electrical resistivity tomography (ERT) is a non-destructive directcurrent geophysical method that uses multi-electrode systems together with the computer control unit and metallic electrodes to measure ground potential (e.g. Reynolds, 2011;Lowrie, 2007;Loke et al., 2013). For a 2D survey a number of evenly spaced electrodes is placed into the ground surface in a straight line. The measurements are made by introducing an electrical current into the ground through two current electrodes and measuring the difference in the resulting voltage at two potential electrodes. By using different electrode spacings at different locations along the cable, a 2D profile of the subsurface is obtained with the selected electrode arrangement. The choice of electrode spacing with a known number of available electrodes defines the resolution and depth range of the ERT method for a particular electrode arrangement. The current and voltage measurements are then converted into an apparent resistivity, taking into account the geometric factor that depends on the configuration of the electrodes, which in turn depends on true subsurface resistivity distribution . True resistivity distribution in the investigated medium is estimated using an inversion procedure based on minimizing a suitable function (Loke and Barker, 1996;Loke and Dahlin, 2002).
The electrical conductivity of rocks and sediments generally depends on the amount of water in the medium, the conductivity of the water, and the water distribution (porosity, the degree of saturation, cementation factor, fracturing) (Kowalczyk et al., 2014;Glover, 2015;Merritt et al., 2016), while also taking into account the lithology (Dobrin and Savit, 1988;Telford et al., 1990). According to previous studies, characteristic resistivity values for peat range from 200 to 500 Ωm, occasionally less (e.g. Legchenko et al., 2011;Comas et al. 2015); clay up to 100 Ωm (e.g. Telford, et al., 1990;Chambers et al., 2013;Kowalczyk et al., 2017), till sediments up to 1200 Ωm (e.g. Pellicer and Gibson, 2011;Pellicer et al., 2012;Dietrich and Krautblatter, 2017) and carbonate bedrock up to several thousand Ωm (e.g. Telford et al., 1990;Pellicer et al., 2012Pellicer et al., , Ž ebre et al. 2019. Successful applications of the ERT method for the characterization of peat bogs include determination of the thickness of the peat layer above the sediments (Comas et al, 2004;Sass et al., 2010;Comas et al. 2011;Kowalczyk et al., 2017;Illés et al., 2019) and the saturation within the peat layers (Legchenko et al., 2011), while at the same time resolving the underlying lake sediment and glacio-marine clay thickness, as well as variability in depth to glacial till, eskers and the underlying bedrock (Comas et al. 2004;Sass et al., 2010;Comas et al., 2011).
The ERT measurements in the Š ijec bog were conducted at the end of February, when the temperature reached 8 • C and the upper 5 cm of the peat was partly thawed, which enabled setting of the electrodes. In order to estimate the thickness of the clayey sediments below the peat layer  and to obtain information on the deeper geological composition of the Š ijec bog we have measured two 2D ERT profiles (Fig. 2), for which we used ARES equipment with 48 electrodes. The dipole-dipole array (DD in continuation), which is more sensitive in detecting vertical structures such as depressions and yields better resolution images compared to other conventional arrays (Dahlin and Zhou, 2004) was chosen as the primary data collecting array at both profiles. To achieve the desired investigation depths, we applied 2 m electrode spacing at profile ERT-1, resulting in 94 m in length and an investigation depth of approx. 17 m, and larger 3 m electrode spacing at profile ERT-2, resulting in 189 m (using the roll-along technique) and an investigation depth of approx. 27 m. Considering the challenging conditions for ERT measurements (partly frozen, partly thawed surface, and areas covered with bushes and related organic litter), a Wenner-alpha array (WA in continuation) was also used for control measurements. WA is much less noise-sensitive compared to DD, and better at resolving horizontal layers (Dahlin and Zhou, 2004), while reaching shallower depths. The WA array was applied with the same electrode spacing as DD at both profiles, for the whole length of the ERT-1 profile and the last 141 m of the profile ERT-2.
The resulting data was processed using Res2DInv inversion software (Loke, 2017). Before processing, we integrated topography data and manually removed the bad data points in the dipole-dipole data sets. In ERT-1 we removed data with standard deviation above 5%, and in ERT-2 data with a standard deviation above 8%. All pseudosections were inverted using the l1 norm smoothness-constrained Gauss-Newton leastsquares optimization method, where the absolute difference (or the first power) between the measured and calculated apparent resistivity values is minimized (Claerbout and Muir, 1973), also known as the blocky inversion method (Loke et al., 2003).

Peat basin
The selected profiles and their main interpreted features are shown in Figs. 4-6. The profiles in each figure are presented as two-way travel time sections and a depth scale. Velocities in peat derived by fitting hyperbolas are in the range 0.034-0.038 m/ns, and velocity used for converting time to depth was the average velocity 0.035 m/ns. Comparison of GPR and peat probe data is represented on radargrams (Figs. 4-6) and in graph in Fig. 7. Some hyperbolas that occurred in depressions and at the undulating reflector at 1-3 m depth yielded higher velocities, between 0.06 and 0.09 m/ns. These velocities belong to a mixture of peat, clay, silt and sand. However, some error with respect to depth is inevitable due to changes in water content in the peat and heterogeneous material within depressions.
The GPR results of 13 profiles reveal the heterogenous morphology of the peat deposit. A dense network of GPR profiles (Fig. 2) reveal four deeper bowl-shaped depressions. Strongly contrasting differences in moisture content at the peat/lake sediment interface provide a distinguishable reflector on all radargrams (Figs. 4-6). At the peat/lake sediments boundary, the signal is rapidly attenuated due to the presence of clay, as determined in the borehole in the middle of the bog (Budnar-Tregubov, 1958;Maja Andrič, personal communication, 2019). In all profiles, this interface forms a conspicuous reflector that varies in depth from about 2-10 m and represents the basal surface of the peat layer. No reflections occurred below the peat/clay boundary. Similar results were also obtained in studies where clayey sediments ("gyttja") occur beneath the peat and the reflector is clearly expressed and signal is rapidly attenuated in clay layer (Slater and Reeve, 2002;Sass et al., 2010;Plado et al., 2011;Comas et al., 2015).
Depressions are marked with numbers 1, 2, 3 and 4 in the radargrams. They were determined by combining GPR results and shown as a 3D model (Fig. 8). The location of the ridges separating depressions corresponds to older moraines deposited after earlier glaciation periods (Šifrer, 1952), indicating that the depressions may have formed prior to the last glaciation. Similarly, glacial deposits within peatlands have been found in other peatlands (e.g. Comas et al., 2004;Chen et al., 2020). The deeper depressions are of different sizes: two larger ones (160-190 m in diameter) in the western part (depressions 1 and 2 in Figs. 4 and 5), and two smaller ones (60-100 m in diameter) in the eastern part (depressions 3 and 4 in Fig. 6). Between them are areas with 4-6 m of peat, which are also areas with denser vegetation and higher bushes. Besides the strong reflections at the interface, some stronger reflections are visible within the depressions, which could be related to increased mineral matter, which consists of sandy, silty and clayey material deposited with surface runoff from the surrounding area during earlier stages of peatland development. The shallower discontinuous undulating reflector at depths of 1.5-3.5 m is also visible in some profiles (Figs. 4 and 5). Some hyperbolas occur at this reflector as well at the bottom of depressions. The velocities evaluated from the hyperbolae fitting are in range from 0.06 to 0.09 m/ns, which corresponds to different materials such as wet clay, sand and silt in the literature (Jol, 2009). This may be due to increased sediment input from the surrounding area or climatic changes and related vegetation and water content changes during a global dry event. It may also represent the boundary between minerotrophic and ombrotrophic peat (when the moss grew a dome above the surrounding soils and the surface water runoff lost its connection) or the effects of peat erosion. A reflection-free and weaker reflection areas are probably related to the homogeneous parts of the peat layer. Some diffractions and ringing effects are visible in the upper parts of the peat due to frozen water in open pools (Profile 1 at distance 188-192 m) and Sphagnum pools (Profile 2 at distance 50-60 m, Profile 6 at distance 172-182 m and Profile 9 at distance 80-90 m).
Two cores (III, VI) presented by Budnar-Tregubov (1958) can be used to better describe the main stratigraphic layers. The locations of the cores are not located exactly on the GPR profile, but in its vicinity (Fig. 2). The core described in Budnar-Tregubov (1958) and shown in Profile 1 (Fig. 4) is 872 cm deep and reaches a grey sand which they described as resedimented moraine sediments. Described core is located above the ridge in Profile 1 (marked III). Peat reaches 4.7 m and is followed by peat mixed with green clay, which transits to only green clay at 5 m depth. This correlates well with the lowest reflections in GPR profile. Grey clay, resembling lake sediments as describes in Budnar-Tregubov (1958), is found at 6 m depth. GPR was unable to detect differences between the two clay types at this location. The other core (marked VI) was not described in as much detail but gives information that peat is also present to a depth of 4.7 m and is later mixed with clay, which lower transits to grey clay. On the GPR profile this change can correlate with stronger reflections within the depression. This indicates that the reflections within the depressions are indeed associated with higher clay content. It is not clear whether sand was also found in the core VI. However, depth at which sand occurred in core III would also correlate with the depth of the medium resistivity body at the edge of the peatland in ERT 1 (Fig. 9).
The results of the peat probing were compared with geophysical data and correlations were established between the depth values obtained from the probing and the GPR reflections. On the elevated areas and edges, the depth determined using the peat probing and GPR is similar, with deviations of approximately 10%. In depressions, the depth values obtained from GPR data are deeper than the probing depths (Fig. 7). The differences between measurements in depressions ranges from 15 to 40%. Variations are greater where the probe only reached reflections that belong to heterogenous material (depression 1 in Profile 1, depression 2 in Profiles 2 and 13; see Figs. 4 and 5). The peat probe did not reach the maximum peat depth within depression, which we explain by the fact that at certain depths there is a mixture of peat, clay, silt and sand that disabled the penetration to the clay layer (interface between peat and lake sediments derived from GPR data). Due to the greater resistance associated with peat probing, the reflections are probably caused by a higher clay and sand content in the peat. At the boundary the transition is gradual, which is apparent from increased probe resistance as the clay content increases. Results of the peat probing indicate that an appropriate velocity for peat (v = 0.035 m/ns) was selected. The results of both approaches, conventional and geophysical, are consistent, as has been reported in some studies (Parsekian et al., 2012;Parry et al., 2014). The results of the study showed that GPR can be applied to determine peat thickness with great precision and can produce even more comprehensive data in areas with heterogeneous material where peat probing is not possible. However, it is important to consider the change of material when selecting the velocity, as part of deviations is also a consequence of the uncertainty regarding the applying single wave velocity used in GPR data.
The thinnest parts of the peat layer are located in the areas outside the depressions and mostly represent the edge zones of the peat bog (Profile 10, Fig. 6). The GPR results reveal the irregular basal surface of the peat, followed by clayey sediments that rapidly attenuate the signal. Location of the depressions and ridges between them, corresponding with glacial deposits from previous glacial periods, indicate that it is possible that four separate smaller lakes existed during the last glacial period and were merged at the end of the period. However, it would also be possible that lake formed after the last glaciation and had four deeper centres separated by older glacial deposits. In this regard we conducted ERT survey to determine the thickness and morphology of the lake sediments below the peat. Similarly, Comas et al. (2004) described a model of peatland formation in which the basin was divided by an esker, with pools forming above the esker. In our case, pools, although much smaller, also formed in the vicinity of the glacial deposit ridge, while the deposit itself represent the driest areas covered with bushes.
Comparison of the surface morphology (Figs. 1 and 2) with subsurface morphology of the peat basin (Fig. 8) shows a good correlation between glacial deposit ridge and drier shrub-covered areas in the central part of the peatland, as well as drier areas with thinner peat layer. The northwest area of the peatland, where the surface is wettest, is located above the deepest depression. Contrary to the western peatland area, where there are deeper depressions and wetter surface areas, the eastern surface area is slightly elevated and drier despite the underlying depression. This suggests that the peatland is unlikely to be uniform in nature, with the eastern area presumably being ombrotrophic, while northwest area is probably minerotrophic. Additional research is required to properly determine the development stage of the peatland.

Stratigraphy below the peat/clay boundary
Although the surface conditions during the measuring were challenging, the ERT survey provided sufficiently good data to explain the stratigraphy of the bog. The comparison of the inverse resistivity models shows minimal variations in the resistivity distribution obtained with dipole-dipole (DD) and Wenner-alpha (WA) arrays (Figs. 9 and 10) and is consistent with the sensitivities of the arrays described in Section 3.2.
The resistivity distribution of the ERT profiles ( Figs. 9 and 10) shows   three well-defined resistivity bodies that can be assigned to different geological layers. A sharp resistivity contrast is present between the uppermost peat layer (P) and underlying clayey deposits (C, D1, D2) at both ERT locations. The resistivity contrast is less sharp and occasionally gradual (particularly in the lateral direction) between the clayey deposits and the medium resistivity body (MRB). The uppermost peat deposits (P) are characterized by high resistivity values (200-600 Ωm) on both ERT profiles. The peat layer gradually thickens with a slight surface subsidence, from 2 to 6 m in an easterly direction in profile ERT-1 and from 3 to 9 m in a north-westerly direction in profile ERT-2, forming two peat depressions (1 and 2 in Fig. 8), which are also interpreted using GPR profiles .
Underlying, predominantly clayey deposits (C, D1, D2) with significantly lower resistivity values (20-70 Ωm) form two deep clayey depressions (D1 and D2 in Figs. 9 and 10), indicating the existence of the deeper lakes in the past (greater than10 m in ERT-1 and greater than 17 m in ERT-2), positioned below the two peat depressions. D1 and D2 might well be connected to one larger clayey depression, which cannot be sufficiently explained based on the two surveyed ERT profiles. Considering the low resistivity values, D1 and D2 generally contain saturated clayey sediments, probably with a certain amount of silty and sandy fractions. The peat, clay, and sandy clay deposits of lacustrine origin of the Pleistocene and Holocene periods were identified in previous research (Budnar-Tregubov, 1958;Jurkovšek, 1987;Maja Andrič, personal communication, 2019). Low resistivity clayey sediments (C) are also present between the peat layer (P) and the medium resistivity body (MRB), reaching a thickness of 4-7 m with a depth of up to ~10 m in profile ERT-1 and a thickness of 5-11 m with a depth of up to ~15 m in profile ERT-2, suggesting the later transition from deeper lakes (or a single lake) to a shallower, extended lake. Given the considerable thickness of the clayey deposits in D1 and D2, it is very likely that the deep lakes (or lake) were formed sometime in the Pleistocene. However, it is not possible to distinguish between Pleistocene and Holocene clayey deposits based on the ERT results, and further coring and dating is required.
The body with medium resistivity values (MRB: 100-200 Ωm) is found at a depth of ~10 m in profile ERT-1 and ~15 m in profile ERT-2 and corresponds with the location of the ridges between the peat depressions observed on the GPR results (Fig. 11). According to previous geological (Buser, 1975;Jurkovšek, 1987) and geomorphological mapping (Šifrer, 1952) of the Pokljuka plateau, glacial deposits (moraines) of different ages and Triassic limestones are present on the surface in the vicinity of the Š ijec bog. The resistivity values of the medium resistivity body (MRB) range from 100 to 200 Ωm, which, according to some authors, can be attributed to fractured and saturated limestones (e.g. Telford et al., 1990;Reynolds, 2011;Redhaounia et al., 2016). However, far more often typical resistivity values for various limestones are significantly higher in Slovenia and elsewhere, ranging from at least 500 Ωm to several thousand Ωm, with an average of more than 1000 Ωm (e. g. Pellicer et al., 2012;Ndougsa-Mbarga et al., 2014;Horn et al., 2018;Ž ebre et al., 2019). Resistivity values within the interval 100-200 Ωm are often associated with unconsolidated and to some degree saturated Quaternary sediments such as till, diamicton, sand and gravel (e.g. Beresnev et al., 2002;Kilner et al., 2005;Chambers et al., 2013;Kowalczyk et al., 2017). The latter is also supported by the gradual lateral transition from medium resistivity values observed in the MRB to the lower resistivity values observed in D1 and D2 present in both ERT profiles, which may reflect a gradual transition from more coarse grained sediments (i.e. glacial deposits: 100-200 Ωm) to finer (sands, silts: 70-100 Ωm) and finest (clayey: below 70 Ωm) sediments.
The morphology of peat deposit in both GPR and ERT results is similar. Both methods present a gradual thickening of peat deposit in the depressions (Fig. 11). However, the exact depth of the peat/clay interface slightly differs. The low resistivity within the lowest parts of the peat deposit is very likely connected to the higher amount of wet clay present within the lowest part of the peat layer, which significantly lowered the resistivity values and can also be explained with a decrease of ERT resolution with depth and a gradual change in electrical conductivity but was still detected by the GPR. This is consistent with peat probing results (see Figs. 4 and 5, Profile 1 and 2) and also confirmed with cores described in Budnar-Tregubov (1958), where the peat reaches 4.7 m depth and is followed by a mixture of peat and clay.

The Š ijec bog development
Based on integrated results (GPR, ERT, manual peat probing) and the description of the previously published core data (Budnar-Tregubov, 1958), we propose a model describing the development of the Š ijec bog (Fig. 12). The model present two separate basins in the western part of the peatland. Similar development model where two basins coalesced into a single peat landform was already presented in Comas et al. (2004). Although ERT survey did not reach sufficient investigation depth to Fig. 11. GPR profiles 1 and 2 plotted as a semi-transparent overlay on a coincident ERT resistivity models. The dashed line separates the peat (P) and clay layer (C) detected on GPR data. V. Pezdir et al. collect data below the clayey depressions, and we are not able to observe whether glacial deposits are present beneath them, we consider that forming a lake directly on carbonate bedrock would be mostly impossible. Therefore, till and glacial deposit material (underlain with carbonate bedrock) represent the basis for the lake sediments deposition. Presumably, lakes first formed sometime before the last glaciation. Retreat of the glacier caused till and glacial deposits that retained the water and caused multiple lake basins (Fig. 12a). After the last glaciation, lakes merged into one larger lake and new glacial deposits formed which acted as dams for the water (Fig. 12b). The older lake sediments were covered with younger sediments. In Holocene the peat started forming on the edges and shallower areas (Fig. 12c) and continued to grow towards the centres of depressions, following typical peat forming processes of terrestrialization. This is confirmed by the results (GPR, ERT and manual probing), as the depressions are filled with a thick layer of peat and clay mixture, while this layer is thinner or absent on the edges and ridges. Today, the area is completely covered with peat and the lake is no longer present (Fig. 12d). Some areas already show domed structure and drier areas (edges and ridges) are covered with shrub vegetation.

Conclusions
GPR and ERT methods proved to be a highly efficient complementary geophysical approach for delineating peat morphology and its underlying stratigraphy, which, together with peat probing, indicate the development of the studied bog. All methods depend on the weather conditions, most important being whether the ground is thawed or frozen. On the basis of the presented results, the chosen conditions proved to be appropriate. Temperature has an important role, as it determines the state of the ground, while precipitation influences the water content of the peat. GPR results provided detailed data about the depth of the peat deposit and delineated peatland surface morphology, which showed that the bog was not initially connected to a single entity; instead, four depressions existed in the form of smaller lakes, which probably merged into a single lake after the last glaciation (Fig. 12). The lakes were almost certainly formed by glacial deposit dams, as the carbonates do not have the capability to retain water on surface. The largest is depression 2 (Fig. 8) with a surface of at least 15,000 m 2 and a depth of 9 m, followed by depression 1 (15,300 m 2 , depth 8 m), 3 (4900 m 2 , depth 7 m) and 4 (3400 m 2 , depth 6 m), while peat depth on the edges Fig. 12. Theoretical schematic model of the peatland formation in a S-N direction (cross-section can be compared to GPR Profile 13). On the carbonate bedrock, till and glacial deposit material represent the basis for the lake sediments deposition. (a) In the first phase, smaller lakes form, separated with a glacial deposit material; (b) In the second phase, new glacial deposits formed after the last glaciation and smaller lakes merged into a larger lake. The older lake sediments were covered with younger sediments; (c) The peat started forming on the edges and shallower areas; (d) Present state of the peatland with forest at the edges and shrub vegetation at the edges and areas inside the bog. ranges from 2 to 4 m. Stronger heterogeneous reflections on the GPR images are generally observed in the deeper parts of the peat depressions which indicate the presence of mixed layers of peat, silt, sand and clay, supported also with ERT results (Fig. 11). Similar conclusions can be made on the basis of manual peat probing, as we often encountered strong resistance from sandy and clayey layers within the peat at greater depths, coinciding with stronger reflections within peat in GPR images. Future coring and core analysis will provide information about the sand and clay layers within the peat and thus explain the possible source of the radargram reflections. The larger amounts of sediments present within the earlier bog depressions can be explained by variations in climate conditions, where the entire area is occasionally exposed to flooding/fluvial sedimentation (Huang et al., 2011).
The peat layer is relatively well defined on ERT results, while the mixture of peat and clay at the base of the peat layer significantly reduced the resistivity values, which could lead to the misinterpretation of its thickness. This proves that the use of different geophysical methods together with direct observations is very important for a more correct interpretation. The distribution of the clayey sediments indicates the existence of deep, narrow lakes below peat depressions 1 and 2, reaching depths of more than 17 m on profile ERT-1 and more than 27 m on profile ERT-2, which were formed sometimes in Pleistocene. Later, however, in the late-Glacial and early-Holocene, the lake expanded across the entire area of today's bog, as shown on schematic model (Fig. 12).
Glacial deposits (MRB, Figs. 9 and 10) are found at an average depth of 9 m at ERT-1 and around 15 m at ERT-2. Considering the thicknesses of the clayey deposits above the glacial deposits and previously mentioned radiocarbon dating of the peat (Maja Andrič, personal communication, 2019), they can probably be attributed to older glacial deposits.
The peat depressions and underlying deeper clayey depressions could also be preconditioned by the carbonate bedrock, meaning that the depressions could have existed in the bedrock as well. Since Pokljuka is considered a karst plateau (Šifrer, 1952), the existence of the dolines in the carbonate bedrock is highly possible below the observed clayey and peat depressions. However, the karst phenomena in the form of dolines etc. were reshaped and covered with glacial, fluvioglacial and lake sediments during the Pleistocene and Holocene (Šifrer, 1952).
The integrated analysis of the geophysical results will serve as a basis for further geochemical studies, which will provide more details on the type and formation of the Š ijec bog. Geophysical methods and peat probing provide data on the morphology, thickness of the peat and stratigraphy of the basin. Information on the peat, the vegetation changes within the peat, sand and clay content, chemical parameters and composition, and their dependence on the peat basin morphology will be determined on the basis of further coring and research.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.