Hydrogeological Model of the Forefield Drainage System of Werenskioldbreen, Svalbard

Abstract

The significant recession of Arctic glaciers caused by climate warming is expanding their proglacial zones. Thus, their importance for the hydrology of glacierised basins is growing. In contrast to the surface waters in such areas, the role of groundwater in the hydrological balance of Svalbard catchments is poorly known. This paper presents the hydrogeological conditions and groundwater flow within the permafrost active layer in the forefield of the Werenskioldbreen glacier basin (44.1 km²), 61% of which is glacierised. Based on field studies of groundwater in the 2017 ablation season and laboratory analyses of the hydrogeological properties of proglacial sediments, a three-dimensional groundwater flow model (FEFLOW) for part of the glacier forefield (4.8 km²) was developed. The main results show the components and characteristics of the groundwater balance and indicate the preferential groundwater flow paths. The volume of water retained in the sediments of the marginal zone is 1.0073 mln m³. The maximum potential free pore space that could be filled by water is 2.0689 mln m³. The calculated groundwater discharge for average conditions is 6076.9 m³ d⁻¹, which is about 2% of the total seasonal catchment runoff from the main glacial river. The results of the spatial analysis for the groundwater depth and the groundwater flow directions are also presented. There need to be further detailed studies of hydrogeological processes in glacial basins in Svalbard in order to develop existing knowledge.

1. Introduction

Ongoing climate changes [1,2,3,4] are responsible for the increased melting and recession of Arctic glaciers [5,6,7,8,9,10,11] and the more intensive thawing of permafrost soils [12,13,14,15,16]. As a result of glacier recession, new areas of land are being uncovered, and glacier forefields are being enlarged [17,18]; changes are also taking place in the hydrological balance and in the flow regime [19,20,21], which may lead to significant changes in the drainage system of glaciated catchment areas.

A thicker active layer and a longer period of thawed sediments allow water to infiltrate more deeply into the ground [16,22,23]. Such infiltration and the formation of a saturation zone above the frost table can activate groundwater flow [24]. Any expansion of an unglaciated area, combined with the deepening of the active layer, generates a larger soil or rock volume potential for the infiltration, flow and storage of water in glacial sediments, which also influences the surface water level in the basin. In the context of the intensively warming Arctic climate, the current importance and potential further role of subsurface water flow and storage capability by expanding proglacial areas have not yet been considered for Svalbard [4,19,20]. Thus, our results may be a helpful example showing how modelling investigations might be used in such areas to estimate groundwater discharge, groundwater balance and groundwater paths in glacierised basins, now and in a warmer future.

In studies of glacierised Arctic basins, conspicuously little research has been devoted to the role of subsurface hydrology in an ice-free area, where the flow of groundwater forms during the short summer [25,26,27,28,29,30,31]. Therefore, we focus here on investigating the drainage system in the subsurface part of the proglacial zone. Particular emphasis is placed on identifying the hydrogeological conditions prevailing in an area with glacial sediments that enable groundwater flow and storage. The research took place in the well-studied Werenskioldbreen basin, which already has a multi-year series of meteorological, hydrological and glaciological data [32,33,34,35,36,37]. The first attempt at hydrogeological modelling was undertaken for the Werenskioldbreen basin to illustrate the subglacial drainage in the ground moraine layer and at the top of the bedrock [38]. The main emphasis in that work was on the glacierised part of the area but with a substantial generalisation of the hydrogeological conditions of the various sediments covering the forefield.

In the present research, we decided to explore the proglacial zone in much greater detail. Our aim was to gain an understanding of the interaction between surface water and groundwater in order to calculate the maximum possible capacity for water storage in glacial sediments and to discover the proportion of the subsurface part in the hydrological balance and catchment area drainage system of the proglacial part.

Groundwater flow models are valuable digital tools supporting our understanding of subsurface flow, enabling us to analyse and predict how an aquifer system will behave under various conditions. We applied groundwater flow modelling in order to understand the role of shallow subsurface processes that govern how groundwater flow is routed through the research area. The first step of this study involved identifying the hydrogeological properties of glacial sediments in order to gather the required input data. A three-dimensional groundwater flow model was developed in the FEFLOW tool for steady-state conditions based on data from the summer ablation season in 2017. The in situ observations of the hydraulic head in temporary observation wells were used for model calibration and validation.

There are several tools on the market that are used for hydrogeological or cryohydrogeological modelling research in cold regions [39]. The FEFLOW model was chosen for its 3D flexible-mesh groundwater flow computation capabilities and its wide scope for future expansion with new modules for heat and/or mass transfer. The model enabled us to calculate the groundwater balance for the whole basin (or a selected part of it) and to simulate the groundwater flow with spatial analysis of the results. This tool is the state-of-the-art model for representing such an environment [39]. FEFLOW was also already successfully tested by Piechota et al. [38], who applied the tool for the first time in Svalbard. No other similar studies have been reported to date.

2. Study Area

Werenskioldbreen is a polythermal glacier in recession situated in Wedel-Jarsberg Land, southwest Spitsbergen (Figure 1) [36,40,41]. Between 1957 and 2019, the glacier retreated by c. 1200 m (25 m/year) [42], and it covers about 61% of the catchment area (of 44.1 km²) [37]. The flat ground moraine studied here (area = 4.8 km²) is part of the forefield. The catchment area boundaries of the forefield are very clearly outlined by the glacier and the ice-cored terminal and lateral moraines.

The Werenskioldbreen basin lies within the Proterozoic Hecla Hoek Succession, which was formed from metamorphic rocks and other rocks of marine origin, with a total thickness (stratigraphic column) of 15–17 km [43]. The formations in this Succession are tectonically disturbed, exhibit numerous overthrusts and are crossed by faults formed during the Caledonian orogeny and modified during the Alpine orogeny [44].

The catchment’s mountain ranges are built of Precambrian metamorphic rocks belonging to the Eimfjellet group (southern part of the glacier), consisting of amphibolite, quartzite and chlorite schists, and also the Deilegga groups (eastern part of the glacier), consisting of phyllites with quartzite and silt intercalations, as well as calcareous and chlorite schists [45]. The north-western part of the basin belongs to the Sofiebogen group, which is composed of green schists, muscovite-carbonate-quartz and carbonate-chlorite-quartz schists [46,47]. The proglacial zone belongs to the Sofiebogen and Eimifjellet groups, represented by mica-carbonate-quartz and grey calcite marbles [46].

Accordingly, the largest part of the proglacial zone is ground moraine (Figure 2) covered by a boulder lag composed of all rock types from alimentation areas or fine-grained outwash deposits of proglacial rivers [44,49]. A significant part of the proglacial zone consists of sandurs and outwash plains. At the glacier front, moraines up to several metres high are formed from recession.

