4.1. Physical description

The physical system in which frost heaving occurs is shown in Figure 1. Top-down freezing occurs over a time-scale of weeks to months. The upper, frozen region is separated from the lower, unfrozen region by the frozen fringe, a partially frozen transition region containing ice, liquid water and soil. The thickness of the frozen fringe in seasonally frozen ground is on the order of 1–10 mm. The frozen region is characterized by regular bands of pure ice called ice lenses. It is the existence of ice lenses that gives rise to surface heave greater than the 9% associated with the volume change of water upon freezing (Reference WashburnWashburn, 1980). In typical field situations where frost heave occurs, permafrost underlies the unfrozen region (Reference Williams and SmithWilliams and Smith, 1989).

Fig. 1.

Diagram of saturated soil undergoing top-down freezing. The region between z_l and z_f is the frozen fringe. Ice lenses occur in the frozen region, giving rise to frost heave.

Because the frozen fringe is relatively thin, it can be reduced mathematically to a plane across which jump boundary conditions are applied (Reference Fowler and KrantzFowler and Krantz, 1994). This plane is located at z _f in Figure 1. Freezing of saturated soils is fairly slow, as can be inferred from the large Stefan number associated with this process. Because the time-scale for freezing is much longer than that for heat conduction, the temperature distributions in the unfrozen and frozen regions are quasi-steady state. The thermal gradients in the unfrozen and frozen region are obtained by solving for the temperature distributions in each region.

Because ice lenses occur in a discrete manner, the frost-heave process is inherently a discontinuous process. However, Reference Fowler and KrantzFowler and Krantz (1994) demonstrate using scaling arguments that the time-scales for equilibration of the pressure and temperature profiles within the frozen fringe are much shorter than the time-scale for the freezing process. Thus, the entire frost-heave process can be modeled as semi-continuous where some parameters are evaluated at the base of the warmest (lowest) ice lens, z _l.

For one-dimensional frost heave, the surfaces z _s and z _f are parallel planes. The frozen region including ice lenses moves upwards relative to z _b but does not deform. In this instance, the rheology of the frozen material is irrelevant. In DFH, however the frozen region deforms as it moves, and the rheology of the material must be specified so the mechanical problem can be solved in conjunction with the thermal and frost-heave problems.

4.2. Thermal problem

Due to the large Stefan number for freezing soil, the thermal problems in both the frozen and unfrozen regions above and below the frozen fringe, respectively, are greatly simplified. The two-dimensional, steady-state temperature distributions in both regions are obtained from the solution to the Laplace equation:

(1)

(2)

The boundary condition at the ground surface is Newton’s law of cooling. For both problems, the boundary condition at z _f specifies the temperature as the standard temperature and pressure (STP) freezing point of water. For this paper, the temperature at z _b will also be specified.

(3)

(4)

(5)

(6)

where H is a heat-transfer coefficient, T _air is the dimensionless air temperature at the ground surface and T _b is the dimensionless temperature of the underlying permafrost. In the absence of significant solutes, T _b = 0, corresponding to 0°C.

4.3. Mechanical problem

The choice of an elastic rheology can introduce several conceptual and modeling difficulties. To address these difficulties, we adopt several assumptions about the geometry of the frozen region as well as its mechanical properties. The frozen fringe is, mathematically, reduced to a plane across which jump boundary conditions are applied (Reference Fowler and KrantzFowler and Krantz, 1994). This plane is the lower boundary of the frozen region, z _f . However, because it will be perturbed and distorted, albeit infinitesimally, and will no longer be a plane in a strict mathematical context, we refer to it as an interface.

Figure 2 shows a schematic of the frozen fringe before and after inception of DFH due to initiation of a new ice lens. In fact, the location of the newly forming ice lens is very close to the base of the lowest ice lens relative to the thickness of the frozen fringe. There is a boundary layer in the frozen fringe near the lowest ice lens wherein most of the pressure drop causing upward permeation occurs. The new ice lens will form within this boundary layer (Reference Fowler and KrantzFowler and Krantz, 1994).

