2.1 Model overview and thermodynamic parameterization

Our simulations were carried out with a version of the Regional Ocean Modeling System (ROMS 3.6; Shchepetkin and McWilliams, 2009) that has been modified to include pressure, friction, and surface fluxes of heat and salt at the ice shelf base (Dinniman et al., 2007, 2011; McPhee et al., 2008; Mueller et al., 2012). ROMS is a hydrostatic, 3-D primitive equation model with a terrain-following (σ-level) coordinate system and Arakawa-C staggered grid. Our model domain (Fig. 1) covers a portion of the southern Weddell Sea, Antarctica, including FRIS. The grid spacing is 5 km with N = 24 vertical levels. A full description of model parameter choices and processing options is given in the Supplement.

Two model geometries were used in our set of simulations, one representing the modern state (standard geometry) and the other representing a possible future state (modified geometry). Model geometry consists of a land mask (including grounded ice sheet), seabed bathymetry (h), and ice draft (z_ice). These grids are described in Sect. 2.2.

Our simulations were initialized with a homogeneous, stationary ocean that has a potential temperature of either ${\mathit{\theta }}_{\mathrm{init}}=-\mathrm{1.9}$ ^∘C (cold case) or ${\mathit{\theta }}_{\mathrm{init}}=-\mathrm{1.4}$ ^∘C (warm case). Initial salinity was set as S_init= 34.65 for all cases. Model hydrography was restored to θ_init and S_init at the boundaries using a mixed radiation and nudging condition (Marchesiello et al., 2001) with a 20-day timescale. Values of θ_init and S_init for the standard geometry cold case were chosen to approximate the primary water mass entering the ice shelf cavity at present (Foldvik et al., 2001; Nicholls et al., 2001, 2009). The consequences to FRIS cavity circulation in choosing a uniform θ_init and S_init are discussed in Sect. 4. The warm case represents a moderate ocean warming scenario with an increase of 0.5 ^∘C in the temperature of water entering the FRIS cavity. This change is much smaller than the 2 ^∘C temperature increase in the inflowing water by the end of this century predicted by Hellmer et al. (2012) but was chosen to investigate whether initial feedbacks due to melt-induced changes in cavity shape from initial warming might be positive or negative. Our idealized simulations do not include wind forcing, frazil ice, or sea-ice formation.

Table 1An overview of the 12 model runs presented in this paper and the case name that is used to reference them. The three runs that were used as the spin-up solutions are referenced by their number in “restart” runs (displayed in bold font) as restart #1, #3, and #9. The runs that include passive dye tracers are flagged as “dye” runs. All run intervals provide steady state solutions as determined by the transient solutions of shelf-averaged basal melt. n/a indicates where initial temperature is not applicable.

Download Print Version | Download XLSX

Circulation develops through buoyancy forcing caused by thermodynamic exchange at the base of the ice shelves and, for tide-forced cases, by boundary conditions of tidal depth-integrated velocity and sea surface height. The thermodynamically driven component of the circulation was introduced by scalar fluxes at the ice–ocean interface beneath FRIS. These fluxes are based on a simplified version of the three-equation parameterization (Hellmer and Olbers, 1989; Holland and Jenkins, 1999) that includes the assumption that the heat flux through the ice shelf is negligible:

In Eq. (1), the surface heat flux ($Q_{\mathrm{T}}^{\mathrm{o}})$ is determined by the combined effect of thermal forcing and turbulent heat exchange. The thermal forcing is represented as $\mathrm{\Delta }T=\left(T_{\mathrm{b}}-T_{\mathrm{o}}\right)$, where T_b is the temperature at the ice–ocean interface and T_o is the temperature of the ocean mixed layer under the ice base. The value of T_b is assumed to be the freezing point temperature, T_f, and depends on the salinity at the ice–ocean interface, S_b, as well as the ice draft, z_ice < 0 (Foldvik and Kvinge, 1974; Dinniman et al., 2007). For T_o, we follow a common approach of using the temperature of the surface σ layer in place of mixed layer values, with the thickness of the surface σ layer beneath the ice shelf cavity in our standard grid ranging from 2–24 m and with 72 % of points between 5 and 15 m. The turbulent heat exchange at the ice–ocean interface is represented by a thermal transfer coefficient, α_h, scaled by a friction velocity, u_∗. This turbulent heat exchange is then adjusted by a meltwater advection term, m, that corrects the scalar fluxes for a computational drift that is introduced as an artifact of assumptions made in the numerical representation of the ice shelf as a material boundary (Jenkins et al., 2001). We define $m=-{\mathit{\alpha }}_{\mathrm{s}}{u}_{\ast }\left(\mathrm{1}-S_{\mathrm{o}}/S_{\mathrm{b}}\right)$ where S_b < 5 and m= 0 elsewhere. The friction velocity is calculated from the surface quadratic stress of the upper sigma level as $u_{\ast }=C_{\mathrm{d}}^{\mathrm{1}/\mathrm{2}}\left|u\right|$, where $C_{\mathrm{d}}=\mathrm{2.5}\times {\mathrm{10}}^{-\mathrm{3}}$ is a constant drag coefficient and $\left|u\right|$ is the magnitude of the surface-layer current. The potential density of seawater, ${\mathit{\rho }}_{\mathrm{o}}(x,y,z,t)$, is evaluated for the uppermost layer, with the heat capacity of the ocean, c_po, assumed constant at 3985 J kg^−1 ∘C⁻¹. ΔS is the salinity equivalent to ΔT and is defined as $\mathrm{\Delta }S=\left(S_{\mathrm{b}}-S_{\mathrm{o}}\right)$, with S_b solved by a quadratic formula combining Eqs. (1)–(3), without the meltwater advection term (m), and with S_o representing the salinity of the surface σ layer.