The surface sediments are composed mainly of poorly sorted sand and gravel and moraine till [50]. These forms are continuously being transformed owing to the highly dynamic thermo-erosive degradation of buried ice and co-occurring fluvioglacial processes.

The surroundings of the study area consist of mountainous terrain with sharp peaks (Figure 1), widespread superficial debris cover and talus slopes based on lateral moraines [44]. Hillslope processes, such as debris flows, rockfall and slopewash, occur on the steep mountain and moraine slopes [51]. The ground is unstable and covered with rock cobbles, making it difficult for soil to form or for vegetation to expand. However, in the more stable zones in the southern lateral moraines, in their flat central plateau and in the foothills, soil formation is taking place, specifically of Cryosols (Leptic Cryosol, Cambic Glacic Cryosols and Leptic Gleysols) [52]. In addition, a dry lichen-moss tundra is developing here. The forefield area is morphologically active, dominated by glacial and fluvioglacial landforms and modified by melt-out processes [44]. The well-developed river network and seasonally occurring flooded areas make surface water flow and sediment transport possible. Therefore, in the northern and central parts of the forefield, within the sandur fans and plains and the sandur with ramparts, such vegetation that exists occurs locally in small patches. The evolution of tundra is perceptible on the extensive ground moraine and around lakes. A unique habitat for flora and fauna surrounds the proglacial Mewie Lake (Figure 2), where soils have developed from stratified limnoglacial sediments [52]. Here, the better-rooted lichen-moss tundra is expanding faster, and gull flocks fertilise it in a natural way. Larger animals rarely appear here; it is only a transit area for reindeer or foxes (occasionally polar bears). This terrain is generally unsuitable for both cultivation and habitation.

Because the root systems of the local and polygenic vegetation are poorly formed, the modelling study did not take them into consideration. The proglacial part of the basin is supplied by rainfall and the glacier’s meltwater. In the frontal zone, the glacier is directly frozen to the ground, resulting in the formation of high-pressure outflows. The ablation water flows out by point exsurgent and pressurised springs or Röthlisberger (R) subglacial outflows at the glacier front and recharges the proglacial zone, giving rise to a number of surface watercourses and lakes [32,53,54]. Outflows of pressurised water also occur in winter, when so-called naleds are formed [55,56]. The main subglacial outflow is situated at the Kvisla gate (Figure 2), in the northern part of the glacier; it drains 60–75% of the proglacial waters [38]. The remaining subglacial outflows include the outflow on the southern side of the central moraine and the southernmost outflow draining the Angellkroken area.

To the west of the central moraine, there is an area with a large accumulation of buried glacial ice, the melting of which contributes to the degradation of the inner outwash plain, where artesian outflows of glacial water occur. The outflow from the Kvisla gate and several other streams form the Breelva River, which drains the whole basin in a south-westerly direction. The most southerly outflows flow into Lake Mewie (Makevatnet) and, when combined with the Breelva, the waters break the end moraine at only one gauging point, approximately 2.5 km from the glacier front, where a permanent hydrometric station is located [34]. Another important source of water supply to the forefield is runoff from the supraglacial drainage system, which supplies the watercourses flowing at the glacier front or the ephemeral pools forming there. This region is characterised by continuous permafrost [57,58]. In the non-glaciated area, permafrost commonly occurs at an average depth from 1 m to 2.5 m [16,22,31]. However, places where the permafrost is deeper may appear locally and are usually related to areas where heat transfer is more intensive [51], for example, in areas where water bodies occur or even in river channels where talik can occur [59,60]. Unfortunately, the accurate identification of such locations would require highly detailed, spatially condensed ground temperature investigations.

The permafrost stops the infiltration and movement of groundwater. The subsurface layer, where water movement is possible, is the active layer thickness (ALT), which thaws during the summer season. Going into greater detail, we can also point to the so-called transient layer between the permafrost and the active layer, which can freeze and thaw on a long-term scale of decades or centuries [51]. The presence of permafrost and the active layer has a significant impact on geomorphological, hydrological and hydrogeological processes [24,51].

In spring, when the snow cover disappears, allowing heat to transfer to the ground, the ground layer in the Werenkioldbreen basin thaws from above [51], activating the water infiltration process. The aquifer is formed in this period (generally July–September), and groundwater flow is possible. In autumn, when the air temperature drops below 0 °C, the ground temperature decreases from above. However, the ALT is also affected by the permafrost temperature (<0 °C) from below and is reduced [51]. Water movement eventually stops (generally in October). The water that is present changes its physical state and is stored until the next thawing season.

In summer, the possibility of water infiltration in the catchment under analysis here changes owing to the thawing of the active layer, and consequently, the potential volume of the aquifer increases. In the Werenskioldbreen basin, the conditions (ALT) become more stable after the end of July.

The average annual air temperature in 1979–2016 at the Polish Polar Station (PPS) in Hornsund (16 km away from the Werenskioldbreen basin) was −3.6 °C [61,62]. In 2017, the mean annual temperature was −1.6 °C [62], and for the summer period (June–September), the daily average air temperature ranged from 1.7 °C to 7.5 °C (Figure 3). The total precipitation in June–September was 226.6 mm [63], and the highest amount of precipitation fell in late September (max. 73 mm per day, 18 September). There is a clear relationship between the groundwater depth and river discharge.

3. Materials and Methods

3.1. Sedimentological Analyses

In order to identify the various glacial sediments and their hydrogeological properties, 44 sediment samples were analysed. They were collected at 18 sites from depths of 0.30 to 1.5 m by manual drilling and digging within the forefield, as well as in the locations of the naturally exposed soil layers in river channels (Figure 2). The sediment samples were analysed at the Soil and Sediment Analysis Laboratory of the University of Silesia in Katowice. The purpose of the study was to identify and classify the sediment types based on the PN-EN ISO 14688-2:2006 norm [64]. Next, the hydraulic conductivity was calculated from the grain-size distribution [65] using the MULTI H_K program written in Visual Basic for Applications (made publicly available by the authors) [65]. By means of a special macro in MS Excel, this makes it possible to calculate the hydraulic conductivity using, according to the stated assumptions, a compilation of 20 different literature formulas. The input data are the percentages obtained from sieve analysis and sample porosity; the latter was calculated using the method described in Rosas et al. (2014) [65].

3.2. Other Model Input Data

3.2.1. DEM and Maps

