## Abstract

Bathymetric sills are important features in the ocean-filled cavities beneath a few fast-retreating ice shelves in West Antarctica and northern Greenland. The sills can be high enough to obstruct the cavity circulation and thereby modulate glacial melt rates. This study focuses on the idealized problem of diabatically driven, sill-constrained overturning circulation in a cavity. The circulation beneath fast-melting ice shelves can generally be characterized by an inflow of relatively warm dense water (with temperatures of a few degrees Celsius above the local freezing point) at depth and cold, less-dense, outflowing water, which exhibits an approximately two-layer structure in observations. We use a two-layer isopycnal hydrostatic model to study the cross-sill exchange of these waters in ice shelf cavities wide enough to be rotationally dominated. A quasigeostrophic constraint is determined for the transport imposed by the stratification. Relative to this constraint, the key parameters controlling the transport and its variability are the sill height relative to the bottom layer thickness and the strength of the friction relative to the potential vorticity (PV) gradient imposed by the sill. By varying these two key parameters, we simulate a diversity of flow phenomena. For a given meridional pressure gradient, the cross-sill transport is controlled by sill height beyond a critical threshold in the eddy-permitting, low-friction regime, while it is insensitive to friction in both the low-friction and high-friction regimes. We present theoretical ideas to explain the flow characteristics: a Stommel boundary layer for the friction-dominated regime; mean–eddy PV balances and energy conversion in the low-friction, low-sill regime; and hydraulic control in the low-friction, high-sill regime, with various estimates for transport in each of these regimes.

## 1. Introduction

In recent decades, warm ocean circulation within ice shelf cavities has led to the accelerated melting of glaciers in West Antarctica and northern Greenland (Schodlok et al. 2012; Mayer et al. 2000). This circulation contributes to glacial melting through the inward transport of warm, salty water masses and amplification of the ice–ocean heat flux through greater flow velocities and temperature gradients (Holland et al. 2008). Recent observations and modeling have also shown that as glaciers retreat past bathymetric maxima, the bathymetry has a leading-order effect on the glacial melt rates as subglacial cavities are shaped by oceanic circulations on the undersides of floating glaciers (Gudmundsson et al. 2012). This rapid retreat, postulated to be caused by the retrograde bathymetric slopes (increasing elevation in the direction of ice flow), has been observed using satellite radar interferometry to varying degrees in West Antarctic glaciers such as the Pine Island, Thwaites, Smith, and Kohler glaciers (Rignot et al. 2014). Autosub cavity-transect measurements in the Pine Island Glacier (PIG) revealed a cavity that hosts the most prominent bathymetric sill yet observed and exhibits the fastest melt rate in the region (Jenkins et al. 2010; Kimura et al. 2016). Floating glaciers in northern Greenland such as the 79 North Glacier (whose sill is external to the cavity) and Petermann Glacier have similar geometries and oceanic circulations (Cai et al. 2017; Schaffer 2017).

Previous studies have suggested that grounding lines retreating on a descending slope will ultimately lead to marine ice sheet instability, a dynamically driven, runaway retreat that can result in complete discharge of the ice stream (Schoof 2007; Gudmundsson 2013; Joughin et al. 2014). The bathymetric influence of many of the glaciers in the Amundsen Sea in West Antarctica, with the PIG as a prime example, are critical in understanding why this region is the fastest melting sector of Antarctica. Here, the melting occurs primarily at the ice–ocean interface and is predominantly forced by warm ocean circulation, as opposed to calving processes or surface melting leading to subglacial discharge (Mouginot et al. 2014; Konrad et al. 2017). The sill under the PIG is speculated to have a controlling effect on melt rates due to its modulation of intruding warm, salty Circumpolar Deep Water (CDW) into Pine Island Bay (De Rydt et al. 2014). Data and observations inside the cavity have been sparse, but diabatic processes deep inside such an ice shelf cavity lead to persistent differences between the water mass properties in the interior compared to the open ocean (Jacobs et al. 2011; Dutrieux et al. 2014).

Previous studies have investigated the dynamics using simplified models with idealized ice shelf geometries without sills (Little et al. 2008) and with sills (De Rydt et al. 2014; De Rydt and Gudmundsson 2016), as well as using more comprehensive regional ocean models with realistic bathymetry (Schodlok et al. 2012; Nakayama et al. 2014; St-Laurent et al. 2015; Seroussi et al. 2017). In these studies, the bathymetric sill was consistently found to act as a topographic barrier to the inflow of warm, salty CDW, provided that the CDW layer was sufficiently thin. These models all predict enhanced friction at the ice–ocean interface as a result of the sill, which leads to greater melt rates.

The problem of pressure-driven flow across a bathymetric obstruction has been widely studied in other oceanographic contexts, such as hydraulically controlled flows (Whitehead et al. 1974; Gill 1977). However, this characterization is also appropriate for glaciers like the PIG as bathymetrically modulated exchange flows. The establishment of hydraulic control is a potential explanation for the sill’s apparent role as an obstruction to CDW inflow mentioned in previous studies (De Rydt et al. 2014; Dutrieux et al. 2014) and is discussed in section 7.

Under an ice shelf, a simple model for the leading-order dynamics is an exchange flow forced by negative buoyancy on the open-ocean side and a positive buoyancy in the far interior. This represents the common situation in West Antarctic and northern Greenland glaciers where a denser bottom layer is thicker on the open-ocean side than in the cavity, forcing an inflow at depth via a pressure head (Dutrieux et al. 2014). The oceanic buoyancy is then increased by freshening of the water mass through a transformation at the ice–ocean boundary occurring due to subglacial melting at the ice–ocean boundary in the far interior of the cavity, which establishes an overall isopycnal tilt along the length of the cavity.

Previous studies have placed less emphasis on the dynamical mechanisms via which the sill constrains the circulation. As a result, relatively little is understood about the flow regimes that manifest in cavity circulation and the important physical parameters that define these regimes and control the cross-sill transport. Using a high-resolution, minimal model, we simulate the circulation patterns over varying flow parameters to provide a qualitative and quantitative understanding of the transport and forms of variability that can be expected in real ice shelf cavities.

The outline of the paper is as follows. In section 2, we discuss the setup of a simplified two-layer dynamical model and the details of the numerical methods used to study this problem. We also discuss the considerations for a robust posing of the idealized problem. In section 3, we discuss the reasoning for partitioning the parameter space into three regimes based on the nondimensionalized sill height and friction velocity. In section 4, we present solutions in the high-friction (HF) regime and explain this theoretically as a friction-dominated Stommel balance regime with equations for the boundary layer and a prediction for the flow structure. In section 5, we present solutions in the low-friction, low-sill (LFLS) regime and discuss the emergence of gyres and eddies using potential vorticity (PV) and energy budgets. In section 6, we present solutions in the low-friction, high-sill (LFHS) regime and discuss the phenomena of shocks, which are sharp interface gradients that appear due to wave propagation in critical flow, and layer-grounding, which occurs when the bottom layer thickness reaches zero. In section 7, we derive theoretical estimates of the cross-sill transport based on rotating hydraulic control theory with uniform PV, discuss their limitations, and compare with numerical results. We classify and discuss the analytical and numerical boundaries between the regimes as a function of friction and sill height. In section 8, we summarize our findings and discuss the implications and caveats of the study. In section 9, we conclude and discuss future research directions.

## 2. Idealized cavity flow

### a. A two-layer model

We pose the problem of sill-controlled circulation under an ice shelf as a two-layer exchange flow with bathymetry (see Fig. 1 for reference geometry). Our shallow water momentum and continuity equations for *n* = 1, 2 on an *f* plane are

