Methods for Estimating Flow Discharge in Ice-Covered Channels

1. Introduction

Ice-covered channels are widespread in high-latitude and high-altitude regions, and their freezing and thawing processes have a significant impact on local hydrological processes, ecosystems, and human activities1. In winter, the river surface freezes, while water still flows beneath, creating complex hydrodynamic conditions. The appearance of ice jams causes a significant rise in water levels upstream, leading to flooding disasters, destroying local infrastructure, and posing a serious threat to human life and property.

In non-frozen rivers, water levels are often continuously monitored by instruments at hydrological stations, and the monitoring data is then converted into flow discharge using the stage-discharge curve established in open-water flow[2,3]. However, the presence of ice cover alters the river's flow structure, causing changes in the stage-discharge relationship curve, a phenomenon that has garnered global attention[4,5,6,7]. After the river surface freezes, the combined effect of the roughness at the bottom of the ice cover and the bed roughness significantly reduces the flow velocity, leading to an increase in water depth. In practical applications, the flow capacity during the ice period is reduced by 30%-37% compared to the open-water conditions, significantly lowering the transport capacity[8,9]. Meanwhile, the position of the maximum flow velocity shifts downward compared to open-water flow, and as the ice cover roughness increases, the location of the maximum velocity gradually moves closer to the bed10. The velocity distribution and along-flow resistance, which directly affects the transport efficiency of water under the ice, making this research highly important11. To determine river flow affected by ice conditions, the velocity-area method is often employed12. This manual measurement method demands significant labor and may expose operators to safety risks. It is extremely necessary to develop a simple and effective method for estimating flow during freezing periods.

In hydraulic methods for calculating winter flow, the Manning equation has the inherent advantage of not needing flow velocity measurements, which are challenging to obtain in frozen conditions. Hydraulic calculations for ice-covered channels are concentrated on deriving the comprehensive resistance coefficient, which incorporates boundary friction effects. Since 1931, Pavlovskiy13 proposed the first calculation formula for ice-covered channels, and famous scholars such as Lotter14, Einstein15, Larsen16, and Sabaneev17 have continued to deepen the research and propose comprehensive roughness calculation formulas. Uzuner18 analyzed the research of these scholars and concluded that the Larsen formula is the best method, although it has fewer equations than unknowns and cannot provide a method for calculating partition depths. With the use of current methods, we can establish the stage-discharge relationship for rivers impacted by ice conditions in cold winter regions. Since it does not require the collection of region-specific ice data, this method allows for a more convenient improvement of flow estimation in ice-covered channels. Therefore, the aim of this research is to establish a physical formula that combines both comprehensiveness and simplicity for estimating flow in ice-covered channels.

The main structure of the research is as follows: (1) Based on the boundary conditions, we establish the relationship between comprehensive roughness, ice cover roughness, bed roughness, and the hydraulic radius of each section, and propose a general calculation method for comprehensive roughness applicable to ice-covered channels; (2) propose methods for calculating physical quantities; (3) compare the accuracy of the methods in this study with the Sabaneev, Larsen, and Pavlovskiy formulas in estimating flow in ice-covered channels and discuss the causes of errors in each formula; (4) develop two simplified formulas for calculating comprehensive roughness. Finally, the research is dedicated to improve the flow estimation method for ice-covered channels and provide a reference for winter river flow discharge measurement.

2. Methodology

2.1. Analytical Method for Ice-Covered Channel Resistance Calculation

When calculating the resistance of partially or completely frozen flows, the Manning equation or the Darcy-Weisbach equation is commonly used. In the case of the Manning equation, the transport capacity of frozen rivers can be expressed as[19,20]

$$Q=AV={\displaystyle \frac{1}{n}}\chi R^{5/3}\sqrt{S}$$

(1)

Where $V$ is the mean velocity, $R$ is the hydraulic radius, $n$ is the comprehensive roughness of the total flow section, $S$ is the energy slope. Empirical formulas based on hydraulic radius are applied to calculate various coverage conditions from open-water to fully frozen channels. For rectangular cross-section channels, the formula can be expressed as21

$$R={\displaystyle \frac{BH}{B\left(1+a/100\mathrm{\%}\right)+2H}}$$

(2)

Where $B$ is the channel width, $H$ is the flow depth, $a$ is the percentage of the cover.

Accurate calculation of the Manning coefficient for the entire flow section is crucial for using the Manning equation to estimate flow under ice cover conditions. In the study of ice-covered channels, the two-layer assumption proposed by Einstein15 is widely used, adopting a resistance partitioning approach. In this assumption, the location of the maximum flow velocity is used as the boundary, dividing the ice-covered flow area into two equivalent open-water flow layers: the ice cover section and the bed section. The two layers are assumed not to interfere with each other and are only influenced by the roughness of the ice roughness and the bed roughness. The fractal theory, introduced by French mathematician B.B. Mandelbrot in the 1970s, is used to describe self-similar and affine shapes in nature. This theory has found extensive applications in river dynamics22.

