Evolution of ice-shelf channels in Antarctic ice shelves

Ice shelves buttress the continental ice flux and mediate ice–ocean interactions. They are often traversed by channels in which basal melting is enhanced, impacting iceshelf stability. Here, channel evolution is investigated using a transient, three-dimensional full Stokes model and geophysical data collected on the Roi Baudouin Ice Shelf (RBIS), Antarctica. The modeling confirms basal melting as a feasible mechanism for channel creation, although channels may also advect without melting for many tens of kilometers. Channels can be out of hydrostatic equilibrium depending on their width and the upstream melt history. Inverting surface elevation for ice thickness using hydrostatic equilibrium in those areas is erroneous, and corresponding observational evidence is presented at RBIS by comparing the hydrostatically inverted ice thickness with radar measurements. The model shows that channelized melting imprints the flow field characteristically, which can result in enhanced horizontal shearing across channels. This is exemplified for a channel at RBIS using observed surface velocities and opens up the possibility to classify channelized melting from space, an important step towards incorporating these effects in ice–ocean models.


Introduction
Almost three-quarters of the Antarctic ice-sheet boundary is in contact with the ocean (Bindschadler et al., 2011) where the majority of ice-mass loss occurs (Rignot et al., 2013;Depoorter et al., 2013).Floating ice shelves extend from the continental ice seawards, providing an interface for melting and refreezing processes at the ice-shelf base.Confined ice shelves transmit stresses landwards, controlling the grounding-line position and the continental mass flux (Gudmundsson, 2013;Pattyn et al., 2013).In the extreme case of an ice-shelf collapse (as observed on the Antarctic Peninsula), the regulating effect is entirely lost and leads to an acceleration of the grounded ice flow, hence sea level rise (Scambos et al., 2004).In some cases in West Antarctica, reduced ice-shelf buttressing has already led to irreversible mass loss (Joughin et al., 2014;Mouginot et al., 2014).
The ice-shelf integrity is influenced by fractures and rifts (e.g., Borstad et al., 2013), as well as by the formation of ice mélange and marine ice with different rheological properties than meteoric ice (Khazendar et al., 2009;Pattyn et al., 2012;Dierckx and Tison, 2013;Kulessa et al., 2014).Melting and refreezing processes at the ice-shelf base are governed by the presence of ice-shelf water.The latter is formed inside the ice-shelf cavity by mixing of fresh water with other water sources such as high-salinity shelf water or circumpolar deep water.The fresh water originates from melting at the ice-ocean interface and/or from basal water of the ice sheet's interior exiting at the grounding line (Jenkins, 2011).Regardless of the specific origin, ice-shelf water is part of a buoyancy-driven circulation described as the "ice pump" (Lewis and Perkin, 1986) which typically leads to high basal melt rates near the grounding line and also to accretion of marine ice elsewhere.These ice-ocean interactions trigger a complex response of the continental ice sheet; knowledge of the basal mass exchange, and its spatial variability in particular, is important for predicting the ice-sheet evolution (Gagliardini et al., 2010).
A distinct feature of many ice shelves is longitudinal channels.They often start near the grounding line, extend many tens of kilometers towards the ice-shelf front (Fig. 1a), and vary between a few hundred meters to a few kilometers in width.Channels may form spontaneously from transverse variability in ice thickness (Sergienko, 2013), and/or originate from subglacial water outlets of in-flowing glaciers (Le Brocq et al., 2013).In both cases, current understand-R.Drews: Evolution of ice-shelf channels ing suggests that the initial thickness perturbations near the grounding line are amplified farther downstream through a buoyancy-driven plume with enhanced basal melting within the channels.Channelized melting, with significantly higher melt rates inside than outside the channels, has been reported for Petermann Glacier Ice Shelf in Greenland (Rignot and Steffen, 2008;Dutrieux et al., 2014), and for Pine Island Ice Shelf (Vaughan et al., 2012;Mankoff et al., 2012;Stanton et al., 2013;Dutrieux et al., 2013Dutrieux et al., , 2014) ) and for Fimbul Ice Shelf (Langley et al., 2014) in Antarctica.Channels may weaken ice shelves structurally through crevasse formation (Vaughan et al., 2012), or by breaking up entirely (Rignot and Steffen, 2008).Channels can also stabilize ice shelves, by preventing area-wide basal melting (Gladish et al., 2012;Millgate et al., 2013).In both cases, the basal mass balance inside the channels is key.
Despite the importance of ice-shelf channels for ice-ocean interaction and for ice-shelf stability, little is known about their evolution from a glaciological perspective.This is because ice shelves are typically simulated in a computationally efficient manner by applying a number of simplifications to the Stokes equations.Most common is the hydrostatic approximation (e.g., Greve and Blatter, 2009, p. 117) which assumes in the vertical momentum balance that the shear-stress gradients ∂σ xz ∂x and ∂σ yz ∂y play no role in balancing the gravitational force (given by the ice density ρ i and the gravitational acceleration g acting in the vertical direction z).In that case, the vertical normal stress σ zz is purely determined by the weight of the overburden ice column.Further approximations neglect vertical gradients of the horizontal velocities, leading to a vertically integrated set of equations with a linear profile of the vertical velocity (the shallow-shelf approximation, Morland, 1987;MacAyeal, 1989).While those approximations are wellsuited to describe the large-scale flow of ice shelves, they fail to capture the full dynamics in ice-shelf channels where ice thickness and stresses change substantially over horizontal distances which are comparable to the mean ice thickness.For example, the negligence of the shear-stress gradients entails that ice is locally balanced by the water pressure (i.e. it is in hydrostatic equilibrium).However, observations indicate (Vaughan et al., 2012;Dutrieux et al., 2013) that some channels are out of hydrostatic equilibrium due to lateral stress transfer along the channel walls, an effect referred to as bridging (van der Veen, 2013, p. 57).
In this study we use a three-dimensional, transient, full Stokes model (Elmer/Ice, Gagliardini et al., 2013) pursuing three main questions: (1) how do channel amplitudes evolve as a function of basal melting and horizontal advection?, (2) under which conditions are channels out of hydrostatic equilibrium and how does the imbalance develop downstream?, and (3) what is the imprint of the channel formation on the flow pattern of the surrounding ice shelf?Once a channel is formed, ice starts to converge into the channels (Dutrieux et al., 2013), presumably causing a reduction of the channel's amplitude while it is being advected downstream.Simulating those effects helps to pinpoint the channel state (active channelized melting or simple advection) as a function of its longitudinal position in the ice shelf.The second question is needed for deriving the basal mass balance in the channels using the continuity equation either in an Eulerian (e.g., Rignot and Steffen, 2008), or a Lagrangian (e.g., Dutrieux et al., 2013) framework.Both approaches require the ice thickness in the channels which is typically inferred from surface elevation using hydrostatic equilibrium.This assumption fails if bridging stresses prevent a full adjustment of the ice surface to hydrostatic equilibrium.Thirdly, observations (Dutrieux et al., 2013) show that melting channels imprint the surrounding flow field in a characteristic fashion.Simulating these patterns in an idealized geometry helps to better understand the mechanisms involved, and opens up the possibility to exploit the observed velocity anomalies for the classification of channels from high-resolution surface velocities.
We approach questions (1)-( 3) using the full Stokes model and by comparing the results with radar and GPS data collected in 2012 and 2013 on the Roi Baudouin Ice Shelf (RBIS), Dronning Maud Land, Antarctica (Fig. 1).Section 2 details the methodology regarding surface elevation and ice flow (Sect.2.1), the hydrostatic thickness (Sect.2.2), the radar thickness (Sect.2.3), and the model setup (Sect.2.4).We then investigate the evolution of channel amplitudes (Sect.3.1), hydrostacy (Sect.3.2), and the imprint of channels on the surrounding ice flow (Sect.3.3).We discuss the results in Sect. 4 and close by drawing conclusions regarding the role and characteristics of channels in ice shelves in Sect. 5.

