Evolution of ice rises in the Fimbul Ice Shelf, Dronning Maud Land, over the last millennium

Abstract

We investigate two ice rises, Kupol Moskovskij and Kupol Ciolkovskogo, in the Fimbul Ice Shelf, East Antarctica, situated ~60 km from each other but differing in their glaciological settings. We apply a thermo-mechanically coupled Elmer/Ice model to profiles going across these ice rises and use it to investigate their past evolution covering present to several millennia ago. We constrain the model results using field measurements, including surface-velocity measurements, and surface mass balance estimated by isochronous radar stratigraphy dated with firn cores. We find that the ice rises are thickening at present (2012–2014), which started only in recent decades. Investigation of deeper radar reflectors suggests a stronger upwind-downwind contrast in surface mass balance in the past for both ice rises, with varying details. This result matches what was previously found on Blåskimen Island ice rise, which is also in the Fimbul Ice Shelf. Moreover, Kupol Moskovskij, situated at a shear margin, shows signs of recent changes in its ice-divide position, while Kupol Ciolkovskogo shows a more stable divide position. This study highlights the long-term influence of surface mass balance on ice rises, as well as the strong influence of local glaciological settings on their evolution.

Introduction

Ice rises, with their local ice flow regime, can act as local archives of the Antarctic coastal climate (e.g. Philippe et al. 2016, Winstrup et al. 2019) and ice dynamics (e.g. Conway et al. 1999, Wearing & Kingslake 2019). Ice rises can vary significantly in their shape and size as well as in their influence on the stability of the coast, which can vary greatly from a net stabilizing to net destabilizing influence (Matsuoka et al. 2015, Furst et al. 2016, Still et al. 2019, Goel et al. 2020). The information retrieved from these ice rises can also vary greatly from climatological to ice dynamical, covering recent decades to many millennia (Mulvaney et al. 2002, Kingslake et al. 2016).

Typically, the ice flow over an ice rise is slower than the ice-shelf flow around it by one to two orders of magnitude, allowing ice rises to store local glaciological information of millennial timescales within their ice (Matsuoka et al. 2015). Some of this stored information can be retrieved by examination of their englacial radar stratigraphy consisting of pseudo-parallel layers, which are widely accepted as isochrones (Kohler et al. 1997, Richardson et al. 1997, Fujita et al. 1999, Spikes et al. 2004). The shape of these layers can be influenced by surface mass balance (SMB), ice-flow pattern, bed topography, basal melting and ice rheology (Raymond, 1983, Parrenin et al. 2006, Sinisalo et al. 2013, Hudleston 2015). An ice-flow model can be used to interpret the englacial stratigraphy and constrain the past evolution of an ice rise and its surroundings. This method of elucidating the past glaciological information of an ice rise has been applied to several ice rises around Antarctica, such as Roosevelt Island (Conway et al. 1999), Berkner Island and Fletcher Promontory (Hindmarsh et al. 2011), Kealy Ice Rise (Martín et al. 2014), Halvfarryggen Ice Rrise (Drews et al. 2013), Derwael Ice Rise (Drews et al. 2015, Callens et al. 2016), Blåskimen Island (BI; Goel et al. 2018), Kroff Ice Rise (Kingslake et al. 2016) and Henry Ice Rise (Wearing & Kingslake 2019).

Dronning Maud Land (DML) has ice rises almost every 100 km along its coast. However, in an earlier study, Goel et al. (2020) examined the presence of surface lineations parallel to the topographical ridges of ice rises in DML and inferred that only a small fraction of ice rises seems stable over the last few millennia. Such lineations are surface manifestations of well-developed englacial features that require a stable ice-divide position to develop (Goodwin & Vaughan 1995, Martin et al. 2009a). Only 4 of 33 ice rises in DML show such lineations in satellite imagery, indicating that the majority of ice rises underwent changes in their divide position over the last few millennia. Investigations of the recent mass balance over the past few decades of five ice rises in DML have shown that changes in SMB represent the primary driver of changes in ice-rise mass balance (Goel et al. 2022). Yet, is this true over longer timescales? A study of the past evolution of BI ice rise in the Fimbul Ice Shelf, DML, showed that the ice rise remained in balance over recent millennia, whereas the SMB pattern over the ice rise evolved with time (Goel et al. 2018). Can the same be said for neighbouring ice rises with different glaciological settings?

Here, we investigate two ice rises, Kupol Moskovskij (KM) and Kupol Ciolkovskogo (KC), east of the Jutulstraumen Glacier outlet in the Fimbul Ice Shelf, DML (Fig. 1a). KM is situated at the calving front, whereas KC is entirely surrounded by the ice shelf. These ice rises, although only ~60 km away from each other, differ significantly in their settings, including their location relative to the calving front, grounding line, shear margins and morphology. KM (area: 508 km²; Moholdt & Matsuoka 2015) is an isle-type ice rise located at the calving front of the Fimbul Ice Shelf, with the Jutulstraumen shear zone to its west (Fig. 1a). It features a distinct ridge running southwards from its summit ~230 m above the surrounding ice shelf, ending in a tail. It is grounded ~200 m below sea level on a bed that slopes gently to the west and more steeply to the east (Fig. 2a,b). KC (area: 450 km²; Moholdt & Matsuoka 2015) is also an isle-type ice rise encircled by a considerably slower-moving ice shelf than KM. It is grounded on a shallower bed (~80 m below sea level) despite being similarly elevated above the ice shelf (Figs 1a,c & 2c,d).

