The Importance of Solving Subglaciar Hydrology in Modeling Glacier Retreat: A Case of Study for Hansbreen, Svalbard

1. Introduction

Global warming and Arctic atlantification are driving significant changes in the polar regions [1,2], with implications for the stability of Arctic tidewater glaciers [3,4]. The retreat of these glaciers serves as a critical climate indicator, reflecting the impacts of rising atmospheric and oceanic temperatures [5,6,7,8,9].

The primary mechanisms of glacier mass loss include surface melting, iceberg calving and submarine melting [10,11,12], all of which interact through different processes to influence glacier dynamics (Figure 1). During the melt season, positive degree days, characterized by atmospheric temperatures above freezing, lead to surface melting in Arctic tidewater glaciers [9]. As the glacier surface melts, the resulting meltwater percolates through the ice and reaches the subglacial bedrock, flowing toward the glacier terminus. This meltwater is subglacially discharged into the denser fjord waters, then forming a buoyant plume, which entrains surrounding ambient fjord waters and carries sediments up to the fjord surface [13,14,15]. The rising plume generates overturning circulation within the fjord, effectively drawing warmer out-fjord waters toward the glacier terminus [16,17,18]. The influx of these warmer waters along with the flux of subglacial discharges enhances submarine melting at the glacier terminus [19,20,21,22,23], which in turn amplifies calving events [7,17,24,25,26,27] (see schematics in Figure 1).

The complex interplay between surface melting, subglacial flow, and fjord dynamics, along with their dependence on ambient factors such as atmospheric conditions, ocean temperatures, and bedrock topography, highlights the challenges in identifying the mechanisms driving glacier retreat in response to changing conditions. In this context, we pay special attention to the subglacial hydrology system, which is crucial for determining the number, width, and flow rate of discharging channels [28,29,30], thereby playing a key role in frontal ablation [31]. Moreover, the flow of meltwater beneath glaciers plays a critical role in influencing ice dynamics, including ice flow velocity and stability [32]. As the atmosphere warms during the melt season, increased surface melting generates more meltwater, which can accumulate at the glacier base [33]. This meltwater can lubricate the bedrock, accelerating ice flow toward the terminus [34,35,36], thus favoring calving by shear stress.

The formation of subglacial channels is a dynamic process that evolves throughout the melt season [37,38]. As meltwater accumulates on the glacier’s surface, it flows to the bedrock through crevasses and moulins, creating a drainage system [36,39,40,41]. Initially, water is distributed in thin films or small passageways, but as more meltwater enters, it concentrates along specific flow paths influenced by bedrock topography and pressure gradients [42]. This flow causes thermal erosion of the ice, creating a positive feedback loop: faster water flow generates heat, melts more ice, and enlarges the channel, leading to the development of larger, more efficient drainage channels over time [36,43,44].

However, despite its potential contribution to glacier retreat, subglacial hydrology represents one of the most significant yet poorly understood components of glacier dynamics. This complexity arises from the challenging access to subglacial environments, which are typically located beneath thick ice masses, making direct observations and measurements exceedingly difficult. The harsh conditions, including extreme cold, high pressure, and the dynamic nature of the ice above, pose substantial logistical challenges to study these hidden systems. The scarcity of observations in subglacial environments has led to a reliance on models to infer the behavior of these systems [41,42].

Modeling glacier dynamics is essential for enhancing our understanding of Arctic glacier-fjord interactions by simulating the complex physical processes that govern these systems [27,42]. Recent studies have incorporated information of the subglacial system into glacier-ocean coupled models [45,46,47,48], recognizing its importance in glacier dynamics. However, when it comes to tidewater mass loss, it is crucial to quantify the trade-off between the increased computational cost of solving a subglacial hydrology model and the results obtained from models that simplify or parameterize the subglacial component.

In this study, we aim to assess the impact of subglacial hydrology on glacier front retreat simulations. We focus on the Hansbreen glacier, which offers a rich dataset for model testing, making it an ideal natural laboratory for our research. We use a 3D glacier dynamics-fjord model which accounts for surface mass balance, calving, and submarine melting. We then compare the outputs from two versions of the model: one that resolves the subglacial hydrology model and another that operates under the assumption of a single discharge channel at the flow line. By quantifying the differences in annual front retreat, submarine melting, and calving between these two models, we evaluate the role of subglacial hydrology in understanding glacier dynamics, which can enhance our ability to make more accurate predictions of future changes in response to ongoing atmospheric and oceanic warming.