Measuring of surface elevation and ice flow
Geodetic GPS data were collected to position radar profiles (kinematic mode) and for measuring ice flow (static mode).All data were processed differentially relative to a base station on a grounded pinning point with velocities of less than 3 cm per day (Fig. 1b).Daily base-station coordinates stem from precise point positioning of the Canadian Geodetic Survey.The solutions agree within centimeters with results from the Atomium software (Defraigne et al., 2008).
The kinematic elevations were corrected for the tidal displacement using the CATS02.01model from Padman et al. (2002).The model was verified by comparing the simulated vertical displacements with the observed vertical displacements of GPS markers which measured continuously over The Cryosphere, 9, 1169-1181, 2015 an 8-day period.Differences between the model and the observations are 0.01 ± 14 cm (mean and standard deviation).The de-tided kinematic profiles differ in height at cross-over points with 0.08 ± 0.7 m.For measuring ice flow, seven 3 m long stakes were installed in 2013 along a 4 km long cross section, traversing a channel (profile O-O' in Fig. 1).The GPS antennas were mounted on top of the stakes, and measured for at least 30 min.About 3-4 stakes were occupied simultaneously, in order to position the stakes in a network approach using closed loops of relative baselines.The processing was done using the GAMIT, GLOBK/TRACK v.10.5 (Herring et al., 2013) software package, and followed the procedure described in Bergeot et al. (2009) including GPS dual frequency observables, precise GPS orbits, and absolute phase center corrections for the ground and satellites antennas.Positioning errors reported from the software are within a few centimeters.The measurements were repeated after 7 days to measure the marker's displacement.Because this time interval is comparatively short to infer velocity differences between the markers, the same stake array was revisited in 2014.In order to cross-check the results with a different processing technique, we used precise point positioning of the Canadian Geodetic Survey to infer the yearly averaged marker velocities.Both sets of measurements show spatially the same pattern, and differ in magnitude by less than 3 %.Regardless of whether this constant offset reflects processing uncertainties and/or a real change in flow behavior between the weekly and the yearly averages, the uncertainties are small enough to use the GPS-derived velocities as ground truth for satellite-based velocities.

Determination of the hydrostatic thickness
The hydrostatic thickness H of a freely floating ice shelf can be derived from the freeboard height h asl (i.e. the ice shelf's height relative to the ocean surface) using the buoyancy formula with the densities of sea water (ρ w = 1027 kg m −3 ), ice (ρ i = 918 kg m −3 ), and air (ρ a = 2 kg m −3 ): To account for lower-density firn, the ice shelf is decomposed into layers of ice (H i ) and air (H a ), so that ρ(z)H ≡ ρ i H i + ρ a H a holds for the depth-averaged density ρ(z).In areas where the ice shelf is not freely floating (i.e., where the ice column is not only balanced by the water pressure, but also by stresses within the ice as is the case near the grounding line and potentially in ice-shelf channels) this equation is not valid.This will be investigated in more detail in Sects.3 and 4.
The freeboard height follows from the GPS height (z, referenced to the WGS84 ellipsoid) by correcting for the geoid (G) and the mean dynamic topography (MDT): h asl = z − (G + MDT).Here, the EIGEN-GL04C geoid (Förste et al., 2008) is used which, at RBIS, varies between 18 and 19 m and deviates less than 0.3 m from the EGM2008 ellipsoid (Pavlis et al., 2012).However, geodetic measurements in nearby Breid Bay (Shibuya et al., 1999) result in a geoid height of approximately 22 m, a value which is close to the presumably outdated EGM96 geoid (developed by NASA, NIMA and Ohio State University).The MDT corrects for long-term differences between geoid and ocean surface.At RBIS it is estimated with −1.2 m (Andersen and Knudsen, 2009), but uncertainties of the MDT are large for www.the-cryosphere.net/9/1169/2015/The Cryosphere, 9, 1169-1181, 2015 the Antarctic perimeter (O.Andersen, personal communication, 2014).
The equivalent air content (H a ) is 15 m, based on icecore measurements (Hubbard et al., 2013), and approximately 12 m based on atmospheric modeling (Ligtenberg et al., 2014).Differences may partially reflect spatial variability, given that the ice core's location, close to both the ice-shelf front and a rift system, is not fully representative of the research grid (Fig. 1b).
The dominating uncertainties for obtaining the hydrostatic thickness are hence rooted in the determination of h asl (mainly due to uncertainties in G, MDT, and cross-over errors of the kinematic GPS profiles) and H a (reflecting uncertainties of the depth-averaged density).Uncertainties of the geoid and the MDT have little spatial variation within the research grid, and biases carry over to the hydrostatic thickness as constant offsets without spatial dependence.To represent the combined uncertainty, we calculate a lower and an upper estimate (the lower estimate assuming h asl − 2.0 m and H a = 15 m, the upper estimate assuming h asl + 2.0 m and H a = 12 m).The mean difference between those two cases is 50 m.

Determination of the radar thickness and internal layering
The thickness of ice can be determined with radar by measuring the two-way travel time of a radar pulse which is emitted at the surface and reflected at the base.For ice shelves, the radar thickness is ambiguous because basal reflections may either originate from a meteoric ice-ocean interface or from a meteoric ice-marine ice interface (e.g., Blindow, 1994;Kulessa et al., 2014).For the latter case, the radar thickness is smaller than the hydrostatic thickness, and the difference corresponds to the thickness of the marine ice layer.For a reliable detection, the marine ice layer should have a thickness larger than 50 m (measured in the ice equivalent density of meteoric ice), given the uncertainties derived in the previous section.
The radar thickness was measured using a 10 MHz radar with resistively loaded dipole antennas (e.g., Matsuoka et al., 2012).Kirchhoff depth migration was applied to account for off-nadir reflections from slanted channel walls.For the examples presented here, the non-migrated data show reflection hyperbolas from crevasse tops.The travel time to depth conversion uses the pure ice velocity (v i = 168 m µs −1 ) and corrects for faster wave propagation in firn by adding 6.5 m.The correction is based on the ice-core density and an empirical density-permittivity model (Kovacs et al., 1995).The previously discussed uncertainty of the depth-density profile (corresponding to the uncertainty in H a ) changes the radar thickness by less than 2 m.Post-processing for visualizing the deeper internal layers (> 50 m depth) included low-cut filtering, bandpass filtering, and a depth-variable gain function.The shallow layering (< 50 m depth) was measured with a 400 MHz radar (GSSI:SIR 3000) which occupied the same profile lines as the 10 MHz radar.The post-processing is that used for the 10 MHz data.

Model setup
The ice-shelf model is based on the ones used in previous studies (e.g., Durand et al., 2009;Favier et al., 2012Favier et al., , 2014)), but simplified by excluding grounding-line dynamics and thermal effects.It is three-dimensional, transient, and full Stokes (based on Elmer/Ice, Gagliardini et al., 2013) with a Glen rheology (Cuffey and Paterson, 2010, p. 51), and has a constant surface mass balance (0.3 m a −1 ) and uniform temperature (−10 • C).Ocean pressure is applied at the front and at the bottom.The buoyancy pressure at time t is determined by the unknown ice-shelf geometry of this time step.Because small hydrostatic imbalances cause large vertical velocities, the geometry cannot easily be approximated with the previous time step.The system is therefore stabilized by introducing a time-dependent scheme for the bottom interface evolution (a.k.a.viscous spring), which is explained in more detail in Durand et al. (2009, Sect. 3.4).
At the lateral boundaries the normal velocities are zero, and for the unconstrained case, horizontal shearing parallel to the boundary is equally zero.The constrained case includes horizontal shearing by using a constant slip coefficient with a linear sliding law.The surface and the bottom can move freely.The landward boundary is in hydrostatic equilibrium; initial inflow velocities and ice thickness are 100 m a −1 and 500 m, respectively.Lateral and vertical velocities are initialized with 0 m a −1 and not further imposed.The thickness at the landward boundary changes during the evolution, and the inflow velocities are adapted so that the mass flux remains constant.The ice shelf is evolved without basal melting until it reaches a steady-state geometry (for the unconstrained cases the steady-state landward thickness and inflow velocities are 436 and 115 m a −1 , respectively).Channels are then initiated through localized melting near the upstream boundary, following the suggestions of Jenkins (2011) and Le Brocq et al. (2013).The prescribed melt rates have the following form: using a coordinate system where x is along flow, y across flow, and z vertically upwards.Across flow, M is a Gaussian curve with root-mean-square width σ y and peak A. The curve is extended along flow with a determining the extension length, and b the steepness of M along x to reach A; x 0 and y 0 describe the center of melting in x and y direction, respectively.Table 1 and Fig. 2 summarize the characteristics and parameters of the individual simulations discussed below.The forced hydrostatic equilibrium at the landward boundary saves computation time, but is unphysical because grounding lines are typically not in hydrostatic equilibrium (e.g., Lestringant, 1994;Durand et al., 2009;Bindschadler et al., 2011).The applied simplifications cause small hydrostatic imbalances at the lateral boundaries which, however, do not imprint the ice-shelf center where the channel evolution is investigated.Because the landward boundary is freely floating, channels are carved into the steady-state ice shelf with an approximately 4 km seaward offset.This hampers the investigation of how channels behave directly at the grounding line, but does not impact the channel evolution farther downstream.The latter is the focus of this study, as field data have been collected relatively close to the ice-shelf front (Fig. 1).
The synthetic ice-shelf geometry applied here hampers a one-to-one comparison with the field data and excludes quantification of basal melt rates.In turn, the simplified geometry allows distinguishing more easily between the different mechanisms acting during the channel creation, advection, and decay.Laterally unconstrained cases are used to study the channel amplitudes and bridging effects.A laterally constrained case, with reduced longitudinal stretching, is used for investigating the channel's imprint on the surrounding flow field.This is done to provide a scenario for an effect seen in the observed surface velocities (i.e.increased horizontal shearing across channels) which is not apparent in the simulations when longitudinal stretching is too dominant.

Development of channel amplitudes
Figure 3 displays the along-flow development of channel amplitudes in an unconfined ice shelf for three melt scenarios differing in the peak magnitude and longitudinal extent (A = 14, 2.3, and 1.6 m a −1 over a = 1, 6.5, and 8.5 km, respectively; MS 1-3, Table 1 and Fig. 2a-c).The channel amplitudes A ch are defined as with H in , and H out for the ice thickness at the channel trough and the channel keel, respectively.The melt functions M(x, y) are chosen so that at the downstream end of the longest prescribed melting (here at distance L = 20 km from the grounding line), the cumulative melt rate is equal for all cases.Along the central flow line (y = 0), this means 3a).
The cumulative effective melt (M eff., Fig. 3b) for each scenario is different, because the residence time in which ice is subject to melt also depends on the modeled along-flow velocities v x : Figure 3c show that channel amplitudes increase in areas where channelized melting is sustained and that they decay when melting ceases.The amplitudes at the downstream end differ less than 10 % for all scenarios.