In Error! Reference source not found., $H$ is the total flow depth of the ice-covered channel, $y_b$ is the flow depth of the bed section, $y_i$ is the flow depth of the ice cover section, $u$ is the flow velocity, ${\chi }_i$ is the wetted perimeter of the ice cover section, ${\chi }_b$ is the wetted perimeter in the bed section, $A_i$ is the cross-sectional area of the ice cover section,$A_b$ is the cross-sectional area of the bed section.

Figure 1. (a) Ice-covered channel's ice cover and bed sections (b) Velocity distribution in the ice-covered channel along the flow depth direction.

Figure 1. (a) Ice-covered channel's ice cover and bed sections (b) Velocity distribution in the ice-covered channel along the flow depth direction.

2.2. General Resistance Formula for Ice-Covered Channels

According to the law of mass conservation, the total flow discharge in the channel equals the sum of the flow discharge in the ice cover section and the bed section

$$Q=Q_b+Q_i$$

(3)

Where $Q$ is the total flow discharge, $Q_i$ is the flow discharge in the ice cover section, ${\mathrm{Q}}_b$ is the flow discharge in the bed section.

The hydraulic radius in the ice cover section and the bed section satisfies the relationship

$$A_i/{\chi }_i=R_i$$

(4)

$$A_b/{\chi }_b=R_b$$

(5)

Based on boundary conditions, it is known

$$\chi ={\chi }_i+{\chi }_b$$

(6)

$$A=A_i+A_b$$

(7)

By using Eq.(1) to express the flow discharge, we can derive

$$Q_i={\displaystyle \frac{1}{n_i}}{\chi }_i{R}_i^{5/3}\sqrt{S_i}$$

(8)

$$Q_b={\displaystyle \frac{1}{n_b}}{\chi }_b{R}_b^{5/3}\sqrt{S_b}$$

(9)

The total hydraulic radius can be expressed as

$$R={\displaystyle \frac{A}{\chi }}={\displaystyle \frac{A_b+A_i}{{\chi }_b+{\chi }_i}}={\displaystyle \frac{R_b+P{R}_i}{1+P}}$$

(10)

Where $P={\chi }_i/{\chi }_b$, is the ratio of the wetted perimeter in the bed section to that in the ice cover section, a dimensionless parameter.

Assuming $S=S_i=S_b$, this assumption holds perfectly only when the flow is uniform. However, the flow in a frozen river may not be stable and uniform in fact. Despite this, many scholars still favor this assumption[23,24], and it has been extensively used in calculating the resistance of ice-covered channels.

From Eq.(3), Eq.(8) and Eq.(9), we can derive

$${\displaystyle \frac{1}{n}}=\left({\displaystyle \frac{1}{n_b}}{\displaystyle \frac{{\chi }_b}{\chi }}+{\displaystyle \frac{1}{n_i}}{\displaystyle \frac{{\chi }_i}{\chi }}{\displaystyle \frac{R_i^{5/3}}{R_b^{5/3}}}\right){\displaystyle \frac{R_b^{5/3}}{R^{5/3}}}$$

(11)

By substituting Eq.(10) into Eq.(11), we get

$$n={\displaystyle \frac{{\left(1+{PR}_r\right)}^{5/3}}{{\left(1+P\right)}^{2/3}\left[1+{{P\left(R_r\right)}^{5/3}/R}_n\right]}}{n}_b$$

(12)

Where $R_r=R_i/R_b$, is the ratio of the hydraulic radius between the bed section and the ice cover section, a dimensionless parameter, $R_n=n_i/n_b$, is the ratio of the Manning roughness between the ice cover section and the bed section, a dimensionless parameter.

If we assume $R_i=R_b=R$ i.e., the hydraulic radii of all sections are equal, by substituting into Eq.(12), we can obtain Lotter’s formulation14

$$n={\displaystyle \frac{1+P}{1+{P/R}_n}}{n}_b$$

(13)

If we assume $V_i=V_b=V$, meaning the flow velocities of all sections are equal. Eq.(1) can be used to obtain

$$R_n={\left(R_r\right)}^{2/3}$$

(14)

By substituting Eq.(14) into Eq.(12), we can derive Sabaneev's formulation17

$$n={\left[{\displaystyle \frac{1+P{\left(R_n\right)}^{3/2}}{1+P}}\right]}^{2/3}{n}_b$$

(15)

If both assumptions of equal flow velocity and equal hydraulic radius are used, despite the serious issue that this leads to equal Manning’s roughness coefficients in the ice cover and bed sections, we can obtain Pavlovskiy’s formulation13

$$n={\left[{\displaystyle \frac{1+P{\left(R_n\right)}^2}{1+P}}\right]}^{1/2}{n}_b$$

(16)

Based on the assumption of equal velocities, by adding the assumption that ${\chi }_i={\chi }_b=\chi /2$, we can obtain Sun Z’s formulation25

