Freshwater produced by the surface melting of ice sheets is commonly discharged into ocean fjords from the bottom of deep fjord-terminating glaciers. The discharge of the freshwater forms upwelling plumes in front of the glacier calving face. This study simulates the meltwater plumes emanated into an unstratified environment using a nonhydrostatic ocean model with an unstructured mesh and subgrid-scale mixing calibrated by comparison to established plume theory. The presence of an ice face reduces the entrainment of seawater into the meltwater plumes, so the plumes remain attached to the ice front, in contrast to previous simple models. Ice melting increases with height above the discharge, also in contrast to some simple models, and the authors speculate that this “overcutting” may contribute to the tendency of icebergs to topple inwards toward the ice face upon calving. The overall melt rate is found to increase with discharge flux only up to a critical value, which depends on the channel size. The melt rate is not a simple function of the subglacial discharge flux, as assumed by many previous studies. For a given discharge flux, the geometry of the plume source also significantly affects the melting, with higher melt rates obtained for a thinner, wider source. In a wider channel, two plumes are emanated near the source and these plumes eventually coalesce. Such merged meltwater plumes ascend faster and increase the maximum melt rate near the center of the channel. The melt rate per unit discharge decreases as the subglacial system becomes more channelized.
The Greenland Ice Sheet’s contribution to sea level rise has increased significantly in the past decade (Rignot and Kanagaratnam 2006; Velicogna 2009; Shepherd et al. 2012). Much of this change is due to accelerated flow of outlet glaciers into deep fjords in western and southeastern Greenland (Howat et al. 2007; Joughin et al. 2004; Stearns and Hamilton 2007). In Greenland, the melting of the ice sheet surface produces freshwater that percolates down to the bed through moulins, imparting a seasonal variability to the motion of ice in the ablation zone (Zwally et al. 2002; Das et al. 2008; Joughin et al. 2008). The percolated freshwater forms a hydrological network beneath the ice sheet (Zwally et al. 2002; Bartholomaus et al. 2008; Fountain and Walder 1998; Bartholomew et al. 2010) and eventually drains to the ocean beneath the termini of the glacier, which is typically a vertical ice face in one of the Greenland fjords. Chu et al. (2009) observed the appearance of meltwater plumes at the termini of the glacier in Kangerlussuaq Fjord soon after the initial onset of surface melting, indicating a rapid coupling between the ice sheet and subglacial discharge. The subglacial discharge has strong seasonal variation, with peak discharge during the melting season, and can flux more freshwater than the melting at the ice–ocean interface (Matthews and Quinlan 1975; Walters et al. 1988; Svendsen et al. 2002; Motyka et al. 2003).
Motyka et al. (2003) concluded that the thermal driving by the ocean alone could not explain the rapid melting of LeConte Glacier in Alaska and argued that the meltwater plume must be enhancing the heat transfer. In Greenland, the water properties in the fjords near the glaciers suggest that the water originated from subglacial discharge and Atlantic Water, a relatively warm water mass that intrudes to Greenland’s western and eastern shelves (Holland et al. 2008; Motyka et al. 2011; Straneo et al. 2012). These observations lend support to the idea that the meltwater plume entrains the Atlantic Water toward the ice face, increasing submarine melting at the termini of the glacier. The discharged water often contains high concentrations of glacially eroded basal sediments and iron from the mechanical and chemical weathering beneath glaciers (Drewry 1986). The discharge of sediments and iron play an important role in shaping the landscape (Powell 1990) and marine ecosystem (Hartley and Dunbar 1938; Apollonio 1973; Arimitsu et al. 2012) in the fjords.
Earlier studies speculated that the trajectory of meltwater plumes is determined by the relative strength of buoyancy and inertia of the discharged water (Syvitski 1989; Powell 1990). Under sufficiently large discharge, Syvitski (1989) argued that horizontal inertia could push the trajectory of a meltwater plume away from the glacier, and the plume would surface some distance away from the terminus. In the case of lower discharge, the buoyancy force was thought to dominate the inertia, and the meltwater plume is expected to rise along the ice face (Powell 1990; Cowan and Powell 1990).
Applications of plume theory in inclined coordinate systems have provided many insights into the overturning circulation beneath ice shelves in Antarctica, rooted in the pressure dependence of the freezing point, and its relationship to the melting of the ice (e.g., MacAyeal 1985; Jenkins 1991). The plume theory consists of a simple one-dimensional model of convection driven from a point source by casting volume, momentum, and buoyancy conservation into three ordinary differential equations (e.g., Morton et al. 1956). A common assumption in the subice shelf plume theory is that the only source of buoyancy that acts to stratify the water column and maintain the overturning circulation is the generation of meltwater at the ice–ocean interface. In the presence of subglacial discharge, Jenkins (2011) predicted that the melt rate scales with the cube root of the freshwater flux, implying slightly more than a doubling of the melt rate for an order of magnitude increase in flux. The first geophysical-scale simulation of a meltwater plume in an idealized two-dimensional fjord by Salcedo-Castro et al. (2011) found that the plume trajectory is attached to the ice face, regardless of the inflow velocity of subglacial discharge and in the absence of the generation of meltwater at the ice–ocean interface. The buoyancy feedback due to ice melting is included in the two-dimensional simulation of Xu et al. (2012) and Sciascia et al. (2013). Both studies confirmed that the trajectory of the meltwater plume is attached to the ice and found the cube root dependence of the melt rate on the subglacial flux predicted by Jenkins (2011). Xu et al. (2013) performed three-dimensional simulations using temperature and salinity profiles from Store Glacier in west Greenland as initial conditions. Their simulations use one channel size and show that the melt rate increases proportionally to thermal forcing and subglacial discharge flux to the power of 1.2–1.6 and 0.5–0.9, respectively, suggesting the role of a warming ocean water in melting the glacier.
The objective of this study is to extend our understanding of meltwater plumes and their role in melting glaciers. As the ice sheet melting season progresses, the subglacial system becomes more channelized. So far no studies have considered the sensitivity of the melt rate to channel geometry or the effect of multiple interacting channels. We also consider the question of detachment of plumes from the ice face. We employ a nonhydrostatic, unstructured, mesh ocean model, Fluidity (Piggott et al. 2008). The model configuration is described in section 2. Section 3 describes the results, including the model uncertainty in entrainment, the variation in melting with respect to source discharge and geometry, and a comparison with established plume theory. Finally, our discussion and conclusions are in section 4.
2. Model configuration
We employ Fluidity (Piggott et al. 2008) to integrate the nonrotating Boussinesq equations, and advection–diffusion equations for temperature and salinity, cast in a Cartesian coordinate system (x, y, z). The relative importance of the ambient rotation rate Ω to the time taken by a meltwater plume particle to cover distance L at speed U can be quantified by the Rossby number, defined as R = 2πU/ΩL, where Ω = 7.29 × 10−5 s−1. If R is on the order of or less than unity, the rotation is important relative to the time scale of plume evolution. In our simulation, the meltwater plume emanates ~500 m below the sea surface with a minimum speed of 0.1 m s−1. This gives us R ~ 17, much larger than unity. We therefore expect that the rotation of Earth has a minor effect on the evolution of meltwater plume at this scale, and the Coriolis force is not included in the equations. An implicit time-stepping scheme is used with a time step of 1 s. This short time step combined with an expensive nonhydrostatic configuration to simulate meltwater plumes limits our computational ability to include other forcing that acts on longer time scales, such as remote buoyancy or wind forcing. Velocity and pressure are discretized using a combination of first-order discontinuous Galerkin (DG) (P1DG) and second-order continuous Galerkin (P2) function spaces (a P1DG–P2 finite element pair), as described by Cotter et al. (2009). The advection–diffusion equations for temperature and salinity are discretized using P1DG. The equation of state is represented using an algorithm described in McDougall et al. (2003).
We consider three different idealized domain configurations A, B, and C (Fig. 1). All the configurations have a width of 1000 m in the x direction, vertical z length of 500 m, and a channel to prescribe the inflow at y = 0. Configurations A and B are of length 4400 m in the y direction, and the southern boundary is located away from the channel (y = −2000 m). Configuration C is 2400 m in length and contains a vertical ice face at y = 0.
b. Ice–ocean thermodynamics
Glacial ice was incorporated into Fluidity by Kimura et al. (2013). The model has been extended to simulate frazil ice formation and deposition to better represent the freezing process (Jordan et al. 2014). The ice melt rate is calculated for the configuration C by applying thermodynamic boundary conditions at the ice–ocean interface. Fully resolving physical processes near the ice–ocean interface while simulating regional oceanic circulation is not computationally feasible, so we rely on a parameterization based on Kader and Yaglom (1972) to estimate turbulent heat and salt fluxes toward the ice–ocean interface in order to compute a melt rate. The melt rate is calculated by the “three-equation” formulation, which links the local freezing relation and balance of heat and salt fluxes at the ice–ocean interface (e.g., Hellmer and Olbers 1989; Holland and Jenkins 2001; Losch 2008). The local freezing relation constrains the temperature Tb and salinity Sb at the ice–ocean interface:
where a = −0.0573°C, b = 0.0832°C, and c = −7.53 × 10−8 °C Pa−1, and P is the local hydrostatic pressure. The balances of heat and salt fluxes between the ice and ocean are
where c0 = 3974 J kg−1 °C−1 and cI = 2009 J kg−1 °C−1 are the specific heat capacities of seawater and ice, respectively. The velocity of the ocean in the direction normal to the ice–ocean interface is represented by m′, and the melt rate of ice is m = ρ0m′/ρice, where ρ0 and ρice are the density of the ocean and ice, respectively. The variable L = 3.35 × 105 J kg−1 represents the latent heat of ice fusion. The far-field internal temperature of ice is assumed to be TI = −25°C.
Nondimensional coefficients of the heat and salt transfer through the boundary layer are represented by γT and γS, where these numbers are a function of the molecular Prandtl number (the ratio of viscosity to thermal diffusivity) and Schmidt number (the ratio of viscosity to saline diffusivity), respectively (Holland and Jenkins 1999). These values are independent of the thickness of the ice–ocean boundary layer, since they are scaled with the far-field velocity. The precise values of γT and γS are uncertain, such that the precise melt rates predicted by the model are subject to uncertainty (Dansereau et al. 2014). We have assumed γT = 1.05 × 10−3 and γS = 3.97 × 10−5 in accordance with the formula presented in Holland and Jenkins (1999).
The terms on the right-hand sides of (2) and (3) are a parameterization of the mixing of heat and salt toward the ice through the oceanic boundary layer. The “far-field” ocean temperature and salinity are represented by T∞ and S∞. The variable u∞ represents the speed of ocean flow oriented parallel to the ice, which is taken to be the source of turbulence that drives the mixing of heat and salt toward the ice. We assume the presence of a minimum background flow of 10−4 m s−1 in the calculation of melt rate, so u∞ = 10−4 m s−1 when ugrid < 10−4 m s−1, otherwise u∞ = ugrid. The variable ugrid is the far-field speed calculated by the model. The melting continues at a background rate whenever the flow rate drops below 10−4 m s−1. For given T∞, S∞, and u∞, the three unknowns, Tb, Sb, and m′, are solved by combining (1)–(3) to produce a quadratic equation, of which one solution of Sb is positive definite (e.g., Hellmer and Olbers 1989; Holland and Jenkins 2001; Losch 2008). In a finite element approach, the prognostic variables vary within elements according to the choice of basis function (here discontinuous linear for velocity and scalars), so the far-field values can be obtained at any location irrespective of the mesh (Kimura et al. 2013). In this study, the far-field properties are taken 1 m away from the ice–ocean interface; that is, the melt rate depends entirely on the pressure, salinity, temperature, and velocity at 1 m from the ice. Melting introduces cool, freshwater into the ocean, referred to here as “meltwater feedback.” When the meltwater feedback is enabled, the temperature and salinity of the ocean are modified by applying inhomogeneous Neumann boundary conditions in accordance with Jenkins et al. (2001). Otherwise, insulated (zero flux) boundary conditions are applied for temperature and salinity along the ice face. The remaining boundaries have insulated boundary conditions for temperature and salinity.
c. Domain configuration
All the configurations are subjected to inflow from a subglacial channel with volume flux Q (Fig. 1). The channel discharges vertically in configuration A and horizontally in configurations B and C. Subglacial discharge water is fresh (Sd = 0) and constrained to be at the local freezing temperature Td = b + cPd. The variable Pd represents the hydrostatic pressure at the depth of the channel: Pd = ρ0gLz, where Lz is the height of the domain. The freezing relationship is linearized to calculate the melt rate in the model, so we employ the same linearization to prescribe the subglacial discharge; if we used a different linearization then unphysical freezing could occur at the inlet. Configuration A allows a simulation of a plume emanated vertically in the absence of obstacles, akin to the environment in which laboratory experiments of plumes are typically performed, and we will use this setup to calibrate our model diffusivity with plume theory. The plume is emanated horizontally in configurations B and C from the bottom of the ocean. An ice face is present directly above the source in configuration C. These configurations are used to identify the effect of an ice face on the evolution of meltwater plumes. The rest of the simulations are performed using configuration C.
Temperature and salinity are restored to the initial conditions of Ta and Sa at the lateral boundary opposite the ice face, with a frequency decreasing linearly between 0.01 s−1 at y = 2400 m and 0 at y = 2000 m (no restoring). At the bottom of the glacier, freshwater with temperature Td enters the ocean, and the ambient temperature and salinity are assumed to be uniform in space with Ta = 2°C and Sa = 33 psu. We impose free-surface boundary conditions on the upper surface. All other boundaries can either take free-slip or no-slip boundary conditions depending on the experiment. Boundary conditions are applied with respect to the discretized function space in which solutions are approximated; that is, the boundary conditions are applied in a “weak” form.
d. Subgrid-scale mixing
The turbulent entrainment processes are subgrid scale and therefore parameterized as a diffusive flux, with the choice of diffusivity dictating the degree of entrainment. We apply Prandtl’s mixing length theory to represent the diffusivity of momentum K. The theory states that the diffusion of momentum is proportional to a typical scale of the fluctuating velocity and the mixing length, defined as the distance traveled by a fluid particle before it loses its momentum anomaly (Taylor 1915; Prandtl 1925). Since our mesh has a variable resolution, we expect the level of turbulence that is subgrid scale to vary spatially; as the mesh gets coarser, the turbulent motions are less well resolved. As a first assumption, we assume that the unresolved turbulent diffusion scales with the ratio between the grid scale and the mixing length. If the mixing length is regarded as a fixed quantity, we may then scale the subgrid-scale mixing by the grid scale. This gives K a spatially varying form:
where u* is the turbulent velocity scale (regarded as fixed) and dr is the length of separation between computational nodes. Our computational domain is represented by tetrahedral elements. Each tetrahedron is composed of four nodes and faces. By changing dr, we can change the distribution of nodes, and the use of unstructured meshes allows total freedom in the placement of nodes. We represent dr as a linear function of y:
where drg and dro represent the spatial resolution along the ice face and open ocean, respectively. The variable Ly represents the width of the domain in the y direction. We allocate finer resolution near the ice face (drg = 5 m) than the open ocean (dro = 100 m) (Fig. 1). We assume that the horizontal and vertical diffusion of momentum, temperature, and salinity are represented by an isotropic diffusivity K.
a. Calibration of subgrid-scale mixing
To calibrate u*, we simulate a meltwater plume in a homogeneous ambient fluid, with the strength of the entrainment of ambient fluid into the plume dictated by our choice of u*. We assess the uncertainty in the entrainment by comparing our simulations having different u* with plume theory, which predicts changes in the plume width with respect to height. When u* = 0, the numerical diffusion is the only source of diffusion in the model, and the plume water remains less than 1.8°C, with individual element boundaries visible in the solution (Fig. 2a). As we increase the physical diffusion, the plume develops a coherent rising trajectory (Figs. 2b,c), the plume water becomes warmer, and the width of the plume increases. When diffusion dominates, the plume trajectory is very different (Fig. 2d).
Classical plume theory assumes that the horizontal profiles of mean vertical velocity and buoyancy are of a similar form at all heights, when a buoyant fluid is released vertically into a quiescent environment (e.g., Morton et al. 1956). Morton et al. (1956) defines the effective radius of the plume R as the distance from the central axis to points at which the amplitude of the vertical velocity is 1/e of its axial value. In the absence of ambient stratification, the buoyancy flux of the plume remains constant, and the slope of R with respect to z is proportional to the entrainment constant αe (Morton et al. 1956):
The entrainment constant αe relates the radial entrainment velocity υe to the mean vertical velocity we in the plume by υe = αewe. We tune u* by searching for an appropriate dR/dz for the case of a “pure” plume in a uniform ambient fluid (Ta = 20°C and Sa = 33 psu) with no ice face (configuration A). The pure plume emanates a buoyancy flux but without momentum or volume fluxes from a virtual source, which is located below the actual source. At the actual source, the rising plume has picked up the momentum and volume fluxes. Hunt and Kaye (2005) categorized plumes into three different regimes by a parameter Γ that characterizes the local balance of momentum, buoyancy, and volume fluxes at the actual source:
where M and F are the momentum and buoyancy fluxes, respectively. The plume is classified to be lazy (Γ > 1), pure (Γ = 1), or forced (Γ < 1). A source of freshwater, constrained to the local freezing point, is released from an upward-looking square channel of 5 m in width (see Fig. 1a). As a result, the inflow velocity of the pure plume is 4 m s−1 (Q = 100 m3 s−1), ensuring that Γ = 1. The theory of Morton et al. (1956) assumes that the plume is emanated from a point source. In contrast, the plume has a finite width in our simulations, so we expect that the comparison to the theory near the source to be less accurate than away from the source.
Laboratory experiments of plumes show that αe varies between 0.07 and 0.16 (e.g., Morton et al. 1956; Linden 2002). Increasing u* results in increasing αe, yielding 0.04, 0.06, 0.08, and 0.3 for u* = 0, u* = 0.05 m s−1, u* = 0.1 m s−1, and u* = 1 m s−1, respectively (Figs. 3a,b). When u* is too high (u* = 1 m s−1), diffusion dominates and erodes the ambient stratification, and the plume is no longer rising in a uniform environment; the plume rises in a stratified environment, and dR/dz is not constant (Fig. 3a). The entrainment constant of 0.08 from u* = 0.1 m s−1 is consistent with laboratory experiments, so we choose u* = 0.1 m s−1 for all subsequent simulations. Increasing u* from 0 to 0.1 m s−1 results in a doubling of the entrainment rate, which suggests that the numerical diffusion is at most 50%.
The volume flux used in our later experiments ranges between 2.5 (Γ ~ 2000) and 200 m3 s−1 (Γ ~ 0.25), and αe is expected to be different in these cases. For example, List and Imberger (1973) argue that αe = 0.056 for jets (Γ = 0), smaller than the pure plume. Kaminski et al. (2005) proposed a model of αe for forced plumes (Γ < 1), and their model suggests an increase in αe with increasing Γ. For plumes, it remains difficult to cite a definitive value for αe as Q varies, as the different studies give notably divergent results, and there is no universal function to describe αe with respect to Γ. We use u* = 0.1 m s−1 throughout our study.
b. Base cases
Most insight into plume dynamics originates from laboratory experiments, where buoyant fluid is released vertically into an open environment (e.g., Priestley and Ball 1955; Morton et al. 1956). In contrast, in our simulations buoyant fluid is discharged horizontally into the ocean from the bottom of the vertical ice front (Figs. 1b,c). Laboratory experiments of buoyant jets in a confined space find that, as the momentum of the jet decreases, buoyancy forces start to dominate and the jet transforms into a plume (e.g., List 1982; Turner 1986). The base case has a square channel of width d = 10 m and height h = 10 m at the base of the ice face and variable inflow velocity of υd. These physical parameters can be combined to form a Froude number, which characterizes the flow at t = 0:
The Froude number represents the ratio of inertial to buoyancy forces. Laboratory experiments and modeling studies show that the trajectory of the jet depends on Fr (e.g., Jirka and Harleman 1979; Kuang and Lee 2001; Sobey et al. 1988; Arakeri et al. 2000; He et al. 2002). A jet with higher Fr has a higher inertial force and therefore a flatter trajectory than to a lower Fr jet for a given buoyancy forcing (Arakeri et al. 2000). In a set of base cases, we study four cases of varying discharge rate, with a range of values of Fr (Table 1).
The subglacial discharge creates a plume attached to the ice face (Fig. 4a). An umbrella of buoyant fluid is formed at the leading edge of the ascending fluid, followed by a steady plume beneath. As the plume ascends, the umbrella becomes wider, and the buoyant fluid within it loses its buoyancy due to the entrainment of the ambient fluid (Fig. 4b). When the plume reaches the surface, the cool, freshwater spreads laterally away from the plume as a gravity current, with a sharp front at its leading edge (Fig. 4c).
The ascending motion of the meltwater plume entrains the surrounding fluid and increases the melt the ice near the trajectory of the plume. As the plume rises, the melt rate along the trajectory of the plume increases. Melting is diagnosed by a vertically averaged melt rate along the center of the ice face above the subglacial discharge and a melt rate averaged over the entire ice face 〈m〉 defined as
where m(x, z, t) and Lx are instantaneous melt rate on the ice face and the width of the domain in the x direction, respectively. All the simulations are run with the same Lx and Lz to keep the and 〈m〉 consistent among different simulations. The trajectory of the plume is confined above the channel and so is the region of high melt rate. The region outside of the plume trajectory does not have high enough melt rate to significantly influence the melt rate averaged over the ice face, so increasing increases the melt rate averaged over the ice face. After at most t = 1000 s, the adjusts to a steady state because the meltwater plume reaches the surface (Fig. 4d). We refer to the melt rate at t = 2000 s as the steady-state melt rate. Because the meltwater plume rises faster with increasing discharge Q, the adjusts to a steady state more quickly, and the faster flow increases the steady-state melt rate. The steady monotonically increases up to Q = 150 m3 s−1.
The spatial pattern of melt rate resembles the shape of the plume, with a band of a high melt rate concentrated above the source (Fig. 5a). The plume is strongly divergent at the surface directly above the source, so it decelerates, and this reduces melting. This is largely unphysical, as it results from the need for a large-scale flow to drive the turbulence in the melting parameterization. When the discharge is low (Q = 30 m3 s−1), the maximum w is located along the ice face (y = 0) and decays away from the ice (Figs. 5b–d). The w along the ice face increases with increasing Q up to Q = 150 m3 s−1. As the discharge increases to Q = 200 m3 s−1, the increase in horizontal momentum pushes the location of maximum w away from the ice, and the w along the ice face decreases (Fig. 5d). Unsurprisingly, the transition to this new behavior occurs at around Fr = 1 (Table 1). The melt rate depends on both the temperature and speed of the plume. The decrease in the vertical velocity along the ice face (Fig. 6a) reduces the melt rate (Fig. 6c), and therefore the temperature inside the plume in Q = 200 m3 s−1 is higher than Q = 150 m3 s−1 (Fig. 6b) due to the reduced meltwater supply.
For illustrative purposes, we also run a case with no discharge (Q = 0). In this case, the temperature and salinity at the plume source are set as Dirichlet boundary conditions, but there is no inflow imposed through the boundary. Despite the lack of an explicit mass flux, the existence of fixed temperature and salinity values causes a diffusive flux, which initiates a significant plume. A diffusive velocity of u* = 0.1 m s−1 combined with a source area of width d = 10 m and height h = 10 m is expected to create a diffusive flux equivalent to an inflow of Q = 10 m3 s−1, and indeed the results are consistent with a flux of this magnitude (Figs. 4–6). It is important to note that this flux exists in all of the model results presented here and indeed in other studies. This diffusive flux can become dominant in cases with a low inflow velocity, that is, a large source area or low discharge (see below).
c. Comparison with meltwater plume theory
Our results are compared to the melt rate predicted by two different formulations of meltwater plume theory (Fig. 6c). The first case shown is for the “line plume” theory detailed in Jenkins (2011), using Q = 200 m3 s−1 and αe = 0.08 and all other parameters as used by Jenkins (2011) or detailed for the simulations above. However, our plumes are initiated from a point source and solved in three dimensions, so, unlike in previous two-dimensional studies, this line-source theory is not strictly applicable. For comparative purposes, we also show the results of an axisymmetric plume equivalent of the model of Jenkins (2011). In the axisymmetric formulation, ice shelf melting is applied to ordinary differential equations for the conservation of heat and salt based on the classical axisymmetric plume model of Morton et al. (1956), rather than the line plume model of Ellison and Turner (1959). The basic difference between the two is that the line plume theory is formulated to conserve momentum, heat, and salt per unit length, whereas these quantities are conserved per unit area in the axisymmetric theory. As a result, in the absence of ambient stratification and for a pure plume, the line plume area increases linearly with height, while the axisymmetric plume area increases quadratically, since in both cases the entrainment assumption dictates that the plume radius increases linearly. The pure axisymmetric plume therefore decelerates in order to conserve its buoyancy flux, while the pure line plume has a velocity constant with height.
In the cases shown, the plume source is very lazy, even with this high value of Q, and so over the first ~150 m the plume is accelerating to adjust to a pure plume balance. After this, the difference in velocity variation between axisymmetric and line plume becomes apparent, with the melt rate decreasing with height in the axisymmetric plume (Fig. 6d), and the line plume melt rate adjusting toward a constant rate as it rises.
The melt rate from the three-dimensional simulations falls between the axisymmetric and line plume results. The shape of the melt rate profiles largely has the characteristics of a line plume, where the melt rate gradually increases with height, despite the use of a point source. The simulated meltwater plume is not a full axisymmetric plume in nature because the wall boundary conditions disrupt the symmetry. This decrease in the entrainment by the wall keeps the pure vertical velocity profile relatively constant in height, and therefore the continual adjustment from the lazy source means that the upper part of the ice melts more than the bottom part (Fig. 6c).
d. Effect of an ice face on a plume
Oceanographic control of sediment transport from subglacial discharge has been the object of many studies (e.g., Syvitski 1989; Powell 1990; Mugford and Dowdeswell 2011). The concentration of sediments suspended in the glacial meltwater is not high enough to overcome its buoyancy, and the sediment-rich water rises as a meltwater plume (Mulder and Syvitski 1995). The trajectory of the meltwater plume dictates the transport of sediments and plays an important role in shaping the sea floor near the glaciers (Syvitski 1989; Powell 1990). Based on plume theory, a number of previous studies predicted that the horizontal momentum of the discharge deflects the plume trajectory away from the source and that the plume surfaces some distance away from the glacier under sufficiently large discharge (e.g., Syvitski 1989; Powell 1990; Mugford and Dowdeswell 2011). The plume theory employed in these studies relies on adding conservation of horizontal momentum to the set of plume equations introduced by Morton et al. (1956) and has been verified in laboratory experiments (e.g., Lane-Serff et al. 1993; Lane-Serff and Moran 2005). Importantly, however, this theory does not consider the influence of a wall above the mouth of the discharge.
The effect of a wall on a jet was first discovered by Young (1800) and later applied to aircraft development by Henri Coandă, known as the Coandă effect (Wille and Fernholz 1965). The Coandă effect is the tendency of a jet to be attracted toward a nearby surface. The effects of vertical walls on plume trajectory have been investigated in laboratory experiments (e.g., Pera and Gebhart 1975; Zukoski et al. 1981). Pera and Gebhart (1975) showed that a plume has a tendency to deflect toward nearby vertical walls. When a semicircular buoyant source was placed next to a vertical wall, the plume attached to the wall (Zukoski et al. 1981). The presence of a wall prevents the entrainment of the surrounding fluid toward the jet on one side, and the jet is attracted toward the wall as a result (Wille and Fernholz 1965).
The manifestation of the Coandă effect is present in our simulations (Fig. 7), suggesting that the Coandă effect influences the trajectory of the meltwater plume. In the presence of a wall, the meltwater plume is attached to the ice face as it rises, even in the relatively high discharge of 200 m3 s−1 (Fig. 7a). In contrast, the horizontal momentum in the discharge deflects the plume trajectory away from the ice without a wall, and the trajectory is closer to the prediction by Lane-Serff et al. (1993) (Fig. 7b). Nevertheless, in neither case does the plume trajectory follow the trajectory predicted by the simple theory of Lane-Serff et al. (1993). Two possible reasons for this are the local buoyancy forcing by the meltwater feedback or the diffusion in the model. The meltwater feedback continuously supplies lighter water along the ice face, which attracts the plume to the wall. The presence of diffusion in the model causes increased mixing and loss of horizontal momentum relative to the plume theory, which has no explicit mixing other than entrainment. However, an additional case with no meltwater feedback or explicit diffusion still follows a path very close to the wall (Fig. 7c). Therefore, we conclude that the Coandă effect is the primary cause of the plume path deviation from the simple theory.
The horizontal momentum and buoyancy forces at the discharge are almost in balance in these simulations, Fr = 1.3 (see Table 1). The two-dimensional simulations of Salcedo-Castro et al. (2011) used Fr ranging between 0 and 3 and found that the plume trajectory was attached to the wall, as in our simulations. However, none of the previous studies have considered simulations with and without a wall, and therefore it has not been clear if the attachment of the plume is due to the presence of a wall or the underlying path preferred by the plume. We conclude that subglacial plumes do not detach from a vertical ice face for realistic levels of discharge, which is critically important to the overall melt rate of the glaciers. However, the shape of the ice face may steer the plume to detach from the wall or attach more firmly. If the ice is slightly undercut at the bottom, the meltwater plume may not have enough horizontal momentum left at the exit of the cavity, and the meltwater plume would rise vertically along the ice face. On the other hand, if the ice is slightly overcut, there may be a critical angle below which the Coandă effect can take place.
The effect of the wall also depends upon the velocity boundary conditions applied at the wall, which are critically important because it is the velocity close to the wall that determines the melt rate. We now consider the difference between boundary conditions of free slip, no slip (applied weakly), and a quadratic drag with a coefficient of 0.0025 (Holland and Feltham 2006). To isolate the effect of the boundary condition on the dynamics, we disable the meltwater feedback in this study. Our default boundary condition in this study is one of free slip (Fig. 8a), because in models where the boundary layer turbulence is not fully resolved it is not justified to apply no-slip conditions. In all cases, the profile of the vertical velocity spreads away from the ice face (Figs. 8c–e). In the case of an ice face with a no-slip boundary condition, the maximum vertical velocity is separated from the ice face at all depths (Fig. 8b). The reduction of the vertical velocity near the ice face leads to a decrease in the melt rate. In the case of a drag boundary condition, the effect is very similar to the free-slip case (Figs. 8c–e).
e. Different channel sizes
We now seek to understand the effect of different source channel geometries. To our knowledge, there are no direct underwater observations of subglacial channel size or volume flux, so we rely on speculated channel size and discharge velocity derived from observations and theory. Motyka et al. (2003) estimated the subglacial discharge to vary from low levels in winter to around 435 m3 s−1 in summer at LeConte Glacier in Alaska. Xu et al. (2012) estimated the discharge of freshwater from each channel in Store Glacier, Greenland, to average 100 m3 s−1 during the May–September melting season, with peak values approaching 250 m3 s−1 in midsummer, using surface runoff from the Regional Atmospheric Climate Model (RACMO) (Ettema et al. 2009) combined with their field observations. Andersen et al. (2010) estimated an average surface runoff of ~174 m3 s−1 during the summers of 2007 and 2008 at Helheim Glacier, east Greenland.
Much of the subglacial water is transported by channels. The channelized system tends to draw in surrounding water and form small numbers of large channels, aligned along the flow of the glacier. The size of a channel in steady state is thought to be determined by balance between the inward deformation of ice and melting of the ice due to turbulent dissipation of flow inside the channel (Röthlisberger 1972; Weertman 1972; Nye 1976). Xu et al. (2012) employed the theoretical framework of Röthlisberger (1972) and estimated the channel area for a given discharge for their numerical experiment. They came up with the area of channels and volume flux to be in the order of 50 m2 and 100 m3 s−1. For low rates of discharge, the subglacial hydrological network can be described as a distributed system, where the water beneath the glacier is spread over a large part of the bed rather than being localized into channels. The distributed system consists of a network of cavities, which are formed by ice flow over bumps in the bed (Kamb 1987). The cavities are interconnected, allowing the flow of water between the cavities and eventually to the terminus of the glacier.
Previous simulations used a two-dimensional domain (e.g., Salcedo-Castro et al. 2011; Xu et al. 2012; Sciascia et al. 2013), so the channel sizes were altered by changing their heights. Here, we explore the effects of changing the height and width of the channel. Four different channel sizes are considered in this study (Table 1). The widest channel, d × h = 100 m × 5 m, is our best approximation of the distributed system, while the rest of the channels are designed to represent the channelized system.
For each channel we consider a range of Q up to 200 m3 s−1. The minimum value of Q depends on the channel size and imposed diffusion. Since we assume the turbulent velocity scale of u* = 0.1 m s−1, the effect of any inflow velocity below 0.1 m s−1 is not distinguishable from the effect of diffusion. Our results are therefore valid only when the inflow velocity is larger than 0.1 m s−1, so the minimum values of Q are 2.5, 10, 20, and 50 m3 s−1 for d × h = 5 m × 5 m, d × h = 10 m × 10 m, d × h = 20 m × 10 m, and d × h = 100 m × 5 m, respectively (Table 1). The ranges of inflow velocity and Q are the same order of magnitude as previous studies (e.g., Xu et al. 2012; Salcedo-Castro et al. 2011; Sciascia et al. 2013; Xu et al. 2013).
The steady-state and 〈m〉 increase with increasing Q, but the proportionality is dependent on channel shape (Figs. 9a,b). Because the melt rate outside of the plume trajectory is lower, 〈m〉 is lower than . For a given Q, the largest and 〈m〉 values are found with d × h = 100 m × 5 m. The meltwater plume theory of Jenkins (2011) predicts that the melt rate is proportional to . As discussed previously, the from smaller channels flattens out after reaching a critical Q because the high velocity of the inflow forces the fast plume flow away from the ice face. The critical Q increases with increasing channel size, Q = 125 m3 s−1 and Q = 150 m3 s−1 for d × h = 5 m × 5 m and d × h = 10 m × 10 m, because the Froude number depends upon the source velocity, not discharge. The experiment d × h = 20 m × 10 m does not reach the critical Q in the range of Q used in our experiments.
The widest channel (d × h = 100 m × 5 m) has the largest melt rate for given Q (Fig. 9). A channel of size d × h = 100 m × 5 m discharging Q = 100 m3 s−1 is the equivalent of 20 adjacent channels, with each channel having d × h = 5 m × 5 m and discharging Q = 5 m3 s−1 (Fig. 10). However, the melt rate is only 2.5 times larger than a single channel of this small source size and discharge. When the channel is small enough, the temperature contours spread out horizontally with increasing z, akin to the axisymmetric plume (Figs. 10a,c). In the case of d × h = 100 m × 5 m, the width of the temperature contours decreases between z = −500 m and z = −450 m and the contours stay nearly vertical afterward (Fig. 10b).
At the bottom, the temperature profiles have a single minimum at the center, corresponding to the temperature of discharged water, and the temperature minimum spans the width of the channel (Fig. 10g). As the plume ascends, the temperature profile retains the same Gaussian shape with increasing plume width in the d × h = 5 m × 5 m case (Fig. 10f). In contrast, the profile has two local minima in the d × h = 100 m × 5 m case, as if there are two plumes ascending next to each other (Fig. 10f). The source has become wider than is sustainable for a single coherent plume of this given flux, so the plume geometry is unstable and partially separates into two cores. Higher up, these two plumes coalesce back into a single feature (Fig. 10e). Near the bottom, temperatures at the center of the plume are the same in d × h = 100 m × 5 m and d × h = 5 m × 5 m cases. As the plume rises, the core of the plume in d × h = 100 m × 5 m is colder than d × h = 5 m × 5 m, and two local minima in the d × h = 100 m × 5 m merge into a single minimum. The coalesced plume is colder and more buoyant than the single plume (Fig. 10d), so it ascends faster and fluxes more heat to the ice base. This results in an increase in the melt rate along the center of the plume for d × h = 100 m × 5 m relative to d × h = 5 m × 5 m cases. In the next section, we examine such plume coalescing in more detail.
f. Coalescing plumes from multiple channels
Observations of meltwater plumes near the terminus of glaciers rely on visual inspection, satellite imagery, or water samples (e.g., Dowdeswell and Cromack 1991; Horne 1985; Greisman 1979; Syvitski 1989; Chu et al. 2009). The surface signature of the plumes is characterized by a turbulent and turbid mass of fluid near the terminus of glaciers. All the modeling studies assume that the meltwater plume is emanated from a single channel (e.g., Xu et al. 2012; Salcedo-Castro et al. 2011; Sciascia et al. 2013; Salcedo-Castro et al. 2013; Jenkins 2011; Mugford and Dowdeswell 2011), but as the subglacial hydrological system develops over the melt season, it is entirely possible that multiple channels would discharge into a fjord. There have been no attempts to study the effects of interacting meltwater plumes and their effect on the melt rate. Two plumes in close proximity are known to coalesce and form a single plume. Pera and Gebhart (1975) studied the interactions of two line plumes, merging to form a single plume. The merging process is a result of the restriction of entrainment into each plume by the presence of the other. If the distance between the plumes is close enough, the two plumes merge and form a single plume (Pera and Gebhart 1975). This type of merging can also happen for pairs of axisymmetric plumes (Kaye and Linden 2004). In the case of two plumes having unequal strength, the stronger plume remains unaffected and attracts the weaker plume (Pera and Gebhart 1975; Kaye and Linden 2004).
We test the effect of coalescing plumes by simulating a subglacial discharge from two different channels separated by a distance dx. The edge of each channel is located dx/2 away from the center. We are considering the multiple-channel stage of subglacial development, when it is very unlikely that each channel supports a large fraction of the total catchment drainage. We therefore choose a low discharge flux estimate of 40 m3 s−1 (e.g., Xu et al. 2012), with each channel discharging 20 m3 s−1 from a channel size of d × h = 5 m × 5 m. The dx has a significant effect on the meltwater plumes and melt rate (Fig. 11).
In the case of dx = 0, a high melt rate is concentrated above a channel of d × h = 10 m × 5 m, discharging Q = 40 m3 s−1, located at the center (Fig. 11a). The melt rate is theoretically proportional to in the presence of a meltwater plume (Jenkins 2011), so a multiple-channel system should be more efficient in melting ice than a single channel for a given total discharge, and this is borne out by the results for dx > 0 (Figs. 11–13). However, the system is not straightforward in the presence of plumes separated by only a short distance, since these plumes may coalesce into a single plume. In the case of dx = 10 m, a confined band of high melting is concentrated above the center, despite the two channels being from separated sources (Fig. 11b). Near the bottom, both dx = 0 and dx = 10 m have similar melt rate, but, as the plume ascends, the melt rate from dx = 10 m becomes higher than dx = 0 due to a merging of the plumes (Fig. 12a). The signals of the plumes in the vertical velocity and temperature are evident near the channels (Figs. 12d,g). As the plume rises, the two extrema of the vertical velocity and temperature merge into one (Figs. 12c,f), and the merged plume is colder and ascends faster than dx = 0 (Figs. 12b,e). The total transport in the dx = 10 m case equals that of the dx = 0 case, but the merged plume is more concentrated spatially, since adjacent plumes suppress the lateral entrainment into each plume. As a result, the pressure between the two plumes becomes relatively low compared to the outer edges, and the two plumes are attracted toward the center, forming a single plume. The merging event makes it harder to define the core of the plume, so we use the maximum melt rate to quantify the effect of coalescing, and we use the area-averaged melt rate to diagnose the overall effect. The maximum melt rate is the largest time-averaged melt rate between t = 1000 and 2000 s. The area-averaged melt rate is the space- and time-averaged melt rate over the ice face.
Kaye and Linden (2004) showed that Γ > 1 at the point in which two plumes with equal volume fluxes merge, that is, the merged plume is lazy. The lazy plume is known to contract during the adjustment phase from lazy to pure, and the plume suffers appreciable acceleration (Morton and Middleton 1973). This acceleration enhances the heat flux toward the ice base, and therefore the local melt rate of the merged plume is higher than a single plume (Figs. 11b, 13a). However, the contraction of the plume makes the merged plume confined to the center and less effective in melting the entire ice face, so the area-averaged melt rate is lower than a single plume (Fig. 13b). This argument lends support to the fact that the mean centerline melt rate from d × h = 100 m × 5 m, which has a merging lazy plume, is higher than d × h = 5 m × 5 m (Fig. 9).
The coalescing becomes less obvious when the plumes are separated further (Fig. 11c). A separation of dx = 50 m is not far enough for individual plumes to ascend on their own (Figs. 12b–d). Bands of ascending fluid appear above each channel with a weak coalescing region in the center. Because coalescing becomes weaker as the separation increases, the maximum melt rate decreases (Fig. 13a). However, two separate plumes can occupy a larger area than a coalesced plume, increasing the melt rate averaged over an ice face (Fig. 13b). At a sufficient distance of separation (e.g., dx = 200 m), the plumes do not coalesce, and the individual plumes of equal strength remain separated (Fig. 11d). Lateral circulation is generated once the plumes surface. This lateral circulation results in shallow overturning circulation between two plumes, so the melt rate of the plumes near the surface is slightly affected.
4. Discussion and conclusions
The dynamics of subglacial meltwater emerging beneath a glacier face in a fjord plays an important role in governing the ocean melting of the glacier and the distribution of sediment deposits in the fjord. In this study, we investigate various features of this scenario using a nonhydrostatic ocean model with a mesh that is unstructured in all three dimensions. Although our simulations are intended to shed light on the role of the meltwater plume in melting a vertical ice face, the following caveats should be noted:
Subgrid-scale mixing, which influences the entrainment of surrounding fluid toward the ice, is calibrated by comparison to a pure plume (Γ = 1). The Γ in our numerical experiments ranges from 0.25 to 2000, depending on the volume flux. To the best of our knowledge, there is no universal function to describe the entrainment constant with respect to Γ.
We have assumed that the ambient temperature and salinity are uniform in space. In the presence of strong stratification, the meltwater plume can intrude into the interior of the ocean and never reach the surface, as in the simulations of Sciascia et al. (2013).
The basic results are similar to previous studies; increasing the discharge rate increases the melting of the glacier face by accelerating the plume and thus increasing the turbulent mixing of heat across the ice–ocean boundary layer. Though the subglacial meltwater is injected horizontally into the ocean, we find that the resulting plume always adheres to the ice face, in contradiction to the simple plume theory employed in several previous studies. The character and melt rate of the plume is critically dependent upon the momentum boundary conditions applied at the wall.
Much of our understanding of the behavior of meltwater plumes comes from one-dimensional models (e.g., MacAyeal 1985; Jenkins 1991, 2011) based on the theory of line plumes (Ellison and Turner 1959). The classical result of a pure line plume is that the plume velocity, and hence the predicted melt rate in the meltwater case, is constant with height. However, in reality subglacial discharges are typically from a channel, akin to a point source, and therefore such plumes might be expected to have an axisymmetric character, in which the classical pure result would be for the velocity, and hence the melt rate, to decrease with height. Our simulations show that, in fact, the predicted melt rates increase with height, and the melt rate falls between the theoretical predictions of line and single-source plumes. We attribute this to the fact that even for high discharge rates, the plumes are found to be “lazy,” in the sense that the momentum of the source is far less than its buoyancy flux would generate in a pure plume momentum balance (Hunt and Kaye 2005). The plumes are continually accelerating toward such a balance over nearly the entire height of the ice face, and this causes an increase in melting with height. This spatial pattern of melting will “overcut” the ice, making the bottom part of the glacier more buoyant than the top. This could play a significant role in the calving of these glacier faces, which is commonly observed to occur through bergs toppling inwards toward the glacier at the surface.
We also test the effect of changing the source discharge and geometry. Different source geometries can give a variety of different melt rates for any given discharge, with higher melting found for a thinner, wider source. Therefore, there is no single relation between discharge and melting, as proposed in many previous studies. For a given channel geometry, however, the melt rate is approximately proportional to discharge to the one-third power, as proposed previously. When two subglacial channels are in close proximity to each other, meltwater plumes can coalesce and form a single plume, which ascends faster and melts more ice within the core of the meltwater plume than the two plumes, but is inefficient in melting a wide area of the glacier. All of these results imply that glacier melting will vary in a highly complex manner over a seasonal cycle, as the subglacial meltwater source changes in discharge, and the geometry of the source concurrently changes as the subglacial system changes from a distributed system to an evolving channelized flow.
Comments from two anonymous reviewers improved the manuscript. The work was supported by the U.K. Natural Environment Research Council under Grants NE/H02333X/1 and NE/J005746/1.