Channels and hydrostatic equilibrium
In order to check whether the simulated channels are in hydrostatic equilibrium, we invert the simulated surface topography for ice thickness using Eq. ( 2) (with H a = 0 because firn is excluded in the simulations) and compare this hydrostatic thickness with the modeled thickness.Channels are in hydrostatic equilibrium if both thicknesses are equal.The channels investigated in the previous section have a keel-tokeel width of approximately 2 km (σ y1 = 500 m) and are essentially in hydrostatic equilibrium everywhere.Some channels at RBIS are narrower, and we investigate bridging by considering two channels with σ y1 = 100 m and σ y2 = 75 m, using otherwise the same boundary conditions as above (MS4, Table 1).This setup is comparable to channels observed at RBIS.In the simulations, we observe no significant interaction between the two channels, and the varying channel width is the main difference regarding the channel evolution.Figure 4a shows that, unlike the wider channel, the narrower channel is out of hydrostatic equilibrium for tens of kilometers along flow.The imbalance increases along flow in areas of channelized melting and decreases farther downstream.Bridging stresses are larger in the narrower channel than in the wider one (Fig. 4b), the surface is less depressed in the narrower channels (Fig. 4c), and the hydrostatic thickness deviates from the modeled thickness inside the narrower channel (Fig. 4d).