Figure 1. Ice rises and field data in the Fimbul Ice Shelf, Dronning Maud Land. a. The ice rises are highlighted using blue outlines (Moholdt & Matsuoka 2015), and ice-flow speed in the region is depicted using the colour scale (Rignot et al. 2011). The background satellite image is a moderate-resolution imaging spectroradiometer (MODIS) mosaic (Haran et al. 2005) with hill shade from the Reference Elevation Model of Antarctica (REMA; Howat et al. 2019). The map's coverage is shown in the inset. Black boxes show the coverage of panels b. and c. Panels b. and c. show close-up views of Kupol Moskovskij (KM) and Kupol Ciolkovskogo (KC) ice rises, respectively. Contours at 20 m intervals are used to depict ice-surface topography. Radar profiles are indicated by green and yellow curves. The dots represent the locations of Global Navigation Satellite System (GNSS) stakes that were used to measure ice velocity, which is shown as black arrow markers and with the magnitude depicted as colour fill in the dots. Blåskimen Island (BI) was investigated by Goel et al. (2017, 2018), and those results are compared with the results from KM and KC in this paper. The figure was made using Quantarctica (Matsuoka et al. 2021).

Figure 2. Radargrams (2 MHz) across (a. & b.) Kupol Moskovskij and (c. & d.) Kupol Ciolkovskogo ice rises. Panels a. and c. show Profile A and panels b. and d. show Profile B of both ice rises. The locations of the radargrams are shown in Fig. 1. The data are shifted using Global Navigation Satellite System (GNSS) elevation measurements along the radargrams. The surface and bed are highlighted in blue. Tracked reflectors are shown in red. m a.s.l. = metres above sea level.

In this study, we adopt the method employed in Goel et al. (2018) to study BI nearby and simulate the englacial stratigraphy of these ice rises across multiple profiles to elucidate their evolution over the past several millennia. These ice rises were surveyed in the 2012–2013 and 2013–2014 seasons. All data used for modelling in the paper were already reported, so here we briefly describe the data and refer the reader to Goel et al. (2022) for a complete description.

Data

To measure the local ice flow and SMB over the ice rises, 3 m-long, hollow aluminium stakes were installed at 80 sites over KM and 56 sites over KC (Fig. 1b,c). These stakes were installed along the flow lines determined after preliminary analysis of the kinematic Global Navigation Satellite System (GNSS) data on site. Each stake's location was established using static GNSS measurements taken over 20 min intervals together with tape measurements of the stake's height above the snow surface. The measurements were repeated during the following year to determine each stake's lateral displacement and surface height difference for inferring SMB. Stake-derived flow speed across the ice divides on the ice rises shows near-stagnant ice (< 2 m a^-1) near the divide, with faster ice flow towards the ice-rise margin (Fig. 1b,c). The ice-flow pattern across the divide shows a 20% higher flow speed on the upwind side for KM (10 vs 12 m a^-1) and a more significant 40% difference on KC (5 vs 7 m a^-1) at ~7 km from the divide towards the flank (Fig. 3b,d).

Figure 3. Surface mass balance (SMB) and flow speeds along Profile A of Kupol Moskovskij (KM) and Kupol Ciolkovskogo (KC). Panels a. and c. show SMB used as the model input (Goel et al. 2022). The red curve depicts the SMB estimated using dated radar reflectors, whereas the green curve shows the SMB derived using stake height measurements. Panels b. and d. show observed flow speeds (markers), together with modelled flow speeds for different scenarios (dashed and solid curves).

This flow-speed contrast (and consequent ice-flow flux contrast) is consistent with the observed SMB contrasts and probably maintains the mass balance of the ice rises. Both ice rises exhibit higher SMB in the upwind slope than the downwind slope, as can be seen in Fig. 3a,c along profiles on both sides of the ice divide. This contrast is more substantial on KC. An anomalous drop in SMB observed at the ice divides is probably associated with wind erosion in this area, which has also been observed on other ice rises (Drews et al. 2015, Goel et al. 2017, Kausch et al. 2020). SMB estimates using all repeat stake height measurements show a mean SMB of 0.52 m of ice equivalent per year (m_ieq a^-1) for KM for the 2012–2013 period and 0.36 m_ieq a^-1 for KC for the 2012–2013 period.

