Diffuse debris entrainment in glacier, lab and model environments

Abstract

Small quantities of liquid water lining triple junctions in polycrystalline glacier ice form connected vein networks that enable material exchange with underlying basal environments. Diffuse debris concentrations commonly observed in ice marginal regions might be attributed to this mechanism. Following recent cryogenic ring-shear experiments, we observed emplacement along grain boundaries of loess particles several tens of microns in size. Here, we describe an idealized model of vein liquid flow to elucidate conditions favoring such particle transport. Gradients in liquid potential drive flow toward colder temperatures and lower solute concentrations, while deviations of the ice stress state from hydrostatic balance produce additional suction toward anomalously low ice pressures. Our model predicts particle entrainment following both modest warming along the basal interface resulting from anticipated natural changes in effective stress, and the interior relaxation of temperature and solute concentration imposed by our experimental protocols. Comparisons with experimental observations are encouraging, but suggest that liquid flow rates are somewhat higher and/or more effective at dragging larger particles than predicted by our idealized model with nominal parameter choices. Diffuse debris entrainment extending several meters above the glacier bed likely requires a more sophisticated treatment that incorporates effects of ice deformation or other processes.

1. Introduction

Entrained debris muddy glacier soles, producing distinctive basal layers that influence ice and landscape dynamics (e.g. Knight, 1997; Hubbard and others, 2009; Waller and others, 2009). Several different entrainment mechanisms have been identified, with supporting observational evidence derived from exposed margins (e.g. Sugden and others, 1987; Hubbard and Sharp, 1995; Lawson and others, 1998, see Fig. 1a), remote segments of glacier interiors (e.g. Jansson and others, 1996; Carsey and others, 2002), the surfaces of capsized icebergs (e.g. Anderson and others, 1980; Pierce and others, 2021), and controlled laboratory experiments (e.g. Knight and Knight, 1994, 1999; Hansen and Zoet, 2022). Direct access to the basal environment is difficult, preventing in situ observations of basal layer formation. Model treatments designed to elucidate the physical and chemical fingerprints of entrainment processes hold promise for providing rare and valuable insights into the subglacial conditions that produce these enigmatic features (e.g. Lliboutry, 1993; Alley and others, 1997; Cook and others, 2006; Rempel and others, 2022).

Figure 1.

(a) Photograph showing diffuse debris in basal ice near Kangerlussuaq, Greenland (reproduced with permission from Knight and Knight, 1994). (b) Photograph of four veins meeting at a node to form a connected network through polycrystalline ice (provided by Heidy Mader; see Mader, 1992a). Note that the center of curvature for each of the vein walls is inside one of the four ice grains depicted so that vein cross-sections more closely resemble triangles than circles, with optical interference producing brightening along the middle (third) edges of the two sub-horizontal veins and darkening on the middle edges of the two sub-vertical veins shown here.

[Figure omitted. See PDF]

During recent cryogenic ring-shear experiments, we slid temperate ice over several unconsolidated substrates and observed different patterns of diffuse debris entrainment. Following ice slip over Horicon till (median particle diameter ~350 μm), a sharp boundary between debris-free ice and a fringe of ice-infiltrated till was maintained, whereas slip over Iowa loess (median diameter ~48 μm) led to diffuse entrainment above the fringe characterized by the deposition of fine particles (median diameter ~20 μm) along the veins that line three-grain junctions in polycrystalline ice, and during the early stages of slip over a bed of glass beads (median diameter ~160 μm) we observed particles moving up into the ice, but infer that these particles had escaped before the end of the ~20 day experimental run. To investigate, here we develop an idealized continuum model of vein liquid flow (see Fig. 1b) that builds upon empirically backed formulations of vein evolution by Mader (1992b, 1992) and Nye (1991), and an idealized treatment of water exchange due to vein flow above subglacial lakes (Rempel, 2005), to predict how flow characteristics respond to changing conditions (e.g. along the basal boundary).

Our analysis suggests a prominent role for the following physical elements in governing the diffuse entrainment we observed. Gradients in liquid potential favor flow toward colder temperatures, lower solute concentrations and reduced ice pressures, so transient temperature, solute, and pressure fields produce changes in flow rate and the maximum sizes of mineral particles that can be entrained and transported. In the idealized case where the ice pressure gradient is hydrostatic, changes in the liquid potential are proportional to a linear combination of changes in the temperature and solute concentration, which are coupled through their influence on the liquid volume fraction. This coupling causes the immediate response to an increase in the local temperature, for example, to be mirrored by an opposing decrease in the solute concentration that leaves the liquid potential nearly unchanged, limiting the magnitudes of the potential gradients that drive liquid flow over the relatively short time scales that are associated with thermal diffusion. Over the longer time scales associated with compositional diffusion, differing rates of heat and solute transport lead inevitably to changes in the local liquid fraction that are accompanied by steepening of potential gradients, driving more rapid liquid flow that can transport larger particles. Moreover, following an abrupt perturbation that initially causes vein liquid transport upward into the ice, we find that the flow direction can reverse so that subsequent vein liquid transport is directed downward toward the bed. Our model predicts that liquid flow rate magnitudes tend to diminish with distance above the glacier bed so that particle transport is confined to the lowermost basal ice, in qualitative agreement with experimental observations and field evidence.

We present our continuum model of vein liquid flow next and develop scaled governing equations that are applicable when the ice pressure distribution can be approximated as hydrostatic. With this framework in place, we illustrate aspects of the system behavior by examining how vein liquid flow can transport particles in response to changes in basal temperature of the magnitude expected during sliding over a heterogeneous bed that contains both macroscopic drainage elements and areas of close contact between basal ice and till or bedrock. This leads into an analysis of a scenario in which vein flow is forced by changes in conditions above the basal interface, as may have been responsible for the diffuse debris entrainment we observed during the ring-shear experiment over loess mentioned above. Measured particle size distributions and concentrations following our loess experiment reveal that particle transport is somewhat more effective than predicted using nominal parameter estimates in our idealized model, suggesting that some combination of higher vein permeabilities and/or the increases to particle drag that are associated with flow in confined spaces might be responsible for causing these discrepancies. To explore an additional degree of complexity that is expected in many basal environments, we next outline how our treatment can be modified slightly to address cases where the ice pressure distribution deviates from hydrostatic. We provide a brief discussion of model implications and limitations before offering concluding remarks.

2. Continuum model of vein flow

