Development of a river ice jam by a combined heat loss and hydraulic model

This paper discusses and shows the heat loss theory and the hydraulic theory for the analysis of the development of wide channel ice jams. The heat loss theory has been used in Iceland for a long time, while the hydraulic theory largely follows the classical ice-jam build-up theories used in known CFD models. The results are combined in a 5 new method to calculate the maximum thickness and the extent of an ice jam. The results compare favorably to the HEC-RAS model for the development of a very large ice jam in Thjorsa River in Iceland. They are also in good agreement with historical data, collected where a hydroelectric dam project, Urridafoss, is being planned in the Thjorsa River. 10


Introduction
Ice jams are among the most dramatic natural events that occur in a river.Understanding of ice jam formation and break up is very important in river engineering, especially dams and water diversion works.As a rule, water levels are greatly increased when an ice jam forms in a river section.Ice jams often lead to potentially unwanted situations for the human activities along the banks of the river.Other major difficulties are reduced flow during the formation of an ice jam and surges of water and ice fragments during break-ups.Uzuner and Kennedy (1976) developed the hydraulic equation system and in Beltaos (1993) a model is applied to three case studies of ice jam events and the results found to compare well with observations.The various model coefficients fall within the expected ranges, with only one exception.A thorough description of the formation and evolution of ice jams is given in Beltaos (2008).It builds on 1990 state of the art, Beltaos (1995) and a large number of publications Correspondence to: J. Eliasson (jonase@hi.is)exist from other authors and institutions about various details as well, see e.g.Beltaos and Burrell (2006) on temperature differences between water and ice jam.Here, the Cold Regions Research and Engineering Laboratory (CRREL) is an important source.It is therefore a reason to believe that CFD models can handle the hydraulic behavior of ice jams correctly, this is demonstrated in Beltaos (2008).In this paper we do not present CFD results, but one-dimensional (e.g.RIVICE, RIVJAM, ICEJAM, HEC-RAS, etc.) and twodimensional (2-D; www.river2d.ca)public-domain ice-jam models are now available Morse and Hicks (2005).Instead we describe the force balance that is used to predict the thickness and shape of the freeze-up jams (Grondal, 2003) and combine it with a heat loss model that has been known for considerable time.The heat loss model can only predict formation of ice mass in the river, and the force balance model can only describe the ice jam thickness that is in equilibrium with the river flow.It is shown that these models can be combined through a single equation.The results are compared with field data from Urridafoss (Fig. 1) in Thjorsa River in Southern Iceland.

Freeze-up ice jam at Urridafoss in Thjorsa
Thjorsa River originates at Hofsjokull glacier in Central Iceland and flows to the South-West approximately 230 km where it discharges into the North Atlantic Ocean, see Fig. 1.
The river system has a large hydropower potential that has been developed quite extensively in the last four decades, but the development has been concentrated in the upper reaches.The freeze-up jam under discussion in this article forms in the relatively flat section just downstream Urridafoss waterfall, as a consequence of frazil ice production in the approximately 50 km long river section downstream of the power plant at Burfell (Fig. 1).The Urridafoss ice jam is formed almost every winter.It typically extends through the lower part of the Urridafoss gorge down to the flood plain, in all a distance of about 3-4 km.The width of the jam in the gorge is approximately 100-400 m, and expands to roughly 700 m on the flood plain.Water level increases up to about 18 m have been observed (Rist, 1962).The formation and evolution of the jam was first described by Rist (1962).Further investigations of the ice conditions are in Eliasson and Gröndal (2006).Some of these investigations are planned to obtain the necessary design data for a dam in the Thjorsa River at the Urridafoss site and a hydroelectric power plant associated to it, but there is also an investigation by Gröndal (2003) reported, where the HEC-RAS model, Brunner (2001) was applied.It proved quite easy to construct an ice jam model in HEC-RAS and get results that compare very well to observations when a roughness coefficient of n c =0.025-0.035was used

