On the evolution of an ice shelf melt channel at the base of Filchner Ice Shelf, from observations and viscoelastic modeling

, Abstract. Ice shelves play a key role in the stability of the Antarctic Ice Sheet due to their buttressing effect. A loss of buttressing as a result of increased basal melting or ice shelf disintegration will lead to increased ice discharge. Some ice shelves exhibit channels at the base that are not yet fully understood. In this study, we present in-situ melt rates of a channel which is up to 330m high and located at the southern Filchner Ice Shelf. Maximum observed melt rates are 2 . 3ma − 1 . Melt rates decline inside the channel along ﬂow and turn into freezing 55km downstream of the grounding line. While closer to 5 the grounding line melt rates are higher within the channel than outside, this reverses further downstream. Comparing the evolution of this channel under present-day climate conditions over 250 years with its present geometry reveals a mismatch. This mismatch indicates melt rates two times higher were necessary over the past 250 years to form today’s channel geometry. In contrast, forcing the model with present-day melt rates results in a closure of the channel, which contradicts observations. Time series of melt rate measurements show strong tidally-induced variability in vertical strain-rates. Author contributions. AH has designed the study, conducted the ﬁeld study together with DS and wrote the manuscript together with OZ, JC, KWN, VH and LH. OZ processed the pRES/ApRES data and analyzed the melt rates together with AH. JC performed the viscoelastic simulations together with TS and with contributions from OZ and AH. JC, AH and OZ analyzed the results together with RM. VH and OZ processed the GPS data. LH and MW performed the inverse modeling. CH performed the seismic measurements and supported the discussions. NN calculated the TDX-DEM. KWN supported the ﬁeld study and contributed to melt rate analysis and its discussion together 525 with HFJC. All authors helped to improve writing.