2. Materials and Methods

2.1. The Glacier-Fjord System of Study: Hansbreen-Hansbukta

The glacier-fjord system Hansbreen-Hansbukta is located in one of the branches of Hornsund fjord in South West Spitsbergen, Svalbard, at ∼77°N, ∼15.6°E (Figure 2). Hansbreen is a polythermal tidewater glacier flowing southward that covers an area of ∼57 km². It is about 16 km long with a low mean surface slope of around 1.8° on average along the central flowline [49]. Its calving front is 1.5 km wide with a vertical face ∼100 m thick at the central flowline, of which 50 to 60 m are below the sea level. The seasonal retreat of Hansbreen usually starts in June/July and lasts until late autumn/early winter, and the average summer and winter fluctuations amount to $-125$ and 79 m, respectively [50]. As for Hansbukta, it is a ∼2 km long bay, with a maximum depth of ∼77 m. Temperature and salinity in Hansbukta experiences strong seasonal variability, ranging from $-1.8$ to 3°C and from 34.6 to 31.8 psu between April and August, respectively.

2.2. Observations

There are a set of observations needed to constrain the glacer-fjord model: glacier velocities, surface meltwater, bedrock topography, glacier front position, and fjord properties.

We use the gridded surface velocity data of [45], which were obtained by applying Bayesian Kriging techniques [52] on daily horizontal velocities measured at a set of stakes located along the glacier (see blue dots in Figure 2).

Surface mass balance (SMB) and surface meltwater (SMW) (Figure 3) were obtained from downscaled European Arctic Reanalysis data at 2 km horizontal and hourly temporal resolutions, constrained by automatic weather stations and stake observations [53]. The surface elevation came from the SPIRIT digital elevation model for gentle slopes, with a 30 m RMS absolute horizontal precision and 40 m resolution. Bedrock topography was inferred from ground-penetrating radar data [49,54].

We take the Hansbreen weekly front positions obtained from time-lapse camera images processed and described in [55]. We also analyzed all the imagery searching for plume-activity locations. In some of the winter images, plume activity was observed as a free-of-ice patch of the fjord surface attached to the glacier terminus. In summer periods, the imprint of the plume was recognised by a sediment patch at the fjord surface close to the glacier front (Figure 4).

Hydrography of Hansbukta has been monitored from 2010 up to date by the IOPAS (Institute of Oceanology Polish Academy of Science) and the data used in this study were already published in [17,51], as well as used in [45]. We give here some information of the fjord data used in this study. CTD (conductivity, temperature and depth profiler) data were collected at different casts in Hansbukta (yellow points in Fig. Figure 1). Temperature and salinity measurements were vertically averaged every 1 dbar (1 kPa), giving a vertical resolution of 1 m.

2.3. The Models

We use the coupled glacier-fjord model of [45], which was built on the open-source, full-Stokes, finite element, ice flow model Elmer/Ice [56]. [45] coupled the different processes of the system: glacier dynamics, calving, and subglacial hydrology models for the glacier component, as well as line plume and submarine melt models for the fjord component. The different models coupled are briefly explained in the following subsections, though we refer to [45] for exhaustive information.

2.3.1. Glacier Dynamics + Calving

The model solves the full-Stokes equations for ice flow, with rheology defined by Glen’s flow law (e.g., [57]) and uses the calving implementation described by [46,47,58,59]. This implementation was initially postulated by [60,61], and later improved by [62,63] for use in a 3D framework. This calving criterion has been chosen because Benn et al. [64] demonstrated that the crevasse depth (CD) calving law reflects the glaciological controls on calving and exhibits considerable skill in simulating seasonal fluctuations. CD are calculated by

$$\begin{array}{c}\hfill {\sigma }_n=2{\tau }_{e}sgn\left({\tau }_{xx}\right)-{\rho }_{i}gd+P_w\end{array}$$

(1)