The Digital Elevation Model and orthophotomap were based on 117 aerial photos taken by a Mavic PRO Unmanned Aerial Vehicle (UAV) fitted with an FC220 camera (focal length: 35 mm; charge-coupled device: 4000 × 3000 pixels, DJI, Shenzhen, China) on 15 September 2017. The final product was an orthophotomap with a 0.15 × 0.15 m raster field resolution and a DEM with a 1 × 1 m raster field resolution, corrected by a 3GCP (ground control point) measured during the aerial survey with accurate GNSS receivers (LEICA 1230). These products were mapped onto the UTM for the zone 33N and on the ellipsoid WGS 84, with an accuracy of x = 1.52 m, y = 0.95 m and z = 0.65 m. The results were used to prepare the model structure and maps in this paper.

In addition, electronic maps from the Norwegian Polar Institute Map Data and Services and Polish Polar Database and prepared by L. Kolondra [66] and Karczewski et al. [49] were used when drafting the figures in order to verify the products and for the spatial analysis of the research area.

3.2.2. Groundwater Level and Rainfall

Groundwater levels were measured in a temporary monitoring network consisting of six piezometers. In July–August 2017, all of the piezometers were measured every 2–3 days. These measurements provided the data used for calibrating the groundwater flow model. The rainfall data were taken from the WMO (World Meteorological Organisation) meteorological station No. 01003 at the Polish Polar Station (PPS) in Hornsund.

3.3. Model Assumption

In this study, the FEFLOW (Finite Element subsurface FLOW system) modelling package, version 7.0 [67], was used as the engine for simulating groundwater flow.

The mathematical formula for groundwater flow in a porous medium is based on two fundamental principles: Darcy’s Law and the principle of continuity. The general flow equation representing three-dimensional (3D) transient groundwater flow for unconfined, heterogeneous and anisotropic conditions is given by Equation (1):

$$K_x\frac{{\partial }^{2}h}{\partial x^2}+K_y\frac{{\partial }^{2}h}{\partial y^2}+K_z\frac{{\partial }^{2}h}{\partial z^2}=Sy\frac{\partial h}{\partial t}-R$$

(1)

where K_x, K_y and K_z are the respective principal components of the hydraulic conductivity along the x, y and z axes (LT-1), h is the hydraulic head (L), S_y is the specific yield (dimensionless), R is a volumetric sink/source term, defined as the volume of water injected per unit time and unit volume of the porous medium (T − 1), and t is time (T). As the flow model of the Werenskioldbreen forefield represents the steady-state conditions in the 2017 ablation season, Equation (2) can be simplified ($\frac{\partial h}{\partial t}=0$):

$$K_x\frac{{\partial }^{2}h}{\partial x^2}+K_y\frac{{\partial }^{2}h}{\partial y^2}+K_z\frac{{\partial }^{2}h}{\partial z^2}+R=0$$

(2)

4. Groundwater Flow Model Parameterisation

The steady-state FEFLOW model was developed for average parameters reflecting the conditions present in the 2017 ablation season. To assess the contribution of the groundwater component in the total runoff from the forefield area, we present a model for the condition when the ALT potentially has its maximum thickness. The monitoring data were collected to provide input (averaged parameters) and validate the modelling results. The model extent includes the glacier marginal zone, which was determined by drawing the elevation boundary for 25–35 m a.s.l. under the steep mountain slopes and terminal and lateral moraines. The area where surface water runoff dominates and groundwater storage is impossible was omitted. In addition, the glaciated area was rejected in this study, so the model boundary was delineated along the glacier front. The total model area was 4.8 km².

4.1. Geometric Model and Finite Element Mesh

To reflect the hydrogeological model framework, relevant input data such as the DEM, the orthophotomap, and the top and bottom of the layers were prepared in ArcGIS software. Finally, three model layers were created with thicknesses of 2.0, 0.3 and 0.2 m. In addition, five geomorphological polygons were defined for model generalisation; for each one, the three sediment layers were diversified using appropriate hydrogeological parameters based on the results of the laboratory analysis. The model was simulated as an unconfined aquifer. The model assumes the possibility of groundwater flow for a maximum thickness of the active layer up to 2.5 m, which lies within the ranges given in the literature for the Hornsund area [16,23,31,68].

4.2. Boundary Conditions

The boundary conditions (BCs) define the system conditions and fluxes into and out of the model domain [39,67]. They are used at the aquifer model’s borders as well as within it, in particular, to represent the relationship between surface and groundwater. Because the catchment in the study area is bounded by ice-cored moraines and a glacier, it was assumed that no exchange of water at these external boundaries of the model was possible (fluid flux boundary conditions = 0). The only place where water can flow out beyond the model’s external boundaries is the gorge of the Breelva River, in the extreme south-western part of the model. Except for this place, the external boundaries of the model constitute an impermeable barrier. The boundary conditions inside the model (surface waters) were located on the basis of an analysis of the orthophotomap and DEM. The stream lines, lake boundaries and polygons were delineated in ArcMap software and imported as *.shp files to FEFLOW. For the Breelva, the main watercourse draining the glacier forefield, the 1st BC (Dirchlet-type, hydraulic-head BC) was set up for each of the model’s three layers. For the other surface waters with limited possibilities of exchange with groundwater, a 3rd BC (Cauchy-type fluid-transfer BC) was established for the 1st and 2nd layers.

4.3. Recharge (Water Input)

In the summer period, when the active layer within the research area has begun to thaw, “rainwater infiltration is influenced by the type of soil, the amount of fissures in the soil or bedrock, rainfall intensity and antecedent soil moisture content” [24].

In the model, the subsurface was recharged mainly by precipitation and, to a smaller extent, by surface water, particularly in a narrow zone along the glacier front (an area c. 100 m wide). Here, a higher level of infiltration was estimated to simulate meltwater inflow, which occurs at the glacier front in contact with the ground. The initial assumed values for subsurface recharge were based on rainfall data (557.2 mm rainfall/year in 2017) and soil permeability data [69,70]:

• Glacial sands and gravels: 110 mm/year (19.7% rainfall);
• Glacial sands and gravels along the glacier front: 420 mm/year (75% rainfall);
• Outcrops: 90 mm/year (16% rainfall).

4.4. Hydraulic Conductivity

Based on the field survey, laboratory results and literature data [52,71], it was found that the foreland material was represented by significant amounts of different types of sediments. The ability of water to flow through sediments was reconstructed in the model using hydraulic conductivity. This parameter was calibrated in the model and lies within the ranges of the variability shown in Table 1 for unconsolidated sediments, the hydraulic conductivity of which was determined in the laboratory. Within the boundaries of the model, no samples were taken for examination in the region where strongly cracked substrate rocks outcrop. The hydraulic conductivity of these rocks was assumed on the basis of literature data [72,73,74].

4.5. Calibration and Validation

The calibration tested a number of sets of model parameters to find the best fits for field measurements. Therefore, this process is crucial, indicates and provides information, and helps understand the modelled system’s behaviour and responsiveness.