Introduction
Melt channels in ice shelves have been hypothesized to destabilize ice shelves and were often linked to enhanced basal melt.
This triggered a variety of observational studies (Le Brocq et al., 2013;Langley et al., 2014;Drews, 2015;Marsh et al., 2016; models, which we will also engage below. In this study, we present in-situ melt rates of a large melt channel feature of the southern Filchner Ice Shelf at the inflow from Support Force Glacier (SFG). Field measurement and satellite-borne data enable us to investigate how this feature evolves using numerical modeling. In addition to the spatial distribution of basal melt, we are analyzing the temporal evolution of melt rates. We split this manuscript into two main parts, starting with observations followed by a modeling section. We present the 70 methodology and the results in each part separately. A synthesis is then focusing on the evolution of the melt channel.

Data acquisition
We acquired data at a melt channel on the southern Filchner Ice Shelf under the framework of the Filchner Ice Shelf Project (FISP). We performed 44 phase-sensitive radar (pRES) measurements (locations are shown in Fig. 1) in the season 2015/16, 75 that have been repeated in 2016/17 as Lagrangian-type measurements. These measurements were taken in 13 cross-sections ranging from 14 to 61 km downstream the grounding line (Fig. 1). This allows us to investigate the spatial variability of basal melt rates. At each cross-section, up to four measurements were performed at different locations: at the steepest western flank (SW), at the lowest surface elevation (L), at the steepest eastern flank (SE) and outside east of the channel (OE; Fig. 1b). In order to achieve an all-year time series, one autonomous pRES (ApRES) station was installed (Fig. 1b). This instrument performed 80 autonomous measurements every two hours resulting in 4342 measurements between 10 January 2017 and 6 January 2018. A GPS station was also in operation at this point from December 24, 2015 to May 5, 2016. To distinguish the single-repeated measurements from the autonomous measurements, we refer to them as pRES and ApRES measurements.  (Morlighem, 2020;Morlighem et al., 2020)). The study area near the Support Force Glacier (SFG) is marked with a black box. (b) Study area with pRES-derived basal melt rates at 13 cross-sections of the melt channel. The different symbols indicate the position relative to the channel, as shown in (c). For each cross-section, the distance from the grounding line and the duration of ice flow from the location furthest upstream are given. The location of an ApRES/GPS station is shown by a star. The seismic I, IV and V lines mark the location of active seismic profiles (Hofstede et al., 2021a, b). The background is a hillshade of the Reference Elevation Model of Antarctica (Howat et al., 2018(Howat et al., , 2019 overlaid by the ice flow velocity (Hofstede et al., 2021a). (c) Sketch of a cross-section of the channel with measurement locations on the steepest western surface flank (SW), at the lowest surface elevation (L), on the steepest eastern surface flank (SE) and outside east of the channel (OE).

pRES device and processing 85
The pRES device is a low-power, ground-based radar that allows for estimating displacement of layers from repeated measurements with a precision of millimeters (Brennan et al., 2014). This accuracy enables investigating even small basal melt rates, taking snow accumulation together with firn compaction and strain in vertical direction into account (Corr et al., 2002;Jenkins et al., 2006). The pRES is a frequency modulated continuous wave (FMCW) radar that transmits a sweep, called chirp, over a period of one second with a center frequency of 300 MHz and bandwidth of 200 MHz (Nicholls et al., 2015). For a better signal to noise ratio, the single-repeated measurements were performed with 100 chirps per measurement and the measurements of the time series with 20 chirps due to memory and power limitations. After collecting the data, each chirp was correlated with every other chirp in order to reject those which had a low correlation coefficient on average. The remaining chirps were stacked.
We followed Brennan et al. (2014) and Stewart et al. (2019) for data processing to get amplitude-and phase-depth profiles.
The final profile that contains the amplitude and phase information as a function of two-way travel time was received from 95 a Fourier transformation. To convert two-way travel time into depth, the propagation velocity of the radar wave is computed following Kovacs et al. (1995). For this the density is required. Here we use a model described by Herron and Langway (1980). As input parameters, accumulation rate and mean annual temperature is needed, for which we use data from the regional climate model RACMO 2.3/ANT (van Wessem et al., 2014(van Wessem et al., , multi-annual mean 1979(van Wessem et al., -2011. Despite the correction of higher propagation velocities in the firn, the uncertainty of the velocity and thus of the depth is 1% (Fujita et al., 2000).

Basal melt rates from repeated pRES measurements
The method for determining basal melting rates, previously described by e.g. Nicholls et al. (2015) and Stewart et al. (2019), is based on the ice thickness evolution equation. The change in ice thickness over time ∂H/∂t consists of a component arising from deformation and accumulation/ablation at both interfaces (e.g. . As our observations are discrete in time, the change of ice shelf thickness ∆H within the time period ∆t, that is caused by changes at the surface and 105 in the firn ∆H s (e.g. snow accumulation/ablation and firn compaction), by strain in vertical direction ∆H ε and by thickness changes due to basal melt ∆H b is considered: (Vaňková et al., 2020;. In order to obtain the basal melt rate, the change in ice thickness must be adjusted for the other contributions. Snow accumulation/ablation, firn compaction but also changes in radar hardware (and settings) can cause an vertical offset near the surface that cannot be distinguished from one another. Following Jenkins et al.

110
(2006), we aligned both measurements below the firn-ice transition. To this end, we compute the depth of pore closure h pc takes place, i.e. the depth at which a density of 830 kg m −3 is reached. To this end, we apply the densification model (Herron and Langway, 1980) and mean annual accumulation rate and temperature from the multi-year mean RACMO2.3 product (van Wessem et al., 2014). In our study area, h pc varies between 62 m and 71 m. The actual alignment is based on a correlation of the amplitudes for a window of 6 m around h pc . No reliable alignment could be obtained from the correlation for nine stations, 115 since the correlation values were not unambiguous. As a consequence, these stations were not considered.
After the alignment, the change in the ice thickness H i below the depth of the pore close h pc is only affected by vertical strain and basal melt. Thus the basal melt rate a b (positive for melting, negative for freezing) is with ∆H ε being the thickness change due to vertical strain ε zz . ∆H ε is derived from integrating ε zz from the aligned reflector Here, h b denotes the averaged depth of the ice base of the measurements. The vertical strain is defined as with the displacement in vertical direction u z .
In order to determine u z , we followed the method described by Stewart et al. (2019). We divided the first measurement in segments of 6 m width with 3 m overlap from a depth of 20 m below the surface to 20 m above the ice base. To determine vertical displacements, we cross-correlated each segment of the first measurement with the repeated measurement. The lag of 125 the largest amplitude correlation coefficient was used to find the correct minimum phase difference, from which we derived the vertical displacement. Since noise prevents the reliable estimation of the vertical displacement from a certain depth on, we calculated the depth at which the averaged correlation of unstacked chirps undercuts the empirical value of 0.65. We name this the noise-level depth limit h nl , which is 743 m on average in this study area. Only those segments located below h pc and above h nl were used to avoid densification processes and noise to influence the strain estimation. A linear regression was calculated 130 from the shifts of the remaining segments, assuming a constant vertical strain distribution over depth as the overall trend.
However, at six stations, all in the hinge zone where the ice is bended by tides, we observed a slight deviation from a linear trend at deeper layers (Fig. A1a). The segments that indicate a non-linear distribution are located below h nl and are hence not taken into account for the regression. Nevertheless, we want to provide a lower limit considering other forms of strain-depth relations. For this purpose, we use a strain model that is decreasing linearly from half the ice thickness (approximately h nl ) to 135 the depth of at which ε zz = 0 (Fig. A1b). This serves as a lower limit of the displacement, whereas a linear ε zz (z) gives the upper limit. The average of both gives ∆H ε and the difference the uncertainty.
In order to derive ∆H i , we used a wider segment of 10 m around the basal return, which was identified by a strong increase in amplitude. Its upper limit is located 9 m above the basal return, while the lower limit is defined 1 m below the basal return.
The vertical displacement of the ice base and thus the change in ice thickness was obtained from the cross-correlation of the 140 basal segment.
The uncertainty of the melt rate results mainly from the alignment of the repeated measurement and the uncertainty of ∆H ε .
This leads to uncertainties in the melt rate of more than 0.2 m a −1 for locations in the hinge zone, while at other locations the uncertainty is predominantly in the range of < 0.05 m a −1 .
In order to classify how representative the melt rates are for the past, we reconstructed the ice thickness based on the values velocity is constant. Next, we treat the change in ice thickness as a transport equation. To this end, we compute the advection of the ice thickness under present day climate conditions (H PDadv ). For this we use interpolated functions of a b (t), ∆H ε (t) and 150 ∆H s (t). The expected ice thickness at H PDadv is then the thickness at t 0 = 0 a plus the cumulative change in ice thickness: We can turn this around and calculate a synthetic melt rate a syn b (t) that reconstructs the ice thickness H: Descriptions of the symbols are given in Tab. A1.

Basal melting from ApRES time series
The processing of the autonomous measured time series differs slightly from the single-repeated measurements. For the ApRES 155 time series, the instrument was located below the surface, thus snow accumulation had no influence on the measured ice thickness and an alignment of the measurements is not necessary. This gives the possibility to determine the firn compaction ∆H f . Without the alignment, thickness change due to strain needs to be considered for the whole ice thickness H For processing, we followed the method described by , which differs slightly from the processing applied by Vaňková et al. (2020). Similar to processing of the single-repeated measurements, we divided the first measurement 160 into the same segments and calculated the cross-correlation of the first measurement (t 1 ) with each repeated measurement (t i ).
The displacement was obtained by the lag of the minimum phase difference. To avoid half-wavelength ambiguity due to phase wrapping, we limited the range of expected lag based on the displacement derived for the period t 1t i−1 .
The estimation of the vertical strain for the period t 1t i is based on a regression analysis of the vertical displacements for chosen segments. Only those segments located below a depth of 70 m and above the noise-level depth limit of h ≈ 600 m were 165 used to avoid densification processes and noise to influence the strain estimation. Assuming constant strain over depth (which is a first guess only), the regression analysis gives the vertical strain and the cumulative displacement u z (z) is where the intercept at the surface is the firn compaction ∆H f . By increasing the time period, the cumulative melt of the ApRES time series is derived.
In order to investigate if the basal melt is affected by tides, we first de-trended the cumulative melt time series and computed 170 the frequency spectrum afterwards. Subsequently, we used frequencies up to the solar annual constituent as input for a harmonic fit of ∆H(t). We then de-tided ∆H(t) by subtracting the harmonic fit and calculated the thinning rate. Assuming, that basal melt causes changes on short time scales of several days, we attribute abrupt increases in the thinning rate to basal melt anomalies.

175
The GPS processing is similar to the method used by Christmann et al. (2021). With the Waypoint GravNav 8.8 processing software, we applied a kinematic precise point positioning (PPP) processing for the GPS data that were stored in daily files.
We merged three successive daily solutions to enable full day overlaps avoiding jumps between individual files. Afterwards, we combined the files in the middle of each 1-day overlap using relative point to point distances and removed outliers. The data has been low-pass filtered for frequencies higher than 1/3600 Hz. For tidal analysis, we calculated the power spectrum of 180 the vertical displacement.

Digital Elevation Model (DEM)
We use the TanDEM-X PolarDEM 90m Digital Elevation data product provided by the German Aerospace Center (DLR) as reference elevation model (DLR, 2020). As the elevation values represent ellipsoidal heights relative to the WGS84 ellipsoid we refer the PolarDEM to the EIGEN-6C4 Geoid (Foerste et al., 2014). In the following, we refer the DEM heights above 185 Geoid as observed surface elevation h TDX . The absolute vertical height accuracy of the PolarDEM is validated against ICESat data and given to be < 10 m (Rizzoli et al., 2017). For our region of interest the accuracy is given to be < 5 m as shown in three times larger than those outside the channel or at the steepest eastern flank (SE) at similar draft.
The distribution of ∆H ε shows a significant thickening of more than 1 m a −1 at the most upstream cross-section at L and OE 205 (Fig. A2). In ice flow direction, ∆H ε declines, reaching about zero above the channel at the cross-section furthest downstream.
In contrast, outside the channel, strain-thinning occurred 30 km from the grounding line on. The change in ice thickness due to firn compaction and accumulation is close to zero in the entire study area (Fig. A2). However, the measurements only show a snapshot, as the variability on longer time scales is unknown. Based on the interpolated melt rates, ∆H ε and ∆H s along the channel (solid lines in Fig. 3a and A2), we computed the advected ice thickness under 210 present day climate conditions H PDadv (solid lines in Fig. 3b). The comparison of H PDadv with the measured ice thickness (dashed lines) shows large differences of up to 185 m above the channel. While the observed ice thickness decreases rapidly above the channel, H PDadv remains almost constant. In contrast, no significant differences between the observed ice thickness and H PDadv can be identified outside the channel. If the present day melt rates were representative for a longer period of time, the channel would be closed within 250 years, as the difference in H PDadv above and outside the channel reaches zero.
the channel must have been higher in the past. How large the melt rates must have been on average can be deduced from the reconstruction of the existing ice thickness. The resulting synthetic average melt rate in the channel is about twice as high as the observed ones, reaching 3.5 m a −1 in the upstream area (yellow dashed line in Fig. 3a). Assuming a steady state ice thickness upstream of the study area (supported by low elevation change found in (Helm et al., 2014)) and constant vertical strain and 220 accumulation in the past, this indicates that melt rates in the last 250 years have been significantly higher than observed now.
In addition to the observations we have presented in this section, we also conducted measurements of the vertical profile of the vertical displacement, that we present below together with simulations.

Time series of basal melting
The ApRES time series outside the melt channel reveals an average melt rate of 0.23 m a −1 (Fig. 4a). A look at the monthly 225 mean melt rates shows increased melt during the summer months (January, February and November, December) in comparison with the winter season. In these months the melt rates show values from more than 0.3 m a −1 up to 0.62 m a −1 . The unfiltered time series of the cumulative melt shows a tidal signal with amplitudes of ∼ 1 cm within 12 h around the low-pass filtered cumulative melt. The spectral analysis shows all main diurnal and semi-diurnal constituents, which is in accordance to the frequencies observed from the GPS station (Fig. A3). The analysis of melt events from the de-tided thinning rate shows several 230 melt anomalies distributed over the entire measurement period (Fig. 4a). These events lasted from a several hours to a few days and melted up to 1.5 cm of ice.  We found evidence for a clear accordance of the strain in the upper ice column with the tidal signal as recorded by GPS measurements. Unfortunately, we are lacking vertical strain in the lower column of the ice due to the noise, which permits to extract the temporal variation of basal melt rates on tidal time scales. As the tidal variation of ∆H/∆t is by far lower than 235 the observed ∆H ε /∆t, either deformation in the upper and lower parts compensates each other or basal melt/freeze takes this role. We can exclude freezing, as we do not find jumps in the amplitude of the basal return in the ApRES signal (Vaňková et al., 2021) over tidal time scales. Consequently, we infer that strain in the lower part compensates the one in the upper part and there is only a small variation of basal melt on tidal time scales.
As our location is close to two hinge zones, upstream and west of the melt channel, only a full three-dimensional model could 240 shade light into the vertical strain in the lower part of the ice column. This is numerically costly for the required non-linear strain theory and not in reach. With melt channels being located (or initiated) in the hinge zone, any kind of ApRES time series performed at ice columns with a thickness of more than 1000 m is affected by the unclear strain-depth profile in the lower part of the ice column. This may be overcome by a radar device with higher transmission power, that allows to detect the vertical displacement of layers down to the base. The observed tidal dependency of the vertical strain is consistent to the finding from 245 other ApRES locations at the Filchner-Ronne Ice Shelf by Vaňková et al. (2020). They found the strongest dependency, even of the basal melt rate at some stations, on the semidiurnal (M 2 ) constituent. Beside depth-independent tidal vertical strain, (Vaňková et al., 2020) found tidal deformation from elastic bending at ApRES stations located near grounded ice.

Viscoelastic modeling
To obtain a more profound understanding of the evolution of the channel, we conduct transient simulations and analyze the 250 change in geometry of a 2D cross-sections over time, as well as the simulated strain-field. The simulations are forced with the basal melt rates (both interpolated and synthetic) obtained in this study (Fig. 3). We transform distance to time in along flow direction of the ice shelf ( Fig. 1) using present day velocities. This enables us to study under which conditions the channel is stable or vanishes.
Ideally, we would have observations of ice geometry and basal melt rates from the transition from inland onward, but our 255 first cross-section with observations is located 14 km downstream of the grounding line (Fig. 1). The initial elastic response of the grounded ice becoming afloat had faded away. Further elastic contributions to the deformation originates from in-situ melt at the base and accumulation at the surface. To initialize our simulations adequately, we therefore conduct a spin-up.

Model
The model comprises non-linear strain theory accounting for finite deformations, as there is no justification to expect a priori 260 the deformation to be small for simulation times of more than 200 a (e.g. Haupt, 2002). We treat the ice as a viscoelastic fluid and solve the system of equations for displacements using the commercial finite element software COMSOL (Christmann et al., 2019). The constitutive relation corresponds to a Maxwell material with an elastic response on short time scales and viscous response on long time scales. For homogeneous, isotropic ice, two elastic material parameters exist (Young's modulus and Poisson's ratio). We conduct all viscoelastic simulations with commonly used values for ice for Young's modulus of 1 GPa and Poisson's ratio of 0.325 (Christmann et al., 2019). Another material parameter of the viscoelastic Maxwell material is the viscosity. It controls the viscous flow of ice. We use a constant viscosity of 5 × 10 15 Pa s and discuss the influence of this material parameter later on. This constant viscosity is at the upper limit of the viscosity distribution derived by an inversion of the rate factor in the floating part of Filchner-Ronne Ice Shelf (Appendix Sec. B1 and Fig. B1). This inversion has been conducted using the Ice Sheet and Sea-Level System Model (ISSM) (Larour et al., 2012) in higher-order Blatter-
The model geometry represents a cross-section through the melt channel ( Fig. 5) with the x-direction being across channel and resembling the seismic IV profile (Fig. 1). By assuming plane strain, it is virtually infinite in the y-direction. The com-275 putational domain is discretized by an unstructured mesh using prisms with a triangular basis involving a refined resolution near the channel. We use the direct MUMPS solver and backward differentiation formula with automatic time step control and quadratic Lagrange polynomials as shape functions for the displacements. The viscous strain is an additional internal variable in the Maxwell model and we use shape functions of linear discontinuous Lagrange type to save computational effort. In some cases, the geometry evolution leads to degraded mesh elements, which requires automated remeshing from time to time.

280
In this study, the ice density is 910 kg m −3 and the seawater density is 1028 kg m −3 . At the upper and lower boundaries, we apply stress boundary conditions: for the ice-ocean interface, a traction boundary condition specifies the water pressure by a Robin-type condition. The ice-atmosphere interface is traction-free. Laterally, we apply displacement boundary conditions. As we take a plane strain approach, we can neglect deformation in the downstream direction. To obtain reasonable lateral boundary conditions, we transfer observed vertical strain and hence, vertical displacements, in horizontal displacements assuming 285 incompressibility ε zz = −(ε xx + ε yy ), so that u x becomes with W the width of the simulated cross-section (Fig. 5). We assume that the horizontal displacements are constant in vertical direction at lateral boundaries, resulting in a compression or elongation perpendicular to the channel (Fig. B2a).
The climate forcing consists of the SMB and basal melt rate. Technically, both are applied by changing the geometry of 290 the reference configuration with the respective cumulative quantities (Fig. B2b,c). For the SMB, we used multi-year mean we conduct individual experiments that are based on our observed melt rates and its variations. As this data is spatially sparse, we need to interpolate those values in across-channel (x) direction. To this end, we assume a smooth cosine-shaped transition.
As we conduct Lagrangian experiments, distance equals time. We define t 0 = 0 a at the pRES measurements furthest upstream ( Fig. 1) that is also the location of the seismic observations by Hofstede et al. (2021a, b). To evaluate our simulations, we compare the simulated surface topography and ice thickness, as well as u z (z) with the observed one. We performed a spin-up to avoid model shocks, introduced by the transient behavior of a Maxwell material, that could be falsely interpreted as the response. The main goal is here is to have the geometry after spin-up fit reasonably to the geometry measured at the seismic IV line (see Fig. 1) that we denote as time t 0 . The spin-up covers 75 a, which corresponds to the time from the grounding line to that profile under present day flow speeds. To this end, we take a constant melt rate equal to the melt rate at t 0 and adjust the geometry at the grounding line to match the geometry at t 0 of the seismic IV profile reasonably well.  Fig. 3). The second experiment aims to derive the best match between simulated and observed geometry. For this experiment, we use synthetic melt rates (dashed lines in Fig. 3a).

First experiment: pRES-derived melt rate
The spin-up for this experiment is starts with a manually adjusted geometry (including the channel at the base) for t = −75 a 320 to fit seismic IV profile for t 0 at the base. We applied a constant melt rate of 1.5 m a −1 at L and 0.5 m a −1 at OE. This forcing enlarges the melt channel during the spin-up as the ice thickness OE increases due to the prescribed displacement at the lateral boundaries. The general shape of the base matches the seismic profile IV reasonably well ( Fig. 1 and Fig. B3). In the experiment, we force the base with a b (solid line in Fig. 3a).
The results of this experiment are displayed in Fig. 6. For both locations, L and OE, the simulated and observed geometry 325 differ significantly. While the simulated ice thickness above the channel only declines by 21 m or 1.7% in 250 a, the observed one is a factor of 9, or 191 m thinner. However, the simulated trend outside the channel shows thinning. This thinning sets in after 50 a, whereas we find continuous thinning in the observations. This delayed onset of thinning is also represented in the simulated surface topography. Most notably is the match between simulated H sim and advected H PDadv ice thickness under present day climate conditions at L. This match confirms that present day melt rates would not lead to the observed channel 330 evolution over 250 a.

Second experiment: Synthetic melt rate
The spin-up for the second experiment also starts with a geometry that has been manually adjusted for t = −75 a to fit seismic IV profile for t 0 at the base. In the second simulation experiment, we force the base with the synthetic melt rate (Fig. 3a).
Again, the melt rate has been kept constant over the spin-up with a syn b (t 0 ). By disregarding the additional melt of the spin-up, 335 the synthetic melt rate leads to a cumulative melt after 250 a of 290 m (Fig. B2a). With that 184 m more ice is melted at L than in the first experiment.
The modeled geometry of this experiment is presented in Fig. 7. The simulated ice thickness at L is in very good agreement with H pRES . There is some mismatch at OE, but the simulated trend of thinning is synchronous to the observation. After 250 a the deviation from the observed ice thickness at OE reaches 53 m. The simulated evolution of the base for the second 340 experiment shows a persistent basal channel (Fig. B4). The mismatch of the surface elevation at L and OE is reversing over time: while the simulated surface topography at OE is first too low, it is too high in the second half of the transient simulation ( Fig. 7). However, the trend of the observed h TDX and simulated h sim elevation behave similarly. Above the channel, the surface elevation is first overestimated by 4 m at the end of the spin-up. After 57 a, it turns from an over-to underestimation that results in an 8 m lower h sim than the observed h TDX after 250 a. To assess if the ice is in buoyancy equilibrium, we take another approach and estimate the mean density under the assumption of buoyancy equilibrium: at t 0 this corresponds to 901 kg m −3 and after 250 a to 896 kg m −3 . As more ice is melted from below and with higher snow accumulation at L, the density decreases, which is to be expected.

350
After 250 a, the simulated freeboard at OE is 1 m higher than the surface elevation of 138 m inferred by buoyancy equilibrium using an ice density of 910 kg m −3 . Similar considerations regarding OE lead to a 3 m higher h TDX than 132 m out of buoyancy equilibrium. Overall, we see convergence to equilibrium state at OE and the simulated surface elevation at L. At the end of the simulation, only h TDX above the channel does not reach buoyancy equilibrium, which leads to the justifiable assumption that the mean ice density at L is lower than OE.

355
At the position of the furthest upstream pRES observations we can see from the seismic IV profile that the influence of the grounding line has not completely vanished. The assumption of buoyant equilibrium is therefore likely to be flawed. At the end of the simulation, the geometry should be close to buoyancy equilibrium despite melting and a 50% higher SMB at L than OE. Hence, simulations carried out using a higher SMB within the channel would result in better agreement with the observed values of h TDX . Next, we exploit the variation of the vertical displacement over depth. The results are presented in Fig. 8. For this purpose, we calculated the cumulative vertical displacement in one year over depth. For comparability, the vertical displacements due to accumulation and snow compaction were removed from the observed distributions.
Most notably, we move from a vertically extensive regime into a compressive by increasing distance to the grounding line.

365
Given the complexity of the problem, the simulations show a reasonable agreement with the observations. The best match is reached at OE, which is not that surprising. The generally good agreement of the simulated displacements outside the channel comes from tuning u x at the lateral boundary to match u z from the pRES measurements at OE. Both simulated and observed vertical displacement distributions show that the strain decreases from L to OE. The only exception here is t = 57 a, where the vertical strain at SE is larger than the one at L. While at 0 a and 26 a the deviation of the simulated displacements between L and 370 OE is small, it increases afterwards. From 105 a, the simulated vertical displacements agree very well with those of the pRESmeasurements, where a displacement distribution was derivable at L and OE. The same comparison for the first experiment ( Fig. B5) shows similar results, with significantly less pronounced differences between L and OE. Hence, the mismatch to the observed vertical displacements for this experiment using the measured melt rates is higher than for the second experiment with the synthetic melt rates. The simulated strain evolution of the cross-section confirms that the viscoelastic model needs to 375 account for finite deformation as the strain exceeds 10% (Fig. B6).
As the last point of this second experiment, we consider the influence of the viscosity on the evolution of the melt channel (Fig. B7). To reach the ice thickness of seismic IV, the simulation applying the smallest viscosity needs a higher initial channel and hence, less ice thickness at L at the beginning of the spin-up (Appendix Sec. B2). The channel thickness of the pRESmeasurement is modeled best using a viscosity of 5 × 10 15 Pa s. A two times higher viscosity leads to an ice thickness in the 380 channel that is 42 m smaller after 250 a, while a five times lower viscosity results in 116 m thicker ice above the channel due to more viscous flow into the channel. The simulated ice thickness OE is similar for all three different viscosities. study. In general, the channel height declines, so the channel fades out. The channel diminishes by melt rates inside the channel falling below those outside the channel. The trend in vertical strain has only a minor contribution to this evolution. We thus do not find any evidence that such channels are a cause for instabilities of ice shelves as suggested by Dow et al. (2018).
One of the main findings of our study is that the present geometry can only be formed with considerable higher melt rates in the past (see Fig. 3). This finding is based on the assumption that the strain-rates were in the past similar to present day. This is 400 justified, as significant changes in strain would require a change in the system that would cause other characteristics to change, like the main flow direction, for which we do not find any indication.
The pRES melt rate observations covered only one year. As the ocean conditions within the sub-ice shelf cavity are known to respond to the ocean forcing from the ice front (e.g. Nicholls, 1997), they would expect to be subject to significant interannual variability. Underlying any interannual variability, a long-term reduction in basal melt rates would be the expected response to 405 a reduction in production of dense shelf waters north of the ice front, resulting from a reduction in sea ice formation (Nicholls, 1997), resulting in turn from a reduction in the southerly winds that blow freshly produced sea ice to the north.
A decrease in northward motion of sea ice has been observed in the satellite record (e.g. Holland and Kwok, 2012), and modeling experiments by Naughten et al. (2021) find decreasing basal melt rates. This reduction is therefore consistent with higher basal melt rates in the past. However, our model results suggest that the mismatch between the past melt rates needed to explain the observed channel geometry and those that were observed applies only to the channel, and not to the ambient ice. This could be explained by historically higher levels of subglacial outflow at the grounding line, or anomalously low levels during the observation period. Subglacial outflow contributes to the buoyant flow up the basal slope and therefore the shearinduced turbulence that raises warm water from deeper in the water column towards the ice base. Smith et al. (2009) found an active subglacial lake at the transition between Academy Glacier and SFG, and also Humbert et al. (2018) suggest in the 415 upstream area of SFG a subglacial lake. Hofstede et al. (2021a) showed that the subglacial channel appears 7 km upstream of the grounding line increasing its height to 280 m at the grounding line. The location of the channel corresponds with increased subglacial flux found by Humbert et al.
(2018) using a simple routing scheme. The channel formed on the grounded part is most likely the source of a grounding line fan and thus carrying sediments, formed at the seabed under the basal channel Hofstede et al. (2021a). Once this topographic 420 feature reaches the ocean, it will focus on the relatively buoyant flow and enhance shear-driven vertical mixing, bringing heat and salt to the ice base leading to higher basal melt rates.
However, with increasing distance along the channel, the basal gradient, and therefore the speed of the buoyant flow, reduces, as does the entrainment of warm water from beneath. Coupled with the increasing pressure freezing point at the ice base this leads to a gradual reduction of the melt rate in the channel. From Fig. 2a, the melt rate in the channel reduces below that of 425 the ambient ice base by about 30 km distance from the grounding line, suggesting that the effect of basal melting thereafter is to suppress the channel. The cause of the strong melt anomalies identified in the ApRES measurements remains unclear as no direct ocean observation exists near SFG. However, the time scale of the events is consistent with the passage of warm cored eddies. Such features have been observed in the ocean cavity beneath the neighboring Ronne Ice Shelf (Nicholls, 2018).
The channel height is found to increase until 30-35 km downstream of the grounding line. Further downstream, the channel 430 begins to close. Our modeling results show that less viscous ice (1 × 10 15 Pa s) would tend to shut the channel faster than the rate we observe. For the best match between observed and modeled geometry, we need viscosities around 5 × 10 15 Pa s of stiff ice to prevent the closure by deformation. This claim is also confirmed by the inversion of the viscosity to model observed The difference in geometry change, due to different values of the viscosity, manifests stronger inside the channel than outside.
This was to be expected because of the load situation resulting from the prescribed geometry (Fig. B7). The simulated geometry change is mainly due to the elastic response to thinning by basal melt and ice accumulation. Any purely viscous simulation We find a difference (−4 m to 8 m) between simulated and observed surface elevation at L. The elevation difference is most likely caused by the constant density that we used for the simulations, as the ice thickness matches well. For the thinner ice 445 above the channel, this could be achieved by an ice density decreasing from outside to inside and from upstream to downstream the channel. However, one has to keep in mind, that the accuracy of the surface elevation product is 5 m.
In general, we benefited highly from having measurements of vertical strain available. This opens new possibilities to identify weaknesses in the modeling and gave us useful insight into the spatial variation of the vertical strain across such a topographic feature. Although the pRES surveys only about half the ice thickness, the slope of u z (z) in the upper half is distinct for the 450 positions L, SE and OE and greatly varies with distance from the grounding line, also influenced by the embayment of the ice shelf. Simulated u z at L start to match only after about 100 a well with observations, which could result from the first few cross-sections still being influenced by the hinge zone. Tidal bending was not taken into account here, due to the 2D setting.
This could in future be investigated, if repeated pRES measurements would be conducted up to the grounding line covering the entire hinge zone, in which it would also be extremely advantageous to obtain basal melt rates.

455
Our study demonstrates that viscoelastic simulations can be a useful but complex tool to analyze melt channel evolution. In an inverse approach, viscoelastic models could also give more insights into basal melt rates of channel systems of ice shelves in general, given that satellite-borne surface elevation is available in high resolution. However, the fact that large deformations require non-linear strain theory will make this a challenging endeavor. As changes in basal melt rates will inevitably lead to surface elevation changes of channel systems, systematic monitoring of the surface topography from space can serve as an 460 early warning system and trigger further in-situ observation similar to this study.

Conclusions
We find basal melt rates in a melt channel and its surroundings on Filchner Ice Shelf to be up to 2.3 m a −1 . Basal melt rates inside the channel drop with distance down-flow, even turning into freezing 55 km downstream of the grounding line. Close to the grounding line, melt rates are larger inside the channel than outside, while further downstream this relationship reverses.

465
Over distance along flow, the channel dimension decreases from a maximum height of 330 m to below 100 m. The channel diminishes because the reduced melt rate is unable to maintain the channel geometry against viscoelastic deformation. Analysis of the present day ice thickness advection revealed large differences compared to the observed ice thickness above the channel, which indicates that melt rates have been about twice as large in the last 250 a. The viscoelastic simulation confirms this statement and indicates that basal melt channels need high basal melt rates and relatively cold ice to persist. The deformation 470 of the basal melt channel is mainly driven by the elastic response to the basal melt rate. The observed and simulated evolution of this melt channel demonstrates that melt channels of this kind are not a destabilizing element of ice shelves. The ApRES time series showed brief melt anomalies distributed over the entire measurement period and slightly increased melt rates in summer. simulation used for this study is available via AWI's gitlab (https://gitlab.awi.de/jchristm/viscoelastic-finite-defos-meltchannel).
Data availability. Raw data and derived products of the single-repeated pRES measurements, raw data of the ApRES time series (https://doi. done on an unstructured finite element grid with a refined resolution of 2 km at the grounding line, in the shear margins as well as at other regions of faster ice flow. In the melt channel domain we further refine the resolution of the grid to 0.5 km. To generate the geometry of the ice shelf the BedMachine Antarctica v2 data set is used (Morlighem et al., 2020;Morlighem, 2020). For the ice rigidity in the grounded region, as well as an initial guess of ice rigidity in the floating shelf, we assume the results of a long-term thermal spin-up also used in Eisen et al. (2020) based on the geothermal flux from Martos et al. (2017). 495 We constrain ice surface velocities to fit the MEaSUREs data set (Mouginot et al., 2019b, a).
Our optimization approach iteratively infers two parameters -the basal friction parameter in the grounded area and the ice rheology parameter in the floating area. For this purpose two cost functions are built. Each cost function consists of two data misfits, linear and logarithmic, as well as a Tikhonov regularization term: The first term will be most sensitive to velocity observations in fast-flowing areas, the second term will be most sensitive to velocity observations in slow-floawing areas, while the third term penalizes oscillations in the optimization parameter. We performed an L-curve analysis to find suitable weights γ 1 , γ 2 , γ t for both cost functions. With this trade-off curve, we can make sure that we find a regularization term that fits the data well without overfitting noise. For the basal friction inversion, we found best weights γ 1 = 1, γ 2 = 5 × 10 −6 and γ t = 1 × 10 −8 , while for the ice rigidity inversion the optimal weights were 505 γ 1 = 1, γ 2 = 0.8 and γ t = 4 × 10 −17 .
We linearize and solve the optimization problem using the M1QN3 algorithm with an incomplete adjoint (Larour et al., 2012).
For this inversion setting we apply a gradient relative convergence criterion gttol = 10 −6 and two points which are less then dxmin = 10 −4 from each other are considered identical. Besides we used a maximum number of iterations and function evaluations of 1000. We show our best-fit results for ice viscosity in the region around the melt channel in Fig. B1. The range of 510 the viscosity is between 5.0563 × 10 13 and 2.6656 × 10 15 Pa s. Figure B1. Ice viscosity in the melt channel area obtained from inverse modeling. The map extent is the same as in Fig. 1

B2 Sensitivity of experiment 2 on viscosity
To capture the influence of the viscosity, different constant values (one smaller and one higher as in the second experiment) are investigated in a further experiment. The spin-up for each viscosity starts at t = −75 a with an arbitrary basal geometry that should fit seismic IV profile at the end of the spin-up (t 0 ). The melt rate a syn b (t 0 ) is again assumed to be constant over the 515 spin-up for all different viscosity values. We force the base with the synthetic melt rate (Fig. 3a), the same melt rate we already used in the second experiment. The initial base for the middle and high viscosity is nearly the same as 5 × 10 15 Pa s is for ice a rather high value requiring cold ice (Fig. B7). For the smallest viscosity, a deeper channel at the beginning of the spin-up is needed.