Consider a block of polycrystalline ice that is subjected to a prescribed set of changes in temperature, liquid pressure, solute concentration and ice normal stress along its exterior boundaries. When liquid flow is directed toward the interior at sufficient rates, sediment particles can be lofted through the vein network and incorporated within basal ice layers. Averaged over a volume spanning many ice grains, Darcy's law describes the volume flux of liquid through a polycrystalline cross section as (e.g. Nye and Frank, 1973)1\[{\bf q} = -{k\over \mu}\left(\nabla P_l - \rho_l {\bf g}\right),\; \]where the permeability k is assumed here to be isotropic, μ is liquid viscosity, P_l is vein liquid pressure, ρ_l is liquid density, and g is the acceleration of gravity (i.e. for vertical unit vector \(\hat {{\bf z}}\), we can write \({\bf g} = - g\hat {{\bf z}}\), where g = |g|). Eqn. (1) indicates that the potential gradient that determines the direction of flow is controlled by gravity and the gradient in liquid pressure. Below, by enforcing local equilibrium between vein liquid and its bounding ice walls, we express the liquid pressure in terms of temperature T, solute concentration c and ice pressure P_i. With the liquid fraction (i.e. vein volume fraction) denoted by ϕ, the conservation conditions give rise to the following evolution equations for the temperature and the solute concentration in the vein liquid:2\[{\partial T\over \partial t} = {K\over \rho_iC_i} \nabla^2\, T - {L\over C_i} {\partial \phi\over \partial t} \; ,\; \]3\[{\partial c\over \partial t} = {D\over \phi}\nabla\cdot\left(\phi\nabla c\right)-{{\bf q}\over \phi}\cdot \nabla c -{c\over \phi}{\partial \phi\over \partial t}.\]Here, ρ_i is ice density, K is thermal conductivity, C_i is the specific heat capacity of ice, L is latent heat, and D is the effective compositional diffusivity. The two terms on the right side of Eqn. (2) represent thermal diffusion through the ice and the latent heat associated with changes in phase; on the right side of Eqn. (3), the three terms represent solute diffusion through the veins, advective transport of solute with the flowing vein liquid, and the concentrating (or diluting) effect of decreases (or increases) in liquid volume. In this formulation: i) ϕ is assumed to be small enough that the thermal conductivity and heat capacity of the liquid play no role and heat transport by advection can be neglected, ii) heat exchange between crystal interiors and boundaries is assumed rapid enough for the thermal state to be characterized by the local equilibrium vein temperature, iii) the total quantity of dissolved solute within the vein network can evolve only through exchange across the ice base, consistent with expectations when solutes are rejected entirely from the ice lattice so that they are segregated to the liquid veins and iv) we neglect both the contributions to vein flow that result from the small difference between ice and liquid densities as the liquid fraction evolves, and the ice deformation that is needed to compensate for gradients in liquid transport rate.

The ice–liquid surface energy and curvature of vein walls produce differences between the ice and liquid pressures, P_i and P_l. The equilibrium relation along these surfaces can be expressed in terms of the deviation of the local temperature from a reference bulk melting temperature T_m (defined for phase equilibrium at reference pressure P_i = P_l = P_m, e.g. T_m ≈ 273 K for P_m ≈ 10⁵ Pa), the ice pressure, the liquid pressure and the solute concentration c as (c.f. Nye, 1991; Rempel, 2005; Dash and others, 2006)4\[T_m-T = {T_m\over \rho_i L}P_i - {T_m\over \rho_l L} P_l - {T_m\left(\rho_l-\rho_i\right)\over \rho_i \rho_l L} P_m + \Gamma c,\; \]where solute concentrations in the example calculations below will be assumed low enough that the liquidus slope Γ can be treated in the limit of dilute solution theory as proportional to the gas constant R_g, with \(\Gamma \approx R_gT_m^2/( \rho _l L)\). Hence, the potential gradient that drives liquid flow and appears in parentheses on the right side of Eqn. (1) becomes5\[\nabla P_l -\rho_l{\bf g} = {\rho_l L\over T_m}\nabla T + {\rho_l L \Gamma\over T_m}\nabla c + {\rho_l\over \rho_i}\nabla P_i -\rho_l{\bf g},\; \]notably describing the tendency for liquid in a polycrystal to be drawn from warm to cold temperatures and salty to fresh conditions. Following Nye and Frank (1973), we assume that the permeability has a quadratic dependence on liquid fraction and can be expressed as k = k₀(ϕ/ϕ₀)², where k₀ is the permeability at characteristic liquid fraction ϕ₀. Substituting this and Eqn. (5) into Eqn. (1), while approximating the ice pressure as hydrostatic for now (i.e. \(\nabla P_i\approx \rho _i{\bf g}\) – we briefly discuss how nonhydrostatic pressure gradients affect vein transport later) gives the Darcy transport rate as6\[{\bf q} = -{\rho_l Lk_0 \over \mu T_m\phi_0^2}\phi^2\left(\nabla T + \Gamma \nabla c\right),\; \]where the ratio of terms that forms the prefactor on the right defines a constant. Because changes in temperature and solute concentration occur at different rates, as reflected by the differences between the right sides of Eqns. (2) and (3), over time, gradients in temperature and concentration are expected to vary in strength and relative importance for determining the direction of the liquid flow described by Eqn. (6), while the quadratic dependence on ϕ indicates that changes in flow magnitude are particularly sensitive as well to changes in the liquid fraction.

The Laplace pressure difference between the ice and liquid phases that is associated with surface energy γ is P_i − P_l = γ/R_v, where R_v is the radius of curvature of the vein walls (see Fig. 1b). The liquid fraction for an equiaxed polycrystal with grain size d bounded along triple-junctions (i.e. where 3 grains meet) by veins characterized by R_v can be obtained from (e.g. Nye, 1991)7\[\phi\approx 3\alpha \left({R_{v}\over d}\right)^2,\; \]where the coefficient \(\alpha = \sqrt {3}\sin ^2( \pi /6-\psi /2) -3( \pi /6-\psi /2) + ( 3/2) \sin ( \pi /3-\psi )\) varies over a fairly tight range for the span of dihedral angles ψ at the contacts between ice grains that are typically observed (Mader, 1992a). Focusing our attention on conditions under which the grain size d can be treated as constant, with α effectively constant as well, Eqn. (7) implies that changes in liquid fraction must be attributed solely to changes in vein radius and the characteristic liquid fraction ϕ₀ can be expressed in terms of a characteristic vein radius R_v0 so that \(\phi _0 = 3\alpha R_{v0}^2/d^2\). We orient the z coordinate vertically and assign the reference pressure P_m as the ice pressure at location z = 0 so that when the ice pressure distribution can be treated as hydrostatic: P_i = P_m − ρ_igz. Upon substituting this, the Laplace pressure difference, and Eqn. (7) into Eqn. (4), the liquid fraction becomes8\[\phi = {\phi_0\gamma^2\over R_{v0}^2}\left[{\rho_l L\over T_m}\left(T_m-T- \Gamma c\right) + \left(\rho_l-\rho_i\right)g z\right]^{-2},\; \]where ϕ₀ and R_v0 are the reference values of liquid fraction and vein radius introduced above. Not surprisingly, ϕ increases when the temperature warms or the solute concentration is enhanced, while changes in gravitational potential energy cause modest reductions in ϕ with elevation z.

2.1 Scaled model equations

The evolution of T and c described by Eqn. (2) and (3) together with the Darcy transport rate from Eqn. (6) and the description for liquid fraction in Eqn. (8) can be solved together once appropriate boundary and initial conditions are chosen. The results we describe below emerge as particular examples of such model solutions. To aid in subsequent interpretations, however, it is useful to first introduce a set of characteristic scales and group parameters so that the model can be cast in dimensionless form and the relative importance of the various terms and the processes they represent can be assessed. Accordingly, in addition to the characteristic liquid fraction ϕ₀ introduced above, we also identify representative values for temperature T₀, concentration c₀, distance Δz, and time Δt, facilitating the following definitions, in which overlying tildes indicate scaled variables and operators:9\[\tilde \phi = {\phi\over \phi_0},\; \quad\tilde T = {T-T_0\over T_m-T_0},\; \quad \tilde c = {c\over c_0},\; \]10\[\tilde \nabla = \Delta z \nabla,\; \quad \tilde t = {t\over \Delta t} = {tK\over \rho_i C_i\Delta z^2},\; \]11\[{\cal C} = {\Gamma c_0\over T_m-T_0},\; \quad{\cal S} = {L\over C_i\left(T_m-T_0\right)},\; \]12\[L_e = {K\over \rho_i C_i D},\; \quad\beta = {\rho_l L k_0\left(T_m-T_0\right)\rho_i C_i\over \mu T_m K}.\]The four key dimensionless parameters in Eqns. (11) and (12) represent: \({\cal C}\) – the fractional importance of the undercooling caused by solutal effects (the remaining \(1-{\cal C}\) is attributed to surface energy along curved vein interfaces); \({\cal S}\) – the heat required to change phase relative to that required to change temperature by the characteristic amount T_m − T₀; L_e – the rate of thermal diffusion relative to compositional diffusion; and β – the ratio of time scales for diffusion of latent heat and liquid flow.