SMB estimated using shallow radar reflectors dated using chronology developed by Vega et al. (2016) shows that the spatial pattern in SMB remained the same for both KM and KC for the last couple of decades. For KM, this method results in an average SMB of 0.92 m_ieq a^-1 for 16 years (1998–2014). For KC, the SMB is estimated to be ~0.36 m_ieq a^-1 for 30 years (1984–2013). For KC, the magnitude and spatial pattern of stake- and radar-derived SMB match well (Fig. 3c). However, in the case of KM, the magnitude of the radar-derived estimate is nearly twice the magnitude of the stake-derived estimate, whereas both estimates show similar spatial patterns (Fig. 3a). Goel et al. (2022) used these SMB estimates for KC and KM to estimate their mass balances using an input-output method and then compared them with the geodetic mass balances of these ice rises derived using satellite altimetry. They found that, for KM, the input-output estimate with stake-derived SMB (0.21–0.13 m_ieq a^-1 for 2013–2014) shows slight thickening, being close to the geodetic mass balance around this period (0.12–0.06 m_ieq a⁻¹ for 2012–2019 and 0.08–0.04 m_ieq a⁻¹ for 2011–2013). However, the input-output estimate with the radar-derived SMB suggests significantly stronger thickening (0.47–0.61 m_ieq a^-1 for 2004–2014) than the geodetic estimate. Goel et al. (2022) speculated that the firn-core chronology used for reflector dating is probably inaccurate. For KC, they found good agreement between geodetic mass balance estimates and the input-output mass balance, with both the stake-derived and radar-derived SMB showing slight thickening.

Borehole temperature measurements between 8 and 12 m depths at the ice-rise summits show nearly uniform firn temperatures of -15.9°C for KM and -17.1°C for KC. Firn cores from both ice rises do not show any major melt features (Vega et al. 2016).

Ice-flow modelling

Model setup

In this study, we model two-dimensional (2D) ice flow along profiles going across the ice rises on the steepest descent path (i.e. the direction of ice flow; Fig. 1b,c). The elevation contours along either side of the ice divide for both of the ice rises are nearly straight (Fig. 1b,c), affirming the 2D flow assumption used for modelling. We model ice flow across two near-parallel profiles across the ice rises, as this gives us more confidence in the validity of the model-derived evolution of the ice rises. We applied the same model setup that was used for the adjacent ice rise, BI (Goel et al. 2018). We use a thermo-mechanical full-Stokes ice-flow model Elmer/Ice model (Gagliardini et al. 2013) adapted for 2D ice flow over ice rises by Martín et al. (2009b). Here, we describe the model setup briefly and refer the reader to the appropriate papers for further details.

We assume 2D-plane flow along the profile with no basal sliding and a stress-free upper surface for the ice flow. Under the shallow-ice approximation (e.g. Van Der Veen 1999, section 4.2), we prescribe a depth profile of ice-flow speed at the downstream end of the numerical domain to match theoretical isotropic-laminar flow. The ice-flow flux via these boundaries is specified to simulate the desired mass balance scenario. These fluxes balance SMB over the entire modelled domain for a steady state and are further adjusted to mimic specific rates of thickening or thinning.

Over the modelled domain of both ice rises, the surface elevation decreases by ~20–25 m/km from the summit to the flank, with a total elevation difference of ~160–200 m. We assume the surface temperature over this surface to be homogeneous, with values based on borehole thermistor data at the ice-rise summits: -15.9°C for KM and -17.1°C for KC. We took geothermal flux from An et al. (2015): 55.3 mW m^-2 for KM and 55.7 mW m^-2 for KC. As the bed is assumed to be non-deformable, there is no heat advection through it. The domain boundaries downstream are assumed to have no horizontal temperature gradient.

The thermo-mechanical coupling is applied such that the flow-law parameter is updated at each time step according to a new ice temperature field (eqn (13) in Martin et al. 2009b, which is based on Dahl & Jensen 1989) determined by solving the heat equation that considers heat advection, diffusion and strain heating (see Martin et al. 2009b for details).

For rheology, we employed Glen's flow law with a rheological power index of n = 4.5 for all of our experiments except the sensitivity tests, which used a rheological power index of n = 3. The enhancement factor E was estimated such that its simulated ice-rise surface matched the observed surface well (Table I). E accounts for the various unknowns, such as grain size, impurities and ice fabric, that are not included in the model (Mangeney et al. 1996, Ma et al. 2010). We assume that the E does not evolve with time; thus, we assume that the various factors it accounts for are not evolving.

Table I. Input parameters used in the model.

KC = Kupol Ciolkovskogo; KM = Kupol Moskovskij.

Evaluation of results

For evaluating the model results, we use six prominent radar reflectors tracked along the profiles in the deep-sounding radargrams over KM and four reflectors along the profiles for KC (Fig. 2a–d; profile positions shown in Fig. 1b,c). Additional reflectors were tracked in the divide region, which were used to analyse the change in arch positions relative to the ice divide. Moreover, two additional reflectors from shallow-sounding radargrams were tracked to evaluate more recent events.

To evaluate how well the modelled isochrones match the tracked reflectors, we calculate the difference between them as D = H_r - H_i. Here, H_r and H_i are the (observed) reflector and (modelled) isochrone elevations, respectively. D is calculated at N points (spaced at 10 m) along the reflector where values of H _r are available. Furthermore, to evaluate the match over a reflector, we use mean absolute difference MAD by normalizing D as in Equation 1:

(1)$$MAD = \sum \displaystyle{{\vert D \vert } \over N}.$$

We use this as a measure of the goodness of match of our model simulations, and we try to achieve the lowest value. There is no specific target value of MAD whilst evaluating models run for different cases.

Although ice flow was modelled for the entire domain, we avoid interpreting modelled results within 2 km (four times greater than the local ice thickness) from the ice divide because the model is not fully adequate to evaluate this region where crystal alignments of ice play more significant roles (Pettit et al. 2007, Martin et al. 2009a). We make an exception to this whilst evaluating shallow reflectors, as the effect of divide flow near the surface is minimal.

Additionally, we exclude a region of 2 km from the downstream boundaries. This is because the error caused by the prescribing of the vertical boundary conditions there can propagate into the model domain up to a few times greater than the ice thickness (Hvidberg 1996, Martín et al. 2009a). To summarize, we evaluated the model results between 2 and 6 km from the summit on both sides.

Results

Past ice-divide positions

Clear, prominent Raymond arches are visible in both of the radar profiles (Fig. 2a,b; Raymond 1983). These arches are formed as a result of opposing driving stresses at the ice divide and the non-linear nature of ice deformation, resulting in stiff ice close to the bed. However, the arches at shallower depths (Fig. 3e,f; Goel et al. 2022) are not caused by accumulated anomalous stresses but by anomalous SMB near the divide (Vaughan et al. 1999). Regardless of their causes, these arches have been used as indications of long-term stability of the summit position and ice-thinning rates (Conway et al. 1999, Drews et al. 2015, Kingslake et al. 2016). These arches' presence also depends on other divide-flow characteristics such as stable summit position, minimal along-ridge flow and absence of basal sliding and melting. The absence of more developed ‘double arches’ and their impression on the ice-rise surface as ‘double lineations’ suggest that the ice divide position has been stable for one to four characteristic timescales (= ice thickness/SMB) but not longer (Martin et al. 2009a). KM and KC show a similar characteristic timescale of ~900 years.

We noted the highest position of these arches relative to KM's summit (Fig. 4) and found both profiles to show an eastwards movement of the arch positions with depth. The shallower arches show a more rapid change in the arch position of ~400 m over ~20 m depth. Deeper arches show a more gradual change in arch positions, with Profile B showing a change of ~400 m over ~140 m of depth and Profile A showing ~150 m over ~160 m of depth. While the two profiles agree on the more recent and faster changes in the shallower arches, they differ on the rate of change at deeper depths, suggesting different driving mechanisms.

Figure 4. Position of maxima of arches for different profiles relative to the ice divides of ice rises in the Fimbul Ice Shelf. Results for Blåskimen Island (BI) are derived from radar data from Goel et al. (2017). KC = Kupol Ciolkovskogo; KM = Kupol Moskovskij; Prof = Profile.

The stratigraphy of KC also shows distinct arches beneath the ice divide. Unlike KM, these arches are located closer to the ice divide, with gradual westwards changes in arch positions of ~200 and ~100 m over depths of 300 m depth. The absence of double lineations on the surface (Goel et al. 2020) and the absence of double arches near the bed suggest that the divide has been stable for a similar period of one to four characteristic timescales. These profiles also clearly show the effect of bed topography on reflectors, with reflectors draping over the bump in the bed (3–4 km east).

Recent mass balance

To determine the recent mass balance of the ice rises, we start with the case of steady state. The model is initiated using the current surface topography, and the ice surface is allowed to evolve freely whilst the ice discharge is maintained to balance SMB. The model is run for 25 characteristic timescales, which is ~25 × 900 years for KM and KC. This ensures that our model reaches steady state and that the isochronous pattern does not change in subsequent iterations. For these simulations, we used KM Profile A and KC Profile A (Figs 1b,c & 2a,c). We used both radar-derived and stake-derived SMB along these profiles as inputs to the model.

The steady-state simulations for KM resulted in the basal temperatures being considerably below freezing, with the maximum temperature below the divide being ~-10°C, suggesting neither basal melting nor sliding and affirming the model assumption. For the simulation with radar-derived SMB as the input, we see that the modelled flow speeds are thrice as fast as observed on both flanks (Fig. 3b), indicating that the ice rise is significantly out of a steady state. The modelled flow speeds for the run with stake-derived SMB also do not match the observations, but the difference from the observed flow speeds is much smaller. As the flow speeds from the steady-state simulation are faster than observations, a thickening scenario seems more probable considering mass conservation.

Thus, going with the more probable hypotheses, we next simulate a thickening ice rise. For this experiment, the elevation of the ice rise was prescribed to be lower by dH everywhere (i.e. ice is uniformly thinner than the present-day value by dH). dH is calculated such that the ice rise reaches the present ice surface after the prescribed thickening. After a 25 000 year model run, the model reaches a steady state with the lower ice surface, which serves as the starting point for a thickening experiment. To simulate thickening, the downstream outgoing ice flux is adjusted to have a net positive balance.