These heat and salt fluxes ($Q_{\mathrm{T}}^{\mathrm{o}},Q_{\mathrm{S}}^{\mathrm{o}})$ depend on scalar transfer coefficients (α_h, α_s) that are proportional to each other by a double diffusive parameter, $R={\mathit{\alpha }}_{\mathrm{h}}/{\mathit{\alpha }}_{\mathrm{s}}$ (McPhee et al., 2008). Chapter 2 in Mueller (2014) provides a more detailed explanation of the background and motivation for this parameterization. Here, we used scalar transfer coefficients based on observations of the RIS sub-ice shelf cavity (Jenkins et al., 2010a), with ${\mathit{\alpha }}_{\mathrm{h}}=\mathrm{1.1}\times {\mathrm{10}}^{-\mathrm{2}}$ and R= 35.5. The meltwater-equivalent melt rate term is derived by scaling the heat flux, $Q_{\mathrm{o}}^{\mathrm{T}}$, by latent heat, $L=\mathrm{3.34}\times {\mathrm{10}}^{\mathrm{5}}$ J kg⁻¹, and the density of ice, ρ_i=918 kg m⁻³, such that

Melting ice is indicated by w_b > 0 and represents the thickness of freshwater added to the ocean surface, per second. These equations highlight that basal melting is driven by ocean heat and motion, with the latter being influenced by thermohaline circulation and tides.

Figure 2Bathymetry (h) and ice draft (z_ice) for the standard and modified geometries: (a) h for both the standard and modified cases; (b) z_ice for the standard case; (c) water column thickness (wct) for the standard case; (d) difference between standard z_ice and modified z_ice, where difference > 0 indicates regions of melting and a corresponding decrease in z_ice in the modified geometry when compared to the standard geometry.

Download

A set of 12 model simulations was performed and is summarized in Table 1. A detailed description of the different simulations is given in Sect. 2.3. Each case involves a combination of standard or modified geometry, cold or warm ocean, and tidal forcing switched off or on. For tide-forced simulations, tide heights and barotropic currents were specified along the domain's open-ocean boundaries (see Fig. 1). The tidal boundary conditions were obtained for the four most energetic tidal constituents (K₁, O₁, M₂, and S₂) from the CATS2008 barotropic inverse tide model, an updated version of the circum-Antarctic model described by Padman et al. (2002). These four constituents account for 94 % of the total tidal kinetic energy for this region, based on CATS2008 estimates. Flather (1976) boundary conditions were used for the barotropic velocity with free surface conditions following Chapman (1985). Radiation conditions for baroclinic velocities were applied following Raymond and Kuo (1984). Tracer equations were radiated across the open boundaries (Marchesiello et al., 2001) and nudged to initial conditions over a 20-day timescale.

2.2 Geometries

The grid of seabed bathymetry (h) over the entire domain (Fig. 2a) was derived from the RTOPO-1 gridded dataset (Timmerman et al., 2010). The ice shelf is represented by a non-evolving, although freely floating, surface boundary based on an ice draft (z_ice, Fig. 2b) that was also derived from RTOPO-1. The land mask was adjusted around the ice rises and rumples in the southern RIS to follow the grounding line provided by Moholdt et al. (2014). Values of h and z_ice in regions of the ice sheet that are grounded in the RTOPO-1 mask but floating (i.e., ice shelf) in the mask obtained from the Moholdt et al. (2014) data were computed by linear interpolation and a nearest-neighbor extrapolation to ice shelf points in the original RTOPO-1 grids.