Adopting the scalings and definitions in Eqns. (9)–(12), the conservation equations are written as13\[{\partial \tilde T\over \partial \tilde t} = \tilde \nabla^2 \tilde T - {\cal S}\phi_0 {\partial \tilde \phi\over \partial \tilde t},\; \]14\[{\partial \tilde c\over \partial \tilde t} = { 1\over L_e\tilde \phi}\tilde \nabla\cdot\left(\tilde \phi\tilde \nabla \tilde c\right)-{\tilde{{\bf q}} \over \phi_0\tilde \phi} \cdot \tilde \nabla \tilde c -{\tilde c\over \tilde\phi} {\partial\tilde \phi\over \partial \tilde t} ,\; \]where the dimensionless Darcy flow rate is defined as15\[\tilde{{\bf q}} = {\Delta t\over \Delta z}{\bf q} = -\beta \tilde \phi^2\left(\tilde \nabla \tilde T + {\cal C} \tilde \nabla \tilde c\right)\; .\]Choosing the temperature scale to correspond with the undercooling at z = 0 when c = c₀, the scaled liquid fraction from Eqn. (8) becomes16\[\tilde \phi = \left({1-\tilde T -{\cal C} \tilde c\over 1-{\cal C}} + {R_{v0}\Delta z\over \lambda_c^2}\tilde z\right)^{-2}\; ,\; \]where the capillary length is defined following the standard convention as \(\lambda _c = \sqrt {\gamma /[ ( \rho _l-\rho _i) g] }\). Equation (15) and the time derivative of scaled liquid fraction from Eqn. (16) can now be substituted into Eqns. (13) and (14) to write the evolution of temperature and concentration (following algebraic manipulations) as17\[\eqalign{& \displaystyle{{\partial \tilde{T}} \over {\partial \tilde{t}}} = \bigg \{ \bigg (C\tilde{c} + \displaystyle {{1 - C} \over {2{\tilde{\phi }}^{1/2}}} \bigg ) {\tilde{\nabla }}^2\tilde{T}-CS\phi _0\tilde{\phi } \bigg [ \displaystyle{1 \over {L_e\tilde{\phi }}}\tilde{\nabla } \cdot (\tilde{\phi }\tilde{\nabla }\tilde{c}) \cr & \left. {\quad \quad + \displaystyle{{\beta \tilde{\phi }} \over {{\tilde{\phi }}_0}}(\tilde{\nabla }\tilde{T} + C\tilde{\nabla }\tilde{c})\cdot \tilde{\nabla }\tilde{c}} \right] \bigg \}\times \bigg (S\phi _0\tilde{\phi } + C\tilde{c} + \displaystyle{{1-C} \over {2{\tilde{\phi }}^{1/2}}} \bigg ) ^{-1},} \]18\[\eqalign{& \displaystyle{{\partial \tilde{c}} \over {\partial \tilde{t}}} = -\displaystyle{1 \over C}\displaystyle{{\partial \tilde{T}} \over {\partial \tilde{t}}} \cr & \quad \quad + \displaystyle{{1-C} \over {2{\tilde{\phi }}^{1/2}}}\left[ {\displaystyle{1 \over C}{\tilde{\nabla }}^2\tilde{T} + \displaystyle{1 \over {L_e\tilde{\phi }}}\tilde{\nabla }\cdot \left( {\tilde{\phi }\tilde{\nabla }\tilde{c}} \right) + \displaystyle{{\beta \tilde{\phi }} \over {{\tilde{\phi }}_0}}\left( {\tilde{\nabla }\tilde{T} + C\tilde{\nabla }\tilde{c}} \right)\cdot \tilde{\nabla }\tilde{c}} \right] \cr & \quad \quad \left( {S\phi _0\tilde{\phi } + C\tilde{c} + \displaystyle{{1-C} \over {2{\tilde{\phi }}^{1/2}}}} \right)^{-1}.} \]Upon defining appropriate initial conditions and boundary conditions, Eqns. (17) and (18) can be solved for changes in \(\tilde T\) and \(\tilde c\) while making use of the expression for \(\tilde \phi\) from Eqn. (16), enabling the spatial and temporal evolution of the scaled Darcy flow rate to be evaluated from Eqn. (15). Appendix A describes the semi-discretized numerical scheme used to generate the results presented below.

Further intuition into the controls on vein flow is facilitated by substituting the gradient of \(\tilde \phi\) from Eqn. (16) into Eqn. (15) to obtain19\[\tilde{{\bf q}} = -\beta \,\left(1-{\cal C}\right)\left({{\tilde \phi}^{1/2}\over 2} \tilde \nabla\tilde \phi + \tilde\phi^2{R_{v0}\Delta z\over \lambda_c^2}\hat{{\bf z}}\right)\; ,\; \]which makes clear from the first term in parentheses that flow tends toward lower liquid fractions; the second term in parentheses represents a small correction that exactly offsets the reduction in liquid fraction with elevation that occurs under hydrostatic loading with uniform temperature and concentration. When circumstances produce slight increases in temperature and/or solute concentration toward the ice base and thereby result in increased liquid fractions with depth, Eqn. (19) predicts upward vein flow that has the potential to cause diffuse debris entrainment.

2.2 Parameter choices

In order to evaluate the relative importance of the different terms appearing in the governing equations, we chose scalings that make the dimensionless fields \(\tilde \phi\), \(\tilde T\), and \(\tilde c\), along with their spatial derivatives, of order unity (with bounds \(\tilde \phi \ge 0\), \(\tilde T\le 1-{\cal C}\tilde c\), and \(\tilde c\ge 0\)). We chose the time scale Δt = ρ_i C_iΔz²/K in Eqn. (10) to make the dimensionless temporal and spatial derivatives of temperature balance in Eqn. (13) as well. The values of relevant physical constants that are summarized in Table 1 imply that the Lewis number L_e from Eqn. (12) is very large, as listed together with the nominal parameter values compiled in Table 2. The other dimensionless parameters \({\cal C}\), \({\cal S}\) and β are each sensitive to one or more variables, in particular c₀, T₀, and k₀, that must be estimated for the physical context of interest.

Table 1. Physical constants

Parameter	Value	Units
Ci	2.1 × 103	J kg−1 K−1
D	4 × 10−10	m2 s−1
g	9.8	m s−2
K	2.2	W m−1 K−1
L	3.3 × 105	J kg−1
Rg	8.314	J mol−1 K−1
Tm	273	K
γ	0.029	J m−2
λc	6.1 × 10−3	m
μ	1.8 × 10−3	Pa s
ρi	920	kg m−3
ρl	103	kg m−3
Γ	0.0019	m3 K mol−1