Fig. 2.

Diagram of the frozen fringe showing the initiation of DFH. Relative dimensions are not to scale.The dashed line in the left figure shows the location of z_f which is defined where the ice content is zero. A new ice lens forms within the water-pressure boundary layer below the lowest ice lens. After initiation of a new ice lens, the location of z_f will move either upwards or downwards relative to its previous position. The two possibilities are shown by and .The ice crystals labeled “possible ice” will exist for and not for .The new location of z_f is predicted by the LSA and is not prescribed in the model.

The dashed line labeled z _f is the bottom of the frozen fringe, defined as where the ice content is zero. At the inception of DFH, the location of z _f can move either up or down. Because of the quasi-steady-state assumption, any movement of z _f is relative to the stationary, one-dimensional boundary. Although z _f is always moving downwards on the time-scale of freezing-front penetration, quasi-steady state must be assumed in order to conduct the LSA. While this is an abstraction from reality, it is an expedient and useful practice for LSA of systems with unsteady basic states (Reference FowlerFowler, 1997). Two possible locations for z _f are shown in Figure 2 (right side): one for upward movement and one for downward movement. Which of these two possible outcomes actually occurs is not prescribed in the model but is a result of the LSA. However, linear theory does prescribe that z _f is either in phase or anti-phase with the ground-surface perturbation, z _s. Any relative displacement other than 0 or π is not allowed. The location of z _l in Figure 2 will be in phase with the ground surface at z _s since the lowest ice lens precludes water migration above it. Thus, Figure 2 demonstrates how z _s and z _f can be either in phase or anti-phase due to the formation of a new ice lens.

Mechanically, when non-uniform ice-lens formation begins, the frozen region will begin to deform and bend. We model this as the bending of a thin, elastic plate, though a viscous, viscoelastic or visco-plastic model might also be used. Laboratory experiments by Reference TystovichTsytovich (1975) demonstrate that frozen soil does behave elastically on short time-scales. The necessity for including the rheology of the frozen material in the model is to provide a “resistive force” to deformation of the frozen material. As will be demonstrated later, when the rheology is neglected, a LSA predicts that the most unstable mode has a wavelength approaching zero, which is intuitively incorrect.

A conceptual difficulty is encountered when assuming an elastic rheology in this problem since frozen material is being accreted. This problem can be avoided if we neglect the addition of newly frozen material in the geometry of the elastic material being bent. What this means physically is that when a non-uniform ice lens begins to form, only the material spanning from the ground surface, z _s to the lowest ice lens, z _l, is actually being bent. In the sense of the elastic model, the newly accreted material has an elastic modulus of zero. This approach is convenient because it avoids some inherent complexity, and there is considerable justification for its use. Experimental data for Young’s modulus of frozen soil as a function of temperature show that the modulus approaches zero as the sub-freezing temperature approaches zero (Reference TystovichTsytovich, 1975, p.185). Large amounts of unfrozen water at 0°C are the cause of this “weakness”. Thus, we would expect the newly accreted material to have a very low Young’s modulus, so that neglecting it completely introduces only a small error. In fact, as a first approximation the entire frozen region is assumed to have a constant Young’s modulus. This approximation is crude since there is a non-zero temperature gradient within the frozen region. In reality, the Young’s modulus monotonically increases with decreasing temperature and has a maximum value at the ground surface where it is coldest.

Due to the simplifications discussed above, the mechanical model reduces to the bending of a thin, elastic plate with a constant Young’s modulus. If the deformation is small relative to the thickness of the plate, the deflection of the plate’s mid-plane from equilibrium, w, can be approximated as (Reference Brush and AlmrothBrush and Almroth, 1975):

(7)

(8)

(9)