The ice draft and bathymetry were each smoothed to minimize errors in the baroclinic pressure gradient that arises with the terrain-following coordinate system used in ROMS (Beckmann and Haidvogel, 1993; Haney, 1991). The two parameters used to quantify smoothing are the Beckmann and Haidvogel number, $rx\mathrm{0}=\left|h\right(e)-h(e^{\prime }\left)\right|/\left(h\right(e)+h(e^{\prime })$ (Beckmann and Haidvogel, 1993), and the Haney number, $rx\mathrm{1}=\left|h\right(e,k)-h(e^{\prime },k)+h(e,k-\mathrm{1})-h(e^{\prime },k-\mathrm{1}\left)\right|/\left(h\right(e,k)+h(e^{\prime },k)-h(e,k-\mathrm{1})-h(e^{\prime },k-\mathrm{1}\left)\right)$ (Haney, 1991), where $\mathrm{1}\le k\le N$, the surface σ layer, and e and e^′ represent two horizontally adjacent grid cells. Together, these parameters establish that the surface (ice) and bottom bathymetry slopes are sufficiently small to reduce or eliminate spurious flows due to a horizontal pressure gradient and ensure hydrostatic consistency throughout the water column at adjacent horizontal grid nodes. Our Beckman and Haidvogel number, rx0, is less than 0.045 along both surface and bottom topographies, and our Haney number, rx1, is less than 10 in both surface and bottom levels except for some areas along the ice shelf front, where rx1 is larger and reaches a maximum value of 17.

Our maximum values of rx1 are larger than typically recommended for ROMS. To test whether large values lead to significant circulation from resulting errors in the baroclinic pressure gradient, we ran unforced models for each of the standard and modified grids. We initialized these models with horizontally uniform temperature and salinity fields, using an extreme density profile in order to get an upper bound on the grid errors. We chose a profile from a location just north of the eastern side of Berkner Island. The T and S profiles in this region yielded a relatively strong pycnocline with potential density changing from 1027.67 to 1027.77 kg m⁻³ between 0 and 400 m, a depth range where the change in ice draft along the ice shelf front is most likely to cause spurious flows. This profile was extrapolated in depth and horizontally to yield the uniformly stratified test hydrography. The resulting velocities in these unforced runs are representative of maximum grid errors that may occur in the full simulations. We quantified grid error by comparing the currents generated by these horizontally uniform, stratified, unforced model runs to the 30-day averaged current speeds in the standard warm tide-forced and modified warm tide-forced cases. The maximum fractional error in our velocity fields is 10 % for the standard grid and 5 % for the modified grid. Large errors were limited to a very small region north and northwest of Berkner Island; the relative error over most of the domain is negligible.

In the smoothed standard geometry, the ice draft beneath FRIS ranges from 1537 m at the deepest part of the grounding line to 11 m at the shallowest point of the ice shelf front. Small values of ice draft near the ice front are unrealistic but are a consequence of necessary smoothing in models that use terrain-following vertical coordinates. The region of thinned ice shelf represented by these small values is a narrow band along the ice front (Fig. 2b and c). The water column thickness, wct $=h+z_{\mathrm{ice}}$, ranges from 50 m (a specified minimum value, chosen for numerical stability) to 1111 m under FRIS. In the open ocean, wct =h and has a maximum value of 1211 m in the Filchner Trough.

Using this standard geometry, we conducted simulations for both the cold and warm cases described in Sect. 2.1, with and without tides. The modified geometry was then created from the output of the two 20-year tide-forced simulations of the standard cold and standard warm cases (see Sect. 2.3.2). In creating this grid, we assumed that RIS and FIS are both in steady state under present-day conditions (Rignot et al., 2013; Depoorter et al., 2013; Moholdt et al., 2014) represented by our standard cold case and that the most accurate simulations of basal melting will be those with tidal forcing included (Makinson et al., 2011). Steady state requires that mass input from lateral ice transport across the grounding line plus snowfall onto the ice shelf is balanced by basal melting and iceberg calving that maintains a constant ice-front position. The difference in local melt between the standard warm and standard cold cases, neglecting any ice dynamical feedbacks, would then be equivalent to the rate of change in thickness of the ice shelf, provided the change in melt rate is not offset by changes in mass inputs to the ice shelf.

We applied the melt rate imbalance for a period of 50 years to provide a sufficiently large change in z_ice to significantly alter the general circulation and tidal currents in the cavity. The resulting modified geometry thins z_ice by an average across the ice shelf of 30 m and a median of 14 m. The ice shelf thickens in the freezing (mid-shelf) regions by a maximum of 14 m and thins in the melt regions by a maximum of 453 m, although about 99 % of all thinning values are less than 200 m (Fig. 2c and d). The combined area where the ice shelf thickens is only 0.1 % of the total ice shelf area and is characterized by an average increase of 5 m. The modified-case bathymetry is the same as the standard case bathymetry, so these changes in z_ice cause an equal magnitude change in wct. We chose to only run the modified geometry as a warm case because it is designed to represent the FRIS cavity under warm conditions, with and without the influence of tidal forcing along the open boundaries.

2.3 Model simulations

Three types of simulations were run: (1) tide resolving with time-averaged output every 2 h and no thermodynamic exchange, (2) simulations with ice–ocean thermodynamics, with and without tide forcing, and (3) passive dye tracer simulations to explore circulation patterns. These runs are more fully described in the following sections.

2.3.1 Tide-resolving simulations with no thermodynamic exchange

We performed two 40-day simulations with 2 h averaged output, one with standard geometry and the other with modified geometry, to predict tidal current speeds. These simulations (tides-only cases) did not include thermodynamic interactions at the ice–ocean boundary, so the ocean remained unstratified at its initial homogeneous state. Absent stratification, the resulting currents are barotropic in nature, although some depth dependence arises from the friction at the seabed and ice base (see, e.g., Makinson et al., 2006).

Figure 3(a) Barotropic current ($|u{|}_{\mathrm{tide}}$, Eq. 5) for the standard tides-only run. (b) Same as (a) but for the modified tides-only run. (c) Percent change in wct between the standard and modified cases with positive values indicating where the wct is greater in the modified geometry. (d) Change in $|u{|}_{\mathrm{tide}}$ between (a) and (b) where percent change > 0 indicates locations where the standard case $|u{|}_{\mathrm{tide}}$ is greater than the modified case $|u{|}_{\mathrm{tide}}$. (e) Ratio of barotropic tidal transport ($\mathrm{wct}\times \left|u{|}_{\mathrm{tide}}\right)$ shown here as modified and standard, with values > 0 showing where there is increased transport for the modified case. (f) Seabed depth (as in Fig. 2a) with existing seismic observation locations shown as black dots. Std h refers to standard case h.

Download

The spatial characteristics of time-averaged tidal currents (Fig. 3a, b) were calculated as the time- and depth-averaged tidal current speed, $|u{|}_{\mathrm{tide}}$, given by

$$\begin{array}{}\text{(5)}& {\displaystyle }|u{|}_{\mathrm{tide}}={〈\sqrt{u_{\mathrm{b}}^{\mathrm{2}}+v_{\mathrm{b}}^{\mathrm{2}}}〉}_t[\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}],\end{array}$$

where $u_{\mathrm{b}}(x,y,t)$ and $v_{\mathrm{b}}(x,y,t)$ are orthogonal components of modeled, depth-averaged current and 〈〉_t represents temporal averaging over the last t= 30 days of the model run, which characterizes two cycles of the 15-day spring–neap cycles generated by the M₂, S₂, K₁, and O₁ tides. The maximum tidal speed at spring tides is, typically, about $\mathrm{2}\times |u{|}_{\mathrm{tide}}$.

2.3.3 Passive dye tracer simulations

Three simulations were run with passive dye tracers to investigate the advection and diffusion of water from different regions. These simulations were initialized with the steady state solutions of the standard cold, standard warm, and modified warm tide-forced cases. They were run for 2 years each, with 30-day averaged output. Two types of dyes were used. (1) Passive meltwater dyes were continuously added to the model's surface sigma layer at a rate of 1×10⁴w_b in six regions – the grounding zones of five tributary glaciers plus South Channel, which is the region of ice shelf south of Henry Ice Rise (Fig. 4). (2) A “bulk” dye was added to the open-ocean region shown in Fig. 4. The bulk dye was initialized at a concentration of 100 %, over the entire water column, but was not replenished after these simulations began.

Figure 4Locations of the six, continuous, modeled meltwater dye releases and the open-ocean bulk dye explained in Sect. 2.3.3.

Download

3.1 Tidal currents in tide-resolving simulations