where *u* is the cross-channel velocity (in the *x* direction), *υ* is the along-channel velocity (in the *y* direction), and the layer thicknesses are and , with as the elevation of the interface between the layers. The reference thicknesses of the top and bottom layer are and , and is the bathymetric height. We use a linear friction velocity *r* (m s^{−1}) that is equal for both the top and bottom layers. In numerical calculations, we control grid-scale energy and enstrophy using a thickness-weighted biharmonic eddy viscosity term , for which and , where , (Griffies and Hallberg 2000). The Montgomery potential is defined as , , where is the rigid lid surface pressure, , and .

The layer thicknesses at the north and south boundaries are linearly restored toward prescribed reference thicknesses, providing a simplistic representation of processes that occur outside of the cavity and the interior. The water mass transformation near ice shelves can be significant wherever ocean waters with temperatures above the local freezing point are in contact with the ice shelf. Melting is typically enhanced near grounding lines, which for Antarctic ice shelves are partly due to the greater depths for which the pressure-dependent freezing point is reduced to lower temperatures (Joughin et al. 2012). Also, in some warm water cavities such as the PIG, which has relatively warm (i.e., above surface freezing temperature) intrusions, the location at which the ice is in contact with the warmest water is primarily concentrated at the grounding line (Dutrieux et al. 2014). Our simplification is valid only for glaciers where the sill maximum is not located near the primary source of water mass transformation. For more general lock-exchange applications, this water mass transformation, which effectively fixes the stratification, represents any external processes that do not occur near the topographic sill as long as it leads to a relatively steady across-sill pressure head.

We model the two-layer channel flow problem purely dynamically and do not include interior diabatic mixing or thermal fluxes between layers or to the ice shelf and bottom boundary, except for a simplified representation of water mass transformation at the northern and southern boundaries. This means that although we are motivated by a flow under an ice shelf cavity, our study is not specific to ice–ocean interactions and the results also apply to channel flows in general. Because of this simplified dynamical framework, we do not make any predictions for the PIG melt rate, but a total transport can be used to constrain melt rate estimates via the water mass transformation. Instead, we focus on the circulation patterns that emerge and the bathymetric and geometric constraints on these patterns and the resulting transport. Also, the depth distribution of the water mass transformation/diabatic forcing is an important topic, but not addressed in the present study, and would be better addressed with more complete vertical resolution than in the two-layer model here.

In the posing of this problem, there is a range of possible choices for upstream and downstream boundary conditions including two extremes: a specified water mass transformation or nudging to a fixed stratification. We select the latter, since it allows the transport (*Q*) to be an emergent property of the solution (except in the high-sill limit discussed in section 6), and it also allows the isopycnal heights at the boundaries to be effectively fixed to specified values. This translates into an effective water mass transformation between the two layers at the boundaries via

with each layer *n* being restored to and [where ] at the northern and southern nudging region and , defined as 12-km-wide regions adjacent to the northern and southern boundaries. We choose the nudging time scale to be day at the domain edges, with a nudging strength that decreases linearly to zero in the interior edge of the nudging zone. By construction, specifying the boundary thickness nudging and sets the diapycnal velocities in the two layers to be equal and opposite, . We experimentally selected the nudging time scale to be large enough to minimize spurious high-frequency noise and small enough that the thicknesses at the boundaries remain close to (i.e., 90%–95% of) the nudged values. Since we choose , there is an upwelling in the northern nudged region , which drives an overturning circulation that inflows in the bottom layer and outflows in the top layer from the open-ocean (northern) boundary.

The cavity has dimensions km × 200 km × 700 m, based on an idealization of bathymetry and glacier geometry of the PIG cavity in Dutrieux et al. (2014). We choose a computational domain that is twice as long (in *y*) as the PIG cavity to cleanly separate the northern and southern boundary forcing effects from the flow dynamics over the sill. The flow in the neighborhood of the sill is essentially insensitive to further increases in *L*. The domain coordinate system is defined as = [−25:25 km, −100:100 km, 0:700 m].

In this study, we define the internal baroclinic deformation radius as

for and defined as the forced top and bottom layer thicknesses at the northern boundary. The deformation radius varies throughout the domain and for each run primarily due to differences in bathymetry and the southern boundary forcing, while the northern open-ocean boundary condition allows to be fixed as a constant for all runs. For the purpose of a scaling analysis, we make the simplifying approximation .

We allow the sill height to vary between 0 and 450 m for the various cases. We prescribe the bathymetry as a Gaussian in *y*

for a chosen sill width scale km that is fixed for all cases; is the relevant scale used to approximate the bathymetric slope . The definition of approximates the shape of the bathymetry beneath the PIG (De Rydt et al. 2014). In our tests, results were insensitive to the width of the sill for sills much wider than the deformation radius, . Also, excluding the minor influences of top topography, we use a flat-topped rigid lid for our simulations. Top topography plays a minor role compared to the controlling effects of the sill. For further discussion on this, see appendix A.

Using representative values for the PIG (Jacobs et al. 2011), the Coriolis parameter is and the densities for the two layers are , which corresponds to a reduced gravity *g*′ = 0.0027 m^{2} s^{−1} and determines the deformation radius and in Eq. (4).

### b. Numerical methods

To run simulations of the channel flow problem, we use the Back of Envelope Ocean Model (BEOM), a publicly available FORTRAN code written by Pierre St-Laurent (St-Laurent 2018). The numerical scheme is similar to the Hallberg Isopycnal Model (Hallberg and Rhines 1996) but offers a special treatment of layer-grounding (isopycnal outcropping when layer thicknesses vanish) using a Salmon layer (Salmon 2002). The Salmon layer introduces an artificial term added to the Montgomery potential that prevents numerical instability due to layer-grounding by raising the potential energy of the layer to infinity as its thickness approaches zero. PV is conserved in the continuous equations with the Salmon layer present.

We modify this code to include a rigid lid pressure solver, variable rigid lid elevation, friction against the rigid lid, and biharmonic viscosity. We assume equal top and bottom friction velocities in this study. The model uses the generalized forward–backward scheme (Shchepetkin and McWilliams 2005), which has been modified for compatibility with the rigid lid implementation and is discussed in appendix B. The boundary conditions are free-slip and closed on all four boundaries (no normal flow), with thickness nudging at the north and south boundaries, as mentioned in section 2a. The model uses energy-conserving finite differences (Sadourny 1975) for momentum and third-order upwinding (e.g., Shchepetkin and McWilliams 1998) for thickness advection on an Arakawa C grid of uniform resolution m (unless otherwise specified) chosen based on the convergence of key flow properties (see appendix A for a discussion on resolution convergence). The runs in some cases require hundreds of days of model time to fully spin up due to the generation of eddies, so we choose a run duration of 1000 days to adequately achieve statistical equilibrium, classified as having small trends (less than 5% change in the last 100 days) in the domain-integrated mean and eddy energy reservoirs of kinetic energy and available potential energy [to be defined in Eqs. (27a)–(27d)].

All time-averaged results in the following discussions are derived from 100-day averages at the end of each simulation, which is approximately one residence time scale, defined by the cavity volume divided by the exchange rate ().

### c. Geostrophic transport

With the boundary condition posing in section 2a, there is a fixed isopycnal tilt along the channel if a strong enough nudging is used. Assuming the end points of the stratification are essentially fixed by the nudging, we can make an estimate for the zonal geostrophic transport. Here zonal (cross channel) geostrophic transport is used as a reference scale for the meridional (along channel) geostrophic transport; generally all of the flow crosses the channel from west to east as it flows from north to south due to topographic steering by the sill, so the zonal geostrophic transport and meridional geostrophic transport are approximately equal. This does not mean the cross-channel and along-channel pressure gradients need always be the same, but the reference zonal transport is found to predict the net cross-sill exchange well in most of the experiments discussed in sections 4–7.

Assuming weak drag, viscosity, and tendency, the bottom layer geostrophic transport is

where BFD is the transport reduction due to bathymetric form drag and IFD is that due to interfacial form drag. Mass conservation implies that

