Modeling the impact of glacial runoff on fjord circulation and submarine melt rate using a new subgrid-scale parameterization for glacial plumes Geophysical Research:

The injection at depth of ice sheet runoff into fjords may be an important control on the frontal melt rate of tidewater glaciers. Here we develop a new parameterization for ice marginal plumes within the Massachusetts Institute of Technology General Circulation Model (MITgcm), allowing three-dimensional simulation of large (500 km 2 ) glacial fjords on annual (or longer) time scales. We ﬁnd that for an idealized fjord (without shelf-driven circulation), subglacial runoff produces a thin, strong, and warm down-fjord current in the upper part of the water column, balanced by a thick and slow up-fjord current at greater depth. Although submarine melt rates increase with runoff due to higher melt rates where the plume is in contact with the ice front, we ﬁnd that annual submarine melt rate across the ice front is relatively insensitive to variability in annual runoff. Better knowledge of the spatial distribution of runoff, controls on melt rate in those areas not directly in contact with plumes, and feedback mechanisms linking submarine melting and iceberg calving are necessary to more fully understand the sensitivity of glacier mass balance to runoff-driven fjord circulation.


Introduction
The melting of Greenland's marine terminating glaciers by warm oceanic waters may be a significant cause of dynamic mass loss from the Greenland ice sheet [Holland et al., 2008;Straneo and Heimbach, 2013]. Relatively little is however known about this process, particularly the interaction between the glacier (and subglacial runoff), fjord circulation, and terminus melt rates . For example, although very high (>5 m d 21 ) melt rates have been predicted in areas directly adjacent to upwelling subglacial runoff [Jenkins, 2011;Xu et al., 2012], the importance of this process when considered across the entirety of the glacier front remains unclear. The significance of glacially forced estuarine circulation within the fjords as a means of replenishing the supply of warm water available for melting at the fjord head is also debated [e.g., Motyka et al., 2013;Straneo et al., 2010].
Obtaining field measurements from Greenland's ice-choked fjords has proven extremely difficult, particularly in the critical area close to glacier fronts . As such, data sets of water properties [e.g., Chauch e et al., 2014;Mortensen et al., 2011;Rignot et al., 2010;Straneo et al., 2011Straneo et al., , 2010 are rare and often limited with regard to timeframe and proximity to the glacier, making interpretation challenging. The convection of glacial plumes was modeled by Jenkins [2011], providing a theoretical relationship between meltwater discharge and submarine melt rate, but this approach cannot provide information on the impact of such plumes on the wider fjord circulation. In recent years, numerical ocean models have also been used to examine the impact of glacial runoff on fjord circulation [Salcedo-Castro et al., 2011;Sciascia et al., 2013;Xu et al., 2013Xu et al., , 2012. Such studies have however been limited by the computational expense of both resolving the fine-scale, highly nonhydrostatic dynamics of the turbulent glacial outflow plume and incorporating the wider circulation of fjords that may be hundreds of square kilometers in area. This constraint has resulted in simulations being limited to two-dimensional fjords [Sciascia et al., 2013;Xu et al., 2012] or small (<0.1 km 2 ) three-dimensional domains [Xu et al., 2013], and short (hours to days) simulation times.
Here we aim to overcome these constraints by parameterizing the glacial plume. Parameterization of vertical motion is a technique commonly employed in ocean modeling, allowing incorporation of processes occurring at a scale finer than the grid resolution [e.g., Campin and Goosse, 1999;Paluszkiewicz and Romea, 1997;Price and Yang, 1998]. We start by developing a plume model based on an input of subglacial runoff from a point source (i.e., a major subglacial channel), as is expected based on the observation of discrete plumes adjacent to tidewater glacier fronts [e.g., Chauch e et al., 2014;Sole et al., 2011]. This is then used to calculate the vertical movement of water, heat, and salt at each time step at a specified location within an ocean model. Doing so eliminates the need to resolve the glacial plume, thus permitting much coarser model resolution. This has the advantage of greatly reducing the simulation time while still allowing the plume to evolve as runoff varies or conditions in the fjord change.
This model is used to examine ice-ocean interaction on annual time scales in an idealized fjord, similar in scale to a large Greenlandic fjord such as Sermilik or Kangerdlugssuaq. Greenland's glacial fjords are complex environments, driven by a range of highly variable terrestrial, atmospheric, and oceanic forcings . We do not attempt to realistically recreate or assess the relative importance of the multiple components of fjord circulation. Instead we single out the glacial forcing, in particular assessing the impact of glacial melting and runoff on fjord circulation and the sensitivity of submarine melt rate to glacial runoff. In doing so, we seek to aid the interpretation of existing oceanographic data and assess the significance of glacially forced fjord circulation for the future stability of Greenland's tidewater glaciers. Furthermore, this idealized study provides an initial framework in which to develop and test the new plume parameterization, which we envisage applying to complex, realistic fjord systems in future work.