The maps of $|u{|}_{\mathrm{tide}}$ defined by Eq. (5) (Fig. 3a and b) highlight the spatial variability of tidal currents beneath FRIS. In particular, they show negligible tidal currents in the inlets of the major ice streams that feed into RIS and FIS, as well as local maxima along the northeastern RIS front and South Channel.

The maximum tidal currents along the northeastern ice shelf frontal zone of RIS are consistent with previous tide models (e.g., Robertson et al., 1998; Makinson and Nicholls, 1999). This region has a relatively small wct (Fig. 2c), so a larger tidal current here is expected. The melt-induced geometry change in the modified case (Fig. 3c) has the overall effect of increasing the water depth in this region and reducing these tidal currents (Fig. 3d).

Tidal currents in South Channel are not as strong as those in the northeastern ice shelf frontal zone, but the melt-induced change in wct and $|u{|}_{\mathrm{tide}}$ is larger than in that region (Fig. 3c and d). The melt-induced change in wct is also large along the southeastern grounding line of FIS (Fig. 3c). These changes in wct and $|u{|}_{\mathrm{tide}}$ enhance barotropic tidal transport $\left(|u{|}_{\mathrm{tide}}\times \mathrm{wct}\right)$ along the southern edge of South Channel, the eastern grounding line of FIS, and along the northeast boundary of RIS with Berkner Island (Fig. 3e). These comparisons show that melt-induced ice shelf thinning generally reduces the local tidal currents (Fig. 3c, d). However, larger-scale reorganization of the barotropic tidal energy fluxes under FRIS also occurs (see, e.g., Rosier et al., 2014 and Padman et al., 2018), so simple scaling of modern tidal currents by the change in wct is not possible.

3.2 Spatial pattern of melt rates (w_b) in the base simulations

We analyzed the six base simulations described in Sect. 2.3.2 to quantify the sensitivity of the spatial pattern of basal melt rates (w_b) to θ_init, tides, and geometry.

Figure 5Melt rates averaged over 1 year of steady state solutions for (top) no-tides runs, (middle) tide-forced runs, and (bottom) the difference between the tide-forced and no-tides melt solutions. Positive values in the bottom subplots show where there is more basal melting in the solutions that include tidal forcing (middle subplots). The left column (a, d, g) shows results for the standard cold case, the middle column (b, e, h) shows results for the standard warm case, and the right column (c, f, i) shows results for the modified warm case.

Download

3.3 Sensitivity of w_b in the base simulations to tides, θ_init, and geometry

We summarize the effect of θ_init and geometry through values averaged over the ice shelf area. Net mass change (M_b: Gt a⁻¹) and averaged values of w_b were calculated for three regions: (1) all of FRIS, (2) areas for which melting conditions are predicted (w_b > 0), and (3) areas where freezing conditions are predicted (w_b < 0). The ratio of M_b and w_b for the freeze-only and melt-only calculations of mass change are not constant, since M_b is influenced by the extent of melting and freezing regions as well as by the mean magnitude of w_b.

Figure 6(a) Integrated mass flux over total FRIS area, melt-only regions, and freeze-only regions for both no-tide (⊘) and tide-forced cases. (b) Same regions and runs as in (a) but showing FRIS-averaged basal melting.

Download

Ocean temperature is the dominant control on total ice-shelf-integrated M_b and ice-shelf-averaged w_b (Fig. 6; Table 2). Regardless of tides and geometry, warming the ocean inflow by 0.5 ^∘C increases net mass loss by a factor of 3 to 5. The integrated mass gain due to freezing (marine ice formation) is insensitive to temperature for the no-tide cases but sensitive to temperature for the tide-forced cases (Fig. 6a), suggesting that accurate predictions of marine ice formation requires accurate representation of tidal currents in simulations.

Tides change the effect of ocean heat on distributions of M_b and w_b. In particular, the addition of tides to the standard cold case (our approximation to the modern state) increases net freezing by a factor of 4 – almost exactly offsetting a factor-of-2 increase in mass loss in melt-only regions. For the warm cases, tides increase total mass loss in melting regions by about 20–40 %, with most of this increase occurring in South Channel and under FIS (see Fig. 5). In contrast to the cold case, the increases in net freezing in the warm cases is small compared with the increased mass loss, so the total basal mass loss for FRIS increases significantly when tides are added to a warm ocean.

Changing geometry has a much smaller but still significant effect on ice-shelf-integrated M_b and ice-shelf-averaged w_b. In the no-tides simulations, the change from standard to modified geometry causes net mass loss to increase slightly, while the tide-forced cases show a slight decrease (Fig. 6). This change in sign in the mass loss anomaly is driven primarily by the anomalous behavior of South Channel and will be discussed in greater detail in Sect. 3.4 and 4.2.

Figure 7Melt rates averaged over the last 12 months of steady state solutions in the standard cold, standard warm, and modified warm cases for some of the regions shown in Fig. 4 and for both no-tides (⊘) and tide-forced simulations. “All FRIS” duplicates the information shown by “total” in Fig. 6b.

Download

3.4 Regional sensitivity of w_b in the base simulations

Regional averages of w_b (Fig. 7) indicate that basal mass loss near the grounding lines of major inflowing ice streams varies by an order of magnitude within a given simulation. The Support Force, Foundation, and Rutford ice streams show the largest values, exceeding 2 m a⁻¹ for the standard cold cases with and without tide forcing. In contrast, modeled melt rates for the Möller and Institute ice streams are in the range of 0.28–0.44 m a⁻¹ for the standard cold case runs.

For each ice-stream inlet, θ_init is the primary control on melt rate near the grounding line, with mean melt rate approximately doubling from the standard cold case to the standard warm case. Tidal forcing is a secondary control that leads to either an increase or decrease in w_b near grounding lines. For the Support Force, Foundation, and Rutford ice streams, adding tides reduces area-averaged melt rates, with the relative change being larger for the warm cases with both standard and modified geometry. The largest fractional change in melt rate due to tides occurs near the Rutford Ice Stream grounding line in the modified-geometry warm case, where adding tides reduces mean melt rate by 40 % from 7.7 to 4.6 m a⁻¹.