which allows Eq. (6) to be rewritten as

where we define

The simplified transport estimate is chosen such that in the quasigeostrophic (QG) limit, for as defined in Eq. (4). The approximation in Eq. (9) is valid if and , which is a reasonable approximation for all of our cases. The QG approximation (e.g., Vallis 2006) is valid for small Rossby number (Ro) and small thickness perturbations , which is a reasonable assumption for most cases except the ones with the tallest sills or lowest friction velocities, for which BFD and IFD become nonnegligible.

## 3. Regime partitioning

For both our analytical and numerical study, we primarily examine the dynamics in a parameter regime relevant to the PIG, but our results also apply to the general class of nearly geostrophic, wide-channel flows, in which the channel width is much larger than the deformation radius , and the channel length is larger than the scale sill width . For the PIG, typical length scales are km, km, and km.

We show a reference case in Fig. 1, which illustrates the geometry of the model and the interface elevation *η* for a case with sill height = 400 m, isopycnal tilt (representing a pressure head) m, and north/upstream boundary forcing height = 550 m. This case approximately matches the properties of the PIG (Jacobs et al. 2011). The imposed stratification in this case (and all our cases) has an internal baroclinic deformation radius km. The regime classifications, which will be discussed in sections 4–6, are also directly applicable for any stratification as long as the deformation radius is less than half the width of the channel (i.e., the boundary currents on opposite sides do not interact). The majority of the following figures will show the dynamics of the bottom layer mapped onto the *x–y* plane because the effects of bathymetry are stronger in the bottom layer.

We now determine the key parameters that control the dynamics of cross-sill exchange. The relevant dimensional parameters are *W*, *L*, *H*, , , *r*, , , , *f*, , and *υ*. We choose and *υ* such that the solutions are insensitive to modest changes, that is, within a factor of 2. Furthermore, we choose *W*, *L*, *H*, , and to resemble the PIG values, but in principle, we can vary and study these parameters for more general cases. Excluding these, only three parameters remain to be independently varied in this study: friction velocity *r*, sill height , and pressure head , which we nondimensionalize as

where is an estimate for the bottom layer thickness at the sill maximum, assuming a linear isopycnal slope. We summarize these parameters in Table 1.

The chosen dimensionless is based on a bulk estimate of friction PV modification in a cross-sill Stommel boundary layer, which we discuss in section 4a. Recent microstructure friction velocity estimates from the PIG lie in the range m s^{−1} (Kimura et al. 2016). Therefore, using for velocities in the range 0.1 < *u* < 2 m s^{−1}, we estimate the linear friction velocity to be in the range from *r* = 1 × 10^{−6} to 1.4 × 10^{−4} m s^{−1} for the PIG.

To understand variations of the parameters controlling the circulation, we study the effects of varying nondimensionalized friction, sill height, and pressure head from Eqs. (10a)–(10c). We plot the mean and root-mean-square deviation (RMSD) of the transport as a function of these three variables in Fig. 2 in which each parameter is varied separately, relative to a reference state defined by m, m s^{−1}, and m (or named as with the *S* subscript denoting and the *F* subscript denoting ). This reference state was chosen such that deviations from this case qualitatively illustrate all of the sensitivities of the transport and its variability. We use a fixed dimensional friction for the sill height and pressure head parameter sweeps, since the mean transport is found to be insensitive to friction (discussed in section 7). The RMSD is measured based on the last 100 days of the run using high temporal resolution calculations of the transport (*dt* = 1/20 day).

In these results, variations in friction cause a 100% change to the variability over the studied range, while sill height causes transport to drop to 25% of QG prediction at . Mean transport varies by less than 5% and variability varies by less than 30% when we vary the pressure head. Therefore, once nondimensionalized, and capture the important changes in mean transport and variability. These results suggest a reasonable partition of the parameters into three regimes: high-friction (HF) (including both low and high sills); low-friction, low-sill (LFLS); and low-friction, high-sill (LFHS). We do not distinguish between low and high sill in the HF regime because in HF, even for very high sills, transport decreases to no less than 80% of the geostrophic estimate.

In the next three sections, we discuss each of these regimes and describe the phenomena that emerge, both qualitatively and quantitatively.

## 4. High-friction regime []

### a. PV balance

In this section, we discuss the numerical solutions in the high-friction regime. We present a solution in Fig. 3 for a high-friction case , where we define the streamfunction and PV for each layer as

There is a western boundary current (i.e., at the edge of the domain toward which topographic Rossby waves would propagate) that narrows and strengthens as it approaches the steepest part of the sill, and a subsequent broadening near the sill maximum due to smaller bathymetric gradients. At the maximum, the bottom layer flow crosses the channel since topographic beta, defined as , changes sign due to a reversal of the bottom slope. Compared to the bottom layer, the sill exhibits a much weaker influence on the top layer. Overall, this solution is representative of the high-friction regime with a circulation that is steady in time and can be interpreted via a time-mean vorticity balance resembling Stommel boundary layer theory (Stommel 1948).

The sill imposes a PV barrier that must be surmounted to allow water mass exchange. We use the PV budget to identify processes that facilitate transport across the PV gradient over bathymetry. The shallow water absolute vorticity equation (or equivalently, the thickness-weighted PV equation for steady flow) for a homogeneous fluid layer with friction, but neglecting viscosity, is

For further discussion of the thickness-weighted PV equation, see appendix C. In the high-friction limit, we can assume a steady state, small Rossby number approximation (where Ro ), , to the thickness-weighted PV equation. Along with the definition in Eq. (3), we can simplify Eq. (12) as the layerwise vorticity balance:

The first two terms represent the traditional Stommel boundary layer balance (Stommel 1948), but with the topographic in place of the planetary *β*. The third term is the torque due to lateral variations in the effective friction associated with variations in the layer thickness. The last term represents the torque imposed by diabatic stretching in the nudging regions. We show the term-by-term vorticity balance in Fig. 4 for the case (same case as Fig. 3), which illustrates the magnitude of each of the terms in Eq. (13). The first two terms are the dominant ones in the HF regime, the third term provides a small contribution to the bottom layer vorticity budget near the sill maximum, and the last term represents the torque imposed by diabatic stretching and is only important in the nudging regions.

### b. Stommel boundary layer

We now extend Stommel’s theory to make an explicit prediction for the structure of the boundary layer, and derive a suitable nondimensionalization for the friction coefficient based on boundary layer modification of PV. We make a QG approximation to Eq. (13) in the bottom layer (*n* = 2) using the assumption of small thickness perturbations , which allows the leading-order balance to be written as

This also implicitly assumes that the variations in the bottom layer thickness are due solely to bathymetry instead of the interface elevation, which is justified when the bottom slope dominates the isopycnal thickness gradient. We assume that the flow is spatially slowly varying in the along-channel direction, , and make the QG approximation, . Thus, Eq. (14) leads to the approximation

the solution of which is

where is a length scale described in further detail below. We can determine the constant using *y* invariance of the transport set by the nudging zone

where is area of the northern nudging region. This assumes all of the along-channel transport occurs in the boundary layer.

The velocity profile in Eq. (16) establishes a Stommel boundary width scale

which approximately matches the boundary structure in Fig. 4 and is shown in the bottom layer streamfunction panel in Fig. 3 as a dotted red line. This demonstrates the agreement of boundary layer width in the numerical results with this simple theory. In our runs, the minimum varies significantly from 20 m to 30 km from case (LFHS) to (HF).

### c. Scaling for the importance of bottom friction