Model 2.1.1. Ocean Model
The Massachusetts Institute of Technology General Circulation Model (MITgcm, http://mitgcm.org) [Adcroft et al., 2004;Marshall et al., 1997aMarshall et al., , 1997b] is a versatile general circulation model that solves the incompressible Navier-Stokes equations using finite-volume methods on an orthogonal curvilinear grid. The model has a nonhydrostatic capability that makes it suitable for the study of areas of complex bathymetry and buoyant upwellings and includes options to incorporate horizontal ice shelves [Losch, 2008] and vertical ice faces [Xu et al., 2012]. To date, it has been used in several studies of glacial ice-ocean interaction both beneath Antarctic ice shelves [Heimbach and Losch, 2012;Losch, 2008] and in Greenland's fjords [Sciascia et al., 2014[Sciascia et al., , 2013Xu et al., 2013Xu et al., , 2012. We use the model in nonhydrostatic mode, but the plume parameterization (described in the following sections) is equally suitable for incorporating glacial plumes in relatively coarse hydrostatic simulations.

Plume Model
Theoretical models of turbulent plumes have been developed and utilized for over 50 years [e.g., Morton et al., 1956] and have proven highly effective at representing a wide range of oceanic and atmospheric phenomena [Kaye, 2008]. The plume model that forms the basis of the parameterization is motivated by that of Jenkins [2011], which coupled a buoyant plume model to equations describing heat exchange at the ice-ocean interface. Jenkins's [2011] model does not include plume expansion or entrainment in an across-fjord direction, and so is only appropriate in scenarios in which the input of runoff is uniform across the width of the glacier. At tidewater glaciers in Greenland, localized turbid plumes have been observed at the fjord surface in summer [e.g., Chauch e et al., 2014;Sole et al., 2011;Tedstone and Arnold, 2012], suggesting the presence of large subglacial channels providing a localized input of runoff to the fjord at the glacier grounding line. As such, an alternative plume geometry, permitting a discrete source of freshwater input, is required (M. O'Leary, Frontal processes on tidewater glaciers, University of Cambridge, unpublished thesis, 2011). In these cases, a half-conical plume geometry, with flat side against the glacier, is more appropriate (Figure 1). We therefore construct a model for an ice marginal plume of this geometry utilizing well-established plume theory [Kaye, 2008]. This model is described here, with symbols and parameter values defined in Table 1. A vertically issuing, conical plume model has previously been described by Morton et al. [1956]. This model was written in terms of density; we modify the equations therein to split out temperature and salinity and to take into account the half-conical geometry. Plume density q p is then given as a function of plume temperature T p and plume salinity S p using the modified UNESCO equation of state of Jackett and MacDougall [1995]. Following Morton et al. [1956], entrainment of ambient water into the plume is proportional to the plume velocity with rate a. We employ the standard Boussinesq approximation, which introduces a reference density q 0 into the equations. With plume radius b and vertical velocity u, the defining equations (1-4) state conservation of mass, momentum, heat, and salt: As in Jenkins [2011], the last term in (1) and the last two terms in (3) and (4) describe the melting of the ice, and the last term in (2) is a reduction in momentum due to friction, with additional multipliers to take into account the plume geometry. T a , S a , and q a are the ambient fjord conditions. The submarine melt rate _ m is related to the plume variables T p , S p , and u and the ice-ocean boundary layer temperature T b and salinity S b using a threeequation formulation which has frequently been used in this setting [e.g., Holland and Jenkins, 1999]: An analytical solution to the plume equations can be obtained using a number of simplifications (supporting information). For example, neglecting the melt and friction terms in (1)-(4) (which is a good approximation for a high initial discharge) and assuming a linear equation of state and uniform ambient conditions, equations (1-4) have the solution described in Straneo and Cenedese [2015] for a point source plume. This analytical solution can provide physical intuition, but for the more general cases of finite source size,  Figure 1. Schematic of the coupled plume-ocean model used in this paper. T a , S a , q a 5 ambient (fjord) temperature, salinity, and density, T p , S p , q p 5 plume temperature, salinity, and density. (1) At each time step, ambient temperature and salinity profiles are sent from the ocean model to the plume model.
(2) These are used in conjunction with subglacial discharge to calculate the vertical fluxes of water, heat, and salt.
(3) These fluxes are then used to calculate source and sink terms in the ocean model. stratified ambient conditions, a nonlinear equation of state, and inclusion of the melt and friction terms, the equations must be solved numerically. In this study, we therefore solve the equations numerically using the DLSODE differential equation solver [Hindmarsh, 1983].
In common with the model of Jenkins [2011], the plume emerges vertically into the domain. This is unlikely to occur in reality, and thus it is not clear how to choose the initial plume radius b i and velocity u i for a given discharge Q, related in our model by Q5pb 2 i u i =2. We find however that the model is insensitive to this choice, in that solutions converge quickly on one another for various combinations of (b i u i ) so long as Q remains fixed. Accordingly, we maintain a constant u i of 1 m s 21 and vary b i to provide the required discharge.