When analysing changes in the distributions of piezometric pressure, we changed the input hydraulic conductivity obtained from calibration for each layer of the model several times. The calibration was based on the results of field surveys for particular types of soil (Table 1).

Hydraulic head parameters were calibrated on the basis of real ranges from monitoring network measurements. The depth to the water table in piezometers varied from 0.05 m to 0.15 m in the model. The model was also validated against data that were not used directly in developing the model, owing to the insufficient precision of the measurements or an insufficient number of such measurements. This type of validation provided some information regarding the possible ranges of variables. The data were used to both calibrate and validate the model against the conceptual model. An example of such data is the river flow measurements taken by the author with a hydrometric current meter, as well as earlier archived data. This enabled us to determine the sections of rivers where infiltration to the water-bearing horizon might occur.

In addition, a detailed spatial analysis of the correctness of the model results was undertaken based on the author’s two-month-long direct observations of the drainage system’s behaviour.

The results of the hydraulic head distribution in model simulation correlate very well with the field monitoring data. The measured and modelled depth to groundwater in piezometers (groundwater level) showed that our calculations and simulations were correct (Section 5.3). The statistical measures obtained were satisfactory: mean error—0.07 (m) (groundwater deviation); standard deviation (σ)—0.096; and mean square error—0.088 (RMS).

The water mass balance data also show that the model was properly designed and calibrated. Differences between the water input and output may appear owing to the numerical approximations of the groundwater equation, but the error should be no greater than one percent [75]. Larger errors would indicate an inconsistency or an error in the model. For the resulting model, this error is 0.007% (Section 5.5.2, rate budget imbalance).

5. Results and Discussion

5.1. Sediment Properties and Hydraulic Conductivity

The granulometric composition of 44 sediment samples was obtained from a laboratory analysis, from which the hydraulic conductivities were calculated, and nine principal types of sediments present in the study area were identified. The largest percentage of these sediments consisted of highly permeable gravelly sand (31.8%), sandy gravel (22.7%) and silty sand (18.1%) (Table 1).

The laboratory analyses showed that the post-glacial sediments (Figure 1), especially the sandur fans and outwash plains, as well as sandurs with ramparts and cones, consisted of layers containing highly permeable mixed sands and gravels with a variety of grain sizes. The samples representing moraine sediments are sandy clayey silt with overlying sandy gravel, with a low to very low permeability.

Table 1. The hydraulic conductivity calculated for the 9 sediment types identified, representing all samples collected in the study area.

ID	Sediment Type (or Bedrock)	Percentage of All Samples (%)	Results for Hydraulic Conductivity (k) by Methodology Rosas et al. (2014) [65]
			Average (m s−1)	Max (m s−1)	Min (m s−1)
1	Gravelly sand	31.8	3.60 × 10−4	6.90 × 10−4	3.00 × 10−5
2	Sandy gravel	22.7	1.05 × 10−3	3.62 × 10−3	5.30 × 10−4
3	Silty sand	18.1	1.65 × 10−5	3.01 × 10−5	1.03 × 10−5
4	Gravelly silty sand	6.8	2.12 × 10−5	2.71 × 10−5	1.29 × 10−5
5	Silty sandy gravel	4.6	5.30 × 10−4	8.20 × 10−4	2.50 × 10−4
6	Sand	4.6	4.00 × 10−5	6.00 × 10−5	2.00 × 10−5
7	Sandy clayey silt	6.8	1.64 × 10−8	4.11 × 10−8	1.70 × 10−9
8	Sandy gravelly silt	2.3	1.28 × 10−5	2.88 × 10−5	4.37 × 10−6
9	Silt	2.3	x	1.03 × 10−5	x

Table 1. The hydraulic conductivity calculated for the 9 sediment types identified, representing all samples collected in the study area.

ID	Sediment Type (or Bedrock)	Percentage of All Samples (%)	Results for Hydraulic Conductivity (k) by Methodology Rosas et al. (2014) [65]
			Average (m s−1)	Max (m s−1)	Min (m s−1)
1	Gravelly sand	31.8	3.60 × 10−4	6.90 × 10−4	3.00 × 10−5
2	Sandy gravel	22.7	1.05 × 10−3	3.62 × 10−3	5.30 × 10−4
3	Silty sand	18.1	1.65 × 10−5	3.01 × 10−5	1.03 × 10−5
4	Gravelly silty sand	6.8	2.12 × 10−5	2.71 × 10−5	1.29 × 10−5
5	Silty sandy gravel	4.6	5.30 × 10−4	8.20 × 10−4	2.50 × 10−4
6	Sand	4.6	4.00 × 10−5	6.00 × 10−5	2.00 × 10−5
7	Sandy clayey silt	6.8	1.64 × 10−8	4.11 × 10−8	1.70 × 10−9
8	Sandy gravelly silt	2.3	1.28 × 10−5	2.88 × 10−5	4.37 × 10−6
9	Silt	2.3	x	1.03 × 10−5	x

The hydraulic conductivity of the samples was sometimes lower than the values found in the literature [74,76]. This is most probably due to their poor sorting, as they often contain mixed finer fractions, which are the remnants of till outwash, for example, the washed-out parts of the middle moraine. This is a characteristic of glacial deposits in many different parts of the world [77,78].

Model calibration yielded the final values of the hydraulic conductivity: 3 × 10⁻⁵–1.05 × 10⁻³ m s⁻¹ for gravelly sand, 1.64 × 10⁻⁸ m s⁻¹ for moraine sandy clayey silt and 8 × 10⁻⁶–7 × 10⁻⁶ m s⁻¹ for bedrock [71,72]. This variability is due to the considerable lithological diversity of the sediments in the Werenskioldbreen marginal zone. Table 2 lists these results in detail. When comparing the results in Table 1 and Table 2, it is well to remember that the differences in the hydraulic conductivity are due, among other things, to the scale effect. The model results represent more highly aggregated hydraulic conductivity, which is an averaged value for a part of the whole of a geomorphological form or lithological formation.

The calculated hydraulic conductivities fall within the range of values determined for Scandinavian tills (sandy-silty tills), i.e., from 5 × 10⁻⁹ m s⁻¹ to 5 × 10⁻⁴ m s⁻¹ [77]. The hydraulic conductivity also accords with the parameters estimated using the PARAMEX method in piezometers deployed in the Ebba River basin (Billefjorden, Spitsbergen), where the range of parameters obtained was from 8.88 × 10⁻⁵ m s⁻¹ to 2 × 10⁻³ m s⁻¹ [28].

Table 2. Structure of the model layers and their final hydraulic conductivity. The five main sediment types represent the generalised result for the model layers within each type of geomorphological polygon.