The transition from high-friction to low-friction cases is determined by friction being in the vorticity budget in Eq. (12); when friction becomes sufficiently weak, the flow must surmount the sill’s PV gradient via advective mechanisms instead. Our solutions (further discussed in section 6) and previous literature (Pratt and Whitehead 2007) indicate that frictionless boundary currents have widths of approximately , so the transition may be expected to occur when becomes smaller than the deformation radius . However, this criterion is ambiguous because typically varies by an order of magnitude along the length of the channel (see Fig. 3), while varies much less (from to about at the sill maximum for the tallest sills). We therefore define an integral Stommel boundary width , which estimates the meridionally integrated influence of friction and motivates the nondimensionalization of in Eq. (10a). Our integral Stommel width also relaxes the assumption of QG flow used to derive .

To find when the friction-induced change in PV within the boundary layer of width from the northern boundary to the sill maximum is , we start with an integral of Eq. (12) over a friction-dominated boundary layer:

For a steady-state semigeostrophic flow for (later discussed in section 6a), with no cross-channel flow across the boundary layer, , this yields

For a boundary layer of width , the terms scale according to the dimensional estimates , , (since the meridional transport ), and . Thus, Eq. (20) yields the scaling relationship

so that the meridional PV difference due to friction scales as . To find how much friction accounts for in the total PV change from the northern boundary to the sill maximum, we find the ratio , which is when friction is dominant, where is the upstream PV, is the bottom layer PV, and is the layer thickness at the sill maximum . This allows us to define the integral Stommel width via

We define as the bottom layer thickness at the sill maximum based on a simple average of the north and south interface heights.

Therefore, the integral Stommel boundary width is

and is an estimate for the bulk frictional boundary layer width. For our runs, ranges from 20 m to 120 km from the LFLS to HF cases. Using , we can equivalently express our nondimensionalized friction parameter from Eq. (10a) as the ratio of the integral Stommel boundary layer width to the deformation radius width

The simplifying approximation is made since this allows to be defined as a fixed value for each run and is a reasonable bulk estimate of the domain. However, it should be noted that at the sill maximum of the tallest sills in our runs. We will use this definition of to define the boundary between low-friction and high-friction regimes and quantify the success of this classification in section 7.

## 5. Low-friction, low-sill regime (, )

### a. Transition to low friction: Gyres and eddies

In the HF cases, the numerical and analytical solutions predict a meridional antisymmetry about the domain center (due to a symmetric sill with equal magnitude and opposite sign in ). As we decrease friction, the flow develops meridional asymmetry about the sill and intensification of the western boundary current. These features are visible in Fig. 5, which shows snapshots of the interface elevation *η* for various values of the dimensionless friction velocity . For weak friction, the boundary layer does not widen and weaken over the sill maximum, as it does in the HF cases. The flows in this regime also exhibit significant variability in time, and generate eddies that intensify as friction is reduced.

Figure 6 shows the bottom layer relative vorticity for several sill heights . These snapshots demonstrate the effect of decreased friction and increased sill heights on the flow variability, most notably in the boundary layer. Not shown is the instantaneous vorticity over a range of , and the interface elevation and velocities over a range of . The vorticity for decreasing friction velocity is consistent with the narrowing and strengthening of the western boundary current and cross-sill separation location (also seen in Fig. 5), as well as the emergence of eddies primarily shed from the boundary current. For greater sill heights, there are greater velocity magnitudes near the sill maximum along with larger isopycnal tilts. The evidence for the sill impedance of the flow is in the lower calculated transports for greater sill heights, as shown in Fig. 2. Qualitatively, the cross-sill separation region occurs progressively closer to the sill maximum (partly observed in the boundary current separation location in Fig. 6).

We show the time-averaged bottom layer PV for various values of in Fig. 7, which filters out the transport variability associated with eddy formation. This figure illustrates the steady-state narrowing and intensification of the western boundary current, a gyre-like circulation around the sill, and the along-channel asymmetry that develop as the friction is reduced. In the LFLS cases, the PV in the boundary current is observably less modified as the flow crosses the sill compared to the HF cases. Upstream of the sill, PV is particularly well conserved within the boundary current, which suggests that flow across mean PV contours primarily occurs in the eddying region downstream of the sill.

To quantify the role of eddies in facilitating cross-sill exchange, we use the time-averaged thickness-weighted PV balance

Here, we define the thickness-weighted average and deviations from that average via

as discussed in appendix C. The mean and eddy components are separated on time scales of (100) days, which is about a residence time scale, as previously stated. The terms on the left-hand side of Eq. (25) together represent the PV advection decomposed into the tendency, mean, and eddy contributions. This is balanced by the vorticity forcing due to friction, viscosity, and vortex stretching due to nudging at the boundaries.

We plot the terms of the PV balance in Fig. 8 for the LFLS case . The dominant balance is between mean and eddy advection, particularly south of the sill where the boundary current separates from the western wall (see leftmost panel in Fig. 7), with nonnegligible viscosity and friction contributions along the western boundary. Despite the prominent eddy PV fluxes, the eddy component of the transport is consistently much smaller than the transport by the time-mean flow [the mean flow has velocities of (1 m s^{−1}) while the eddying velocity is (0.1) m s^{−1}]; it contributes up to 15% of the transport in the lowest friction cases, as seen in Fig. 2. Even though the eddies do not directly influence the cross-sill transport (whose controls are discussed in section 7), the eddies mix and homogenize PV, and the resulting eddy PV flux convergence allows the flow to cross mean PV contours in the absence of friction. There is also a nonnegligible modification of the PV by viscosity as the flow crosses the sill.

### b. Eddy energetics

To understand the mechanisms of eddy generation and their locations, we now examine the energy budget. Here, we consider kinetic energy (KE) and available potential energy (APE) reservoirs separately and split them into mean and eddy components using thickness-weighted averaging according to the following isopycnal eddy-mean decomposition,

We partition the time derivatives of these four quantities into transport and conversion terms following Aiki et al. (2016).

The barotropic and baroclinic mean-to-eddy conversion terms in isopycnal coordinates are

and correspond to the forces exerted on the mean flow by Reynolds stress convergence and isopycnal form drag convergence. These terms decelerate the mean flow by converting the mean kinetic energy to eddy kinetic energy for positive values and extract energy from the eddy kinetic energy to accelerate the mean flow for negative conversion values. These terms approach the barotropic and baroclinic analogs in the QG limit of Lorenz (1955), and quantify the degree of local horizontal shear and baroclinic instability processes that generate eddies in our model. We show these terms in Fig. 9 for the LFLS case and the LFHS case , with a fixed friction , so the effects of friction are controlled for a fair comparison between two different sill heights. Note the axes of this and the next two figures (Figs. 9–11), which we show in a zoom view of the domain.

The LFLS panels in Fig. 9 show that both the barotropic and baroclinic eddy energy sources are concentrated close to the western boundary sill maximum. The horizontally integrated energy production is approximately 40% barotropic and 60% baroclinic, with a small region of intense negative barotropic conversion just before a localized peak in baroclinic conversion. Because of the sill, the bottom layer domain-integrated dissipation [calculated from the nonconservative terms in the thickness-weighted momentum equations from Aiki et al. (2016)] accounts for ⅔ of the total energy dissipated, while the top layer accounts for the remaining ⅓, since the bottom layer is more strongly eddying and the EKE dissipation by friction is roughly proportional to the EKE itself.

## 6. Low-friction, high-sill regime [, ]

### a. Varying sill height

Increasing the sill height causes the boundary current velocities to increase, which strengthens the eddy intensity, as shown in Fig. 6. However, the bottom layer thickness decreases, so the vertically integrated transport and energy conversions generally decrease. The dynamics of the separation region of the western boundary current near the sill maximum change in structure and as seen in Fig. 2, the transport variability increases for intermediate sill height, and then decreases for the tallest sills. Figure 6 shows the bottom layer relative vorticity snapshots for varying sill height. Although it is not clearly observed in Fig. 6 why transport variability first increases, then decreases as a function of increasing sill height, the western boundary current reaches the farthest extent (south) past the sill maximum in the case and exhibits the greatest magnitude oscillations in the boundary current separation location. The tallest sill cases exhibit similar vorticity peaks (with eddies that scale with ), but compared with intermediate sill heights, the eddy energy generation is reduced. In the rightmost two panels of Fig. 9, both the barotropic and baroclinic eddy energy production rates decrease by about 20% as the sill height increases from 200 to 400 m. The conversion partition between and stays relatively constant over this range of sill heights, but is more dominant for sill heights below 150 m.

