How Accurately Should We Model Ice Shelf Melt Rates?

Assessment of ocean‐forced ice sheet loss requires that ocean models be able to represent sub‐ice shelf melt rates. However, spatial accuracy of modeled melt is not well investigated, and neither is the level of accuracy required to assess ice sheet loss. Focusing on a fast‐thinning region of West Antarctica, we calculate spatially resolved ice‐shelf melt from satellite altimetry and compare against results from an ocean model with varying representations of cavity geometry and ocean physics. Then, we use an ice‐flow model to assess the impact of the results on grounded ice. We find that a number of factors influence model‐data agreement of melt rates, with bathymetry being the leading factor; but this agreement is only important in isolated regions under the ice shelves, such as shear margins and grounding lines. To improve ice sheet forecasts, both modeling and observations of ice‐ocean interactions must be improved in these critical regions.


Introduction
In certain locations along the Antarctic coastline (Arneborg et al., 2012;Dutrieux et al., 2014;Greenbaum et al., 2015), warm Circumpolar Deep Water (CDW) exists on the continental shelf as a result of Ekman upwelling, weaker sea ice growth, and deep oceanic troughs (Jenkins et al., 2016;Petty et al., 2013;Walker et al., 2013), leading to high ice-shelf basal melt rates. In recent years, this melt has led to a large reduction in ice-shelf mass, particularly in the Amundsen Sea region (Paolo et al., 2015;Pritchard et al., 2012). This reduction lessens buttressing of the ice sheet, increasing ice sheets' contribution to sea levels (Jacobs et al., 2012;Joughin et al., 2014;Shepherd et al., 2004;Thomas, 1979).
Estimates of melt rates under Amundsen ice shelves have typically been area-averaged or area-integrated; either because estimates are based on hydrographic measurements (e.g., Jacobs et al., 2011;Miles et al., 2016;Randall-Goodwin et al., 2015;Rignot et al., 2013) or because the spacing of satellite altimetry tracks does not allow for spatially resolved measurement (Paolo et al., 2015;Pritchard et al., 2012). However, a number of studies have found spatially resolved measurements through high-resolution remote sensing methods (Berger et al., 2017;Dutrieux et al., 2013;Gourmelen et al., 2017), showing that melt rates can differ widely from their areal average at spatial scales on the order of kilometers.
Meanwhile, there has been a great deal of effort in the modeling of ice-ocean interactions in the Amundsen (e.g., Dutrieux et al., 2014;Kimura et al., 2017;Nakayama et al., 2017;Payne et al., 2007;Robertson, 2013;St-Laurent et al., 2015). While regional ocean models have been successful in reproducing ocean circulation 10.1029/2018GL080383 and its link to bulk ice-shelf melt, ice modeling suggest that the location of ice removal from an ice shelf, in addition to its bulk value, may impact its buttressing capacity (Arthern & Williams, 2017;Goldberg & Heimbach, 2013;Goldberg et al., 2012;Seroussi et al., 2017). The extent to which ocean models reproduce this spatial variability is unclear, and there is a need to strengthen the link between ocean and ice modeling if assessments of ice-sheet response to ocean forcing are to be made.
In this study, we employ a high-resolution ocean model with newly derived bathymetric data, validated against high-resolution satellite observations of melt, to better constrain the spatial variations in ice-shelf melt rates and evaluate their effect on ice-sheet stability using an adjoint-modeling approach. Focusing on Dotson and Crosson ice shelves, both situated in the Amundsen Sea and subject to strong CDW forcing, we examine the effects of different representations of bathymetry, ice-shelf draft, and physics of the ice-ocean boundary layer upon both melt rates and impact to grounded ice. We find that a number of factors are important to reproducing the observed spatial melt variability; but that capturing this variability is more important in some locations than others, at least where ice-sheet response is of interest.