$$n={(1/2)}^{2/3}{\left[1+{\left(R_n\right)}^{3/2}\right]}^{2/3}{n}_b$$

(17)

If we assume the channel is wide, meaning ${\chi }_i={\chi }_b$. and that the hydraulic radius of each section equals the flow depth of the respective section, substituting these assumptions into Eq.(12) yields Larsen’s formulation16

$$n={\displaystyle \frac{{\left(1/2\right)}^{2/3}{\left(y_i/y_b+1\right)}^{5/3}}{{\left(y_i/y_b\right)}^{5/3}/R_n+1}}{n}_b$$

(18)

As discussed above, many resistance prediction formulas for ice-covered channels are simplifications of the general formula proposed in this study under specific assumptions. Eq.(12) provides a deeper explanation of the factors affecting resistance in ice-covered channels, which is significant for studying flow discharge in ice-covered channels.

2.3. Division of Regions in Ice-Covered Channels

The velocity of the channel is an important basis for studying sediment transport, river evolution, flow discharge, and other related issues26, By integrating the average flow velocity of the river, the transport capacity under specific flow depths can be computed. To derive the turbulent velocity distribution, Prandtl L. assumed (1) that the shear stress near the flow equals the wall shear stress, and (2) that the mixing length is proportional to the distance of the flow particle from the wall. In high Reynolds number flow, he derived the logarithmic velocity distribution formula for turbulent flow. However, this assumption has theoretical limitations, such as the need for fluid particles to pass through the mixing length before exchanging momentum with nearby particles, which contradicts the fact that momentum exchange occurs continuously as particles move. Additionally, when the logarithmic velocity distribution is applied to ice-covered channels, its inherent properties lead to discontinuities in the location of the maximum velocity point. Research[27,28,29,30] indicates that the logarithmic velocity distribution results in significant errors in the core flow area when describing the velocity distribution under the ice cover. When applying the logarithmic velocity distribution to describe flow beneath the ice cover, around 70% of the velocity distribution in the ice cover and bed sections aligns with the logarithmic distribution, but in about 30% of the area near the maximum velocity point, there is a significant discrepancy where the actual velocity is lower than the predicted value.

To address the issue of discontinuous velocity distribution in the theoretical model, Tsai and Ettema29, based on Odggard's research31, proposed using the two-power law to predict the velocity distribution under the ice cover, with two exponents representing the resistance in the ice cover and bed sections. The larger the value of the exponent, the smoother the boundary, and the smaller the flow resistance. Research by Chee and Haggag27, Teal29, and Li32 concluded that the two-power law can describe the continuous distribution of velocity with a single expression, avoiding the inherent flaws of the logarithmic velocity distribution law. Experiments and field observations indicate that the two-power law matches measurement data better[32,33] and facilitates the determination of the maximum velocity location, which has led to its application in research. The two-power law velocity distribution equation is

$$u=K_0{\left(y/H\right)}^{1/m_b}{\left(1-y/H\right)}^{1/m_i}$$

(19)

Where $K_0$ is the constant for the given flow velocity, $m_b$ and $m_i$ are the dimensionless parameters related to the bed and ice cover, and $y$ is the vertical depth. Teal et al. used over 22,300 vertical velocity profiles collected from 13 river flow measurement stations by the U.S. Geological Survey during the winter of 1988-1989, applying nonlinear regression methods. They estimated the exponent values to range from 1.5 to 8.5, with most values between 3 and 4. For open-channel flow, the exponent values range from 6 to 7.

The velocity gradient under the two-power law velocity distribution can be obtained from Eq.(19) as

$${\displaystyle \frac{du}{dH}}={\displaystyle \frac{K_0}{h}}\left[{\displaystyle \frac{1}{m_b}}{\left({\displaystyle \frac{y}{H}}\right)}^{{\displaystyle \frac{1}{m_b}}-1}{\left(1-{\displaystyle \frac{y}{H}}\right)}^{{\displaystyle \frac{1}{m_i}}}-{\displaystyle \frac{1}{m_i}}{\left({\displaystyle \frac{y}{H}}\right)}^{{\displaystyle \frac{1}{m_b}}}{\left(1-{\displaystyle \frac{y}{H}}\right)}^{{\displaystyle \frac{1}{m_i}}-1}\right]$$

(20)

Let $du/dh=0$, we can solve for the maximum velocity under the ice cover as

$$u_{max}=K_0{\left(h_m/H\right)}^{1/m_b}{\left(1-h_m/H\right)}^{1/m_i}$$

(21)

Where the maximum velocity occurs

$${\displaystyle \frac{h_m}{H}}={\displaystyle \frac{m_i}{m_b+m_i}}={\displaystyle \frac{R_m}{1+R_m}}$$

(22)

Where $h_m$ is the depth where maximum velocity occurs, and $R_m=m_i/m_b$, is the ratio of the exponents between the ice cover and bed sections.

