Ice cascade growth and decay: a thermodynamic approach

Abstract

The ice volume evolution of a frozen waterfall (or ice cascade) was studied using a thermodynamic model. The model was developed from meteorological data collected in the vicinity of the waterfall and validated from ice volume measurements estimated from terrestrial lidar images. The ice cascade forms over a 45 m high rock wall located in northern Gaspésie, Québec, Canada. Two stages of formation were identified. During the first stage, the growth is mainly controlled by air convection around the flowing and free-falling water. The ice cascade growth rate increases with decreasing air temperature below 0°C and when the water flow reaches its lowest level. During the second stage, the ice cascade covers the entire rock-wall surface, water flow is isolated from the outside environment and ice volume increases asymptotically. Heat is evacuated from the water flow through the ice cover by conduction. The growth is controlled mainly by the conductive heat loss through the ice cover but also by the longwave radiation emitted at the ice surface during the night. In spring, melting of the ice cascade is dependent on the air convection over the ice surface but also on the sensible heat carried by the increasing water flow and the solar radiation received during the day.

List of Symbols

Introduction

In northern Gaspésie, eastern Canada, the regional road is enclosed between a series of relatively high rock walls and the Saint Lawrence estuary. On those cliffs, many waterfalls lead to the formation of ice cascades during the cold season. Each spring there are a number of ice-block falls, which present a traffic hazard, as reported by the Quebec Department of Transport (MTQ) (Girard and Hétu, 1994; Hétu and others, 1994; Gauthier, 2008; Gauthier and others, 2012). Additionally, with an increasing interest in mountaineering and ice climbing in North America, the number of traumatisms reported since 1951 has increased (American Alpine Club and Alpine Club of Canada, 1999). Of the 200–400 accidents reported each year, 20% were caused during ice climbing. Although the mechanical (in)stability of vertical ice cascades has been studied recently in the Alps (Weiss and others, 2011), a better understanding of this natural hazard is still needed.

Recent studies show that air temperature is certainly the main variable controlling the growth and decay of ice cascades (Gauthier, 2008; Montagnat and others, 2010; Gauthier and others, 2012). The statistical models developed by Gauthier (2008) and Gauthier and others (2012) to predict ice-block falls along the northern Gaspésie road show that rising air temperature above the melting point is the main variable controlling the melt and collapse processes of the ice cascades. The model based on melting degree-days calculation yields acceptable correlation coefficients with the occurrence of ice-block falls. Montagnat and others (2010) proposed a simple conductive heat flux model to quantify the evolution of a frozen waterfall in the French Alps. Their model is based on the crude assumption that the water is flowing slowly between the growing ice cover and the rock-wall surface. However, observations show that the growth dynamics is much more complex. Ice cascade formation begins with the growth of a multitude of independent ice stalactites (Gauthier, 2008; Montagnat and others, 2010; Gauthier and others, 2012) and with the formation of ice accretions caused by the freezing of water spray around the main water flow (Gauthier, 2008; Gauthier and others, 2012). These observations show that the growing and melting processes result from a complex mix of processes common to ice stalactites (e.g. Laudise and Barns, 1979; Knight, 1980; Makkonen, 1988; Maeno and others, 1994; Short and others, 2006), growth of icings (aufeis) (e.g. Schohl and Ettema, 1986, 1990; Hu and others, 1999), ice accretions similar to freezing rain accumulation (e.g. Jones, 1996; Makkonen, 2000) and frazil-ice formation (e.g. Hanley and Michel, 1977; Osterkamp and Gosink, 1983; Andres, 1995; Ye and Doering, 2004). Energy-balance assessments of these superficial ice formations have already been undertaken by the authors cited, but no specific model has yet been proposed for assessing parameters and simulating the growth and decay of ice cascades.

In this paper, we focus on hydrometeorological data collected around an ice cascade and we propose a two-stage growth-and-decay model based on observations and data analysis. We are not considering the mechanical destabilization that frequently occurs during large and rapid temperature changes as the result of thermal stresses (Weiss and others, 2011). Instead, we propose an energy-balance model based on, and developed from, a large number of hydrometeorological data collected around a typical north-facing ice cascade. The model is calibrated from ice volume measurements made from a series of terrestrial lidar images. The model clarifies the role of the different heat transfer components involved in ice cascade growth and decay. The following questions are addressed: Is the conductive heat flux coming from the rock wall a major source of heat? Does latent heat flux (evaporation and sublimation) play a significant role in ice cascade thermodynamics? What is the amount of ice produced during the fall of the water down the cliff involved in the formation of the ice cascade? Is the sensible heat brought by the flowing water a major energy flux? Does the radiation heat budget play a significant role in the ice cascade growth-and-decay process? The main objectives of this study are to identify the growth-anddecay mechanisms, to evaluate the importance of the different heat fluxes controlling the ice cascade growth and decay and to improve our general understanding of this complex system.

Location

Fieldwork was conducted during the 2010/11 winter on an ice cascade named ‘Le Voile de la Mariée’ (G-Voile) located along the north shore of the Gaspé Peninsula, eastern Canada (Fig. 1). The climate of the region is cold temperate with a strong maritime influence. The average daily minimum and maximum temperatures in January are −16.1°C and −7.2°C, respectively (Gagnon, 1970; Environnement Canada, 2000). The maritime influence of the Gulf of Saint Lawrence provides an average rainfall of ∼1000 mm spread evenly throughout the year. One-third of this precipitation falls as snow from mid-November to late April. Liquid precipitation during the fall is also relatively high, with ∼200 mm for October and November. Along the coast, the rock walls are strongly affected by the wind, with only 14% of calm days per year and with a wind speed frequently exceeding 50 km h⁻¹ (Gagnon, 1970; Hétu and Vandelac, 1989).

Fig. 1. Location of the Voile de la Mariée ice cascade (waterfall) and its catchment (CA).

The ice cascade forms on a 45 m high north-facing (354°) and almost vertical rock wall (angled at ∼80°) composed mainly of stratified shale, limestone and sandstone (Brisebois and Nadeau, 2003). The ice cascade does not show freestanding structure such as a suspended ice structure or ice column. The water comes from a small stream draining a 1.8 km² catchment. The mean winter discharge of the stream is ∼14 L s⁻¹ but drops to 6–7 L s⁻¹ in January and February. The peak event arises during the melt season (>250 L s⁻¹ in May 2011).