Table 2. Nominal parameters, dimensionless variables and scales

Parameter1	Value	Range	units
\(^{\dagger} c_{B0}\)	0.05	0.01 − 10	mol m−3
c0	17	1 − 103	mol m−3
\(^{\dagger} d\)	3	1 − 10	mm
k0	5.4 × 10−14	10−17 − 10−12	m2
\(^{\dagger} R_{v0}\)	0.3	0.05 − 0.5	mm
Tm − T0	0.031	0.0019 − 1.9	K
\(^{\dagger} \psi\)	25	0 − 35	degrees
α	0.10	0.061 − 0.16	[]
ϕ0	0.003	10−5 − 0.05	[]
\({\cal C}\)	0.9975	0.75 − 1	[]
Le	2.8 × 103	n.a.	[]
\({\cal S}\)	5.0 × 103	80 − 8 × 104	[]
β	1.0	0.008 − 1.2 × 106	[]
ρp	2650	2000 − 3000	kg m−3
Δz	0.12	0.074 − 0.74	m
Δt	1.3 × 104	4800 − 4.8 × 105	s

¹ Indicates a control variable that was assigned to determine values for some of the other parameters, as described in the text.

The approach followed here to obtain self-consistent estimates of the control variables involves first assigning values for the ice grain size d and the radius of curvature of vein walls R_v0, so that the characteristic liquid fraction ϕ₀ can be obtained as noted above following Eqn. (7). Our interest in particle entrainment motivates a focus on vein radii that are sufficiently large to accommodate particles of the observed size range, with typical maximum diameters of approximately 2R ≈ 10 − 100 μm; identifying the radius of the largest circle to fit within a vein cross section as R_fit, geometrical considerations (e.g. Mader, 1992a) yield \(R_{v0}\approx ( {6\sqrt {6}/\alpha }) ^{1/2} R_{\rm fit}\approx 12 R_{\rm fit}\). For comparison, in detailed observations on ice from 7–60 m depth beneath the surface of Blue Glacier, Raymond and Harrison (1975) found that average vein cross-sectional areas measured within an hour of sample recovery implied R_v0 values near 100 μm; enhanced solute concentrations in basal ice would be consistent with larger values and accordingly we assign a nominal R_v0 = 300 μm, which is similar to the vein sizes measured in the ice permeability study by Fowler and Iverson (2022). Guidance in estimating the permeability scale k₀ comes from theoretical considerations for polycrystalline intergranular flow described by Frank (1968) suggesting that \(k_0 = \phi _0^2 d^2/\chi\), with the constant χ ≈ 2000 according to Nye and Frank (1973). It should be noted, however, that recent experimental measurements indicate a greater sensitivity of k₀ to grain size d (Fowler and Iverson, 2022) suggesting that the uncertainty in assigning permeability values that is reflected in the broad range listed in Table 2 may be somewhat conservative. Fowler and Iverson (2022) also advocate for a modestly reduced value of χ ≈ 1500 that we employed to assign the nominal value of k₀ in Table 2.

Ice recovered from natural glaciers and produced for laboratory experiments is often quite pure, as gauged on the basis of the weight of soluble impurities relative to the total. For reference, distilled water typically has less than 0.5 ppm by weight total dissolved solids, while municipal tap water might be expected to exceed this concentration by two or three orders of magnitude. An analysis of tap water from the laboratory in Madison, WI, for example, found 8 ppm Ca, 4 ppm Mg, and 4 ppm Na; whereas impurity levels in the distilled water that was frozen for our ring-shear experiments were below the detection limit of 0.02 ppm for each of these species. Impurity levels measured in polar ice cores vary, but are typically similar to those for distilled water except in the bottom few meters of basal ice, which can have impurity concentrations that are several orders of magnitude higher (e.g. Hallet and others, 1978; Knight, 1997). Defining the bulk concentration of mobile impurities on a molar basis as c_B0, we expect the dissolved concentration in the vein liquid to be given by c₀ = c_B0/ϕ₀, where the tendency for c_B0 and ϕ₀ to covary limits the range reported in Table 2. The nominal value of c_B0 that we report in Table 2 is an order of magnitude lower than the tap water concentrations noted above, yet much less pure than the distilled water that was introduced to the loess bed before sitting and interacting with particle surfaces for days prior to the start of the ring-shear experiment.

Using the equilibrium relation from Eqn. (4), while accounting for the effects of surface energy γ and vein-wall curvature in setting P_i = P_l + γ/R_v0, and taking P_i = P_m leads to an estimate for the temperature scale of T_m − T₀ = T_mγ/(ρ_l L R_v0) + Γc₀. We note that reasonable choices for the length scale Δz might be based on factors like the typical size of bedrock obstacles or the dimensions of ice samples inserted within an experimental apparatus; a natural choice in the absence of such considerations could involve setting \(\Delta z = \lambda _c^2/R_{v0}\) (approximately 0.1 m for the nominal R_v0 listed in Table 2 and spanning from 0.07 − 0.7 m for the listed range of R_v0 values), which simplifies the final term in Eqn. (16) to more clearly indicate that this is the length scale over which changes in gravitational potential significantly affect vein geometry. The nominal choice of Δz made in this way leads to a thermal diffusion time scale Δt of about 3.5 hours.

In order for vein cross-sections to be large enough to convey debris particles that are tens of microns in size, dissolved impurity levels must generally be much more effective than surface energy in generating vein liquid and hence the nominal value of \({\cal C}\) listed in Table 2 is very close to unity, as is its typical range. For these estimates of \({\cal C}\), Eqn. (15) implies that temperature gradients and concentration gradients are of comparable importance for driving liquid flow. The product \({\cal S} \phi _0\) in Eqn. (13) is of order unity for the nominal parameter values, suggesting that the sensible heat associated with changes in temperature is similar in scale to the latent heat associated with changes in phase. Note, however, that the broad range of typical values listed for ϕ₀ and \({\cal S}\) in Table 2, allow for conditions where the product \({\cal S}\phi _0\) is either large or small compared to unity, in the former case describing a regime in which temperature changes are primarily limited by the removal of latent heat, whereas in the latter case sensible heat is more important. Meanwhile, the enormous Lewis number ensures that \(L_e^{-1}\) is always very small, suggesting that compositional diffusion is comparatively unimportant in Eqn. (14). Note however that the relative importance of advective and diffusive compositional transport depends not only on the Lewis number, but also on the magnitude of the transport rate; we show below that the counteracting effects of temperature and concentration gradients often conspire to make \(\tilde {{\bf q}}\) so small that diffusion actually dominates compositional transport despite the large value of L_e (this is consistent with the a priori assumption made by Nye, 1991, in his idealized treatment of vein evolution). The nominal value listed for the transport factor β is of order unity, though the potential variations in the permeability and undercooling allow for a particularly broad range of values for this parameter. While it is tempting to jettison the small terms, for example by expanding the governing equations and retaining only the leading order dependence on the small parameter \(1-{\cal C}\), this procedure doesn't yield an obvious numerical advantage so we retain everything for now. We observe, however, that in the small \(1-{\cal C}\) limit, Eqn. (18) implies that changes in the scaled concentration are very nearly opposite to changes in the scaled temperature, with notable consequences for limiting the efficacy of vein flow, as discussed below.