2.4. New Predictor for Frozen River Estimation

By substituting Eq.(12) ino Eq.(1), we obtain

$$V={\displaystyle \frac{{\left(1+P\right)}^{2/3}\left[1+{{P\left(R_r\right)}^{5/3}/R}_n\right]}{{\left(1+{PR}_r\right)}^{5/3}}}{\displaystyle \frac{1}{n_b}}{R}^{2/3}\sqrt{S}$$

(23)

The research focuses on fully covered flow($a=1$), thus

$$R={\displaystyle \frac{BH}{2B+2H}}={\displaystyle \frac{\beta }{2\left(1+\beta \right)}}H$$

(24)

Where $\beta =B/H$ is the width-to-depth ratio of the channel. If $\beta \to \mathrm{\infty }$, then $R\to H/2$, which implies that the hydraulic radius of the flow section is half of the flow depth, a common assumption in wide channels.

Eq.(23) can be written as

$$V=K{\displaystyle \frac{1}{n_b}}{R}^{2/3}\sqrt{S}$$

(25)

$$K={\displaystyle \frac{{\left(1+P\right)}^{2/3}\left[1+{{P\left(R_r\right)}^{5/3}/R}_n\right]}{{\left(1+{PR}_r\right)}^{5/3}}}$$

(26)

Where $K$is a physically based coefficient that indicates the correlation between the flow discharge of an ice−affected channel and that of an open−water channel. When the flow is converted to open−water flow, we easily find that $K=1$. In this case, Eq.(25) simplifies to the Manning equation, indicating that the Manning equation is a special case with zero cover roughness34. By calculating the value of $K$, the flow discharge of the ice-affected channel can be estimated using the stage-discharge curve established for open-water flow.

2.5. Estimation of Parameters

The Manning formula for steady, uniform flow in open channels is related to the Darcy-Weisbach coefficient through the following equation

$${\displaystyle \frac{R^{1/6}}{n\sqrt{g}}}=\sqrt{8/f}$$

(27)

Where $f$ is the Darcy-Weisbach resistance coefficient for the total flow cross-section,，and $g$ is the local gravitational acceleration. It is assumed that each section also satisfies a similar resistance relationship35.

$$n_b=R_b^{1/6}\sqrt{f_b/\left(8g\right)}$$

(28)

$$n_i=R_i^{1/6}\sqrt{f_i/\left(8g\right)}$$

(29)

Where $f_b$ and $f_i$ are the Darcy-Weisbach resistance coefficients for the bed and ice cover sections.

The relationship between the exponents and the Darcy-Weisbach resistance coefficients is described by36

$$m_b=\mathsf{\kappa }{\left(8/f_b\right)}^{1/2}$$

(30)

$$m_i=\mathsf{\kappa }{\left(8/f_i\right)}^{1/2}$$

(31)

where κ is the von Karman constant.

Therefore, we can obtain

$${n_i=R}_i^{1/6}\cdot {\displaystyle \frac{\kappa }{m_i\sqrt{g}}}$$

(32)

$$n_b=R_b^{1/6}\cdot {\displaystyle \frac{\kappa }{m_b\sqrt{g}}}$$

(33)

$$R_n={\left(R_r\right)}^{1/6}/R_m$$

(34)

From Eq.(8) and Eq.(9), we can obtain

$$\alpha ={\displaystyle \frac{{1/n}_i}{{1/n}_b}}{\displaystyle \frac{R_i^{2/3}}{R_b^{2/3}}}=R_m{\left(R_r\right)}^{1/2}$$

(35)

Where $\alpha$ is the velocity ratio between the ice cover section and the bed section.

From section2.3, we can get

$$\alpha ={\displaystyle \frac{1/\left(h-h_m\right){\int }_{h_m}^{h}u\hspace{0.17em}dy}{1/h_m{\int }_0^{h_m}u\hspace{0.17em}dy}}={\displaystyle \frac{h_m}{\left(h-h_m\right)}}{\displaystyle \frac{{\int }_{\frac{h_m}{h}}^1{t}^{\frac{1}{m_b}}{\left(1-t\right)}^{\frac{1}{m_i}}\hspace{0.17em}dt}{{\int }_0^{\frac{h_m}{h}}{t}^{\frac{1}{m_b}}{\left(1-t\right)}^{\frac{1}{m_i}}\hspace{0.17em}dt}}$$

(36)

From Eq.(35) and Eq.(36), we can obtain

$$R_r={\left({\int }_{{\displaystyle \frac{R_m}{1+R_m}}}^1{t}^{{\displaystyle \frac{1}{m_b}}}{\left(1-t\right)}^{{\displaystyle \frac{1}{m_i}}}\hspace{0.17em}dt/{\int }_0^{{\displaystyle \frac{R_m}{1+R_m}}}{t}^{{\displaystyle \frac{1}{m_b}}}{\left(1-t\right)}^{{\displaystyle \frac{1}{m_i}}}\hspace{0.17em}dt\right)}^2$$