Polygon Symbol	Type of Polygon Layer (Geomorphology)	Number of Layers: 1—2 m 2—0.30 m 3—0.20 m	Sediment	Final Hydraulic Conductivity k (m s−1) in Model
P1	Rock outcrops	1	Rock outcrops	8 × 10−6
		2
		3
P2	Ground moraine and bedrock	1	Gravelly sand	3.6 × 10−4 (or 4.6 × 10−4)
		2	Weathered crystalline	8 × 10−6
		3		7 × 10−6
P3	Sandur with ramparts and cones on dead ice	1	Sandy gravel	1.05 × 10−3
		2	Sandy gravel or gravelly sand	1.05 × 10−3 (or 3.6 × 10−4 and 4.6 × 10−4)
		3	Sandy gravel	1.05 × 10−3
P4	Parts of terminal moraine	1	Gravelly sand (or parts of silty sand)	3.6 × 10−4 (or 3 × 10−5)
		2	Sandy clayey silt	1.64 × 10−8
		3
P5	Sandur fans and plains (areas flooded by rivers)	1	Gravelly sand or sandy gravel (locally)	6.9 × 10−4 (or 1.6 × 10−3 locally)
		2
		3	Weathered crystalline rock	7 × 10−6

Table 2. Structure of the model layers and their final hydraulic conductivity. The five main sediment types represent the generalised result for the model layers within each type of geomorphological polygon.

Polygon Symbol	Type of Polygon Layer (Geomorphology)	Number of Layers: 1—2 m 2—0.30 m 3—0.20 m	Sediment	Final Hydraulic Conductivity k (m s−1) in Model
P1	Rock outcrops	1	Rock outcrops	8 × 10−6
		2
		3
P2	Ground moraine and bedrock	1	Gravelly sand	3.6 × 10−4 (or 4.6 × 10−4)
		2	Weathered crystalline	8 × 10−6
		3		7 × 10−6
P3	Sandur with ramparts and cones on dead ice	1	Sandy gravel	1.05 × 10−3
		2	Sandy gravel or gravelly sand	1.05 × 10−3 (or 3.6 × 10−4 and 4.6 × 10−4)
		3	Sandy gravel	1.05 × 10−3
P4	Parts of terminal moraine	1	Gravelly sand (or parts of silty sand)	3.6 × 10−4 (or 3 × 10−5)
		2	Sandy clayey silt	1.64 × 10−8
		3
P5	Sandur fans and plains (areas flooded by rivers)	1	Gravelly sand or sandy gravel (locally)	6.9 × 10−4 (or 1.6 × 10−3 locally)
		2
		3	Weathered crystalline rock	7 × 10−6

5.2. Model Structure

Detailed examination of the granulometric composition of the sediments from the study area enabled the ranges of the parameters representing them in the model (hydraulic conductivity) to be narrowed down. As a result, there was no need to test broad ranges of literature data during the calibration. The final structure of the model layers diversified by the hydraulic conductivity is presented above in Table 2.

The final maximum thickness of the model layers did not exceed the maximum thickness of the active layer in this area (2.5 m). A long series of data relating to ground thaw monitoring around the Polish Polar Station in Hornsund is available [15,16,79,80,81]. The average active layer thickness (ALT) in Hornsund in the coastal zone in 1981–2000 was 1.75 m (max 2.2 m) [82]. New data (2014–2017) for the Fuglebekken basin (Hornsund) [31,83] show that the maximum average depth for ALT was 2.4 m. In recent years, other locations near the Werenskioldbreen basin have also been studied. The results from a few studies conducted within the coastal zone and the neighbouring Brateggdalen valley indicate that the active layer thickness did not exceed 2.5 m, while for sorted coastal material (10 m a.s.l.), it was 1.0–1.5 m, and for unsorted tundra material (10 m a.s.l.), it was 2.0–2.5 m [23,68]. With respect to the rate of soil thaw in the Hornsund area, our model assumed the most relevant maximum depth of thawing to be 2.5 m, as a rule of thumb for the model’s construction, in order to determine the groundwater retention potential for the applied model domain.

5.3. Hydraulic Head (Observed and Modelled)

Piezometric monitoring of the groundwater showed the amplitude of their levels in July–August 2017. Deployed in the marginal zone of the glacier, the piezometers recorded the mean groundwater level as 0.13–0.37 m in depth. The average amplitude of the groundwater level was 0.15 m (range 0.08–0.33 m). The correctness of the model’s performance was verified on the basis of the above monitoring data (

Table 3. Groundwater level measured against modelled values.

Field Data (Monitoring)						Modelled Data
			Depth to Groundwater
ID	Elevation (m a.s.l.)	(A) Groundwater Level Average (m a.s.l.)	Average (m)	Max (m)	Min (m)	(B) Hydraulic Head (m a.s.l.)	Depth to Groundwater (m)	Groundwater Deviation (B − A) (m)
P1	14.71	14.43	0.28	0.22	0.30	14.52	0.20	0.08
P3	14.92	14.79	0.13	0.00	0.33	14.72	0.20	−0.07
P4	14.71	14.34	0.37	0.22	0.44	14.50	0.22	0.16
P5	15.28	14.90	0.37	0.29	0.44	14.90	0.38	−0.01
P7	28.38	28.07	0.31	0.24	0.33	27.98	0.40	−0.09
P8	26.64	26.34	0.30	0.24	0.31	26.32	0.33	−0.02

). Comparison of the mean groundwater level measured in the field with that calculated using the model showed that the latter had been properly and reliably calibrated. The mean error between the modelled water level and that measured in the field was 0.07 m.

Modelling the hydraulic head yielded the distribution of groundwater levels for the whole of the glacier’s forefield (Figure 4). Values were the highest in the eastern part of the model (36–48 m a.s.l.), diminishing towards the west. The lowest values were obtained in the Breelva River gorge at an altitude of 13 m a.s.l.

The modelled groundwater level largely tallied with the area’s morphology, which is characteristic of groundwater occurring just below the ground surface. In areas no more than 20 m a.s.l., the groundwater level varied from a dozen or so centimetres to 0.5 m below the surface, while at altitudes no higher than 35 m a.s.l., the level was at 0.5–1 m below the surface.

Modelling highlighted the areas where the groundwater level exceeded the ordinate of altitude, which accords with what was observed during the fieldwork. Pools of water, no more than 1 m deep, occur in these areas. The modelling results are reflected in the terrain, where they coincide with periodically observed flooded areas and lakes (Figure 5).

5.4. Directions of Groundwater Flow

The geomorphological and geological forms of the Werenskioldbreen forefield, varying in space and with depth, are largely responsible for the different directions of groundwater flow and the distribution of drainage zones. The intensity of water flow is the greatest in the first and second layers of the model down to a depth of 2.30 m, which, in these layers, is due to the predominance of formations with good or very good permeability.