Integration of the Plume and Ocean Models
To allow the simulation of plumes at subgrid scale, we embed the plume model as a new module within the MITgcm (Figure 1). In grid locations where subglacial runoff is specified, ambient conditions (T a , S a , q a ) are passed at each time step from the MITgcm to the plume model. These are used, in conjunction with the specified runoff (b i , u i , T i , S i ), to calculate plume radius, velocity, temperature, and salinity as described above. The vertical extent of the plume is also calculated, with the plume terminating upon reaching neutral buoyancy or the fjord surface. Water, heat, and salt are then removed from MITgcm cells at the rate predicted by the plume model in those vertical layers where the plume is entraining, and put into the cell at the depth at which the plume is predicted to terminate. Further detail on the model implementation is provided in the supporting information.
The primary advantage of this approach over models which seek to resolve the plume is that it permits a relatively coarse spatial and temporal resolution to be used throughout the domain, greatly increasing computational efficiency. Xu et al. [2013], for example, required grid cells of 1 m 3 at the ice front to resolve the plume, necessitating limiting their domain to a section of fjord 500 m long, 500 m deep, and 150 m wide. Such high spatial resolution also requires a very small time step, particularly since vertical plume velocities may exceed 2 m s 21 [Xu et al., 2013], which again increases the computational demands. Because the plume parameterization presented here undertakes vertical transport outside of the MITgcm framework, the time constraint imposed by high vertical velocities within the model domain is avoided. An additional benefit of this method is that it may improve the prediction of plume dynamics relative to intermediate resolution models (10-20 m horizontal resolution) which seek to resolve the plume, as such models remain highly dependent on viscosity and diffusivity parameterization terms within the ocean model that are difficult to constrain [Sciascia et al., 2013;Xu et al., 2012].
Along the remainder of the ice front, there is melting but no-runoff input. The melt rate in these cells is calculated using the same formulation as in the plume model (equations (5-7)) using the temperature, salinity, and velocity as determined by the ocean model. The impact of this melting is taken into account by generating a virtual salinity and heat flux out of the domain, thus cooling and freshening the cells, based on the approach used to incorporate the melting of glacial ice in previous studies using MITgcm [e.g., Losch, 2008;Xu et al., 2012]. In reality, this melting would likely establish convective cells against the ice front on spatial scales too small to feasibly resolve in a fjord model [Wells and Worster, 2008]. Although these convective cells are expected to be much weaker than those forced by subglacial runoff (with a vertical extent on the order of 10 m in a stratified ocean environment [Huppert and Josberger, 1980;Jacobs et al., 1981]), it is important that these currents are taken into consideration. This is because they control the melt rate on parts of the calving front not in contact with a runoff-driven plume, which may be the majority of the ice front. As such, in some experiments we specify a minimum velocity in equations (5-7) (similar to Xu et al. [2013]). Experiments were undertaken using a minimum velocity of 0.01, 0.06, and 0.10 m s 21 , based on a range of modeling and field experiments, to assess the influence of this value on the processes simulated in this study (supporting information).

Initial Fjord Conditions
We utilize temperature and salinity profiles taken in front of Store Glacier, west Greenland, during August 2010 [Chauch e et al., 2014] as initial and ocean-boundary conditions for all experiments presented in this paper ( Figure 2a). These records, chosen because they provide a detailed survey of temperature and salinity within a few kilometers of the calving front of a major Greenlandic outlet glacier, were averaged to provide representative near-glacier conditions. Although in reality conditions down-fjord may vary from those recorded near to the calving front, we specify horizontally uniform initial conditions. This means that no model spin-up period is required, with initial velocities throughout the domain zero and any subsequent circulation resulting from the evolution of water properties over the course of the simulation.

Plume Model
As a first step, we examine the plume model independently of the ocean model. Plumes are calculated based on subglacial discharges of 1, 10, and 100 m 3 s 21 , using an initial plume temperature of 0 C and salinity of 0 psu and ambient fjord conditions as described in section 2.2. We compare the results to that of Jenkins's [2011] model, run using the same initial plume conditions, to assess the difference in the model solution generated by the different plume geometries. The horizontal extent of the plume is described in terms of the plume radius in the half-conical model and the thickness of the plume perpendicular to the ice front in Jenkins's [2011] model (because of the differing geometry); however, these two variables are directly comparable in terms of the expansion of the plume as it rises through the water column.
The results of the half-conical plume model, along with those from the model by Jenkins [2011], are shown in Figure 3. The two models produce similar plume behavior: the plumes expand as they rise, entraining salty water until they become denser than the ambient waters (due to the salinity stratification), which occurs at greater depth for a smaller initial discharge. Because the water column is warmest at depth, both models predict the plumes to be warmer than the ambient conditions above 350 m depth. The principal difference between the two models is that our model entrains more rapidly. This is because the half-conical form gives the plume a greater surface area across which to entrain relative to the sheet form of

Ambient
10 m main implication of this therefore is that plumes are likely to reach neutral buoyancy somewhat lower in the water column than predicted by Jenkins [2011].

Coupled Plume-Ocean Model
The next step is to examine the representation of one of these plumes (subglacial discharge 5 100 m 3 s 21 ) in the coupled plume-ocean model. The focus of this experiment is across-fjord variation in conditions at the fjord-glacier boundary. To do this, we use a rectangular domain 6.7 km long, 1.8 km wide, and 500 m deep. Grid resolution is 200 m in the cross-fjord direction (chosen such that the width of a cell is comparable to the predicted maximum width of the plume) and 10 m in the vertical. In the along-fjord direction, model resolution is 200 m for the 2 km of the domain closest to the glacier, decreasing linearly to 1500 m at the fjord mouth. The sides of the domain are closed, while a glacier extends the full width of the domain at the fjord head, with the plume in the center. A relaxation zone is used to restore properties to the initial conditions at the ocean boundary. The model is run for a simulation time of 1 h, during which time currents generated at the glacier do not reach the ocean boundary.
Initial conditions throughout the domain are set based on the temperature and salinity profiles from Store Glacier (Figure 2a). Diffusivity and viscosity parameters are chosen to be as low as possible without generating numerical noise [Kantha and Clayson, 2000] and are constant between experiments (Table 1; see supporting information for assessment of the sensitivity to these and other parameters). In order to assess the impact of cross-fjord circulatory patterns on melt rate, we do not specify a minimum background velocity in this experiment (section 2.1.3).
The interaction between the plume and fjord is best seen by looking at the modification of water properties across the ice front ( Figure 4). As in the plume model output (Figure 3), the plume forms a zone of relatively fresh, warm water rising up the center of the ice front (Figures 4a and 4b). As it rises, the plume entrains water from the fjord, generating relatively weak (0.01 m s 21 ) currents flowing toward the center of the ice front. Where the plume reaches neutral buoyancy, this pattern is reversed, with somewhat stronger currents (up to 0.07 m s 21 ) flowing away from the center of the ice front ( Figure 4d). These currents are however dwarfed by the vertical motion occurring within the plume itself, where velocities exceed 2 m s 21 . Consequently, while melt rate does increase toward the center of the fjord, this variation is extremely small compared to the melting taking place where the plume contacts the ice front (Figures 4e and 4f). Furthermore, these patterns are only visible because there is no specified minimum velocity used in the melt calculation in this scenario, and so melt rate is negligible except where currents are stimulated by the plume. If in reality other processes such as melt-driven convection, tides, and wind-forced circulation generate additional currents parallel to the ice front, then it is likely that the significance of plume-induced cross-fjord circulation on melt rate will be further diminished. Although the relatively coarse resolution of the model may dampen currents in the area around the plume, we note similar findings from higher-resolution simulations, which suggest extremely low melt rates except where the glacier is in direct contact with the plume [Xu et al., 2013].