3. Vein flow forced by changes in basal temperature

The generalized Clapeyron equation that describes equilibrium between ice and liquid within the veins also applies at the glacier bed, here assigned location z = 0. Along the ice-walled surfaces overlying macroscopic drainage elements (e.g. channels, cavities), the ice and liquid pressures must equal each other and with P_i = P_l = P_m Eqn. (4) implies that the scaled temperature increases above its nominal value of \(\tilde T = 0\) to reach \(\tilde T = 1-{\cal C} \tilde c\). (We note that cases with P_i = P_l < P_m along z = 0 must commonly occur as a result of stress bridging that allows P_i above macroscopic drainage elements to be lower than the overburden stress (e.g. Rempel and others, 2022); the warmer values of \(\tilde T\) that represent such circumstances can be determined using Eqn. (4), and the corresponding changes in ice pressure distribution above the bed can be incorporated into a generalized treatment of changes in liquid fraction, with appropriate modifications to Eqn. (16).) In regions where the basal ice approaches near enough to enable stress transmission across thin premelted liquid films to the underlying till or bedrock, however, the ice pressure must exceed the liquid pressure (e.g. Dash and others, 2006; Rempel, 2008), much as it does along the boundaries of veins, where surface energy is responsible for maintaining the ice–liquid pressure difference. This suggests an idealized model scenario in which ice that is initially at a uniform temperature, with \(\tilde T( \tilde t = 0,\; \, \tilde z) = 0\), containing veins with a uniform composition, \(\tilde c( \tilde t = 0,\; \, \tilde z) = 1\), encounters an abrupt temperature change along its basal boundary so that \(\tilde T( \tilde t,\; \, \tilde z = 0) = \nu ( 1-{\cal C}\tilde c)\) for a chosen ν ≤ 1, while the boundary solute concentration is held fixed at \(\tilde c( \tilde t,\; \, \tilde z = 0) = 1\). With ν = 1 this might represent the conditions obtained when ice slides over a subglacial cavity, for example, whereas values of ν < 1 could represent conditions resulting from changes in effective stress that are caused by transient behavior in the basal hydrological system. Far-field boundary conditions of \(\tilde T( \tilde t,\; \, \tilde z_{\rm max}) = 0\) and \(\tilde c( \tilde t,\; \, \tilde z_{\rm max}) = 1\) are applied at some distance \(\tilde z_{\rm max}\) beyond the range at which thermal and compositional changes are able to propagate in the idealized one-dimensional results described below.

Figure 2 summarizes some of the model predictions for the near-basal effects of an abrupt change in bed temperature. Immediately afterward, the ice above the bed is warmed and this change in temperature diffuses upward as time progresses, resulting in the different profiles shown with black lines in Fig. 2a. Coupling of the thermal and compositional fields through their mutual control on the changing liquid fraction cause the solute concentration (red lines) to evolve in response to changes in temperature. In particular, slight increases in \(\tilde \phi\) that take place with increases in \(\tilde T\) cause local freshening of vein liquid and corresponding drops in \(\tilde c\) (displayed as \(\tilde c-1\) in the profiles of Fig. 2a and the time histories of Fig. 2b to facilitate comparison with the simultaneous changes in \(\tilde T\) on the same axes). The tendency for the time evolution of \(\tilde T\) and \(\tilde c\) to mirror each other makes the contributions of \(\tilde \nabla \tilde T\) and \(\tilde \nabla \tilde c\) to driving liquid flow in Eqn. (15) nearly cancel. However, as described by Eqn. (18), the compositional and thermal evolution are not perfectly matched, and compositional diffusion toward the ice interior from the comparatively solute-rich fluid near the bed does eventually allow the concentration to return toward its boundary value, as shown by the time evolution of \(\tilde c-1\) for the lowest layer (solid red) depicted in Fig. 2b and the spatial variation of \(\tilde c-1\) at the longest times in Fig. 2a. As might have been anticipated, Fig. 2a shows that a concentration profile sampled at a particular time exhibits an abrupt decrease with elevation close to the bed (at a location controlled by the rate of compositional diffusion), which is separated from an abrupt increase with elevation higher up (at a location that is controlled by the rate of thermal diffusion); the distance between these changes scales approximately as \(\sqrt {L_e}\approx 50\). Similarly, the time history shown in Fig. 2b for the lowest displayed location, 0.5 cm above the bed, exhibits a drop in concentration that is associated with processes of thermal diffusion that precedes a later rise in concentration that is associated with compositional diffusion; the time delay between these changes scales approximately as L_e ≈ 2800.

Figure 2.

Perturbed near-bed conditions following an abrupt temperature rise from \(\tilde T = 0\) to \(\tilde T = 0.99( 1-{\cal C})\) along \(\tilde z = 0\), modeled using the nominal parameter values listed in Table 2. (a) Profiles of scaled temperature \(\tilde T\) (black) and change in dimensionless concentration \(\tilde c-1\) (red) as a function of height above the bed at the different labeled times. (b) Evolution of scaled temperature \(\tilde T\) (black) and change in dimensionless concentration \(\tilde c-1\) (red) at the different labeled elevations. Note that the axis limits of \(\tilde T = \pm 2\times 10^{-3}\) represent excursions in T away from \(T_0\approx -0.031\, ^\circ\)C of magnitude \(6.3\times 10^{-5}\, ^\circ\)C, while the same axis limits on \(\tilde c-1\) represent deviations in c from c₀ ≈ 17 mol m⁻³ of magnitude 0.03 mol m⁻³.

[Figure omitted. See PDF]

The nominal vein radius R_v chosen for these calculations is of sufficient size to permit the passage of silt-sized particles (as noted above, Mader, 1992b, 1992, found that \(R_v\approx \sqrt {6\sqrt {6}/\alpha }\, R_{\rm fit}\approx 12 R_{\rm fit}\), so that particles a factor of 12 smaller than the radius of curvature R_v can fit within veins; for example R_fit ≈ 25 μm for the nominal R_v0 and α). Whether such particles can be dragged by the flowing vein liquid and lofted above the bed depends on the speed of liquid flow. Recognizing that the average velocity of water molecules is given by the ratio of the Darcy flow rate to the liquid fraction, we approximate the radius R of the largest particle of density ρ_p that can be lifted by this flow as corresponding with the Stokes’ settling velocity, yielding20\[R\approx \sqrt{{9 \mu \vert {\bf q}\vert \over 2\left(\rho_p-\rho_l \right)g \phi}} = \sqrt{{9 k_0\left(\rho_l-\rho_i\right)\tilde \phi\left\vert \tilde \nabla \tilde T + {\cal C}\tilde \nabla \tilde c\right\vert \over 2\phi_0\left(\rho_p-\rho_l\right)\left(1-{\cal C}\right)}}\; .\]We note that Eqn. (20) represents a conservative estimate of the largest entrained particles, since it is based on the expression for particle drag that determines the steady speed of a particle relative to an unbounded fluid (i.e. for R ≪ R_fit), whereas drag is enhanced during flow in confined geometries (e.g. Humphrey and Murata, 1992); we revisit this complication later. Under the idealized conditions where Eqn. (20) applies, Fig. 3 shows the radii of particles that can be carried by the liquid flow through veins, with upward transport indicated in black and downward transport in red, which we model as also following Eqn. (20). Immediately after the basal temperature perturbation, rates of liquid flow are very low because changes in temperature and solute concentration combine to keep liquid potential gradients small. Over time, the differing rates of heat and solute transport facilitate increases to both the driving potential gradient and the liquid fraction so that the Darcy flow rate increases. The faster liquid flow is able to mobilize larger particles, as shown by the profiles in Fig. 3a and the increases in R over time at particular elevations in Fig. 3b. This late-stage particle transport is delayed until compositional diffusion is sufficiently effective at transporting solute from the basal interface to diminish concentration gradients appreciably, enabling residual temperature gradients to more effectively drive liquid flow. Comparison against the corresponding values of R_fit (not shown) confirm that particle entrainment for this scenario is limited by the rate of liquid flow rather than the physical dimensions of the veins themselves.

