**1. Introduction**

Most rivers at high altitude in cold northern regions always freeze in the winter and form ice sheets. The ice sheet in some Canadian rivers is at least 0.6 m thick and lasts for at least 4 months [1]. The wetted perimeter of the cross section and flow resistance in ice-covered flows increases with the presence of the ice sheet, which significantly affects the hydraulic characteristics of river and topographical features and greatly change the flow velocity distribution, flow transport capacity and sediment transport rate [2–11]. Therefore, it is necessary to study the ice-covered flows.

Unlike the open channel flow, flows in the ice-covered channel have asymmetric forms. The presence of the ice sheet makes the maximum streamwise velocity appear at the inner center of the flow, and the location of the maximum velocity is generally considered the division point in the asymmetrically distributed flow [12,13]. The asymmetric distribution mainly depends on the roughness of the ice sheet and the riverbed. The division point of the velocity tends to be away from the rougher surface [14]. Previous studies have shown that the main effects of the ice sheets on alluvial channels can be summarized as: they increase the water level (compared to open channels at the same flow rate), reduce

**Citation:** Zhang, J.; Wang, W.; Li, Z.; Li, Q.; Zhong, Y.; Xia, Z.; Qiu, H. Analytical Models of Velocity, Reynolds Stress and Turbulence Intensity in Ice-Covered Channels. *Water* **2021**, *13*, 1107. https:// doi.org/10.3390/w13081107

Academic Editors: Paweł M. Rowi ´nski and Jueyi Sui

Received: 28 February 2021 Accepted: 14 April 2021 Published: 17 April 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

the average flow velocity, increase the channel drag force and reduce the bed sediment transport rate [15].

Previous investigations on the flow characteristics of ice-covered flows are mainly obtained through experiments [16–18]. Many laboratory experimental data and field observations show that the vertical distribution of streamwise velocity forms a double-layer, which is characterized by the plane of maximum velocity. Shen and Harden (1978) [19] and Lau and Krishnappan (1981) [20] applied this double-layer theory in the ice-covered flows. They divided the ice-covered flows into two separate layers: upper ice layer and lower channel bed. Parthasarathy and Muste (1994) [13] found that the zero shear stress plane in the ice-covered flows was inconsistent with the maximum velocity plane. Chen et al. (2015) [21] proposed that the horizontal plane of zero shear stress should be determined as the dividing plane of sublayers and modified the double-layer assumption.

Some researchers have studied the vertical profiles of the longitudinal velocity by numerical simulation methods and proposed two-dimensional and three-dimensional models [22,23], which may require sensitive hydraulic parameters. However, these hydraulic parameters do not have explicit expressions to be calculated, and they have uncertainty that cannot be ignored. Except the numerical methods, previous researchers also attempted to obtain the vertical distribution of streamwise velocity by theoretical methods [24,25]. Uzuner (1975) [26] separated the ice-covered flow into two layers based on the position of the maximum velocity and assumed that the velocity distributions in the upper ice layer and lower bed layer were consistent with the logarithmic distribution law, and Manning formula could be independently applied to each layer. In addition, a two-power-law function was adopted to calculate the vertical distribution of the longitudinal velocity. Teal et al. (1994) [27] reported a reasonable fit of streamwise velocity data to the two-power-law function. Compared to the two-power-law function, the logarithmic law appears to overestimate velocities near the location of maximum velocity. The two-power-law function is a reasonable extension of the power-law expression that is used to describe velocity profiles of open channel flow. The advantage for the two-power-law function is that it describes the entire flow with a single continuous curve. Hence, the two-power-law function deserves further study.

The objective of this study is to obtain the vertical profiles of longitudinal velocity, Reynolds stress and turbulence intensity in ice-covered channels. Hence, this study focuses on (1) adopting a two-power-law function to calculate the vertical profile of the longitudinal velocity, (2) simplifying the time-averaged momentum equation to analyze the vertical distributions of Reynolds stress, (3) improving the exponential model to calculate the turbulence intensity and (4) employing experimental data to validate the theoretical models and giving detail discussion on the coefficients in the theoretical models.

### **2. Material and Methods**