where P _b is the pressure, E is Young’s modulus, h ₁ is the frozen-region thickness and ν is Poisson’s ratio. In this analysis, ν = 0.5.

Boundary conditions for Equation (7) must be specified in the x and y directions. Because the problem is unbounded in these directions, periodic boundary conditions are used. For the x direction,

(10)

(11)

(12)

(13)

There is a corresponding set for the y direction (however, only the two-dimensional problem will be considered here). The mechanical problem is coupled to the frost-heave problem through the load pressure at z _l. In the one-dimensional basic state, the frozen region is planar and there is no bending (P _b = 0). In the two-dimensional state, the non-zero P _b is a function of x. The load pressure is the sum of the weight of the overlying soil and the pressure due to bending.

(14)

4.4. Non-dimensionalized frost-heave equations

The Fowler and Krantz frost-heave equations are two coupled vector equations that describe the velocity of both the upward-moving ice at the top of the lowest ice lens, v_i , and the downward-moving freezing front, V_f . Specifically, the equations are:

(15)

(16)

where G_i and G_f are the thermal gradients in the unfrozen and frozen regions, respectively, φ is the soil porosity, W _l is the unfrozen-water volume fraction at the top of the frozen fringe (base of the lowest ice lens), and k _f and k _u are the frozen and unfrozen thermal conductivities, respectively. Vectors are denoted in bold, and the coordinate system is oriented with z pointing up as shown in Figure 1. The value of W _l is determined using the lens-initiation criterion (Reference Fowler and KrantzFowler and Krantz, 1994) that involves the characteristic function.

(17)

The solution to the one-dimensional problem is fairly straightforward.The steady-state, one-dimensional Laplace equation for temperature is solved to determine the relevant temperature gradients G_i and G_f . Work by Reference Fowler and NoonFowler and Noon (1993) and Reference Krantz and AdamsKrantz and Adams (1996) has shown that this model predicts the one-dimensional frost-heave behavior observed in the laboratory quite well.

6.1. Effects of parameter variations

In the following subsections, the effects of varying different parameters from their reference-set values are investigated. The results are presented by plotting the dimensionless growth rate, Γ, as a function of the dimensionless wavenumber, α (see Equation (19)). The scaling factors for these parameters are

(34)

(35)

For comparison purposes, the results will be presented in dimensionless form. The stability criterion, Equation (33), is quadratic, and therefore there are two roots Γ₊, Γ₋, with Γ₊ > Γ₋. In this analysis, the lower root Γ₋ was always negative. The positive value, Γ₊, is plotted for the reference set using Chena Silt in Figure 3. In all cases, both roots are real numbers for the range of α values presented, and the principle of exchange of stabilities is valid (Reference Drazin and ReidDrazin and Reid, 1991).

Fig. 3.

Plot of the positive root of the linear stability relation as a function of the wavenumber for the reference set of parameter values.The soil is Chena Silt.Wavenumbers greater than α_ntl are stable and α_max is the wavenumber with the maximum growth rate under these conditions.

There are three values that are of most importance in each plot: Γ_max, α _max and α _ntl, indicated in Figure 3. The maximum value of the growth rate, Γ_max, and the corresponding value of the wavenumber, α _max, indicate which mode grows fastest in the linearized model. Also, α _max indicates the lateral spacing of the most highly amplified, two-dimensional mode through Equation (35). The neutrally stable wavenumber, α _ntl, indicates the maximum wavenumber for which modes can be linearly unstable. These values are not explicitly marked in subsequent figures, to reduce clutter.

6.2. Effect of soil type

The type of soil is the most influential parameter in the model when determining the stability of one-dimensional frost heave. It is important to note that only a small fraction of soils exhibit frost heave (Reference Williams and SmithWilliams and Smith, 1989). A delicate balance between the hydraulic conductivity of the partially frozen soil, particle size and its porosity is necessary for frost heave to occur. Silts and silty clays meet these requirements in general. Furthermore, DFH is not necessarily observed in all frost-susceptible soils. Determining whether this observation is due to the specific soil or its environmental conditions is a major objective of this analysis. Figure 4 plots the growth rate as a function of wavenumber for the three soil types specified in Table 1.

