Geometric Controls of Tidewater Glacier Dynamics

Retreat of marine outlet glaciers often initiates depletion of inland ice through dynamic adjustments of the upstream glacier. The local topography of a fjord may promote or inhibit such retreat, and therefore fjord geometry constitutes a critical control on ice sheet mass balance. To quantify the processes of ice-topography interactions and enhance the understanding of the dynamics involved, we analyze a multitude of topographic fjord settings and scenarios using the Ice-sheet and Sea-level System Model (ISSM). We systematically study glacier retreat through a variety of artificial fjord geometries and quantify the 5 modeled dynamics directly in relation to topographic features. We find that retreat in an upstream widening or deepening fjord does not necessarily promote retreat, as suggested by previous studies. Conversely, it may stabilize a glacier because converging ice flow towards a constriction enhances lateral and basal shear stress gradients. An upstream narrowing or shoaling fjord, in turn, may promote retreat since fjord walls or bed provide little stability to the glacier where ice flow diverges. Furthermore, we identify distinct quantitative relationships directly linking grounding line discharge and retreat rate to fjord topography, 10 and transfer these results to a long-term study of the retreat of Jakobshavn Isbræ. These findings offer new perspectives on ice-topography interactions, and give guidance to an ad-hoc assessment of future topographically induced ice loss based on knowledge of the upstream fjord geometry.

on paleo-time scales, retreat has occurred over large distances, but the temporal resolution of geomorphological studies is limited by the available geological data and key information is missing to discern different drivers of glacier retreat (Briner et al., 2009;Åkesson et al., 2020). Meanwhile, numerical studies that can address these issues have so far mostly used widthand depth-integrated flowline models, which carry many assumptions that do not hold in some settings (Nick et al., 2009;60 Åkesson et al., 2018b;Steiger et al., 2018). In particular, they parameterize or do not account for factors that are thought to be instrumental to explain ice-topography interactions, such as lateral drag, across-flow variations in glacier characteristics and viscosity changes due to variations in ice temperature. The latter, for instance, was found to be key to explain Jakobshavn Isbrae's recent retreat (Bondzio et al., 2017).
Here, we use a numerical ice-flow model resolving two horizontal dimensions, we include a suite of 21 experiments and 65 present a systematic approach to compare the relative importance of basal and lateral fjord topography. This setup allows to assess how fjord topography controls glacier retreat on interannual to centennial time scales. We hypothesize that there are quantifiable relationships between glacier retreat and topography that apply to a wide range of glaciological settings. Such general relationships would yield substantial predictive power for a broad assessment of expected future outlet glacier retreat.
We create a large ensemble of artificial fjords that include geometric features (referred to as 'perturbations' in the following) 70 typically found in glacial fjords, such as sills and overdeepenings (referred to as 'bumps ' and 'depressions', respectively; together 'basal perturbations') as well as narrow and wide fjord sections (referred to as 'bottlenecks' and 'embayments', respectively; together 'lateral perturbations'). We then force synthetic glaciers to retreat through this variety of fjords by increasing ocean-induced melt rates, and assess key retreat metrics such as the grounding line retreat rate. The ice dynamics of each simulation are compared and quantitatively linked to the characteristics of the respective fjord geometry. Finally, we 75 investigate whether the retreat dynamics and ice-topography interactions are transferable from the idealized setup to a longterm study on Jakobshavn Isbrae in western Greenland.