Measurements

All instruments were installed during summer 2009 and tested during winter 2009/10. The final set-up was completed during summer 2010 and the measurements were conducted during the 2010/11 cold season from early October to the end of June. Instantaneous measurements were sampled every 10 min but we used hourly mean values for data analysis and to develop the model.

Instrumentation

Two independent and automatic weather stations (Onset: HOBO Micro Station Data Logger, H21-002) were set up around the ice cascade. The first station (WS1) was located at the top of the waterfall (Fig. 2). The air temperature, relative humidity, air pressure, precipitation and solar radiation were measured with a 12-bit temperature/RH smart sensor (S-THB-M002), a barometric pressure smart sensor (S-BPB-CM50), a 0.2 mm rainfall smart sensor (S-RGB-M002) and a solar radiation sensor (Silicon Pyranometer, S-LIB-M003). The pyranometer was set up at the same vertical angle and orientation as the ice cascade. Its measurement range is from 0 to 1280 W m ⁻² across a spectral range of 300–1100 nm.

Fig. 2. Ice cascade evolution between 12 December 2010 and 6 May 2011. The positioning of the weather stations is shown (WS1 and WS2).

The second station (WS2) was located ∼50 m east of the middle of the waterfall and at ∼4.5 m from the rock wall. The wind speed and rock-wall temperature were measured with a wind-speed and -direction smart sensor (S-WCA-M003) and three 12-bit temperature smart sensors (S-TMB-M017). The temperature sensors were fixed on the rock-wall surface and embedded 38 and 80 cm deep into the rock.

Water pressure and temperature were measured at the top and bottom of the waterfall with two HOBO U20 water-level data loggers (U20-001-04). Unfortunately, the pressure sensor at the foot of the waterfall did not produce reliable data because it was displaced by a storm surge in December 2010. The water temperature measurements were still valid, but the sensor was located ∼10 m away from the foot of the waterfall. Air pressure was subtracted from the water pressure measurements. As proposed by the Federal Office of Metrology (METAS), we took into account an altitude air pressure correction factor of the form

(1)

where P _wll is the atmospheric pressure at the water-level height (Pa), P is the atmospheric pressure at the barometer (Pa), T is the air temperature (°C) and H is the height difference between the water-level logger and the barometer. The corrected water pressures were then used to develop a discharge calibration curve.

Discharge measurements

Discharge measurements were made using the slug salt injection method (Moore, 2005). We used 1 kg of salt diluted in 900 mL of water. The solution was dumped ∼50 m upstream of the cascade. Less than 1 m from the top of the waterfall, we measured the electrical conductivity of the stream water using a YSI handheld multiparameter instrument (YSI 556 MPS). All measurements were calibrated to determine the relation between the relative salt concentration and the water electrical conductivity. Moore (2004a,b, 2005) provides a more detailed description of the method.

Under suitable conditions, streamflow measurements made by slug injection can be accurate at about ±15% (Day, 1976; Moore, 2005). Measurements were made two to three times a month depending on stream conditions (e.g. ice formation) and discharge fluctuations, and in total 22 measurements were made during the study (between October 2010 and June 2011). Four measurements had to be eliminated because mistakes were made during handling. The 16 reliable discharge measurements ranging from 2.0 to 215.1 L s⁻¹ were kept to build the calibration curve (Fig. 3). A second-order polynomial relation provides the best fit to these data. It is important to note that above the maximum discharge value (215.1 L s⁻¹), discharges are extrapolated and could be overestimated. However, as the peak discharge comes after the complete melting of the ice cascade, this extrapolation does not affect the results presented in this paper.

Fig. 3. Discharge calibration curve.

Ice volume extraction

Ice-volume changes due to accretion and decay were estimated using terrestrial lidar measurements. The scans were performed using a Leica Scan Station 2. We first scanned the rock-wall (waterfall) surface in September 2010 when the water discharge was at its lowest level. Subsequently, 11 scans were performed during the ice cascade growing and melting season. The scanner works up to 300 m for a surface albedo of 0.9, and 134 m for a surface albedo of 0.18. The scans were performed at a distance of 30 m from the base of the ice cascade and ∼60 m from the top. The laser beam did not penetrate the dry, relatively porous and irregular ice cascade surface as shown by the absence of scatter in the point measurements at the ice surface. The laser beam diameter is ∼4 mm and the maximum sampling density is <1 mm. The accuracy of the lidar is 4 mm in distance and 6 mm in positioning if located at <50 m from the scanned surface. The scans were performed using a 5 mm resolution at a distance of 45 m.

The next step consisted of developing high-resolution triangulated irregular network models (HRTINMs) from the lidar measurements. The HRTINMs were developed using the Leica Geosystems software Cyclone 7.1. The registration function in Cyclone was used to merge the point cloud into the same field of coordinates. The use of five fixed survey targets allowed an overlapping precision of <5 mm. For better consistency between HRTINMs, the point clouds were resampled at 20 mm resolution. We were also forced to trim the HRTINMs to eliminate the bordering vegetation and the snow accumulation at the foot of the ice cascade. Doing so, the calculated ice volume does not represent the total ice accumulation but rather a lower estimate (probably between 80 and 90%). The same surface area was extracted from each HRTINM. At the end, we extracted the HRTINMs volume using the function Mesh Volume in Cyclone 7.1. The given volume is the difference between the HRTINMs and a reference plan (Fig. 4). The HRTINMs generated from the September 2010 lidar measurements were used as the reference surface topography. The extracted volume from this reference surface was then subtracted from the other calculated volumes. The differences represent the ice volume accumulation at the time of the scan. Rabatel and others (2008) and James and others (2006) provide a more detailed description of the method.

Fig. 4. Computation examples of the volume of the ice cascade during winter 2010/11. The volume between the surface topography and the reference plan is computed for each survey date using the function Mesh Volume in Cyclone 7.1. The difference between the volume calculation across surveys and the volume of the reference surface topography (9 September 2010) represents the change in ice volume.

The total uncertainty (σ _tot) can be estimated (Baltsavias, 1999; Rabatel and others, 2008) using

(2)