South Channel is an exception to the general result that regional sensitivity of w_b is more strongly affected by changes in θ_init than tidal forcing. In this region, w_b increases by roughly an order of magnitude between the no-tide and tide-forced simulations, whereas the fractional change due to θ_init is much smaller, about a factor of 2 (Fig. 7). South Channel also experiences a large reduction (∼ 30 %) in regionally averaged w_b in the modified warm tide-forced run compared with the standard warm tide-forced run. We attribute this change in w_b to the reduction in tidal currents as geometry is changed (Fig. 3).

3.5 Ocean circulation within the FRIS cavity in the passive dye tracer simulations

General patterns of water mass circulation into and under FRIS are demonstrated by output from the 2-year simulations with passive-tracer dyes (see Sect. 2.3.3). We focus only on tide-forced simulations because, as discussed in Sect. 3.3 and by Makinson et al. (2011), tidal currents are known to be critical to patterns of basal melting beneath FRIS.

Figure 8Distribution of the dye concentration from the last time step and the upper model layer of the standard cold, tide-forced case described in Sect. 2.3.3. These distributions are from 2 years after the initiation of the model dye releases following 20 years of model circulation spin-up time. (a) Bulk dye representing penetration of water initially north of the FRIS ice front. (b–d) Meltwater dyes with sources in the Foundation, South Channel, and Rutford regions, respectively (see Fig. 4 for dye release regions). The white line across South Channel represents the location of the transects in Fig. 13.

Download

3.5.1 Dye tracer circulation in the standard cold tide-forced case

The concentration of the open-ocean bulk dye tracer in the upper σ layer (Fig. 8a) reveals that the FRIS cavity has two different sources of heat inflow. The FIS and innermost RIS cavities are flooded by a southward transport of the open continental shelf waters across the FIS front, whereas the cavity circulation in the northeastern portion of RIS is dominated by incursions of water from across the RIS front. The latter inflow does not penetrate deep into the central and southern RIS within 2 years of simulation, although this is likely to be an artifact of the omission of the high-salinity shelf water that is known to form at the RIS front and is assumed to fuel gravity currents that reach the deep western grounding lines of RIS. In our simulations, the water entering through FIS circulates clockwise along the deep grounding line. After 2 years, some dye has reached as far west as Carlson inlet; however, very little of this dye is found under the central RIS north of the ice rises and rumples. The most direct impact of changes in open-ocean inflow is observed in the Support Force and Foundation inlets, as well as in South Channel.

Meltwater produced near the Foundation Ice Stream grounding line (Fig. 8b) reveals similar clockwise circulation, which is qualitatively similar to the plume flow described by Holland et al. (2007). This water reaches the western RIS ice front in about 2 years. Meltwater from Foundation inlet flows into all ice-stream inlets to the west of Foundation. The flow of this meltwater through South Channel is limited to the southern side of the channel.

Meltwater produced in South Channel also reaches all of the western RIS within 2 years (Fig. 8c), including much of the central region where refreezing occurs. Meltwater produced in Rutford inlet flows northward to the west of Korff Ice Rise (Fig. 8d).

These dye maps demonstrate that water found in contact with the ice shelf base in the uppermost layer, in a specific ice-stream inlet, is a mixture of the incoming high-salinity ocean water and meltwater that was produced at other inlets further upstream. As an important consequence, changes in meltwater production in different regions will alter the meltwater plume characteristics (e.g., temperature) experienced by downstream ice-stream grounding zones. In the following, it will be shown that the interactions within different grounding zone regions lead to non-local feedbacks of the melting response to changes in ocean temperatures and ice shelf geometry.

Figure 9A more expanded distribution of dye concentration than that shown in Fig. 8 to include all meltwater dyes in the three tide-forced base simulations. As in Fig. 8, dye concentrations are from the last 30-day average of upper model layer solutions from the runs described in Sect. 2.3.3. These distributions are from 2 years after initializing the dye release following 20 years of model circulation spin-up time. The left-hand column of this graphic includes the same four regional plots as shown in Fig. 8.

Download

3.5.3 Regional meltwater dye comparison for tide-forced cases

Regional meltwater dye production and advection is evaluated from the surface levels of the Foundation, Möller, and South Channel inlet regions (as in Fig. 4). Figure 10 shows the integrated values over these regions for the standard cold, standard warm, and modified warm tide-forced cases. As described in Sect. 2.3.3, the meltwater dye from a particular region is a scaled quantity of w_b that reflects the magnitude of meltwater produced in that region. The resulting passive dye is then transported through the domain through a combination of advection and mixing and acts as a proxy for the meltwater plume. In this section, we use the quantity of these meltwater tracers to demonstrate how meltwater circulation is affected by changes in θ_init and geometry.

Figure 10Integrated meltwater dye (Sect. 2.3.3) by region for the tide-forced and base simulations, showing (a) the Foundation region, (b) the Möller region, (c) the South Channel region, (d) the Institute region, and (e) all of FRIS. Regions are shown in Fig. 4.

Download

Foundation inlet shows an increase in the integrated meltwater dye with the 0.5 ^∘C increase in θ_init between the standard cold and standard warm cases (Fig. 10a). This increase in the integrated meltwater dye is not sustained with the change in cavity geometry. Instead, the net amount of the dye is reduced in the modified warm case such that the value of the integrated dye in the surface level more closely matches that of the standard cold case. This reduction in the integrated meltwater dye between the standard warm and modified warm case carries forward into the Möller region, where the reduction in the Foundation dye between the two cases is even greater than in Foundation inlet (compare Fig. 10a and b). At the same time, the reduction in the Foundation dye in the Möller region is somewhat compensated for by the Möller meltwater dye, which is consistent between the two warm cases (Fig. 10b). Overall, the Möller region appears to be less affected by the change in geometry than the Foundation region. Within South Channel, the influence of geometry on the quantity of the meltwater dye is compensated for by changes in circulation that allow for more Foundation dye in the surface level of South Channel in the modified warm case than the standard warm case (Fig. 10c). This increase in the surface level Foundation dye in South Channel is caused by changes in circulation that distribute the dye more evenly across South Channel in the modified warm case than in the standard warm case (Fig. 9).

