## Abstract

Fjord-scale circulation forced by rising turbulent plumes of subglacial meltwater has been identified as one possible mechanism of oceanic heat transfer to marine-terminating outlet glaciers. This study uses buoyant plume theory and a nonhydrostatic, three-dimensional ocean–ice model of a typical outlet glacier fjord in west Greenland to investigate the sensitivity of meltwater plume dynamics and fjord-scale circulation to subglacial discharge rates, ambient stratification, turbulent diffusivity, and subglacial conduit geometry. The terminal level of a rising plume depends on the cumulative turbulent entrainment and ambient stratification. Plumes with large vertical velocities penetrate to the free surface near the ice face; however, midcolumn stratification maxima create a barrier that can trap plumes at depth as they flow downstream. Subglacial discharge is varied from 1–750 m^{3} s^{−1}; large discharges result in plumes with positive temperature and salinity anomalies in the upper water column. For these flows, turbulent entrainment along the ice face acts as a mechanism to vertically transport heat and salt. These results suggest that plumes intruding into stratified outlet glacier fjords do not always retain the cold, fresh signature of meltwater but may appear as warm, salty anomalies. Fjord-scale circulation is sensitive to subglacial conduit geometry; multiple point source and line plumes result in stronger return flows of warm water toward the glacier. Classic plume theory provides a useful estimate of the plume’s outflow depth; however, more complex models are needed to resolve the fjord-scale circulation and melt rates at the ice face.

## 1. Introduction

Convective motions are ubiquitous in the ocean and atmosphere, arising from statically unstable density differences between a source fluid and its environment under the influence of gravity. Buoyant plumes are one example of turbulent gravitational convection, where an isolated source of buoyancy drives the mean flow (Turner 1973). These plumes, typically characterized by turbulent dynamics with a high Reynolds number, occur over an enormous range of scales and are forced by both natural and anthropogenic sources, such as volcanic eruptions, hydrothermal vents at ocean ridges, fires, and the discharge of pollutants (Woods 2010).

The primary goal of this paper is to apply the general model of turbulent plumes to a timely environmental problem: the subglacial discharge of meltwater into Greenland’s outlet glacier fjords. The rate of mass loss from the Greenland Ice Sheet (GrIS) quadrupled over the last two decades and currently accounts for one-quarter of the global sea level (GSL) rise (Cazenave and Llovel 2010; Straneo and Heimbach 2013; Enderlin et al. 2014). Submarine melting due to warm ocean waters has been increasingly implicated as a major factor in controlling the stability and acceleration of outlet glaciers worldwide (Joughin et al. 2012; Straneo and Heimbach 2013). With a projected GSL rise of 0.5–1.2 m by 2100 under the IPCC “business as usual” scenario, understanding how a warming ocean contributes to dynamic mass loss from the GrIS is critical to predict and mitigate future climate change (Kopp et al. 2014).

Numerous studies indicate that subglacial discharge at the grounding line of Greenland’s glaciers drives a turbulent plume that rises along or near the ice face (Motyka et al. 2003; Salcedo-Castro et al. 2011; Xu et al. 2012, 2013; Sciascia et al. 2013). The turbulent plume entrains warm Atlantic-origin water (AW) at depth, providing a mechanism for delivering heat to the ice face. In a steady-state balance, the outflowing plume is balanced by a return flow of AW at depth. The resultant, turbulent plume–driven flow has been proposed as a mechanism for setting up a fjord circulation pattern that transports heat to Greenland’s outlet glaciers (Rignot et al. 2010; Straneo et al. 2012; Motyka et al. 2013); however, in situ evidence that this mode of circulation governs heat transport over any time scale is lacking because of the difficulty of making sustained measurements near the glacier terminus in these remote, ice-choked fjords.

Although we lack sufficient observations to test the plume-driven circulation hypothesis, progress has been made in understanding plume dynamics near the ice–ocean interface using numerical methods. Current models include 1D theory-based models (Jenkins 2011), 2D numerical simulations (Salcedo-Castro et al. 2011), more complex general circulation models (GCMs) (Xu et al. 2012, 2013; Sciascia et al. 2013, 2014), and finite-element methods (FEMs) (Kimura et al. 2014). Recent studies have parameterized plumes in coarser fjord-scale models (Cowton et al. 2015) and investigated the effect of subglacial hydrology on melt rates (Slater et al. 2015). Both 2D and 3D conduit geometries have been modeled with subglacial discharge entering the fjord either through a continuous crack along the grounding line (i.e., a line source) (Fig. 1a) or from a number of discrete subglacial conduits (Fig. 1b). To date, an eddy-resolving simulation with appropriate handling of the turbulent energy cascade has yet to be run, and instead models have relied on simplified parameterizations for the turbulent entrainment and diffusivity.