where ${\sigma }_n$ is the net stress (positive for extension and negative for compression). The terms on the right-hand side represents the balance of forces: the first corresponds to the opening force of longitudinal stretching, where ${\tau }_e$ represents the effective stress, ${\tau }_e^2={\tau }_{xx}^2+{\tau }_{zx}^2$ and the sign function ($sgn$) ensures that crevasses opening is only produced under longitudinal extension; the second term corresponds to the ice overburden pressure, which leads to creep closure, where ${\rho }_i$ is the ice density, g is the gravity acceleration and d the crevasse depth. $P_w$ stands for the water pressure inside crevasses, which is here considered to be zero for surface crevasses and controlled by the subglacial hydrological system for basal crevasses, and at the calving front can be expressed as:

$$\begin{array}{c}\hfill P_w=(Z_{sl}-Z){\rho }_{w}g\end{array}$$

(2)

being ${\rho }_w$ the density of water at the calving front and Z the elevation with respect to sea level. The sea-level elevation, $Z_{sl}$, is set to 0 m. When calving occurs, the model calculates a calving vector which is normal to the calving front, maps pre-calving to post-calving node positions, and remesh the main mesh.

2.3.2. Buoyant Plume + Submarine Melting

Due to the density differences between meltwater and fjord water, subglacially discharged water rises in contact with the calving front, mixing turbulently with the surrounding water and producing melting at the ice–water interface.

We rely on the line plume model of [21,65] to simulate the buoyant plume generated by sublacial discharge of meltwater entering the fjord. This line plume model has shown to accurately reproduce this phenomenon in tidewater glacier terminus [28,66], and is implemented in Elmer/Ice [46,47] as a continuous sheet-style ’line’ plume, split into coterminous segments across the calving front. To compute the submarine melt rates at the glacier terminus, the thermodynamical equilibrium at the ice-ocean interface [67] is solved (see [17] for further details).

The plume model is initialized by the subglacial discharge at each node of the grounding-line, where, depending on the experiment (Section 2.4), the subglacial discharge values are obtained either as a solution of the subglacial hydrology model (Section 2.3.3), or within a single imposed channel centered with the flow line. The calculated melt rates are then applied to modify the geometry of the submerged part of the calving front.

2.3.3. Subglacial Hydrology

The subglacial hydrology of Hansbreen is simulated by using the Glacier Drainage System module (GlaDS) of Elmer/Ice [68]. GlaDS simulates both inefficient distributed drainage, represented by a sheet of water that covers the whole area of the glacier, and efficient channelized drainage, represented by a series of channels generated on the edges of the mesh elements of the domain (see more detail in [69]). The main parameters of the model are set out in Table 1.

We use GladS to obtain subglacial discharge estimates at the grounding line. Water is not permitted to flow through the lateral boundaries, and we set the hydraulic potential, $\varphi$, to zero at the grounding line. From the combination of Eq. (2) and the hydraulic potential, we get

$$\begin{array}{c}\hfill \varphi ={\rho }_{w}gZ+P_w\end{array}$$

(3)

Water entering the hydrological system is mainly derived from surface meltwater (Figure 3). We account also for a small portion of basal meltwater, which is assumed to be spread and is calculated using a geothermal heat flux of 63 mW m⁻² [68].

2.4. Experiments

We simulate the glacier-fiord system of Hansbreen-Hansbukta from April to December 2010. Two distinct models were employed, both incorporating the same components: glacier dynamics, calving, buoyant plume, and submarine melting. The primary difference between the models lies in the treatment of subglacial hydrology.

• Model 1 is coupled with the subglacial hydrology model, allowing meltwater to flow freely through the discharge channels derived from the simulation. This approach enables a more realistic representation of the interactions between the glacier dynamics and the underlying hydrological system.
• Model 2 does not resolve the subglacial hydrology component. Instead, it imposes a single discharge channel, centered along the glacier’s flow line, with an approximate width of 200 meters, as described by [17,51]. With this ’conventional’ configutration, we test the impact of simplifying subglacial hydrology on the overall dynamics of the glacier-fjord system.

By comparing the results from both models, we aim to gain insights into the role of subglacial hydrology in glacier front retrat and the dynamics of Hansbreen-Hansbukta system.

3. Results

3.1. Discharging Channels from the suBglacial Hydrology Model