Ice sheet model
We use the Ice-sheet and Sea-level System Model (ISSM; Larour et al., 2012) with the shallow-shelf approximation (SSA; 80 Morland, 1987;MacAyeal, 1989). A discussion on the appropriateness of this approximation compared to a Full-Stokes model is provided in sect. 4.2. Our domain is rectangular (80 km×10 km) with x and y representing the along-flow and acrossflow coordinates, respectively (Fig. 1a). We create an unstructured mesh with a fixed resolution of 100 m close to the GL, comparable to other high-resolution studies of Greenlandic fjords (e.g. Morlighem et al., 2016Morlighem et al., , 2019. The temporal resolution is set to ∆t = 0.01 yrs (3.65 days) to satisfy the Courant-Friedrichs-Levy condition (Courant et al., 1928). We apply a sub-85 element GL migration scheme (Seroussi et al., 2014) and enable a moving calving front with the level-set method . We use a thermal model to account for the feedback between frictional heating and ice viscosity. The spin-up ice temperature is -5 • C, representative of conditions in South/Central Greenland and the Southern Fennoscandian ice sheet (Nick et al., 2013;Åkesson et al., 2018a). To simulate calving, a von-Mises law is used , where the calving rate c is given as: where v is the velocity vector,σ a scalar quantity representing the tensile stress at the ice front, and σ max is a stress threshold.
This formulation, demonstrated to perform well in Greenlandic fjords (Choi et al., 2018), means that calving occurs when the tensile stress at the glacier front exceeds a fixed threshold. σ max generally needs to be determined experimentally for studies on real-world glaciers. Here we fix σ max to 1 MPa for grounded ice and 200 kPa for floating ice. These values yield 95 a representative setup and are within the range suggested for Greenland outlet glaciers Choi et al., 2018Choi et al., , 2021. Note that the choice of the calving law is an important, yet poorly constrained, control on the grounding line dynamics and ice front behaviour (Schoof et al., 2017;Haseloff and Sergienko, 2018). In the absence of a universal calving law, a reasonable assumption has to be made, and we justify the choice of the von-Mises law through its relatively good performance in real-world applications (Choi et al., 2018).

100
Basal sliding is parameterized with a Budd type friction law (Budd et al., 1984) of the form where τ b is basal drag, k is a friction parameter, v b is the the basal velocity, and N is the effective pressure. N is given as N = ρ i gH − ρ w g max(0, −z B ) where ρ i and ρ w are the density of ice and salt water, respectively, g is gravitational acceleration, H is ice thickness and z B is bed elevation with respect to sea level. N is thus the difference between the ice overburden and water 105 pressure assuming perfect connectivity between the subglacial water layer and the ocean. The effective pressure in the friction law induces an elevation-dependence of the basal resistance to flow. This dependence is motivated through the assumption that weak sediments are more likely to be present in low-lying areas. Implications for our results are discussed in sect. 4.2. We set k spatially uniform to isolate the topographic signal of retreat in our results and thus to reduce the number of degrees of freedom for the interpretation. We fix k = 40 (Pa yr m −1 ) 1 2 , which is mid-range among values typically found in glaciological settings 110 resembling ours (Bondzio et al., 2017;Haubner et al., 2018;Åkesson et al., 2018a).
Melting under floating ice is parameterized through prescribed melt rates that are invariant of depth. The reference forcing used for model spin-up is a subshelf melt rate of 30 m yr −1 and a frontal rate of undercutting of 200 m yr −1 , where the latter acts on the vertical ice front if the terminus is grounded. Both values are on the lower end of observed melt rates (Motyka et al., 2003;Xu et al., 2013), thus reflecting a configuration prior to recent warming when glaciers were 115 more in balance with the ambient climate than today (King et al., 2020).
As part of the idealized setup, the surface mass balance (SMB) is fixed to zero, except in the uppermost 10 km of the domain ( Fig. 1a), where we set an accumulation rate of 55 m yr −1 . This creates a realistic fixed upstream ice flux and is not meant to represent local SMB found on real glaciers. Additionally, mass is added to the model domain by fixing an ice velocity of v x = 50 m yr −1 at the upper domain boundary, which creates an influx as a function of the glacier thickness. These two 120 approaches to adding mass represent a total accumulation (A) which can be expressed as where C = 10 × 10 3 2 × 55 = 5.5 × 10 9 m 3 yr −1 is a constant accumulation. Through the thickness-dependent influx represented by the second term on the right-hand side, we account for a reduction in accumulation for a shrinking glacier, thus parameterizing the SMB -altitude feedback (Harrison et al., 2001;Oerlemans and Nick, 2005). 125 We impose free-slip boundaries at the lateral margins of the domain (where y = 0 and y = y max ), meaning that no mass can leave the system laterally. In summary, the only mass source is at the upstream end of the domain, while the only mass wasting occurs where the glacier is in contact with the ocean (either through calving or melting).

Fjord Geometries
Our reference geometry is a fjord sloping linearly towards the ocean with a wide section in the upstream area from which ice is 130 funneled towards a 5 km wide outlet channel with parallel walls (Fig. 1a). The fjord topography is given by and

135
where The parameter values and descriptions are found in Table 1. This formulation is inspired by the MISMIP setup (Gudmundsson et al., 2012;Asay-Davis et al., 2016), but adapted to our purpose.
To insert basal or lateral perturbations in the outlet channel and thus alter the fjords' depth or width in specific areas, we modify the parameter Θ(x) in either eq. 4 or eq. 5 such that 140 Altering Θ(x) only in one of the terms on the right-hand side of eq. 5 allows to produce fjords with one-sided lateral perturbations, thus making them asymmetric. Combined, these equations reproduce the typical U-shape of fjords (Fig. 1b) and yield a setting that is representative of a wide range of outlet glaciers.
The metric used to quantitatively link fjord shape with glacier retreat is the wetted area S: the submerged cross-sectional area of the fjord (Fig. 1b), which can be calculated at every point along an outlet channel. This metric combines information 145 about the width and depth of a fjord and is thus a comprehensible parameter for comparing basal and lateral perturbations. Furthermore, it is straightforward to calculate its first derivative, dS, which yields information on the along-fjord change in width or depth.
Besides our reference setup, we test 20 fjord geometries ( Table 2), each of which contains either a small, medium or large geometric perturbation. The magnitude of the perturbation is defined by how much the width or depth of the fjord deviates from 150 the reference fjord. Our 'core experiment', which the results will focus on, comprises 12 fjords, each of which features one of the four perturbation types (depression, bump, bottleneck, embayment) of one of the three magnitude classes (small, medium, large). The depressions and embayments in each magnitude class increase the wetted area S at the center of the perturbation x C by the same amount, while the bottlenecks and bumps in each magnitude class reduce S at x C by the same amount. The along-flow horizontal extent of all perturbations in the core experiment is 20 km (Fig. 2).

155
In the eight simulations outside our core experiments, we test asymmetric and longer perturbations to verify if the results from our core experiment can be transferred to a wider range of settings. We test two asymmetric embayments, which have the same S at x C as the small and medium embayments and depressions, as well as two asymmetric bottlenecks which have the

Reference Glacier
All experiments start from an artificial reference glacier, which is produced by relaxing a rectangular block of ice in the reference fjord. The spin-up is run until the relative ice volume ((dV /yr)/V << 0.05%) and GL position are stable. The length of the spun-up reference glacier is 82 km, its GL is at x = 73 km and the velocity at the GL along the central flow line of the glacier v GL = 2.5 km yr −1 . At steady-state, the total mass gain is ∼ 6.1 km 3 yr −1 (≈ 5.6 Gt yr −1 ), which is balanced 165 by mass losses through melting at the ice -ocean interface (∼ 0.9 km 3 yr −1 ) and calving (∼ 5.2 km 3 yr −1 ). In a sensitivity experiment with doubled ocean melt rates (400 m yr −1 undercutting and 60 m yr −1 subshelf melting), mass loss through ocean melt increases to ∼ 2.1 km 3 yr −1 , while calving reduces to ∼ 4 km 3 yr −1 . The GL remains largely stable, indicating that the reference glacier is not very sensitive to ocean forcing due to compensating effects in the mass wasting processes.
The setup represents a medium sized fjord-glacier system, which has similar dimensions and dynamics as, for example, the 170 present-day Alison glacier in NW Greenland, where the fjord width is about 5 km, water depth is around 500 m and observed ice discharge has increased from ∼ 4 to ∼ 8 Gt yr −1 in the past 20 yrs (Mouginot et al., 2019). It is furthermore broadly  representative of outlet glaciers from the Fennoscandian ice sheet during the last glacial, such as the Hardangerfjorden glacier (Mangerud et al., 2013;Åkesson et al., 2020).

175
We slightly modify the reference glacier to match the new fjord geometry when introducing geometric perturbations. For embayments, we extrapolate the glacier surface laterally to fill the newly introduced lateral cavities. For depressions, we fill the new basal cavity with ice, but keep the glacier surface the same. For bumps or bottlenecks, we remove ice while keeping the glacier surface unaltered. Subsequently, we relax the glacier in each geometry for 50 yrs, resulting in an ice volume change (dV /yr)/V < 0.5% at the end of relaxation for every setup tested and a stable GL.

180
After relaxation, we increase the ocean forcing to trigger a retreat. We aim to force the GL to retreat through the entire geometric perturbations. The ocean melt rates required to induce such a retreat depend on the fjord geometry, which we elaborate on in the results section. To determine what melt rates are needed to force this complete retreat in a particular fjord, we strengthen the ocean forcing using multiples of the reference forcing (200 m  forcing (e.g. 20 times the reference forcing) did not trigger complete retreat, suggesting that glaciers in these geometries are not sensitive to ocean melt.
Since we want to explore the response of outlet glaciers to melting at the ice-ocean interface, we keep the SMB constant with time, and let the upstream ice flux vary only through our parameterized SMB -altitude feedback (Eq. 3).
We assess 16 glacier metrics during the retreat, which we expect to show a response to local topography ( Table 3). All of 190 these can be observed in-situ or via remote sensing techniques (e.g. Mouginot et al., 2019;King et al., 2020), which means that our results are readily transferable to real-world settings. The GL position (x GL ) and its derivative, the GL retreat rate (dGL), the front position (x F r ) and its derivative, the frontal retreat rate (dF r), but also the velocity at the GL (v GL ) and the shelf length (L S ) are measured along the central flow-line of the glacier.

195
We want to verify to what degree the dynamics seen in our experiments are also prevalent in real-world settings. This is challenging since we investigate decadal to centennial time scales. Specifically, we would need observations with high temporal resolution on glacier metrics (Table 3) for a glacier that has retreated over tens of kilometers, through a fjord with variable and known topography. There are perhaps only a handful of glaciers worldwide that may fulfil these requirements, and even so, acquiring the necessary data is difficult and outside the scope of the present study.

200
To test our idealized results in a real-world setting, we instead turn to model simulations from Jakobshavns Isbrae's (JI) evolution in the Holocene (Kajanto et al., 2020). We focus on a simulation of the retreat of JI from a sill at the fjord mouth of Jakobshavn Isfjord, to a point inland of today's GL position. This model retreat is forced using a step reduction in the equilibrium line altitude early in the Holocene (experiment SE_CM in Kajanto et al. (2020)). While this is a sensitivity exper-

Maximum Ice Thickness
Hmax m iment not meant to reflect the actual evolution of JI (Kajanto et al., 2020), it is convenient for our purposes since it produces a 205 long-lasting, dynamic retreat.
Just like in our idealized experiments, we calculate S and dS for Jakobshavn Isfjord. We then asses whether the relationships found in our idealized settings are also prevalent in JI's dynamics (sect. 3.5).
3 Results In summary, stagnant GL positions are found where the fjord widens/deepens in the direction of glacier retreat (positive dS, Fig. 2) and rapid retreat occurs through areas where the fjord becomes narrower/shallower (negative dS). Thus, the alongfjord change in width or depth (dS) is a key control on GL retreat. However, glaciers in narrower or shallower fjords than the reference fjord (bottlenecks and bumps) can also temporarily stagnate where S is small. This shows that the wetted area 235 constitutes an additional important control on GL retreat. The experiments with asymmetric and longer perturbations confirm these findings (see Fig. SA2).

Forcings and timings of retreat
Now, we investigate how retreat from stagnant and ephemeral GL positions is correlated with fjord topography (i.e., S and dS). Two parameters are important in this context ( Fig. 4): First, the amplitude of the forcings needed to induce complete 240 retreat through the different geometric perturbations. As mentioned previously, distinctions exist both between the different perturbation types (bumps, depressions, bottlenecks, embayments) as well as the different magnitude classes (small, medium, large) for a given geometry type. Second, the approximate residence time of the GL in a stagnant position. The stronger the GL is stabilized by a particular geometric perturbation, the longer it will be stagnating.
All glaciers in embayments require the same increase in forcing to retreat completely (6 times the spin-up forcing). This 245 increase is larger than what is needed to induce retreat through the linear reference fjord (4 times the spin-up forcing). The residence time of the GL in the stagnant positions at the downstream end of the embayments are such that the glacier in the smallest embayment is the earliest to retreat (after 61 yrs of GL stagnation), and the one in the largest the latest (after 173 yrs) ( Fig. 4). This implies that the larger the embayment, the more stability it provides to the glacier at its downstream end, before retreat through the entire perturbation is possible. A larger embayment means a locally larger along-flow change in wetted area  dS at its downstream end (Fig. 2). Thus, there is a positive correlation between GL stability and dS. This indicates that dS not only determines the location of stagnant GL positions in embayments, as shown before, but that it also quantitatively impacts how stagnant the GL is.
The glaciers in fjords with depressions require different forcings to retreat completely (small: 6× the reference forcing, medium: 5×, large: 4×). The residence time also varies; retreat over small depressions occurs ∼65 yrs later than over medium 255 and large depressions, which retreat after about the same time (after 169 and 170 yrs, respectively). These findings indicate that the stabilizing effect of a depression declines the deeper it is (Fig. 4). Thus, there is a negative correlation between GL stagnation and S. In sect. 3.1, ungrounding in the central part of depressions was identified as the trigger of rapid retreat from the temporary stillstands. Such GL retreat occurs more easily in a deeper fjord (larger S). Therefore, it is consequential that a deeper depression is less stabilizing. Note, however, that the fjord depth several kilometers upstream of the GL determines for 260 how long the GL stagnates. There is no direct correlation between S or dS at the GL and the stability provided to the glacier by the fjord in our settings with depressions.
The glacier in a fjord with a 'small' bottleneck required a four-fold increase in oceanic melt rates, and retreated from its stagnant position after 126 yrs of stagnation. This is a weaker forcing than for the glaciers in the embayments, as well as for the medium and small depressions, and thus suggests that this bottleneck provides less stability than these geometries. This 265 contrasts with the common pattern, where a small S and a positive dS should stabilize the glacier strongly. It is unclear why this is the case here. We hypothesize that it might be related to a combination of high driving stresses due to a steepened surface inside the bottleneck in conjunction with high modeled calving rates (not shown). The two experiments with glaciers in geometries with narrower bottlenecks ('medium' and 'large' bottleneck) did not retreat through the entire perturbation.
This, in turn, aligns well with the general notion of a confined (low S) and downstream narrowing (positive dS) fjord yielding 270 strong stability to the glacier. Likewise, none of the glaciers in fjords with bumps retreated completely, which follows the same concept. However, the strong stability that bumps provide to the glacier may also be related to the choice of model parameters (sect. 4.2).

Stress balance response to fjord geometry
To assess the underlying mechanisms behind the geometric controls described before, we now analyze the stress regimes across 275 the studied geometries. We focus on lateral shear stress gradients for lateral perturbations, and longitudinal stress gradients for basal perturbations as given by the SSA in x-direction by where σ ′ is the deviatoric stress and τ b is basal drag. We interpret the second and third term on the right-hand side as longitudinal stress gradient and lateral shear stress gradient, respectively. With the imposed spatially uniform friction coefficient, 280 variations in the investigated stress fields are largely caused by variable fjord topography, and are hence convenient to investigate for our purpose.
For embayments and bottlenecks, variations in lateral shear stress gradients can be seen along the fjord walls (Fig. 5a,b).
Specifically, strongly negative shear stresses are found where ice is funneled in a downstream narrowing fjord. This occurs, for example, where 55 km < x <65 km in embayments (Fig. 5a) and at 45 km < x <55 km in bottlenecks (Fig. 5b). This indicates 285 enhanced resistance to flow for the glacier originating from the fjord walls, which provides stability to the glacier. Where ice flow diverges in a widening section of the fjord (in the upstream half of embayments and the downstream half of bottlenecks, Fig 5a,b), lateral shear stress gradients are comparatively weak. This indicates that the glacier -fjord wall contact is reduced here, and that the fjord walls provide little support to the glacier in these areas.
For depressions and bumps, we see clear variations in longitudinal stress gradient along the glacier bed (Fig. 5c,d). In 290 depressions, a band of negative values stretches across the full width of the outlet channel where the bed turns from being prograde to retrograde (at x = 55 km, Fig. 5c), indicating that ice flow is being blocked here. Likewise, a marked reduction in positive longitudinal stresses is seen at the onset of the bump where the bed slope switches sign (at x = 45 km, Fig. 5d).
Together, the stress regimes in basal perturbations demonstrate that a retrograde glacier bed, tilted against the direction of flow, reduces longitudinal stress gradients considerably as it increases the basal resistance to flow, which ultimately stabilizes the 295 glacier.
In summary, the stress analysis above suggests that increased lateral shear stress gradients or negative longitudinal stress gradients are found wherever ice flow is forced to converge, either horizontally or vertically, towards a narrowing or shoaling area downstream. Simulations using asymmetric as well as longer perturbations confirm that these findings are robust (see Fig. SA3). Through the convergent flow, the contact between the glacier and the fjord is enhanced, leading to increased resis-300 tance to flow. Overall, along-flow change in fjord width or depth (i.e., dS) is found to define areas of increased lateral shear gradients or negative longitudinal stress gradients, and thus GL stagnation.

A universal quantitative relationship for ice-topography interaction
We hypothesize that there is a quantitative relationship between fjord geometry and glacier retreat, valid across a range of different geometries. To test this, we correlate a variety of metrics indicative of glacier retreat (Table 3) against relevant 305 metrics of fjord geometry, that is, the submerged cross-sectional area (S) and its derivative (dS). We restrict the data to those instances when the GL is located within a geometric perturbation (gray-shaded areas in Fig. 3). Among all combinations of retreat and geometry metrics tested, including those of asymmetric and 'longer' perturbations (Table 2), the clearest and most universal relationship found is a negative, close to linear correlation between the ratio of the GL flux and the submerged crosssectional area Q GL /S over the change in submerged cross-sectional area dS (Fig. 6). This relation expresses that a widening 310 or deepening fjord in the downstream direction (negative dS) promotes a high GL flux per wetted area (Q GL /S). Conversely, a glacier retreating in a fjord that becomes narrower or shallower downstream (positive dS), will have a reduced Q GL /S. Note that the ratio Q GL /S for basal geometry perturbations (grey to black and green colors in Fig. 6) is on average lower than for lateral geometry perturbations. This means that basal perturbations generally inhibit ice flux across the GL more efficiently than lateral perturbations (note that this may be influenced by our modelling choices (sect. 4.2)). Also, note that the GL flux 315 is the product of the velocity v GL and the flux gate area at the GL A GL , that is Q GL = v GL × A GL . The ratio Q GL /S is thus proportional to v GL when there is hydrostatic equilibrium at the GL (because in that case, S = 0.9 × A GL ), and so we find a comparable, negative linear relationship between v GL and dS (Fig. SA4).
We find an additional, yet less distinct, negative relationship between the wetted area S and the GL retreat rate dGL (Fig. 7a).
This shows that a wider or deeper fjord promotes faster GL retreat, while a narrower or shallower geometry stabilizes the 320 glacier. This relation is not as universal as the previous one since one value for dGL is not uniquely linked to one value for S across different geometries. Furthermore, it is not linear, but rather such that for a range of low S values, dGL does not vary noticeably. Only above a certain threshold in S, the GL retreats markedly faster (Fig. 7b). This threshold varies between different fjord geometries. However, we find that it is always associated with the location of GL stagnation (Fig 3.1). This means that a relationship between GL retreat rate dGL and S only unfolds if a local stability position is passed. These stagnant 325 positions can be either where S is low, or where dS is high, as shown previously. For instance, dGL does not increase as the GL retreats very slowly at the stagnant position in the downstream half of embayments. Only once it has retreated passed this point of GL stagnation, a correlation between dGL and S can be seen.
For depressions we do not see a distinct relationship between dGL and S (Fig. SA5). This is because we measure the GL position x GL and therefore also dGL as the furthest downstream grounded point along the central flow-line of the glacier.

330
When the glacier ungrounds in the center of a depression, where the fjord is deepest, the dynamics of retreat are triggered several kilometers upstream of the GL, as mentioned in sect. 3.1. Therefore, there is a correlation between fjord depth and GL retreat in depressions. However, it is not reflected when only considering processes at the GL. Not finding a dGL over S relation for depressions is hence expected by construction of our methodology, and not an actual feature.

335
Given our previous results, we now aim to assess whether our principal geometric relationship Q GL /S over dS can be found for Jakobshavn Isbrae. To this end, we calculate the wetted area S along the topography of Jakobshavn Isfjord as used in Kajanto et al. (2020) which depicts overall higher values and larger along-fjord changes in dS than our idealized settings (Fig. 8b).
Plotting all available data points for Q GL /S over dS at Jakobshavn, we do not find the aforementioned geometric relation- ship. This may have many reasons related to the complex dynamics of Jakobshavn Isbrae (Bondzio et al., 2017), but most  , 2018). This alters the stress balance at the GL compared to our experiments where the glacier is always closely confined between fjord walls. Specifically, it means that the rigid glacier-wall interface in our experiments is replaced by a changing ice-ice contact. This has implications for the lateral friction that the fast-flowing ice in the main channel experiences and for the processes transferring stabilizing back stress from the sides to the center of flow. 345 We thus expect our findings to be more easily transferable to settings where Jakobshavn Isbrae is enclosed by fjord walls.
This is only the case in one part of the outlet channel, upstream of the present-day front (Fig. 8c). Indeed, values for Q GL /S are inversely related to dS in a qualitative way here, such that an increase in dS is generally associated with a decrease in Q GL /S and vice versa (Fig. 8d,e), consistent with our findings for synthetic geometries above (Fig. 6). Even though this relationship is only qualitative, meaning that one value for dS is not uniquely associated with one value for Q GL /S, we find these results 350 encouraging given the complexity of Jakobshavn Isbrae's dynamics. For settings resembling our setup more closely, such as medium-sized outlet glaciers found in, e.g. Greenland (Carr et al., 2014;Bunce et al., 2018;Catania et al., 2018), Svalbard (Schuler et al., 2020) and Novaya Zemlya (Hill et al., 2017), we expect an even stronger imprint of topography on retreat dynamics.

Mechanisms behind geometric controls of glacier dynamics
The current study offers new quantitative insights into how topography influences the evolution of marine outlet glaciers, and their response to ocean warming. We demonstrate that two topographic metrics, the wetted area S and its derivative dS, jointly control the dynamics and retreat of glaciers constrained by fjord walls. Together, these metrics largely explain variations in grounding line mass flux Q GL , which is important in the context of sea-level rise, and the grounding line retreat rate dGL.

360
Based on our stress analysis and physical principles, we propose the following physical interpretation for these results: First, a downstream narrowing or shoaling fjord (positive dS) stabilizes the glacier, as ice flow is funneled through the constriction enhancing the glacier-fjord contact (Fig. 5a,c). This increases the basal or lateral resistance to flow, which stabilizes the glacier.
Conversely, a downstream widening/deepening fjord (negative dS) provides little support to the glacier as glacier-fjord contact is reduced (Fig. 5b,d). Second, a narrow fjord (low S) stabilizes the glacier, because the distance between the lateral ice margins, 365 where friction with the fjord walls is high, and the center of flow is small. This means that the part of the glacier where ice flow is largely undisturbed is reduced (Raymond, 1996;van der Veen, 2013). Third, a shallow fjord (low S) stabilizes the glacier because the glacier is further away from flotation, and thus grounding line retreat is less likely to occur with a given amount of thinning (Pfeffer, 2007;. In our experiments, the area exposed to ocean melt does not have a large effect on retreat dynamics. Even high oceanic melt rates, which could compensate for a small ice-ocean interface, do not trigger 370 retreat through geometric perturbations where S is low. For a particular fjord geometry, the relative importance of S or dS in providing stability to the glacier may vary. This will be discussed with two examples from our results: In embayments, S is larger than the reference fjord. Therefore, if S was the dominant control for glacier dynamics here, retreat through embayments should occur more easily than through the reference fjord. However, we find that the opposite is true; the grounding line stagnates at the downstream end of embayments (Fig. 3a), 375 while it retreats steadily through the reference fjord (Fig. SA1). This indeed confirms that S alone does not explain glacier retreat. Rather, dS controls glacier dynamics, because the point of grounding line stagnation is where the fjord changes from wide to narrow in the direction of ice flow. For bed bumps in our experiments, the picture is different. Our model glaciers stagnate on or near the crest of the bumps, where dS is close to 0 or negative (Fig. 3d). This should not be an obstacle for retreat if dS was the dominant control on glacier dynamics. Therefore, it must be the shallowness of the fjord at this point (indicated 380 by low S) which governs the dynamics here.
Given these disparities between different settings, it is all the more compelling that we find the geometric relationship Q GL /S over dS universal to all our tested fjords. It implies that given the current grounding line mass flux Q GL and the upstream subglacial topography of a particular glacier, a well-founded estimate of the topographically induced future contribution to sea-level rise can be made. To the authors knowledge, this type of universal quantitative link between fjord topography 385 and glacier response has not been established before, going beyond the qualitative descriptions of ice-topography interaction offered in previous studies Carr et al., 2014;Bunce et al., 2018;Catania et al., 2018;Åkesson et al., 2018b). For projections of future sea-level rise, this direct coupling between topography and ice discharge is highly relevant, as it enables an ad-hoc assessment of the expected future sea-level contribution of a glacier on decadal to centennial time scales.
Q GL is readily available for glaciers where the velocity and bathymetry is well-known. However, the physical interpretation 390 of the relationship between Q GL /S and dS is not straight-forward. Since Q GL /S is proportional to v GL when there is hydrostatic equilibrium at the grounding line, the expression can be thought of as relating grounding line velocities to along-fjord changes in fjord topography, through the mechanisms described in our stress analysis (sec. 3.3). Accordingly, our results show that velocity evolution at the grounding line over time is also a good proxy for the dynamic response of a glacier to fjord topography (Fig. SA4). This may be specifically useful for less well-studied glaciers with unknown bathymetry. Notably, our 395 geometric relationship is distinct from a typical mass-conservation argument, which simply states that velocities must increase for a decreased flux gate, and vice versa, to maintain the same grounding line discharge. This is because in such an argument, velocities are related to absolute values of fjord width or depth, and not the along-fjord change of fjord geometry. While massconservation mechanisms certainly play a role in our simulations, we do not find that such a relationship alone is sufficient to fully explain the dynamics we observe.

400
Our second quantitative relationship between dGL and S confirms the widely accepted concept that a wide or deep fjord promotes fast grounding line retreat (e.g. Warren and Glasser, 1992;Carr et al., 2013;Bunce et al., 2018;Catania et al., 2018;Åkesson et al., 2018b). However, we highlight that this relationship may not hold if the fjord is narrowing or shoaling downstream (positive dS). In practical terms, this means that retreat from a uniformly wide and deep channel into an upstream widening or deepening section does not automatically imply that retreat has to accelerate. Rather, we suggest that 405 a glacier will sit at the downstream end of the section, where S is increasing upstream, for a considerable time, or may not even retreat further, because it is particularly stagnant here. Only if this position is abandoned, fast retreat will occur (Fig. 5a).
This fast retreat is in fact facilitated by the long residence time of the grounding line, since concurrent upstream thinning preconditions the glacier for fast retreat.
In contrast to most previous studies, we emphasize the role of along-fjord change in fjord topography (dS) to explain geo-410 metric controls of outlet glaciers. Along-fjord change in fjord depth is also key in the context of the marine ice sheet instability (MISI) theory, according to which retrograde beds promote retreat (Schoof, 2007;Gudmundsson et al., 2012;Gudmundsson, 2013). However, even though we do have retrograde beds in our fjords with basal perturbations, we do not see any influence of the MISI on glacier dynamics. This is simply because our tested glaciers never retreat into an area where the bed slope is strongly negative, and where the MISI effect would be expected to occur. On bumps, the glaciers stop to retreat on the 415 downstream side where the bed is prograde (Fig. 3d). In depressions, the grounding line stagnates where the bed slope is only slightly negative (Fig. 3c), which is not enough to trigger a MISI feedback loop. Retreat off this stagnant position occurs through ungrounding several kilometers upstream of the grounding line. This process is not related to typical MISI dynamics.
Besides that, our quantitative relation Q GL /S over dS may seem contradictory to widely accepted concepts of glacier dynamics, because we project high Q GL /S for prograde beds (i.e., negative dS in our study). This may give the impression that 420 prograde beds should lead to accelerating ice discharge. However, we emphasize that we assess ice discharge per area, not absolute values of ice discharge. A glacier retreating on a prograde bed will experience a reducing wetted area as it recedes, and thus the ratio Q GL /S may increase, but not the absolute grounding line flux. This is exemplified by our experiments with bumps, where the glacier stagnates on a prograde bed even though Q GL /S is relatively high (Fig. 6).
Only a few studies have considered the influence of along-fjord changes, rather than absolute values, in fjord width on 425 glacier dynamics, and available observations are limited in time (Carr et al., 2014;Bunce et al., 2018). The main consensus is that a fjord widening in the direction of glacier retreat promotes fast grounding line recession, while a narrowing fjord reduces retreat rates. Furthermore, retreat onto a pinning point can stabilize the grounding line. This is related to our findings in that we also see accelerating retreat the further the grounding line moves into a wider fjord upstream (c.f. grounding line positions in the downstream half of the embayment (55 km < x < 65 km) in Fig. 3a). However, in our results, this only occurs after a 430 phase of grounding line stagnation at the downstream end of such fjord sections. We do not find conclusive evidence in the observational record whether these points of grounding line stagnation are a relevant phenomenon in real-world settings or not (Carr et al., 2014;Bunce et al., 2018;Catania et al., 2018). Further research analysing a range of fjord geometries and glacier retreat histories is required to test this result. We do see some signs in our experiments that retreat slows down the further the grounding line recedes into a narrower fjord upstream. Overall though, retreat in upstream narrowing fjords is markedly faster 435 than if the fjord is upstream widening (compare grounding line positions in the downstream half of bottlenecks (55 km < x < 65 km) with the ones in the upstream half (45 km < x < 55 km) in Fig. 3b). This we explain with the aforementioned enhancement (reduction) in fjord-glacier contact for an upstream widening (upstream narrowing) fjord (Fig. 5a,b). Thus, we confirm that retreat slows down in an upstream narrowing fjord, but in the context of a retreat cycle through both upstream widening and narrowing fjord sections, overall faster retreat occurs through upstream narrowing fjords. Therefore, the observational records 440 of glacier retreat in Greenland may be too short and the fjord-width variations too small to testament similar dynamics as we observe in the model (Carr et al., 2014;Bunce et al., 2018;Catania et al., 2018). However, in line with the observational record, we can reproduce the stability that lateral pinning points offer. This is clearly demonstrated by the strong stability that narrow bottlenecks provide in our experiments.
Our experiments are set-up to simulate grounding line retreat, and hence we do not offer any insights on ice-topography 445 interactions for advancing glaciers. Previous studies, however, suggest that fjord geometry induces hysteresis in the retreatadvance cycle of a glacier, meaning that a reversal to colder conditions after a phase of climate warming does not allow the grounding line to advance to the same position it occupied initially if the fjord is widening or deepening in front of the glacier (Brinkerhoff et al., 2017;Åkesson et al., 2018b). We expect this also to hold for our experiments if we had simulated ocean cooling following the warming scenarios tested.

Study limitations
As in all numerical studies, our results have limitations related to the choice of model parameters. In particular, there are three aspects that warrant further discussion. First, the SSA used here is not a full representation of the stress regime in a glacier which especially bears relevance near the grounding line. For weak beds and fast-flowing outlet glaciers, as we aim to mimic in our synthetic setup, the SSA is a reasonable approximation and widely used in the glaciological literature. It may fail, however, 455 on steep bed slopes, such as the ones in our simulations of basal perturbations. In particular, the effect of basal high-friction points may be overestimated which might explain why basal perturbations generally yield lower Q GL /S values than lateral perturbations in our experiments (Fig. 6), and why bumps were found to be excessively stabilizing. However, we do not expect this effect to compromise our conclusions in general. For instance, Favier et al. (2012), using a setup similar to ours to simulate the effect of a basal pinning point on ice dynamics with a Full-Stokes model, present stress distribution patterns that agree 460 favorably with the ones shown here. The influence of using the SSA as compared to a Full-Stokes model on our results is therefore expected to be limited.
Second, choosing a Budd-type friction law which introduces an elevation-dependence of the basal friction through the effective pressure adds complexity to the interpretation of our results. Specifically, it means that bed bumps and depressions alter the basal resistance to flow purely through their elevation difference to the surroundings. This may be another reason commonly used ones, previous studies have also shown that different friction laws can lead to substantial differences in transient ice dynamics and steady-states (Brondex et al., 2017;Åkesson et al., 2021). A comparison between different friction laws is outside the scope of this study, and thus the impact of this modeling choice is difficult to estimate. Third, the calving parameterization chosen is an important control on the simulated dynamics. Other idealized studies have applied a calving 470 law with a prescribed ice front position or a prescribed ice shelf length (Schoof et al., 2017;Haseloff and Sergienko, 2018), which both have the disadvantage of being unknown in general, whereas we opt for the von-Mises parameterization due to its relatively good performance when applied to real-world glaciers (Choi et al., 2018). Again, in the absence of a universal calving law, the potential effect of this modeling choice on our results is hard to assess. Overall, it can be assumed that the geometric relationship found here is not as distinct in complex real-world settings. To what degree our results are transferable to a specific 475 glacier will depend on the degree to which the stress regime at the glacier front and the grounding line is comparable to our setup. This is what we demonstrated with the example of Jakobshavn Isbrae. Nonetheless, since we use a state-of-the-art model with parameterizations that have been calibrated against typical values for marine outlet glaciers, we expect our findings to be applicable to a wide range of settings.
Many of the medium-sized glacier catchments in Greenland may yield a long-term contribution to sea-level rise as thinning at 480 the outlet glaciers may propagate far upstream (Felikson et al., 2021). Here, we offer a quantitative perspective on the processes at the grounding line, and highlight the importance of assessing both the wetted area and the along-fjord change in wetted area in order to accurately describe ice-topography interactions. These two parameters together determine the geometrically induced ice discharge to the ocean, which is crucial for sea-level rise, and the expected future retreat of marine outlet glaciers.

485
The shape of a fjord can promote or inhibit glacier retreat in response to climate change. Here, we use a numerical model to study such ice-topography interactions in a synthetic setup under idealized conditions. We find that variable fjord topography induces gradients in lateral or basal shear stresses which then influence glacier dynamics. Increased shear at the ice-fjord interface, which stabilizes the model grounding line, is caused by converging flow towards a downstream constriction, because such flow enhances the glacier-fjord contact. Conversely, areas of reduced shear, which promote fast retreat, are found where ice 490 flow diverges, because glacier-fjord contact is reduced. In practical terms, this means that retreat of a glacier into an upstream widening or deepening fjord does not necessarily promote retreat, but may in fact stabilize the glacier. We also confirm that rapid retreat is more likely to occur through deep and wide fjords, while slower retreat is expected for narrow and shallow topographies.
Furthermore, using the concept of the wetted area, which is the submerged cross-sectional area of a fjord, and its along-fjord 495 change, we quantitatively link grounding line discharge and retreat rate to fjord topography. Specifically, we postulate that given the current grounding line flux and the upstream subglacial topography of a particular glacier, an ad-hoc estimate of the topographically induced component of future mass loss can be made. For less well-studied glaciers, we demonstrate that the velocity evolution over time is a promising proxy for the dynamic response of a glacier to fjord topography.
We expect that the quantitative relationships between topography and retreat dynamics are most likely to be transferable to 500 real-world glaciers confined by fjord walls, while being less relevant for ice streams where lateral ice flow influences dynamics considerably. We demonstrate this using the example of Jakobshavn Isbrae.
Future studies should aim to verify our findings using real-world observations. Particularly valuable in this context would be long-term observations of (sub-)annual grounding line positions, calving fronts and velocity changes, combined with detailed bathymetric maps for glaciers confined by fjord walls in Greenland and the Arctic.