Increasing the sill height also decreases the water column thicknesses at the sill maximum substantially, and most of this decrease occurs in the bottom layer. Therefore, the minimum width of boundary currents, which scales as at the sill maximum, decreases from approximately 4 km (low sills) to 2 km (tall sills). To achieve the same amount of geostrophic transport, velocities in the boundary current must increase. This acceleration of the boundary current proceeds as we increase the sill height (decrease the bottom layer thickness) until shocks form.

### b. Shock formation

Shocks (or hydraulic jumps) are sharp interface gradients that arise when wave steepening due to nonlinear propagation occurs faster than any counteracting wave dispersion or dissipative processes (Helfrich et al. 1999). We observe the existence of standing shocks in our LFHS cases (defined here as the isopycnal steepness in a localized region of width exceeding the average isopycnal tilt by an order of magnitude). These conditions occur for 0.45 (threshold discussed in section 7a), which coincides with velocities that exceed the baroclinic gravity wave speed and lead to a reduction in transport (also discussed in section 7a).

For the highest sills 0.73, the amplitude of the standing shocks can be on the same order as the bottom layer water column thickness, leading to grounding of the interface between the two density layers. In the model, this manifests as a Salmon layer with thickness . The layer-grounded region is generally localized and occurs along the bottom layer western boundary current just downstream (or upstream for high enough sills) of the sill maximum in the highest sill cases. We show an example of this in Fig. 10.

Shocks (and, as a subset, layer-grounding) occur when the flow is critical at the sill maximum. We define criticality using the composite Froude number *G* (Pratt and Whitehead 2007) as

where

There are additional criteria for the criticality of flows with generalized PV (Stern 1974, Pratt and Whitehead 2007), which have an associated necessary, but not sufficient condition that *G* must be unity for some value in the range . However, this condition is incomplete since this only indicates that at least one, and not necessarily all, of the wave modes is arrested. Since the dynamics in the wide channel case are primarily localized to a boundary layer of width , the integral relation to determine the critical condition should be calculated within this boundary layer. Furthermore, since the relevant waves here are Kelvin waves, which are entirely arrested at criticality within this boundary layer (with a dominant contribution from the bottom layer), an appropriate measure in our case for criticality is the maximum value of *G* since the transport, waves, and shocks are localized to the boundary layer where the velocity reaches a maximum. For grounded layers, the internal wave speed changes due to the modification of the potential energy by the Salmon layer term (Salmon 2002). However, this is still a reasonable approximation for the criticality of the boundary current flowing around the grounded region.

In Fig. 11, we plot *G* in the vicinity of the sill maximum for three different values of . In these runs, the standing shocks appear along the western boundary and have widths of order with background current velocities similar to internal gravity wave propagation ~1 m s^{−1}, as shown in Fig. 10. This suggests they are Kelvin wave shocks, that is, shocks formed from Kelvin waves that are arrested and reach criticality within the boundary layer (Fedorov and Melville 1996; Hogg et al. 2011). Kelvin wave shocks have been observed in previous numerical simulations (Pratt et al. 2000) and attempts to produce them in laboratory conditions have been made (Pratt 1987).

In our LFHS cases, we find that the meridional velocities are much larger than the zonal velocities near the shock region, similar to the dynamics in the western boundary layer before reaching the shock. The velocity field near the shock is approximately semigeostrophic (not shown); that is, meridional velocities across the shock are in near-geostrophic balance (with small frictional and Salmon layer contributions), while the zonal momentum balance has a leading-order advective contribution.

Even though the boundary current mainly flows around the separated region for these wide channel cases, we find that the bathymetric form drag of the shock acts to reduce the overall transport. The influence of form drag on the transport is further discussed in section 7a.

The small transport contribution within the grounded region is due to the Salmon layer, discussed in section 2b. The Salmon layer term in the region where the bottom layer is important, otherwise the bottom layer would be grounded. The role of the modified potential energy of the Salmon layer can be estimated from the momentum budget using geostrophic balance (accounting for bottom friction). The Salmon layer term is largest for the region of separation seen in Fig. 11. Also, the integrated effect of viscosity on energy is insignificant, but can be important locally in small regions near layer-grounding, shocks, and western boundary currents for the lowest friction cases. However, the velocity fields remain smooth across the shocks and the transport within the Salmon layer is not significant.

The positive vorticity region south of the sill coincides with the location of separation of the western boundary current and elevated barotropic conversion shown in Fig. 9. The location of the shock oscillates meridionally along the western boundary with period of days and distance (not shown) accompanied by an increase in amplitude for shocks that are farther south of the sill. The amplitude of the shock and its meridional oscillation distance both increase as friction decreases. The strength of the oscillatory motion decreases for sills higher and lower than , which is consistent with the peak in RMSD transport shown in Fig. 2 and combined EKE generation in section 5b for intermediate sills.

Since a steady-state viewpoint of the existence of Kelvin wave shocks appearing under critical conditions necessarily involves hydraulic control, we discuss the theory involved and how this influences the cross-sill exchange in the next section.

## 7. Constraining cross-sill exchange

In Fig. 2, transport is shown to be insensitive to friction, while sill height significantly reduces transport for . The quasigeostrophic transport in Eq. (9) accurately predicts the transport for the LFLS and to a lesser extent, the HF cases. The small difference between our observed results and the geostrophic prediction is primarily due to the observed pressure head deviating slightly from the prescribed forcing , which is affected by the choice of nudging time scale .

Additionally, the central gyre-like recirculation strength (which varies proportionally with the western boundary current) increases for low-friction cases. The recirculation strength reaches 3 times the total transport for the lowest friction case (shown in Fig. 7). There are also weaker gyre-like recirculations that develop in the northern and southern regions of the domain in response to the central gyre and connect it to the nudged regions in Fig. 7. The existence and strength of the gyre depend on a cross-isobath PV flux (unknown a priori), which is balanced by frictional vorticity destruction. However, the overall transport is approximately determined by only the forcing strength and the sill height, which leads to the difference in recirculation and overall transport.

The time scale needed to establish the total transport is on the order of tens of days (a few recirculation time scales, which is defined by the approximate time the current takes to cross from the northern to the southern boundary), while the western boundary current/recirculation strength takes hundreds of days to achieve equilibrium (a few residence time scales; not shown). In ice shelf cavities, the boundary layer transport influences the gyre-like recirculation strength, while the overall transport exerts a more direct control on the heat flux brought to the ice sheet (ignoring mixing processes, which may be important).

### a. Hydraulic control theory

In this section, we test three different theoretical constraints for the cross-sill transport against the numerically diagnosed transports. These theoretical predictions rely on a steady-state analysis in the limit of no friction, often used for hydraulically controlled flows. Our aim is to predict the LFHS cross-sill exchange, which is important since the net cross-sill exchange most directly relates to basal melt.

The theory of hydraulically controlled oceanic flows was established on the grounds of comparing oceanic abyssal flows between deep basins separated by bathymetric ridges to similarly controlled flow states in dams and reservoirs (Whitehead et al. 1974). A key characteristic of hydraulically controlled exchange flows is that the flow dynamics can be expressed in terms of one controlling parameter (e.g., bathymetry), using appropriate approximations to derive an analytical solution for the steady-state flow. Here we use classical rotating hydraulic control theory under simplifying assumptions (i.e., one-layer zero PV and two-layer uniform PV) (Whitehead et al. 1974; Gill 1977; Pratt and Armi 1990) to understand the dynamics of idealized glacial cavity circulation. For the simplest case of one-layer flow with a free surface, assuming slow variations of the flow in the along-channel direction (Pratt and Whitehead 2007) leads to a system of three equations [geostrophy in the cross-channel direction (i.e., semigeostrophy), uniform PV, and Bernoulli equations):