where σ ₁–σ ₄ represent the different independent processing errors. σ ₁ is the error inherent to the lidar itself. σ ₂ is the error in the development of the HRTINM. It corresponds to the point-cloud resolution chosen to generate the HRTINM: 20 mm in our case. σ ₃ is the error in the point-cloud merging (overlapping) process (5 mm). σ ₄ represents the errors related to the scan surface itself and to the environment. Because of the good climatic conditions during the scanning acquisition process, except for scans on 10 December 2010, the error related to the environment is certainly negligible. On 10 December 2010, light snowfall reduced the quality of the acquired point-cloud image. Many points representing snowflakes had to be removed before the generation of the HRTINM. The few remaining points, too close to the ice surface to be removed, caused a slight overestimation of the ice volume. Furthermore, the error associated with the scan surface is difficult to estimate. Even if the ice cascade surface was generally dry, some limited zones were lubricated by localized flow on the ice surface. The point-cloud resolution over these wet zones was not as high as on the remaining surface. It affected the size of the mesh locally, but we were unable to estimate the exact error associated with this phenomenon. This leads to a minimum uncertainty of 22 mm over the trimmed surface area, representing an absolute error of 6.5 m³ on the calculated ice volumes or 115% of the smallest calculated ice volume (5.5 m³ on 26 November 2010) to 2.8% of the largest calculated ice volume (231 m³ on 30 March 2011).

The lidar measurement only leads to a precise reconstruction of the external shell of the ice cascade. It is not possible to estimate the volume of the voids between the ice cascade and the rock wall, especially if ice melting occurs at the ice– rock-wall interface due to a positive sensible heat flux carried by the flowing water.

Hydrometeorological Conditions and Ice Cascade Evolution

We focus here on the local hydrometeorological conditions in the vicinity of the ice cascade. We also relate the evolving conditions to the ice cascade volume evolution. Figure 5 shows the hydrometeorological data recorded during the study commencing on 15 November 2010 (day 0). Because the ice cascade is facing north, only diffuse radiation reaches the ice surface during winter. There is still a daily fluctuation, but its order of magnitude is small (Fig. 5a). At the end of January (∼day 65), the incoming diffuse radiation slowly rises. The waterfall was submitted to a short period of solar radiation during late afternoon at the beginning of May, but at that time the ice cascade had almost entirely disappeared. Unsurprisingly, the relative humidity remains rather high during the cold months of winter and tends to decrease slightly when the temperature increases in spring (Fig. 5b). The wind conditions are almost never completely calm and the frequency of wind events with hourly mean values close to, or over, 10 m s⁻¹ (36 km h⁻¹) is high (Fig. 5b). We assume that the wind plays an important role in the sensible heat transfer at the water–air and ice–air interfaces. We also assume that relatively low humidity favors evaporation/sublimation and heat transfer to the air. Also, presumably, the slowly increasing incoming solar radiation will transfer energy to the system.

Fig. 5. (a) Global solar radiation; (b) wind speed and relative humidity; (c) air and water temperature; (d) precipitation and water discharge; and (e) freezing degree-days and measured ice cascade volume. Period of measurements is from 15 November 2010 (day 0) to 31 May 2011 (day 197). The dashed gray arrow represents a more probable ice volume evolution.

The first sustained period with air temperature below 0°C seems to control the beginning of the ice cascade growth (Fig. 5c). For example, a period of 4 days (day 4 to day 8) with air temperature below 0°C was sufficient to initiate the formation of ice structures on each side of the waterfall (Fig. 5c; lidar image of 26 November 2010 in Fig. 4). On the other hand, a drastic increase in temperature at approximately days 28–29 resulted in the melting of the young ice cascade. These observations suggest the calculation of a freezing or melting degree-hours parameter (FDH):

(3)

where T _f is the freezing point of water and T _a is the air temperature. The FDH was calculated at an hourly time-step (t) starting on 15 November 2010 at 00:00 with a minimum value of 0. Similar calculations but at a daily time-step were previously performed to model the growth of ice cascades in the French Alps (Montagnat and others, 2010) and also to predict the collapse of specific ice walls in northern Gaspésie (Gauthier, 2008; Gauthier and others, 2012). At the beginning of the growing process, FDH shows a significant lag compared with the measured ice volume accumulation rate (Fig. 5e). FDH also responds rapidly to air temperature variations. We think that the fragmentation of the water flow during its free fall increases convective heat transfer with the atmosphere. Under low temperature (<0°C), ice formed on each side of the waterfall due to cold air convection around the free-falling water droplets, thus producing important ice accretion. When air temperature increases (>0°C), heating is controlled by air convection around the static ice formations, thus limiting the melting of ice. This imbalance between ice formation and melting could explain why the FDH does not adequately reproduce the initial growth-and-decay dynamics.

Although the water temperature is affected by air temperature variations, we cannot at this stage ascertain whether the ice volume evolution is also dependent on these fluctuations. Most likely, the sensible heat flux carried by the flowing water will play a role in the global evolution of the ice cascade. Three ad hoc water temperature measurements made at the beginning of the growing process show that the water at the foot of the waterfall was in a supercooled state (between −0.01°C and −0.06°C ± 0.015°C) even when the air temperature was just below 0°C. These measurements are not in agreement with the bottom water temperature presented in Figure 5c. We think that the logger was located too far from the foot of the ice cascade to record this information. In steep stream systems, supercooled temperatures (usually with maximum value of ∼−0.04°C) are observed only in frazil production zones such as rapids or riffles (e.g. Stickler and Alfredsen, 2009). At the foot of the 45 m high waterfall, a huge amount of frazil ice was observed but it was almost entirely evacuated by the current. Only a small amount of ice seems to accumulate on the branches and little on the rock structures located around the main water flow path. Ice accretions (similar to freezing rain accumulation; Makkonen, 2000) and ice stalactites grow slowly on and under those structures (Fig. 2).

Even if the discharge slowly decreases during this period, the limited amount of liquid precipitation received continues to affect the streamflow, bringing more sensible heat to the system (Fig. 5d). But because the ice is formed mainly on the edges of the main waterfall flow, the slowly growing ice structures are not significantly affected by these fluctuations. On the other hand, snowfall events (e.g. 26 cm received on 26 January 2011, day 72) can presumably favour the growth process by lowering the water temperature. In this scenario, the melting snowflakes in the water flow might act as ice nuclei.