Fig. 4.

The linear growth coefficient for three soil types using the reference set of parameter values.

All three soils are considered frost-susceptible. However, for the reference set of parameter values, only Chena Silt has a propensity for DFH (i.e. positive value of Γ₊). It should be noted that these results do not indicate that Calgary Silt and Illite Clay will never heave differentially. Other parameters can have a significant effect on the stability predictions, and under the correct conditions the soils may heave differentially.

The utility of a model such as this one hinges on whether it can predict observations made in the field. Unfortunately, there is a limited number of soil types to analyze due to a dearth of adequate characterization data including hydraulic conductivity and sub-freezing temperature as functions of unfrozen water content. In fact, the work by Reference Horiguchi and MillerHoriguchi and Miller (1983) is one of the few comprehensive datasets for frost-susceptible soils available in the open literature. Fortunately, the predictions shown in Figure 4 do agree with observations made in the field. In almost all instances, hummocks are observed in silty-clayey materials. The value of α _max for Chena Silt is of the order unity. This wavenumber corresponds to a dimensional wavelength on the order of meters, which is also observed in the field (Reference Williams and SmithWilliams and Smith, 1989). The value of Γ_max is several orders of magnitude less than unity, indicating that DFH will take (dimensionally) several years to develop. It is commonly believed that most types of patterned ground form on the 10–100 year time-scale (Reference MackayMackay, 1980; Reference Hallet, Nicolis and NicolisHallet, 1987).

6.3. Effect of the elastic constant

In this and subsequent subsections, we investigate the effects of parameter variations away from the reference set for Chena Silt, which has the greatest propensity for DFH. Although calculating the effect of the elastic constant on linear stability is not difficult, interpreting the results can be confusing due to some of the assumptions made in the mechanical model. The complication arises due to the difficulty in determining an appropriate value for the Young’s modulus, E. While there are several references that address the Young’s modulus for frozen soils (Reference TystovichTsytovich, 1975; Reference Puswewala and RajapaksePuswewala and Rajapankse, 1993; Reference Andersland and LadanyiAndersland and Landanyi, 1994; Reference Yuanlin, Ping, Jiayi and JianmingYuanlin and others 1998), there is a large range of reported values. Reference Andersland and LadanyiAndersland and Ladanyi (1994) report that for frozen silts, sands and clays, a Young’s modulus of about 5 MPa is typical for the temperature range of interest: −10 → 0°C. The effect of Young’s modulus on the growth-rate predictions is shown in Figure 5.

Fig. 5.

The linear growth coefficient, Γ₊, as a function of wavenumber, α, for three values of the Young’s modulus, E.The soil type is Chena Silt.

As might be expected, the range of unstable modes decreases with increasing Young’s modulus. The value of α _max decreases, indicating that a “stiffer” soil results in less closely spaced patterning. Because Γ_max decreases with increasing E, it would also take longer for patterns to become expressed in stiffer soils. It appears that increasing elasticity is a stabilizing mechanism. As it becomes more difficult to “bend” the frozen region, the propensity for differential heave is reduced.

It is worth noting that only in the limit of E → ∞ does the propensity for DFH at some value of α disappear. Both α _ntl and α _max approach zero only when E → ∞. One major cause for larger values of E for a particular soil would be colder temperatures. However, colder temperatures also affect the thermal regime in the frozen region.This effect is explored separately in section 6.4 and 6.5.

6.4. Effect of the ground surface heat-transfer coefficient