Study Area
Smith, Pope, and Kohler Glaciers are three narrow interconnected ice streams in the Amundsen sector of West Antarctica, which drain into Crosson and Dotson Ice shelves. For purpose of discussion, we adopt terminology from Khazendar et al. (2016) and Gourmelen et al. (2017) and refer to them (in east-to-west order) as Pope, Smith,Kohler East,and Kohler West (Figure 3a). Although their contribution to ice flux from the continent is ∼7-8 times smaller than that of Thwaites and Pine Island Glaciers (Shepherd et al., 2002), their observed thinning rates are even larger than that of these bigger ice streams (McMillan, Shepherd, Sundal et al., 2014). They have exhibited significant grounding line retreat in recent years, with the Smith grounding line retreating at rates upward of 2 km/a . Ice-sheet modeling suggests that this retreat may have been induced by a decrease in buttressing from the Crosson and Dotson Ice Shelves (Goldberg et al., 2015), consistent with observations of increased velocities close to the grounding line of these ice streams (Lilien et al., 2018;Mouginot et al., 2014).
This drop in buttressing may be related to submarine melt-induced thinning, which can decrease buttressing (e.g., Shepherd et al., 2004). High melt rates have been observed for both Dotson and Crosson in recent years (Depoorter et al., 2013;Gourmelen et al., 2017;Lilien et al., 2018;Miles et al., 2016;Randall-Goodwin et al., 2015;Rignot et al., 2013). Between 2003 and 2008, Dotson and Crosson had net average thinning rates of 3.1 and 6.5 m/a, respectively ; and both have had strong thinning trends for the last two decades (Paolo et al., 2015).
Previously, numerical modeling of ice-ocean interactions under these ice shelves has been challenging due to inaccurate bathymetric information (Schodlok et al., 2012). A previous estimate of bathymetry, RTOPO (Timmermann et al., 2010), was constructed from a series of bathymetric soundings. However, the data set contains little information underneath Crosson and Dotson. A recent study (Millan et al., 2017) used gravity data from Operation IceBridge to generate a far more detailed bathymetric map of the region, revealing a significant cavity beneath Crosson Ice Shelf as well as a substantial oceanographic connection between Crosson and Dotson. The findings raise questions of whether models require accurate bathymetry to assess oceanographic influence on ice sheets.

Melt Rates From Remote Sensing
We generate swath elevation of Dotson and Crosson from CryoSat-2 between 2010 and 2015 (Gourmelen et al., 2018) and, to avoid interference of advecting ice-shelf topography, solve for the Lagrangian rate of surface elevation change on a 500 by 500 m grid . The Lagrangian rate of change is performed using Sentinel-1 derived velocities (McMillan, Shepherd, Gourmelen, et al., 2014). The melt rate is assessed through the following (Jenkins & Doake, 1991): where m is basal melt rate,ȧ is the surface mass balance (van Wessem et al., 2016), i is ice density of 917 kg/m 3 , w nominal ocean density of 1,028 kg/m 3 , u is ice velocity, and s is surface elevation from the digital elevation model (DEM), corrected for a 1.5-m penetration bias.

Ocean Cavity Modeling
We use the Massachusetts Institute of Technology general circulation model (MITgcm; Marshall et al., 1997) to model the circulation and melt rates underneath Dotson and Crosson Ice Shelves. The ocean model uses a stereographic polar projected grid and is restricted to a small domain ( Figure 1) which includes the ice-shelf cavities. External ocean boundary conditions are imposed from the output of a regional ocean simulation of the Amundsen Sea and shelf break . The Kimura simulation was forced by atmospheric reanalysis and agrees well with available observations and can be considered a reliable product for conditions at our domain boundaries. Monthly averages of temperature, salinity, and velocity for 2010-2014 are interpolated to our domain boundaries. The model is spun-up for 2 years with 2010 forcing. No sea-ice or ocean surface forcing is included in the model. Sub-ice shelf melt rates are calculated with a viscous sublayer model, which parameterizes turbulent fluxes of heat and salt just beneath the ice (Losch, 2008). These fluxes are determined by turbulent exchange coefficients (Holland & Jenkins, 1999). While some studies assume constant exchange coefficients (e.g., Losch, 2008;Seroussi et al., 2017), MITgcm explicitly represents their dependency on near-ice velocities (Dansereau et al., 2014). We carry out simulations with both velocity-dependent and nonvelocity dependent parameterizations. In the velocity-dependent runs, the frictional drag coefficient c D in the formulation (where u 2 * is normalized interfacial drag and U is near-ice velocity) is chosen to give area-average modeled melt similar to that of the observations for Dotson and Crosson. In the nonvelocity dependent run, the temperature exchange coefficient ( T ) is chosen to achieve the same (with S , the salt exchange coefficient, held to a fixed ratio). Experiments are summarized in Table 1, and other relevant ocean model parameters are given in Table S1 of the supporting information. Note. Melt values in parentheses indicate an alternative method of filtering the observations. The final column represents a linear estimate of VAF loss relative to the CONTROL run, calculated via equation (3). VAF = Volume Above Floatation.