Previous investigators (Sciascia et al. 2013) have used line plume theory (Ellison and Turner 1959) as a basis for parameterizing turbulent entrainment in 2D models. However, this 2D line plume theory is not applicable to meltwater discharged from discrete subglacial conduits. Other 3D modeling efforts have attempted to resolve the turbulent plume (i.e., ≤1-m grid spacing) from an individual channel, but these models are too computationally intensive to investigate the far-field plume (Xu et al. 2013). Additionally, coarser-resolution models that resolve the fjord often assume idealized homogenous or two-layer ambient stratification (Sciascia et al. 2013; Kimura et al. 2014). Although these previous efforts have been useful, the necessary limitations imposed by numerical models combined with the expansive parameter space encompassed by the 200+ GrIS outlet glaciers highlight the need for further GCM simulations.

In this paper, we develop a general framework for determining circulation patterns forced by meltwater plume dynamics, using buoyant plume theory and a 3D, nonhydrostatic version of the Massachusetts Institute of Technology general circulation model (MITgcm) (Marshall et al. 1997). This study differs from previous work in that we shift the focus to downstream of the ice face, investigating how basic model parameters relevant to all outlet glacier systems (e.g., subglacial discharge rates, conduit geometry, turbulent diffusivity, and fjord stratification) determine fjord-scale circulation. This approach allows us to investigate the effect of the spatial distribution of submarine melting along the glacier terminus, identify the transition between point source and line plume regimes, and quantify the downstream properties and mixing of the plume. To correctly estimate the magnitude of ocean heat transport toward the glacier in future fjord-scale models, we need to resolve the vertical tracer structure and terminal level of the outflowing plume.

## 2. Background

### a. Physical setting

For the case studies considered here, the ambient temperature and salinity is based on data collected from Rink Fjord, west Greenland (Fig. 2). This deep, tidewater glacier is exposed to a complex density stratification that is typical of outlet glacier fjords in Greenland (Straneo et al. 2011; Chauché et al. 2014), providing an ideal physical setting to investigate meltwater plumes rising through strong vertical density gradients. Warm, salty AW occupies the bottom layer, overlaid by relatively cold, fresh polar water (PW). A thin layer of seasonal surface water (SW), consisting of runoff and ice melt warmed by solar radiation, is present near the surface. At many locations within the fjord, conductivity–temperature–depth (CTD) measurements were taken with a 6-Hz XR-620 RBR sensor. From these temperature and salinity profiles, we calculate a mean summer buoyancy frequency *(z)* profile, where the overbar represents an average of profiles taken during our synoptic survey (Fig. 2c). In the summer, stratification reaches a maximum at 30-m depth, at the strong density interface between the PW and SW layer.

Estimates of subglacial discharge in GrIS fjords, such as the Rink Fjord, are not well constrained and depend on unknown parameters such as ablation rate, catchment surface area, and the efficiency of the subglacial channel network (Chu 2014). Previous investigators have estimated summer subglacial discharges of ~200–500 m^{3} s^{−1} at Store Glacier (Xu et al. 2012, 2013), while estimates of surface meltwater entering Sermilik and Jakobshavn Fjord are ~174 m^{3} s^{−1} (Andersen et al. 2010) and 750–1500 m^{3} s^{−1} (Echelmeyer et al. 1991), respectively. Chauché et al. (2014) estimate a summer subglacial discharge into Rink Fjord of 1000 ± 300 m^{3} s^{−1} in 2009 and 1500 ± 450 m^{3} s^{−1} in 2010. Given the wide range and large uncertainties for estimates of subglacial discharge in GrIS fjords, we use a subglacial discharge flux of 1, 10, 75, 150, 300, 500, and 750 m^{3} s^{−1} to simulate a variety of possible summer discharge scenarios. The discharges tested were adequate to produce a range of subsurface and surface-trapped plumes (see section 4a for details).

### b. Buoyant plume theory

To develop insight into the interaction between subglacial buoyancy flux and stratification, we consider an idealized model of a nonrotating, axisymmetric turbulent plume discharged from a point source of buoyancy into a stratified ambient fluid. We build upon the classical model of point source plumes in a uniformly stratified environment (Morton et al. 1956, hereinafter referred to as MTT56). The MTT56 model has been widely used in modeling turbulent plumes for the last half-century and has been verified extensively using both laboratory experiments and observations over a diverse parameter space (Turner 1973; List 1982; Woods 2010). To apply this model to Greenland’s outlet glacier fjords, we assume that the point source plume is bisected into a half-plume by the vertical ice face, the buoyancy flux of the plume is due to steady subglacial discharge (i.e., we assume submarine melting has a second-order effect), the plume initially has no horizontal momentum, and the half-plume rises in a continuously stratified ambient fluid (Caulfield and Woods 1998). To represent a half-plume, we scale the prescribed subglacial discharge flux by a factor of 2 [(5)].

The MTT56 model is based upon three fundamental assumptions. First, the rate of horizontal entrainment at the edge of the plume *u*_{e} is linearly proportional to the vertical velocity *w* at that height:

with the dimensionless entrainment constant *α* (Fig. 1). The edges of a turbulent plume entrain quiescent water from outside, progressively increasing the mass flux of the plume as it rises. At a high Reynolds number, this process is independent of molecular viscosity. We assume that the plume velocity and buoyancy have a “top-hat” profile, where the value is constant inside the plume and zero outside. The value for the entrainment constant depends on the choice of the plume profile; we adopt a value of *α* = 0.13, representative of a pure plume forced by buoyancy alone (Linden 2000).