Figure 3.

Approximate radii R of the largest particles that can be carried upward by modeled rates of liquid flow through the vein network, treating the fluid drag on particles as corresponding with that for Stokes’ settling velocity and evaluated using Eqn. (20) for the nominal parameters listed in Table 2. Black lines represent upward flow, whereas red lines are for downward flow toward the bed. (a) Profiles of R at the different labeled times. (b) The evolution of R in the lowermost basal ice. Note that Eqn. (20) implies that the average velocity of vein liquid |q|/ϕ is quadratic in R such that |q|/ϕ ≈ 2.0 × 10⁻⁴ μm/s along the lower axis limit of R = 10⁻² μm, rising to |q|/ϕ ≈ 1.8 mm/s at the upper axis limit of R = 30 μm.

[Figure omitted. See PDF]

4. Vein flow accompanying interior relaxation

Diffuse debris entrainment occurred during a controlled ring-shear experiment that slid a temperate ice ring over a bed of saturated loess. Analysis following the conclusion of the experiment confirmed that diffuse debris entrainment had led to the deposition of loess particles along grain boundaries (see Fig. 4), but the imposed experimental conditions make it unlikely that this behavior was forced by an abrupt basal temperature change. This recognition led us to investigate whether transient evolution of the vein network that began prior to the initiation of slip might have primed the system for particle entrainment.

Figure 4.

Entrained debris following a ring-shear experiment sliding ice over a bed of Iowa loess. (a) The base of the ice ring with a patchy frozen ‘fringe’ developed during sliding. (b) Side view of the near-basal region (flipped so that the \(\hat {{\bf z}}\) direction is downward, as labeled) with a thin fringe layer that was overlain during the experiment by ice that is discolored from diffuse concentrations of entrained debris. (c) Close-up of the fringe and lowermost debris-laden ice. (d) Thin section from near the fringe boundary viewed with cross-polarized light to show debris concentrated along grain boundaries.

[Figure omitted. See PDF]

For context, we provide a brief description of the cryo-ring-shear experiment and refer the reader to Hansen and Zoet (2022) for further detail regarding equipment specifications and protocols. For this experiment, a ~15 cm ring of ice emplaced atop an 11 cm thick bed of saturated Iowa loess (median particle diameter 48 μm, distribution shown in Fig. 5c) was loaded under uniaxial compression. A vertical axial stress (mean ~176 kPa) was applied to the sample chamber by a hydraulic ram for ~20 days, and two pressure transducers installed in the sidewalls 7 cm below the ice–bed interface recorded water pressures throughout the experiment (mean P_w ≈ 89 kPa). Temperature was monitored by four RTD thermistors embedded in the sample chamber base-plate and eight RTDs in the acrylic sidewalls. Heat flux into the sample chamber was regulated with a temperature-controlled glycol–water bath. Twelve drainage lines connected the sample chamber to atmospheric pressure, allowing excess porewater to drain from the system.

Figure 5.

Data from a ring-shear experiment sliding ice over Iowa loess. (a) Frequency distribution of ice grain diameter measured on a vertical thin section containing ~320 grains. (b) Sampled volumetric debris concentration plotted as a function of distance from the sliding surface. Blue line shows approximate location of the fringe upper boundary. (c) Particle size distribution of the Iowa loess (blue, solid) shown with an approximate lognormal fit (black, dashed) indicating a median diameter of exp μ ≈ 48 μm. (d) Median radii of entrained particles found at diffuse concentrations and shown as a function of distance above the sliding interface.

[Figure omitted. See PDF]

At the onset of a ring-shear experiment, the upper surface of the ice ring must be frozen into the notched Delrin® platen that drives shear. To facilitate this coupling, the lower half of the ~15 cm thick ice ring was maintained near the bulk melting temperature through heat exchange with the circulating bath that was lowered to a position ~7.5 cm down the ice ring. Simultaneously, the upper half of the ice ring was exposed to an ambient freezer temperature of \(-4\, ^\circ\)C and an axial force was applied. Once the ice ring was securely affixed to the upper platen, the level in the circulating bath was raised to submerge the full thickness of the sample chamber so that the upper half of the ice ring warmed to approach the bulk melting temperature. The drive was engaged at this point, sliding the temperate ice ring atop the loess bed at a constant centerline velocity of 73 m a⁻¹.

At the conclusion of the experiment, we observed a ~0.5 − 1 cm thick layer of frozen debris (i.e. a ‘frozen fringe’, see Rempel, 2008) adhered to the base of the ice (Fig. 4a), but also a diffuse fraction of sediment entrained above the upper boundary of the fringe (Figs. 4b,c; note that the ice ring was turned over so the fringe is at the top of these images). Thin sections cut from the ice ring in this diffuse zone reveal particles concentrated along grain boundaries (Fig. 4d). To assess how debris concentration varied with distance from the ice–bed interface, a hot wire was used to section off horizontal segments of the ice ring (~15 cm), which were subsequently weighed and melted. Above the upper boundary of the fringe, volumetric debris concentrations y (in %) decreased monotonically (Fig. 5b), well described by a power-law of the form y ≈ 12.3 z^−2.051 (with z in cm measured from the sliding surface along the fringe base). As shown in Fig. 5c, the loess bed had an approximately lognormal particle size distribution, and the median radius of (expμ)/2 ≈ 24 μm was noticeably higher than the median size of entrained particles, which is shown as a function of distance from the sliding interface in Fig. 5d.

The evolution of vein temperature, solute concentration and liquid fraction for this experimental scenario are expected to satisfy the governing equations described above, but with different initial conditions and boundary conditions. For simplicity, in our idealized model treatment we considered one dimensional transport in the vertical direction, assigning the location \(\tilde z = 0\) to the ice base and \(\tilde z = \tilde z_{\max } = 0.15\, {\rm m}/\Delta z\) to the upper boundary. We again used the nominal parameter values from Table 2, basing d on the average observed grain size from Fig. 5a, and imposed a uniform initial solute concentration so that \(\tilde c( \tilde t = 0,\; \, \tilde z) = 1\). To approximate the initial thermal state we set \(\tilde T( \tilde t = 0,\; \, \tilde z< 0.5\tilde z_{\rm max}) = 0\) and \(\tilde T( \tilde t = 0,\; \, \tilde z> 0.5\tilde z_{\rm max}) = -4\, ^\circ {\rm C} /( T_m-T_0)\). On the bottom boundary we imposed constant temperature and concentration conditions with \(\tilde T( \tilde t,\; \, \tilde z = 0) = 0\) and \(\tilde c( \tilde t,\; \, \tilde z = 0) = 1\), whereas we treated the upper boundary at \(\tilde z_{\rm max}\) as thermally insulating so that \(\partial \tilde T/\partial \tilde z = 0\) and imposed the same uniform boundary concentration as before so that \(\tilde c( \tilde t,\; \, \tilde z_{\rm max}) = 1\).