Using streamlines in FEFLOW to show trajectories of particles following the velocity vectors of the flow based on a steady-state flow field, it was possible to analyse the directions and velocities of groundwater flow, from which preferential flow paths could be determined, as well as areas of less intense flow, associated with the zone of poorly permeable sediments in the central part (Figure 5). The modelling results indicate that a watershed is forming in the central part, separating two subsystems of shallow groundwater circulation, i.e., a northern and a southern one.

In general, groundwater flows westwards from the northern and north-eastern zones towards the old, no longer active Breelva outlet point within the northern margin of the end moraine. Until the 1970s, when large volumes of water destroyed the remains of the medial moraine, water from the northern part of the basin used to drain westwards, and water from the southern part drained in a south-westerly direction [53,84]. However, nowadays, the water from this part of the basin ultimately enters the Breelva channel and flows in a south-westerly direction to the only outlet point, within the southern margin of the end moraine (Figure 5). In the northern part of the catchment, the distinctive sandurs and outwash plains constitute extensive areas with preferential groundwater flow paths.

Generalisation of the model in the northwestern part indicates that the aquifer is supplied with water and drained solely by the Breelva River. In fact, however, there is also a network consisting of a great many narrow, shallow channels fed by ablation waters, which form ephemeral pools. The water from these pools ultimately finishes up in the main channel of the Breelva, albeit farther down its course. We did not build an extensive network of either the watercourses or the lakes in this area, because the former were continually changing their courses, and even the diurnal changes in water states associated with the supplies of ablation waters were too great here. As a result of the schematisation of the real conditions, we focussed on the constant, stable hydrographical components of the catchment area, in other words, the main channel of the Breelva. The model results confirmed the localised appearance on the surface of waters from the active layer of the substrate (Figure 5); this was observed from time to time during fieldwork.

The northern part of the catchment area is fed mainly by water from precipitation, the Breelva and River 3 (Figure 5), an artesian outflow that is part of the subglacial drainage regime. On the other hand, water flows out of the model to the Breelva, which makes a very considerable (29% of the total model outflow) contribution to the drainage in this part of the model.

The groundwater in the southern part is supplied mainly by precipitation and by a watercourse formed by the inflow of ablation waters near the slopes of Angellfjellet (Figure 1) (River 1: 13% inflow, 4% outflow); there is also an exchange of waters with respect to flow-through lakes. The largest of these, Lake Mewie (Makevatnet), largely drains the water-bearing layer (15% of the total outflow from the model). In addition, it should be noted that River 4 flows out of this lake, finally draining this part of the basin into the Breelva.

The modelling results show that there is considerable differentiation in the exchange of surface and groundwater. Most surface watercourses are diverse in character and, along certain stretches, feed or drain the aquifer; only the Breelva River predominantly drains, which confirms its role as the main centre of groundwater drainage.

5.5. Groundwater Balance

5.5.1. Aquifer (Active Layer) Storage Capacity

The domain volume is 1.1973 × 10⁷ m³, and the volume of the pore space contained in the model, with possible infiltration (void volume), is 2.0689 × 10⁶ m³. This is the maximum volume of water that can be retained by the sediments in the studied part of the forefield, assuming an active layer thickness of 2.5 m.

The modelled average conditions for which the model was calibrated indicate that the fluid volume in the modelled sediment layers is 1.0073 × 10⁶ m³. This is the volume of water present between the groundwater surface and the bottom of the model.

We also performed a prognostic simulation, assuming a 0.5 m increase in the thickness of the deepest layer in the model, while retaining all of its previous parameters, in particular, the lithological differentiation of the sediments. This enhanced the groundwater storage potential in the sediments, which, given ongoing climate warming, appears to be a more and more likely scenario, i.e., that the active layer will become deeper in the future. As a result, the domain volume of the model domain increased to 1.4368 × 10⁷ m³, and the void volume was 2.4986 × 10⁶ m³. In other words, the potential water retention in the modelled area with a 3 m thick active layer increased by 20.8%.

5.5.2. Groundwater Balance

The dynamic resources of the active layer for the model are 6077 m³ d⁻¹ (Table 4); 18.2% of the groundwater recharge comes from the Breelva River, and 55.3% comes from other watercourses and surface water bodies, while 26.5% comprises infiltrating water from precipitation and surface melting along the glacier front. The substantial dominance of water infiltrating from surface watercourses over recharge due to precipitation is because the surface water state was mapped onto the model from the glacier’s ablation period. High surface water levels cause intensive infiltration of these waters into the aquifer.

The highest groundwater recharge in the model (Figure 4), of the order of 420 mm/year, resulting from the natural infiltration of water from precipitation or ablation, occurs in a 50–200 m wide zone along the glacier front (5% of the model’s area), where there is intensive surface runoff from the glacier and where pools sometimes form. The least recharge—90 mm/year—occurs where precipitation water infiltrates from rocky outcrops, especially by the slopes of Angellfjellet (Figure 1) and in the northern part of the catchment area at the foot of the slopes of Wernerknatten (7% of the model’s area). The recharge of the remaining surface of the model (88%) from glacial sediments is uniform: 110 mm/year, of which 19.7% is rainfall.

Recharge of the aquifer from surface waters is the greatest in the northern and central parts of the catchment area (Table 4), which are covered by sandurs consisting of highly permeable sediments: here, 18.6% of the water in the aquifer comes from River 6 (artesian spring), and 9% comes from the Breelva. In turn, in the southern part of the catchment area, within the ground moraine, where there are numerous water bodies and watercourses flowing through them, 13% of the aquifer’s water infiltrates from River 1, c. 5.6% arrives from Lake 1, and 2.4% comes from Lake 2, into which River 1 flows.

Water flows out of the model into watercourses and surface water bodies. The principal watercourse here is the Breelva, which removes almost 40% of the groundwater from the model; the other watercourses and water bodies drain 60% of the groundwater from the model.

The recharge or drainage in such a shallow groundwater circulation system is highly variable. The dominant drainage stretch in the main river, the Breelva, is in its upper course, as far as the groundwater divide. The northern part of the Breelva’s catchment drains 29% of the aquifer but supplies 9% of its water, whereas the southern part drains 10% of the aquifer but supplies 9% of its water. Mewie Lake plays an important part in draining groundwater from the active layer in the southern part of the modelled area, removing 15% of its water. Overall, the water bodies in the northern part of the catchment are responsible for draining 43% of the water from the aquifer, and in the southern part, they are responsible for draining 57%.

Table 4. Rate budget in FEFLOW. The locations of each object, i.e., rivers (R) and lakes (L), are marked in Figure 5. The border between the northern and southern parts of the catchment areas for summary calculations is the groundwater divide line, as shown in Figure 5.