Ice Sheet-Ice Shelf Modeling
We use the STREAMICE ice flow package of MITgcm (Goldberg & Heimbach, 2013) to model the response and sensitivity of Smith, Pope, and Kohler Glaciers to melt rates under Dotson and Crosson. We use it as a stand-alone model, run in the domain indicated in Figure 1a with 450-m resolution, and a fixed time step of 1 24 years. BEDMAP2 data give bathymetry and initial ice thickness. To address the lack of cavity data in BEDMAP2, we artificially deepen the bed by 50% seaward of its grounding line. While our modification of BEDMAP2 could bias against grounding line advance, the historic trend has been one of thinning and retreat. Still, this highlights the need for more reliable topographic data sets that extend over the entire continent.
In order to assess sensitivities, the model is calibrated to observations, that is, a model inversion is carried out. As described in the section S2.2, we constrain the time-evolving model, which is forced by ocean-modeled melt, to MEaSUREs (450 m) velocities (Rignot et al., 2011) as well as a record of grounded thinning rates (Gourmelen et al., 2018). Basal traction and Glen's flow law coefficient (Cuffey & Paterson, 2010) are used as controls-as in Goldberg et al. (2015), grounded ice stiffness is determined by estimating the thermal steady-state, and Glen's law coefficient is adjusted only in floating ice.
The number of control parameters is roughly 2.5 × 10 5 , so to minimize model-data misfit, an adjoint approach is used (MacAyeal, 1992). We use the Automatic Differentiation tool OpenAD (Utke et al., 2008) which allows adjoint sensitivities of STREAMICE to be generated easily in both time-independent and time-dependent modes (Goldberg et al., 2016).
Finally, calibrated parameters are used to initialize time-dependent model runs. The time-dependent adjoint model is used to assess sensitivity of grounded ice volume to melt rates over 15 years. We do not force our model with surface accumulation as we expect its low values in this region (30-40 cm/year; Arthern et al., 2006) to have minimal dynamic impact over the time scale investigated; however, such forcing would be necessary for century-scale runs.
We stress that our use of thinning observations in our calibration is not meant to reproduce evolution of the system over a specific time window; rather, it is to initialize the model in a dynamic state representative of that of Smith, Pope, and Kohler. The ice model, calibration and initialization processes, and adjoint sensitivity calculation are explained in more detail in Supporting Information S1 (Fürst et al., 2015;Goldberg, 2011;Pattyn et al., 2013).  high-melting regions are near the Smith and Pope grounding lines, as well as an elongated region south of Bear Peninsula, just east of the Dotson-Crosson shear margin. Thinning is evident in this region from the altimetry ( Figure S1).

Remotely Sensed Melt Rates
The results suggest little melt in the southeast portion of Crosson and even localized freezing. Freezing is likely an artifact of our Lagrangian tracking, since Crosson is heavily rifted in these regions, and freezing is unlikely given nearby observed ocean temperatures (Jenkins et al., 2018;Randall-Goodwin et al., 2015). Figures 2b-2d show melt rate results, averaged for each of the simulations over the years 2011-2015. Area-average melt rates (separately for each ice shelf and combined) are given in Table 1. For each model result, the average is over the region where there is circulation beneath an ice shelf. For the satellite-derived melt rates, two values are found: one in which rates are filtered between −100 and +100 ma −1 (from examination of outliers in a melt-rate distribution) and one between 0 and +100 ma −1 . The latter value assumes that the negative melt rates found are artifacts, and the ocean melt-rate parameters c D and T are based on this value.