Figure 6a shows model predictions for profiles of the scaled temperature and concentration at the labeled discrete times as the vein network relaxes from the perturbed initial conditions. As expected, the temperature (black) evolves quickly toward a steady-state \(\tilde T = 0\) that is consistent with the imposed boundary conditions. In the lower half of the ice ring (i.e. below z = 7.5 cm), changes in the solute concentration (red) mirror the changes in temperature, which is similar to the behavior shown earlier in Fig. 2a. The rapid warming that takes place above the midpoint, where the ice was initially cold (\(\tilde T( t = 0,\; \, z> 7.5\, {\rm cm}) \approx -130\), corresponding with \(T = -4^\circ\)C), causes the vein liquid content to increase, thereby diluting the solute and causing \(\tilde c\) to approach its limiting value of zero above z = 7.5 cm. As shown by the evolution of scaled temperature and concentration in Fig. 6b, this region persists as being much fresher than the lower half of the ice ring for many days of model time. In the lower half of the ice ring, Fig. 6b shows how the response of the scaled temperature and concentration to the changing conditions above is both enhanced with proximity to the ring midpoint and increasingly delayed as the ice base is approached.

Figure 6.

Modeled conditions for ring-shear experiments. (a) Profiles of scaled temperature (black) and concentration (red) at the different labeled times. (b) Evolution of scaled temperature (black) and concentration (red) at the different labeled heights above the bed. (c) Maximum particle radius that can be mobilized by flow upward into the ice (black) and downward toward the ice base (red) when particle drag is assumed to correspond with Stokes’ settling velocity, displayed for the different labeled distances above the ice base. (d) Profiles of maximum particle radius that can be mobilized by upward (black) and downward (red) flow at the different labeled times. In (a) and (b) note that the axis limits of \(\tilde T = \pm 40\) represent excursions in T away from \(T_0\approx -0.031\, ^\circ\)C of magnitude \(1.26\, ^\circ\)C, while the same axis limits on \(\tilde c\) represent deviations in c from c₀ ≈ 17 mol m⁻³ of magnitude 667 mol m⁻³. In (c) and (d), the average velocity of vein liquid that entrains particles of the size given by the lower axis limit of R = 10⁻² μm is |q|/ϕ ≈ 2 × 10⁻⁴ μm s⁻¹, while R = 10 μm requires |q|/ϕ ≈ 2 × 10² μm s⁻¹. (Note, in panel (a), the predictions for 1 day and 10 days nearly coincide and give the appearance of vertical long-dashed lines.).

[Figure omitted. See PDF]

We show how the size of the largest particles that can be mobilized by the flowing vein liquid evolves in Fig. 6c and display profiles of this maximum size at different times in Fig. 6d. Below the ring midpoint, the more rapid response to changing conditions nearer to the warming region above z = 7.5 cm causes upward transport to first be capable of carrying micron-sized particles at higher elevations, but within an hour our model predicts that vein transport is capable of mobilizing micron-sized particles up into the ice from the bed (i.e. see the dashed black line in Fig. 6d). The predicted flow direction reverses later so that by the time one day has elapsed, the particle transport direction is downward throughout the lower half of the ice ring. Notably, over the time frame represented in Fig. 6c at each horizon the maximum particle size that can be carried upward at early times is larger than the maximum particle size that can be carried downward later on – implying that the outcome of this vein flow history is capable of depositing particles within the vein network. Above the mid-point, upward transport is maintained and increases in magnitude for the model time displayed in Fig. 6c.

4.1 Model–data comparison

Consistent with the results of earlier experiments (e.g. Knight and Knight, 1999), the median radii of entrained debris particles are smaller than the median radii of bed particles from which they are drawn (see Figs. 6c,d). Nevertheless, the predicted maximum sizes of entrained particles shown in Fig. 6 are smaller than the observed vein particles documented in Fig. 5d. These discrepancies might be explained if vein flow is faster than predicted by our idealized model, for example if the transport factor β is larger than that which corresponds with the adopted nominal parameter values. Both the impurity loading that determines c₀ and controls T_m − T₀, and the permeability scale k₀ are amongst the least well-constrained parameters for describing experimental conditions. Since β in Eqn. (12) is proportional to their product, it is conceivable that higher rates of liquid flow might be explained by the actual values of k₀ and c₀ being higher than the nominal values summarized in Table 2. An increase in β by an order of magnitude, for example, would raise predictions for the maximum entrained particle size by a factor of \(\sqrt {10}\approx 3\) so that particles of the median size reported in Fig. 5d would be carried upward during the early model stages illustrated in Fig. 6c and d. We note that recent experimental results reported by Fowler and Iverson (2022) yielded vein permeabilities that were indeed on the high end of previously published estimates. Moreover, in addition to the possibility that our idealized model may underestimate vein flow rates, we recognize as well that the drag on a particle in an unbounded fluid can be much lower than the drag associated with particle transport in Poisseuille flow. Humphrey and Murata (1992) summarize theoretical and experimental estimates of particle settling velocities in circular tubes. When the particle radius is an appreciable fraction of the tube radius, the Stokes’ settling velocity can be much larger than measured and predicted particle velocities. For example, results compiled by Humphrey and Murata (1992) suggest that a given rate of flow through a circular tube can entrain a spherical particle that is an order of magnitude larger than would be predicted by Eqn. (20) once the particle radius approaches ~85% of the tube radius (a three-fold enhancement is seen when the particle radius is approximately half the tube radius – this effect is roughly equivalent to the behavior noted above for the case when β is raised by an order of magnitude). We are not aware of analogous work examining particle transport velocities in triangular, or vein-shaped tubes, but the implication is clear: particle transport is enhanced when fluid flow is restricted significantly by partial blockage of the fluid conduit, making the Stokes’ settling limit encapsulated in Eqn. (20) a conservative estimate for the maximum size of entrained particles at a given q.

Prior to the loess experiment, we followed the procedure described above to examine fringe growth during shear over a bed of Horicon till. The median grain size of the till was approximately 350 μm, which we infer to have been too large for vein entrainment as the fringe–ice interface remained sharp, with no evidence of diffuse entrainment. After the loess experiment, we used the ring-shear apparatus to conduct yet another experiment with slip over glass beads that had a median grain size of approximately 160 μm. High resolution video evidence obtained through the transparent sample chamber wall clearly showed particle movement upward into the ice early on during this experiment, consistent with the model predictions discussed above. However, upon extraction of the ice ring at the conclusion of this experiment, we observed that the fringe–ice interface was sharp, with negligible concentrations of diffuse debris remaining in the overlying ice. This suggests that the predicted reversal in the direction of vein flow, discussed above, may have been responsible for removing the particles that we observed moving upward during the early stages of the glass-bead experiment.