(37)

If we continue the investigation, many interesting results can be found. Based on the previous discussion, the calculation process becomes relatively simple, and only the results are presented.

$$R_b=AP/\left[{(R}_r+P){\chi }_b\right]$$

(38)

$$R_i=A{R}_r/\left[{(R}_r+P){\chi }_i\right]$$

(39)

2.6. Estimation of $m_i$ and $m_b$

In practical engineering, the flow cross-section of natural rivers is often simplified to a two-dimensional rectangular or trapezoidal shape, with geometric properties like channel width, flow depth, and slope easily measured directly. Based on the above basic information, the research provides three methods to estimate the exponents $m_i$ and $m_b$.

• Given the velocity distribution along the cross-section, the regression method proposed by Attar and Li32 is applied to obtain the values;
• Given the Darcy-Weisbach resistance coefficients f_b and fi, they can be estimated directly through Eq.(30) and Eq.(31);
• Given the Manning roughness coefficients n_b and ni, by assuming an initial value of m_b, iterative solutions are obtained using Eq.(10), Eq.(32), Eq.(33) and Eq.(37).

3. Results and Discussion

3.1. Verification of the Proposed Formula

3.1.1. Data Collection

The applicability of the proposed formula is verified using experimental data from fully covered channels, including 29 sets of laboratory measurements by Parthasarathy and Muste [34], Wei and Huang [35], Engmann [37], Smith and Ettema [38] and Zhang J [39] and 12 sets of field observations in natural river by Attar and Li [32] and Tatinclaux and Göğüs [40]. The geometric parameters and hydraulic conditions of the channel are listed in Error! Reference source not found..

Table 1. Geometric parameters and hydraulic conditions used for validation.

Data Sources	Runs	$\mathit{B}\left(\mathit{m}\right)$	$\mathbf{S}\left(\mathbf{‰}\right)$	$\mathit{H}\left(\mathit{m}\right)$	${\mathit{m}}_{\mathit{b}}$	${\mathit{m}}_{\mathit{i}}$	$\mathit{Q}\left({\mathit{m}}^3/\mathit{s}\right)$
Smith and Ettema (1997)	SE-S2	0.912	1.37	0.181	4.51	7.56	0.0754
	SE-M2	0.912	1.29	0.195	4.79	5.72	0.0760
	SE-R2	0.912	1.33	0.208	4.7	4.56	0.0756
	SE-S4	0.912	1.34	0.182	4.59	7.69	0.0765
	SE-M4	0.912	1.3	0.19	4.9	5.85	0.0755
	SE-R4	0.912	1.3	0.209	4.7	4.56	0.0754
Tatinclaux and Gogus (1983)	TG-C1	425	0.5	5	5.73	2.44	1850
	TG-C2	425	0.5	4	5.47	2.44	1230
	TG-C3	425	0.5	3.5	5.04	1.85	850
Parthasarathy and Muste (1994)	PM-R1	0.912	0.197	0.218	7.02	8.47	0.0501
	PM-R2	0.912	0.197	0.245	6.63	6.35	0.0501
	PM-R3	0.912	0.197	0.29	5.73	4.75	0.0501
Wei and Huang (2002)	WH-Test 1	0.5	0.478	0.24	9.68	8.26	0.0501
	WH-Test 2	0.5	0.459	0.241	9.68	8.26	0.0501
	WH-Test 3	0.5	0.446	0.242	9.68	8.26	0.0501
	WH-Test 4	0.5	1.072	0.218	9.58	8.16	0.0699
	WH-Test 5	0.5	1.124	0.199	9.48	8.08	0.0600
	WH-Test 6	0.5	0.896	0.165	9.26	7.92	0.0401
	WH-Test 7	0.5	0.774	0.144	9.26	7.92	0.0303
	WH-Test 8	0.5	0.66	0.22	9.58	8.15	0.0507
	WH-Test 10	0.5	2.845	0.192	8	3.15	0.0500
	WH-Test 11	0.5	3.058	0.211	8	3.15	0.0602
	WH-Test 12	0.5	2.526	0.201	8	3.15	0.0502
	WH-Test 13	0.5	2.534	0.201	8	3.15	0.0502
	WH-Test 15	0.5	2.572	0.236	3.45	3.04	0.0507
	WH-Test 16	0.5	2.277	0.217	3.45	3.04	0.0412
J Zhang(2021)	Case1	1	1	0.15	6.35	4.84	0.0536
	Case2	1	1	0.185	7.13	5.31	0.0783
Engmann (1977)	EN-101	1.22	0.65	0.0497	3.1	6.46	0.0071
	EN-102	1.22	0.79	0.0649	4.12	8.08	0.0156
	EN-103	1.22	2.49	0.0384	4.6	7.54	0.0127
	EN-104	1.22	1.61	0.0396	4.6	7.54	0.0114