Since the free water surface is covered with ice, the open channel flow changes to a closed conduit flow but retains the flow characteristics of the open channel flow. The vertical distribution of longitudinal velocity significantly changes in the open channel flow with ice-cover, which reflects that the plane of maximum longitudinal velocity shifts from the free water surface to the inner water (see Figure 1). Here, we used a Cartesian coordinate system with the *x* axis in the main flow direction and *z* axis in the water depth direction. Corresponding to the *x* and *z* axes, *u* and *w* are the longitudinal and vertical velocities, respectively.

**Figure 1.** Schematic profile of the longitudinal velocity in an iced-covered open channel. **Figure 1.** Schematic profile of the longitudinal velocity in an iced-covered open channel.

#### *2.1. Vertical Distribution of Longitudinal Velocity 2.1. Vertical Distribution of Longitudinal Velocity*

The presence of the ice cover causes the increase of the wetted perimeter, which increases the composite flow resistance. Owing to the different roughness of ice cover and channel bed, the velocity profile is vertically asymmetric (see Figure 1). In Figure 1, the vertical flow structure can be divided into two independent layers at the plane of maximum velocity, i.e., the upper ice layer and lower bed layer. The flow in the upper ice layer is mainly affected by the ice cover, and the lower flow is primarily affected by the channel bed. The vertical location of the maximum velocity is determined by the roughness of the ice cover and channel bed, i.e., the maximum velocity will not occur at the free water surface but near the middle of the water depth. The presence of the ice cover causes the increase of the wetted perimeter, which increases the composite flow resistance. Owing to the different roughness of ice cover and channel bed, the velocity profile is vertically asymmetric (see Figure 1). In Figure 1, the vertical flow structure can be divided into two independent layers at the plane of maximum velocity, i.e., the upper ice layer and lower bed layer. The flow in the upper ice layer is mainly affected by the ice cover, and the lower flow is primarily affected by the channel bed. The vertical location of the maximum velocity is determined by the roughness of the ice cover and channel bed, i.e., the maximum velocity will not occur at the free water surface but near the middle of the water depth.

Tsai and Ettema (1994) [28] adopted a power-law function to predict the velocity profile in the open channel flow. For the ice-covered flow, a two-power-law function was used to obtain the vertical profile of the longitudinal velocity [27]. The first advantage of the two-power-law function is that it describes the flow just using a single continuous curve. The second is that the analytical solution of velocity is adjustable in accordance with the roughness change of channel bed and ice cover. The two-power-law is shown as Tsai and Ettema (1994) [28] adopted a power-law function to predict the velocity profile in the open channel flow. For the ice-covered flow, a two-power-law function was used to obtain the vertical profile of the longitudinal velocity [27]. The first advantage of the two-power-law function is that it describes the flow just using a single continuous curve. The second is that the analytical solution of velocity is adjustable in accordance with the roughness change of channel bed and ice cover. The two-power-law is shown as

$$
\mu = K\_0 \left(\frac{z}{H}\right)^{1/m\_b} \left(1 - \frac{z}{H}\right)^{1/m\_i} \tag{1}
$$

where is the vertical axis with = 0 at the channel bed; is the water depth; is the flow parameter based on a given flow discharge per unit flow width and and are the parameters corresponding to the boundary roughness of the channel bed and ice cover. When parameter approaches infinity, the velocity profile becomes equivalent to a single power law expression, which indicates that the ice cover disappears. where *z* is the vertical axis with *z* = 0 at the channel bed; *H* is the water depth; *K*<sup>0</sup> is the flow parameter based on a given flow discharge per unit flow width and *m<sup>b</sup>* and *m<sup>i</sup>* are the parameters corresponding to the boundary roughness of the channel bed and ice cover. When parameter *m<sup>i</sup>* approaches infinity, the velocity profile becomes equivalent to a single power law expression, which indicates that the ice cover disappears.

The velocity gradient obtained from Equation (1) is The velocity gradient obtained from Equation (1) is

$$\frac{du}{dz} = \frac{K\_0}{H} \left[ \frac{1}{m\_b} \left( \frac{z}{H} \right)^{\frac{1}{m\_b} - 1} \left( 1 - \frac{z}{H} \right)^{1/m\_i} - \frac{1}{m\_i} \left( \frac{z}{H} \right)^{\frac{1}{m\_b}} \left( 1 - \frac{z}{H} \right)^{\frac{1}{m\_i} - 1} \right] \tag{2}$$

