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, Reference Nye and Frank1973)

(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, Reference Nye1991; Rempel, Reference Rempel2005; Dash and others, Reference Dash, Rempel and Wettlaufer2006)

(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) becomes

(5)$$\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 (Reference Nye and Frank1973), 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 as

(6)$${\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, Reference Nye1991)

(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, Reference Mader1992a). 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 becomes

(8)$$\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 as

(13)$${\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 as

(15)$$\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) becomes

(16)$$\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) as

(17)$$\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 obtain

(19)$$\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

Table 2. Nominal parameters, dimensionless variables and scales

¹ 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, Reference Mader1992a) 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 (Reference Raymond and Harrison1975) 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 (Reference Fowler and Iverson2022). Guidance in estimating the permeability scale k ₀ comes from theoretical considerations for polycrystalline intergranular flow described by Frank (Reference Frank1968) suggesting that $k_0 = \phi _0^2 d^2/\chi$, with the constant χ ≈ 2000 according to Nye and Frank (Reference Nye and Frank1973). It should be noted, however, that recent experimental measurements indicate a greater sensitivity of k ₀ to grain size d (Fowler and Iverson, Reference Fowler and Iverson2022) 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 (Reference Fowler and Iverson2022) 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, Reference Hallet, Lorrain and Souchez1978; Knight, Reference Knight1997). 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, Reference Nye1991, 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, Reference Rempel, Meyer and Riverman2022); 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, Reference Dash, Rempel and Wettlaufer2006; Rempel, Reference Rempel2008), 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⁻³.

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, Reference Mader1992b, Reference Mader1992, 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, yielding

(20)$$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, Reference Humphrey and Murata1992); 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.

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.

For context, we provide a brief description of the cryo-ring-shear experiment and refer the reader to Hansen and Zoet (Reference Hansen and Zoet2022) 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.

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, Reference Rempel2008) 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.).

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.

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 to

(21)$$\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) becomes

(22)$$\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, Reference Schoof and Hewitt2016); we leave the problem of reconciling these different approaches to future work.