Attar and Li (2012)	S.W. Miramichi R., NB	92	0.07	2	3.59	7.39	51
	Burnt R., ON	32	0.04	1.9	3.2	5.48	10
	Pembina R., AB	74	0.13	0.7	3.23	6.25	12
	Halfway R., BC	39	0.8	0.54	2.77	5.96	7.4
	Peace R., NWT	525	0.04	4.5	5.44	9.22	1111
	Yellowknife R., NWT	72	0.01	3	3.55	5.92	24
	Fraser R., BC	95	0.1	1.3	3.25	6.37	32
	Takhini R. YT	46	0.08	1.4	3.12	5.96	14
	Yukon R., YT	145	0.4	2.5	3.69	7.06	246

Table 1. Geometric parameters and hydraulic conditions used for validation.

Data Sources	Runs	$\mathit{B}\left(\mathit{m}\right)$	$\mathbf{S}\left(\mathbf{‰}\right)$	$\mathit{H}\left(\mathit{m}\right)$	${\mathit{m}}_{\mathit{b}}$	${\mathit{m}}_{\mathit{i}}$	$\mathit{Q}\left({\mathit{m}}^3/\mathit{s}\right)$
Smith and Ettema (1997)	SE-S2	0.912	1.37	0.181	4.51	7.56	0.0754
	SE-M2	0.912	1.29	0.195	4.79	5.72	0.0760
	SE-R2	0.912	1.33	0.208	4.7	4.56	0.0756
	SE-S4	0.912	1.34	0.182	4.59	7.69	0.0765
	SE-M4	0.912	1.3	0.19	4.9	5.85	0.0755
	SE-R4	0.912	1.3	0.209	4.7	4.56	0.0754
Tatinclaux and Gogus (1983)	TG-C1	425	0.5	5	5.73	2.44	1850
	TG-C2	425	0.5	4	5.47	2.44	1230
	TG-C3	425	0.5	3.5	5.04	1.85	850
Parthasarathy and Muste (1994)	PM-R1	0.912	0.197	0.218	7.02	8.47	0.0501
	PM-R2	0.912	0.197	0.245	6.63	6.35	0.0501
	PM-R3	0.912	0.197	0.29	5.73	4.75	0.0501
Wei and Huang (2002)	WH-Test 1	0.5	0.478	0.24	9.68	8.26	0.0501
	WH-Test 2	0.5	0.459	0.241	9.68	8.26	0.0501
	WH-Test 3	0.5	0.446	0.242	9.68	8.26	0.0501
	WH-Test 4	0.5	1.072	0.218	9.58	8.16	0.0699
	WH-Test 5	0.5	1.124	0.199	9.48	8.08	0.0600
	WH-Test 6	0.5	0.896	0.165	9.26	7.92	0.0401
	WH-Test 7	0.5	0.774	0.144	9.26	7.92	0.0303
	WH-Test 8	0.5	0.66	0.22	9.58	8.15	0.0507
	WH-Test 10	0.5	2.845	0.192	8	3.15	0.0500
	WH-Test 11	0.5	3.058	0.211	8	3.15	0.0602
	WH-Test 12	0.5	2.526	0.201	8	3.15	0.0502
	WH-Test 13	0.5	2.534	0.201	8	3.15	0.0502
	WH-Test 15	0.5	2.572	0.236	3.45	3.04	0.0507
	WH-Test 16	0.5	2.277	0.217	3.45	3.04	0.0412
J Zhang(2021)	Case1	1	1	0.15	6.35	4.84	0.0536
	Case2	1	1	0.185	7.13	5.31	0.0783
Engmann (1977)	EN-101	1.22	0.65	0.0497	3.1	6.46	0.0071
	EN-102	1.22	0.79	0.0649	4.12	8.08	0.0156
	EN-103	1.22	2.49	0.0384	4.6	7.54	0.0127
	EN-104	1.22	1.61	0.0396	4.6	7.54	0.0114
Attar and Li (2012)	S.W. Miramichi R., NB	92	0.07	2	3.59	7.39	51
	Burnt R., ON	32	0.04	1.9	3.2	5.48	10
	Pembina R., AB	74	0.13	0.7	3.23	6.25	12
	Halfway R., BC	39	0.8	0.54	2.77	5.96	7.4
	Peace R., NWT	525	0.04	4.5	5.44	9.22	1111
	Yellowknife R., NWT	72	0.01	3	3.55	5.92	24
	Fraser R., BC	95	0.1	1.3	3.25	6.37	32
	Takhini R. YT	46	0.08	1.4	3.12	5.96	14
	Yukon R., YT	145	0.4	2.5	3.69	7.06	246

3.1.2. Model Verifications