The three unknowns are , and is the layer thickness infinitely far upstream (assumed to be infinite for zero PV). The Bernoulli equation is derived from the Bernoulli function being conserved along streamlines (constant) for a steady flow.

Assuming most of the flow is contained within the boundary layer current, we estimate the cross-sill exchange based on a boundary layer of fixed width . We find the prescribed transport that leads to the flow reaching criticality at the sill maximum, which may be used as an estimate for the transport necessary to establish control as a function of . The analytically predicted critical transport according to one-layer, zero-PV, rotating hydraulic control theory (Pratt and Whitehead 2007) is

where *w* is the width of the flow, yielding a simple expression that captures the first-order dynamics of critical flow.

We also derive a more comprehensive critical transport estimate for a two-layer rotating uniform PV assumption within a fixed width boundary layer and transport that occurs entirely within the lateral boundary regions. We present a derivation of the rotating two-layer uniform PV solution in appendix D.

For the highest sills, we can also predict the transport required for layer-grounding to occur, by considering a one-layer uniform-PV boundary layer of dynamic width . This provides an analytical solution for the boundary layer width as a function of bathymetric height and is discussed further in appendix E. For the high-sill cases, this predicts the transport required for layer-grounding.

We compare the geostrophic, rotating one-layer zero PV, rotating two-layer uniform PV, and dynamic width transport estimates to the numerical results in Fig. 12. This shows that even the rotating one-layer zero PV version of the theory is reasonably accurate, particularly in the LFHS regime, since the Froude number in the bottom layer is generally much larger than the top-layer Froude number for higher sills . Both the rotating two-layer uniform PV and the rotating one-layer zero PV transport estimates assume no friction, so they are well suited to the LFHS regime.

For low sills, the transport matches the quasigeostrophic scaling [Eq. (9)] and decreases when critical conditions are reached. The quasigeostrophic scaling regime applies when the sill does not strongly affect the circulation strength and is also an upper bound that slightly overestimates the transport. One factor that accounts for this overestimation is the boundary nudging, which establishes a north–south interface gradient slightly smaller than the prescribed . Even though the flow profiles in low-sill cases may be somewhat modified by the sill, the sills do not control the transport imposed by the boundary conditions.

For intermediate sill heights where , the observed transport decreases with the trends predicted by one-layer and two-layer theories. We use a spread of in the boundary layer width used to calculate the solution to the one-layer zero PV case, since the width of the boundary layer changes by a factor of ~2 with varying sill height. Given that the actual flow is unsteady and a uniform PV assumption is made, we do not expect these solutions to match the predictions perfectly. However, the predicted sill height necessary for critical flow (as defined in section 6b) is for the one-layer zero PV theory and for the two-layer uniform PV theory, which is consistent with the numerical result showing . The theoretical prediction for the transition between critical and eddying regimes is defined as the sill height for which the geostrophic transport is predicted to become critical from Eq. (31). While the theory captures the overall drop-off in transport beyond = 0.45, there are still substantial differences, such as up to a factor of 2 difference between theory and simulation, and different functional dependencies on .

An alternative way of understanding the transport reduction is that for LFHS cases, the flow becomes asymmetric across the sill and produces appreciable form drag terms, which were previously assumed to be small in Eq. (6). The bathymetric and interfacial form drags contribute to the balance of the along-channel pressure force due to the imposed horizontal stratification, and thereby reduce the zonal (cross-sill) transport. The interface has approximate rotational symmetry in the HF cases, but does not in the LFLS cases. In both these cases, the interface does not exhibit large vertical excursions unless a shock forms, which is necessary for an appreciable form drag.

We calculate the relative contribution of transport reduction by IFD and BFD. The gray, dotted line in Fig. 12 shows , for and BFD defined in Eq. (6), which demonstrates that the pressure head transport term is primarily balanced by the transport reduction due to BFD for the LFHS cases, while IFD and other nonconservative terms reduce transport by approximately 15% for all cases. A zonally averaged breakdown of the terms in Eq. (6) is also provided in Fig. 12 for a LFLS and LFHS case and shows that the pressure head transport is the dominant term for lower sill heights and both pressure head transport and BFD terms are dominant, partially canceling terms for taller sills.

### b. A regime diagram for cross-sill circulation

Based on the dynamical regimes discussed in the preceding sections, we now characterize the regimes of cross-sill flow over the entire parameter space of nondimensionalized friction and sill height parameter space in Fig. 13, which includes representative snapshots of the instantaneous bottom layer thickness. This diagram summarizes the three regimes that have been discussed in sections 4–6 and the analytical predictions and numerical results that define the regime boundaries.

We empirically classify each individual simulation, which is represented as a point on the regime diagram, as critical if the maximum composite Froude number satisfies max() = max() ≥ 1. We choose this boundary based on the success of hydraulic control theory in predicting the transition between geostrophically constrained and hydraulically controlled transport in section 7a. The analytical boundary between critical and eddying regimes (red dotted line in Fig. 13) is defined based on the theory in section 7a, which is consistent with the numerical classification based on the maximum composite Froude number.

Similarly, we classify individual simulations as eddying or friction-dominated based on the thickness-weighted PV Eq. (25). When the ratio of the friction to eddy terms spatially integrated and time averaged over the last 100 days of the run exceeds 1 {}, we classify the point as friction-dominated based on our findings in section 4. We define the theoretical boundary between friction and eddying regimes as = 1/5 (blue dotted line), which we choose empirically. The theoretical and numerical classifications for the friction and eddying regimes agree approximately, but exhibits some disagreement at lower sill heights. This is due to limitations of Stommel theory for intermediate friction, which is plausible since the theory in section 4 primarily applies for asymptotically large . The constant value of 1/5 validates the assumptions made in deriving the scaling argument in Eq. (23) since is constant over varying sill height and .

We note that there is no explicit prediction for the boundary between the critical regime and the friction-dominated regime since friction does not have a strong effect on the transport, as discussed in section 3. However, we have extended the = 1/5 boundary line between the critical and friction-dominated regimes based on experimental evidence and the definition for , which takes into account the varying sill height.

## 8. Summary and discussion

In this study, we model the dynamics of ocean circulation in sill-impeded subshelf cavities, motivated by the fast-melting PIG (Jenkins et al. 2010; Dutrieux et al. 2014). To investigate the dynamics, we simplified the problem to a two-layer diabatically forced circulation with a simple channel geometry overlying a bathymetric sill. Our numerical solutions exhibit a variety of phenomena, including Stommel boundary layers, gyre-like circulations, eddies, hydraulic shocks, and layer-grounding.

This study demonstrates that friction and sill height have important controlling effects on the flow structure and variability of the flow in a simple buoyancy-forced flow. Increasing the friction does not strongly influence the transport, but decreases the variability. The transport under ice shelves is approximately geostrophic and generally controlled by pressure head and local water mass transformations due to the interaction with the ice base, except when the sill height begins to impede the flow. In the case of the PIG, this occurs for pycnocline depths m and pressure head m, which are comparable to values observed in the PIG [ = 500–600 m, = 100–200 m, based on Dutrieux et al. (2014)]. The transport and variability as a function of friction, sill height, and pressure head suggests a division of the parameter space into three separate regimes: HF, LFLS, and LFHS.

The theory for the HF regime provides an analytical solution for the flow profile and western boundary current over the sill, which narrows with larger , and exhibits sill crossing at the sill maximum due to a sign change in .