Impact of Runoff on Fjord Circulation
In this experiment, we seek to examine the impact of glacial melt and runoff on the circulation of a fjord over the course of a year. For simplicity, we use a rectangular fjord of constant width and depth. Model resolution is 200 m in the cross-fjord direction and 10 m in the vertical, with a domain width and depth of 5 km and 500 m, respectively (chosen as appropriate dimensions for a large Greenlandic tidewater glacier). Resolution in the along-fjord direction is 200 m at the ice front, decreasing linearly to 2400 m at the fjord mouth, with a domain length of 130 km. The fjord is considered to comprise 100 km of this length, while the remaining 30 km represents a coastal strip. In this coastal zone, a relaxation zone restores water properties to the initial conditions-in this way, coastal properties remain constant, with all variability driven at the ice front [e.g., Sciascia et al., 2013]. The Coriolis parameter is set to 1.37 3 10 24 , appropriate for a typical Greenlandic latitude of 70 N [Kantha and Clayson, 2000].
Three experiments are undertaken, running for 12 months from January to January. In the first, the head of the fjord is simply a rock wall, with no melting and no input of runoff. In the second, there is an ice wall at the head of the fjord, allowing melting, but there is no runoff. In the third experiment, runoff enters from the base of the glacier in the center of the fjord. We use an idealized runoff input, interpolating linearly between monthly values. These values are based on a normal distribution centered on a mean July value of 300 m 3 s 21 (Figures 2b and 2c). While runoff is expected to vary significantly between glaciers and years (meaning there is no representative value for a typical Greenlandic fjord), this is chosen as an appropriate peak runoff as predicted for outlet glaciers such as Helheim Glacier [e.g., Mernild and Liston, 2012;Sciascia et al., 2013] and Store Glacier [e.g., Xu et al., 2012]. The effect of varying the duration and intensity of the melt season, and hence total annual runoff, is examined in the following experiment (section 3.4).
Each of these scenarios was run using three different values of background velocity (0.01, 0.06, and 0.10 m s 21 ). The choice of background velocity affects the predicted submarine melt rate (as described in section 3.4), but has little impact on the overall circulation (supporting information). As such, only those scenarios with the moderate background velocity (0.06 m s 21 ) are described in this section (with additional details in section 3.4). The velocity fields generated in these experiments are approximately symmetrical across the fjord (with along-fjord velocities increasing toward the fjord center) and so presentation and discussion of results is focused on properties along the fjord centerline (across-fjord variation in velocity is illustrated in the supporting information). Figure 5 shows the along-fjord velocity patterns produced in these simulations. In the absence of a glacier, there is little movement: only weak currents near the ocean boundary generated by horizontal salinity gradients forming in response to gradual vertical diffusion (Figure 5ai). The addition of an ice front, but no runoff, generates a series of vertical circulation cells, driven by melting at the fjord head (Figure 5bi). These currents remain very weak, however, with maximum predicted velocities of the order of 10 23 m s 21 . The addition of runoff stimulates a much stronger circulation, with a clear seasonality (Figures 5ci-fi). As runoff increases during the spring, a down-fjord current forms at the depth at which the plume reaches neutral buoyancy, compensated by a weaker but thicker up-fjord current below (Figure 5ci). The depth of the down-fjord current decreases and its strength increases as runoff increases, until by peak runoff it occupies the upper 100 m of the water column with a maximum velocity of the order of 0.1 m s 21 (Figure 5di). At this stage, the fjord resembles a single-cell estuarine circulation. The reverse of this process happens during autumn, with the down-fjord current increasing in depth and decreasing in strength as runoff decreases (Figure 5ei).
In the no-glacier scenario, temperature change is limited to the vertical smoothing of the temperature profile by diffusion (Figure 5aii). The addition of the ice front generates cooling near the fjord head, most notably in the lower part of the water column where the warmer waters are cooled by 0.1 to 0.2 C (Figure 5bii). A similar cooling occurs in the runoff scenario, but greater change to the temperature structure occurs due to the vertical movement of water within the plume (Figures 5cii-5fii). In keeping with the results of the plume model experiment (Figure 3), the waters exiting the plume at the depth of neutral buoyancy tend to be warmer than the ambient conditions at that depth, particularly when the plume terminates in the cold band between approximately 80 and 180 m depth. This generates a general warming of the upper part of the water column, with the warmest layer in the upper part of the water column often coinciding with the depth of the downfjord current at that time ( Figure 5).
Changes to the salinity distribution in the fjord are less notable (Figures 5aii-5fii). Outflow from the plume inputs a large quantity of water of constant salinity at some depth in the water column. This water may pileup near the ice front, generating a slight along-fjord salinity gradient. As runoff increases and the depth of the plume outflow decreases, the up-fjord currents generated by plume entrainment serve to restore the original salinity distribution.