As soon as the air temperature is maintained under the freezing point of water and discharge reaches its base flow (∼6–12 L s⁻¹), the ice cascade forms rapidly (Figs 2 and 5). The ice structures on the edges of the main flow gradually grow from the periphery inward to the center of the decreasing flow. Eventually, the ice covers the entire rock-wall or waterfall surface. Although the FDH index shows a significant lag (∼10 days) compared with the ice volume evolution, it still follows the fast-growing trend of the ice cover.

When the ice covers the entire waterfall area, the main water flow is isolated from the outer environment. There is still some external flow on the ice cascade, mainly by exfiltration, but it remains marginal and localized. Observations show that, between 12 January (day 58) and 5 February 2011 (day 82), ice formed over the entire waterfall (Figs 2 and 4). The transition from a regime controlled by the air–water interaction to the complete isolation of the main water flow is fast (<6 days). Between 18 and 23 January 2011 (days 54–69), the water temperature at the foot of the waterfall increased by nearly 0.4°C from a supercooled state to +0.34°C (Fig. 6). From 21 January 2011 (day 67), the water temperature continued to rise even though the air temperature remained far below 0°C (Fig. 6). This is a sign that the ice outer shell is completely closed; the isolated water was no longer directly affected by the air temperature. It is important to mention that at the end of December a deep snow cover had already covered the stream between the foot of the ice cascade and the location of our water-level and temperature logger. Moreover, the only important snowfalls in January occurred on 22 January (day 68) with an accumulation of 12 mm w.e. and 26 January (day 72) with 26 mm w.e. We believe that the only condition responsible for this temperature change is that ice completely covered the waterfall.

Fig. 6. Air and water temperature evolution during the growing transition phase: from 1 January 2011 (day 47) to 10 February 2011 (day 87).

Because the ice cover is thin at the beginning of this transition stage, the ice volume accumulation remains substantial, but it rapidly decreases at the beginning of February. This asymptotic behavior was also observed by Montagnat and others (2010) on an ice cascade in the French Alps. These observations suggest a transition between a growth dynamics controlled by convective heat exchange between water and the atmosphere to a dynamics where the heat is removed from the flow through the ice cover.

On 12, 16 and 17 March 2011 the air temperature rose above 0°C and liquid precipitations induced an increase of the discharge above the base level. The effect was not visible on the lidar measurement, but we noticed the formation of multiple voids between the ice shell and the rock wall. The ice cascade surface was still intact, but the ice volume suffered a clear and drastic decrease (gray arrow, Fig. 5). Numerous small holes opened in the ice cascade (ice cascade pictures taken on 30 March and 10 April, Fig. 2). After these events, the air temperature fluctuated around the melting point and the water discharge increased. The ice cascade rapidly melted. Finally, the major rain event of 5 April 2011 resulted in a rapid increase in the stream discharge. This sequence of events resulted in a radical dismantling and collapse of the ice cover between 10 and 14 April 2011.

Modelling

The ice cascade growth-and-decay process is controlled by a series of complex and highly fluctuating factors. The atmospheric conditions (air temperature, radiation and evaporation) affect heat transfers in many directions: from the flowing water and the ice surface to the ambient air, and from the rock wall to the forming ice cascade. Additionally, water temperature, precipitations, discharge, streamflow velocity and their spatial distribution over the rock wall all play a role in the formation of the ice cascade. To better understand the role played by each of these variables, a thermodynamic model (model 1) was designed based on the crude assumption that the water is flowing over an ice cascade fixed to a vertical rock wall (Fig. 7). The different components of the energy balance are transferring heat to or from the water flow at a temperature close to the freezing point. As shown below, this model strongly overestimates the amount of ice produced over the cold season. Therefore, we developed a two-stage model (model 2) where the transition from a regime controlled by the air–water interaction to a regime where the isolation effect of the ice cover is taken into account. Graphic inspections and coefficients of determination (R ²) were used to assess model performances and determine good empirical constants (e.g. Gottlieb, 1980; Hock, 1999).

Fig. 7. Schematic illustration of an ice cascade growth-and-decay model (model 1).

Model 1

This model is based on energy budget and applies theoretical considerations inspired mainly by ice stalactite growth models (e.g. Makkonen, 1988; Maeno and others, 1994), river ice and icing growth models (e.g. Schohl and Ettema, 1986; Ashton, 1989; Hu and others, 1999), ice accretion accumulation models (e.g. Jones, 1996; Makkonen, 2000) and frazil-ice production models (e.g. Hanley and Michel, 1977; Ye and Doering, 2004). The various heat fluxes associated with water flowing down a vertical and simplified plane (the rock wall) are presented in Figure 7. At any time during the freezing or melting process the energy balance can be written as

(4)

where Q _cv and Q _evap are the sensible and latent heat flux, respectively, at the air–water interface, Q _rad is the radiative heat budget, Q _cc is the heat brought by conduction from the rock wall to the ice cascade, Q _w is the sensible heat carried by the flowing water and Q _ice is the latent heat released during the ice production process (phase change). In this model, we assume that the flowing water around the forming ice is at the freezing point.

Sensible and latent heat flux

As the water falls on and over the rock wall, air convection around the flow and the falling droplets evacuate or bring energy to the water. The convective heat transfer between the air and the water can be expressed as

(5)

where h _wa is a heat transfer coefficient that depends on meteorological conditions, especially the wind speed. The difficulty in determining Q _cv is, in fact, in determining h _wa. The option often taken is to use a bulk aerodynamics approach (e.g. Oke, 1987; Hock, 1998, 2005; Hu and others, 1999; Wagnon and others, 2003) where h _wa can be expressed as

(6)

in which ρ _a is the air density, C _a is the specific heat of the air under constant pressure, U is the wind speed and K _s is the bulk-exchange coefficient. Generally, K _s is determined by calculating the relative effect of air buoyancy on the mechanical forces of eddies over a horizontal surface. Such an approach is inapplicable in our case because the heat budget is calculated for water flowing over a vertical plane. We decided to opt for a simplified approach where K _s is adjusted empirically.

We also applied this approach to determine the latent heat flux (e.g. Makkonen, 1988; Maeno and others, 1994):

(7)