With radar-derived SMB as the input and a thickening rate of 0.5 m a^-1, we find that the surface flow speeds are still twice the observations (Fig. 3b). To reproduce the flow speed using the radar-derived SMB, thickening rates of ~1 m a^-1 or higher are required, which are very different from the satellite-based mass balance observations (0.12–0.06 m_ieq a^-1 thickening for 2012–2019 and 0.08–0.04 m_ieq a^-1 thickening for 2011–2013; Goel et al. 2022). Using the stake-derived SMB, we find that the matches between the modelled flow speed and observations are asymmetric. So next, as has been successfully modelled for the adjacent BI (Goel et al. 2018), we use two different thickening rates for each side of the summit. We find that rates of 0.25 m_ieq a⁻¹ on the eastern side and 0.15 m_ieq a⁻¹ on the western side result in the best match with the observed flow speeds, and these thickening rates are much closer to the satellite observations for recent years.

Both radar- and stake-derived SMB estimates show nearly identical spatial distributions. As radar-derived SMB seems to produce results that significantly differ from other field- and satellite-based observations, we only use the stake-derived SMB for KM, achieving closer results. We speculate that the firn-core chronology used for radar-derived SMB is not accurate.

For KC, unlike KM, the radar-derived and stake-derived SMB estimates are very similar in magnitude and spatial pattern (Fig. 3c). Moreover, for KC, the mass balance estimates obtained using these estimates agree with each other and the satellite altimeter measurements (Goel et al. 2022). As radar-derived SMB has higher spatial resolution and is representative of a longer period, we use the radar-derived SMB for the KC model runs. From the steady-state run, we find the bed to be frozen at -10.6°C, and the surface flow speeds are ~50% higher than the observations (Fig. 3d). Thus, KC is also not in a steady state and is probably thickening. Next, similar to KM, we simulated uniform thickening at different rates of 0.1 and 0.2 m_ieq a⁻¹ and found the matches to be asymmetric. Eventually, with the differential thickening rates of 0.12 m_ieq a⁻¹ on the south-eastern side and 0.08 m_ieq a⁻¹ on the north-western side, the model reproduced the observed flow speeds the best.

Onset of ongoing thickening

To determine the onset of the current thickening of KM, we ran the model at the determined thickening rate for different periods of 20, 40, 60, 100 and 200 years. We then tried to match the modelled isochrones to two reflectors tracked on shallow-sounding radargram at depths of 8 and 21 m near the summit (top two reflectors from the surface in Fig. 2a). We found MAD to be < 3 m between the isochrone and these two shallow reflectors (Fig. 5a,c). For the second shallow reflector from the surface, the best match (< 2 m) is observed at the thickening of 60 years and then worsens. For the first shallow reflector from the surface, however, the match only worsens with increasing thickening period. Both shallow reflectors show the worst match ~3–6 km from the divide.

Figure 5. Goodness of match between radar reflectors and modelled isochrones for steady-state and recent thickening scenarios. Panels a. and b. show matches between modelled isochrones (orange) and observed reflectors (blue) for Kupol Moskovskij (KM) Profile A and Kupol Ciolkovskogo (KC) Profile A. Panel c. shows the goodness of match measured as the mean absolute difference (MAD, Equation 1) for a range of thickening periods.

Similarly, for KC at steady state, the two shallow reflectors at depths of 6 and 18 m near the divide (Fig. 2c) showed MAD to be < 2.5 m, with a larger mismatch at -3 to -6 km from the divide (Fig. 5b,c). This match does not improve with increasing thickening period.

Although KM shows slight improvement with thickening over recent decades, KC shows no improvement. Thus, it can be said that the observed thickening has not sustained more than a few decades for both ice rises, with thickening over KC being more recent. Moreover, it is also possible that the ice rises have gone through several phases of thickening and thinning over this period, with insignificant net change.

SMB changes over the last millennium

To evaluate changes farther back in time, we now evaluate the deeper reflectors tracked from deep-sounding radargrams. We modelled two near-parallel radar profiles across the ridge for both ice rises. These two profiles were used to define flow bands for mass balance estimates using the input-output method (Goel et al. 2022). Our radar-mapped SMB patterns and GNSS-measured ice velocities are very similar between the two profiles (Figs 1b,c, 6 & 7). Bed topography is also similar between the profiles (Fig. 2a–d). These similar settings give us an opportunity to validate independent 2D models with each other and to increase confidence in the ice-rise evolution inferred with our modelling work.

Figure 6. Comparison of modelled isochrones to radar reflectors over two profiles of Kupol Moskovskij. In panels a. and d., the black curves show stake-derived surface mass balance (SMB) used as an input to the model for the two profiles, A and B. Red curves show potential SMB in the past (ΔSMB_West = -10%). Panels b. and e. show observed radar reflectors (blue) and modelled isochrones for the steady-state scenario (orange). Black curves show the surface and bed. Red labels show the ages of the modelled isochrones in years BP (estimated using SMB in panels a. and d.). Red labels show the ages of the modelled isochrones in years BP (estimated using SMB in panels a. and d.). Panels e. and f., are the same as b. and e. with differential SMB of ΔSMB_West = -10% and ΔSMB_East = 0%. m a.s.l. = metres above sea level.