These results highlight the fact that the regional sensitivities of the meltwater dye to θ_init and geometry may influence but not necessarily determine the quantity and temperature of meltwater in downstream regions. This result reveals the degree to which θ_init and cavity shape precondition the quantity and origin of meltwater in any given region. For example, the FRIS-integrated surface dye quantity (Fig. 10e) for Foundation and Support Force is equivalent between the standard warm and modified warm cases, even though there are strong regional variations in these cases (Fig. 10a–c). In addition, the FRIS-integrated values of meltwater dyes from the ice-front regions (RIS west, RIS east, and FIS) are similar among all cases while they differ among regions, showing greater amounts of dye in the Foundation, Möller, and South Channel regions for the warm cases than in the cold case. These regional and integrated changes demonstrate that the FRIS meltwater product is a result of regional feedbacks that are affected by a combination of production, mixing, and advection.

Table 2Values for integrated mass transport (M_b, Gt a⁻¹) and FRIS-averaged basal melt rate (w_b, m a⁻¹) for the six runs shown in Fig. 6: standard cold no tides (S −1.9), standard cold tide forced (S −1.9), standard warm no tides (S −1.4), standard warm tide forced (S −1.4), modified warm no tides (M −1.4), and modified warm tide forced (M −1.4). The symbol “⊘” is used here to denote “no tides.”

Download Print Version | Download XLSX

Table 3Publication sources and abbreviations used in Fig. 11. The methods for calculating melt rates are summarized as follows, following the nomenclature used in Table S2 of the Rignot et al. (2013) Supplement: ocean model (OM), standard glaciological method (GM), ocean observation (OO), and geophysical tracer (GT). The time period of observation(s) or forcing files is listed together with source of data or model output. n/a indicates studies in which time period is not applicable.

Download Print Version | Download XLSX

4.1 Comparison of modeled, ice-shelf-averaged w_b estimates with observations

The melt rate averaged over the area of an ice shelf is a common metric for evaluating ice shelf mass balance (e.g., Rignot et al., 2013; Depoorter et al., 2013). Our estimate of melt rate averaged over FRIS for the standard case is 0.14 m a⁻¹ freshwater equivalent, which is equal to ∼ 48 Gt a⁻¹ of net mass loss (Fig. 6 and Table 2). The range of values reported by other studies extends from 0.03 m a⁻¹, the lower bound in Depoorter et al. (2013), to 0.55 m a⁻¹ from the first oceanographically derived estimates (see Fig. 11 and Table 3; Jenkins, 1991; Jacobs et al., 1992). The range in estimates of w_b is a result of variations in observation type and model choices. Estimating w_b from observations typically requires averaging other ice shelf mass budget terms, derived from satellite observations and atmospheric models, over several years. Estimates from models are affected by the model setup.

Compared with the three most recent satellite-constrained estimates, our value is near the central estimate of 0.12 m a⁻¹ of Depoorter et al. (2013) and near the lower limit of the ranges reported by Rignot et al. (2013) and Moholdt et al. (2014). Given the model simplifications discussed above, we regard the general agreement between our integrated mass loss and prior studies as evidence that our simulations are sufficiently realistic for further sensitivity studies and interpretation of the role of tides in FRIS evolution for future climate states.

Figure 11FRIS-averaged basal melt rate comparison between this study (TS) and others. Model results are shown as black dots while observations are shown as grey dots. Error bars for remote sensing observations are shown as black lines for M14 (Moholdt et al., 2014), R13 (Rignot et al., 2013) and D13 (Depoorter et al., 2013). Min and max values are connected by thick, solid, black lines to show the range of values reported by N03 (Nicholls et al., 2003) and G94 (Gammelsrød et al., 1994). A summary of the studies presented here and their abbreviations is provided in Table 2.

Download

4.2 Sensitivity of w_b to ΔT and surface currents

We explore the regional variations of thermal forcing (ΔT) and turbulent exchange (u_∗) on w_b using the 30-day averaged values from the end of each simulation to calculate ΔT and $\left|u\right|$. Note that w_b in the 30-day averaged model output is based on the average of instantaneous heat fluxes and, therefore, includes the model's knowledge of covariances between ΔT and $\left|u\right|$ on much shorter timescales. In contrast, $\left|u\right|$ is calculated from 30-day averaged u- and v-velocity components. We use a linear combination of non-tidal and tidal currents, U, given by

where $|u{|}_{\mathrm{tide}}$ is from Eq. (5), to represent the local forcing for turbulent exchange. For the no-tides cases, $|u{|}_{\mathrm{tide}}$ is zero and $U_{\mathrm{no}\mathrm{tides}}=\left|u\right|$. We include $|u{|}_{\mathrm{tide}}$ in the tide-forced cases to more closely approximate the non-time-averaged relationship described by Eq. (1), because the 30-day average removes the tidal signal in ΔT and $\left|u\right|$ in the tide-forced cases.

Figure 12Regional influences of ΔT and current speed on melt rates. Tide-forced cases are plotted using a solid marker style, ▪, and no-tide cases are plotted using an open marker style, □. (a) Current speed (U, Eq. 6) vs. thermal forcing (ΔT), color-coded according to melt rate (w_b). Black contours follow $\mathrm{\Delta }T=c/U$ (with c being a set of different scalars), along which constant values of w_b (as in Sect. 3.2) are expected to be found. (b) ΔT vs. w_b for each region. (c) U vs. w_b for each region. (d) ΔT difference between no-tides and tide-forced cases such that positive values show where thermal forcing is stronger in the no-tides cases, (e) current speed difference between tide-forced (U_tides) and no-tides (U_no tides) cases, and (f) w_b difference between tide-forced and no-tides cases.