Impact of Runoff on Submarine Melt Rate
Theoretical and numerical modeling experiments have shown that submarine melt rate increases with the discharge of subglacial runoff from a single channel [e.g., Jenkins, 2011;Sciascia et al., 2013;Xu et al., 2013Xu et al., , 2012. However, the significance of this relationship when considered across the entire ice front and on annual and interannual time scales is not well understood. To investigate this, we conducted a series of 1 year experiments in which the magnitude of the summer melt season was varied. Taking the seasonal runoff distribution used in section 3.3 (7 months long, peaking at 300 m 3 s 21 in July) as a base case ''standard'' year, we conducted a series of experiments with total annual runoff scaled to 0.5, 1.5, and 2 times this value, representing a generous allowance for interannual variability. (For comparison, Mernild et al. [2010] estimate interannual variability in total annual runoff from Helheim Glacier to reach a maximum of approximately 625% of the mean over the period 1999-2008.) In the first set of experiments, we achieve this by varying the intensity of melt season, such that its duration remains constant but the magnitude of runoff is scaled accordingly. In the second set of experiments, we keep the peak July runoff at 300 m 3 s 21 , and instead vary the duration of the melt season, maintaining a normal distribution. The resulting runoff time series are shown in Figures 2b and 2c.
As seen in Figure 4, the average ice front melt rate can be broken into two components: ''plume melting,'' where the plume comes into contact with the glacier, and ''background melting'' where it does not. The former component scales with discharge, causing average melt rate to rise toward the middle of the melt season then fall again (Figures 6a and 6b). If only the central 200 m of the ice front are considered, plume melting dominates and we find a similar power law relation between melt and discharge to that identified in previous studies [ when the entire width of the ice front is considered, the background melt rate becomes much more significant. As in Figure 4, we do not find that the plume generates sufficient cross-fjord circulation to strongly influence this background melt rate, although the existence of relatively warm glacially modified waters leads to a slightly higher melt rate in the autumn relative to spring (Figures 6a and 6b). Instead, we find that the background melt rate is strongly dependent on the minimum velocity value prescribed in equation (6) (see also supporting information). This is problematic, as this value remains somewhat artificial and is difficult to constrain.
Increasing the intensity of runoff results in higher peak summer melt rates (Figures 6a and 6b). However, as identified here (Figure 6c) and in previous studies [Jenkins, 2011;Sciascia et al., 2013;Xu et al., 2012], the sensitivity of plume melt rate to discharge decreases as discharge rises (Figure 6c), meaning the difference between scenarios is somewhat subdued. For example, doubling peak runoff from 300 to 600 m 3 s 21 results in a 25% increase in peak melt rate relative to the background melt rate (Figure 6a). Increasing the duration of the melt season has a greater effect for the same increase in annual runoff, because plume melt rate is more sensitive to variation in discharge at lower discharge values (Figure 6c). This effect can be seen most clearly when the duration of the melt season is altered while the total annual runoff is kept constant (Figure 7). Although the more intense melt seasons generate higher peak melt rates, this is more than compensated by the greater duration of enhanced melt rates in the longer melt seasons (Figure 7).
The significance of annual and interannual variability in runoff on total melt rate depends on the background melt rate, which in turn depends on the prescribed background velocity (Figures 6d-6f). In the scenario with the smallest background melt rate (Figure 6d), doubling total runoff from the standard scenario results in a 15-45% increase in total annual melting. In the moderate background melt rate scenario, this drops to a 3-12% increase (Figure 6e), while in the highest background melt rate scenario, this falls further to 2-8% (Figure 6f). This behavior is discussed further in section 4.
Following this, the model was run for 10 consecutive years to examine whether there was a cumulative impact on melt rate over this time (for example, due to a sustained cooling near the ice front). Two scenarios were run-in the first, there was no runoff; in the second, the standard runoff forcing was repeated for 10 cycles. In order to reduce the computational time required, the horizontal model resolution was increased to a uniform 1000 m in along and across-fjord directions. The effect of this coarser resolution is assessed in Figure 8a, which shows the submarine melt rate over the course of 1 year, as predicted using the original horizontal model resolution and the coarser grid. During the first part of the year, the predicted melt rates in the two simulations are almost identical, while in the second part of the year, the predicted melt rate is slightly reduced in the coarser domain. The overall similarity of results demonstrates the importance of the plume melting (which is independent of the model resolution) on the total melt rate and justifies the use of the coarser-resolution domain for these experiments (see also supporting information). The results of the 10 year simulations are shown in Figure 8b. There is a small increase in melt rate following the first year in the standard runoff scenario due to the displacement of cold water from adjacent to the glacier by warmer glacially modified water. With the exception of this, there is however little change in melt rate in either scenario over the 10 years. This suggests that, with or without runoff, there is no significant cumulative cooling or warming of the fjord waters surrounding the glacier front on this time scale (assuming no other external change in properties of the fjord water entering the fjord mouth).

Glacial Modification of Fjord Oceanography
Oceanographic survey data from Greenland's fjords are sparse and frequently limited to late summer when ice cover in the fjords is at its annual minimum. The area close to the ice front, often choked with icebergs and sea ice and at risk from calving events, is extremely difficult to survey, making data from this crucial location in which glacially modified water masses are formed particularly rare . This makes interpretation of existing data sets and model results challenging, limiting our understanding of the impact of glaciers on fjord oceanography.
Nevertheless, several authors have used temperature and salinity characteristics to argue for the existence of glacially modified layers in fjord survey data. At Sermilik Fjord, east Greenland, Straneo et al. [2011] identified a relatively warm layer extending over 50 km from the ice front at a depth of between 100 and 200 m. They argued that this layer was formed by subglacial runoff and entrained warm Atlantic waters, rising upward until it reached neutral buoyancy at the steep halocline between the Atlantic waters and the overlying cooler, fresher Polar waters. Similarly at Store Glacier, Chauch e et al. [2014] reported the presence of warm, glacially modified layers in the upper part of the water column. In August 2010, this layer occurred at the surface, above the coldest water, while during the cooler August of 2009 the glacially modified waters were found both at the surface and lower in the water column. This occurred, they argued, because reduced runoff meant that the plume reached neutral buoyancy lower in the water column during the cooler year.
Our findings support and add further detail to these hypotheses. We find that as runoff increases through the spring, the plume grows in strength and rises higher in the water column. As it does so, the waters near the ice front are replaced with a mixture of entrained fjord waters, glacial runoff, and the products of submarine melt. Of these, the fjord waters are dominant and the submarine meltwater is negligible. For example, by the point of neutral buoyancy in the highest discharge scenario in Figure 3, the plume is transporting >6000 m 3 s 21 of entrained fjord water, 100 m 3 s 21 of glacial runoff, and 1.5 m 3 s 21 of freshwater from submarine melting. If the fjord was vertically uniform in temperature, this mixture would be slightly cooler than the ambient waters. However, because deep silled Greenlandic fjords are characterized by a warmer, saltier layer of Atlantic waters underlying a cooler, fresher layer of Polar waters [Straneo and Heimbach, 2013], this glacially modified water is typically warmer than the ambient waters at the depth at which it reaches neutral buoyancy.
The plume outflow sets up a circulatory cell, with a thinner, stronger down-fjord current overlying a thicker, weaker, up-fjord current. Maximum modeled along-fjord velocities are on the order of 0.01-0.1 m s 21 , in agreement with the residual runoff-driven circulation identified by Sutherland and Straneo [2012] in Sermilik Fjord, east Greenland. Glacially modified waters are reentrained by this cell as the outlet point of the plume rises during the spring, and so in summer there may be little evidence of this modification below the depth of the main plume outflow (e.g., Figure 5d). This is likely to be the scenario observed by Chauch e et al.
[2014] during August 2010, when the only runoff-modified waters are in the surface layer. As runoff decreases in the autumn, the plume loses strength and drops down through the water column again (Figure 5e). Glacially modified waters become trapped in the relatively slow flowing zone above the level of the main circulatory cell as the depth of the down-fjord current increases during this phase. These glacially modified waters may form a warm zone around the glacier front, as observed by Chauch e et al. Not all fjord survey data fit this picture. In Kangerdlugssuaq Fjord, east Greenland, Inall et al. [2014] report cold glacially modified waters, while warm surface waters originate on the shelf, the opposite of the scenario reported at Store Glacier [Chauch e et al., 2014]. This may be because an extensive ice m elange exists in front of Kangerdlugssuaq Glacier, acting as a substantial additional heat sink [Enderlin and Hamilton, 2014]. This demonstrates the complexity of the glacial fjord system, with multiple freshwater inputs (including subglacial and surface runoff, as well as meltwater from icebergs and glaciers) interacting with complex and evolving water masses originating beyond the fjord mouth. In this study, we are focusing on a simplified system, in which only subglacial runoff and melting from a single glacier act upon the fjord. As such, the features that we identify here will in reality be but one component of more complex and variable circulation [Sutherland and Straneo, 2012].
A further level of complexity that may be present in real-world scenarios is across-fjord variability in circulation and fjord water properties, introduced by both topographic complexity and the effects of rotation. The Rossby radius of deformation of our idealized fjord is 9 km, which is large compared to the fjord width of 5 km (Kelvin number 5 0.55) [Garvine, 1995]. As such, the modeled fjord is dynamically narrow and rotation is not important. Large Greenlandic fjords such as Sermilik Fjord and Kangerdlugssuaq Fjord are typically slightly wider than the 5 km width used here, and the Kelvin number can equal or exceed unity in places . Sutherland et al. [2014] find geostrophic currents to be of importance offshore of the fjord mouths, but report only limited evidence of lateral gradients within the fjords themselves. Contrastingly, Inall et al. [2014] report significant across-fjord gradients midway along Kangerdlugssuaq Fjord, arguing that circulation within the fjord is in geostrophic balance. It is therefore possible that in Greenland's largest (and in particular broadest) fjords, rotation may become significant, generating greater across-fjord variation in circulation patterns than is apparent in the experiments undertaken here (supporting information).

Interaction of Runoff, Fjord Circulation, and Submarine Melt Rate
There are several ways in which variability in runoff could influence glacier melt rate. The first, as investigated in previous studies [e.g., Jenkins, 2011;Xu et al., 2012], is that increased runoff generates a more vigorous upwelling, and as such generates higher melt rates where this plume comes into contact with the ice front. The second mechanism would be if this plume generated stronger currents across a wider part of the ice front, causing an increase in melt rate extending beyond the area in contact with the plume. The third mechanism is that the along-fjord circulation driven by the runoff could help to supply warm water to the fjord head.
In keeping with previous studies [Jenkins, 2011;Sciascia et al., 2013;Xu et al., 2013Xu et al., , 2012, we find that melt rate is strongly dependent on runoff, when measured across the region immediately surrounding the plume (Figure 6c). A single plume covers only a small proportion of the ice front however, and when averaged across the remainder of the glacier, the effect is reduced accordingly. We do not find that runoff generates significant cross-fjord circulation away from the plume (Figure 4). We also do not observe a strong difference in along-fjord temperature gradient between the runoff and no-runoff scenarios, with melt-driven circulation sufficient to prevent sustained cooling of waters adjacent to the ice front even in the no-runoff scenario ( Figure 5). A similar cooling at depth close to the ice front is observed in the runoff and no-runoff scenarios even after 10 years, although there is a modest (<3%) increase in melt rate in the runoff scenario due to the displacement of cold water masses in the upper part of the water column by glacially modified waters (Figure 8). This situation may be different if the temperature of the water mass at the fjord mouth changed over time, as the glacially forced circulation could serve to accelerate the rate at which these changes affected the properties of the water in contact with the glacier. Further work is required to address this uncertainty, in particular examining the relative importance of glacially and coastally forced fjord shelf exchange in a range of geographical settings [e.g., Jackson et al., 2014;Sciascia et al., 2014;Sutherland et al., 2014].
The cumulative effect of these processes is that average melt rate increases in the summer due to more rapid melting where the ice is in contact with the plume and remains slightly elevated the following winter due to the displacement of cold water masses ( Figure 6). Whether or not these processes result in a significant increase in the total mass loss through submarine melting depends on the background melt rate, which remains comparatively steady across large sections of the ice front over the course of the year. Melt rate is strongly dependent on the velocity at the ice-ocean interface [Jenkins et al., 2010], and so cannot be properly predicted without knowledge of these currents (supporting information). However, resolving buoyant convection due to melting on this boundary would require an extremely fine mesh in this zone, which is computationally impractical. It may be possible to include this process through the use of an additional coupled model similar to our parameterization of the main runoff-forced plume, but this would require the development of an appropriate melt-driven convection model, which lies beyond the scope of this paper. Alternatively, a simple parameterization based on observation of submarine melt rates at tidewater glaciers would be invaluable, but obtaining the data required to formulate such a parameterization remains extremely challenging in this hostile environment .
If the background melt rate is in reality very low, then variability in runoff may be responsible for large relative changes in melt rate (Figure 6d). However, it is difficult to reconcile low overall melt rates with those derived from field-based heat flux estimates which place typical submarine melt rates at Greenland's glaciers in the range of 1-10 m d 21 [Inall et al., 2014;Rignot et al., 2010;Sutherland and Straneo, 2012], which would imply melt rate away from the main plume is high and contributes significantly to the total melt rate. If this is the case, then runoff variability, particularly on the scale of interannual variability, may not result in substantial variability in submarine melt rate. There are however large uncertainties in the melt calculations used in both modeling and observational studies. The melt parameterizations used here and in other models of submarine melting [Sciascia et al., 2013;Xu et al., 2013Xu et al., , 2012 remain poorly validated for the vertical faces and strong currents of tidewater glacier fronts, having been developed and tested based on heat transfers beneath ice shelves and sea ice [Holland and Jenkins, 1999]. Field-based melt rate estimations [Inall et al., 2014;Rignot et al., 2010;Sutherland and Straneo, 2012] on the other hand, are often based on heat flux estimates from relatively sparse observations acquired at some distance from the ice front, and do not account for the potentially considerable loss of heat through melting of the ice m elange [Enderlin and Hamilton, 2014]. As such, comparing field and model-based melt rate estimates remains problematic.
A further consideration is that runoff may not be concentrated into a single large channel when it enters the fjord. Although field observations seldom report more than one large plume [e.g., Chauch e et al., 2014; Sole et al., 2011], there may be numerous smaller plumes that never reach the fjord surface and so remain undetected. If runoff is shared between several smaller channels rather than (or in addition to) emerging from a single main source, the individual plumes will be weaker, and so the melting associated with each individual plume smaller. Nevertheless, average melt rate across the entire ice front may be increased as a greater proportion of the ice is affected by plumes [Kimura et al., 2014]. The basis of this result lies in the sublinear relationship between freshwater discharge and average plume melt rate, with melting scaling with discharge raised to the power of 1/3 to 1/2 [Jenkins, 2011;Sciascia et al., 2013;Xu et al., 2012]. This means that for a given total runoff, concentrating this runoff leads to a decrease in total melting. This reasoning also applies with respect to time, with a longer melt season generating greater total melting for a given total runoff (Figure 7).
The distribution of runoff has the potential to influence the seasonal variability in melt rate. A greater number of channels would result in a greater proportion of the ice front being influenced by seasonal runoff input, potentially increasing the submarine melt rate in summer relative to winter. Conversely, a distributed Journal of Geophysical Research: Oceans 10.1002/2014JC010324 drainage system may persist during the winter months (fed by basal frictional melting), reducing seasonal variability in submarine melt rate. It is clear that a more complete understanding of tidewater glacier hydrology is required to better constrain submarine melt rate, and how this varies in space and time .
A further factor that is likely to be important in this respect is the feedback between submarine melting and calving as the morphology of the calving front is modified by spatially uneven melting [O'Leary and Christoffersen, 2013]. If submarine melt is concentrated around a large channel then this may cut back into the glacier, causing enhanced calving in the surrounding region that spreads out across the ice front [Chauch e et al., 2014]. Alternatively, if runoff is more widely distributed across the grounding line, greater melting at depth may undercut the ice front [O'Leary and Christoffersen, 2013]. Better understanding of the links between submarine melting and calving mechanics is therefore necessary in order to assess whether the average ice front melt rate, or the spatial distribution of submarine melting, is more important in controlling tidewater glacier mass balance.

Comparison With Previous Numerical Modeling Studies
This paper represents a first attempt to model the impact of glacial runoff on fjord oceanography by coupling a theoretical plume routine to an ocean circulation model. In recent years however, several studies have used MITgcm to model the formation of plumes arising from subglacial discharge by using a highresolution set up in which the plume turbulence is partially resolved [Sciascia et al., 2013;Xu et al., 2013Xu et al., , 2012. Several key elements are consistent between this study, previous numerical modeling studies and field surveys: subglacial runoff generates a plume that rises adjacent to the glacier, producing a glacially modified down-fjord current underlain by an up-fjord current. One feature of field surveys that modeling studies have sought to reproduce is the observation that glacially modified currents may form at an intermediate depth between the surface and bed [Chauch e et al., 2014;Straneo et al., 2011], in contrast to a classic estuarine circulation in which the outflow occurs at the surface [e.g., Motyka et al., 2013]. While previous modeling studies have succeeded in generating glacial plumes which reach neutral buoyancy and flow down-fjord at some intermediate depth, this has only occurred at discharges less than 5 m 3 s 21 [Sciascia et al., 2013;Xu et al., 2012] to 30 m 3 s 21 [Xu et al., 2013]. These values are very low compared to the expected summer runoff from these glaciers (hundreds of cubic meters per second), indicating that runoffmodified waters should be present at the surface throughout much of the year.
In the case of Xu et al. [2012] and Sciascia et al. [2013], this discrepancy occurs principally because the domains used in these earlier studies were two dimensional. This means that the surface area across which fjord water can be entrained into the plume is greatly reduced (i.e., there is no entrainment in an acrossfjord direction), allowing the plume to remain buoyant higher into the water column. This case is similar to the plume model of Jenkins [2011], in which entrainment is only permitted on the down-fjord face of the plume, allowing the plume to rise higher into the water column relative to the half-conical plume geometry that we employ in this study (Figure 3).
In the three-dimensional domain of Xu et al. [2013] (which uses a similar Store Glacier temperature and salinity stratigraphy to this paper), the plume reaches the surface at a discharge of less than 30 m 3 s 21 . However, the glacially modified waters discharged by the plume remain saline relative to the surface waters and quickly sink back to a depth in excess of 100 m [Xu et al., 2013, Figure 2c]. Similarly, we find that if momentum is allowed to carry the plume past the point of neutral buoyancy in our model, then the plume will reach the surface at a subglacial discharge of 50 m 3 s 21 (supporting information). As found by Xu et al. [2013], this water would however be denser than its surroundings and sink back to the depth of neutral buoyancy (below 100 m) as it moves away from the upwelling. To avoid this density inversion, we terminate the plume at the depth of neutral buoyancy. This means that while in reality the plume may break the surface before sinking back to a depth of neutral buoyancy (a process observed at Store Glacier by Chauch e et al. [2014]), our simulated plume is not visible at the surface until the glacially modified waters discharged by the plume are of density equal to or less than the surface waters. This requires a subglacial discharge of several hundred cubic meters per second ( Figure 5 and supporting information). A wider implication of this mechanism is that while plumes may become visible at the surface at relatively low discharges (tens of cubic meters per second), much greater discharges (hundreds of cubic meters per second) may be required to generate runoff-modified waters that flow down-fjord at the fjord surface.

Conclusions
In this paper, we have used a plume model to parameterize the buoyant upwelling of subglacial runoff in an ocean circulation model. This greatly improves the computational efficiency of modeling glacial fjords, permitting simulation of larger domains (hundreds of square kilometers in area) over longer time periods (years to decades) than has previously been possible. This model has been used to examine the effect of runoff on fjord circulation and submarine melt rate independent of coastal forcing. In these simulations, runoff generates a relatively thin, strong, down-fjord current overlying a relatively thick, slow up-fjord current. The strength and depth of the down-fjord current varies with runoff on annual time scales, being strongest, and forming highest in the water column, when runoff is greatest. Because this current contains entrained Atlantic waters from the lower part of the water column, it is typically warmer than the ambient conditions at the depth at which it forms. This finding supports the identification of anomalously warm, glacially modified waters in some Greenlandic fjords [Chauch e et al., 2014;Straneo et al., 2011].
Submarine melt rate increases strongly with runoff in areas where ice is in contact with a buoyant runoffdriven plume. This contact area may be small however relative to the width of the tidewater glacier margin, and the impact is greatly reduced when averaged across the entire ice front. We find the submarine melt rate to be only weakly sensitive to runoff when considered on the scale of interannual variability and across the full width of the glacier. The significance of the relationship between runoff and submarine melt rate is however dependent on the distribution of runoff across the grounding line and the strength of the background melt rate in those areas not in contact with the plume, both of which remain poorly constrained. In recent years, there has been much interest in the relationship between melt rate and discharge in the vicinity of ice marginal plumes [e.g., Jenkins, 2011;Sciascia et al., 2013;Xu et al., 2013]; further work is needed to better constrain both the proportion of the ice front directly affected by runoff-driven plumes and the submarine melt rate across the remainder in order to improve the prediction of melt rate across the full width of the glacier front. Furthermore, this must be combined with a better knowledge of how the rate and distribution of melting influences calving mechanics before the significance of submarine melting for the mass balance of tidewater glaciers can be accurately assessed.