The heat loss model
There is considerable experience in heat loss calculations, Carstens (1970) and Freysteinsson and Benediktsson (1994) report both experimental results and field observations.In the heat loss model that was used to estimate the volume of the Urridafoss ice jam two equations are solved, namely a heat transport equation and an ice transport equation: and distance along longitudinal axis T water temperature in cross section C ice concentration in cross section, m 3 /m 3 V flow velocity S heat loss from water column y depth of flow ρ w density of water ρ i density of ice c p specific heat of water L latent heat of fusion of water According to Eq. ( 1) the temperature of the water decreases when there is net heat loss from the water surface.As soon as the temperature of the water has dropped to the freezing point of the water, the temperature decrease stops.Instead, frazil ice begins to form at the rate corresponding to the heat loss, according to Eq. ( 2).Thus, by solving Eq. ( 1) and ( 2) in combination, one can find the total ice produced in a river section, given that the heat loss, S, can be determined.Heat loss from the river is governed by 1. Rate of heat exchange with the atmosphere 2. Rate of heat exchange with the river bed 3. Heat transfer via groundwater inflow

Frictional heating
In Thjorsa, term 1 (Rate of heat exchange with the atmosphere) (Carstens, 1970b) is the dominating one, and the other terms can be neglected without serious error.Net rate of heat exchange with the atmosphere in W/m 2 , is a sum of the effects of terrestrial or long wave radiation, heat transfer due to evaporation or condensation of water, sensible heat transfer due to convection and heat transfer due to precipitation, minus the effects of incoming solar or short wave radiation.Grondal (2003) discusses methods that can be used to quantify heat loss caused by these processes.
Figure 2 shows the result of the calculations of ice volume in the winters 1958/59 to 1963/64 and 1998/99 to 2001/02.Calculations are done for average winter flow, but the actual river discharge does not affect the ice production directly, it is a function of the size of open river surface and the weather parameters.The size of the open water area is taken as a constant value obtained from measurements made during winter conditions when shore fast ice occupies the borders of the river.According to the heat loss model about 35 to 40 mil.m 3 of solid ice are produced on the average each winter.In mid winter accumulated volume is often about 20 mil.m 3 .At this time there is often a large ice jam at Urridafoss (Fig. 1). Figure 3 gives an idea how this production is distributed throughout the winter.

Forces in an ice jam, hydraulic theory
The external forces acting on the jam arise from the following factors: Friction between ice cover and flowing water, backwater pressure, the longitudinal component of the ice and pore water weight.They are balanced by internal normal stresses and boundary shear stresses at the riverbanks.There are slightly different methods to formulate this so a brief description of the method applied will be given here.
As the jam lengthens upstream and thickens, the forces acting on the jam increase, until internal stresses in the jam become too large.At that point the ice jam lengthening process stops, which may lead to shoving, i.e. consolidation and thickening of the jam.Broadly speaking, this process then repeats while the supply of ice from the river upstream continues.Here we follow Beltaos (1995) and Uzuner and Kennedy (1976), they derive the one dimensional force balance equation for floating ice jams.Useful explanatory figures are Fig. 5.42 in Ashto (1986) andFigs. 3.13 and4.4 in Beltaos (1995).Their analysis leads to the following equation for the thickness of the jam: (3) slope of water surface m/m K x =tan(π/4+φ/2) equivalent Rankine passive pressure coefficient k 0 =tan φ angle of internal friction in jam k 1 coefficient of lateral thrust

C i cohesion in jam
shear stress between water and underside of jam density of ice ρ w density of water g gravity constant α water surface slope angel in radians Now it is assumed that the cohesion C i can be neglected, Eq. ( 3) then reduces to: All the above mentioned authors use quasi-steady momentum and energy equations for the flow as local acceleration in natural rivers is very low because of slow changes in the flow.For steady state flow, the energy equation is used to calculate the water surface profile in the jam.The reader is referred to Fig. 5.42 in Ashton (1986) 2.
Y 1 , Y 2 water depth at two cross sections 1 and 2 Z 1 , Z 2 elevation to channel bottom V 1 , V 2 average velocities α 1 , α 2 velocity weighting coefficients H energy head loss 2 Ice jam thickness and extent