The second assumption is that mean vertical velocity and buoyancy exhibit self-similarity at all heights, with a conic shape emanating from the point source (Fig. 1b). Mass and momentum fluxes are defined as integrals of the mean values taken across the width of the plume. The third and final assumption is that local variations of density in the plume *ρ*′ are small compared to the background density *ρ*_{0} (where *ρ* = *ρ*_{0} + *ρ*′). Under this assumption, we invoke the Boussinesq approximation (i.e., fluxes of mass can be considered fluxes of volume).

The MTT56 formulation yields a model for three key characteristics of the plume as a function of height *z* above the source: the radius *b*, vertical velocity *w*, and reduced gravity *g*′ = *g* (*ρ*_{p} − *ρ*_{a})/*ρ*_{0}, where *g* is the gravitational acceleration, *ρ*_{p} is the density of the plume, and *ρ*_{a} is the ambient density. Fluxes of volume *Q*, momentum *M*, and buoyancy *F* are expressed as

where *r* is the radial distance. Averaging over the cross-sectional area of the plume leads to a system of three prognostic governing equations:

Following Caulfield and Woods (1998), we extend the original MTT56 formulation to have the buoyancy frequency *N* as a function of height above the point source instead of a constant:

This system of ordinary differential equations is integrated numerically using a fourth-order Runge–Kutta method (Burden and Faires 2011). Boundary conditions for volume *Q* and momentum *M* fluxes are zero at the grounding line depth. The buoyancy flux for the half-plume is prescribed at the boundary as

where *Q*_{sg} is the subglacial discharge flux. By scaling the solutions, we recover the diagnostic variables *w*, *b*, and *g*′:

Depending on the strength of the prescribed buoyancy flux and ambient stratification, the plume reaches neutral buoyancy at depth or at the free surface. We define the plume outflow depth as the level of neutral buoyancy, where *ρ*_{p} = *ρ*_{a} (i.e., where the buoyancy flux *F* = 0). As the plume reaches its maximum height (defined to be where the momentum flux *M* = 0), the assumptions of self-similarity and linear entrainment no longer hold and the integration is stopped. More complex models are required to resolve the terminal level and downstream properties of the plume.

## 3. Buoyant plumes in ocean GCMs

The extended MTT56 model provides a useful foundation for investigating point source plumes in a stratified environment; however, it only resolves plume properties along the upward centerline trajectory, assumes a constant entrainment coefficient and linear model for the entrainment rate, and does not incorporate a buoyancy flux from the melting ice face. To overcome some of these limitations, we use the MITgcm in a 3D, high-resolution, nonhydrostatic configuration to simulate line and point source meltwater plumes rising along a melting glacier terminus. The MITgcm is a developed version of Marshall et al. (1997), which solves the primitive Boussinesq form of the Navier–Stokes equations on an Arakawa staggered C grid (Arakawa and Lamb 1977). The MITgcm is useful for modeling buoyant plumes because of its nonhydrostatic capabilities (Marshall et al. 1998) and has been used to simulate nonhydrostatic dynamics in glacier environments (Xu et al. 2012, 2013; Sciascia et al. 2013, 2014; Gladish et al. 2015; Cowton et al. 2015; Slater et al. 2015).

### a. MITgcm configuration

We use the MITgcm to investigate the sensitivity of near-glacier and downstream plumes to variations in horizontal eddy diffusivity *κ _{H}*, subglacial discharge

*Q*

_{sg}, and the number of subglacial conduits (Table 1). We define downstream as oriented in the direction of the outflowing plume (i.e., away from the ice face in the along-fjord direction). The model domain is an idealized representation based on Rink Fjord, with a length of 100 km and depth of 850 m. The along-fjord horizontal resolution Δ

*x*is 10 m within 1 km of the terminus, telescoping to 10 km at the open western boundary (fjord mouth). A fjord width of 1 km is used, with a uniform across-fjord resolution Δ

*y*of 10 m and periodic lateral boundary conditions. For the range of

*Q*

_{sg}tested, the fjord width was adequate to prevent the plume from reaching the lateral boundaries in the near-glacier field. A uniform vertical resolution Δ

*z*of 10 m with 85 levels was used. Temperature and salinity are restored to the prescribed initial conditions at the open western boundary, which includes a 50-km restoration region to prevent internal waves reflecting back into the near-glacier field.

Near the glacier, the internal Rossby radius of deformation is larger than the fjord width, and we assume that rotational effects have a second-order effect (i.e., *f* is set to zero in the model simulations). The equation of state (JMD95Z) follows Jackett and Mcdougall (1995). The model has a nonlinear free surface and no-slip conditions enforced at the seabed and ice face. Additional simulations with a free-slip condition on the ice face resulted in slightly fresher plumes with larger vertical velocities near the grounding line. Bottom drag is parameterized with a quadratic drag law coefficient *C*_{D} = 2.5 × 10^{−3}. We use a volume-conserving “virtual” salt flux for the subglacial discharge, where the freshwater flux is implemented as a salt flux (Huang 1993; Sciascia et al. 2013) and has no associated mass flux in the continuity equation. The Rink Isbræ glacier terminus is represented as a 850-m-deep solid boundary, with a uniform ice temperature of −2°C. Submarine melting and freezing processes on the ice wall are parameterized with a system of equations that represent the conservation of heat and salt combined with a linear equation for the freezing point of seawater (Hellmer and Olbers 1989; Holland and Jenkins 1999; Losch 2008; Xu et al. 2012).

We perform a range of experiments with subglacial discharge values of 10, 75, 150, 300, 500, and 750 m^{3} s^{−1} injected at the grounding line. The subglacial discharge has a salinity of 0 psu with temperature set to the salinity pressure–dependent freezing point at the grounding line depth (−0.628°C). The efflux velocity is linearly restored over a 1-day time scale to increase numerical stability.

The relative importance of buoyancy and inertial forces in a meltwater plume can be described by the Froude number (Syvitski 1989; Mugford and Dowdeswell 2011). The Froude number can also be written as a ratio of the Reynolds and Grashof numbers (Arakeri et al. 2000; Salcedo-Castro et al. 2011):

The Reynolds number is defined as Re = *uh*/*ν*, where *u* is the horizontal velocity in the subglacial conduit, *h* is the height of the subglacial conduit (30 m), and *ν* is the kinematic viscosity of seawater. The Grashof number characterizes the buoyancy flux:

For our upper-bound *Q*_{sg} of 750 m^{3} s^{−1}, Re ~ 2.5 × 10^{7} at the subglacial conduit and Gr ~ 7.3 × 10^{15}, ensuring a Fr < 1. Therefore, our simulations correspond to a plume regime dominated by buoyancy forces (i.e., the injected subglacial discharge does not result in a jet), allowing for comparison with the MTT56 model.

### b. Parameterization of turbulent entrainment

Turbulent entrainment processes in the MITgcm are not resolved and are parameterized here by a constant eddy diffusivity and viscosity. Sciascia et al. (2013) demonstrated that for 10-m grid resolution, the choice of horizontal eddy diffusivity *κ*_{H} was critical to correctly represent the dynamics of a 2D line plume modeled in MITgcm. The choice of *κ*_{H} determines the rate of turbulent entrainment, constraining the reduced gravity and terminal level of the plume. Xu et al. (2012) used a Laplacian vertical viscosity of 0.1 m^{2} s^{−1} and a biharmonic *κ _{H}* of 300 m

^{4}s

^{−1}for a horizontal grid resolution of 20 m. Finer-scale models have used a Laplacian

*κ*

_{H}of 0.25 (Sciascia et al. 2013, 2014) and 0.01 m

^{2}s

^{−1}(Xu et al. 2013) for 10- and 1-m grid resolution, respectively.

For our experiments, we vary a biharmonic *κ*_{H} to evaluate the sensitivity of the downstream plume to this parameter. The horizontal eddy diffusivity and viscosity are set to be equal and range from 0.125, 0.25, and 0.5 m^{4} s^{−1}. The eddy diffusivities are chosen to be large enough to suppress numerical instabilities, yet small enough to allow the model to reproduce sharp gradients and eddies. Following Sciascia et al. (2013), we use a constant Laplacian vertical eddy viscosity of 10^{−3} m^{2} s^{−1}. A third-order, direct-space, time flux–limited advection scheme is used to increase the accuracy of the plume front and eliminate extrema in the tracer field. The simulations are integrated for 5 days, which is sufficient time for the plume to span the 1-km, 10-m resolution subdomain.

### c. Conduit geometry and plume interaction

To simulate 3D point source plumes, subglacial discharge is injected horizontally into the fjord through a single 30 m × 30 m conduit at the grounding line depth. Additionally, we simulate a line plume and multiple point sources to evaluate how subglacial conduit spacing affects the downstream properties of the plume in a 3D domain. For the line plume, subglacial discharge is injected uniformly as a buoyancy source across the entire width of the glacier through a 1000 m × 30 m conduit at the grounding line depth with an efflux velocity of 3.33 × 10^{−4} m s^{−1}. To investigate the transition between coalescing point source and line plume regimes, we perform experiments with multiple 30 m × 30 m conduits at 100-, 200-, and 400-m spacing distributed across the 1-km glacier width (10, 5, and 3 conduits). The efflux velocities for the 100-, 200-, and 400-m conduit spacing are 1.11 × 10^{−3}, 2.2 × 10^{−3}, and 3.7 × 10^{−3} m s^{−1}, respectively.

## 4. Results

### a. Extended MTT56 experiments

Numerically integrating the MTT56 governing equations with the continuous stratification profile from Rink Fjord gives solutions for the vertical structure of volume, momentum, and buoyancy flux (Fig. 3). The theory predicts a range of surface and subsurface plumes with ambient stratification representative of summer. For *Q*_{sg} = 75 and 150 m^{3} s^{−1}, the plumes reach maximum height in the PW and SW layers (<50 m). The lower-bound *Q*_{sg} of 1 m^{3} s^{−1} results in a subsurface plume that reaches a maximum height of 400 m. For *Q*_{sg} = 300 m^{3} s^{−1}, the plume penetrates the meltwater pycnocline and reaches the free surface (Fig. 3b); however, the level of neutral buoyancy remains near 100 m within the PW region.