Download

Comparisons of the six base simulations show that w_b generally follows the expected functional dependence on ΔT and U (Fig. 12a): in all six cases, values of w_b increase with stronger currents and more thermal forcing, with values roughly falling along lines of constant ΔT⋅U. In general, our values are in range of those shown by Holland et al. (2008), (compare their Fig. 1 with our Fig. 12b), although regional differences can be seen in the bivariate relationships between w_b and either ΔT or U (Fig. 12b and c). Most ice-stream inlet averages show a similar increase in w_b with respect to ΔT (Fig. 12b), suggesting that reasonable estimates of melt rate in the ice-stream inlets could be obtained from variability of ΔT and a constant, assumed low, value of U. South Channel and, to some degree, Institute inlet diverge from this relationship, demonstrating a larger variability in w_b in relation to ΔT than is seen in other inlets (Fig. 12b). This larger variability in w_b in South Channel arises because changes in modeled melt in this area are controlled primarily by changes in U (Fig. 12b).

Comparisons of the ratios for ΔT, U, and w_b at each site between simulations without and with tides (Fig. 12d–f) show how each region responds to the combined effects of tide-induced changes in ocean conditions. With the exception of South Channel, adding tides always cools (decreases ΔT) the upper layer of ocean water adjacent to the ice base (Fig. 12d). On average, the largest reductions occur for RIS ice-stream inlets. We attribute this result to the cooling of water entering RIS inlets by inclusion of meltwater from upstream freshwater sources, with RIS inlets being influenced by rapid melting in the Support Force and Foundation inlets, as well as in South Channel (Figs. 8 and 9).

Figure 13Transects of dyes (a) and potential temperature (θ, b) across South Channel at the western tip of Henry Ice Rise. The upper four panels are for tide-forced runs only. Transect location is shown in Fig. 8.

Download

The differences between the tide-forced and no-tide cases show up more strongly in the regionally averaged comparison of U (Eq. 6, Fig. 12e). In all regions, the effect of adding tides is greater for the cold standard cases than for warm standard cases. Since the value of $|u{|}_{\mathrm{tide}}$ in Eq. (5) is the same for the standard geometry runs, this difference represents the increase in the thermohaline-driven $\left|u\right|$ from the cold to warm cases.

The largest differences in $\left|u\right|$ amongst all three model runs are in the Möller, South Channel, and Institute inlets. For the warm cases, modifying the geometry increases the ratio of $U_{\mathrm{tide}-\mathrm{forced}}/U_{\mathrm{no}\mathrm{tides}}$ for these three regions even though tidal currents decrease (Fig. 3) as wct increases. This response implies that U_no tides also declines in the modified-geometry case. A decline in U_no tides in the modified geometry is consistent with a reduction in z_ice in the inlet regions, which would reduce the thermal forcing and, hence, reduce the ice pump circulation. If true, this feedback is an artifact of our model geometry, which excludes the possibility of deeper ice that could be exposed when the grounding line migrates due to the imposed thinning. Corollary evidence for this reduction in ice pump circulation is seen in the top row of Fig. 9.

The role of South Channel melt on cooling downstream ice-stream inlets, its sensitivity to tides, and tidal sensitivity to changing z_ice suggest that reliable predictions of change in modeled w_b in the southern RIS ice streams for future climate scenarios depend on the correct representation of changes to South Channel geometry.

4.3 Role of advection through South Channel

As the maps of surface-layer dye tracers (Figs. 8 and 9) show, most water entering the FRIS cavity in our simulations flows southward under the FIS front and then circulates clockwise around the FIS and RIS grounding lines. A water parcel takes about 2 years to travel from the FIS front to the southwestern RIS region of Rutford inlet. During that time, each water parcel is subjected to mixing with meltwater, so the properties of water entering each inlet depend on the processes along the entire upstream path. Meltwater produced along each flow path continues to circulate clockwise with the inflowing warmer and saltier water that originated north of the ice front; see, e.g., dye distribution representing meltwater from Foundation inlet (Fig. 8b).

Water near the ice base in South Channel contains significant freshwater contributions from the Support Force and Foundation inlets, with a smaller contribution from Möller inlet (Fig. 9). A north–south transect across the western end of South Channel (location shown in Fig. 8 and transects in Fig. 13) shows that all water within that transect is colder than θ_init. Although not shown, this water is also fresher, i.e., some meltwater from upstream is present at all depths. Furthermore, distributions of meltwater contributions from individual regional sources vary in the vertical and horizontal (Fig. 13), with the path-integrated buoyancy fluxes determining the meltwater fractions that drive stratification. The standard geometry simulations show a core of open water dye along the bottom and northeastern slope of the trough under South Channel. The Support Force dye is concentrated near the ice base, toward the southwestern end of the transect. The Foundation dye appears in both the surface and deep model layers, concentrated on the southern side of South Channel (see also Fig. 8b). The bifurcation of the Foundation dye in South Channel reflects two sources: one in which fairly pure inflow water melts ice in Foundation inlet and then flows directly into deeper portions of South Channel and the other in which the Foundation dye flows into Möller inlet and is mixed down to the bottom of Möller inlet, which shares a similar shoaling of bathymetry (and wct) as in South Channel (Fig. 2). As expected, the South Channel dye has the highest concentration in the surface waters of this transect.

The spatial pattern in dye distribution is fairly consistent between the warm and cold standard geometry cases, although much more open water dye from north of the ice front is present in the warm case. This quantitative difference is consistent with the overall understanding that warmer θ_init drives a stronger thermohaline circulation that enhances cavity circulation and leads to a shorter residence time (Fig. 9).

More qualitative differences between simulations arise from the change in cavity shape. Except for the South Channel dye, the modified geometry shows more laterally uniform dye concentrations across the channel. Dye distribution remains vertically stratified in all three cases with the depth of the upper layer being similar in the standard warm and modified warm cases. However, even though the averaged w_b in South Channel is similar between the standard cold and modified warm cases (Fig. 7), the South Channel meltwater product does not mix down as far in the modified warm case as it does in the standard cold case (Fig. 13), which we attribute to the much weaker tidal currents in this region (Fig. 3) for modified geometry.