where L _e is the latent heat of evaporation at 0°C, C _w is the specific heat of water, P is the air pressure, R is the relative humidity and and are the saturation water-vapor pressures over water at T _a and T _f.

Because ρ _a, C _a, L _e, C _w and T _f are constant and T _a, U, P, R and are measurable variables, the only empirical value that remains is K _s. With this approach, we assume that h _wa varies (increases) linearly with wind speed.

Net all-wave radiation

The net radiation of the ice cascade surface is the difference between the incoming and outgoing radiation absorbed and emitted by the surface:

(8)

where G is the global (direct and diffuse) solar radiation, a is the albedo of ice and IR is the longwave radiation emitted by the ice surface. We are aware that Q _rad is probably slightly underestimated by our measurements as the spectral range of measurement of the pyranometer used in the study does not include the incoming longwave radiation. Furthermore, an albedo value for this type of ice is unavailable. The lack of literature on the subject suggests that such measurements have not yet been realized. The albedo of natural ice is highly variable: it can have a value of 0.1 for dirty glacier ice, 0.3 for newly formed sea ice, 0.4–0.5 for slush ice or freezing ice and even ∼0.9 for ice covered with fresh snow (e.g. Bolsenga, 1969; Oke, 1987; Hu and others, 1999; Hock, 2005). During the ice cascade growth and decay, the albedo of the ice–water system changes from a bluish ice formed around and over the dark waterfall surface to a nearly white and slightly translucent ice cascade that completely covers the rock wall (Figs 2 and 4). Consequently, we used a value of 0.5, which represents a qualitative estimation and a mean value for different types of natural ice.

Finally, IR can be written

(9)

where σ is the Stefan–Boltzmann constant and a is the radiation linearization constant (8.1 × 10⁷ K³). In Eqn (9), it is assumed that the emissivities of the ice cascade surface and the environment are both unity (Makkonen, 1988, 2000).

Rock-wall conductive heat flux

As shown in Figure 6, the conductive heat flux from the rock wall was taken into account:

(10)

where k _rw is the thermal conductivity of the rock. In this case, the rock wall is a highly fractured assembly of bedded shale, limestone and sandstone (Brisebois and Nadeau, 2003). The thermal conductivity of this type of rock varies from 1.2 W m⁻¹ K⁻¹ for a highly hydrated shale to 2.3 W m⁻¹ K⁻¹ for a dry mudstone–sandstone (Bloomer, 1981). It essentially depends on water saturation, rock type and bedding planes. In our case, we assume that the rock mass is highly hydrated and we choose a value of 1.5 W m⁻¹K⁻¹. T _rw corresponds to the mean annual temperature in northern Gaspésie (3.65°C) (Environnement Canada, 2000), T _f is the temperature of the rock–ice interface (0°C), X _rs is the rock–ice interface depth (0 m) and X _rw is the depth at which the rock temperature does not vary from the mean annual temperature (2.1 m). To determine this depth, we measured the rock-wall temperature from the surface to 80cm underneath the surface. After a 1 year measurement period, we were able to extrapolate the depth where no temperature change occurs during the whole year. At the end, a constant value of 2.6 W m ⁻² was calculated for the rock-wall conductive heat flux.

Energy flux from flowing water

The fifth term in the energy-balance equation represents the heat flux carried by the flowing water. This heat flux is very important because it also represents the volume of water available for the production of ice. It can be written

(11)

where C _w is the specific heat of water, D _in is the incoming discharge, T _wt and T _wb are, respectively, the water temperature at the top and bottom of the waterfall and S is the mean surface area covered by the ice cascade (or the waterfall) on the rock wall.

Energy flux released by the freezing water

When the temperature drops below 0°C and a sufficient amount of heat is lost through the environment, ice forms in the flowing water. The latent heat released during the freezing process is

(12)

where L _f is the latent heat of freezing (or melting) and dM/dt is the ice production rate per unit mass. The energy balance (Eqn (4)) can then be used to estimate the ice production:

(13)

where dM/dt is in kg s⁻¹. Finally, the volume rate production is given by

(14)

where ρ _i is the density of the ice. The ice density was determined from different ice structures (ice stalactites, ice accretions) sampled directly on the ice cascade and measured by tomodensitometry

(e.g. Hounsfield, 1973; Knoll, 1989; Boespflug and others, 1994; Calmels and Allard, 2004). A mean value of 880 kg m⁻³ was used for calculations. Finally, K _s remains the most critical adjustable variable in the model (Eqn (6)).

Model 1: results

Figure 8 shows different ice production curves with different values of K _s ranging from 0.0005 to 0.2. By giving K _s a high value, the overall ice production in the cascade follows the general shape of the FDH curve presented in Figure 5e. In this situation, the sensible heat flux and the latent heat flux are the dominant terms in the heat balance, which minimizes the effect of other terms used in the model. This agrees with what can be expected from the hydrometeorological data analysis where the air temperature evolution represented by the FDH appears to be the main variable controlling the ice cascade growth. Using a small K _s means reducing the effect of Q _cv and Q _evap in the ice production process, hence leading to a heat budget mainly controlled by the other terms, especially the sensible heat carried by the flowing water, Q _w. The increase in ice production after day 164 in this case (Fig. 8) is a spurious effect of _{Q w.} During this period, the air temperature, which remained above the melting point, transfers heat to the water, resulting in a water temperature higher at the bottom than at the top of the waterfall. This leads to a negative value of Q _w. In fact, it represents the energy transferred to the water by Q _cv and Q _evap that leaves the system along with the flowing water. This effect is hidden by the high value of Q _cv and Q _evap when using a high value of K _s.

Fig. 8. Calculated ice production (model 1) with different values of K _s: from 15 November 2010 (day 0) to 31 May 2011 (day 197).

When running the model, we realized that the linear effect of the wind on h _wa (Eqn (6)) was incoherent with our observations. During cold but calm days we observed that the ice cascade was actually forming very rapidly. In the model, a very low value of wind speed tends to lower the h _wa value near 0. The result is an underestimated value of Q _cv and Q _evap during these calm periods. To adjust this effect into the model, we were forced to increase the value of K _s. To best fit the changes in measured ice volume, a value of K _s between 0.01 and 0.2 must be used.