Figure 7. Comparison of modelled isochrones to radar reflectors over two profiles of Kupol Ciolkovskogo. In panels a. and d., black curves show radar-derived surface mass balance (SMB) used as an input to the model for the two profiles, A and B. Red curves show potential SMB in the past (ΔSMB_West = -10%). Panels b. and e. show radar reflectors (blue) and modelled isochrones (red) for the steady-state scenario. Black curves show surface and bed topography. Red labels show the ages of the modelled isochrones in years BP (estimated using SMB in panels a. and d.). Panels e. and e. are the same as b. and e. with differential SMB of ΔSMB_West = -10% and ΔSMB_East = +40%. m a.s.l. = metres above sea level.

Because the firn cores used for this study are only tens of metres long, covering the past few decades, we do not know the exact age of the deeper reflectors used for this work. To estimate the ages of these reflectors, we associate them with the ages of the modelled isochrones that match them the best (Figs 6b,e & 7b,e). We use the present SMB values to date these reflectors.

We begin with the simulation of KM Profile A. As discussed earlier, the stake-derived SMB estimate is the most representative SMB estimate for KM that we have, and we use it for our simulations. As the thickening is very recent and its effect on deeper isochrones is minimal, we ignore it for simplicity and assume a steady state when deeper reflectors are interpreted. For steady-state simulation, isochrones on the eastern side match very well with the observed reflectors, with MAD = 2–4 m (Fig. 6b), whereas the match on the western side is somewhat variable: the top three reflectors match very well with MAD = 2–4 m, while the bottom three match to a lesser degree (MAD = 5–9 m). Two possible scenarios could explain the differences in match across the divide: 1) differential thinning and 2) differential SMB changes across the ice divide. Both scenarios, with thinning on the western side or lower SMB on the western side, would result in shallower isochrones on the western side and thus match better with the reflectors. However, we can rule out option 1) since a sustained thinning on one side would result in a divide migration of scale, which is not visible on KM. Moreover, such a scenario would result in an eastwards tilt in the arches, which is the opposite of the observations. We explore the likelihood of option 2) by trying out different input SMB scenarios and determining which comes closest to the observations.

Thus, we ran a steady-state experiment with shifted SMB values: ΔSMB_west = (SMB_past - SMB_present)/SMB_present = -10% on the western side. The results showed improved matches for the bottom three reflectors, with MAD = 3–7 m (Fig. 6c), and an increased mismatch for the top three reflectors (MAD = 4–8 m). The model, however, could not reproduce the curved shape of the reflectors on the western side, which is why MAD < 3 is not achieved.

The steady-state run for Profile B of KM shows similar features to Profile A. Overall, the modelled isochrones match very well with the radar reflectors on the eastern side (MAD = 1–2 m) but show significant mismatch on the western side (MAD = 5–14 m; Fig. 6e). We ran a steady-state experiment with ΔSMB_west = -10% on the western side, as we did for Profile A. We find that the top three isochrones match very well (MAD = 0–1 m; Fig. 6f). The fourth and fifth reflectors show almost no changes, whereas the deepest reflector shows a near-perfect match, with MAD = ~0 m.

For KC Profile A, with the radar-derived SMB as the input, we ran a steady-state simulation. We find that, unlike KM, there is a significant mismatch on both sides of the divide on this rise, except for the small region within 3 km east of the divide (Fig. 7b). This region corresponds well with the less steep region next to the divide due to the presence of a bump in the bed. We find the mismatch on the eastern side to be approximately the same for all of the reflectors (MAD = 8–11 m), whereas, on the western side, the amount of mismatch varies more (MAD = 6–15 m). Considering the variation in mismatch along the profile, we prescribe input SMB with ΔSMB_west = -10% for the western side and ΔSMB_east = +40% between 3 and 8 km on the eastern side. The results show a very good match for the reflectors on the eastern side, with MAD = 2–4 m (Fig. 7c). On the western side, the top three reflectors show a very good match (MAD = 1–3 m), whereas the deepest reflector does not match well (MAD = 13 m). We ran another experiment with ΔSMB_west = -30% and ΔSMB_east = +40%. The results show a very good match for the deepest reflector (MAD = 2 m), whereas the mismatch for the shallower reflectors increases (MAD > 15 m).