In all the previous results, we have assumed a constant temperature of −10°C at the ground surface. However, a constant-temperature condition is unrealistic for prolonged periods of time (i.e. > 1–2 days). A lumped parameter boundary condition with an overall heat-transfer coefficient is more reasonablefor long-term analysis.The constant-temperature condition is obtained from the heat-transfer-coefficient condition by allowing H → ∞. Figure 6 plots Γ₊ as a function of α for three values of the dimensionless heat-transfer coefficient. A dimensional value of the heat-transfer coefficient, h, can be obtained by reversing the scaling .

Fig. 6.

The linear growth coefficient, Γ₊, as a function of wavenumber, α, for three values of the dimensionless heat-transfer coefficient . The value of Γ_max for H = 10 is two orders of magnitude greater than for a constant-temperature boundary condition (H → ∞).

It is readily apparent that the “insulating” effect represented by smaller values of H causes a greater range of unstable modes. Furthermore, as H decreases from ∞ to 10, Γ_max increases two orders of magnitude. Insulation obviously has a destabilizing effect.

Both α _max and the corresponding value of Γ_max increase as H decreases. There are several scenarios in the field in which small values of H might be applicable. The most obvious instance is that of increased ground cover in the form of grasses, mosses and shrubs. It is noteworthy that the extensive survey of earth hummocks byReference Tarnocai and ZoltaiTarnocai and Zoltai (1978) supports this prediction. In essentially all places where earth hummocks are found, there is a moderate amount of organic ground cover. The fact that the tops of hummocks are sometimes more barren than the trough regions is most likely the result of changes in soil moisture after the hummock shape has matured. It is reasonable to assume that the differential height between the crest and the trough would cause water drainage from the top that collects in the inter-hummock spaces, thus propagating further organic growth in the depressions.

6.5. Effect of ground-surface temperature

Snow cover, which is almost always present to some degree in arctic and subarctic regions, also has an insulating effect. However, snow is more likely to maintain a constant temperature at the ground surface due to its ability to function as a “heat sink” because of the relatively large latent heat associated with freezing. To analyze the effect of snow cover, it is more useful to return to the constant-temperature boundary condition and vary the value of T _air. Since H → ∞ for the constant-temperature condition, the ground-surface temperature is equivalent to T _air. A smaller (in magnitude) value of T _air corresponds to a thicker and/or less thermally conductive snow layer. In fact, the only limitation in this case is T _air < 0. The effect of the ground-surface temperature is shown in Figure 7.

Fig. 7.

The effect of warmer constant-temperature conditions on the growth of differential modes for the reference set of parameter values.The growth rate, Γ is plotted in dimensional form.

In Figure 7 it can be observed that as the ground-surface temperature warms, Γ_max decreases. Furthermore, in the limit of T _air → 0, the growth coefficient for all unstable modes tends to zero. These results indicate that the occurrence of DFH should be more prevalent in less snow-covered areas than those where snow accumulates, and in fact this is observed. A good example is the Toolik Lakes (Alaska, U.S.A.) field site where D. A. Walker (personal communication, 2000) has observed that frost boils tend to be more prevalent in the wind-blown areas. Frost boils are small mounds of soil material, ∼1 m in diameter, presumed to have been formed by frost action (NSIDC, 1998). Ground-surface temperature measurements in both regions confirm that surface temperatures are indeed colder in the barren regions.

6.6. Effect of freezing depth

It is appropriate to investigate the effect of freezing depth, h ₁ ≡ z _s − z _f, on the stability theory predictions.The non-dimensionalization establishes that h ₁ is of the order one.When ground freezing commences, h ₁ begins at zero and increases (approximately with the square root of time (Reference Fowler and NoonFowler and Noon, 1993)). By examining the behavior of α _max and Γ_max for increasing values of h ₁, we are able to observe how the stability of the system changes as ground freezing progresses. Figure 8 shows the effect of freezing depth on the stability analysis predictions. As one might expect, the largest range of unstable modes occurs at the smallest values of the freezing depth and slowly decreases as h ₁ increases. However, contrary to previous results where unstable modes existed for all values of a parameter, stabilization occurs at a finite value of h ₁. For the reference set of parameter values shown here, that value is h ₁ ≈ 0.23, which corresponds to a dimensional value of 23 cm.