The amount of ice predicted with this model is higher than the actual ice volume (Figs 5 and 8). If the model is well balanced, then the ice production is well estimated. This implies that in reality the system must evacuate most of the ice produced. Only a small fraction of the ice remains fixed to the rock wall. This finding is not surprising and is in close agreement with the production of frazil ice reported in fluvial systems (e.g. Martin, 1981; Osterkamp and Gosink, 1983; Turcotte and Morse, 2011). Frazil ice is produced in rapids with highly turbulent flow and carried away by the current at the foot of the riffles. Frazil-ice concentration can reach 10⁷ ice discs per m³ of water and the production rate can reach 100 × 10⁶ m³ d⁻¹ in a large open-water river (Michel, 1971; Martin, 1981; Osterkamp and Gosink, 1983). Sometimes the ice pellets are fixed to the gravel bed as anchor ice, but they usually form floating frazil flocs, frazil pans or slush, which accumulate against booms or static ice covers in low-velocity zones. We observed during low-temperature periods that the water temperature at the bottom of the waterfall was in a supercooled state and that a large amount of frazil ice was leaving the system. Moreover, by comparing the average amount of ice formed between lidar measurements with the discharge, we calculate that no more than 0.5% of the discharge is retained in the system to form the ice cascade (before 21 January 2011). After the complete closure of the ice cover over the rock wall (after 21 January 2011), we found that 0–5% of the discharge is retained in the system to form the ice cascade. When we consider the calculated value during the initial growth stage (before 21 January 2011), the energy evacuating the system might be sufficiently high to cause the entire flow to freeze. This behavior may seem unrealistic, but it is in agreement with the observations: a large amount of the supercooled slush-like moving water flows away. After the ice completely covers the rock wall, the discharge becomes very low and the flowing water is isolated from the environment by the thickening ice wall on the cascade face. To adapt the model to those issues, an empirical variable representing the percentage of ice formed (%IF) and remaining in the system must be used.

Model 2

As mentioned above, during the initial stage of the growing process the water is constantly in contact with the air and the ice only forms at the edges of the water flow. From 21 January 2011 (day 67), an ice cover forms over the whole rock-wall surface and the growth dynamics changes (Figs 2 and 6). The ice cover isolates the flowing water from the outer environment. A convective sensible heat transfer at the water–air interface does not apply to this new situation. To improve the efficiency of the model, we now propose a two-stage model. Figure 9 presents a schematic illustration of this second model.

Model 2: stage 1

During the first stage (days 0–66), the flowing water is directly in contact with the atmosphere (Fig. 9). The same heat budget used in model 1 is applied (Eqn (4)), but this time we modified the calculation of the heat transfer coefficient:

Fig. 9. Schematic illustration of a two-stage ice cascade growthand-decay model (model 2).

(15)

where U _eff represents the effective contribution of the wind:

(16)

where U _d is the droplet (water flow) mean falling speed. The ‘aeration effect’ that occurs during the fall of the water increases the area of the air–water interface because the water flows consist of sprays, droplets and broken streams (e.g. Chanson and Cummings, 1996; Chanson, 1997; Zhang and others, 2001; Cheng and others, 2004). The free-falling speed of a raindrop ranges between 1 and 9 m s⁻¹ depending on the size of the droplet (e.g. Gunn and Kinzer, 1949; Atlas and others, 1973). After 1 m of drop, a large droplet (0.6 cm) can reach the maximum speed of free fall. Since the droplet size and the proportion of the water flow that fragments and falls freely are unknown, we assumed a mean value of 5 m s⁻¹ for a 45 m waterfall. Such a constant will reduce the tendency to minimize the role of air convection around the water flow during calm days. Indeed, the wind generated during the free fall of the water flow favors convective heat exchanges with the air. Adding a mean value of the droplet (water flow) falling speed results in a minimum amount of wind speed effect in the calculation of h _wa. To adapt the model to this improvement, a value of 0.048 was given to K _s to best fit the measured ice accumulation trend (graphic inspection) and obtain the highest coefficient of determination (R ²) between the measured and calculated ice volumes (Fig. 10a).

Fig. 10. (a) Calculated and measured ice cascade volume evolution. (b) Calculated versus measured ice cascade volume. The thin dotted line is the model 1 tendency with a %IF of 0.2% (K _s = 0.2). The thin black line is the first stage of the second model from day 0 to day 66 (15 November 2010 to 21 January 2011) (%IF = 0.2%; K _s = 0.048). The gray line is the continuation of the second model’s first stage after day 66 (%IF = 0.2%). The thick black line is the second stage of the second model from day 66 to the end of the modelling (21 January 2011 to 31 May 2011) (%IF = 1.5%; K _s = 0.9). The dashed gray arrow represents a more probable ice volume evolution.

The other variables were not modified. However, we added an empirical constant representing the percentage of ice formed and remaining in the system (%IF = 0.002 or 0.2%) to adjust the model to the ice volume measurement:

(17)

Model 2: stage 2

During the second stage (from day 67 until the end of the observations, day 197), the ice covers the whole rock wall and different physical processes take place. The ice isolates the water flow from the atmosphere. In the model, all the heat is evacuated from the isolated water flow even if occasional exfiltration over the ice cascade is observed and involved in the growing process. The heat budget can now be written as

(18)

As with the model proposed by Montagnat and others (2010), we suggest a conductive heat transfer through the ice to the water, keeping in mind that the outside surface temperature of the ice cover (T _is) tends to T _a and the temperature at the ice–water interface is T _f. In this case, we assume a linear gradient of the temperature within the ice cover and the conductive heat transfer flux through the ice (Q _ic) can be described by (e.g. Incropera and DeWitt, 2007)

(19)

where I _t is the ice-cover thickness and K _ice is the thermal conductivity of ice (2.22 W m⁻¹ K⁻¹ at 0°C). We calculate I _t at each time-step by dividing the calculated ice volume by the mean surface area covered by the ice cascade (S). Q _ic takes effect when ice formation occurs (T _a <0°C). When T _a is above the melting point of ice, air convection at the ice–air interface favors heat transfer to the ice cascade. The ice surface temperature (T _is) is at T _f and Q _cv can be written

(20)

where h _ia is the heat transfer coefficient at the ice–air interface, again calculated as a bulk formulation:

(21)