The subglacial hydrology model is coupled with the glacier dynamics and mass balance of Hansbreen. The subglacial water drainage system evolves dynamically over time, influenced by factors such as surface meltwater input, bedrock topography, and water-induced erosion within the ice. Figure 5 presents the modeled configuration of the subglacial drainage system in September, at the end of the melt season, when the drainage system is considered to be well established and more efficient.

The subglacial hydrology model identifies three distinct channels, each exhibiting unique characteristics. The largest channel, approximately 100 m in width, is situated in the deepest region ($\sim 70$ m depth) on the central-western side of the Hansbreen terminus (highlighted in orange and labeled as channel 2 in Figure 5). The other two channels, each approximately 50 m wide and shallower ($\sim 30$ m depth), are located closer to the glacier margins (identified as channels 1 and 3 in Figure 5). This three-channel configuration predicted by the model corresponds closely with observational data. Time-lapse imagery (Figure 4) shows distinct plume patches at the glacier terminus, which are consistent with the positions of the modeled channels, indicating the presence of subglacial water discharge.

3.2. Performance of the Two Models Against Observed Glacier Terminus

After running Model 1 and Model 2 from April to December 2010, we obtain the various positions of the Hansbreen front over time. Figure 6 graphically represents the observed and modeled monthly positions for a qualitative analysis between Models 1 and 2. On the other hand, Table 2 quantifies the longitudinal differences of each model against observations. This quantification results from calculating the area difference between observed and modeled, divided by the observed width of the glacier front. Negative values indicate that observed front is more advanced than modeled one.

The results in Figure 6 show a very good approximation of both models to the observations during the initial months of the simulation (May and June), during which there are minimal changes in the position and morphology of the glacier front. The longitudinal difference in April is so small ($<1\mathrm{m}$) that it is not graphically represented (see values in Table 2). For the months of May and June, both models exhibit a longitudinal difference of $<5\mathrm{m}$ compared to the observed front. However, it is noteworthy that Model 1 presents a glacier front that is a bit more retreated than observed.

Starting in July, this situation changes, and the front begins to retreat more significantly until October, thereby increasing the differences between both models and the observations. These differences become remarkable in the eastern margin from July, and at the center of the glacier terminus from August, and maintained in both cases until October. The longitudinal differences of Model 1 (Model 2) against observations increase from ∼13 m (∼32 m) m in July, to ∼ 76 m (∼122 m) in October (Table 2). Although neither model exactly matches the observed front position, Model 1, which is coupled with the sophisticated subglacial hydrology model, provides results that are much closer to the observations, reducing the longitudinal difference presented by Model 2 by up to 50%.

3.3. Comparison of Frontal Ablation Between the Two Models

Frontal ablation refers to the mass loss of a glacier at its terminus. As previously described in the introduction (Section 1) and illustrated in Figure 1, the front of a tidewater glacier can experience mass loss primarily through two mechanisms: calving, which is associated with the detachment of icebergs, and submarine melting. The latter is influenced not only by the temperature of the fjord but is also enhanced by buoyant plumes generated by subglacial meltwater discharges. In this section, we quantify the frontal ablation associated with these two main mechanisms: calving and submarine melting, for the two models employed in the simulations.

3.3.1. Calving

In general terms, the simulation results indicate that the volume of mass loss due to calving increases throughout the melt season, peaking in August (Figure 7a). Following September, the mass loss from calving begins to decline. This trend is consistent with the pattern observed for surface meltwater (Figure 3), highlighting the interrelationship between these two mass-loss mechanisms.

Focusing on the calving volume presented in Figure 7a, the minimum monthly mass loss rates, approximately $3\times {10}^4{\mathrm{m}}^3{\mathrm{month}}^{-1}$, occur in May for both models. In April and June, Model 2, which assumes a single discharge channel, exhibits a greater calving volume than Model 1. However, this situation reverses in the subsequent months (July, August, and September), when the mass loss from calving is more pronounced. In July and September, the calving volume for Model 2 is around $4\times {10}^5{\mathrm{m}}^3{\mathrm{month}}^{-1}$, while Model 1 reaches calving volumes of up to ${10}^6{\mathrm{m}}^3{\mathrm{month}}^{-1}$. This suggests that the resolution and coupling of the subglacial hydrology model with the glacier dynamics model (Model 1) facilitate more than double the frontal ablation due to calving during July and September, which are months characterized by significant surface melting (see Figure 3). In August, both models exhibit their highest monthly calving volumes, with Model 1 reaching $3\times {10}^6{\mathrm{m}}^3{\mathrm{month}}^{-1}$ and Model 2 reaching $2\times {10}^6{\mathrm{m}}^3{\mathrm{month}}^{-1}$.