Profile B for KC shows a very similar isochrone pattern to Profile A. The mismatches on both sides are very similar to what is observed for Profile A. On Profile B, the bed bump is located farther to the east, ~4 km from the divide. For the steady-state run, we see that the region with the good match also extends up to the bump from the divide region (Fig. 7e). The reflectors east of this region show mismatches of MAD = 7–10 m, whereas for the region west of the divide MAD = 4–8 m. We attempted another experiment: ΔSMB_west = -10% and ΔSMB_east = +40%, which gave the best result for Profile A. This run results in improved matches on either side of the divide. On the east, we find very good matches for the middle two reflectors (MAD < 2 m; Fig. 7f). The MAD does not change much for the top reflector, although D changes its sign, whereas the match shows improvement for the bottom reflector (MAD = 7 m; improved from 10 m for the steady-state case). On the western side, all reflectors show good matches (MAD < 3 m), with a near-perfect match (MAD = ~0 m) for the third reflector from the top.

Overall, for KC, the steady-state run shows mismatches on both sides, excluding the region between the ice divide and the bed bump to the east. The mismatch east of the bump is improved for the scenario ΔSMB_east = +40%, whereas on the western side different scenarios between ΔSMB_west = -10% and -40% result in improved matches for the different reflectors. Using the present SMB to date these reflectors, we infer that, on the eastern side, SMB was higher than present by 40% between 2100 and 400 BP. Since 400 BP, the SMB on the eastern side has decreased to the present values. On the western side, the SMB was lower by up to 40% from the present SMB at 2100 BP and has been increasing since. At 1100 BP, it had increased by 30% and has since increased by the remaining 10% to its present value. Profile B shows similar changes with slightly different ages.

Model sensitivity

Due to the nature of the ice flow at ice divides, the dominant effect of non-linear rheology is limited to the divide region of ice rises (Pettit et al. 2007, Martin et al. 2009a). To see how sensitive our results are to the n, we re-simulated one of our primary results with n = 3 and compared the results. For KM, we simulated Profile A with ΔSMB_west = -10%. Although we find significant differences in the divide region, with smaller arches, the isochrones outside remain nearly the same. For the top three reflectors, we find a similar match, with MAD = 4–8 m (MAD = 4–8 m with n = 4.5). For the bottom three reflectors, MAD = 3–6 m (MAD = 3–7 m with n = 4.5). Essentially, the MAD values did not change for the choice of n; however, since lower n values result in softer ice, we find slight changes (5, 10, 15, 20, 30 and 45 years older) in the ages of the isochrones producing these matches.

Similarly, for KC, we re-simulated the Profile A steady-state run with ΔSMB_west = -10% and ΔSMB_east = +40% and n = 3. On the western side, we find similar matches of the top three reflectors, with MAD = 1–3 m (MAD =1–3 m for n = 4.5) and MAD = 13 m (MAD = 13 m for n = 4.5) for the deepest reflector, whereas on the eastern side, we found MAD = 1–4 m (MAD = 2–4 m for n = 4.5) for all four reflectors. The resulting changes in the ages of the four isochrones producing these matches are 15, 60, 85 and 150 years older.

Discussion

Uncertainty

Our model assumes 2D plane-strain flow, which can be inaccurate in the presence of a more complex geometry. An along-ridge flow on an ice rise can reduce near-bed viscosity and thus result in arches being smaller than the case without along-ridge flow or result in differences in isochrone depths (Martin et al. 2009b). Both KM and KC are ridge-shaped, and the profiles were made to follow the steepest decent paths on either side of the ridges. Being close to the summit, the along-ridge slopes between the two profiles for both of the ice rises is < 0.4°. Furthermore, the ratio of flow components orthogonal to a profile is < 1–2% for all four profiles based on ice-flow measurement. As a result, we do not expect any significant effects of either along-ridge flow or flow divergence on our results. Moreover, our comparison of model results along the two profiles, A and B, shows nearly identical results for both KM and KC. This suggests an absence of a significant along-ridge flow perpendicular to our profiles. The slight differences in the results are more likely to be attributed to differences in individual profiles, such as bed topography and input SMB.

For KM, the stake-derived single-year SMB and radar-derived 16 year-average SMB show similar patterns (Fig. 3a). This suggests that the spatial distribution of SMB has remained consistent over the 16 year period. For KC, we use radar-derived SMB as the input, which is representative of 30 years. Our long-term simulations suggest that the SMB distribution has probably changed over the last millennium in such a way that the upwind-downwind contrast was higher than at present. Our method, however, cannot determine the changes in the magnitude of SMB if the spatial distribution stays the same. This is because we have no ice-core chronology at these depths, and thus uniformly higher or lower SMB could give older and younger ages of ice at a given depth.

SMB changes: possible regional hypotheses

To discuss the regional evolution of the Fimbul Ice Shelf, in this section we evaluate our results together with the prior findings from BI (Goel et al. 2018).

The comparison of flow-speed measurements with the modelled flow speeds indicates that both KM and KC are thickening at present. However, further thickening simulations and comparison with shallow reflectors suggest that this thickening has not been sustained for more than a few decades. This is consistent with the prior findings on BI. Moreover, recent mass balance investigations over the three ice rises in the Fimbul Ice Shelf show that the mass balance trends over the last two to three decades are closely related to the SMB trends (Goel et al. 2022). However, we cannot determine whether this relationship has been persistent over longer periods.