where K _s is the bulk-exchange coefficient between the air and the ice surface. An empirical value of 0.09 returns the best fit with the measured ice volume values. Because there is now no more direct interaction between the water and the outside environment, Q _evap is replaced by a sublimation heat flux:

where L _s is the latent heat of sublimation at 0°C, C _i is the specific heat of ice and is the saturation water-vapor pressure over the ice at T _is. Because we assume that T _is = T _a when the temperature is below 0°C, Q _sub remains null during the freezing period. When the temperature is above 0°C, T _is = T _f and sublimation occurs.

Finally, we consider that the other fluxes remain identical as in the first stage. An empirical variable representing the amount of ice retained in the system and participating in the ice cascade growth process is required (%IF). As mentioned previously, no more than 5% of the water mass flux is retained as ice in the system during this period. Hence an empirical value of 0.015 (1.5%) was used for %IF in the ice volume calculation (same forms as in Eqn (17)).

Model 2: results

Figure 10 presents the calculated ice cascade volume evolution given by model 2. In Figure 10, model 1 is also presented, but with the use of the same %IF as in the first stage of model 2 (0.2%) and with a K _s value of 0.2. Figure 11 shows the normalized evolution of the different fluxes for the entire duration of the model run and Figure 12 presents representative examples of an energy balance during freezing and melting periods for the two different stages of the model.

Fig. 11. Comparison of normalized heat fluxes during the first and second stages. Each normalized heat flux is the ratio between fluxes and the sum of fluxes of the same sign (positive or negative).

Fig. 12. Examples of 24 hour energy balance during freezing and melting periods for the two-stage modelling. 13 December 2010 was a partly cloudy day with an average air temperature of 9.9°C, a minimum of 1.9°C and a maximum of 16.7°C. 7 January 2011 was a partly sunny day with an average air temperature of −10.3°C, a minimum of −15.2°C and a maximumm of −4.3°C. 27 February 2011 was a sunny day with an average air temperature of −16.0°C, a minimum of −19.1 °C and a maximum of −13.9°C. 27 March 2011 was a sunny day with an average air temperature of −6.3°C, a minimum of −8.0°C and a maximum of −4.5°C. 10 April 2011 was a partly sunny day with an average air temperature of 3.0°C, a minimum of −1.8°C and a maximum of 6.2°C. 16 April 2011 was a sunny day with an average air temperature of −2.2°C, a minimum of −7.5°C and a maximum of 1.4°C.

The results highlight differences between the two models. First of all, the use of a constant representing the water free-falling speed (U _d) during the early stage of formation has forced the readjustment of K _s from 0.2 to 0.048. The resulting value is now closer to values used to calculate energy exchanges over a melting glacier (e.g. Hay and Fitzharris, 1988). Except for the ice volume measurement of 10 December 2010, the calculated values are very close to measured volumes (Fig. 10a and b). After the complete insulation of the water flow by the ice cover on 21 January 2011, both models (model 1 and model 2: stage 1) overestimate the ice cascade volume. At this stage of formation, the transition between a regime ruled by air convection around the water flow and a regime ruled by conductive heat transfer through the ice cover (model 2: stage 2) has led to a better representation of the bulk evolution of the ice volume. As reported by Montagnat and others (2010), ice cascades show an asymptotic growth behavior when they reach a given thickness. In our case, the ice reached a mean thickness value of 8 cm before such behavior was observed (∼day 80). The model also depicts this tendency. Following this slow-growing stage, a small increase of the discharge after a 15 mm rain event (between day 116 and day 120) accelerated the melting of the ice cover (Fig. 5d and e). The model shows a small decrease during this event, but is not able to predict the rapid decrease in ice volume associated with the partial collapse of the ice structure (Fig. 10a).

During stage 1, Q _cv is the dominant flux controlling the ice cascade growth and melt process (Figs 11 and 12). Evaporation (Q _evap) also carries heat out of the water, but its contribution remains small. The diminishing energy carried by the flowing water (Q _w) is of great importance at the beginning of the growth process. If the discharge remains high, the ice cascade will only form on the edges of the flow and the other heat fluxes will not be efficient enough to freeze the water in the main flow path. The radiation heat flux, mostly the solar diffuse radiation, carries a little heat to the system but it seems to contribute only when the air temperature is close to 0°C. The amount of heat coming from the rock wall appears negligible in comparison with the other fluxes.

The second stage of model 2 follows the asymptotic tendency of the measured ice volume. The use of a conductive heat transfer flux into the ice (Q _ic) yields a better approximation of the heat exchange between the water and the outside environment. Q _ic and Q _cv (at the air– ice interface) remain energy fluxes of primary importance (Figs 11 and 12). Right after the transition phase, the discharge is still at its lowest level and the diffuse solar radiation received at the ice surface is very low. The ice cascade can pursue its growing process rapidly. Even if the heat conduction into the ice is low, the ice is thin and allows a relatively good heat exchange between the thin liquid film and the outside cold air. But as the ice cascade gets thicker, the heat exchange (Q _ic) also decreases. During this period, the ice cascade is in a delicate equilibrium state. The growth of the ice cascade is now also controlled by the radiative energy balance (Q _rad), but more specifically by the increasing contribution of direct solar radiation from the end of February (Figs 11 and 12). Even if the direct solar radiation reaches the surface with a high incident angle, the radiative energy budget tends towards positive values favoring heat accumulation at the surface of the ice. This leads to a decrease in growth rate during cold days, but favors melting when air temperature tends toward 0°C (Fig. 12). The rising discharge in mid-March increases the amount of sensible heat carried by the flowing film of water (Q _w). Our observations show that the ice cascade melts rapidly at the ice–rock interface. As mentioned earlier, we observed massive voids opening between the ice and the rock wall on 30 March and 10 April 2011 (Fig. 2). Actually, the increase in discharge began on 17 March 2011 after a 12.5 mm rain event (Fig. 5d). But the ice cascade melts and collapses rapidly only after the significant increase of the discharge on 11 April 2011 (Figs 10a and 12). The ice volume measured on 15 April 2010 was composed of small residual ice structures standing on the edges of the main water flow. The model trend cannot reproduce the thermoerosive action of the flowing water and the collapsing process occurring during this fast melting period. At this point, model 2, stage 2 becomes incompatible with the ice formation and melting process. During this last period, some ice may form around the main flow during a cold spring night but rapidly melts over the following days. The growth and melt processes then return to a regime that is more compatible with the first stage of model 2.

