## Abstract

Thermodynamic flowline and plume models for the ice shelf–ocean system simplify the ice and ocean dynamics sufficiently to allow extensive exploration of parameters affecting ice-sheet stability while including key physical processes. Comparison between geophysically and laboratory-based treatments of ice–ocean interface thermodynamics shows reasonable agreement between calculated melt rates, except where steep basal slopes and relatively high ocean temperatures are present. Results are especially sensitive to the poorly known drag coefficient, highlighting the need for additional field experiments to constrain its value. These experiments also suggest that if the ice–ocean interface near the grounding line is steeper than some threshold, further steepening of the slope may drive higher entrainment that limits buoyancy, slowing the plume and reducing melting; if confirmed, this will provide a stabilizing feedback on ice sheets under some circumstances.

## 1. Introduction

Ice shelves, the floating extensions of the outlet glaciers and ice streams that drain inland ice sheets, cover approximately 40% of the Antarctic continental shelf (Williams et al. 1998). As the underlying ocean is effectively isolated from any atmospheric influence, circulation beneath ice shelves is primarily thermohaline in nature, driven by melting and freezing at the shelf base, with tides also contributing to vertical mixing (MacAyeal 1984). This local circulation is of global interest for two reasons. First, the outflow of cooled, freshened Ice Shelf Water (ISW) formed by sub-shelf melting contributes to the formation of Antarctic Bottom Water (AABW), a key driver of thermohaline circulation. Second, the presence of warm, dense Circumpolar Deep Water (CDW) in the Amundsen Sea has led to high melt rates (tens of meters per year), thinning and retreat of ice shelves in this region of West Antarctica, and acceleration of inland ice streams (Rignot 1998; Rignot and Jacobs 2002; Payne et al. 2004; Shepherd et al. 2004). Although melting of floating ice has little direct effect on sea level (Jenkins and Holland 2007), loss of the buttressing provided by ice shelves causes increased flux of grounded ice into the ocean, drawing down the interior ice sheet and contributing to sea level rise (e.g., Dupont and Alley 2005, 2006). Uncertainty about this process was a key factor in the most recent Intergovernmental Panel on Climate Change (IPCC) report providing projections of sea level rise “excluding future rapid dynamical changes in ice flow” (Solomon et al. 2007). While recent community modeling efforts such as Sea-level Response to Ice Sheet Evolution (SeaRISE; http://websrv.cs.umt.edu/isis/index.php/SeaRISE_Assessment; Bindschadler et al. 2013) and Ice2Sea (http://www.ice2sea.eu/) have provided more realistic bounds on glacially driven sea level rise, several studies (Walker et al. 2008; Gagliardini et al. 2010) have shown that both the magnitude and spatial distribution of basal melting strongly affect ice shelf buttressing, suggesting that ice-sheet models will require accurate oceanic forcing to project sea level rise with greater precision.

However, modeling of sub–ice shelf circulation remains subject to considerable uncertainty, largely because of the inaccessibility of these regions and the resulting scarcity of observations. Until the recent development and deployment of autonomous submersibles (Nicholls et al. 2006; Jenkins et al. 2010a, 2012), observations were limited to ship-based measurements in the open ocean near ice fronts (Jacobs et al. 1996; Smethie and Jacobs 2005; Jacobs et al. 2011) and a handful of studies in which instruments were lowered through holes drilled in ice shelves (e.g., Nicholls et al. 2009). Several key elements in the modeling of sub–ice shelf circulation are borrowed from analogous oceanic phenomena that have been far more extensively observed. Dynamical models of buoyant ISW plumes (e.g., Jenkins 1991; Holland and Feltham 2006) are based upon treatments of dense overflows down slopes in laboratory experiments (Ellison and Turner 1959, 1973) and in the ocean (e.g., Killworth 1977; Bo Pederson 1980; Stigebrandt 1987; Arneborg et al. 2007). Causing greater uncertainty, theories of thermodynamics at the ice shelf–ocean interface (e.g., Holland and Jenkins 1999) rely upon the assumption that observations of melting and vertical mixing below sea ice (McPhee 1990, 1992; McPhee et al. 1999; McPhee 2008) are applicable on larger scales. While these sea ice observations are both extensive and highly detailed, there are not yet enough sub–ice shelf measurements to confirm this hypothesis, and it is uncertain whether basal roughness and ocean currents are sufficiently similar. Recent work by Jenkins et al. (2010b) shows that melt rates observed at one site beneath the Ronne Ice Shelf are consistent with multiple parameterizations of the turbulent boundary layer, highlighting the need for further observations under a wide range of conditions to narrow the uncertainties.

In this study, we examine several uncertainties in interface thermodynamics and their likely impact on ice dynamics. We begin by assessing whether three plausible versions of the thermodynamical equations produce melt rates sufficiently different to significantly affect the ice–ocean system. Our analysis of the results suggests a significant role for the drag coefficient, motivating a sensitivity study on this parameter. We then apply a simplified, fully coupled ocean–ice shelf–ice stream model to assess the range of grounding-line retreat resulting from uncertainty in the drag coefficient. Finally, because ice shelf geometry strongly affects basal melting (Jenkins 1991; Walker and Holland 2007; Little et al. 2009), we repeat some of our fixed-shelf experiments for multiple ice shelf depth profiles, in some cases finding an interesting relationship between melt rates and basal slope.

## 2. Ocean model

### a. Ocean plume dynamics

This study uses a slightly modified version of the one-dimensional ocean plume model introduced by Jenkins (1991):

where *U*, *D*, *T*, *S*, and *ρ* indicate plume velocity, thickness, temperature, salinity, and density, respectively; *ρ*_{0} is the reference density; *g* is the acceleration due to gravity; *θ* is the slope of the ice shelf base; *C _{d}* is the drag coefficient; and

*T*,

_{w}*S*, and

_{w}*ρ*are the temperature, salinity, and density of the ambient seawater outside the plume, respectively. We modify the original momentum equation by including a geostrophic factor

_{w}*ɛ*in (2) to parameterize the otherwise neglected Coriolis force (cf. Wright and Stocker 1991; Walker et al. 2009; Parizek and Walker 2010). While this model cannot fully capture the effect of water column thickness on depth-integrated barotropic flow, studies using three-dimensional isopycnic-coordinate models with explicit mixed layers (e.g., Little et al. 2009) have shown that basal slope strongly affects melt rates, providing some justification for the plume approximation. We also modify the calculation of the entrainment rate

by setting *e*_{0} = 1.2 as in Stigebrandt (1987), rather than using constant *e*_{0} as in the original model. This parameterization allows entrainment to depend on the drag coefficient, as in Bo Pederson (1980) and Arneborg et al. (2007), while remaining simple enough for use with a reduced model.

### b. Laboratory-based interface thermodynamics

The melt rate , ice shelf basal temperature and salinity *T _{b}* and

*S*, and turbulent exchange velocities

_{b}*γ*and

_{T}*γ*are calculated by a separate set of equations for thermodynamics at the ice–ocean interface. We typically use a three-equation formulation (Hellmer and Olbers 1989; Holland and Jenkins 1999), in which heat and salt balances are calculated across a thin sub-layer assumed to be at the local freezing point:

_{S}where *λ _{i}* are coefficients in the linearization of the freezing point (Millero 1978),

*p*is the pressure at the ice shelf base,

_{B}*c*

_{pi}and

*c*

_{pw}are the specific heat capacities of ice and seawater,

*L*is the latent heat of fusion, and

_{f}*T*

_{surf}is the ice shelf surface temperature. The turbulent exchange velocities are calculated using the equations derived by Jenkins (1991) from the analyses of laboratory studies reported by Kader and Yaglom (1972, 1977):

### c. Geophysically based interface thermodynamics

Several studies based on field observations have led to alternative formulations of ice–ocean interface thermodynamics. While these methods lead to simpler sets of equations, their primary motivation is to describe interface thermodynamics in terms of variables that can be measured in a geophysical setting.

Jenkins et al. (2010b) found that melt rates observed at the base of Ronne Ice Shelf by phase-sensitive radar could be accurately calculated by a suitably tuned version of (9), but that their data could be matched equally well using a simpler parameterization in which the Stanton numbers () are constants. Noting that the complexity of expressions like (9) has not been shown necessary on geophysical scales, and that in practice this formula gives nearly constant values because of the dominance of molecular diffusion in the interface sublayer, they recommended the use of constant Stanton numbers with (6)–(8).

McPhee (1992) and McPhee et al. (1999) recommended (at least under sea ice) that (8) could be dropped by assuming the interface salinity to be equal to the plume salinity and parameterizing the rate limiting process of salt diffusion through the molecular sublayer by a suitable choice of effective heat transfer velocity, leaving the two-equation formulation

which agrees well with measurements of heat fluxes beneath sea ice of widely varying roughness. Jenkins et al. (2010b) found that this formulation could also match their observations when used with constant Stanton number (), although they cautioned that it is likely less applicable across a broader range of oceanographic conditions than three-equation thermodynamics.

## 3. Ocean experiments

Our experiments using the plume model beneath a fixed ice shelf emulate the study of Holland et al. (2008) using a three-dimensional isopycnic-coordinate ocean model. We consider shelves of lengths 275 and 550 km, with depth ranging from 600 m at the grounding line to 200 m at the ice front according to the formula

where *α _{i}* are chosen to give the desired depths at grounding line and ice front, and the exponent

*n*allows the ice shelf profile (Fig. 1) to vary from linear (

*n*= −1) to highly concave (

*n*= 4). The grounding-line depth has been reduced from 1000 m in the earlier study to avoid situations involving low melt rates (at lower ocean temperatures) and steep basal slopes, which may cause the plume to reach neutral buoyancy and separate from the ice shelf base, terminating the simulation. This shallower grounding line allows us to run experiments with deep ocean temperatures ranging from −1.8° to 2.0°C at 0.2°C intervals. We apply the temperature–salinity profiles used by Holland et al. (2008), setting the vertical salinity gradient and deep ocean temperature constant for depths greater than 200 m.

In presenting these experiments, we should emphasize that, unlike Jenkins et al. (2010b), we do not have sufficient observations to evaluate the accuracy of each thermodynamical scheme. Instead, we will assess how closely the geophysically based two- and three-equation constant Stanton number methods and the laboratory-based three-equation formulation agree. We will also examine which parameters are primarily responsible for differences in model results, with an eye toward guiding further observations.

### a. Comparing geophysically and laboratory-based thermodynamics

In our first set of experiments, all parameters are as determined by Jenkins et al. (2010b) through tuning the various methods to match their observations. The laboratory three-equation method uses *C _{d}* = 0.0062, while the geophysical methods use a somewhat higher value [

*C*= 0.0097, with = 3.1 × 10

_{d}^{−5}, = 0.0011, and () = 5.9 × 10

^{−4}]. For the flattest shelf, the 550-km linear profile, maximum and mean melt rates from the geophysical methods are 9%–11% higher than laboratory. For the steepest shelf, the 275-km nonlinear profile with

*n*= 4, the two-equation constant Stanton number method produces the lowest melt rates, while the three-equation constant Stanton method closely matches the laboratory method (Fig. 2).

In an effort to separate the effects of thermodynamical scheme from those of parameterization, we also consider experiments in which the laboratory three-equation method uses the same drag coefficient (*C _{d}* = 0.0097) as the constant Stanton number methods. These runs produce better agreement between methods for the 550-km linear ice shelf by increasing laboratory melt rates to less than 3% higher than either geophysical method. However, agreement between methods worsens for the 275-km nonlinear shelf (Fig. 2, bottom).

The differences between the preceding sets of runs can be explained by examining the three roles that the drag coefficient *C _{d}* plays in the model. First, for the laboratory three-equation method, is a factor in calculating the turbulent exchange velocities using (9), so that increasing the drag coefficient produces more vertical mixing (and thus more melting) at a given plume velocity. (For the other methods, the drag coefficient is already included as part of the constant Stanton numbers.) Second, a factor of also appears in our entrainment parameterization (5), so that at a given plume velocity, increasing the drag coefficient will increase entrainment of warmer, saltier ambient water into the plume, producing more melting. Third, the drag coefficient appears in the momentum equation (2), where a higher value will decrease plume velocity, working against vertical mixing (at the ice–ocean interface) and entrainment (at the plume–ambient water interface). It is thus not immediately clear whether increasing the drag coefficient will increase or decrease basal melting, but we will find later that the influence of drag coefficient on entrainment is the strongest of these effects.

For the 550-km linear shelf, the lower drag coefficient (*C _{d}* = 0.0062) laboratory runs produce a slightly faster-flowing plume than the (three equation) constant Stanton number runs. Despite this advantage in velocity, the lower drag coefficient leads to a smaller thermal Stanton number [mean = 9.7 × 10

^{−4}from (9), versus constant value of 1.1 × 10

^{−3}] and a smaller thermal exchange velocity

*γ*. The lower

_{T}*C*also leads to a slightly lower entrainment rate from (5), resulting in a smaller temperature contrast between the plume and the interface sublayer. As solving (6)–(8) for melt rate gives, to leading order, , the constant Stanton number runs produce approximately 10% greater melt. When the laboratory runs are repeated with

_{d}*C*= 0.0097, the higher drag coefficient more than offsets a decrease in velocity, leading to increased values of thermal exchange velocity and entrainment rate and approximately 2% greater melt than for the constant Stanton number runs.

_{d}For the 275-km nonlinear shelf, the extra velocity gained by the laboratory (*C _{d}* = 0.0062) runs compensates for the lower drag coefficient, leading to thermal exchange velocities and entrainment rates only slightly lower than those of the constant Stanton number runs. The two methods thus agree rather closely for this case. In contrast to the linear shelf case, increasing the drag coefficient for the laboratory runs worsens the agreement significantly. The plume slows slightly, but a high thermal Stanton number leads to an increased turbulent exchange velocity. With the entrainment rate for the laboratory runs actually slightly exceeding that of the constant Stanton number runs, the laboratory runs now produce greater melting than the constant Stanton number runs.

Overall, these results indicate that three-equation thermodynamics with constant Stanton numbers and the laboratory-based formulation match rather well, if calculated melt rates are of greatest interest. Two-equation thermodynamics can also provide a reasonable match, except in cases involving both high ocean temperatures and high ice shelf basal slopes (cf. Fig. 2). Differences between the two three-equation methods can likely be further reduced by suitable choices of model parameters, although determination of the drag coefficient presents its own difficulties.

It must be noted that Jenkins et al. (2010b) urge caution in the use of their relatively high values for the drag coefficient, pointing out that their data do not allow independent evaluation of both drag and turbulent transfer coefficients. While their values for Stanton numbers are well constrained, deriving a value for the drag coefficient requires the assumption that the turbulent transfer coefficient found by McPhee (1992) for sea ice is valid beneath ice shelves. Uncertainties in the measurement of temperature and velocity also contribute to uncertainty in the drag coefficient.

### b. Assessing sensitivity to drag coefficient

The difficulty of precisely determining the drag coefficient even where observations exist motivates sensitivity studies on this parameter. While we have seen that the three-equation methods are comparable, we use laboratory thermodynamics because evaluation of (9) is more straightforward than determining the proper constant Stanton numbers for each value of *C _{d}*; the extra computational expense is negligible for our one-dimensional model. We repeat the experiments already described with three commonly used values of

*C*(0.0015, 0.0025, and 0.0035) and the two values suggested by Jenkins et al. (2010b) (0.0062 and 0.0097). Both the 550-km linear and 275-km nonlinear (Fig. 3) shelf profiles show significantly higher melt rates as

_{d}*C*increases, indicating that the effects of enhanced vertical mixing and entrainment are stronger than the effect of reduced velocity. Compared to

_{d}*C*= 0.0025, which has been a standard estimate since MacAyeal (1985),

_{d}*C*= 0.0097 produces 36%–40% higher mean and maximum melt rates for the 550-km linear shelf; for the 275-km nonlinear shelf, with faster plume flow, mean and maximum melt rates are up to 46% higher. Typical values of the diffusive Stanton number nearly double (from 2.20 × 10

_{d}^{−5}to 4.33 × 10

^{−5}) and the entrainment rate increases by nearly half (from 1.32 × 10

^{−5}to 1.93 × 10

^{−5}m s

^{−1}for the 275-km nonlinear shelf at 0°C); together, these effects more than make up for the plume slowing (from 8.34 to 6.05 cm s

^{−1}for the 275-km nonlinear shelf at 0°C).

To test the relative importance of drag-dependent processes, we repeat the above experiments with the entrainment parameter *e*_{0} set to the constant value used in Jenkins (1991). When entrainment no longer depends on the drag coefficient, we find that the effect of reduced velocity overcomes increased Stanton numbers and mean melt rates drop as *C _{d}* increases. The relative change is largest for linear ice shelves, which lack the steep near-grounding-line slopes needed to establish plume velocity against higher drag coefficients. For the 275-km nonlinear ice shelf, maximum melt rates for

*C*= 0.0097 can be up to 15% higher than for

_{d}*C*= 0.0025, even though mean melt rates are lower. This effect is driven by steep near-grounding-line slopes that allow plume velocities fast enough to take advantage of the increase in Stanton numbers at high

_{d}*C*; the relative difference is largest at low temperatures, where the slope is most important in driving melting, and decreases as the ocean warms. Still, the overall result is decreasing melt as the drag coefficient is increased, the opposite of what is found in experiments with drag-dependent entrainment. We thus conclude that entrainment is the most significant of the drag-dependent processes in our model.

_{d}## 4. Coupled ice–ocean experiments

Having seen that uncertainty in the poorly known drag coefficient can lead to a broad range of calculated melt rates, we will now examine the potential impact of this uncertainty on the flow of grounded ice. This investigation will require coupling the plume model to a reduced-dimensional model of ice stream and ice shelf flow.

### a. Ice model and coupling approach

We use a version of the Dupont and Alley (2005) flowline model, which is a one-dimensional “shelfy stream” model (MacAyeal 1989). This type of model is the simplest that includes the longitudinal stresses responsible for ice shelf buttressing of inland ice flow.

The momentum equation is given by

where *h*, *u*, and *ρ _{i}* are the ice thickness, velocity, and density, respectively;

*z*is the elevation of the ice base; and

_{b}*L*is the half width of the ice stream. The effective viscosity is defined as

_{y}where *B* = *A*^{−1/3} is the ice hardness parameter, and *A* is the ice softness parameter in Glen's law. Lateral drag is parameterized by

and basal drag by

where *τ _{y}* and

*τ*are constant coefficients, and

_{b}*h*is the hydrostatically determined flotation thickness. At the downstream end (

_{f}*x*=

*L*), we impose a boundary condition consistent with hydrostatic pressure against the ice front,

_{x}while at the upstream end (*x* = 0) we set the velocity to be consistent with the imposed influx of ice. Advection of ice is determined from the continuity equation

in which we neglect surface accumulation so that the basal melt rate is the only forcing. We apply a constant flux boundary condition at the upstream end, while allowing free outflux of ice at the downstream boundary.

The model is discretized using linear finite elements, with a Petrov–Galerkin upwinding method applied to (18). Nodal spacing is 250 m, with refinement to 10 m for at least 2 km up- and downstream of the grounding line when solving (13). This relatively high resolution is necessary for smooth grounding-line migration without numerical artifacts similar to those observed by Vieli and Payne (2005), despite our use of grounding-line interpolation and partial-element basal drag as in Parizek et al. (2010). Time discretization is fully implicit, with 24 time steps per year. Coupling between ocean and ice models is handled as in Parizek and Walker (2010), with the ice shelf depth profile passed to the plume model and melt rates at ice shelf nodes returned. The full output of the plume model is also saved for analysis.

Because the plume model is steady state, it can be called on time scales consistent with ice shelf evolution, avoiding the difficulty of allowing continuously changing domain geometry and the expense of running a time-dependent model at relatively fast ocean time scales. In contrast with earlier simplified coupled models (Walker and Holland 2007; Walker et al. 2009), the time required to solve an ocean model consisting only of ordinary differential equations (ODEs) is negligible, allowing more computational resources to be dedicated to improved resolution of grounding-line migration. A similarly efficient approach is taken by Gladstone et al. (2012), who couple a flowline ice model with a box model of ocean circulation that is also a set of ODEs. In the experiments described below, the time between ocean model calls is one year.

### b. Assessing sensitivity to (ocean) drag coefficient

In our coupled experiments, the plume model uses laboratory three-equation thermodynamics for each of the five values of the drag coefficient used in section 3; as discussed there, the uncertainty in melt rates due to this parameter can be much greater than that due to differing thermodynamical schemes. To clearly assess the effects of ocean parameters on the coupled system, it is necessary to begin with an ice shelf–ice stream that is a steady state of the ice model without oceanic forcing. Attempting to use the arbitrary ice shelf profiles given by (12) would introduce potentially large ice flow transients, as seen by Parizek and Walker (2010), significantly complicating analysis. The initial ice configuration is identical to that used in Walker et al. (2008), with a roughly 99-km-long ice stream flowing up an inland-deepening bed (slope 3 × 10^{−3}) and an ice shelf occupying the remainder of the 150-km domain (Fig. 4). The original study found that applying basal melting averaging 10 m yr^{−1} would cause grounding-line retreat to a new stable position if the melt was relatively evenly distributed across the shelf, but would cause unstable retreat (i.e., flotation of all ice in the domain) if the melt was sufficiently concentrated toward the grounding line. Basal melting averaging 15 m yr^{−1} was strong enough to cause unstable retreat in all cases.

As might be expected, experiments with our coupled model at an ocean temperature of −1.8°C show only minor differences in final steady state as *C _{d}* is varied. Steady mean and maximum melt rates are 1.50 and 1.79 m yr

^{−1}for

*C*= 0.0015, increasing to 2.43 and 2.94 m yr

_{d}^{−1}for

*C*= 0.0097. While this is a difference of over 60%, melt rates in all experiments remain small enough that grounding-line retreat ranges only from 1.41 to 2.65 km.

_{d}The effect of *C _{d}* becomes slightly more pronounced as the ocean temperature is increased to −1.2°C. Steady mean and maximum melt rates now range from 4.51 and 4.74 m yr

^{−1}for

*C*= 0.0015 to 7.62 and 8.37 m yr

_{d}^{−1}for

*C*= 0.0097. Melt rates are still insufficient to cause complete flotation, but grounding-line retreat increases, ranging from 5.44 to 12.06 km.

_{d}The most dramatic effect of varying *C _{d}* is seen at an ocean temperature of −0.6°C, where the drag coefficient determines stability (Fig. 5). For experiments reaching a new steady state, final mean and maximum melt rates range from 9.04 and 9.46 m yr

^{−1}at

*C*= 0.0015 to 11.79 and 13.16 m yr

_{d}^{−1}at

*C*= 0.0035. While the latter run does have melt exceeding the 10 m yr

_{d}^{−1}threshold at which unstable retreat was possible in Walker et al. (2008), its spatial distribution is only weakly concentrated near the grounding line, resulting in a long (41 km in 1260 years) but stable retreat. Both of the larger drag coefficient values introduced by Jenkins et al. (2010b) produce unstable retreat, with

*C*= 0.0062 leading to flotation of the entire ice stream in 349 years, and

_{d}*C*= 0.0097 in 214 years. Mean and maximum melt rates at the end of these experiments are 14.18 and 17.77 m yr

_{d}^{−1}for

*C*= 0.0062, and 16.50 and 20.80 m yr

_{d}^{−1}for

*C*= 0.0097.

_{d}Finally, warming the ocean to −0.2°C is sufficient to cause unstable retreat for all our values of *C _{d}*. Final mean and maximum melt rates range from 13.52 and 15.58 m yr

^{−1}for

*C*= 0.0015 up to 26.13 and 31.38 m yr

_{d}^{−1}for

*C*= 0.0097. This doubling of melt rate produces significant differences in the rate of retreat, with complete flotation requiring as long as 554 years or as little as 95 years. We note that the apparently slight difference between

_{d}*C*= 0.0015 and 0.0025 is enough to cause retreat times to vary by more than a factor of 2.

_{d}In all experiments, melting increases rapidly in the first few decades as the ice shelf base steepens, demonstrating the feedback seen by Walker and Holland (2007). Once the ice stream has time to respond dynamically to the onset of forcing, advection of ice moderates slopes and leads to melt rates increasing more slowly (for unstable runs) or reaching a steady state (Fig. 6). As the grounding line retreats, small variations (on the order of several centimeters per year) are superimposed on the overall trend in melt rate. When the ice at a model node thins to flotation, the slope of the first few kilometers of the ice shelf base increases slightly, leading to greater entrainment and a small but sudden increase in melting. As the newly floating ice adjusts to oceanic forcing, this increase fades over decades, until the ungrounding of another node restarts the cycle or a final grounding-line position is reached. Because the variations are generally less than 1% of the final melt rate, we expect that our asynchronous coupling method has not significantly affected our results. Rather, we point out subtle variations in calculated melt rates as evidence of the sensitivity of plume models to near-grounding-line slope.

## 5. Conclusions

The broad range of reduced-model experiments conducted in this study leads to three principal conclusions regarding ice–ocean thermodynamics. First, the constant Stanton number method proposed by Jenkins et al. (2010b) and the laboratory-based transfer velocities of Holland and Jenkins (1999) produce similar melt rates across a broad range of shelf geometries and ocean temperatures when both schemes are appropriately parameterized. The two-equation method (McPhee 1992) is also in reasonable agreement when basal slopes and ocean temperatures remain relatively low. The uncertainty in melt rates resulting from the choice of thermodynamics can often be smaller than that resulting from uncertainty in a single parameter, the drag coefficient at the ice–ocean interface, which influences the entrainment of warmer ambient water into the plume. Second, grounding-line retreat in a coupled ocean–ice shelf–ice stream model is highly sensitive to the drag coefficient, again because of the influence of this parameter on entrainment. Third, as will be shown in the appendix, ice shelves with very steep near-grounding-line slopes can experience lower melt rates than flatter shelves when high entrainment in this region limits buoyancy and slows the plume.

Our use of reduced models allows exploration of a broad parameter space, but requires some caution when interpreting the results. While the plume model produces a reasonable near-shelf boundary layer that couples well with the interface thermodynamics, it also simplifies sub–ice shelf circulation by ignoring horizontal variations transverse to ice shelf flow and deep-ocean dynamics. The parameterization of entrainment, while based on observations, is simpler than the vertical mixing schemes of most ocean general circulation models, so that effects attributed here to the drag coefficient may be sensitive to multiple parameters in a more complex model. In addition, ice shelf grounding-line dynamics can be more complicated in nature than can be fully captured by a reduced-dimensional model on an idealized domain. Nevertheless, our results strongly suggest that uncertainty about the effects of vertical mixing on ice–ocean thermodynamics limits the ability of coupled models to make accurate projections of grounding-line retreat.

Given the scarcity of available data on sub–ice shelf ocean circulation, nearly any new observations would be useful for improving projections. Autonomous submersibles (e.g., Jenkins et al. 2010a) can map bathymetry and water properties, while the combination of phase-sensitive radar (Corr 2002; Jenkins et al. 2006) and thermistors deployed through boreholes (e.g., Nicholls et al. 2009) allows limited validation of parameterizations of ice–ocean thermodynamics. However, as noted by Jenkins et al. (2010b), such observations constrain only the Stanton numbers and cannot be used to determine the drag or turbulent exchange coefficients. Direct measurements of turbulent transfer at the ice shelf base will be necessary to produce a clear picture of thermodynamical processes in the ice–ocean boundary layer and fully validate model parameterizations. Such measurements and their analysis are currently in progress for Pine Island, Larsen C (Nicholls et al. 2012), and George VI ice shelves, and the results, when available, are expected to inform projections of sub–ice shelf ocean circulation and its role in forcing grounding-line retreat and sea level rise.

## Acknowledgments

RTW, BRP, and RBA were supported by NASA through Grant NNX10I04G. RTW was also supported by NASA Grant NNX12AD03A to ESSIC. BRP and RBA were also supported by NSF through Grant 0909335 and through Grant 0424589 to the Center for Remote Sensing of Ice Sheets (CReSIS). DMH was supported by NSF through Grants ARC–0806393, ANT–0732869, and ANT–1043395, and by NYU Abu Dhabi Grant G1204. We thank the editor and two anonymous reviewers for their contributions to an improved final manuscript.

### APPENDIX

#### Effect of Steep Basal Slopes near the Grounding Line

Although issues with plume separation prevent us from using the Holland et al. (2008) ice shelf profiles across a broad range of ocean temperature, we run some experiments with a deeper (1000 m) grounding line to examine the effect of relatively steep basal slopes. While high slopes near the grounding line are generally expected to enhance melting, we find cases in which melt rates for the steepest (*n* = 4) shelf are lower than melt rates for the *n* = 3 and 2 profiles. This results from near-grounding-line entrainment as well as the overall shape of the profiles, all of which have the same mean slope.

These experiments span the full range of fixed ice shelf profiles described in section 3, with laboratory thermodynamics using all five previously considered values of *C _{d}*. For the highly concave (

*n*= 4) profile, plume separation occurs within the first 50 km at lower ocean temperatures (≤0°C for 275-km length and ≤−1.2°C for 550-km length) regardless of

*C*At higher ocean temperatures, at least one of the lower-exponent 275-km-long shelves experiences greater melt than the highest-concavity shelf in all experiments. The 550-km shelves, which have lower concavity, display a slightly more complex dependence on ocean temperature and drag coefficient. For these profiles, once the ocean is warm enough to drive sufficient melting to prevent plume separation, there is a temperature range for which at least one lower-concavity shelf experiences greater melt than the highest-concavity shelf. (This range depends on the drag coefficient, extending 0.4°C warmer for

_{d}.*C*= 0.0015 than 0.0097.) When the ocean is warmed further, melting becomes sufficiently strong relative to entrainment that melt rates increase with concavity.

_{d}As an example, we analyze runs with 275-km-long shelves, *C _{d}* = 0.0025, and ocean temperature of 0.8°C (Fig. A1). For this geometry, the

*n*= 4 shelf has a maximum slope of nearly 25°, roughly 3 times greater than for

*n*= 3 (8.3°) and 10 times greater than for

*n*= 2 (2.5°). The entrainment rate (5), which depends on the sine of the basal slope, follows roughly the same proportion in the first several hundred meters of these three experiments, resulting in a much thicker

*n*= 4 plume at 1 km along the flow (6.4 m, versus 2.6 m and 0.9 m, respectively). This close to the origin of the plume, melting has had little chance to contribute to buoyancy, so the momentum equation (2) requires a rapidly thickening plume to slow down; the

*n*= 4 plume's velocity falls below that of the

*n*= 3 plume by 200 m along the flow and below that of the

*n*= 2 plume by 900 m along the flow. Because the

*n*= 4 profile remains steeper than the

*n*= 3 and 2 profiles for less than 2.5 and 5 km, respectively, the

*n*= 4 plume never fully recovers from this early deceleration, reaching a peak velocity of only 17.3 cm s

^{−1}versus 21.6 and 22.9 cm s

^{−1}for the

*n*= 3 and 2 plumes. Although the

*n*= 2 plume reaches the highest velocity, it does so farther downstream and thus shallower than the

*n*= 3 plume; the greater contrast at depth between the local freezing point and ambient ocean temperature results in the

*n*= 3 plume producing the highest maximum melt rate. The dependence of melt rate on plume velocity [through the turbulent exchange velocities (9)] then leads to the second- and third-steepest ice shelf profiles producing higher melt rates than the steepest.

Should this result hold in nature, rather than being an artifact of a particular model, it presents the possibility of a limiting mechanism for melting of a coupled ocean–ice shelf system. With an early coupled model, Walker and Holland (2007) found a positive feedback between basal melting and ice shelf basal slope. At the relatively low slopes used in their experiments, it appeared that slope–melt feedback would grow indefinitely, until halted by reaching equilibrium with ice flux into the shelf. The results of the present study suggest that in some cases the feedback could be slowed by the ice shelf base reaching a sufficiently steep slope, thus complicating the process by which the coupled system reaches equilibrium, and possibly reducing grounding-line retreat. However, we caution that this preliminary result requires confirmation by observations or by a three-dimensional ocean model capable of resolving the sub-shelf boundary layer.

## REFERENCES

*Polar Oceanography. Part A: Physical Science,*W. O. Smith Jr., Ed., Academic Press, 287–334.

*Climate Change 2007: The Physical Science Basis.*Cambridge University Press, 996 pp.