Mathematical Symbols
1 Introduction
Freshly exposed bedrock in front of retreating glaciers often exhibits features indicative of subglacial chemical alteration. In regions composed of limestone, a prominent example is the set of well-developed surface furrows and carbonate deposits, occurring in close association with undulations in the bed topography (Fig. 1). The appearance of these furrows and deposits is most strikingly revealed in low-angle sunlight.
Decimetre-scale bumps of the deglaciated limestone bedrock surface infront of Blackfoot Glacier, Montana, U.S.A., adorned with carbonate deposits (light patches) and sets of solutional furrows (arcuate patterns, concave down-glacier). Photograph taken in the direction of former ice flow.
Following earlier work by Reference HalletHallet (1976, Reference Hallet1979), such features have generally been attributed to the subglacial diagenesis of calcite. The active agent is water in the regelation film flow that is sustained as basal ice moves past (small) obstacles in the bed. On the stoss side of a bump, ice melting produces chemically aggressive water, which dissolves the bed to form furrows. Freezing of this water on the bump’s lee side elevates its solute concentration to the point that CaCO3 precipitates. This description accounts for the observed pattern well, especially the gross locations of dissolution and precipitation. However, the morphological characteristics of the assemblage remain intriguing, and merit explanation.
It is, for instance, common to find solutional furrows organized into distinct sets of arcs, with characteristic spacing and geometry. Moreover, the calcite deposits adopt a fluted surface parallel to the former direction of ice flow, and in some cases they even develop pin-shaped protrusions (“spicules”). The regularity and abundance of such forms are suggestive of specific mechanisms that underlie their formation. Their origins are considered in this paper. We address two particular questions: Why do solutional furrows appear in the first place, and why is the surface of the deposits fluted? Drawing on results of recent laboratory experiments (Reference Killawee, Fairchild, Tison, Janssens and LorrainKillawee and others, 1998), we argue that solutional furrows are intimately linked to the exsolution of gaseous CO2 and heterogeneous inclusion of it into basal ice during the freezing of regelation film water. We then propose an explanation for the fluted surface of the deposits in terms of an interfacial instability. We formulate quantitative models to describe these mechanisms.
In section 2 the relevant background material is reviewed, and we begin by developing a simplified description of how the bed surface evolves when dissolution and deposition are coupled to regelation physics. The generation mechanisms of furrows and flutes, respectively, are then treated in sections 3 and 4, where we also consider spicules as a special case of the depositional instability. Apart from the intrinsic fascination of these landforms, the current study raises the possibility of inferring the conditions that favour their formation, those which may be representative of the former environment at the ice—bed interface. We consider related issues and wider glaciological implications towards the end.
2 Background
21 Surface morphology
Figure 2 illustrates in plan view a limestone specimen from the foreground of Castleguard Glacier in the Canadian Rockies. Its surface pattern is typical of those observed elsewhere where calcareous bedrock has undergone substantial chemical modification (e.g. Reference Ford, Fuller and DrakeFord and others, 1970; Reference Sharp, Tison and FierensSharp and others, 1990). The assemblage, consisting of furrows and deposits adjacent to a bump, represents a basic unit that is often repeated over large bedrock areas (Reference Walder and HalletWalder and Hallet, 1979; Reference Hallet and AndersonHallet and Anderson, 1981) (Fig. 1).
The limestone specimen of Reference HalletHallet (1979) showing solutional furrows (A) orientated transversely to the former ice-flow direction (which is from left to right), and fluted calcite deposits (light areas, B) orientated parallel to it. A1, transverse furrows; A2, oblique furrows.
Solutional furrows
Solutional furrows are elongated sub-millimetre to millimetre-sized channels etched into the rock surface on the up-glacier side of bedrock bumps (Fig. 2, region A). They occur in sets with a millimetre-scale spacing and are oriented transversely to the former ice-flow direction, curving round the sides of bumps to form an arcuate pattern, generally concave down-glacier. This pattern apparently reflects the two-dimensional water flow field during regelation. The alignment and varying tortuosity of the solutional furrows strongly suggest an origin from chemical dissolution of the substrate, in contrast with abrasion.
Neighbouring furrows generally parallel one another, but truncations and intersections are common. Often, directly upstream of bumps, individual furrows are difficult to identify and the surface appears “scalloped” (Fig. 2, region A1). Where solutional furrows turn to point downstream, they tend to straighten out into distinct grooves (Fig. 2, region A2). In section 3, we highlight additional observations to support our hypothesis for their genesis.
Depositional flutes and spicules
On the lee side of bed obstacles, a light-coloured crusty deposit is found that contrasts conspicuously with the substrate below (Fig. 2, region B). It is made up predominantly of crystalline calcite (CaCO3) precipitate and can be up to a few centimetres in thickness. Such deposit also outlines the up-glacier side of some solutional furrows, and where extensive, it covers the furrow pattern.
The calcite deposit is usually laminated. Its spatial distribution, and variations in structure, chemistry, impurity content and isotopic composition have been reported by many authors (Reference Hanshaw and HalletHanshaw and Hallet, 1978; Reference Lemmens, Lorrain and HarenLemmens and others, 1983; Reference Souchez and LemmensSouchez and Lemmens, 1985; Reference Sharp, Tison and FierensSharp and others, 1990; Reference Frisia and BorsatoFrisia and Borsato, 1994; Reference Hubbard and HubbardHubbard and Hubbard, 1998). Based on such variations, Reference Sharp, Tison and FierensSharp and others (1990) have identified two forms of the deposit: sparite (sparry carbonate deposit), which is relatively thin, macrocrystalline, and occurs in patches of several square centimetres (Fig. 3a); and micrite, which consists of microcrystals in clear laminations, and which covers larger areas and is darker than sparite (Fig. 3b). Flutes sub-parallel to the local ice-flow direction are ubiquitous on the deposit (Reference Ford, Fuller and DrakeFord and others, 1970; Reference HalletHallet, 1976, Reference Hallet1979; Reference Peterson and MoresbyPeterson and Moresby, 1981). They are found on both sparite and micrite, and have surface morphology ranging from sub-millimetre scale ridges- and-grooves and columnar spicules (Fig. 3a) to well-defined flutes (Fig. 3b). Ice-flow direction is clearly important in determining the orientation of the flutes, but abrasion by fine particles cannot adequately explain their formation, since they do not resemble striae. Evidence from electron micrographs of the deposit cross-sections supports this, showing internal laminae that undulate with the surface, with widespread truncation being absent (Reference HalletHallet, 1979; Reference Sharp, Tison and FierensSharp and others, 1990).
Common examples of the surface morphology of calcite deposits on bedrock: (a) spicules on a sparite deposit, Blackfoot Glacier, Montana; (b) flutes on a micrite deposit, Grinnell Glacier, Montana; (c) stalactite-like spicules protruding from a steep lee side (specimen from Castleguard Glacier, Canadian Rockies). The inferred direction of former ice flow is (a) right to left; (b) top to bottom; (c) approximately out-of-page.
Where the obstacle’s lee side is steep to vertical, conspicuous spicules may develop, resembling stalactites (Fig. 3c), but instead of being vertical, they appear to parallel the local ice-flow direction. Their growth seems to be fuelled partly by bed seepage; and along inferred seepage lines where they develop into rows, the spicule separation is extremely regular (≈1 mm in Fig. 3c). The spatial organization of both spicules and flutes, and their similar spacing, suggest they are morphological expressions of the same instability. We examine this idea in section 4.
22 Regelation over a soluble bed
Temperate glaciers can move past bed obstacles by a combination of regelation and viscous ice deformation. The bumps of interest to us are small and close enough in separation (roughly ≾ 0.1m) for regelation to be dominant, so we neglect viscous ice flow in most of the discussion here.
Regelation occurs by pressure melting of basal ice on the stoss side of obstacles and refreezing of water on their lee side. The classical theory by Reference NyeNye (1969) and Reference KambKamb (1970) details the conductive heat transfer (from lee side to stoss side) and stress concentrations (high pressure at stoss, low pressure at lee) that enable its operation. For simplicity, here we consider two-dimensional geometry, illustrated in Figure 4a. Melt-water flows in a thin film from the stoss side to the lee side, where it refreezes, releasing latent heat. This water will tend to dissolve the bed initially because of its low solute concentration, and in this process it becomes enriched in solute.
(a) Schematic diagram of the basal regelation process. Dashed line indicates the cross-sectional view of Figure 7. (b) Sinusoidal bed profile z(x) with amplitude a and wavenumber k (thus, wavelength λ = 2π/k) and the associated film water flux q(x), where x denotes distance, (c) Dimensionless plot of film solute concentration 
against distance X (= kx) for three values of α (0.1, 1, 10), as given by the analytic solution in Equation (4), based on v = 0, β = 1. Dashed lines mark the saturation level (the curves become inapplicable above 
) and are equivalent to the composite model solutions in (d) at the limit αp → ∞. (d) Dimensionless numerical solutions 
of the composite dissolution/precipitation model for v = 0.01, α = β = β
p
= 1 and two values of α
p
(1, 10). The case α = 1 in (c) is included as dashed line for comparison. (e) Dimensionless dissolution rates of the composite model for v = 0.01, α
p
= β = β
p
= 1 and three values of α (0.1, 1, 10). (f) Dimensionless precipitation rates of the composite model with the same parameter values as in (e).
To make this idea precise, let us refer to the dissolution of calcite. If the concentration of Ca2+ ions in the film is c(x), and the film thickness and water flux per unit width are h and q(x), respectively, a basic equation to describe the solute concentration profile at steady state is
where x denotes distance down-glacier. The first term on the righthand side has been found to be a reasonable first-order model for the rate of calcite dissolution in laboratory experiments (e.g. Reference Brown, Tranter and SharpBrown and others, 1996) and karst landform studies (Reference DreybrodtDreybrodt, 1990); usually β ≈ 1. K is a rate constant and cs denotes the saturation value of the solute concentrationFootnote 1 . The last term in Equation (1) describes solute diffusion in water, for which a typical value of the diffusivity D is 10−9 m2 s−1. The Nye–Kamb theory (Reference NyeNye, 1969; Reference KambKamb, 1970) predicts the film thickness h to be several microns. Reference HalletHallet (1979) gave an independent estimate, h ∼ 10 μm, based on size analysis of debris fragments within the carbonate deposit, consistent with the theoretical value.
The behaviour of the solution to Equation (1) may be examined by first neglecting the diffusion term. Such an approximation remains valid until concentration gradients become high. Where c < c s, the solute flux (magnitude of q × c) increases in the flow direction, i.e. dissolution always adds to the solute flux. Variation in c may therefore be deduced by dividing the solute flux by q. Along lee sides, q(x) is a decreasing function due to freezing, so c increases, and particularly we anticipate c becoming large where q diminishes (e.g. q = 0 at the midpoint of the lee side in Figure 4a and b would suggest that c becomes singular there). Of course, this does not happen without limit because of diffusive effects. But whether or not diffusion is neglected, Equation (1) breaks down at saturation, when c = c s. Beyond this, the dissolution model must be replaced by one that describes precipitation. We see below that in the adjoining region, c is somewhat above c s, constituting finite oversaturation; in that region, continued freezing ensures the precipitation of calcite.
It is instructive to use an exact result for q(x) obtained from Nye’s (Reference Nye1969, Reference Nye1973) theory in the solution of Equation (1). We take β = 1, and assume ice sliding at velocity U over the bed profile in Figure 4b, consisting of a single sinusoid z = a sin kx (a denotes amplitude and k the wavenumber). Thus, x = 0 denotes the midpoint of the stoss side (“stoss midpoint’’), and the distance between successive bumps is λ = 2π/k. As in the linearized Nye–Kamb theory, we take ak ≪ 1 such that the bump height is much less than its length and bed slopes are small. If sliding is due entirely to regelation, the melt rate is approximatelyFootnote 2 equal to U dz/dx, and since the water flux is given by integrating melt rate over distance (together with a density correction due to the phase change), we find
This is depicted also in Figure 4b using a different vertical axis for the water flux.
The unsaturated region (c < c s)
By using the result in Equation (2), Equation (1) may be rewritten in the form
in which we have defined a dimensionless dissolution rate, 
. The parameter v = ρ
w
hDk/ρ
i
aU is a dimensionless measure of diffusion. It is typically small and is inversely proportional to the flow Péclet number, Pe = qλ/hD, which is large (λ is the natural length scale). If v = 0, only one boundary condition for the differential equation is necessary: dc/dx = 0 at x = 0, owing to problem symmetry (and the initial solute flux is zero). Then, for β = 1, it is straightforward to solve Equation (3) by using the method of integrating factor, and we obtain
as the approximate solution for the unsaturated region. Here we use the dimensionless variable X (= kx) to represent distance. This concentration profile is qualitatively representative of other cases where β ≈ 1, and we show it in Figure 4c for several α values. For comparison, taking a nominal value for the rate constant K = 2.5 × 10−7 m s−1 (from Reference DreybrodtDreybrodt, 1990) would give α ≈ 1.5 for sliding at 10 m a−1 over bed bumps 5 mm in height, spaced 50 mm apart. In this case, v is indeed small, ≾ 10−2 (and Pe is large) for typical film thickness h ∼ 10 μm (Reference HalletHallet, 1979). Note that despite “fresh” water being created at the stoss midpoint, the solute concentration is non-zero there, c(0) = αc s/(1 + α), because it is a depth-averaged value.
Combining with saturated region
For c ≥ c
s, precipitation may be modelled by writing 
instead of the dissolution term in Equation (1), where Kp
and βp
are the corresponding rate constant and exponent. On non-dimensionalizing we derive an equation similar to (3), except that the first term on its righthand side is 
(the parameter α
p is the counterpart of α). This equation, applicable for c ≥ c
s, has to be solved together with Equation (3). A necessary condition is that the solute fluxes in both the dissolution and precipitation regions are equal at the point of saturation. This allows the position of the saturation point to be determined as part of the solution. Because both q and qc are smooth across the transition, at which c = c
s, this leads to the concentration gradient being continuous there.
We present a composite numerical solution in Figure 4d, taking into account diffusion, with v = 0.01, and α p = 1 or 10, β p = 1, for the case α = β = 1. We have used an iterative algorithm (Newton’s) and imposed boundary conditions based on symmetry, as well as the condition of solute flux continuity mentioned earlier. The asymptotic solution in Equation (4) (shown in Fig. 4c, reproduced as a dashed line in Fig. 4d) approximates well for c < c s, because v ≪ 1 and concentration gradients are not high. When c ≥ c s, a limiting behaviour is that the degree of oversaturation is minimized if precipitation is fast enough to accommodate solute flux changes due to freezing. This is illustrated by the difference between the solutions for α p = 1 and α p = 10. If α p ≫ 1, we have c/c s ≈ 1. Locally, the precipitation rate may then be taken as −c s × dq/dx, where dq/dx < 0.
The position where saturation is attained should depend on how much dissolution has occurred already. In Figure 4e and f, we show the dimensionless dissolution and precipitation rates for α = 0.1, 1, 10, when α p = β = β p = 1. As expected, the saturation point, marking the switch from dissolution to precipitation, varies with α. Putting aside the small effect of diffusion, this point lies down-glacier of the bump crest (π/2 < X < π) and never on stoss sides; this is because on stoss sides, melt dilution, while q is still increasing, offsets dissolution when c is sufficiently close to saturation. Such a prediction agrees well with the field observed pattern, in that stoss sides are normally free of precipitates. On the other hand, lee-side calcite deposits on most specimens extend back close to the bump crest, suggesting relatively high values of α and dissolution rates.
The carbonate system
Equation (1) and its counterpart drastically simplify the dissolution/precipitation process, which involves various species and a number of reactions. Calcite–water systems are usually studied by means of the equilibrium chemical equationFootnote 3
The protons may be derived from a number of sources, notably dissociation of dissolved CO2:
The reaction kinetics of pure calcite–water systems have been investigated by Reference Plummer, Wigley and ParkhurstPlummer and others (1978). For detailed modelling, two other factors need to be considered: (i) the geometry of the system (the regelation film), which determines the distribution of chemical species as a result of transport and reaction within it, and (ii) whether the system is “open” or “closed” with respect to CO2 and water, which determines the pCO2 of the solution (Reference Tranter, Brown, Raiswell, Sharp and GurnellTranter and others, 1993). Through its role in regulating the solution’s pH value via the reactions in Equation (6), pCO2 controls both the equilibrium concentrations (specifically, c s) and rate constants (K, K p), so that these may vary with distance along the film. On the stoss side of a bump, CO2 is continuously replenished by gas bubbles within melting ice, approximating an open system, so we may assume [CO2]water,stoss ≈ [CO2]ice and anticipate K and c s depending directly on the CO2 content of the source ice. For reaction in micron-thick films, such as in the current case, Reference Buhmann and DreybrodtBuhmann and Dreybrodt (1985) showed that K depends also on the film thickness.
For our purpose, there is another motivation for considering pCO2 variations, concerning the situation encountered along lee surfaces. Because of freezing there, gaseous CO2 is not derived from the ice, and the system may be treated as closed if we move with the water flow. pCO2 of a water packet is expected to decrease over a short distance downstream of the bump crest, due to dissolution (if α is not very large). But beyond the saturation point (see Fig. 4e and f), progressive freezing will increase pCO2, because dissolved CO2 is rejected by the freezing water as the reactions in Equation (6) are reversed and calcite precipitation takes place. With continued freezing, CO2 gas bubbles will nucleate in the film and they may be incorporated into the regelation ice. Although gas bubbles associated with the regelation process are not always present in basal ice, they have been observed before (Reference Kamb and LaChapelleKamb and LaChapelle, 1964). As we argue in section 3, their existence can explain the formation of solutional furrows.
23 Model for bed evolution
The dissolution/precipitation rates such as those derived from Figure 4e and f allow us to describe how the bed surface changes with time. First, notice the implication of Equation (2) that q(x) ∝ Uz(x) applies not only to a single sinusoid, but to an arbitrarily shaped periodic bed with small bed slopes, neglecting viscous deformation. This geometrical property holds for the closed system envisaged here — where there is no through-flow and all the film water derives from regelation — as long as z(x) averages to zero over one period 
. In other words, one can deduce the water flux directly from the bed profile (their shapes would be identical) if we choose the vertical datum at the level of the average bed height.
In this case, a compact model representing combined regelation and bed evolution is given, dimensionlessly, by the equations
(following Equation (3)), and
We have represented the bed by writing z(x, t) = af (X, t), f being its shape, now function of time t. 
is the normalized solute concentration c/c
s, and f
m denotes the mean of f over one period. The distance X is as defined before. In Equation (8), c
d and c
b are Ca concentrations within the calcite deposit and the bed substrate, respectively; we are referring to the Ca that resides in carbonate minerals.
Knowing f at any time, Equation (7) may be solved for 
; the bed profile then evolves according to Equation (8). f
m appears in the model for the reason given earlier of representing water flux correctly (f
m becomes non-zero in the event of net aggradation or degradation, which results when c
b ≠ c
d). In arriving at these dimensionless equations, we have chosen a time-scale based on dissolution, t
s = ρ
w
c
b/ρ
i
Ukc
s. Generally c
b ≫ c
s, so t
s is large, reflecting the fact that bed evolution is a very slow process. For example, taking c
b for a limestone with 30 wt% −Ca and c
s for a 10−3 M solution leads to c
b/c
s ∼ 104 . Together with the typical values U = 10 m a−1 and k = 2π/0.05 m−1 used previously, this gives a time-scale of 10 years.
The preceding model may seem rather elaborate, but we have formulated it not so much for simulation purposes as to illustrate two specific points. First, it brings out how the linearity in Reference NyeNye’s (1969) theory is lost in the coupling of water-transport and chemical processes, even in the case of constant model parameters and first-order kinetics (β, β
p = 1). As a result of solute advection (the product 
in Equation (7)) and the switch-over at saturation, a sinusoidal bed surface does not remain sinusoidal, and it is not possible to predict bed changes by considering separate Fourier components.
This leads us to the second point. We shall argue in section 4 that depositional spicules may be the manifestation of a spatial instability at the lee sides of bumps. It is tantalizing to explain the solutional furrows in a similar fashion, by demonstrating that small-wavelength bed perturbations on stoss sides are unstable. So far, such an approach (not presented here) has yielded no instability. For this purpose, Equations (7) and (8) constitute essentially a non-linear wave equation, of the form 
, from which it may be inferred that (solutional) bedforms propagate down-glacier, typically at velocity (kt
s)−1 ≈ Uc
s/c
b (≪ U). A furrow-forming mechanism based on instability cannot be ruled out, however, because our wave equation describes only two dimensions (an extension to three dimensions necessitates a calculation for the vector film water flux over the bed, which is not straightforward), and also we have not considered possible feedbacks in dissolution kinetics. On the other hand, we believe there is a much simpler explanation for the furrows, as discussed next.
3 A Hypothesis for Solutional Furrow Formation
The corrosive origin of solutional furrows is unequivocal, but their morphology has not been explained. Apart from the lack of an obvious instability, their location and orientation suggest that the localized dissolution recorded by them is unrelated to water flow in Nye channels or Röthlisberger channels. Neither can they be explained by incorporating chemical dissolution into the channel-forming instability examined by Reference WalderWalder (1982). At stoss sides, channels at the observed scale would close rapidly due to insufficient heat dissipation. Further, channels incising into ice and transversely aligned to the ice flow are advected down-glacier, and cannot support prolonged dissolution at the same place on the bed. Calcite deposits found on the up-glacier face of individual solutional furrows often develop the flutes described in section 2.1, and this supports Reference Walder and HalletHallet’s (1979) argument that the furrows form in intimate ice–bed contact.
We put forward a new interpretation here, hinging on the fate of the carbon dioxide exsolved at the lee side of bumps. We mentioned CO2 builds up in the water film during freezing because it is rejected by growing ice, at the same time as calcite precipitates. Experiments on the freezing of dilute solutions containing calcium and bicarbonate ions, reported recently by Reference Killawee, Fairchild, Tison, Janssens and LorrainKillawee and others (1998), show that bubbles that nucleate at the freezing front are incorporated into the ice as gas inclusions (this is, provided the ice is on top of the solution). Importantly, these inclusions were invariably found to be rich in CO2, with higher concentrations than that in the surrounding ice. Based on these observations, we envisage the following mechanism: once CO2-rich bubbles form and are incorporated into regelation ice on a lee side, they are advected down-glacier, and when they reach the stoss side of the next bump, they locally enhance dissolution of the bed surface by maintaining high pCO2 there, initiating furrows (Fig. 5). Film flow has the secondary effect of advecting chemically aggressive water along streamlines, leading to streaks of enhanced dissolution; this is manifested in the observed furrow pattern.
Cartoon showing how CO2-rich gas inclusions in the ice are advected down-glacier, leading to localized bed corrosion when they arrive at the next stoss surface (inset). Dashed line marks the elevation threshold for the resulting solutional furrows.
Several features, which are particularly evident on bedrock in front of Blackfoot Glacier, Rocky Mountains, Montana, support the proposed mechanism:
-
(i) On stoss sides, the domain occupied by solutional furrows most often terminates abruptly, rather than gradually, before reaching the bump crest (Fig. 2). Remarkably, the upper limit of calcite deposits lying immediately up-glacier, on extension along the (inferred) ice-flow direction, seems to coincide with this elevation threshold (dashed line in Fig. 5). This may be verified visually by sighting down-glacier along the bed surface, and suggests a link between the deposits and the furrows. Although the agreement is not perfect, the geometrical connection is compelling. We attribute varying degrees of agreement to viscous deformation, which would have some effect on local ice flow at such length scale (≾ 0.1 m), even though regelation may be the dominant sliding mechanism.
-
(ii) Small pits resembling pock-marks are found on stoss sides alongside solutional furrows. They extend transversely and may coalesce. We interpret them as furrows at an early stage of development, and attribute “scalloped” surfaces to their coalescence as a result of corrosion at neighbouring sites by CO2-rich bubbles. Small-scale or shallow transverse channels that may be classified as “incipient” solutional furrows have not been observed.
-
(iii) The solutional furrow width is comparable to the diameter of gas inclusions in the experiments of Reference Killawee, Fairchild, Tison, Janssens and LorrainKillawee and others (1998), which, regardless of the inclusion shape, is of the order of 1 mm or less.
In summary, the empirical evidence is consistent with solutional furrows reflecting recurrent impingement of isolated (CO2-rich) gas inclusions on stoss sides. Above a certain elevation threshold, solutional furrows are rare and dissolution occurs more uniformly due to the lack of inclusions advected from CaCO3 precipitation zones up-glacier. This mechanism accounts for the observed dissolution pits, furrow tortuosity and intersection, and overlapping of deposits with furrows. Trains of spherical and cylindrical bubbles have been identified in the basal regelation ice of Blue Glacier, Washington, U.S.A., by Reference Kamb and LaChapelleKamb and LaChapelle (1964, fig. 4), although their chemical composition was not reported. Well-defined pits in the solutional furrow morphology may be the result of such bubble trains, originating from repeated nucleation at preferred sites on the deposit (as described by Reference Killawee, Fairchild, Tison, Janssens and LorrainKillawee and others, 1998). As yet, our field observations do not indicate that these sites are necessarily located at the tips of the precipitate spicules. We point out in addition that Figure 5 is only schematic. Further investigation is required to establish what controls the spatial distribution of bubbles in the ice, and whether or not a repeated delivery of CO2 from bubble trains is needed to generate the observed furrow depth.
The observed pattern of solutional furrows may be used to test our hypothesis. Three factors control the course of CO2 release from a gas inclusion that is brought in contact with the bed: (i) regelation water flow advects dissolved CO2; (ii) ice motion carries the inclusion some distance along the bed surface if the bed slope is low, and then the CO2 source itself is “mobile” and its trajectory may leave behind a dissolution groove on the bed; (iii) CO2 diffuses. In Figure 6, we plot the pressure and water-flow fields (Fig. 6b) over a hypothetical bed consisting of a two-dimensional array of small bumps, z = a sin kx sin ωy (Fig. 6a), where ak, aω ≪ 1. This model is only a crude representation of the bed topography encountered in reality, but is easy to analyze. The water-flow direction indicators have been calculated using Reference NyeNye’s (1969) vectorial formulation — via −∇p w, where p w is film water pressure, ∝ cos kx sin ωy —assuming only water transfer associated with melting and refreezing (closed system). With regard to ice motion, we bear in mind that it will diverge slightly around the bumps if viscous flow is taken into account, but this effect will not be significant given the (assumed) small separation between bumps, and ice velocities will be close to (u,v) = (U, 0), or exactly so in the case of pure regelation.
(a) Shaded relief and contours of a two-dimensional hypothetical bed topography with the form z ∝ sin kx sin ωy and of low amplitude; light (dark) grey denotes high (low) elevations, (b) Water-pressure field p w resulting from pure regelation over the topography in (a); light (dark) grey denotes high- (low-)pressure areas. Arrows represent the vector gradient −∇p w and point in the local direction of regelation water flow. Regions A1 and A2 are discussed in the text.
(a) Transverse furrows
Near the stoss midpoint, relatively high bed slope ensures that inclusions melt out locally and thus the CO2 sources should have low mobility. The water-flow field is then the dominant factor, and we see the way water streamlines emanate from the stoss midpoint and diverge ahead of the bump (Fig. 6b, region A1) reflected in the furrow pattern of region A1 in Figure 2. However, the streamlines traverse the bump crest whereas the observed furrows do not. This discrepancy may be due to the fact that the streamlines are modified as the furrows develop to have large amplitudes, not accounted for by our current calculation (moreover, our model topography is an idealization).
(b) Oblique furrows
Along the sides of bumps (Fig. 2, region A2), the furrows are again well represented by the local water-flow direction (Fig. 6b, region A2), sub-parallel to the ice flow. If, instead, this is in reality a region of dividing water streamlines and low water fluxes (such as towards the bottom half of box A2 in Fig. 6b), the bed slope is low, so the effect of ice motion will dominate and cause the occurrence of groove-like furrows parallel to it, on the arrival of inclusions. In either case, furrows with the observed orientation are predicted.
(c) Péclet number
Localized CO2 release will not cause the well-defined streaks of dissolution as observed if lateral diffusion of dissolved CO2 within the film is fast compared to advection by water flow. The relevant parameter to consider is a Péclet number, Pe = qδ/hD, the ratio of advective and diffusive effects, and diffusion is negligible if Pe ≫ 1. δ denotes a length scale, and we use the notations defined previously, taking D ∼ 10−9 m2 s−1. In section 2.2 we put δ = λ, the bump separation. For the current calculation, we are interested in how a CO2-rich packet of solution evolves in plan view over the bed surface, after having been released over an area comparable to the inclusion cross-section. Therefore, the appropriate choice for δ is the diameter of the inclusion (∼ 0.1mm, from Reference Killawee, Fairchild, Tison, Janssens and LorrainKillawee and others (1998)). With the example given in section 2, where U = 10 m a−1 and a = 5 mm, the typical water flux is q ≈ 2 × 10 m2 s−1 from Equation (2). A top-end estimate of h = 100 μm then leads to Pe ≈ 2. This is a conservative result, because Pe ∝ h −1 and the film thickens to 100 μm only occasionally, h ∼ 1−10 μm being common (Reference HalletHallet, 1979). We therefore anticipate Pe ∼ 10−102 (≫ 1) in general, consistent with the observed manifestation of furrows.
4 Depositional Instability
In this section we show that the coupled processes at the lee side are spatially unstable. Consider a cut made there in the bed-normal direction as indicated by the dashed line in Figure 4a. In the transverse plane containing the cut (Fig. 7), the ice separates from the bed at the same time as regelative water flows into the gap to freeze — the freezing rate is such that the regelation film thickness remains constant. We include calcite precipitation in the following consideration, but the kinematic balance between deposition and separation is unaltered.
Definition diagram for our mathematical model (section 4) applied to the region near the midpoint on a lee side. Box-arrows indicate separation velocities of the ice and the bedrock (including calcite precipitate) relative to the interface (near z = 0) where freezing and precipitation occur.
Let us assume freezing at the planar ice–bed interface occurs at rate V (water volume per unit time per unit area), equal to −dq/dx using the notations of section 2. The latent heat released is conducted to the stoss sides through ice and rock. The heat flow is therefore three-dimensional, but our model focuses on what happens close to the centre of the lee side, and this enables the use of a local approximation: when length scales much smaller than the bump are considered, the heat flux far from the interface may practically be taken to occur in a direction normal to the interface (even though the lines of heat flow ultimately curve up- and down-glacier). This assumption greatly simplifies the boundary conditions of the temperature calculation in section 4.1 below.
Referring to Figure 7, take x and y, respectively, to be the down-glacier and transverse coordinates, and define z to be the bed-normal coordinate relative to the moving interface. At steady motion, the vertical ice velocity is V/r, where r = ρ i/ρ w. As calcite is actively precipitating from a saturated film, the bed has a relative velocity −Vc s Strictly, c s is the volume of precipitate that results from freezing a unit volume of saturated film water, corrected for density change, solute inclusion in the ice, and oversaturation. Typically c s ≪ 1, so downward motion of the bed is slow.
The question is whether or not the precipitate surface remains stable when infinitesimal perturbations are imposed on it. The perturbations will grow in an unstable manner if, through positive feedback, they promote an increase in the rate of freezing — thus, proportionally, of calcite precipitation — at the crests of the bed surface, and a decrease at the troughs. This picture is motivated by what happens in the Stefan problem for a binary alloy, in which a freezing front advances into a liquid melt, and where instability leads to the formation of dendrites (Reference Mullins and SekerkaMullins and Sekerka, 1964). Our mechanism turns out to be unique, however, as it involves pressure and flow effects in the regelation film and within the ice, in addition to freezing.
Variations in freezing rate along the interface originate from a change in the heat loss from it, caused by perturbations. These perturbations are inherent in the irregularity of any real bed surface, but those that exist in the film pressure and film thickness are also relevant and have to be considered (sections 4.1 and 4.2). Variations in film pressure and thickness drive ice flow out of, and water into, areas where the freezing rate is enhanced; their magnitudes are tied to the freezing rate through flow processes. On the other hand, we mentioned that film pressure and thickness play a part in controlling the freezing rate. Consequently, our deduction below of how the bed surface evolves through freezing/precipitation involves the simultaneous solution of several problems (section 4.3). As we shall see, the perturbative effects of film pressure and thickness are respectively stabilizing and destabilizing, and the outcome of their competition is that flutes develop at the observed length scale (section 4.4). To examine linear stability of the system, various ingredients are first derived.
41 The temperature problem
As assumed in the Nye–Kamb theory, the regelation film is thin, so the temperature difference across it may be neglected. We shall allow film thickness, h, to vary, however. We are concerned with perturbation wavelengths that are large compared to h (which is of the order of microns), because the observed flute spacing (of the order of millimetres) is much greater than this.
When the ice–bed interface is perturbed from its position at z = 0 to, say, z = ∊f(x, y), where f has zero mean and ∊ ≪ 1, the conductive heat losses into the ice and bed are modified. Specifically, their sum will depart from the average (V times latent heat of freezing, L), resulting in non-uniform freezing and precipitation along the interface. Variations in water pressure p w in the film will contribute to this effect through altering the melting point.
To establish the quantities involved, notice that within the ice and bed the temperature θ satisfies Laplace’s equation,
Based on the local approximation introduced at the beginning of this section, we prescribe heat losses far from the interface in the z direction only. Hence, appropriate boundary conditions for Equation (9) are (from Fourier’s law)
in which we use subscripts i and b to refer to ice and bed, respectively, k i,b to denote the thermal conductivities (k b assumed equal for rock and precipitate), and q i,b to denote the (constant) heat-loss values. Additionally, both materials share the same interfacial temperature, equal to the pressure-melting point. If θ and p w are measured relative to their unperturbed (and average) values, we can write
at the boundary, where C is the Clausius–Clapeyron constant. p w is derived later by considering processes in the regelation film.
When f ≡ 0 (the unperturbed situation), θ is a function of z only and is linear in each material. This solution of Equation (9) is shown in Figure 8a. We have deliberately chosen
to ensure that the ice and bed temperatures are symmetrical about z = 0, as in the classical basal sliding theory (see Reference NyeNye, 1969). The temperature solution is θ(z) = −VL|z|/(k i + k b). (Our far-field temperatures then match smoothly with the large-scale three-dimensional temperature distribution at the lee side, observed at the length scale of the bump.) In this case, the interfacial heat losses are the same as the heat fluxes within the ice and the bed, q i and q b (and q i + q b = VL). But they are not equal to each other when k i ≠ k b, despite the temperature symmetry. As we shall see, this is a necessary condition for the lateral instability.
The vertical temperature profile θ(z) (right column) at various horizontal positions calculated for the following three situations of the ice–bed interface (left column): (a) an unperturbed interface; (b) a sinusoidal perturbation in regelation film pressure only; (c) a sinusoidal perturbation in interface position only, with the regelation film thickness and pressure unperturbed. The linear profiles in (a) are shown as dashed lines in (b) and (c) for comparison.
We now calculate the effect of perturbing the interface to first-order accuracy (in ∊). Anticipating that film-pressure fluctuation is a first-order effect, let
By defining T(x, y, z) = θ(x, y, z) − q b z/k b, the problem for bed temperature reduces to
(b.c.’s are boundary conditions), and this motivates an asymptotic expansion also for the temperature,
Following Reference NyeNye (1969), Fourier transform is employed to solve the first-order problem. We do this in two dimensions and use the notation 
to denote the transform of T
i, i.e.
taking k and ω as wavenumbers for the x and y directions, respectively. The solution is found to be
and this implies
An identical method applied to the ice (where we define T = θ + q i z/k i) leads to
We retain these results in the transformed notation for use in the stability analysis below.
The physical meaning of Equations (17–20) is illustrated in Figure 8b and c for perturbations at a fixed value of k or ω. Each term within the square brackets indicates the sign and magnitude of the temperature or temperature-gradient response for given perturbations in film pressure (p
w1) and interface position (f). Far from the interface, temperatures are unperturbed 
because the effects of the “ups and downs” at the interface cancel out.
An elevation in p w depresses the melting point and reduces the temperature lapse rate into both materials, thereby reducing the heat loss and the rates of freezing and precipitation at the interface locally (and vice versa) (see Fig. 8b). The sign of this effect is independent of the magnitudes of k i and k b.
On the other hand, where the interface has been lifted (above z = 0), the temperature lapse rate into the ice is increased and the lapse rate into the bed reduced — by the same amount, since in Equations (18) and (20) q b /k b ≡ q i/k i (see Fig. 8c). The temperature field is distorted in this way because the (relatively) warm interface is now situated where it was originally colder. Locally, the heat loss to one material is enhanced while the heat loss to the other is reduced, and the net effect depends on which material is more conductive: where f < 0, we have enhanced net heat loss (more freezing and precipitation) if k b > k i, but reduced net heat loss (less freezing and precipitation) if k i > k b. The opposite result is encountered where f > 0. To first-order accuracy in ∊, this effect is absent if k b = k i.
42 Interfacial processes
If a perturbed interface forces a perturbation in the rate of precipitation, modifying the interface over time, a feedback may fuel an instability so that a fluted bed surface develops. Apart from the effect of pressure, yet to be discussed, the lower and upper boundaries of the film are subject to different processes, despite the film being thin. It is therefore necessary to perturb film thickness h as well, and this removes the constraint that the ice “sees” the same interface as the bed, z = ∊f. Formally, we let h = h
0 + ∊h
1(x, y) (+higher-order terms), where h
0 denotes the uniform-state film thickness. From now on, ∊f and ∊ (f+h
1) thus refer respectively to perturbations of the lower (bed) and the upper (ice) interface of the film. Accordingly, we replace 
by 
in results (19) and (20).
Calcite deposition
The total heat loss from the film is given by
In conjunction with Equations (18) and (20), the freezing rate of water (as a volume rate per unit area of the interface), φ, may be calculated from this explicitly, leading to the expression, in transformed notation,
φ − V is the freezing-rate perturbationFootnote 4 (a velocity difference), precisely what we need to evaluate spatial variation in precipitation rate. Since the rate ratio of CaCO3 precipitation to film water freezing is c s, as assumed earlier, an evolution equation for the bed surface is ∂(∊f)/∂t = c s φ − c s V, or equivalently
where 
is taken from Equation (22). Notice now we treat f (and other variables) as a function of time. If precipitation occurs faster (slower) at the crests (troughs) of the bed profile, the perturbations will grow. The goal is to derive their growth rate by relating 
to 
, so we need additional relations for 
and 
to close the model.
The freezing front
Our main assumption is that there is intimate ice–bed contact during flute development on the deposit. More generally, the ice–water interface evolves in time also, and thus there is in principle a separate evolution equation for f + h 1, but we assume it relaxes rapidly and is practically at equilibrium given the very long time-scale of calcite precipitation (caused by c s ≪ 1 in Equation (23)). This allows us to neglect the derivative ∊∂(f + h 1)/∂t which otherwise would appear on the lefthand side of Equation (25) below.
Consequently, our description here is identical to that of Nye and Kamb. Regelation is assumed to be quasi-steady, and non-uniform freezing (and hence, ice production) is everywhere accommodated by the ice sliding motion, directed down-glacier approximately at velocity U, superimposed onto which there is ice deformation in response to spatial pressure variations. If the normal ice velocity at the interface is
where w n denotes the perturbation caused by p w1, then to first order, the balance that is satisfied at the freezing front is
For obtaining w n, Reference NyeNye (1969) has already solved the Stokes flow problem for ice of constant viscosity η i. Extending his analysis to three dimensions yields the relation
As one may expect, this indicates that ice is driven upwards faster at high-pressure areas, and vice versa. Substitution of Equation (26) into the Fourier transform of Equation (25) leads to
Film water flow
The final ingredient comes from water mass conservation because spatial variations in the freezing rate require a lateral water transport in the film; more water has to flow into regions where freezing is more intense. Specifically, 
obtained by solving Equation (27) will drive water flow, but the resulting flux convergence will generally differ from what is needed for the freezing. The deficit is compensated by film-thickness changes, because h controls the water flux for a given pressure gradient.
If the water viscosity is η w and laminar flow is assumed, the water flux is given by
(e.g. see Reference WeertmanWeertman, 1972), and the vector components here satisfy the conservation equation
The factor 1 + c s refers to water + solute. (Again, a time derivative ∂h/∂t on the lefthand side of Equation (29) has been neglected under the quasi-steady assumption.) Equation (28) is a Poiseuille flow approximation which breaks down at perturbation wavelengths comparable to or smaller than h 0, but is adequate for our purpose.
To ascertain the effect of perturbations, Equations (28) and (29) need to be linearized about their basic state. Near the lee midpoint, this is defined by q⊥ = 0, and
which describes water flow into the plane of Figure 7 in the unperturbed situation. Expanding Equations (28) and (29) and using Equation (30) leads to
at first-order accuracy. Since q 0 vanishes as we approach the lee midpoint, we may neglect the last term in this equation, in accordance with the local approximation in section 4.1. In transform notation, Equation (31) then becomes
This provides the value of 
whereby a (perturbation) water source from up- and down-glacier balances both the lateral divergence of water flux and water loss to freezing.
43 Is there an instability?
The dispersion relation
The model for analyzing the stability of the interface is based on Equation (23), supplemented by Equations (22), (27) and (32). The processes being considered are coupled, but we can use the last three equations to solve for 
, 
and 
simultaneously, and then express the freezing-rate perturbation 
in terms of 
alone. To facilitate this procedure, let us define
κ and λ 1… λ 5 being positive, dependent on model constants only (κ is a thermal conductivity ratio). After performing the algebra, we obtain
in which the transfer function is
To keep this expression concise, we have used the shorthand 
to denote the composite wavenumber, and also λijk
… to denote the product λiλjλk
…, etc.
The evolution equation (23) now takes the form 
. If we use the trial solution 
, the complex growth rate, σ, is given by the dispersion relation
Its real part indicates whether the perturbations are stable (Re σ, Re G < 0) or unstable (Re σ, Re G > 0). The dependence on k and ω enables identification of the wavenumber range over which instability is predicted. In particular, one would expect the most unstable wavenumber (at which Re σ is most positive) to correspond to the observed depositional flute separation, assuming that it is identical to the spacing of the fastest-growing infinitesimal flutes.
Lateral instability
To examine G, consider first the transverse direction, putting k = 0 such that ξ ≡ |ω| and G = G(ω) only. G is then real. The second term in the denominator of G, containing λ
5, is negligible compared to the first because η
i is large (below). It follows that the denominator has a zero crossing near 
, due to the term ξ|(1 + c
s)λ
1234 − ξ
2| (remember c
s ≪ 1). Because the numerator is non-zero for κ ≠ 1, G is singular at this wavenumber. This underlies the instability we have been seeking.
In anticipation of this behaviour, our following discussions are facilitated by non-dimensionalizing Equation (35), letting k* = k/ω c , ω* = ω/ω c (and thus ξ* = ξ/ω c), where we define the wavenumber scale ω c to be
Equation (35) then becomes
in which the dimensionless parameters μ 1, μ 2 and μ 3 are given by
These parameters measure the importance of viscous deformation, sliding geometry and freezing rate, respectively. (μ 2 is velocity-independent because V/U is given by the bed slope.) G has the unit of a growth rate, s−1, via λ 3.
Remembering that the growth rate of perturbations is proportional to G, we estimate realistic values for the parameters appearing in Equation (38). Calcareous bedrock and calcite are typically more conductive than ice, k b ≈ 4 W m−1 K−1 compared to k i = 2.1 W m−1 K−1 (e.g. Reference ClarkClark, 1966), so we take κ = 2. With h 0 = 1 μm, η i = 3 × 1012 Pa s and the other model constants as listed, we obtain ω c = 1.8 × 104 m−1, and hence μ 1 ≈ 4 × 10−7 and μ 3 ≈ 500. The effect of μ 2 is considered later. The fact that μ 1 and μ 1 /μ 3 ≪ 1 in Equation (38) implies that σ ∝ −(κ − 1)/[ω ⋆(1 + c s − ω ⋆2)] when k ⋆ = 0 and ω ⋆ = O(1), confirming our earlier indication that a singular growth rate occurs close to ω ⋆ = 1 (or ω ≈ ω c). The corresponding wavelength is 2π/ω c ≈ 0.4 mm. Notice the denominator of Equation (38) has no other zeros, and that σ ∝ −ω ⋆ as ω ⋆ → 0, σ ∝ ω ⋆−3 as ω ⋆ → ∞.Footnote 5 From now on, we refer to the ω-value at the singularity as the “critical wavenumber”.
Figure 9 illustrates the wavenumber-dependent factor of Equation (38), μ 3 G/λ 3 μ 1 (proportional to Re σ, the real part of σ), for k ⋆ = 0. Re σ becomes large and positive as we approach the critical wavenumber (at ω ⋆ ≈ 1) from above. Perturbations with this critical wavenumber will therefore be the fastest-growing, responsible for the observed depositional flutes. Our linear theory predicts that the flute spacing should be of the order of a millimetre if h 0 ∼ 1 μm, and this agrees well with observations. (Using h 0 = 10 μm gives λ = 40 mm, which is still reasonable.) While the growth rate does become singular in this model, in reality it is limited by time evolution of the freezing front and film water flow that occurs in a fast time-scale, hitherto assumed to be “instantaneous” under the quasi-steady assumption. Interestingly, the instability vanishes if μ 1 ≡ 0 (corresponding to the rigid ice limit η i → ∞), but persists as long as the ice has a finite viscosity.
Dispersion diagram calculated for perturbations with a transverse orientation (for which k⋆ = 0), taking h 0 = 1 μm, κ = 2 and other model constants as described in the text. Horizontal axis is the dimensionless wavenumber ω⋆. The wavenumber scale is ωc = 1.8 × 104 m−1. The growth rate of perturbations (Re σ) is proportional to the vertical axis. Inset shows details near ω⋆ = 0, and bubbles indicate the asymptotic behaviour at large and small wavenumbers.
Spicules vs pins
A singular growth rate in the wavenumber domain seems to explain why depositional flutes should be well defined spatially. Equation (38) reveals further properties consistent with observations, regarding stability in two dimensions. For general values of k ⋆ and ω ⋆, the imaginary part in both the numerator and denominator of G is non-zero unless μ 2 ≡ 0, where μ 2 (∝U/V) is the dimensionless ratio of the sliding velocity to the average freezing rate. This has two implications:
-
(i) When U ≠ 0 (and thus μ 2 ≠ 0), the growth rate Re σ becomes singular only for perturbations with a purely transverse orientation, for which k ⋆ = 0. As we deviate from this orientation, the singularity is removed and the amplitude of Re σ is greatly reduced. In other words, the instability is most intense in the transverse direction, but becomes less so in directions that contain a component of the basal sliding velocity U. We believe this is responsible for the strong preferred alignment of flutes and spicules down-glacier. Such directionality is anticipated whenever μ 2 is non-negligible, which is typical for a low-amplitude wavy bed. For instance, taking the earlier example where a = 5 mm (bump height), λ = 50 mm (bump separation), and using the corresponding maximum (lee) slope to deduce U/V ≈ λ/2πa, leads to μ 2 ∼ 5.
-
(ii) μ 2 becomes small when U/V ≪ 1. A special case is encountered where the lee side is vertical and the relative ice motion is normal to the interface (in the z direction only), and U = μ 2 = 0. Then the instability is radial symmetric and our model predicts the formation of cylindrical or pin-shaped spicules, such as the ones in Figure 3c. The preferred wavelength is identical to that for the flutings because G(ξ ⋆)|μ 2 = 0 and G(ω ⋆)| k ⋆=0 have the same form.
44 Physical basis
The issue of stability rests on the phase relationship between the bed surface perturbation f and the freezing-/precipitation-rate perturbation that results from it, ∝(φ − V)/∊. To facilitate the discussion here, we denote (φ − V)/∊ by ψ If 
, the system is stable because 
(from Equation (23)); if 
, the system is unstable because 
. How the phase relationship is determined and varies with wavenumber is complex. To help interpret Figure 9 we provide the following descriptive mechanism, assuming as before that the bed is thermally more conductive than the ice (κ > 1).
Bed perturbations f induce non-uniform freezing by their tendency to increase/decrease the interfacial heat loss relative to its average value at the troughs/crests (section 4.1). If no other perturbations exist, this mechanism is stabilizing because calcite precipitates at a higher rate at the troughs than at the crests; this happens at the long-wave-length limit. However, non-uniform freezing (ψ ≠ 0) is accompanied by film-pressure and -thickness perturbations (p
w1, h
1 ≠ 0), and these control the freezing rate also (Equation (22)). As the wavelength is reduced, the effect of the pressure and thickness perturbations becomes progressively more significant, to the extent that there is in fact a sign change, and transition to instability. The crux lies in understanding the (fast) feedback whereby p
w1 and h
1 that accompany ψ also control ψ, and how this features in the (much slower) feedback between the bed surface 
and its evolution 
.
The “fast” feedback
The transverse or pin instability with k or U = 0 is probably simplest to explain. First, how are p
w1 and h
1 regulated? Since sliding is absent, ice deformation has to accommodate for regions of enhanced freezing by an upward motion, driven by high pressure; thus, Equation (27) indicates 
, where λ
5 is defined in Equation (33). Since the film-pressure perturbation is then in phase with the freezing-rate perturbation 
, it drives water away from regions where water is removed at higher rates by freezing. In equilibrium, both losses are balanced by means of a thickened film, h
1 > 0, where ψ > 0 (and vice versa), because film thickening/thinning leads to greater/less water-flux convergence from up- and down-glacier. This balance is described by Equation (32), which, combined with what we have deduced so far, has the form 
(we neglect c
s). Thus, except at very low wavenumbers, “stiff” ice (1/λ
5 being large) leads to significant pressure variations, and this equation shows that film-thickness adjustments compensate mostly for laterally driven flux deficits, with 
.
We have established that the film is thicker and at higher pressure where freezing is enhanced (and vice versa): the perturbations h 1 and p w1 are in phase with ψ. Crucially, their feedback effects on ψ are opposite, as is apparent from Equation (22), rewritten here as
The reason high film pressures reduce heat loss via melting-point depression has been discussed in relation to Figure 8b, and this is described by the negative sign in front of the pressure term in Equation (40). Film thickening promotes high freezing rate by increasing heat loss up into the ice (and vice versa), and explains the positive sign in front of 
. The interplay of these two processes at different spatial frequencies controls the overall sign of the feedback.
To see the consequence, we write the feedback terms 
explicitly as
If ξ → 0, the feedback vanishes because non-uniform freezing over long distances induces vanishing pressure fluctuations, and Equation (40) implies 
. This is the stable longwave limit discussed before, at which the precipitation feedback involves perturbations in the bed surface only. At the length scales of interest, however, ξ is sufficiently large for the feedback terms to be much greater than
itself because, over short distances, (i) large film-pressure variations are induced to drive ice flow, and (ii) large film-thickness adjustments are required to compensate for the large water flux being driven sideways. The pressure/ thickness feedback on controlling heat loss becomes comparable with that due to 
. In this case, 
adjusts (in a fast time-scale relative to the time-scale of bed evolution) so that the dominant balance in Equation (40) exists between those terms appearing on its righthand side only. Whether 
or 
then depends on the overall sign of the feedback, governed by (−λ
4 + ξ
2/λ
123) in Equation (41), which is wavenumber-dependent.
Bed surface instability
At wavenumbers above (approximately) 
, the effect of film-thickness perturbations overrides that of pressure melting and the feedback is positive (−λ
4
+ ξ
2/λ
123 > 0). Where freezing rates are elevated (ψ > 0), the film thickening that develops (to compensate for pressure-driven water flow) tends to increase heat loss. But in order that the net heat-loss perturbation be consistent with the sign and magnitude of ψ, this region can only be situated at the crests, where f > 0 induces a heat-loss reduction. Denoting the change in heat loss by “CHL’’, we can summarize the heat-loss balance by writing
A similar argument shows why ψ < 0 occurs where f < 0, and thus 
. The positive phase relationship between ψ and f is established in the fast time-scale and gives rise to instability: precipitation rates are higher at the crests, lower at the troughs. Dominance of the film-thickness effect is also responsible for the shortwave-limit (ξ → ∞) behaviour indicated in Figure 9.
When ξ ≾ ω c, pressure melting dominates the feedback (now negative) and this is stabilizing. Where ψ > 0, the (relatively) high film pressure that develops (to drive ice flow upward) overrides the film-thickness effect, and the overall feedback reduces heat loss. This can only happen at the troughs where f < 0 induces heat-loss enhancement. The equivalent of Equation (42) is
In this case, , and the system is stable.
Even though the terms on the righthand side of Equation (40) are dominant, a small correction arises from the left-hand side. Transition to instability therefore occurs not at ω c, but a value very close to it (given by one of the roots of rλ 5 = ξ 2(−λ 4 + ξ 2/λ 123)), at which the feedback equalizes the freezing-rate variation that causes it. In the neighbourhood of this critical wavenumber, the response in freezing rate becomes extremely sensitive to bed surface perturbation, leading to highly unstable growth above, and decay (high stability) below, which is why the growth-rate singularity appears in Figure 9.
In summary, our detailed model of the physics at the precipitating interface indicates that a spatial instability arises, and in particular we expect the interface to become irregular in a way that is consistent with both the spacing and orientation of the observed depositional flutes and spicules. A possible extension of this is given next.
45 Scaling laws?
In the dispersion relation derived in section 4.3 (Equation (38)), because 
and 
remain small as long as h
0 ≾ 20 μm (the upper bound inferred by Reference HalletHallet (1979)), it is reasonable to take the critical wavenumber as effectively given by ω
c. For film thickness within this bound, therefore, we can write
for the wavelength of the instability (or 
, λ
c in mm, h
0 in μm). As noted before, there is good order-of-magnitude agreement of the λ
c predicted by this theory with observation.
Classical basal sliding theory provides an independent estimate of h 0. Taking the case of pure regelation, we have
(Reference NyeNye, 1973), in which λ b denotes the separation between bumps. Combining with the earlier result leads to
Both Equations (44) and (46) take the form of scaling laws, and in this respect they motivate a search for field data on flute spacing and the associated value of h 0 or λ b. A fundamental problem with testing the latter law, however, arises from the difficulty of deciding which observed bump separation should be assigned to which particular group of flutings, since, on a real glacier bed, bumps with different length scales are often nested within one another (e.g. Fig. 1). Direct application of the first law provides another way of testing our theory, since h 0 may be constrained by analyzing the rock fragment size found locally within the calcite deposit. Both approaches are pursued in our ongoing research effort.
5 Conclusions
We have endeavoured to seek out the mechanisms responsible for the extensive patterning found on deglaciated limestone bedrock surfaces which comprises solutional furrows and fluted carbonate deposits. That these features might have resulted from subglacial regelation over bed undulations has long been recognized, but up to now their intricate surface morphology has received much less attention, not least the notion of using it to infer past behaviour.
In this paper, we provide the first quantitative analysis that addresses this issue, based on a coupled theory of film chemical processes and regelation physics. Apart from initiating a mathematical framework that encapsulates the early ideas of Reference HalletHallet (1979) (section 2), we offer explanations for (i) why solutional furrows are collectively organized into arcuate patterns, with characteristic spacing, and why they should appear at all (section 3), and (ii) why a fluted surface texture is ubiquitous on the calcite deposits, and what determines the length scale of those flutes (section 4).
The understanding of (i) hinges on the role of CO2 exsolution in freezing and calcite precipitation at lee sides of bed bumps. We envisage that CO2-rich inclusions carried down-glacier by ice flow enhance bed surface corrosion on arrival at the next stoss faces, leading to solutional furrows. The surface flutings in (ii) are inherently the manifestation of a spatial instability at the depositing interface at lee side. We elucidate the coupled physical processes that give rise to this instability. We also put forward scaling relations between the observed flute spacing, regelation film thickness and the separation between bumps (section 4.5). Tentatively, we can combine one of the scaling relations with general results from basal sliding theory (instead of using Equation (45), which is a special case), in which the film thickness is expected to depend on sliding velocity U. The replacement of Equation (46) would then connect the observed flute spacing to U, and this raises the possibility of corroborating other estimates of former ice velocity.
Our modelling provides precise insight into the behaviour of the water film that lubricates the bed of a temperate glacier, and highlights the complex interactions leading to some of the striking glacier-bed features shaped by subglacial chemical processes. Given our analysis refers to the initial development of these features and explores mainly linear interactions, the problem of predicting finite-amplitude bed features awaits further investigation. However, it raises some interesting questions already. For instance, what is the effect of the through-flow of subglacial water? Does CO2 recycling between ice and regelation film water determine the spatial distribution of these bed features under a glacier? What can we infer from the laminations and isotopic compositions within the deposits? And what do observed differences between glaciers tell us? Given the rather complex coupling of mechanisms proposed here, it seems that an integrated approach with strong guidance from more detailed observations and experiments will be the key in future investigations.
Acknowledgements
F. Ng acknowledges the generous support of a Royal Society-Fulbright Foundation Science Fellowship to visit B. Hallet at the Quaternary Research Center, University of Washington. The authors thank R. C. A. Hindmarsh, B. P. Hubbard and scientific editor J. S. Walder for their careful and constructive reviews of the manuscript.
(with V being constant) introduces a delta function of wave-number which complicates our notation.