An important model limitation concerns valid questions as to whether or not Eqn. (20) really sets the maximum size of particles that can be mobilized by downward flow. As long as a particle's density exceeds that of water, it will have a tendency to sink and rest against the base of an inclined vein, where it might be expected to encounter a flow rate that on average would be slower than |q|/ϕ because of the proximity to vein walls. Perhaps more importantly, the acceleration of gravity is aligned with the prevailing direction of vein flow when it is oriented downward. As a result, particles might be expected to simply sink downward toward the ice base even in the absence of significant downward vein flow. In reality, flow pathways, whether conveying transport that is on average upward or downward, are tortuous and expected to involve a range of vein orientations from vertical through horizontal, with differing local flow speeds that are captured by a Darcy-type flow law only as an overall average transport rate in the direction of the macroscopic hydraulic potential gradient. These considerations suggest that the particle radius described by Eqn. (20) is best regarded as an estimated scale for the typical size of a particle that can be induced to move in the direction of vein flow. Individual particles along particular vein transport paths can be expected to experience conditions that deviate from the average behavior in ways that a continuum model such as ours is not designed to describe. Keeping these complications in mind, we are encouraged, nevertheless by the qualitative agreement between experimental observations and the predictions of our idealized model, which suggest that the controlling physical effects are captured reasonably well by our continuum treatment.

5. Effects of non-hydrostatic ice pressure gradients

In the idealized scenarios described above, the ice pressure was treated as hydrostatic. To address more general circumstances where this approximation does not apply, we can introduce a non-hydrostatic pressure perturbation \({\cal P}\), defined so that the ice pressure is \(P_i = P_m-\rho _i g z + {\cal P}\). The scaled liquid fraction from Eqn. (16) is modified under these perturbed conditions to21\[\tilde \phi = \left({1-\tilde T -{\cal C} \tilde c\over 1-{\cal C}} + {R_{v0}\Delta z\over \lambda_c^2}\tilde z-{R_{v0}{\cal P}\over \rho_i g \lambda_c^2}\right)^{-2},\; \]while the scaled Darcy flow rate from Eqn. (15) becomes22\[\tilde{{\bf q}} = -\beta \tilde \phi^2\left[\tilde \nabla \tilde T + {\cal C} \tilde \nabla \tilde c + {\left(1-{\cal C}\right)\rho_lR_{v0}\over \left(\rho_l-\rho_i\right)\rho_i g \lambda_c^2} \tilde \nabla {\cal P}\right]\; .\]Equation (21) suggests that non-hydrostatic pressure changes have a significant effect on the liquid fraction when they are comparable to \(\rho _i g \lambda _c^2/R_{v0}\), which is approximately 1.1 kPa for the nominal value of R_v0 from Table 2. The tendency for non-hydrostatic pressure gradients to drive liquid flow is made explicit by the final bracketed term in Eqn. (22); however, we note that the direct dependence of the multiplicative coefficient on \(1-{\cal C}\), which is expected to be very small under conditions where veins are enlarged enough by the presence of solutes to accommodate micron-scale particles, limits the importance of non-hydrostatic pressure gradients for driving liquid flow that entrains particles, in comparison to the effects of temperature and concentration gradients. Nevertheless, since non-hydrostatic pressure gradients drive ice deformation and perform work, they can have a significant indirect influence by acting as a heat source that alters the temperature field with consequences for both of the other two gradient terms that combine to determine the bracketed liquid potential gradient that drives flow on the right side of Eqn. (22). A more complete assessment of the consequences of non-hydrostatic ice pressure gradients for driving vein liquid flow and associated debris entrainment could be made with a more sophisticated treatment that is coupled to the ice deformation problem. Previous coupled models have treated vein transport as a consequence of viscous compaction, while ignoring the influence of surface energy and solute content (e.g. Schoof and Hewitt, 2016); we leave the problem of reconciling these different approaches to future work.

6. Discussion and Conclusions

Our idealized continuum model confirms that vein liquid flow can entrain and transport silt-sized particles up from the glacier bed into the lowermost basal ice. Predicted liquid flow velocities are reduced by the tendency for gradients in temperature and solute concentration to counteract each other and limit the overall magnitude of the driving gradient in liquid potential (see Fig. 2). However, because thermal and solute transport rates differ, in response to perturbed bed conditions, temperature and concentration gradients undergo different rates of change; as a result, the vein liquid fraction can grow and cause flow to accelerate so that the radii of the largest particles that can be transported increase with time elapsed since the bed perturbation occurred (see Fig. 3). The predicted growth in vein liquid fraction that helps facilitate increasing flow rates is in qualitative agreement with the asymptotic analysis of changes in vein radii by Nye (1991), who considered an idealized system in which the density difference between ice and liquid water was neglected so there was no liquid flow. Because small changes in basal temperature are an inevitable consequence of anticipated heterogeneous changes in the difference between the ice normal stress and the subglacial liquid pressure, this mechanism can be expected to lead to widespread diffuse debris entrainment that complements other basal mass-exchange processes (e.g. Alley and others, 1997; Rempel and others, 2022).

Examining a scenario motivated by the relaxation of temperature and solute concentration that occurred at the beginning of our ring-shear experiments, our model suggests that vein flow can lead to a pattern of particle entrainment that is in qualitative agreement with detailed observations. For our nominal parameter choices and the conservative assumption that particle drag can be approximated from analysis of the terminal settling velocity in an unbounded viscous fluid, the predicted maximum sizes of entrained particles are smaller than the largest particles entrained during our experiments (see Figs. 5 and 6). These discrepancies could be an indication that the permeability of polycrystalline ice under given conditions is higher than assigned by our treatment, or it could result from our neglect of the enhancement to particle drag that results from requiring flow to intrude between the narrow gaps that separate particle surfaces from vein walls. Ongoing experimental efforts to measure the permeability of polycrystalline glacier ice (e.g. Fowler and Iverson, 2022) are expected to reduce the uncertainty in assigning values for this important parameter. Relatively straightforward model extensions could be made to include the dependence of particle drag on vein size, though it should be noted that detailed observations of particle settling in capillary tubes (e.g. Humphrey and Murata, 1992) suggest additional complications, including nontrivial sensitivities to tube angle and particle position (relative to the tube axis) that might combine to add a stochastic element and further challenge the ability of continuum formulations to make accurate predictions. It is worth noting as well that our model formulation using Darcy's law to describe vein transport is appropriate over length scales that are much larger than the typical diameter of ice grains, and this condition is violated near the basal interface where the most rapid vein liquid transport is expected to occur. Moreover, we recognize that predicted spatial gradients in flow magnitudes, sometimes including reversals in flow directions, must entail corresponding ice deformation in order to conserve mass; for the low liquid fractions and modest flow rates that are typical, we anticipate that the errors associated with neglecting this complication are likely to be small, but note that the non-hydrostatic ice pressure gradients that develop in consequence are of interest for future work.

When the gradient in ice pressure deviates from hydrostatic equilibrium, liquid flow is drawn toward regions of anomalously low ice pressure (i.e. described by Eqn. (22)). Ice deformation takes place as well under such nonhydrostatic stress states and viscous dissipation provides a heat source that can affect temperature gradients, thereby further altering the pattern of vein liquid flow. Hence, a rigorous model treatment of such scenarios requires an analysis that couples the flow of both ice and liquid (in addition to heat and solute). Such thermomechanical coupling is expected to be important not only for driving vein flow and enabling diffuse debris entrainment, but also for helping to control the liquid volume fraction and thereby affecting the ice softness that influences both the form drag that limits basal sliding (e.g. Weertman, 1957; Gudmundsson, 1997) and the deformation that takes place in shear margins (e.g. Schoof and Hewitt, 2016; Minchew and others, 2018; Haseloff and others, 2019).