For maximum *Q*_{sg} = 750 m^{3} s^{−1}, the level of neutral buoyancy reaches the free surface and the plume becomes surface-trapped (not shown). During the initial plume rise, the buoyancy flux is nearly constant because of the vertical homogeneity of the AW properties near the grounding line depth where the fluid is very weakly stratified (Fig. 3c). Once the plumes rise ~250 m above the grounding line depth, stratification increases as the ambient density begins to decrease (Fig. 2), as does the buoyancy flux. At the level of neutral buoyancy, the plume’s buoyancy flux crosses zero and the plume becomes negatively buoyant. The plume gradually loses its vertical momentum as it continues to rise through the AW and PW layers (Fig. 3b). For *Q*_{sg} < 75 m^{3} s^{−1}, the level of neutral buoyancy is reached well below the ~50-m-deep SW layer (Fig. 3c). For *Q*_{sg} ≥ 75 m^{3} s^{−1}, the plume’s residual vertical momentum allows for penetration into near-surface depths. Downstream of the glacier terminus, the terminal level of the plume should be bounded by the maximum height and the level of neutral buoyancy.

### b. MITgcm experiments

For the 3D MITgcm point source simulations, the discharge of buoyant meltwater at the grounding line results in a turbulent plume that rises vertically along the glacier terminus. The plumes remain attached to the vertical glacier terminus until they reach their maximum height, at which point they separate and flow horizontally away from the ice face (Fig. 4). Note that we use the term “plume” to describe the both vertical plumes rising along the ice face and horizontal buoyancy-driven currents that flow away from the glacier as surface gravity currents or subsurface intrusions. As the plumes rise, they entrain dense ambient fluid, increasing the volume flux and expanding radially (Figs. 4, 5). The initial dilution of the plume is determined by the rate of subglacial discharge and the choice of horizontal eddy diffusivity/viscosity. The salinity of the plume in the first wet cell adjacent to the subglacial conduit is 32.98 and 25.46 for a *Q*_{sg} of 10 and 150 m^{3} s^{−1}, corresponding to a reduced gravity *g*′ of 0.01 and 0.07 m s^{−2}, respectively. The reduced gravity for pure meltwater at the grounding line is 0.22 m s^{−2}, demonstrating that the plume is quickly diluted with entrained AW as it exits the subglacial conduit.

The model simulations are in good agreement with the MTT56 level of neutral buoyancy and maximum height, implying that the MITgcm bulk plume properties are consistent with theory. The simulations also reproduce the inertial overshoot that is evident in the MTT56 solutions. The inertial overshoot occurs uniformly in the radial direction, constrained only by the vertical glacier terminus. For large subglacial discharges, the horizontal length scale for a turbulent plume to transition into a steady outflow can exceed 500 m (Fig. 6a). We note that for our largest *Q*_{sg} of 750 m^{3} s^{−1}, the plume did not exhibit an inertial overshoot and remained surface trapped (not shown). At 1 km downstream, the terminal level of the outflowing plumes are centered near the MTT56 level of neutral buoyancy (Figs. 6a,b). A *Q*_{sg} of 10 m^{3} s^{−1} results in a weak outflow, with a maximum horizontal plume velocity of 0.04 m s^{−1} 1 km downstream of the terminus (Fig. 6b). A *Q*_{sg} of 300 m^{3} s^{−1} results in a faster-flowing plume, with a maximum downstream horizontal velocity of 0.20 m s^{−1} (Fig. 6b). For *Q*_{sg} < 150 m^{3} s^{−1}, submarine melting along the glacier terminus results in a weak inflow/outflow below the terminal level of the discharge-driven plume. This small-scale structure is amplified and biased toward inflow as subglacial discharge is increased. The relatively slow downstream plume velocities provide a sharp contrast to the large vertical velocities at the glacier terminus, which are on the order of meters per second. Our results also highlight the importance of the seasonal SW layer in bounding the terminal level of the plume. For large subglacial discharges, vertical momentum enables the plume to penetrate through the stratification maximum at the PW–SW interface. However, for *Q*_{sg} < 750 m^{3} s^{−1}, the level of neutral buoyancy is bounded by the SW layer, resulting in a plume that reaches the free surface at the ice face and is trapped underneath buoyant SW as it flows away from the glacier.

A comparison of plume water properties at the glacier terminus and 1 km downstream is shown in Fig. 7. At the glacier terminus, the potential temperature–salinity (*θ*–*S*) properties of the plume fall along the runoff mixing line, indicating the mixing of fresh subglacial discharge with the deep AW layer. For both low and high subglacial discharge fluxes, the near-glacier and downstream water properties are bounded by the runoff and melt lines. The heat required to melt ice combined with the mixing of meltwater with ambient fluid results in a melt line in *θ*–*S* space, with a so-called Gade slope:

where*T* is temperature; *S* is salinity; *L* is the latent heat of fusion for ice; *c*_{i} and *c*_{p} are specific heat capacities of ice and water; *T*_{i} is the temperature of ice; and *T*_{f} is the freezing point temperature (Gade 1979; Jenkins 1999; Mortensen et al. 2013). For the low discharge case of 10 m^{3} s^{−1}, the SW layer is more buoyant than the initial dilution of the plume as it exits the subglacial conduit (Fig. 7a). Compared to the initial *θ*–*S* conditions, the profile 1 km from the glacier shows warming and a slight increase of salinity at depth, indicating the influence of glacially modified AW from the plume (Fig. 7b). Compared to the low discharge case, a high discharge of 150 m^{3} s^{−1} exhibits drastically different downstream *θ*–*S* properties and ambient water mass modification. At this *Q*_{sg}, the plume reaches a maximum height at ~5 m below the free surface, encompassing the entire water column at the glacier terminus. The dynamics of the fjord are determined by the large subglacial discharge flux, with submarine melting having a second-order effect. As it rises, the plume vigorously entrains heat and salt into the PW/SW layers, increasing the temperature and salinity of the upper water column.

### c. MITgcm sensitivity analysis

Both the vertical velocity and maximum height of the plume at the glacier terminus are robust to the tested range of *κ*_{H} (Fig. 8). An increase of *κ*_{H} by a factor of 2 results in slightly increased vertical velocities (<10%) and turbulent entrainment. For a *κ*_{H} of 0.25 m^{4} s^{−1}, a subglacial discharge flux of 10 m^{3} s^{−1} results in a maximum vertical velocity of 0.78 m s^{−1} (Fig. 8a). Subglacial discharge fluxes of 75 and 150 m^{3} s^{−1} produce a more vigorous plume, with a maximum vertical velocity of 1.46 and 1.79 m s^{−1}, respectively (Fig. 8b,c). Away from the grounding line depth, the MITgcm vertical velocities agree well with the MTT56 solutions (Figs. 8a,b,c). At the grounding line, the MTT56 vertical velocities asymptote to infinity due to the zero volume flux (*Q* = 0) boundary condition. The corresponding MITgcm momentum fluxes (Figs. 8d,e,f) are within a factor of 2 of the MTT56 solutions, with the largest discrepancies occurring at midcolumn depths.

Volume transports (*Q*_{out} and *Q*_{in}), binned by salinity class, indicate the shift in entrainment between model runs (Fig. 9a). Positive transports indicate flow toward the glacier terminus; negative values represent the outflowing plume. For all simulations, the net flow across the open boundaries is balanced (i.e., *Q*_{out} = *Q*_{in}). Volume transports are averaged between 900 m and 1 km downstream of the glacier and computed in salinity bins of 0.05. To provide an estimate of the bulk salinity of the in- and out-fjord transports, we compute volume weighted salinities as

where *u*_{in} and *u*_{out} and *S*_{in} and *S*_{out} are the in- and out-fjord velocities and salinities, respectively. A thin return flow of SW with a salinity of ~32.25 is present in all cases and is not included in the calculations. A subglacial discharge flux of 10 m^{3} s^{−1} results in a multicell circulation in the fjord (Fig. 9a), with a deep cell driven by subglacial discharge and a weaker upper-column cell driven by submarine melt. For higher *Q*_{sg}, the circulation is primarily two-layer, with a single inflow and outflow. An increase of *κ*_{H} results in a weaker, more viscous transport of AW toward the glacier (Fig. 9b). Increasing *κ*_{H} results in a slightly more diluted plume, shifting the downstream outflow to lower salinity classes. For *κ*_{H} of 0.25 m^{4} s^{−1}, this gives = 34.52, 34.49, 34.47, and 34.44 for a corresponding subglacial discharge flux of 10, 75, 150, and 300 m^{3} s^{−1}. The corresponding values are 34.34, 34.06, 33.90, and 33.81.

### d. Comparison of line and point source plumes

To investigate the role of subglacial conduit spacing on the downstream properties of the plume, we simulate multiple point source plumes spaced at 400-, 200-, and 100-m intervals across the glacier terminus (Fig. 10). A *Q*_{sg} of 10 m^{3} s^{−1} is used in these simulations. For a conduit spacing of 400 and 200 m, the plumes remain separated and interact weakly at their lateral boundaries (Figs. 10a,b). When the conduit spacing is reduced to 100 m, the plumes coalesce and begin to approximate a continuous line plume (Fig. 10c). The distribution of the subglacial discharge through multiple conduits is modeled by a lower efflux velocity, decreasing the vertical velocity of the plume.

Volume transport binned by salinity class for 1, 3, 5, and 10 point source plumes and a continuous line plume is shown in Fig. 11a. For both line and point source plumes, the flux of subglacial discharge drives a multicell circulation in the fjord, with multiple inflow and outflows distributed across a range of salinity classes. The maximum positive volume transports are centered between salinity classes of 34.49–34.52, corresponding to the return flow of AW toward the glacier terminus. While the return flows fall within a similar range of salinity classes, the outflows vary depending on the terminal level and plume type (i.e., line or point source). The line plume strongly dilutes the subglacial discharge, resulting in a much larger *Q*_{in} compared to the initial *Q*_{sg} and driving a vigorous buoyancy-driven circulation in the fjord. As the conduit spacing increases, the magnitude of volume transport decreases, reaching a minimum in the case with a single subglacial conduit. Line plumes with subglacial discharge fluxes of 10–150 m^{3} s^{−1} result in deep outflows centered at depths > 300 m (Fig. 11b).