Fig. 8.

The linear growth coefficient, Γ₊, as a function of wavenumber, α, for three depths of freezing, h₁.There are no positive roots (Γ₊, Γ₋ <0) at freezing depths greater than h₁ ≈ 0.23.

The fact that there is a critical value of the freezing depth is important when attempting to gauge what wave-number might commonly be most expressed in the field. In the limit h ₁ → 0, all wavenumbers have a positive growth rate. As freezing progresses and h ₁ increases from zero, the largest wavenumbers are “cut off ” (i.e. Γ₊ < 0) while the smaller wavenumbers continue to grow, albeit at a slower rate. As h ₁ continues to increase, more and more wavenumbers are cut off. Eventually when h ₁ ≈ 0.23, all wavenumbers are cut off and no differential modes are unstable. It can be seen in Figure 8 that very small wavenumbers always have small positive growth rates until cut off. Large wavenumbers have large growth rates but are cut off early. Mid-range wavenumbers have moderate growth rates for intermediate periods of time. Thus, mid-range wavenumbers might become most expressed in the field since they have a moderate length of time to grow and moderate growth rates for most of the time.

Although linear theory is not strictly valid once differential modes grow to a finite amplitude, it is insightful to see what wavenumber the model predicts would become most expressed during the entire period of freeze-up. This assumption is not excessively crude since the dimensionless growth rates are much less than unity in this case, and the perturbations are going to grow relatively slowly. As a simple approximation, we define the normalized growth rate as follows:

(36)

where ζ is the perturbation in the ground-surface location. The magnitude of ± at time t′ is determined by integration:

(37)

(38)

where ζ ₀ is the magnitude of the original perturbation at t = 0.

Figure 9 shows the numeric results of integrating Equation (38) for several values of the wavenumber α, indicated by the dots. The integration is performed until a time t′, at which point all modes become stable (i.e. h ₁ ≈ 0.23). The most expressed wavenumber under these conditions is α = 2.3 which corresponds to a dimensional wavelength of 2.7 m. This result is quite encouraging since the spacing of hummocks is typically 1–3 m (Reference Tarnocai and ZoltaiTarnocai and Zoltai, 1978). Because all modes are stable for times greater than t′ (or h ₁ > 0.23), a further increase in frost depth will result in all modes decreasing in amplitude.

Fig. 9.

Relative amplitude of the ground-surface perturbation, ζ/ζ₀, for several values of the wavenumber, α, at the time t′ when all modes become stable. Soil type is Chena Silt and H = 10. At t = t′, h₁ ≈ 0.23.

To determine the evolution of ζ/ζ ₀ in time as h ₁ increases, Equation (38) is integrated in time up to t′ for a given wavenumber. The results for α = 2.6 are shown in Figure 10. This figure shows the evolution of ζ/ζ ₀ as freezing progresses (i.e. h ₁ increases). The normalized perturbation amplitude reaches a maximum at h ₁ ≈ 0.08 and then begins to decrease. At h ₁ = 0.2, ζ/ζ ₀ ≈ 2.65 as expected from Figure 9. The amplitude returns to its initial value at h ₁ ≈ 0.36 and continues to decrease as h ₁ increases further. These results indicate that an initial perturbation of magnitude ζ ₀ will have increased in magnitude at the end of the freezing season for active-layer depths < 0.36 m.

Fig. 10.

Evolution of the normalized ground-surface perturbation as a function of the freezing depth, h₁.

There is no assumption in this derivation that the perturbation grows as exp[Γ t] as linear theory requires. In this sense, this is a real-time analysis. However, other assumptions have implications that must be examined. Perhaps most importantly, all product terms in perturbed variables are assumed small and neglected. Hence, this real-time analysis is still only applicable when all perturbations are relatively small in magnitude.