By setting the velocity gradient డ௨ డ௭ to zero, we deduced the position of the maximum velocity as By setting the velocity gradient *<sup>∂</sup><sup>u</sup> ∂z* to zero, we deduced the position of the maximum velocity as

$$\left(\frac{z}{H}\right)\_{u\_{\max}} = \frac{h\_m}{H} = \frac{m\_i}{m\_i + m\_b} \tag{3}$$

where ௫ is the maximum longitudinal velocity and ℎ is the height of the maximum velocity from the channel bed. where *umax* is the maximum longitudinal velocity and *h<sup>m</sup>* is the height of the maximum velocity from the channel bed.

To predict the vertical distribution of Reynolds stress in ice-covered channels for a steady uniform flow, the time-averaged momentum equation in the longitudinal direction

*2.2. Vertical Distribution of Reynolds Stress* 

can be simplified to

### *2.2. Vertical Distribution of Reynolds Stress*

To predict the vertical distribution of Reynolds stress in ice-covered channels for a steady uniform flow, the time-averaged momentum equation in the longitudinal direction can be simplified to

$$-\overline{u'w'} + \nu \frac{du}{dz} + \frac{u\_{\ast l}^2 z}{H} = u\_{\ast b}^2 \left(1 - \frac{z}{H}\right) \tag{4}$$

where *u* 0 and *w* <sup>0</sup> are turbulent fluctuations of the longitudinal and vertical velocities, −*u* 0*w*0 is Reynolds stress; *ν* is the flow kinematic viscosity; *u*∗*<sup>i</sup>* is the shear velocity at the ice cover and *u*∗*<sup>b</sup>* is the shear velocity at the channel bed.

Since the viscosity shear stress in the flow is much smaller than the Reynolds stress, the viscosity shear stress can be neglected in Equation (4) [29,30]. Then, there is a linear relationship between the shear stress (*τ*) and the vertical axis (*z*). Equation (4) has another form and is shown as

$$
\pi(z) = \pi\_b - (\pi\_b - \pi\_i)\frac{z}{H} \tag{5}
$$

where *τ*(*z*) is the shear stress at distance *z*; *τ<sup>b</sup>* and *τ<sup>i</sup>* are the Reynolds stresses at the channel bed and ice cover; *z* is the vertical distance from the channel bed and *H* is the water depth.

Following the method of Rowinski and Kubrak (2002) [31], who combine the eddy viscosity and flow conditions by a mixing length concept, the shear stress can be given as

$$
\pi(z) = \rho l^2 \frac{\partial u}{\partial z} \left| \frac{\partial u}{\partial z} \right| \tag{6}
$$

where *ρ* is the flow density and *l* is the proposed mixing length. Chen et al. (2015) [21] demonstrate that the mixing length theory can be applied at a certain distance from the fixed boundaries. Two mixing lengths, *l<sup>b</sup>* and *l<sup>i</sup>* , are proposed corresponding to the two fixed boundaries. The first one is the mixing length of channel bed and the second one is the mixing length of the ice cover. Hence, the linear relationships of two mixing lengths are approximately written as

$$\begin{cases} \begin{array}{c} l\_b = \kappa z, \ 0 \le z \le \delta\_b \\\ l\_i = \kappa (H - z), \ H - \delta\_i \le z \le H \end{array} \end{cases} \tag{7}$$

where *δ<sup>b</sup>* and *δ<sup>i</sup>* are the distances from the channel bed and ice cover, where the mixing length theory is valid. *κ* = 0.41 is the von Karman constant.

By substituting Equation (7) into Equation (6), we updated the expression of shear stress as

$$\begin{cases} \begin{array}{c} \tau(z) = \rho \kappa^2 z^2 \frac{\partial u}{\partial z} \Big| \frac{\partial u}{\partial z} \Big|, \; 0 \le z \le \delta\_b \\\ \tau(z) = \rho \kappa^2 (H - z)^2 \frac{\partial u}{\partial z} \Big| \frac{\partial u}{\partial z} \Big|, \; H - \delta\_i \le z \le H \end{array} \tag{8}$$

The gradient of shear stress from Equation (8) is