Modeled Melt Rates
Both runs with the Millan bathymetry and velocity-dependent melt (Figures 2b and 2c) show a channelized feature along the western margin of Dotson, similar to observations. However, melt is elevated along the entire margin, in contrast to observations. It is worth noting that elevated melt is indicated by the observations along the west margin, just upstream of the grounding line protrusion. Thus, it is possible that these two "tributaries" of the channelized melt region are simply expressed in differing degrees by the model and observations.
Melt rates with the CryoSat draft (Figures 2c and 2d) have a similar pattern to observations along the western margin of Crosson, just south of Bear Peninsula. Here the mixed layer is likely guided by inverted depressions in the ice shelf ( Figure S2), while Coriolis focuses the outflow on the margin. In contrast, the topography of the Millan draft guides the flow northward (Figure 2b).
With a velocity-independent melt parameterization (Figure 2d), melt is actually decreased in the location of the channelized feature and in Crosson's west shear margin, suggesting a velocity-driven mechanism in the channel. On the other hand, there is better agreement with observations near the Pope, Smith, and Kohler East grounding lines. (All models other than the RTOPO model indicate high melt near the Kohler West grounding line.) The low melt rates near the grounding line in the velocity-dependent models are due to low velocities just beneath the shelf. This is in line with idealized models using velocity-dependent melt rates (Little et al., 2009;Snow et al., 2017), which also suggest low melting at the grounding line. The RTOPO model (see Figure S3) does indicate elevated melt rates along Dotson's west margin, but the poor agreement in every other respect is likely due to the incorrect bathymetry.
The time series of melt shows a generally decreasing trend ( Figure S4). This is in line with oceanographic estimates (Jenkins et al., 2018), although a temporary increase in 2013 is seen. As our study focuses on melt rate patterns, this is not detrimental to our aims, but care should be taken when interpreting our modeled melt rate evolution.

Grounded Ice Sensitivity to Melt Rates
Adjoint sensitivities of Volume Above Floatation (VAF; Dupont & Alley, 2005) to melt rates are calculated for Dotson and Crosson Ice Shelves (Figure 3b). Specifically, these are found with respect to a "control run" (CONTROL) forced with time-average melt from Model 1, so chosen due to the close correspondence between the Millan draft and the initial ice draft. VAF is used as it is a measure of potential contribution to sea levels; but it is not the only measure of melt rate impact on grounded ice, as discussed below.
Upon examining the adjoint sensitivities, some interesting patterns emerge. Sensitivities are seen to be small over most of Dotson, aside from the grounding line of Kohler West. Sensitivity is slightly elevated where channelized melt-driven thinning takes place, but this is still small. On Crosson, sensitivities are the largest in the vicinity of ice rumples and along the Pope, Smith, and Kohler East grounding lines. Of note, however, is the high sensitivity along the velocity shear margin of Crosson where it borders Dotson and the southern edge of Bear Peninsula. We note that the results are broadly similar to those of Reese et al. (2018), who examined instantaneous velocity response of a time-independent model to ice-shelf mass removal on a coarse grid.
The calculated adjoint sensitivities can be used to generate linearized responses of VAF to different melt rate perturbations as follows. If m i is the melt rate in an ocean grid cell i, then the incremental VAF response (relative to that of the CONTROL experiment) is found by that is, a summation over all cells i, where m ref i is the melt rate from Model 1 and * m i is the sensitivity of ΔVAF to melt rate in the cell i: evaluated at m ref .
Equation (3) is evaluated for each melt field (modeled and observed), with results given in Table 1. Despite the observed melt pattern having a smaller spatial average than that of Model 1, it yields a larger VAF loss. The reason can be traced to greater melt rates near grounding lines, particularly Kohler West and Kohler East. Still, the ice-sheet impact is relatively similar among the models (aside from the RTOPO model).
It is also informative to consider the melt rate pattern of "maximal impact" from a grounded ice loss perspective-this is a melt rate perturbation which is an exact scaling of melt rate sensitivities: where n is the total cell count and M is a perturbation spatial average. Choosing M = 3 ma −1 (in line with the approximate thinning rate of both Crosson and Dotson over the past two decades; Paolo et al., 2015) leads to a linearly predicted VAF loss of 32.1 km 3 . For reference, a spatially uniform perturbation of 3 ma −1 yields predicted loss of 8.6 km 3 .
The above are linear estimates-a limitation of the adjoint approach. For instance, grounding line retreat leads to loss of backstress from basal traction and can lead to increased grounding line thickness, which cannot be detected by linearizing about a fixed trajectory. We run two additional time-dependent simulations of the same length as CONTROL: one in which melt rate is equal to (m ref + Δm max ) and one in which it is equal to (m ref + M). The former is referred to as the FOCUS run below, while the latter is referred to as CONST. The impact of the perturbations on thinning and ice speed relative to CONTROL are shown in Figures 3c-3f. FOCUS yields considerably higher grounded thinning of the ice streams (up to 70 m over the modeled period in some locations) and also increased grounded speeds (up to 220 ma −1 ), as well as considerable speedup of Crosson. The associated VAF losses in the FOCUS and CONST experiments are 41.3 and 14.0 km 3 /a, respectively. These are higher than the predicted linear responses, likely due to model nonlinearities.