## 5. Discussion

The 3D MITgcm plume simulations presented here are in good agreement with the extended MTT56 model, demonstrating that ocean GCMs can capture the bulk properties of a point source plume based on similarity theory. While the MTT56 model provides a useful tool for exploring the range of plume regimes that exist across outlet glacier parameter space (i.e., subsurface versus surface trapped), 3D models such as the MITgcm are essential to capture the fjord-scale circulation and spatial distribution of submarine melt across the ice face, which has important consequences for grounding line stability and glacial undercutting. The application of the MTT56 model to estimate submarine melt rates is limited, particularly in shallow fjords, because of the vertical velocity tending to infinity at the grounding line. A point source formulation of more complex models such as Jenkins (2011), where the initial vertical velocity is prescribed by assuming a balance between buoyancy and frictional drag, would provide a computationally efficient tool for estimating melt rates in systems with discrete subglacial conduits. We find that the use of realistic stratification profiles and a reasonable choice of horizontal eddy diffusivity are critical to represent the bulk dynamics of the plume, consistent with the results of previous modeling studies (Sciascia et al. 2013; Kimura et al. 2014). We note that the MTT56 model is only valid for a summer regime, where plume dynamics are dominated by subglacial discharge and submarine melting has a second-order effect. The downstream outflow depths of the modeled line plumes are biased high compared to the level of neutral buoyancy calculated using the model from Jenkins (2011). This bias is reduced as the discharge rate is increased (Fig. 11b).

Our study demonstrates that the terminal level and water properties of meltwater plumes can vary dramatically with increasing distance from the glacier terminus. The presence of strong stratification can result in plumes that reach the free surface in the near-glacier field and then transition into subsurface intrusions at roughly the level of neutral buoyancy downstream (Fig. 6). After the initial overshoot of the plume, several dampened vertical oscillations about the level of neutral buoyancy occur, which act to further mix ambient fluid into the plume and modify water properties as the plume flows away from the ice face. In this region, the plume behaves as a fountain, consisting of a turbulent shear flow driven by upward vertical momentum with an opposing negative buoyancy flux (Turner 1966; Ansong et al. 2008; Kaye 2008). We note that the use of constant eddy diffusivity does not allow for representation of internal wave or shear-driven mixing, which may be important in this region.

To investigate a wider range of subglacial discharge fluxes and stratification profiles, we use the extended MTT56 model to explore plume outflow depth in GrIS outlet glacier fjord parameter space (Fig. 12). MTT56 was run with both constant *N* (Fig. 12a) and idealized stratification profiles with varied pycnocline depths (Fig. 12c), representing a range of seasonal stratification profiles in a typical GrIS fjord. For constant *N* < 3 × 10^{−4} s^{−1} and *Q*_{sg} > 150 m^{3} s^{−1}, the level of neutral buoyancy reaches the free surface and the plumes become surface trapped (Fig. 12b). A pycnocline with *N*_{max} of 0.02 s^{−1} bounds the level of neutral buoyancy at depth for the entire range of *Q*_{sg} tested, resulting in subsurface outflows (Fig. 12d). Subglacial discharge of meltwater into strongly stratified outlet glacier fjords with deep grounding lines results in subsurface outflows; for shallow fjords with weaker stratification, we would expect plumes to be surface trapped.

During winter, stratification likely decreases near the surface due to the lack of SW. However, plumes may still be trapped at depth because of the seasonal reduction of *Q*_{sg}. For our 3D point source simulations, low subglacial discharge results in a multicell circulation in the vertical (Fig. 6), consistent with observational results from Sermilik Fjord, east Greenland (Straneo et al. 2011; Sutherland and Straneo 2012). Sciascia et al. (2013) showed that for low subglacial discharge fluxes in a two-layer stratification, large amounts of relatively fresh glacially modified water can be exported at depth at the AW–PW interface. Our results indicate that large discharges of meltwater in deep, strongly stratified fjords can result in plume outflows with positive temperature and salinity anomalies (Figs. 4, 5). In this flow regime, turbulent entrainment along the glacier terminus acts as a mechanism to vertically transport heat and salt upward in the water column. These results are in agreement with recent fjord-scale models (Cowton et al. 2015) and hydrographic observations from Godthåbsfjord, west Greenland (Kjeldsen et al. 2014). During ice-dammed lake drainage events, large amounts of warm, saline water were brought to the near-surface layers, indicating turbulent entrainment of AW at the glacier terminus. Warm plumes may result in increased submarine melting of icebergs and sea ice, resulting in a weakening of ice mélange and additional fluxes of freshwater to the ambient stratification.