$$\begin{cases} \begin{array}{c} \frac{\partial \tau}{\partial z} = 2\rho \kappa^2 z^2 \frac{\partial u}{\partial z} \left( \frac{1}{z} \frac{\partial u}{\partial z} + \frac{\partial^2 u}{\partial z^2} \right), \; 0 \le z \le \delta\_b\\ \frac{\partial \tau}{\partial z} = 2\rho \kappa^2 (H - z)^2 \frac{\partial u}{\partial \overline{z}} \left( \frac{1}{(H - z)} \frac{\partial u}{\partial \overline{z}} - \frac{\partial^2 u}{\partial z^2} \right), \; H - \delta\_l \le z \le H \end{array} \tag{9}$$

By substituting the velocity gradient (Equation (2)) into the Reynolds stress gradient (Equation (9)) and letting *∂τ <sup>∂</sup><sup>z</sup>* = 0, we obtained the two maximum Reynolds stresses, *τbmax* and *τimax*, near the channel bed and ice cover, respectively, which are presented as

$$\begin{cases} \quad \tau\_{b\max} = \rho \kappa^2 z\_{b\max} \frac{2}{\delta z\_{b\max}} \left| \frac{\partial u}{\partial z\_{b\max}} \right|, & 0 \le z \le \delta\_b \\\quad \tau\_{i\max} = \rho \kappa^2 (H - z\_{i\max})^2 \frac{\partial u}{\partial z\_{i\max}} \left| \frac{\partial u}{\partial z\_{i\max}} \right|, & H - \delta\_i \le z \le H \end{cases} \tag{10}$$

where *zbmax* and *zimax* denote the positions where the corresponding maximum shear stresses *τbmax* and *τimax* occurred, and they can be calculated as

$$\begin{cases} \ z\_{b\text{max}} = \frac{-A\_b - \sqrt{A\_b^2 - 4B\_b \mathbb{C}\_b}}{2B\_b} H\\\ z\_{i\text{max}} = \frac{-A\_i - \sqrt{A\_i^2 - 4B\_i \mathbb{C}\_i}}{2B\_i} H \end{cases} \tag{11}$$

where *A<sup>b</sup>* , *B<sup>b</sup>* , *C<sup>b</sup>* , *A<sup>i</sup>* , *B<sup>i</sup>* , and *C<sup>i</sup>* are the constants relevant to the roughness of the channel bed and ice cover. Their expressions are

$$\begin{cases} \begin{cases} & A\_b = -m\_b - 2\frac{m\_i}{m\_b} - 2 \\ & B\_b = (m\_b + m\_i) \left( \frac{1}{m\_b} + \frac{1}{m\_i} \right) \\ & \mathcal{C}\_b = \frac{m\_i}{m\_b} \end{cases} \\\ \begin{cases} & A\_i = -m\_i + 2\frac{m\_i}{m\_b} + 2 \\ & B\_i = -(m\_b + m\_i) \left( \frac{1}{m\_b} + \frac{1}{m\_i} \right) \\ & \mathcal{C}\_i = m\_i - \frac{m\_i}{m\_b} \end{cases} \end{cases} \tag{12}$$

The shear stresses at the fixed boundaries, i.e., the channel bed and ice cover, can be considered the maximum shear stresses in the lower bed layer and upper ice layer. Hence, the linear relationship (Equation (5)) between shear stress and vertical position can be rewritten as

$$
\pi(z) = \tau\_{b\text{max}} - (\tau\_{b\text{max}} - \tau\_{i\text{max}})\frac{z}{H} \tag{13}
$$

### *2.3. Vertical Distribution of Turbulence Intensity*

The flow turbulence characteristics can be presented by the turbulence intensity, which can be calculated by the root mean square of the fluctuating longitudinal velocity, i.e., *urms* = q (*u* 0) 2 . In the open channel flow, the term *urms* reaches its maximum value near the channel bed and it has a linear relationship with the vertical distance. Nezu and Rodi (1986) [32] established an exponential model to describe the vertical distribution of the turbulence intensity. However, in the ice-covered channel flow, it first decreases from the channel bed and reaches its minimum value near the middle flow depth, after which it gradually increases towards to the ice cover. For the profile of the turbulence intensity in the ice-covered channel, we adopted the division scheme to investigate the profile of the turbulence intensity. Based on the location of zero shear stress, its profile is separated into two regions. Then, the exponential model can be applied in each region as follows