Imprint of channels on the surrounding flow field
The imprint of channel formation on the surrounding flow field is depicted for a laterally constrained ice shelf (MS5, Table 1) in Fig. 6a-f.The large-scale pattern of the alongflow velocities and the along-flow strain rates (˙ xx ) is inconspicuous (Fig. 6a and d).Across-flow velocities (which are zero without a channel), however, point towards the channel with magnitudes of a few meters per year (Fig. 6b).This convergence is prominent in the across-flow and verticalstrain rates (˙ yy and ˙ zz , Fig. 6d and e), both showing an increasing anomaly up to the point where channelized melting stops.Without channelized melting, the characteristic patterns decay at greater distances.The otherwise smoothly varying along-flow velocities increase step-wise in across-flow profiles (Fig. 7a) near areas of channelized melting.The velocity step is evident in lateral shearing (˙ xy ), and accompanied by a peak in the effective strain rate . A comparable feature is found in the field data: Fig. 7b displays the measured speed for a line of markers, densely placed across a channel (O-O', Fig. 1).The marker's displacements were measured after 1 week and after 1 year.Both type of measurements show the same tendency: speed increases from east to west, following the general ice-flow pattern in this area.Near the channel's trough, however, the velocity increase is anomalously larger than the underlying east-west tendency.The signal is equally apparent in four different sets of satellite-based velocities derived from different techniques, time intervals, and sensors (1speckle and phase offset tracking in RADARSAT data from 2000 (Callens et al., 2014); 2 -mosaicked velocities by Rignot et al. (2011); 3 -interferometric synthetic aperture radar (SAR) from European Remote Sensing (ERS) satellites in 1996 (S.Berger, personal communication, 2014); 4speckle tracking using data from the Advanced Land Observing Satellite in 2010 (S.Berger, personal communication, 2014)).Strain rates derived from the ERS data surrounding the stake array are depicted in Fig. 8a-c.

Discussion
The simulations presented here confirm that channelized basal melting is a feasible mechanism to transform initially small thickness perturbations near the grounding line into mature channels farther downstream.In that sense, the full Stokes model supports, from a glaciological perspective, the studies of Gladish et al. (2012) and Sergienko (2013), both using simplified ice dynamics to allow for an efficient coupling with ocean models.The along-flow amplification rate of the channel amplitudes is primarily determined by the upstream ratio of basal melting and along-flow advection (determining the effective melt M eff.).The maximum amplitude of MS1 is largest, because all basal melt occurs where ice is slow; the ratios between the maximum effective melt of the different scenarios are equal to the ratios between the corresponding maximum channel amplitudes (Fig. 3b and c).Without basal melting, channels decay.In the unconstrained case considered here, horizontal advection is fast, and the channels still have a considerable amplitude when reaching the ice-shelf front.In the constrained case (not shown here), horizontal advection is reduced resulting in smaller final amplitudes for the same melt scenarios.Because channels can be sustained without basal melting for many tens of kilometers, their mere existence in satellite imagery or radar data is not a sufficient condition to infer channelized melting at this location.Conversely, the disappearance of a channel does not necessarily imply basal accretion, and more data are required to characterize the channel's state.A typical approach to obtain the channel geometry is the hydrostatic inversion of highly resolved elevation models, a method which fails if bridging stresses prevent a full adjustment of the channel surface to hydrostatic equilibrium.
In order to compare the geometries of the simulated channels with the observations at RBIS, we use the ratio (α) between the keel-to-keel width and the ice thickness at the channel keel.Figure 4 illustrates that for steady-state conditions, and for the range of melt parameters considered here, channels with α > 5 are essentially in hydrostatic equilibrium.Narrower channels (α < 5) can deviate from hy-   drostatic equilibrium.In Fig. 4a, the imbalance is largest at the downstream end of channelized melting, and contains a memory of the upstream melt history: melt scenarios with higher peak values over a shorter longitudinal extent (e.g., MS1-type vs. MS4-type, Fig. 2) produce channels with comparable α, but with a smaller imbalance.Regardless of the upstream melt history, a typical feature for simulated channels that are out of hydrostatic equilibrium is that the surface depression is wider and shallower (Fig. 4c) than what would be suggested from hydrostatic equilibrium.Consequently, the hydrostatic thickness is larger at the channel trough, and smaller at the channel flanks.The spatial pattern of the imbalance may be different when melting is concentrated on the channel's walls only (Dutrieux et al., 2013(Dutrieux et al., , 2014)).
To validate whether bridging does play an important role for channels observed at RBIS, we investigate the anomalously large hydrostatic thickness in Ch.-2 on profile line B-B' (Fig. 5d): while Ch.-1 and Ch.-2 are almost equally incised at the ice-shelf bottom, the surface of Ch.-2 is less depressed than that of Ch.-1 (Fig. 5c).This causes a mismatch between radar thickness and hydrostatic thickness, which is significant because, with the exception of Ch.-2 on B-B', the radar thickness is within the error bounds of the hydrostatic thickness (including Ch.-1 and Ch.-2 on profile line A-A').
For the large hydrostatic thickness in Ch.-2 on profile line B-B' three explanations can be put forward: (1) the depthaveraged density varies laterally across channels, (2) bridging stresses prevent a full relaxation of Ch.-2, and/or (3) Ch.-2 contains a 20-60 m thick layer of marine ice which is not detected by radar.Variations of the depth-averaged densities across channels can be caused by changes in the surface mass balance (SMB), which is a function of surface slope (e.g., Lenaerts et al., 2014).Increased SMB inside the channels can explain the observed synclines in the shallow radar layers (Fig. 5a), and accords with observations from Langley et al. (2014) for channels perpendicular to the wind direc-tion (as is the case here).Rearranging Eq. ( 2) allows us to calculate how the equivalent air content H a must change spatially in order to obtain hydrostatic equilibrium in the absence of marine ice (cf.Holland et al., 2011).This scenario necessitates H a ∼ 22 m inside Ch.-2, compared to a baseline of H a = 13-16 m outside the channel (Fig. 5d).Such a large change is only required for Ch.-2 along B-B', not for any other channel.Given that the surface depression of Ch.-2 is smaller than that of Ch.-1 (Fig. 5c), this seems unlikely.Regarding the channel width, Ch.-1 is always wider than Ch.-2, and Ch.-2 along B-B' is narrower (2.5 < α < 3.5) than Ch.-2 along A-A' (4 < α < 5).The observations cannot directly be compared with the simulations from above (due to the unknown melt history, and because of ambiguities in locating the channel keels); however, the narrowing of Ch.-2 and the subsequent mismatch between hydrostatic and radar thickness accord with what would be expected from the simulations.While other mechanisms (i.e.spatially variable density or marine ice) may still be at work, the combined evidence from modeling and observations presented here shows that bridging effects are non-negligible for narrow ice-shelf channels.Studies inverting elevation models to obtain ice thickness inside channels must take this effect into account.
The simulations presented in Fig. 6a-f exemplify characteristic points of channel formation: basal melting reduces the ice thickness inside channels and causes vertical velocities at the ice-shelf bottom to be negative.This increases the vertical strain rates (˙ zz ) inside the channels, causing subsequent lateral convergence (˙ yy ) while longitudinal stretching (˙ xx ) is only slightly affected.The spatial patterns of ˙ yy and ˙ zz mirror the prescribed channelized melting: they increase downstream as long as bottom melting is sustained, and they decrease where bottom melting is absent.Dutrieux et al. (2013) have noted lateral convergence of ice into melting channels, and this idealized simulation provides the theoretical underpinning for these observations.It becomes evident that lateral convergence increases as long as bottom melting is sustained, after bottom melting ceases, the convergence decreases.For the specific setting investigated here (i.e. for a laterally constrained ice shelf vs. a laterally unconstrained ice shelf), changes of |˙ xx | and |˙ yy | result in a local maximum of the effective strain rate inside the channel.Because ice deforms more readily when the effective strain rate is elevated (according to Glen's flow law, Cuffey and Paterson, 2010, p. 51), the simulated velocity changes step-wise inside the channel (Fig. 7a).This feature is not universal for channelized melting.It can be essentially absent when longitudinal stretching is too dominant (e.g., in the unconstrained case in which changes in | ˙ xx | offset changes in | ˙ yy | in the effective strain rate).
The simulations propose that channelized melting imprints the channel's flow field in a characteristic way which can be detected in the surface velocities.The flow of ice shelves, however, is often dominated by other mechanisms (e.g., divergence/convergence through tributary glaciers) which may mask the comparatively small effects of channelized melting.Nevertheless, along profile O-O' (Fig. 7b) a similar velocity step occurs in the observational data.This step is significant, given the coherence of the yearly and weekly averaged GPS velocities, confirming the available satellite data.The latter may partially be biased, because interferometric surface velocities require elevation models (e.g., Neckel et al., 2012), which typically do not fully resolve the channel topography.This effect, however, appears to be negligible here.Figure 8a-c shows that the cross section displayed in Fig. 7b is spatially coherent, and that the channel in this area exhibits a similar behavior as suggested by the simulations above: lateral convergence is clearly evident and accompanied by enhanced horizontal shearing.Longitudinal extension changes little.Although no direct evidence for the melt history of Ch.-1 exists, this is a likely showcase for both, i.e. how channels imprint the surrounding flow field, and how this can be measured with GPS and satellite-based data.This foreshadows a large potential of analyzing high-resolution surface velocities to determine channelized melting from space.
The basic model scenarios studied here can be extended by including merging, diverging, and meandering channels to investigate to what extent ice shelves archive temporal changes in basal hydrology at the grounding line.Because vertical velocities vary across channels, enhanced advection of cold ice from the surface should be accounted for by including thermo-mechanical coupling.From the observational side, it will be important to better understand the surface mass balance anomaly inside channels and to establish a link between the channel orientation and main wind direction.This directly impacts on the quantification of basal melt rates either using mass conservation principles or using dipping internal radar layers in an inversion procedure.The observed and simulated velocity anomaly presented here is restricted to comparatively simple strain regimes.The application of the model on a real-case geometry will guide the detection of velocity anomalies in more complex flow settings, facilitat-ing a larger-scale mapping of melting channels using surface velocities only.

Conclusions
The full Stokes modeling confirms enhanced channelized melting as a feasible mechanism for the formation of iceshelf channels.If melting is not sustained, channels gradually decay but may still persist in ice shelves for many tens of kilometers.Therefore, the mere existence of channels in satellite or radar data does not directly imply channelized melting at the location where the channel is observed.In turn, the disappearance of channels does not necessitate basal ice accretion.
The simulations show that channels can be out of hydrostatic equilibrium, and a corresponding example has been discussed for an ice-shelf channel at RBIS.Assuming that melting peaks at the channel's trough, bridging results in a hydrostatic thickness which is larger at the channel troughs, and smaller at the channel flanks.The imbalance is a function of both the channel width and the upstream melt history.This effect must be taken into account when inverting the surface elevation for ice thickness in order to obtain the basal mass balance using mass conservation.
The channel formation imprints the surrounding flow field characteristically.In areas where longitudinal stretching is not too dominant (e.g., for ice shelves which are sufficiently constrained by embayments or pinning points), this increases the effective strain rate, locally softens the ice, and produces a characteristic velocity anomaly across the channels which has been observed at RBIS in ground-and satellite-based velocities.Independent of the specific flow setting, melting channels produce generic velocity patterns on kilometer scales which are likely suited for identifying channelized melting from space, allowing to pinpoint the important role of channels on ice-ocean interactions and ice-shelf stability on large spatial scales.

Figure 1 .
Figure 1.Overview of the Roi Baudouin Ice Shelf, Dronning Maud Land, Antarctica: (a) ice-shelf channels start near the grounding line and extend into the ice shelf (marked with blue dots in the Radarsat Mosaic, Jezek and RAMP-Product-Team, 2002).The red box marks the close-up in (b) locating radar/GPS profiles of the surveys.Channels (e.g., Ch.-1, Ch.-2) appear in the background image (Landsat 8 in December 2013) as elongated lineations.

Figure 3 .
Figure 3. Colored curves display the along-flow cumulative melt rates for the different melt scenarios (MS1-3, Table1) in (a), the effective melt in (b) and the corresponding channel amplitudes in (c) for a laterally unconstrained ice shelf in steady state.

Figure 4 .
Figure 4. Steady-state results for a unconstrained ice shelf with two channels.The prescribed melt function has the same peak amplitude for both channels but varies in the across-flow width (MS4, Table 1).The deviation from hydrostatic equilibrium is shown in (a).The white curves mark zones of melting, the dashed line depicts the cross section showing bridging stresses in (b), the surface in (c), and ice thickness in (d).

Figure
Figure 5a-d compare this situation with two radar/GPS cross sections (A-A', B-B' Fig. 1) across two channels (Ch.-1,Ch.-2) with variable widths.Hydrostatic and radar thickness agree in most areas within their error bounds, except for the narrower Ch.-2 on profile B-B', where the hydrostatic thickness is anomalously larger.This anomaly is consistently observed in profiles across Ch.-2 farther upstream.

Figure 5 .
Figure 5. Observations for two channels (Ch.-1,Ch.-2) along profiles A-A' and B-B' (Fig. 1) are shown for the 400 MHz data in (a) and for the 10 MHz data in (b).Based on the surface profile (c), radar and hydrostatic thickness in (d) agree within the error bounds except for Ch.-2 along B-B' where the hydrostatic thickness is anomalously larger.The gray curve shows the required changes in the equivalent air content to obtain hydrostatic equilibrium.

Figure 6 .
Figure 6.Steady-state results for a laterally constrained ice shelf (MS5, Table 1) showing the horizontal velocities (a, b), ice thickness (c), and principal strain rates (d-f).Gray contours mark surface speed in (a), black curves in (a)-(f) delineate zones of channelized basal melting.For a cross section see Fig.7a.

Figure 8 .
Figure 8. Strain rates inferred from ERS1/2 surface velocities on a 50 m grid dating from 1996 (S.Berger, personal communication, 2014).All values are calculated in a local coordinate systems (x along flow, y across flow).The ice-shelf channel is evident due to across-flow convergence in (a) and enhanced horizontal shearing in (b), whereas along-flow extension remains inconspicuous (c).The red triangles mark positions of GPS-derived flow velocities (Fig. 7b).