Except for the first three months (April to June), the total mass loss accumulated due to calving is greater in Model 1 along the melt season, which reaches a total of $6\times {10}^6{\mathrm{m}}^3$ (∼6 Mt) by November (Figure 7b). In Model 2, the mass loss from calving amounts to a total of $4\times {10}^6{\mathrm{m}}^3$ (∼4 Mt), which is one-third less than that of Model 1. This highlights the significant importance of incorporating the subglacial hydrology model into the glacier dynamics model when calculating frontal ablation due to calving.

3.3.2. Submarine Melting

Submarine melting refers to the submerged portion of a glacier terminus that melts due to direct contact with seawater or fjord waters. In our simulations of Hansbreen glacier, we computed the monthly submarine ice melt (Figure 8a) and noticed that both models exhibit a similar general pattern to that of calving (Figure 7a) and surface meltwater (Figure 3) throughout the simulation period.

In April, submarine melting is minimal, coinciding with cold temperatures in the fjord and low subglacial discharge. This melting progressively increases, reaching its peak in July, which occurs one month prior to the maximum calving volume (see Figure 7a). Following this peak, submarine melting decreases again until the end of the simulation.

The results indicate that monthly submarine melting rates are consistently higher in Model 1 compared to Model 2 throughout the simulation, with the exception of July, when melting reaches its maximum of 5000 m³ month⁻¹ in both models. Figure 8b illustrates the accumulation of frontal ablation resulting from submarine melting. It is evident that Model 1 generates greater mass loss due to submarine melting than Model 2 for the entire duration of the simulation. By the end of the simulation, the total ice-volume loss attributed to submarine melting is 15000 m³ in Model 1 and 12500 m³ in Model 2.

4. Discussion

4.1. Subglacial Hydrology as an Important Component of the Glacier-Fjord System

Subglacial hydrology plays a crucial role in the dynamics of tidewater glaciers. It is interconnected with various processes, relying on surface meltwater while also influencing ice velocities, subglacial discharge, submarine melting, and iceberg calving [15,23,70,71,72]. Understanding these interactions is essential for accurately modeling glacier behavior [42] and predicting future changes in response to climate variability.

Our study demonstrates that resolving the subglacial hydrology model significantly enhances the accuracy of predicting the position of the glacier front (Figure 6 and Table 2). By resolving the subglacial water flow, we can account for the complex network of drainage channels that develop beneath the glacier (Figure 5). This complexity allows for a more realistic simulation of how water moves and exerts pressure at the glacier base, altering the basal sliding coefficient and ultimately affecting glacier velocity and frontal stability [32,34,35].

These additional channels facilitate greater subglacial water flow, which aligns with sediment plume observations in our study (Figure 4) and has also been reported by [15,71]. This increased flow contributes to higher submarine melting rates (Figure 8), consistent with the findings of [28,72]. Moreover, the presence of multiple drainage channels amplifies mass loss through calving (Figure 7) as a consequence of enhanced net stress differences at the glacier terminus [51].