$$\begin{cases} \frac{\underline{u\_{\rm rms}}}{\underline{u\_{\rm \*}}} = D\_b e^{-E\_b z/h\_{\rm \tau}}, 0 < z < h\_{\rm \tau} \\\frac{\underline{u\_{\rm rms}}}{\underline{u\_{\rm \*}}} = D\_i e^{-E\_i z/h\_{\rm \tau}}, h\_{\rm \tau} < z < H \end{cases} \tag{14}$$

where *D<sup>b</sup>* , *E<sup>b</sup>* , *D<sup>i</sup>* and *E<sup>i</sup>* are constant parameters.

### **3. Experimental Verification**

Experiments were performed in a rectangular glass flume with a total length of 20 m, a width of 1 m and a depth of 0.5 m in the State Key Laboratory of Water Resources and Hydropower Engineering Science in Wuhan University to investigate the longitudinal velocity profile in the ice-covered flow (Figure 2). The flow is circulated by a pump system and adjusted to be steady by a tailgate at the end of the flume. A plastic foam board was used to simulate the ice with the length of 15 m and a width of 1 m [33].

and adjusted to be steady by a tailgate at the end of the flume. A plastic foam board was

used to simulate the ice with the length of 15 m and a width of 1 m [33].

(**b**)

**Figure 2.** (**a**) Schematic of the experimental flume, the dark blue denoting the flow in the flume and (**b**) images of the experimental site. **Figure 2.** (**a**) Schematic of the experimental flume, the dark blue denoting the flow in the flume and (**b**) images of the experimental site.

The rectangular coordinates are as follows: x is the main flow direction with x = 0 at the beginning of the ice cover; y is the lateral direction with y = 0 at the sidewall and z is the vertical direction counting from the channel bed. The velocity was measured at the cross section with x = 9 m, which is sufficiently far from the ice cover entry to form a fully developed flow. In the measured cross section, 14 measured lines are arranged, where the measured points are also arranged with a uniform vertical distance 1 cm. The detail layout is shown in Figure 3. The rectangular coordinates are as follows: x is the main flow direction with x = 0 at the beginning of the ice cover; y is the lateral direction with y = 0 at the sidewall and z is the vertical direction counting from the channel bed. The velocity was measured at the cross section with x = 9 m, which is sufficiently far from the ice cover entry to form a fully developed flow. In the measured cross section, 14 measured lines are arranged, where the measured points are also arranged with a uniform vertical distance 1 cm. The detail layout is shown in Figure 3.

**Figure 3.** Layout of the measured cross section, lines and points. **Figure 3.** Layout of the measured cross section, lines and points.

A 3D acoustic Doppler velocimeter (ADV) with a precision of ±0.25 cm/s was adopted to measure the instantaneous velocity. The maximum sampling frequency was 50 Hz, and the sampling time for every measured point was 120 s, which resulted in 6000 velocity data for each point. Five experimental conditions were considered here, and the variables among these five cases were the cover conditions and water depth *H*, which changes from 15 to 20 cm. The details of the experimental parameters are listed in Table 1. A 3D acoustic Doppler velocimeter (ADV) with a precision of ±0.25 cm/s was adopted to measure the instantaneous velocity. The maximum sampling frequency was 50 Hz, and the sampling time for every measured point was 120 s, which resulted in 6000 velocity data for each point. Five experimental conditions were considered here, and the variables among these five cases were the cover conditions and water depth *H*, which changes from 15 to 20 cm. The details of the experimental parameters are listed in Table 1.

**Table 1.** Basic details of the characteristic parameters. is the width of the flume, is the water depth, is the bed slope and is the Reynolds number. **Table 1.** Basic details of the characteristic parameters. *B* is the width of the flume, *H* is the water depth, *S* is the bed slope and *Re* is the Reynolds number.


### **4. Model Parameters 4. Model Parameters**

To predict the profiles of velocity, Reynolds stress and turbulence intensity in the icecovered flow, the model parameters, i.e., and , and , , etc., should be given first. To predict the profiles of velocity, Reynolds stress and turbulence intensity in the ice-covered flow, the model parameters, i.e., *m<sup>b</sup>* and *m<sup>i</sup>* , *n<sup>b</sup>* and *n<sup>i</sup>* , *K*0, etc., should be given first.