Transects for dye tracers are only provided for the tide-forced cases. However, a comparison of temperature transects for the no-tides and tide-forced cases (Fig. 13) shows that the thermocline is deeper in the standard geometry tide-forced cases, with the thermocline most affected in the standard cold case. As shown in Fig. 12d, tides increase the thermal forcing in South Channel in the standard geometry by a factor of 3, while having a negligible effect on ΔT in the modified geometry. These results suggest that the lowered thermocline in this region is caused by tide-induced mixing rather than advection, so its depth responds to the reduced tidal currents in the modified warm case.

4.4 Implications of regional melt on ice sheet mass balance

Walker et al. (2008) showed that the spatial distribution of ice shelf melt rates was critical to the behavior of the buttressed grounded ice streams; for the same integrated mass loss from an ice shelf, grounded ice loss was significantly faster when the melting was concentrated near the grounding line. Gagliardini et al. (2010) confirmed this analysis and also noted that a grounded ice stream could thicken and its grounding line could advance (even when net melting increased) if the melt rate decreased near the grounding line. In the context of our study, the implication is that – even when the change in the area-averaged melt rate is small – substantial variability in melt rates near ice-stream grounding lines could still have a large impact on loss (or gain) of grounded ice.

In addition to being affected by the spatial distribution of w_b, dynamic mass loss of grounded ice is also affected by bedrock slope and ice sheet topography. These factors introduce additional spatial heterogeneity in the influence of basal melting on overall mass loss from the grounded ice streams flowing into FRIS. Wright et al. (2014) used the BISICLES ice-sheet model to test the sensitivity of the grounded ice sheet to changes in FRIS mass loss at the grounding line. They found that the Institute and Möller ice streams are the most sensitive to changes in basal mass balance that might be caused by a warming ocean inflow. This result was confirmed by Martin et al. (2015) using the Parallel Ice Sheet Model (PISM). These two ice streams rest on top of steep reverse bed slopes with low basal roughness – conditions which have been shown to contribute to grounding line instability and retreat (Schoof, 2007). Furthermore, these ice streams are also sensitive to changes in the buttressing effect from ice shelf mass change around the Henry and Korff ice rises and the associated change in basal sliding over these ice rises. Our results indicate that tides currently exert a strong influence on basal mass balance in the area around the Henry and Korff ice rises (Fig. 5) by increasing melting in South Channel and increasing marine ice accretion north of the ice rises and Doake Ice Rumples.

Möller and Institute are among the lowest meltwater producing regions (Fig. 7) and receive the largest fraction of meltwater product from Foundation inlet basal melt (Fig. 10). The relative quantities of these meltwater products are sensitive to changes in advection and mixing imposed by changes in θ_init and geometry (Sect. 3.5). As shown in Fig. 7, the 0.5 ^∘C increase in θ_init increases the Möller grounding-line region w_b to 1.03 m a⁻¹ (an increase of 134 %) and Institute grounding-line region w_b to 1.30 m a⁻¹ (an increase of 100 %) for the tide-forced cases. The grounding lines in these regions appear to be very sensitive to changes in θ_init and less sensitive to changes in the cavity ocean circulation imposed by a change in model geometry. Even with Möller and Institute sensitivity to θ_init, however, these inlets are buffered from variations in open-ocean heat due to the combined influence of circulation pathways and inflowing meltwater derivatives (Fig. 10).

Of the nine grounding-line regions explored in this study, Foundation inlet has the highest averaged melt rate of 2.76 m a⁻¹ for the standard cold case (Fig. 7), a rate which more than doubles to 6.01 m a⁻¹ when θ_init increases from −1.9 to −1.4 ^∘C. However, grounded ice mass flux from Foundation Ice Stream is less sensitive to changes in basal melting than Möller and Institute (Wright et al., 2014). According to the results presented in Wright et al. (2014), even the higher melt rate with the warmed ocean in our study is insufficient to drive grounding-line retreat and significant acceleration of grounded ice loss through Foundation inlet. Therefore, it is possible that the dominant effect on the grounded ice mass budget of large w_b at Foundation is through the effect of Foundation inlet meltwater on downstream inlets, particularly Möller. As shown in the dye results presented in Sect. 3.5, Möller is somewhat isolated from FIS inflow but flooded by Foundation meltwater.

4.5 Implications of regional freeze conditions on ice sheet mass balance

As described in Sect. 3.2.1 and 3.2.2, refreezing occurs in our simulations throughout a large region of the central RIS. Refreezing in this region is qualitatively consistent with estimates of basal mass balance from satellite-based remote sensing (e.g., Joughin and Padman, 2003; Rignot et al., 2013; Moholdt et al., 2014). The extent of freezing conditions is important because marine ice accretion supports ice shelf stability through its effect on the mechanical properties of ice (Kulessa et al., 2014; McGrath et al., 2014; Li et al., 2017) and by altering the extent of grounding on topographic highs such as ice rises and rumples.

Our standard cold tide-forced case produces local maxima in marine ice growth rates in the northwestern RIS, the region northeast and east of Korff Ice Rise, and the region to the north and west of Henry Ice Rise (Fig. 5d). These regions of freezing are broadly consistent in all our model runs (Fig. 5) and the net mass increase in refreezing regions is increased when tide forcing is added (Fig. 6). Our standard cold tide-forced case has roughly 4 times more mass gain than the standard cold no-tides case, generally consistent with Makinson et al. (2011). In both warm cases, standard and modified geometry, adding tides increases net marine ice formation by a factor of 2. That is, tides will continue to be important for marine ice accretion beneath FRIS if ocean temperatures rise as predicted by Hellmer et al. (2012, 2017) and will, therefore, continue to play a role in FRIS stability.