Model 2: sensitivity analysis

Three constants (K _s, U _d and α) were tested to verify the sensitivity of the model to changes in parameter values. The results are presented in Figure 13. Increasing or lowering the value of the bulk-exchange coefficient (K _s) by 30% affects the magnitude of the ice volume evolution (Fig. 13a). This leads to a difference in the percentage of ice fixed to the rock wall (%IF) of +0.2%. The same kind of effect is observed when changing the value of the mean free-falling speed of the droplets (or water flow) (U _d) down the rock wall (Fig. 13b). Lowering U _d from 5 to 3 m s⁻¹ leads to an increase of 0.13% of retained ice in the system (%IF), and increasing U _d to 7 m s⁻¹ leads to a decrease of %IF of 0.2%. Because these two factors (K _s and U _d) affect the calculation of the heat exchange coefficient (h _wa) and therefore the convective heat flux (Q _cv), their revisions mostly change the ice growth-anddecay curve during the first stage of formation when Q _cv is the main heat transfer flux.

Fig. 13. Calculated and measured ice cascade volume evolution for different values of (a) the bulk-exchange coefficient (K _s), (b) the water flow or droplet free-falling speed (U _d) and (c) the albedo (α).

Since the radiative heat budget (Q _rad) has a significant impact during the second stage of formation, it is relevant to evaluate the sensitivity to the ice surface albedo (α). To verify the effect of α on the ice cascade volume evolution, higher and lower values were tested (0.3 and 0.7) (Fig. 13c). These modifications do not cause any changes during the first stage of formation. However, during the second stage, a lower α value leads to an underestimation of the ice volume fixed to the rock wall (%IF) of nearly 3%. A higher α value leads to an overestimation of %IF by 3%. After complete formation of the ice cover over the rock wall, an α value of 0.3 is clearly insufficient considering the bluish, nearly white color of the ice cascade (Figs 2 and 4). An albedo between 0.5 and 0.7 is certainly much closer to reality.

Even if modification of these three parameters changes the amplitude (increase or decrease) of the ice volume evolution curve, it does not affect its shape and general trend. The amplitude of the variations is not critical in this study because it is adjustable through %IF. However, it is essential that the shape of the curve is not modified in order to identify the predominance of the different heat transfer processes during the various growth phases.

Discussion and Conclusions

Modelling the ice cascade volume evolution has revealed the importance of the different hydrometeorological variables. First, the importance of the convective heat transfer (Q _cv) in the initial growing stage was confirmed. The temperature gradient between the flowing water near 0°C and the air temperature supports the freezing degree-day (FDH) calculation. Even though the FDH curve shows a significant lag compared with the measured ice volume, it shows the same general trend. It might eventually be used to model other kinds of rock-wall icings or to predict their collapse (Gauthier 2008; Gauthier and others, 2012).

The linear dependency between the convective heat transfer coefficient (h _wa) and the wind speed (U) is of utmost importance. Unlike many thermodynamic models developed to simulate the growth of ice stalactites (Makkonen, 1988; Maeno and others, 1994), icing and river ice formation (Schohl and Ettema, 1986, 1990; Ashton, 1989; Hu and others, 1999) or glacier melt models (e.g. Munro, 1989; Hock, 1998, 2005; Wagnon and others, 2003), the present model is adjusted by introducing a variable representing the mean free-falling speed of the water flow (U _d). This adjustment provides a better representation of the physical processes involved in the formation of ice accretions around the main water flow after the freezing of water spray, a general process similar to freezing rain accumulations (Jones, 1996; Makkonen, 2000; Gauthier, 2008; Gauthier and others, 2012). Furthermore, since a very high percentage of the ice produced during the free fall of the water flow is evacuated at the base of the frozen waterfall, an ad hoc constant representing the percentage of ice fixed to the rock wall and contributing to the ice cascade volume evolution (%IF) had to be added to the model. The values of %IF were fixed empirically from the ratio between the discharge and the mean amount of ice produced between each lidar measurement. This constant was proven to be stable and the model is robust to parameter changes.

During the initial growth stage, the model also reveals the relative importance of the latent heat flux (Q _evap). This flux, mainly controlled by the relative humidity, is significant when the air temperature is slightly below zero and the relative humidity is low. Under these conditions, Q _evap greatly favors ice formation and becomes the main heat flux. The low mean annual air temperature in northern Gaspésie leads to a low conductive heat flux (Q _cc) emanating from the rock wall. This flux remains low and trivial in our case, but might play a more relevant role in warmer climates such as the European Alps or western America where ice forms sporadically during winter. The diminishing water flow at the beginning of winter favors the growth of the ice cascade. Conversely, in spring, the increasing discharge becomes a major heat source during the ice cascade melting process. These observations are consistent with the results obtained by Maeno and others (1994). They showed that ice stalactites tend to grow faster under a low-discharge regime. Furthermore, the relevance of the sensible heat flux carried by the flowing water (Q _w) may not be significant under a low-discharge regime as in the case of growth on ice walls, which forms after the gradual freezing of groundwater seepage on a cliff.

Finally, one of our major findings is the critical role of the radiation budget (Q _rad) after the isolation of the water flow by a complete ice cover. While the radiative heat transfer (Q _rad) intake is insignificant during the initial growth stage, it becomes a major source of heat during the second stage. This transition between a system controlled by air convection (Q _cv) and one controlled by a conductive heat transfer flux through the ice (Q _ice) changes the physical growth process. The growth rate is reduced during the day when direct solar radiation is received at the ice surface. As the solar radiation increases in spring, the growing and melting process becomes highly controlled by this intake.

To our knowledge, this paper is the first to describe an energy-balance model for this kind of ice formation. Moreover, the two-stage model proves to be relevant for simulating the volumetric variations of the ice cascade. However, further studies need to be conducted to evaluate the accuracy of the empirical variables and constant values used in existing models: the ice cascade albedo, the heat exchange coefficient expressed in terms of the bulk-exchange coefficient, and the percentage of ice produced and evacuated from the system, to mention just a few. Furthermore, the models need to be validated in different environmental conditions (e.g. slope aspect, low-discharge regime, other climate regimes, etc.).