Properties of the jam thickness equation
When investigating local behavior of h it is natural to assume that convective acceleration plays a minor role compared to gravity so changes in velocity head can be neglected.This makes the friction slope equal to the slope of the water level inside the jam.The water level relation becomes

Maximum jam thickness
When h and y are constant in the variable x, S f =S 0 .Now dh/dx can be zero for two values of h, found by solving (4) after inserting Eq. ( 6) and putting the left side to zero.The resulting quadratic equation has two roots, one negative but the other one is positive This h m is the maximum thickness the jam can reach.Similar result was obtained by Beltaos (1995).In Eq. ( 7) y may be calculated from the Manning equation using Note that Eqs. ( 6)-( 8) assume internal strength on the ice jam to be balanced by water drag, gravity, and bank resistance.
Ice jams therefore move during high flow period but sit still on the banks at low flow periods.As the strength parameters are not time dependent, Q in Eq. ( 8) should therefore be a little higher than the average in the particular period to give a correct picture of the development of the jam.

Change of slope
Equation ( 7) reveals that the h m is inversely proportional to S 0 .The quantity 1/a may be regarded as the length scale of the jam.When we have a slope change from a large S 01 to a small S 02 this length scale is reduced and with it h m .Upstream of the point of slope change we will have an ice jam with increasing thickness in the streamwise direction, h approaching h m1 .Downstream of the point of slope change the maximum thickness will be h m2 <h m1 .Figure 4 shows this development and it will be discussed in Sect.3.4.
3 Jam volume and length

Jam length
If we define K y =s i +y/2h it may be argued that K y is of the order one in thick jams.We use this approximation to put K y constant, insert Eq. ( 7) in ( 4) and get: One may notice that Eq. ( 9) produces almost the same maximum as the more accurate Eqs. ( 4) and ( 7), as long as the assumption K y is of order one holds.Equation ( 9) contains a new constant Equation ( 9) may be integrated Equation ( 11) is valid above a point of slope change.As x is a streamwise coordinate, we see that the jam thickens in the direction of the flow but never quite reaches the maximum value h m1 .If a jam, still under development, extends a distance L below the zero thickness point Eq. ( 11) gives the thickness when x=L is inserted.This is the largest thickness of the jam if it is still under development and the volume in the jam corresponds to the accumulated ice volume produced upstream of x=L.In an ice jam where there are no sudden changes in the channel parameters (slope S 0 or width B) Eq. ( 11) thus combines the heat loss and the hydraulic theory into one equation (Sect.3.3).It can also be used in a piecewise constant channels.
But ice jams normally occur where we have sudden changes in the channel parameters so this situation is considered in the next sections.

Change of width and slope
In Eq. ( 7) change of the width of the river channel, B, has the same effect as change of slope.Large changes in width do, however, usually bring larger changes in water profile than mere changes in slope.Now consider a river profile that suddenly changes from a large slope with maximum jam thickness h m1 to a smaller one with maximum thickness h m2 <h m1 .If actual jam thickness in the slope change point x=0 is h m2 <h L <h m1 , then the thicker upstream jam will be pushed beyond the slope change point x = 0 into the region with the lower maximum jam thickness h m2 and we will have below the slope change point when the ice jam is fully developed.Care must be taken in using Eq. ( 12) as the condition of low convective acceleration may very well not be fulfilled.This condition may very well hold for gradually funneling river channels, but not for abrupt channel changes, e.g. at the end of a gorge or down a waterfall, see Fig. 4. Here, ice sludge is being carried down the waterfall, below it a jam thick enough to drown the waterfall can build up.The ice jam will sit on the bottom until the water level inside the jam is high enough to lift it up, here we have an ice jam flooding situation with flooding levels that will increase until the waterfall is submerged and the ice jam build up can continue in the upstream reach.Provided ice production continues, the upstream jam will build up until it approaches maximum thickness h m1 .The length of this ice jam can be estimated using Eq. ( 11).Below the waterfall the situation is more complicated.The very thick jam is floating on the flood water and sitting on the bottom instead of being supported by the river banks only but Eq. (6) will still be valid for the section of the jam where hydrodynamic forces of the flowing water and internal forces in the dam are in balance.We can therefore consider Eq. ( 12) valid with the thickness of the jam just below the waterfall as the upstream boundary value h L , and the maximum thickness of the downstream section h m2 as the downstream boundary value where the force balance Eq. ( 4) is again active.In between there may very well be a different length scale in the exponential variation between the two values.This is discussed in Sect.3.7.