While comparing deeper reflectors with modelled isochrones, we observe that, unlike KC, BI and KM produce a good match with the assumption of steady-state conditions. This match is further improved by simulating increasing SMB on the western side by ~10% over the last millennium. However, for KC, a good match for the deep reflectors requires SMB changes on both sides of the ice divide. The SMB changes are significant (~30–40%) on both sides; however, on the south-eastern side, the change happened in large steps at 2100 and 450 BP, while on the north-western side, the changes happen more gradually since 2100 BP, similar to BI and KM.

All of these differences in past histories on these sites highlight the importance of the local glaciological conditions in the interpretation of the radar layers. For example, BI and KM are near fast-flowing ice that potentially drives the mass balance of the ice shelf, whereas KM has its main ridge perpendicular to the direction of the prevailing wind direction. The fact that our inference of past changes in SMB agrees on these three sites with such distinctive glaciological settings suggests that these findings are probably applicable over the larger region of the Fimbul Ice Shelf.

Arch positions and ice dynamics

Although we are not evaluating the model results in the divide region to determine its evolutional history, it is still interesting to look at the model results in this region. Similar to the prior section, we evaluate our results together with the prior findings from BI (Goel et al. 2018).

Both of the profiles of KM and KC and the across-ridge profile from BI (Goel et al. 2017) show a tilt in the same direction over a similar period (Fig. 4), suggesting a large-scale forcing acting on all of the ice rises in the Fimbul region. In addition, both Profiles A and B of KM show that the arches shifted east from the ice divide close to the surface. Unlike a gradual tilt in arch positions seen over several ice rises, including KC and BI, this large shift of the arch position over a short depth range suggests a recent and rather rapid change in the ice-divide position. In our simulations, this feature is visible as a lateral offset between the modelled isochrones and the observed reflectors, resulting in overall worsening of the matches. We speculate, that this shift is probably a result of a localized forcing at KM's dynamic western shear margin with the Jutulstraumen Glacier outlet (Fig. 1a).

Overall, with n = 4.5 (as compared to n = 3), our model is able to better reproduce the amplitude of the arches in the divide region but not the width or other details. This value was inferred by Gillet-Chaulet et al. (2011) by looking at changes in along-flow englacial vertical velocity at the divide arches. Our results reinforce the point that using n = 4.5 is a good choice to simulate the divide region using an isotropic model. However, in agreement with Martin et al. (2009a), we find that we are unable to reproduce other details of the radar layers near the divide with our isotropic model, particularly their width or how they evolve under divide migration.

Role of bed topography

The bed underneath the arches is smooth in all four profiles from KM and KC (Fig. 2a–d). This suggests that the divide position and the position of the arches underneath are not controlled by a feature in the bed. However, we do observe the effect of bed topography on the englacial stratigraphy with the draping of reflectors on the bed bump 2–4 km east of the ice divide in KC (Fig. 2c,d), and this reflector pattern appears at the bump location all the way to the surface. As the bump is farther away and larger in Profile B, we see this effect move and become stronger. In this way, this feature in the bed topography also seems to determine the size of the low-sloping region (from the ice divide to the bump) on the eastern side of KC, and in this region we did not observe changes in SMB contrast with time. This further suggests the strong role of local topography in determining isochrone shape.

Conclusion

Using ice-flow models to infer changes in the cryosphere from the stratigraphy of an ice rise is becoming more common. The main difficulty is that stratigraphy is the combined result of all of the changes in the cryosphere (mainly changes in accumulation, geometry or flow in the surrounding area) and local ice flow around the ice divide. The latter is notoriously difficult to model, with current theory suggesting the important role in this played by the evolution of the crystal orientation fabric.

Here, we use a 2D plane-flow isotropic thermal model, but we only optimize our model with data in the region outside the low-strain divide-flow region. We employ this model on profiles across two ice rises in the same climatic regime but differing in their glaciological settings. This guides us in understanding the role of glaciological settings on the ice-rise evolution.

Our main findings on the KM and KC ice rises agree with prior findings from the adjacent BI ice rise. First, we observe current thickening on both of these ice rises, but we conclude that this is a recent occurrence and has not been sustained for longer than a few decades. Second, we infer a change in SMB pattern over the last millennium, with higher accumulation in the upwind slope and lower accumulation in the downwind slope. We note that these changes in SMB on the ice rises over the last millennium could not have been driven by changes in the surrounding ice-shelf flow, as there is no record of eastwards migration of the divide positions over this timescale. A more probable cause would be the variability in the factors driving SMB changes in this region, including katabatic winds, atmospheric westerlies and sea-ice extent. Third, although the observed changes in SMB and arch positions broadly agree, there are differences, probably resulting from the influence of the local glaciological conditions of the ice rises. We find that, unlike the reported changes in the mass balance of the ice rises in the Fimbul region in recent decades, these ice rises have been at steady state over the last several millennia. Nonetheless, their evolution over this period has been primarily influenced by changes in SMB, with the timing and the spatial distribution of these changes being strongly controlled by the local topography.