In the LFLS regime, barotropic/baroclinic eddies are abundant in the cavity, and reshape the flow into gyre-like circulations with strong inertial western boundary currents. We showed in section 5a that the inertial boundary layer and eddy PV fluxes, rather than friction, facilitate flow across the sill-imposed PV gradient. The central gyre-like recirculation does not reach the northern and southern nudged boundaries, but may still be important for diabatic mixing processes in the interior especially over the sill maximum (which may indirectly contribute to the heat transfer into and out of the domain and locally increase the heat exchange at the ice base).

In the LFHS regime discussed in section 6, standing shocks emerge in the numerical solutions when the flow reaches criticality, and either an increase in sill height or a shock amplitude increase due to low friction can lead to layer-grounding. The composite Froude number provides justification of the shock as an arrested Kelvin wave for the lowest friction cases.

The HF and LFLS cases generally exhibit a constant transport predicted by QG, while in the LFHS regime, the transport is sill-constrained and is reduced in accordance with hydraulic control theory. A regime diagram (Fig. 13) summarizes the HF, LFLS, and LFHS regimes. The dynamics of various glaciers in Antarctica and Greenland may be studied with their position in the regime diagram in mind. Specifically, the latest measurements in the PIG suggest that the sill-controlled circulation is in the LFHS critical regime, based on microstructure estimates that suggest , as briefly discussed in section 3 (Kimura et al. 2016). The results from the bottom circulation in a high-resolution model of the PIG cavity (Dutrieux et al. 2014) and a process model with additional dynamics (De Rydt et al. 2014) are also qualitatively similar to our results from low-to-intermediate friction and those with quadratic friction ( appendix A). All of these show an inflow that crosses the sill from west to east close to the sill maximum and exhibit similar gyre recirculations.

For ice shelf cavities, an additional way for flow to cross mean PV contours in the absence of eddies and bottom friction is lateral friction against the sidewalls (e.g., Little et al. 2008), but typically the depth tapers to zero at lateral boundaries in an ice shelf cavity (Kimura et al. 2016) and the top/bottom friction becomes dominant instead. The importance of lateral friction for each scenario can be formally calculated by finding its associated boundary length scale and comparing with and as discussed in section 4c. Also, the evolution of glaciers as they melt may generally allow the behavior of ocean circulation to change across different regimes. For hydraulically controlled circulations, the ice shelf cavity potentially evolves in a way to favor the development of controlled states. These states may be preferred since they provide the greatest exchange of warm water masses carrying heat to the ice shelf for a specified geometry. The existence of bathymetry acts to reduce the transport, so it is generally less than the geostrophic transport.

Given the transport and , this also has implications for the ice–ocean heat flux if the transport is geostrophic. Based on this fixed transport, the total water mass transformation in the interior and ice–ocean heat flux may be calculated using the water mass properties of the exiting shelf water. Using the geostrophic transport as an estimate of ocean circulation in ice shelf cavities, we can calculate the residence time of the water masses as a function of the pressure head to be 3 months for the PIG, based on Sv from Jenkins et al. (2010). If the time-averaged stratification at the ice shelf cavity entrance is small enough to allow the sill to constrain the transport, our hydraulic control estimate predicts a transport of ~0.24 Sv for the choice of forcing m used in our study [based on Jacobs et al. (2011)], which yields months. These estimates do not take recirculations into account (see Borenas and Whitehead 1998), which may be significant for low-friction cases. It is also likely that in the presence of large external variability, the total transport is more likely to be able to respond to external variability. See Holland (2017) for a discussion of the importance of circulation residence times on glacial melt rates.

Although this study is primarily motivated by the hydraulically controlled sill under an ice shelf cavity in the PIG, the regimes observed in Fig. 13 can also be generalized to many more sill-influenced exchange flows. Specifically, the top-layer friction does not strongly influence the regime diagram, since this layer is weakly modified by the sill and exhibits much weaker velocities (the composite Froude number is dominated by the bottom layer unless the top layer is much shallower). Therefore, these results can be directly applicable to scenarios with a free surface in open-ocean, wide-sill overflows with similar stratifications such as the Denmark Strait and the Faroe–Bank Channel (Pratt and Whitehead 2007), among many others.

Given the highly idealized nature of the study, there are numerous caveats associated with our results. In addition to geometric simplifications and quadratic friction (discussed in appendix A), there is no direct representation of ice–ocean thermodynamic interactions, no tidal flows, and no tidal mixing in the cavity. Also, the numerical model implements an advection and spatial discretization scheme that is not optimized for the study of shocks. However, with high horizontal and temporal resolution, and an observed steadiness in shock location, we are able to resolve these features with improved accuracy. Layer-grounding is treated in a special way in our simulations using Salmon layers, which prevent instabilities involving layer-grounding.

## 9. Conclusions and future directions

Our results corroborate well with previous modeling studies with idealized geometry (De Rydt et al. 2014; De Rydt and Gudmundsson 2016), and our idealized posing and high-resolution sweep over the key dimensionless parameters, allowing an exhaustive exploration of the dynamics. Based on this study, there are direct implications for future observations under ice shelves. These predictions may be used to guide measurements of the flow properties in ice shelf cavities and look for western boundary currents and shocks that may cause large diabatic fluxes, observable in western boundary currents downstream of the sill maximum. Such a region is likely to be a location of elevated mixing due to the sharp isopycnal tilt that forms when the bottom layer thins as it flows over the sill. The thinness of this layer and measured transport should have an influence on the magnitude of the shock amplitude and therefore, the turbulent mixing of subglacial water masses. Enhanced mixing near the sill is also likely to cause dynamical feedbacks; that is, when an initial pressure head inducing a circulation bringing heat toward the ice shelf causes increased water mass transformation, this further increases the pressure head. Measurements indicate that melting is concentrated close to the grounding lines (Dutrieux et al. 2014), but that it is nonnegligible elsewhere and the influence of that melt on the circulations we have explored remains to be investigated. Also, the importance of the missing surface forcing is important in nudging the stratification in otherwise more weakly stratified cases but has not been considered in this study.

Future open questions involve understanding how these dynamics interact with the evolving ice shelf and its impact on the cavity geometry and glacial melt rates. Additional considerations, including the variability of tidal and external currents, winds, accurate bathymetry, ice shelf geometry, and continuous stratification are all important next steps for the community to consider.

## Acknowledgments

The authors thank Pierre St-Laurent for allowing the use of the open source code BEOM and insightful comments on the manuscript, Aviv Solodoch for comprehensive suggestions on the manuscript, Janine Schaffer for useful discussions and hydrographic data, Tore Hatterman for extensive feedback, and two anonymous reviewers for helpful comments. The model source code is available at www.nordet.net/beom.html. This material is based in part upon work supported by the National Science Foundation under Grant PLR-1543388.

### APPENDIX A

#### Model Sensitivities

Although we have reduced the scope of our problem to heavily idealized cases, we can determine the sensitivity of the model results with three additional considerations: horizontal grid spacing, top topography, and quadratic friction. We find that in all cases the dynamics remain qualitatively unchanged, and our results may be expected to translate well to more realistic cavity configurations.

We find empirically that the transport, interface, velocity, and sill crossing location are not sensitive to horizontal grid spacings higher than m, indicating that this resolution is sufficient for the statistical flow properties presented in main text. Figure A1 shows a comparison of bottom-layer PV for four different grid spacings. Even at lower resolutions, the qualitative features of the flow are similar, but there are significant differences in boundary current velocity and transport. Both the overall and gyre-like recirculation strength time-averaged over the last 100 days between resolution m and m vary by less than 10% and shows a reasonable degree of convergence for m.