In contrast, simplifying subglacial hydrology by assuming a single discharge channel at the tidewater glacier terminus—an approach commonly found in the scientific literature [e.g., [17,21,22,51,73]—leads to less accurate predictions of the glacier front position, as shown in Figure 6. In our study, the results from the oversimplified model indicate that the glacier front is further advanced than both the comprehensive model and actual observations. This suggests that subglacial hydrology is a crucial component in modeling the entire glacier-fjord system.

This discrepancy arises because the simplified model is limited to a single discharge channel, which fails to account for the multiple drainage pathways present at the glacier front. Consequently, it does not accurately capture the pressure distribution along the entire base of the glacier. This results in smaller differences in glacier velocities at the front, which hinders the destabilization of the glacier front and promotes less calving (Figure 7). Ultimately, this leads to a glacier terminus that is more advanced than what is actually observed (Figure 6).

In summary, the interaction between subglacial water dynamics and frontal ablation is essential for understanding the overall mass balance of Arctic tidewater glaciers.

4.2. Implications of Simplifying the Subglacial Hydrology Component in the Modeling Context

The complex interplay between surface meltwater, subglacial discharge, and iceberg calving highlights the essential requirement for precise modeling of subglacial hydrology [42,46,47]. Our research demonstrates that in months characterized by reduced surface meltwater production (Figure 3), the variance in glacier front position between the hydrology-resolving model and the simplified model is negligible (Figure 6 and Table 2). Conversely, as surface meltwater production increases, the discrepancies between the two modeling approaches become increasingly significant (Figure 6 and Table 2). Importantly, the model that incorporates subglacial hydrology exhibits a closer alignment with observations, indicating its superior accuracy in representing the underlying physical processes.

In our simulation, which encompasses a single year from April to December 2010, we observed a notable difference in mass loss attributed to calving between the two modeling approaches (Figure 7). The coupled model indicated a mass loss of 6 Mt, whereas the simplified model recorded a significantly lower mass loss of 4 Mt. This discrepancy underscores the importance of subglacial hydrology, particularly as surface meltwater increases. Consequently, in the context of rising atmospheric and oceanic temperatures, and the high sensitivity of tidewater glaciers to these changes [6,7,8,50], the reliance on predictive models that oversimplify subglacial hydrological processes may result in a considerable underestimation of glacier retreat rates in future climate change scenarios.

While simplified models present advantages in terms of reduced computational costs when studying glacier dynamics [51,61,74,75], it is crucial to consider these benefits in light of their inherent limitations. To properly adjust the parameterization of the subglacial hydrology, further comparative simulations of both modeling approaches (simplified and resolved) are necessary, using empirical observations from a diverse range of glacier-fjord systems with varying characteristics and under different atmospheric and oceanic conditions. This approach will provide a more comprehensive understanding of how each model responds to specific environmental factors. Furthermore, such studies would facilitate the incorporation of an appropriate error range for the simplified models, thereby addressing their limitations and enabling a more accurate assessment of mass loss and front retreat in Arctic glaciers.

Furthermore, the incorporation of artificial intelligence methodologies into glaciological research is increasingly being recognized as a valuable advancement for mapping glacier basal sliding [76], detecting glacier-snow lines [77], and modeling glacier melt dynamics [78] or surface mass balance (e.g., [79,80,81]). These techniques require training on observational datasets or data generated from existing models. Consequently, to ensure the reliability of predictions derived from these advanced methodologies, it is imperative to produce model results that are both accurate and representative of the real glacier-fjord system. This foundational aspect will significantly enhance the credibility of future forecasts concerning glacier dynamics and their responses to climate change.

5. Conclusions

This study highlights the critical role of subglacial hydrology in the dynamics of tidewater glaciers, emphasizing its influence on glacier front position, submarine melting, and iceberg calving. Our findings demonstrate that coupling the glacier dynamics model with a resolved subglacial hydrology model significantly improves the accuracy of predicting glacier front positions compared to the simplified model that assumes a single discharge channel. The complexity of subglacial drainage networks seems essential for capturing the pressure distribution at the glacier base, which directly affects glacier stability and mass loss through calving.

As surface meltwater production increases, the discrepancies between the comprehensive and simplified models become more pronounced. Within a single melt season, the simplified model underestimates up to 1/3 of the mass loss due to calving, underscoring the necessity of incorporating subglacial hydrology in predictive models. The observed differences in mass loss due to submarine melting and calving further illustrate the importance of accurate modeling in understanding glacier responses to climate change.

In light of rising temperatures and the sensitivity of tidewater glaciers, reliance on oversimplified models may lead to significant underestimations of glacier retreat rates. Future research should focus on comparative simulations across diverse glacier-fjord systems to refine parameterizations of subglacial hydrology and enhance the reliability of predictions. Additionally, properly integrating artificial intelligence methodologies into glaciological research holds promise for reducing computational effort for glacier dynamics modeling.

Our findings highlight that a comprehensive approach is essential for improving predictions of glacier behavior and understanding the broader implications for sea-level rise and climate change.