Discussion
In our experiments, the ocean simulation, which gives the best agreement with observations in terms of reproducing large-scale features (Model 2), nonetheless underestimates melt in key areas such as grounding lines. The results raise questions as to the requirements of ocean cavity models to best predict future impacts of ocean forcing on Antarctica. If the most important aspect of the melt field is near the grounding line, then accurate bathymetry-which determines delivery of dense CDW-becomes crucial.
The importance of melt near the grounding line also highlights the importance of the ocean model's melt-rate parameterization. Although our velocity-independent melt model reproduces the high melt rates observed near the grounding line, this does not necessarily mean such a parameterization is the correct one to use, as it could neglect important processes, such as potential accelerated melt due to runoff (Berger et al., 2017;Smith et al., 2017) or potential ice-shelf collapse due to channelized melt . Furthermore, we do not represent tidal effects, which could potentially be important (Jourdain et al., 2019). Moreover, our analysis assumes the satellite-inferred melt rates to be "truth," but the assumption of hydrostatic floatation could lead to systemic errors, particularly within ∼5 km of the grounding line (Brunt et al., 2010). Thus, improved observations of melt rates in the vicinity of the grounding line are needed, as well as an improved representation of ocean physics in this critical region.
In our analysis, we have assumed submarine melting to be the primary driver of loss of grounded ice. However, there are other processes that can affect ice-shelf buttressing. Ice stiffness (the Glens law parameter) influences ice flow in a similar manner to thickness, and ice-shelf weakening can have a similar effect to melt-induced thinning. In fact, Lilien et al. (2018) infer weakening of the Dotson-Crosson margin from 1996-2011. Adjoint sensitivity to Glen's law parameter (not shown) has a pattern similar to that of melting, and it is possible that observed speedup of Smith, Pope, and Kohler East is due to weakening in this shear margin. Alternatively, thinning in the western shear margin of Crosson could potentially be influencing and accelerating this weakening: as an ice shelf thins in its shear margin, shear stress and strain rates increase. Larger shearing stresses might then lead to higher levels of ice damage (Borstad et al., 2016) and thus further weakening. If such a process were to continue indefinitely, it could lead to an effective separation of Crosson and Dotson ice shelves, as has been observed for Thwaites Ice Tonge and Thwaites Eastern Ice Shelf-an event which has led to a large shift in the grounded velocity of Thwaites Glacier .
The FOCUS ice model experiment leads to far more thinning and speedup than the CONST run. Still, the additional mass loss, ∼3 km 3 /a, is not large relative to the ∼21 km 3 /a currently being lost from the region. Moreover, there is little modeled grounding line retreat, despite extensive retreat observed . The lack of grounding line retreat (which would lead to additional VAF loss) may be because the nature of the experiments precludes melt under newly floating ice; other modeling studies (Arthern & Williams, 2017;Seroussi et al., 2017) suggest that melting of newly exposed shelf near the grounding line has a large impact on retreat. Additionally, the initial model ice thickness could be predisposed against retreat: BEDMAP thicknesses are much higher than initial thickness used in Goldberg et al. (2015) along most of the grounding line ( Figure S7). That study produced large grounding line retreat using the same model at the same resolution. Thus, our experiments show that melt pattern-and not just melt volume-can have an important impact on grounded ice; but other processes are required for extensive retreat.

Conclusions
By comparing high-resolution satellite-inferred observations of ice-shelf melt against ocean cavity models, we have shown that reasonable agreement can be achieved with sufficiently accurate boundary conditions such as ice-shelf draft and ocean bathymetry. However, analysis of sensitivities of an ice sheet-ice shelf model suggests this agreement may only be important in certain locations, if the aim is to model and understand ice-sheet response to ocean forcing. Equivalently, melt rate patterns can be as important as bulk melt in determining grounded ice response to melt.
For small, narrow ice shelves like Crosson and Dotson, these locations of high sensitivity to melt are likely to include those near the grounding lines and regions of high shear. Thus, it is very important that ocean models represent ice-ocean physics accurately in these critical locations. Moreover, it is important that observations of melt in these critical locations be improved-since without this, the veracity of ocean models in these locations, and hence their utility in predicting future ice-sheet response to climate variability and change, cannot be assessed.
In this work, we have utilized an adjoint model to investigate melt sensitivities. Despite its being a linear approximation of nonlinear processes, we would advocate such an approach in future investigations of ocean forcing of ice sheets, as it can identify locations where understanding of ice-ocean processes is crucial.