The proposed Eq.(23) is compared with the flow discharge results calculated using the traditional Lotter, Sabaneev, Larsen, and Pavlovskiy formulas. Since there are significant differences in flow discharge under different conditions, velocity errors are used instead of flow discharge for comparison. The percentage error in velocity is equal to the percentage error in flow discharge, making this substitution method valid. We use the relative error $M$ to evaluate the performance of each formula.

$$M={\displaystyle \frac{1}{N}}\sum _{i=1}^N\left|{\displaystyle \frac{V_{cal}-V_0}{V_0}}\right|$$

Where $N$ is the number of all verified conditions$(N=41$), $V_{cal}$ is the calculated velocity, $V_0$ is the measured velocity.

The comparison of the measured velocity results with the actual values shows that the average error $M$ of the proposed formula Eq.(23) is 3.97%, with a maximum error of 15.08% and a minimum error of 0.31%, shown in Error! Reference source not found.; for the Lotter formula, the average error is 16.35%, the maximum error is 55.55%, and the minimum error is 0.31%, shown in Error! Reference source not found.; for the Sabaneev formula, the average error is 5.23%, the maximum error is 15.94%, and the minimum error is 0.40%, shown in Error! Reference source not found.; for the Larsen formula, the average error is 5.89%, the maximum error is 20.93%, and the minimum error is 0.18%, shown in Error! Reference source not found.; for the Pavlovskiy formula, the average error is 6.82%, the maximum error is 17.97%, and the minimum error is 0.11%, shown in Error! Reference source not found.. Compared to traditional formulas, the proposed formula Eq.(23) not only accurately predicts the flow discharge in ice-covered channels, but also shows a significant improvement in performance.

Figure 2. Comparison of estimated results of Eq.(23) with the measured data.

Figure 2. Comparison of estimated results of Eq.(23) with the measured data.

Figure 3. Comparison of estimated results of the Lotter formula with the measured data.

Figure 3. Comparison of estimated results of the Lotter formula with the measured data.

Figure 4. Comparison of estimated results of the Sabaneev formula with the measured data.

Figure 4. Comparison of estimated results of the Sabaneev formula with the measured data.

Figure 5. Comparison of estimated results of the Larsen formula with the measured data.

Figure 5. Comparison of estimated results of the Larsen formula with the measured data.

Figure 6. Comparison of estimated results of the Pavlovskiy formula with the measured data.

Figure 6. Comparison of estimated results of the Pavlovskiy formula with the measured data.

3.2. Hydraulic Radius

It is widely assumed in past calculations of resistance in ice-covered channels that the hydraulic radius of the sections are equal, i.e., $R=R_i=R_b$. This assumption does not reflect reality, as it leads to the conclusion that the ice cover roughness and bed roughness are nearly equal, which is obviously contradictory. As shown in Error! Reference source not found., in the given channel , there is a negative correlation between the hydraulic radius of the ice cover section and that of the bed section, and the total hydraulic radius is between the two. Studies indicate that in the vast majority of cases, the difference between the hydraulic radius of the ice cover section and the bed section and the total hydraulic radius exceeds 10%, and sometimes even exceeds 30%, making the assumption unreasonable.

Error! Reference source not found. shows the trend of $R_r$ as $R_m$ varies. The result shows that the Pearson correlation coefficient between the hydraulic radius ratio and the exponent ratio is $r=-0.8032$, indicating a strong negative correlation. Based on the data, we use a function to fit $R_r$ and $R_m$ and evaluate the fitting curve with the Coefficient of Determination (COD) and Standard Error (SE) indicators.

$$R_r={R_m}^{-1.694}$$

(40)

Where $\mathrm{C}\mathrm{O}\mathrm{D}=0.9960,\mathrm{S}\mathrm{E}=0.0084$. Eq.(40) provides a simple method for calculating the Manning roughness or the exponent.

Figure 7. (a) Comparison of $R$, $R_i$ and $R_b$ (b) Magnified view of local details.

Figure 7. (a) Comparison of $R$, $R_i$ and $R_b$ (b) Magnified view of local details.

Figure 8. Variations of $R_r$ and $R_m.$

Figure 8. Variations of $R_r$ and $R_m.$

3.3. Velocity

The assumption of equal average velocity across sections is widely favored by researchers studying ice-covered channels[41,42,43]. According to the results of the study, the assumption of equal velocity across sections is not suitable, as the average velocity in each section is greatly influenced by the roughness of that section.

Error! Reference source not found. shows the relationship between $\alpha$ and $R_m$. The results show a positive correlation between $\alpha$ and $R_m$. In the validation data, $\alpha$ values typically range from 0.9 to 1.1, suggesting that the assumption of equal velocity causes a 5% to 10% error in most cases. $\alpha =1$ only when the ice cover roughness is equal to the bed roughness. The research advises against using the assumption of equal velocity to calculate flow discharge during the initial freezing stage of ice-affected channels in winter, as the large difference in roughness between the ice cover and bed would lead to large errors.

Figure 9. Variations of $\alpha$ and $R_m.$