Jam volume
Equations ( 11) and ( 12) make it possible to estimate the total volume of the jam, Eq. ( 11) is integrated over the reach L and the average h a found, now we have As before L is the reach of the jam upstream of the point of slope change and h m1 the maximum thickness of the jam, the last approximation (preceded by ≈) is valid for very large L, that is fully developed jams.This remarkably simple estimate is based on that B and S 0 do not vary so much that the integration of Eq. ( 9) is seriously affected.If they do piecewise integration along the channel may still be possible.
To complete the volume estimate for situations like Fig. 4, L is found by successive approximation using Eq. ( 11) for the reach upstream of the sudden slope change (waterfall) and Eq. ( 13) for the channel reach downstream of it.We must begin by estimating how much ice volume, M * jam there is below the reach L. Ice jams form in the same place from year to year so the start of the downstream reach is known and the volume up to the point x=L can be estimated using h L ≈h m1 as the first guess.When M j is estimated from heat loss calculations M L jam =M j −M * jam and now an L value that satisfies Eq. ( 13) must be found.This is inserted for x in Eq. ( 11), h=h L calculated, that used in Eq. ( 12) to find new M * jam and the procedure continued until the approximation is completed to the sufficient degree of accuracy, which of course depends upon the accuracy of the original data.
This use of Eqs. ( 11)-( 13) combines the two theories, the heat loss theory for calculating volume of ice production, and the hydraulic theory for ice jam thickness for situations like the one in Fig. 4. The heat loss theory gives no information on jam thickness and the hydraulic theory gives no information on ice production.This combination is new theory that provides both.

Flooding because of ice jam building
In theory, the flood from an ice jam can be as high as the water level inside an ice jam of maximum height.The majority of the ice jam thickness will be below the water level, so it is on the safe side to estimate the maximum flood equal to h m in Eq. ( 7) above ice free water level in the river as the ice jam does not get thicker than that.

Building a dam in an ice jam river
When a dam is to be built in a river reach where frazil ice formation and ice jam building is known to take place, it is necessary to make the dam high enough so the water level inside the dam does not reach over the crest in the jam flood.The dam must thus be higher than the maximum thickness for the pond inflow channel h m , if expected ice production is large enough to build ice jams up to that level.Otherwise, Eqs. ( 11)-( 13) have to be used as indicated in Sect.3.3.Rist (1962) and calculated in meters.
Place in Fig. 4 Observ.Eq. ( 11) Upstream of waterfall, max. 9 9.4 Downstream of waterfall, max.18 17.4 4 km downstream, average 8 N/A 4 km downstream, maximum 12 N/A 3.6 Jam measurements in the years 1956-1962 In the winters 1954-1959 there was a great ice jam build-up, and also in the winter 1961/62.The maximum extent of the jam at Urridafoss (Fig. 4) was measured and reported by Rist (1962) Fig. 15.The biggest jams are in December 1958 and February 1961.The average difference in the thickness of these two jams is under 1 m, so bearing in mind the uneven surface of a frozen ice jam these two events produce identical jams, as would be expected from the accumulated ice production in the winters 1958/1959 and 1960/1961.They follow each other closely in the period from mid December to mid January on Fig. 2. Their surface profile, reported by Rist (1962) is shown on Fig. 4. Maximum thickness reported is in Table 1.
In Table 2 are shown the elevation measurements of the large, almost identical, jams in December 1958 and January 1961 (column 2 Jams).These two are still the largest that have been reported.It must be stressed, that the fact that these two are identical does not prove that ice jams in two different years, but at the same location formed by same accumulated amount of ice production are necessarily identical.Both flow discharge and periods with temperatures above freezing have their say.Column x is the distance in Eqs. ( 11) and ( 12), The theory (Sim values in Table 2) are calculated using piecewise integration (Sect.3.3) with actual S 0 values represented by the "River bed" line in Fig. 4 and river discharge Q=300 m 3 /s.The sensitivity of this figure is however small.A double discharge (600) would change the h m figure in distance 21.4 in Table 2 from 9.4 to 10.1 and have no other effect.The effect of division by two is even smaller.
The observations compare very well with theory.However, there are two artificial B values in the table indicated in italics, using these values changes the L scale shown in Table 2.The ice free width is 300-500 m in the respective river bed sections.There is no observed justification for this other than with ice free widths of the river bed the jam height will be 3 and 5 m too low.This is discussed in the next section.