6.7. Mechanism for instability

It has been demonstrated above that there is a range of modes that are unstable for most cases of interest. In fact, if the Young’s modulus is allowed to approach zero (E → 0), all modes are unstable for some soils while no modes are unstable for other soils. In this limit, the bending resistance of the frozen soil is essentially removed (see Equation (7)). Figure 11 shows the LSA prediction for the reference set of parameter values and Chena Silt except that E = 0. It is evident that the elasticity of the frozen layer is providing the stabilizing mechanism for DFH.

Fig. 11.

Stability predictions for the fundamental case when theYoung’s modulus is set to zero (E = 0).

In order to help explain the mechanism that causes the instability, Figure 12 shows a schematic of the in-phase perturbed surfaces. Also sketched are approximate isotherms for the constant-temperature boundary condition at the top surface. The large arrow pointing to the right shows the net direction of differential heat flow that is occurring from a crest region to a trough region. Before the perturbations occur, heat transfer occurs only in the vertical direction. However, once the system is perturbed, heat transfer can occur in both the horizontal and vertical directions.

Fig. 12.

Schematic of the perturbed surfaces showing isotherms and the direction of differential heat flow.The basic state surfaces and are shown by the dashed lines, and the perturbed surfaces are the thick solid lines.

Solving Equation (33) for η/ζ indicates that the magnitude of the bottom surface perturbation, η, is less than the top surface perturbation, ζ, albeit by a small difference. However, linear stability theory predicts perturbations grow exponentially. Thus, although the difference between η and ζ is small initially, that difference will grow exponentially fast (at least until non-linear terms become significant). In fact, it is difficult to speculate about what form the perturbations will take once non-linear terms become significant. Figure 12 shows that the thickness has increased from the basic state value underneath a crest, and decreased underneath a trough. Because the temperatures at the boundaries z _f and z _s are fixed at the basic-state values and the conduction path length has changed, the temperature gradient at the bottom surface, Gi, has also changed. Underneath a crest the gradient has decreased in magnitude and, according to Equation (15), the velocity of the freezing front will also decrease. Continuing this line of reasoning, since the conduction path underneath a trough has decreased, the gradient at z _f has increased and the velocity of the freezing front has correspondingly increased. Since there is a positive feedback mechanism occurring at this point, the perturbations will continue to grow.

In order to understand the increase in Γ with increasing wavenumber observed in Figure 11, it is necessary to consider the magnitude of differential heat flow shown in Figure 12. As the wavenumber increases, the magnitude of heat flux from a crest region to a trough region increases. In order for frost heave to occur, the latent heat of the incoming water at z _l must be removed. Since there is differential flow from a crest to a trough, additional heaving can occur underneath the crest region. This also corresponds to less heave occurring in the trough region.

6.8. Implications of the steady-state assumption

This analysis is based on the quasi-steady-state assumption that perturbations grow much faster than changes in the basic state. However, predictions indicate that the growth rate Γ is actually less than unity in most cases, and thus perturbations grow slowly relative to changes in the basic state. For some time-periodic problems with slow growth rates, Floquet theory has been successfully applied when solving the linear stability problem (Reference Drazin and ReidDrazin and Reid, 1991, p. 354). However, the entire freeze–thaw cycle (or period) must be modeled. This analysis has only addressed the frost-heave process that occurs during freezing, which is only part of a larger cycle. The thawing process involves solifluction, soil consolidation, water percolation and evaporation, and possibly other processes.

The present analysis indicates that the frost-heave process alone is capable of causing spontaneous pattern generation, albeit on a 10–100 year time-scale. Combination of the frost-heave model with a thaw model could possibly yield more conclusive theoretical evidence that patterns are generated by the processes discussed here. It may be that the mature periglacial patterns observed are actually the result of a complex combination of all the processes mentioned above.