Figure 9. Variations of $\alpha$ and $R_m.$

3.4. Manning Resistance Coefficient

Teal29 pointed out that the larger the exponent, the smoother the boundary and the lower the flow resistance. In previous studies44, Eq.(14) was used to estimate flow resistance, and as discussed in Section 2.2, this simplified formula is based on the assumption of equal velocity. Section 3.3 demonstrates that the equal velocity assumption carries some degree of error, especially when the roughness difference between the ice cover and bed sections is significant. Therefore, the research delves further into the relationship between Manning's roughness, the exponents, and hydraulic radius. As shown in Error! Reference source not found., $R_n$ and $R_m$ display a distinct negative correlation, with a Pearson correlation coefficient of $r=-0.8459$. By fitting the curve, we derived

$$R_n={R_m}^{-1.282}$$

(41)

Where $\mathrm{C}\mathrm{O}\mathrm{D}=0.9999,\mathrm{S}\mathrm{E}=0.0052$.

By combining Eq.(40) and Eq.(41), we deduced

$$R_r={R_n}^{1.322}$$

(42)

Figure 10. Variations of $R_n$ and $R_m.$

Figure 10. Variations of $R_n$ and $R_m.$

3.5. Two Simplified Formulas for Flow Prediction

By substituting the results of Section 3.2 and Section 3.4 into Eq.(23), the research derived two simplified formulas for flow prediction.

$$V={\displaystyle \frac{{\left(1+P\right)}^{2/3}\left[1+P{R_n}^{C_2}\right]}{{\left(1+P{\left(R_n\right)}^{C_1}\right)}^{5/3}}}{\displaystyle \frac{1}{n_b}}{R}^{2/3}\sqrt{S}$$

(43)

Where $C_1$ and $C_2$ are empirical parameters, with the recommended values are$C_1=1.32$，$C_2=1.20$.

$$V={\displaystyle \frac{{\left(1+P\right)}^{2/3}\left[1+P{\left(R_m\right)}^{C_4}\right]}{{\left(1+P{\left(R_m\right)}^{C_3}\right)}^{5/3}}}{\displaystyle \frac{1}{n_b}}{R}^{2/3}\sqrt{S}$$

(44)

Where $C_3$ and $C_3$ are empirical parameters, with the recommended values are$C_3=-1.69$ and $C_3=-1.54$.

Compared with measured data, the average error $M$ of the proposed formula Eq.(43) is 4.20%, with a maximum error of 15.17% and a minimum error of 0.23%, shown in Error! Reference source not found.; the average error $M$ of the proposed formula Eq.(44) is 4.16%, with a maximum error of 14.93% and a minimum error of 0.49%, shown in Error! Reference source not found.. The two simplified formulas, while maintaining high computational accuracy, greatly reduce the labor costs and time required for measurements and calculations in practical engineering, providing operators with more options to calculate flow discharge in ice-covered channels. For example, using the channel resistance before freezing to replace the bed resistance after freezing, or using methods to measure roughness height for estimating roughness, or measuring the exponents of velocity distribution, etc. Hence, we highly recommend using Eq.(43) or Eq.(44) to calculate flow discharge in ice-covered channels in practical engineering.

Figure 11. Comparison of estimated results of Eq.(43) with the measured data.

Figure 11. Comparison of estimated results of Eq.(43) with the measured data.

Figure 12. Comparison of estimated results of Eq.(44) with the measured data.

Figure 12. Comparison of estimated results of Eq.(44) with the measured data.

3.6. Shortcomings

The assumption of equal energy slopes in each section, as used by Larsen, P. A16 and Ashton, G. D45, essentially assumes that the flow under the ice cover is uniform, which is also adopted in this study. However, the flow in natural rivers is quite complex, and may not always form a stable uniform flow. The energy slope in the bed section is often smaller than that in the ice cover section46, which may negatively affect the accuracy of the research. While field measurements suggest the results are trustworthy, the data sample is limited, requiring more field data for validation.

Eq.(43) and Eq.(44) rely on the determination of empirical parameters, the research provides optional values, but more field data is needed in the future to determine the range of the parameters.

4. Conclusions

The research derives a flow discharge calculation equation for ice-covered channels based on physical conditions and provides two simple formulas, which can easily estimate the flow discharge in ice-covered channels. The analysis results are in strong agreement with laboratory measurements and field observations, showing that the formulas proposed in this study can be effectively used for flow discharge prediction in ice-covered channels. The main conclusions are as follows

• Assuming equal flow velocity or equal hydraulic radius in each section leads to errors in predicting resistance or flow discharge in ice-covered channels, especially the latter, which may result in unacceptable errors.
• Compared to commonly used traditional formulas, such as the Lotter formula, Sabaneev formula, Larsen formula, and Pavlovskiy formula, the general formula and two simplified formulas proposed in the research exhibit superior performance. It is recommended to use the methods or simplified formulas presented in the research for more accurate flow discharge prediction.