Analysed Element	Inflow		Outflow
	(m3 d−1)	(%)	(m3 d−1)	(%)
The Breelva river (1st BC)	1107.6	18.2	2401.5	39.5
Breelva—North part	544.4	8.9	1758.2	28.9
Breelva—South part	563.2	9.3	643.3	10.6
Other streams and lakes (3rd BC)	3360.2	55.3	3675.4	60.5
River R1 (S)	774.0	12.7	263.8	4.3
River R2 (S)	221.2	3.6	139.6	2.3
River R3 (S, connection with Mewie)	187.3	3.1	263.8	4.3
River R4 (S)	330.0	5.4	283.7	4.7
River R5 (S)	225.8	3.7	336.4	5.5
River R6 (N, artesian spring)	1132.7	18.6	794.9	13.1
Lake L1 (S)	343.2	5.6	427.6	7.0
Lake L2 (S)	143.9	2.4	191.7	3.2
Mewie Lake (S)	2.0	0.0	918.3	15.1
Lake L3 (N, moraine)	0.0	0.0	55.8	0.9
Source—effective rainfall	1609.5	26.5	x	x
Summary North part of the basin	2619.6	43.1	2611.1	43.0
Summary South part of the basin	3457.7	56.9	3465.8	57.0
Summary (total)	6077.3	100.0	6076.9	100.0
Imbalance (Inflow–Outflow)	0.4	0.007	x	0.0

Table 4. Rate budget in FEFLOW. The locations of each object, i.e., rivers (R) and lakes (L), are marked in Figure 5. The border between the northern and southern parts of the catchment areas for summary calculations is the groundwater divide line, as shown in Figure 5.

Analysed Element	Inflow		Outflow
	(m3 d−1)	(%)	(m3 d−1)	(%)
The Breelva river (1st BC)	1107.6	18.2	2401.5	39.5
Breelva—North part	544.4	8.9	1758.2	28.9
Breelva—South part	563.2	9.3	643.3	10.6
Other streams and lakes (3rd BC)	3360.2	55.3	3675.4	60.5
River R1 (S)	774.0	12.7	263.8	4.3
River R2 (S)	221.2	3.6	139.6	2.3
River R3 (S, connection with Mewie)	187.3	3.1	263.8	4.3
River R4 (S)	330.0	5.4	283.7	4.7
River R5 (S)	225.8	3.7	336.4	5.5
River R6 (N, artesian spring)	1132.7	18.6	794.9	13.1
Lake L1 (S)	343.2	5.6	427.6	7.0
Lake L2 (S)	143.9	2.4	191.7	3.2
Mewie Lake (S)	2.0	0.0	918.3	15.1
Lake L3 (N, moraine)	0.0	0.0	55.8	0.9
Source—effective rainfall	1609.5	26.5	x	x
Summary North part of the basin	2619.6	43.1	2611.1	43.0
Summary South part of the basin	3457.7	56.9	3465.8	57.0
Summary (total)	6077.3	100.0	6076.9	100.0
Imbalance (Inflow–Outflow)	0.4	0.007	x	0.0

On the day when the field measurements were made (15 September 2017) to establish the state of surface waters, the Breelva was draining c. 0.31 million m³ d⁻¹ of water from the forefield (Figure 3), the daily mean was 3.6 m³ h⁻¹, and the mean water level was 0.18 m. The modelled groundwater inflow (outflow) was 6077 m³ d⁻¹, which makes up 1.96% of the drainage via the Breelva. The amount of water drained from the catchment area on that day was relatively small, as was the case on the day before (14 September 2017) and on the day after (16 September 2017). During these three days, the mean flow rate was 4.2 m³ h⁻¹, and the mean water level was 0.21 m, whereas for the entire period during which the Breelva was monitored (5 July–27 September, 85 days), the mean volume of flow was 9.3 m³ h⁻¹, and the mean level was 0.40 m. These data also confirm the lack of any significant ablation water supply during this time, which would have been recorded as a peak in the discharge. In addition, there was no rainfall on 14–15 September 2017 and very little (just 12.6 mm) between 3–13 September.

The surface water level state in the Werenskioldbreen forefield on 15 September 2017 can therefore be regarded as a stable mean and representative. The inference from this is that the volume of groundwater flow out of the forefield makes up an average of 2% of the water flowing in the Breelva (31 mm for the entire catchment area). Based on this assumption, it was estimated that for the 85-day monitoring period, the outflow from the catchment area via the Breelva was c. 68 mln m³, and that the volume of outflowing groundwater during this period was c. 1.36 mln m³.

Piechota et al. (2009–2011) [38] obtained similar results when modelling the Werenskioldbreen catchment area. Despite the considerable generalisation of the unglaciated zone, the groundwater flow through the Werenskioldbreen forefield was 3% of recharge from the basin and 8% for the whole catchment area including the glacierised part. It should be added here, however, that the area that those authors modelled was much larger—36.2 km², of which 25.7 km² (in 2017 year) [37] was glacierised—and that the active layer was 2 m thick. Lower generalised hydraulic conductivities were then applied (5 × 10⁻⁵ m s⁻¹ and 1.5 × 10⁻⁵ m s⁻¹), as a result of which a slower mean groundwater flow rate was obtained for the modelled catchment area, i.e., 4624 m^3/d [38]. The extensive network of rivers and surface water bodies within the forefield were not taken into consideration either; hence, the modelling results display differences in the model’s recharge, where 93.4% consists of rainfall and ablation waters, defined as the overall surface recharge, and 6.6% involves water supplied from the river network [38]. In the present model of the forefield, the river network and surface water bodies are the main sources of groundwater recharge in the model (73.5%), whereas in the ablation season, rainfall contributes only 26.5% to the recharge.

Research in the Finstenwalderbreen basin (Spitsbergen, 1999) [25,26] showed that groundwater made up just 1% (11 mm) of the outflow from the catchment area. Our catchment area was much larger than that analysed by R. Cooper et al.’s research team (65.7 km², of which 43.5 km² was glacierised) [25,26]. The annual outflow from their total catchment area was 1073 mm, whereas from our catchment area, it was 1542 mm.

Later research by Marciniak et al. (2007–2009) [28,29] for the Ebba River (Petunia Bay, Central Spitsbergen) highlighted the possible variability in the proportion of groundwater in the outflow from the catchment area in the ablation season. They used a mathematical model for their calculations using Darcy’s law, applying a similar range of hydraulic conductivities, and obtained the following results, showing that the estimated average percentage of groundwater runoff in the Ebba River was:

• 2007—Estimated: 2.0% groundwater; calculated: average = 4.57% (0.11 m³/s), min = 0.98%, max = 15.05%;
• 2008—Estimated: 3.8% groundwater; calculated average = 8.42% (0.14 m³/s), min = 1.36%, max = 27.36%;
• 2009—Estimated: 3.1% groundwater; calculated average = 5.66% (0.15 m³/s), min = 1.06%, max = 13.63%.

Those authors stated unequivocally that according to their calculations, the mean percentage of groundwater in the Ebba River ranged from 2.0% to 3.8%, whereas the maximum could be as much as 15–27% of that river’s runoff [29]. This is quite a high result, bearing in mind that they assumed a similar value for the hydraulic conductivity and took into account the possibility of water flowing within the active layer at depths not exceeding 1 m. That result was taken to agree with periodic high water levels in the river, which were provided on the basis of flow-rate measurements made along a specified stretch of the river.

Summarising, one cannot rule out the possibility that in the Werenskioldbreen basin, the exchange of water with respect to the Breelva River’s long-term, permanent channel may be very much greater (because of possible deeper infiltration) when water levels are high and surface water runoff during the ablation period has reached a maximum. Confirmation of this hypothesis will require a better understanding of the infiltration capabilities and thickness of the active layer in the river channel.

6. Summary and Conclusions

This paper discusses the hydrogeological modelling results for a glacier forefield (Werenskioldbreen, Svalbard). The aim of the groundwater modelling was to determine the subsurface flows and to estimate the components of the water balance in an active layer of glacial sediments to 2.5 m thickness, where water infiltration is possible.

The indirect objective was to identify the hydrogeological properties of the sediments within the forefield. Field and detailed laboratory studies were conducted to analyse the granulometric composition and hydraulic conductivity of the sediments. In addition, because there have been very few modelling studies on the importance of groundwater in the Arctic basin, the feasibility of using FEFLOW software in this environment was tested to facilitate the spatial representation and analysis of the results.

The main findings from this study can be summarised as follows:

• The sedimentological structure of the Werenskioldbreen forefield consists mainly of gravelly sand (31.8%) or sandy gravel (22.7%) and silty sand (18.1%). The calculated hydraulic conductivities are typical of glacial sediments. Nevertheless, detailed analysis of samples sometimes revealed hydraulic conductivities lower than those reported in the literature. This is most probably due to their poor sorting, as they often contain mixed finer fractions, which are the remains of tills from the outwash. The hydraulic conductivities in our sediments were: sandy gravel (1.60 × 10⁻³ m s⁻¹ and 1.05 × 10⁻³ m s⁻¹), gravelly sand (3.6 × 10⁻⁴ m s⁻¹ and 4.6 × 10⁻⁴ m s⁻¹ and 6.9 × 10⁻⁴ m s⁻¹), silty sand (3 × 10⁻⁵ m s⁻¹), sandy clayey silt (1.64 × 10⁻⁸ m s⁻¹) and weathered crystalline rock (8 × 10⁻⁶ m s⁻¹ and 7 × 10⁻⁶ m s⁻¹).
• The groundwater table in the model area is perfectly consistent with the terrain’s morphology, typical of shallow groundwater depths from a few centimetres to 0.5 m in plains. The model’s generalisation did not take into consideration periodic changes in the river network or locally occurring flooded areas when surface water levels were high. On the other hand, these areas are shown as modelling results where the groundwater level is above the ground elevation. This indicates that the model calculations are consistent with the real environment.
• The model’s results indicate the existence of a groundwater divide in the central part of the foreland, i.e., that there are two subsystems of shallow northern and southern groundwater circulations. From the north and north-eastern zones, groundwater flows westwards towards the Breelva channel and south-westwards along the end moraine direction to the single gauging point on its southern margin. In the southern part of the basin, groundwater flows towards smaller rivers and lakes, especially Lake Mewie, which drains a significant proportion of the groundwater aquifer. Finally, groundwater also flows in a south-westerly direction towards the Breelva gauge.
• The dynamic resources of the active layer in the model area are equal to 6077 m³ d⁻¹. The groundwater aquifer is recharged by inflow from the Breelva River (18.2%), by inflow from other watercourses and surface water bodies (55.3%), and by effective infiltration from precipitation (26.5%). Drainage from the model is to rivers and lakes. The main drainage component is the Breelva, responsible for 40% of the total outflow from the model, while the remaining rivers and lakes drain 60%.
• The model for steady-state conditions shows that the groundwater is responsible for about 2% of the total runoff from the catchment area. The authors also estimated, based on the monitoring of catchment runoff, that during the 2017 ablation season, the runoff over 85 days was 68 million m². Assuming a 2% contribution of groundwater to the river runoff for each day, the total runoff was estimated at c. 1.36 mln m².
• Under typical hydroclimatic conditions, the water-filled interstices in the layers of the modelled sediments, i.e., the fluid volume (1.0073 × 10⁶ m³), make up some 50% of the void volume (2.0689 × 10⁶ m³) in the forefield area. This demonstrates the significant capacity of sediments to accumulate and transport liquid water. If the active layer thickness is increased by 0.5 m, the groundwater capacity increases by 21%.
• The implementation of hydrogeological models can be very useful for filling the knowledge gap concerning the future evolution of groundwater flow and storage conditions in polar regions, where intense changes associated with permafrost reduction are occurring as a result of climate warming [39]. Our study contributes to improving the understanding of hydrogeological trends in cold regions. Detailed recognition of the hydrogeological conditions of the active layer of sediments provided interesting results on spatial variations in the velocity and direction of groundwater flow and in the active layer storage capacity. Such a detailed approach is still unique in model studies for the Arctic region owing to the difficulties of data availability and the high cost of such research.
  The hydrogeological model of a typical Svalbard glacier forefield can also be a good example of applying such models to other catchments in cold climates, as well as in glacierised catchments. This study also indicates the potential of applying this model in water management research to identify the groundwater component resource and its variability due to permafrost loss. Some regions of cold climates are also supplied with drinking water from the aquifer that forms above the permafrost [85,86]. The development of studies of hydrogeological processes and changes in the thickness of the active layer or permafrost loss, which integrate cryohydrogeology studies [39], enables better water management in these regions.
• A tool like FEFLOW [67] enables the spatial variation in hydrogeological conditions to be accounted for in a model based on detailed field and laboratory results, which significantly impacts the reliability of the results. Moreover, the mesh flexibility of the FEFLOW model reflected the spatial hydrogeological structure, where layer thicknesses are very thin and the surface morphology is complex. Although our model was only produced for steady-state conditions, the software chosen allows us to quite easily extend the model to transient conditions and to add modules responsible for heat and/or mass transport in the future once appropriate input data from the research area have been collected. The FEFLOW tool belongs to the DHI software family; therefore, it can be integrated with the hydrological Mike HydroRiver model and thus combine the modelling of surface and groundwater hydrology.