Discussion
Inspection of Eq. ( 9) reveals two dimensionless parameters for the development of the ice jam.
Here L (not to be confused with L in Eqs. ( 12) and ( 13)), previously called 1/a in Eq. ( 7), is the length scale of the problem.When the two constants are equal for two different jams they are dynamically similar, i.e. a scale model of each other.The constants contain the width, the bed slope and the coefficients for the mechanical strength of the ice jam.They are the natural dimensionless groups to use in dimensional analysis.
In view of the fact that convective acceleration is neglected in the development of Eq. ( 9) this is the result one would have expected a priori.No objections to the use of Eqs. ( 11)-( 13) can be found in the composition of these dimensionless coefficients.
As the length scale appears as B/L and not in other context, we see that changes in B have the same effect as changes in the length scale.Narrow rivers (small B) can therefore be scale models of wide rivers (large B), provided that the other parameters in the coefficients have values that make the coefficients equal for model and prototype.
A partially floating ice jam, i.e. an ice jam pushed down a waterfall or accumulated around a point of sudden slope change, will partially sit on the shallow bottom near the banks with an effective width of the mid-channel water flow considerably smaller than the ice free width.When ice jams start thawing this middle floating section usually disappears first and exposes the effective middle width section for some time Rist (1962).
The equation system Eqs.( 9)-( 13) thus give a realistic picture of the build-up phase of an ice jam.It is still a problem however, what happens in the break-up phase.Ice jams in Iceland can be a product of repeated weather periods with frost and thaw.As may be seen in Fig. 3 there are many periods of thaw in between the periods of frazil ice run in one winter.Jasek (2003) states, that the interaction between the ice mechanics and unsteady flow leads to results that are often unpredictable with open water unsteady flow models.He also points out considerable differences of opinion on the degree of significance of this water-ice interaction.His conclusions lead to that considerable more experience has to be gained in research and analysis before ice jams resulting from complicated weather history as Fig. 3 presents, can be effectively predicted.A support for this may be seen by comparing the ice discharge figures in Fig. 3

Conclusions
The ice production model combined with solving the force balance equation can be used to predict the size of an ice jam, given that the parameters that appear in the force balance equation can be estimated.In the analysis at hand, assumptions were made that allowed for a relatively simple solution, but nonetheless a reasonably accurate result emerged.By using the heat loss theory to calculate the expected ice mass in an ice jam, Eq. ( 12) can be used to find the thickness and extent of a jam that corresponds to the expected ice production in a river and the results used in designing the storage, dam height and other features of the project.
Repeated periods of thaw will disrupt the process and make the estimate of the extent and volume of the ice jam very difficult.

Fig. 1 .
Fig. 1.The Thjórsá river system with glaciers indicated in grey, scale in km.

Fig. 2 .
Fig. 2. Accumulated solid ice volume produced in the river Thjorsa from Burfell to Urridafoss.Water discharge is set at Burfell Power Station operating flow, 300 m 3 /s.

Fig. 3 .
Fig. 3. Calculated ice discharge at Gauging Station 30 at Krokur.River discharge is taken same as in Fig. 2, 200 m 3 s −1 .Horizontal bars indicate days with ice observed.Light blue bars indicate slush or frazil ice runs in the river.Dark blue bars indicate ice covered river.

Fig. 4 .
Fig. 4. Ice jam at Urridafoss in the Thjorsa river, theory compared to observations in Table2.