In the main text, we restricted our attention to flat-topped rigid lids for simplicity. Therefore, these results are also applicable to the free surface exchange flows since the surface wave effects are negligible for our chosen parameters. However, we can test varying cases of top topographies, as shown in Fig. A2, which does not significantly impact the observable dynamics of the flows, except in extreme cases, which are possible since ice shelf drafts can vary by over hundreds of meters in domains with bathymetry comparable to ours. Top topography is defined using with the top slope starting at the midpoint and is flat north of this, which is a crude approximation for the PIG geometry (Kimura et al. 2016). We found no significant differences as long as the sill’s topographic beta is much larger than the top topographic beta, which is true except for very gentle sills.

We also tested longer domains and bathymetry for varying sill width scale using meridionally stretched domains. These results did not show any qualitative difference when mapped back onto the regular domain except in the HF regime, where the boundary layer widens, which is consistent with our theory that .

Finally, we tested the effect of a realistic quadratic friction with a dimensionless coefficient of 2 × 10^{−3}, in contrast with the range of linear friction coefficients we have explored in the main text. With this quadratic friction, we observe the same dynamical regimes discussed in this study for equivalent values of friction depending on the specified pressure head. For example, with a pressure head of m and m, the quadratic solution exhibits a similar mean and RMSD transport to the linear friction case with , as shown in Fig. 2.

### APPENDIX B

#### Forward–Backward Time-Stepping in a Rigid Lid

We present a method for finding a convergent, accurate pressure field at the top surface for the forward–backward time-stepping scheme. A similar approach can also be taken for the generalized forward–backward time-stepping scheme.

Traditionally, in the forward–backward time-stepping scheme with a free-surface, each layer depth *h* is calculated forward in time based on the velocity fields *u* and *υ* at a time *t* [i.e., ], while the velocities are calculated backward in time based on as

where is the momentum tendency from Eq. (1). The meridional velocity *υ* benefits from a backward Coriolis term due to and, to minimize errors, this process alternates the term that is updated first.

However, with a rigid lid implementation, the velocities are only partially updated via an explicit time step, and are then corrected to include the tendency due to the surface pressure gradient. Specifically, the velocity vector for each layer is updated so that , using the definition for the top-layer Montgomery potential , so all the forces in the momentum balance are included except surface pressure gradient. This results in a force decomposition where

The initial velocity after a partial update, including all forces except pressure gradient, becomes , and after including the surface pressure gradient force becomes . This is calculated for each time step in this order:

In Eq. (B3c), the half-updated velocity can no longer include the backward Coriolis contribution from , since continuity must be satisfied at full steps *t* and and not half steps; that is, and are satisfied, but is not necessarily satisfied. Therefore, cannot be updated using since continuity is satisfied for , calculated based on , while it is not necessarily satisfied for .

The velocities with the fully updated pressures must have zero depth-integrated divergence, satisfying the continuity equation at every grid point in the horizontal plane. Thus, the pressure gradient is constrained in terms of the known partially updated variables.

### APPENDIX C

#### Thickness-Weighted PV Balance

In our problem, the PV balance allows us to identify the dominant mechanisms that facilitate cross-sill exchanges in various parameter regimes.

By taking the curl of the momentum equation, the mean shallow water thickness-weighted PV equation for a homogeneous fluid layer with friction and viscosity is

This equation is time-averaged and advection can be partitioned into mean and eddy components of PV flux with the thickness-weighted identity

where we define the thickness-weighted average and deviations from that average in Eq. (26). This allows Eq. (C1) to be written as

The diabatic terms can be grouped together under the continuity equation in Eq. (2) as

This equation states that in a frame following the thickness-weighted average velocity, the thickness-weighted PV is modified by eddy PV flux convergence, torques imposed by friction and viscosity, and stretching imposed by diabatic velocities.

### APPENDIX D

#### Rotating Two-Layer Uniform PV Solution

In addition to the rotating one-layer zero PV solution provided in section 7, we present a rotating two-layer uniform PV solution following Pratt and Armi (1990). This discussion is restricted to isopycnal layers of inviscid fluid flowing due to pressure and gravity in a channel of slowly varying cross section. The two-layer channel flow solution over bathymetry assumes uniform PV for each layer *n*

where is the upstream layer thickness. The Bernoulli functions for the two-layer case are defined as

and are constant along streamlines. The prescribed integration constants and are defined as fixed potentials infinitely upstream and are chosen based on the interface depth for a quiescent flow. The analog of this for our finite domain is defined using the thicknesses at the northern boundary conditions , .

Using Eqs. (D1)–(D2b) and Eqs. (30a)–(30c), the general forms of solutions for thickness and along-channel velocity are

where are constants that can be determined in terms of and analytically (Pratt and Armi 1990). Our analytical predictions for the critical transport based on this theory were calculated iteratively by finding the Bernoulli constants that exceeded the threshold for critical velocities at the sill maximum.

As prescribed pressure head increases, the analytical solution reaches a critical and eventually a supercritical state, which generally leads to PV adjustment processes that cause the breakdown of semigeostrophic theory and uniform PV becoming a poor assumption (Pratt and Armi 1990; Dalziel 1991). Therefore, there must be nonuniform PV effects, as otherwise the cross-channel velocity gradient would be highly supercritical for channel widths greater than , due to the barotropic cross-channel shear necessary to maintain uniform PV. Specifically, in our numerical solutions, we find that uniform PV is a reasonable approximation in the wide channel limit only within a western boundary layer of scale width , and the flow transitions to a smaller PV in the quiescent interior. This is also supported by analytical nonuniform PV models and laboratory experiments. For instance, Borenas and Whitehead (1998) demonstrated in laboratory experiments that upstream recirculation with flow reversals are necessary for flow separation and maintain criticality due to the recirculatory flow. Either recirculation or nonuniform PV is therefore unsurprising in the wide limit to maintain hydraulically controlled states since highly supercritical flows are unlikely to remain steady.

### APPENDIX E

#### Rotating Uniform PV Layer-Grounding Theory

In the study, layer-grounding occurs for sufficiently high sills with supercritical flow. Here, we derive a theoretical estimate for the boundary current transport necessary to achieve this state using the one-layer Eqs. (30a)–(30c) for a boundary current of uniform PV of width , to be determined, that evolves in the along-stream direction, adjacent to a stagnant interior. We consider an inflow coming from the northern boundary, between the western wall and . Implicitly, at from Eq. (30c). The boundary conditions for the western boundary current are

where is the along-channel transport.

The boundary conditions in Eqs. (E1a) and (E1b) combined with Eqs. (30a)–(30c) allows us to write the solutions for *h* and *υ* as a function of and *x*:

Applying the third boundary condition in Eq. (E1c), the expression involving transport can be written as

which yields two solutions for the boundary current width

which is and changes as a function of bathymetric height, upstream forcing, and transport, reaching a minimum at the sill maximum. Layer-grounding can occur when using Eq. (E2a), but the solution branch can be excluded since the thickness is always negative somewhere.

Therefore, the condition that predicts layer-grounding is when the discriminant in (E4) becomes negative:

This represents a sharp transition of a boundary current of width to a boundary of undefined width (i.e., grounding). This appropriately represents the behavior observed in our simulations since boundary currents have width immediately before layer-grounding, instead of a width that approaches zero.

## REFERENCES

*Indo-Pacific Climate Variability and Predictability*, S. K. Behera and T. Yamagata, Eds., World Scientific, 109–134.

*Rotating Hydraulics: Nonlinear Topographic Effects in the Ocean and Atmosphere*. Springer, 592 pp.

*Atmospheric and Oceanic Fluid Dynamics*. Cambridge University Press, 745 pp.

## Footnotes

Supplemental information related to this paper is available at the Journals Online website: https://doi.org/10.1175/JPO-D-18-0076.s1.

© 2018 American Meteorological Society. For information regarding reuse of this content and general copyright information, consult the AMS Copyright Policy (www.ametsoc.org/PUBSReuseLicenses).