This study also demonstrates the importance of having a priori knowledge of the subglacial hydrology and conduit/ice face geometry. With constraints on fjord stratification, grounding line depth, and conduit geometry, our results suggest that plume properties could be used to provide a back-of-the-envelope calculation of subglacial discharge rates, which are difficult to measure in situ. During summer conditions, variability in surface ablation results in a time-evolving network of subglacial conduits (Cowton et al. 2013), which could lead to a wide range of line and point source plumes rising along the glacier terminus. Previous laboratory and modeling results have shown that when two point source conduits are within close proximity to each other, separate axisymmetric plumes can coalesce into a single plume (Kaye and Linden 2004; Cenedese and Linden 2014). Our simulations with multiple conduits agree with recent studies by Kimura et al. (2014) and Slater et al. (2015), demonstrating that narrow conduit spacing in deep fjords can result in coalescing plumes. In our simulations, coalescing point source and line plumes result in a more vigorous buoyancy-driven fjord circulation (Fig. 11a). For the equivalent *Q*_{sg}, line plumes produce deeper, more diluted outflows than point source plumes. These results suggest that current 2D models may overestimate the strength of the buoyancy-driven circulation in systems where subglacial discharge is injected through discrete conduits. The discharge of meltwater through a distributed system of subglacial conduits would result in weaker efflux velocities at the grounding line, decreasing the vertical velocity of the plumes. However, the increased surface area available for turbulent heat transfer to the glacier face may result in a net increase of melt across the entire glacier terminus, especially for wide glaciers. Recent results from Slater et al. (2015) demonstrate that subglacial hydrology has an important control on submarine melt rates, with the transition from a single conduit to a distributed system increasing melt by a factor of 5. Clearly, further studies are needed to assess the impact of ice face topography, channel shape, and plume regimes on submarine melt rates.

We note that these results rely on the simplified assumption of fixed boundary conditions at the fjord mouth, which provide a constant source of ocean heat in the AW layer. In Greenland, the renewal of basin water may play an important role in the heat budget of fjords with sills. Observations from Godthåbsfjord show that dense inflows contribute to the basin water, lasting 1–3 months, with warm inflows related to pulses in the AW content of the West Greenland Current (Mortensen et al. 2011, 2013). It is unclear if this renewal process is ubiquitous across other Greenland outlet glacier fjords; however, preliminary results from Ilulissat Icefjord suggest that renewal occurs over the shallow sill and toward Jakobshavn Glacier (Gladish et al. 2015).

Dynamics other than buoyancy-driven flows—tides, katabatic winds, and intermediary circulation (Straneo et al. 2010; Jackson et al. 2014; Sciascia et al. 2014; Sutherland et al. 2014)—may also be critical to heat transport and the renewal of water in the fjords. Other limitations of our model include the highly idealized model domain. Many Greenlandic fjords, including Rink, have one or multiple sills, which may strongly constrain the plume-driven circulation and the return flow of AW toward the glacier. The neglected rotational effects and cross-fjord circulation may be important farther downstream of the glacier where the fjord width widens to ~15 km. Additionally, the use of a constant eddy diffusivity does not represent the full range of mixing processes in the overshoot region, which may result in deeper outflows than an eddy-resolving simulation.

## 6. Conclusions

Greenland’s fjords couple the ice sheet margins to the continental shelf, forming the boundary between the accelerating GrIS and the coastal ocean. We use buoyant plume theory and nonhydrostatic, three-dimensional MITgcm simulations to investigate how meltwater plume dynamics and the resultant fjord-scale circulation depend on subglacial discharge, ambient stratification, turbulent diffusivity, and subglacial conduit geometry. While classic plume theory provides a useful estimate of the plume’s outflow depth, 3D ocean–ice models are needed to resolve the fjord-scale circulation and spatial distribution of submarine melting along the ice face. Plumes with large vertical velocities penetrate to the free surface near the ice face; however, subsurface stratification maxima can create a barrier that can trap plumes at depth as they flow away from the glacier. Large subglacial discharges result in vigorous turbulent entrainment of bottom water, transporting heat and salt upward in the water column and producing warm, salty outflows. We find that fjord-scale circulation is highly sensitive to subglacial conduit geometry; multiple point source and line plumes result in stronger return flows of Atlantic water toward the ice face. Our results indicate the need for further observations that constrain the heat and salt budgets between GrIS outlet glacier fjords and the coastal ocean. In situ measurements of turbulent dissipation rates and velocity/tracer profiles in the near-glacier plume are critical.

## Acknowledgments

This work was partially supported by the National Aeronautics and Space Administration Grant NNX12AP50G and the University of Oregon. We thank Captain H. Rink and the crew of the R/V *Sanna* for their support in making oceanographic measurements in Rink Fjord.

## REFERENCES

*Advances in Research and Applications,*J. Chang, Eds., Methods in Computational Physics Series, Vol. 17, Elsevier, 173–265, doi:.

*Perspectives in Fluid Dynamics,*G. K. Batchelor, H. K. Moffatt, and M. G. Worster, Eds., Cambridge University Press, 289–345.

**119,**7078–7098, doi:.

**42,**2861–2868, doi:.