*4.1. m<sup>b</sup> and m<sup>i</sup>*

expressed as

*4.1. and* In the single-power law for the free surface flow, exponent is affected by Darcy– Weisbach resistance coefficient , i.e., =ቀ଼ ቁ .ହ [30]. For the two-power expression, exponents and can also be related to Darcy–Weisbach resistance coefficients of channel bed and ice cover . Specifically, the value of one exponent can be related to the resistance coefficient of one fixed boundary. Hence, exponents and can be In the single-power law for the free surface flow, exponent *m* is affected by Darcy– Weisbach resistance coefficient *f* , i.e., *m* = *κ* 8 *f* 0.5 [30]. For the two-power expression, exponents *m<sup>b</sup>* and *m<sup>i</sup>* can also be related to Darcy–Weisbach resistance coefficients of channel bed *f<sup>b</sup>* and ice cover *f<sup>i</sup>* . Specifically, the value of one exponent *m* can be related to the resistance coefficient of one fixed boundary. Hence, exponents *m<sup>b</sup>* and *m<sup>i</sup>* can be expressed as

$$\begin{cases} \begin{array}{c} m\_b = \kappa \left( \frac{8}{f\_b} \right)^{0.5} \\ m\_i = \kappa \left( \frac{8}{f\_i} \right)^{0.5} \end{array} \tag{15} \\\ m\_i = \kappa \left( \frac{8}{f\_i} \right)^{0.5} \end{cases} \tag{15}$$

⎪ ⎨ =൬<sup>8</sup> ൰ The resistance coefficients can be calculated as [34]

$$\begin{cases} \begin{array}{c} f\_b = \frac{8n\_b^2 g}{R\_b^{1/3}} \\\ f\_i = \frac{8n\_i^2 g}{R\_i^{1/3}} \end{array} \end{cases} \tag{16}$$

where *n* denotes Manning's roughness coefficients; *R* indicates the hydraulic radiuses and the subscript *b* and *i* in them represent the channel bed and ice cover. *g* is the gravity acceleration. Here, the ratio of channel width to water depth was approximately greater than 5, which indicates that the flow could be considered a shallow flow. Hence, the hydraulic radius of each sublayer could be simplified as the corresponding sublayer depth, i.e., *R<sup>b</sup>* = *h<sup>b</sup>* and *R<sup>i</sup>* = *h<sup>i</sup>* . The flow depths of the lower bed layer and upper ice layer have been defined from the location of the maximum velocity and are represented by *h<sup>m</sup>* and *H* − *hm*. Hence, we easily obtained

$$\begin{cases} \quad h\_m = h\_b\\ \quad H - h\_m = h\_i \end{cases} \tag{17}$$

Substituting Equation (16) into Equation (15), we obtained

$$\begin{cases} \begin{aligned} m\_b &= \kappa \left( \frac{R\_b^{1/3}}{n\_b^2 g} \right)^{0.5} \\ m\_i &= \kappa \left( \frac{R\_i^{1/3}}{n\_i^2 g} \right)^{0.5} \end{aligned} \tag{18}$$

Then, we substituted Equations (17) and (18) into Equation (3) ( *<sup>h</sup><sup>m</sup> <sup>H</sup>* = *mi mi*+*m<sup>b</sup>* ) and obtained

$$h\_b = \frac{H n\_b (H - h\_b)^{1/6}}{n\_b (H - h\_b)^{1/6} + n\_l h\_b^{1/6}} \tag{19}$$

Once Manning's coefficients *n<sup>b</sup>* and *n<sup>i</sup>* are known, the depth of the lower bed layer (*h<sup>b</sup>* ) is iteratively solved. Then, exponents *m<sup>b</sup>* and *m<sup>i</sup>* were calculated by Equation (18).

*4.2. n<sup>b</sup> and n<sup>i</sup>*

The basic calculation model of Darcy–Weisbach friction coefficient *f* is

$$f = \frac{8u\_\*^2}{\mathcal{U}^2} \tag{20}$$

where *U*(= *Q*/*BH*) is the averaged velocity in a cross section. Substituting Equation (16) into Equation (20), we obtained

$$m = \frac{\mu\_\*}{\mathcal{U}\sqrt{\mathcal{g}}} \mathcal{R}^{1/6} \tag{21}$$

For the ice-covered flow, the flow is divided as two layers, and two depth-averaged velocities, *U<sup>b</sup>* and *U<sup>i</sup>* , are averaged from the channel bed and ice cover to the position of the maximum velocity. Hence, Manning's coefficient in each sublayer can be expressed as

$$\begin{cases} \begin{array}{c} n\_b = \frac{\mu\_{\ast b}}{\Pi\_b \sqrt{\mathcal{S}}} h\_b^{1/6} \\ n\_i = \frac{\mu\_{\ast i}}{\Pi\_i \sqrt{\mathcal{S}}} h\_i^{1/6} \end{array} \tag{22}$$

These two Manning's parameters can be determined by the measured velocity profile. For the fully developed asymmetric flow, the method of logarithmic law can be adopted to predict the vertical distribution of the longitudinal velocity. However, because the velocity gradient of the maximum velocity is discontinuous, the method of logarithmic law may not be completely feasible for the entire water depth. For the ice-covered flow,

the method of logarithmic law can be applied in each sublayer, i.e., the lower bed layer and upper ice layer. Bonakdari et al. (2008) [35] further demonstrate that the turbulent boundary layer consists of the inner region near the sidewall and outer region far from the sidewall and propose that the flow velocity inside the inner region can better satisfy the distribution form of the logarithmic law. In this study, we define that the inner region inside the lower bed layer starts from the channel bed (*z* = 0) to 0.2*H*, and the inner region inside the upper ice layer is from 0.8*H* to the ice cover (*H*). Therefore, we used the following logarithmic law to describe the velocity profile in the two inner regions [36]

$$\begin{cases} \frac{\underline{lI}\_b}{\underline{u}\_{\ast b}} = \frac{1}{\underline{k}} \ln \left( 30 \frac{\underline{z}}{\underline{k}\_{\text{sb}}} \right) \\\ \frac{\underline{lI}\_i}{\underline{u}\_{\ast i}} = \frac{1}{\underline{k}} \ln \left( 30 \frac{\underline{z}}{\underline{k}\_{\text{si}}} \right) \end{cases} \tag{23}$$

where *ksb* and *ksi* are the roughness heights of the turbulence boundary in the lower bed layer and upper ice layer. The logarithmic law can be simplified to

$$\begin{cases} \, \mathcal{U}\_b = a\_b \ln(z) + F\_b \\ \, \mathcal{U}\_i = a\_i \ln(z) + F\_i \end{cases} \tag{24}$$

where *<sup>a</sup><sup>b</sup>* = *<sup>u</sup>*∗*b*/*κ*, *<sup>a</sup><sup>i</sup>* = *<sup>u</sup>*∗*i*/*κ*, *<sup>F</sup><sup>b</sup>* = *<sup>a</sup><sup>b</sup> ln* <sup>30</sup> *ksb* and *F<sup>i</sup>* = *a<sup>i</sup> ln* <sup>30</sup> *ksi* . These parameters can be obtained by regression analysis based on the measured longitudinal velocity data; then, shear velocities *u*∗*<sup>b</sup>* and *u*∗*<sup>i</sup>* and roughness heights *ksb* and *ksi* can be obtained.

*4.3. K*<sup>0</sup>

According to the previous work [28,37], *K*<sup>0</sup> is the flow model parameter for a known flow rate per unit flow width and is

$$K\_0 = \mathcal{U}/K\_1\tag{25}$$

where *K*<sup>1</sup> is the normalized depth-averaged velocities for total flow and is calculated by

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

### **5. Results and Discussion**

### *5.1. Model Verification*

With the above model parameters, the theoretical model can be applied to the experimental data. Figure 4 compares the predicted velocity from the model with the measured ones. Both profiles followed the assumed one, where the maximum velocity occurs near the middle of the water depth, and the minimum velocity appears near the fixed boundaries. The analytical results were consistent with the experimental velocities, so the proposed model and corresponding parameters are reasonable and reliable to be used to predict the velocity profile in the straight open channel flow covered by the ice cover.

Figure 5 presents the profiles of measured Reynolds stress and analytical Reynolds stress for all experimental cases. The Reynolds stress in the longitudinal direction demonstrates a linear distribution along the water depth. Zero shear stress occurred at the location of *h<sup>τ</sup>* near the middle of the water depth, after which the absolute value of Reynolds stress gradually increased and reached each peak near the fixed boundaries. The measured and predicted data in Figure 5 were basically consistent and had identical trends, which indicates that our proposed model for Reynolds stresses was feasible.

**Figure 4.** Comparison of the measured and analytical velocities. Black squares denote the measured data and red lines denote the analytical ones. **Figure 4.** Comparison of the measured and analytical velocities. Black squares denote the measured data and red lines denote the analytical ones.

Figure 5 presents the profiles of measured Reynolds stress and analytical Reynolds stress for all experimental cases. The Reynolds stress in the longitudinal direction demonstrates a linear distribution along the water depth. Zero shear stress occurred at the location of ℎఛ near the middle of the water depth, after which the absolute value of Reynolds stress gradually increased and reached each peak near the fixed boundaries. The measured and predicted data in Figure 5 were basically consistent and had identical trends,

which indicates that our proposed model for Reynolds stresses was feasible.

**Figure 5.** Comparison of the measured and analytical Reynolds stresses. Black squares denote the measured data and red lines denote the analytical ones. **Figure 5.** Comparison of the measured and analytical Reynolds stresses. Black squares denote the measured data and red lines denote the analytical ones.

Figure 6 shows the results of the turbulence intensity from the experiments and analytical model. The turbulence intensity reached its minimum near the middle water depth,

who demonstrate that the turbulence production near the central region is small, and its diffusion effect is significant, but the turbulence production near the fixed boundaries reaches its maximum value. Overall, the analytical model can catch the general trend of the measured turbulence intensity, although the measured one has some fluctuations, which confirms that the model can be applied to calculate the turbulence intensity.

(**a**) Case 1 (**b**) Case 2

lines denote the analytical ones.

Figure 6 shows the results of the turbulence intensity from the experiments and analytical model. The turbulence intensity reached its minimum near the middle water depth, from which it had an increasing trend towards the channel bed and ice cover. The profiles of turbulence intensity are consistent with the results from Papanicolaou et al. (2007) [38], who demonstrate that the turbulence production near the central region is small, and its diffusion effect is significant, but the turbulence production near the fixed boundaries reaches its maximum value. Overall, the analytical model can catch the general trend of the measured turbulence intensity, although the measured one has some fluctuations, which confirms that the model can be applied to calculate the turbulence intensity.

*Water* **2021**, *13*, x FOR PEER REVIEW 12 of 17

To find the difference between analytical and measured data, the error analysis was conducted. An absolute error is defined as the difference between analytical and measured time-averaged velocities. Hence, the average absolute error is calculated as

$$\overline{\Delta\_a} = \frac{1}{N} \sum\_{1}^{N} \left| (I)\_{analytic} - (I)\_{measured} \right| \tag{27}$$

where *N* is the number of measured points in each measured line for each case, (*I*)*analytical* and (*I*)*measured* are the analytical and measured values and *I* represents the variables, i.e., longitudinal velocity, Reynolds stress and turbulence intensity. (**e**) Case 5

The average relative error is defined as **Figure 5.** Comparison of the measured and analytical Reynolds stresses. Black squares denote the measured data and red

$$\overline{\Delta\_r} = \frac{1}{N} \sum\_{1}^{N} \left| \frac{(I)\_{\text{analytical}} - (I)\_{\text{measured}}}{(I)\_{\text{measured}}} \right| \times 100\% \tag{28}$$

As shown in Table 2, the time-averaged longitudinal velocities obtained from the proposed model were reliable within an accuracy of 0.019 m/s in terms of the average absolute error ∆*a*. The average relative error ∆*<sup>r</sup>* was 2.86–10.97%. The average absolute error and average relative error of the Reynolds stress were within 0.093 m/s and 13.63%, respectively. The average absolute error and relative error of the turbulence intensity were within 0.0019 m/s and 12.54%, respectively. All values in Table 2 were below 20%, which further confirmed that the proposed analytical model was reliable and feasible to predict the velocity and turbulence structure in the ice-covered flow. lytical model. The turbulence intensity reached its minimum near the middle water depth, from which it had an increasing trend towards the channel bed and ice cover. The profiles of turbulence intensity are consistent with the results from Papanicolaou et al. (2007) [38], who demonstrate that the turbulence production near the central region is small, and its diffusion effect is significant, but the turbulence production near the fixed boundaries reaches its maximum value. Overall, the analytical model can catch the general trend of the measured turbulence intensity, although the measured one has some fluctuations, which confirms that the model can be applied to calculate the turbulence intensity.

**Figure 6.** *Cont.*

**Figure 6.** Comparison of the measured and analytical turbulence intensity. Black squares denote the measured data and red lines denote the analytical ones. **Figure 6.** Comparison of the measured and analytical turbulence intensity. Black squares denote the measured data and red lines denote the analytical ones.


To find the difference between analytical and measured data, the error analysis was conducted. An absolute error is defined as the difference between analytical and meas-**Table 2.** Error statistics for the longitudinal velocity, Reynolds stress and turbulence intensity. ∆*<sup>a</sup>* is the average absolute error and ∆*r* is the average relative error.

## *5.2. Discussion*

∆ തതത= <sup>1</sup> ቤ()௬௧ <sup>−</sup> ()௦௨ௗ 5.2.1. Manning's Coefficients *n<sup>b</sup>* and *n<sup>i</sup>*

The average relative error is defined as

()௦௨ௗ <sup>ቤ</sup> ଵ As shown in Table 2, the time-averaged longitudinal velocities obtained from the proposed model were reliable within an accuracy of 0.019 m/s in terms of the average absolute error ∆ തതത. The average relative error ∆ തതത was 2.86–10.97%. The average absolute Table 3 lists the calculated Manning's roughness coefficients of the channel bed and ice cover for each case, i.e., *n<sup>b</sup>* and *n<sup>i</sup>* . Manning's roughness coefficient of the channel bed was 0.012–0.015, and its mean value was 0.0138 with the standard deviation of 0.0012, which verified that *n<sup>b</sup>* hardly changed among all studied cases. Manning's roughness coefficients of the ice cover in all cases were 0.017–0.02, which had small fluctuations and were lightly

× 100% (28)

ே

larger than those of the channel bed. The mean value of *n<sup>i</sup>* was 0.0182 with a standard deviation of 0.0012. For the experimental channel, Manning's roughness coefficients of the channel bed and ice cover should be considered two specific constants. We took the mean Manning's roughness coefficient in each layer as the final Manning's roughness coefficient. Hence, *n<sup>b</sup>* was set to be 0.0138, and *n<sup>i</sup>* was equal to 0.0182.


**Table 3.** Values of the model parameters to predict the velocity and Reynolds stress.

5.2.2. Flow Parameters *m<sup>b</sup>* and *m<sup>i</sup>*

The values of exponents *m<sup>b</sup>* and *m<sup>i</sup>* in this study are shown in Table 3. All values of *m<sup>b</sup>* were 5.4 ≤ *m<sup>b</sup>* ≤ 7.1, and all *m<sup>i</sup>* were 4.3 ≤ *m<sup>i</sup>* ≤ 5.3. Teal et al. (1994) [27] estimated *m<sup>b</sup>* and *m<sup>i</sup>* by nonlinear regression for more than measured 2300 vertical velocity profiles, and they found that these two parameters were 1.5–8.5, which includes the theoretical range (see Table 3). For an open channel flow, only exponent *m<sup>b</sup>* is considered and is approximately 6–7, which is unlike the values used here. Hence, both *m<sup>b</sup>* and *m<sup>i</sup>* in the covered flow are affected by the roughness characteristics of both the channel bed and ice cover.

The shape of the vertical profile of the longitudinal velocity is determined by exponents *m<sup>b</sup>* and *m<sup>i</sup>* . The ratio of them can be given from Equation (17) as

$$\frac{m\_i}{m\_b} = \left(\frac{h\_i}{h\_b}\right)^{1/6} \frac{n\_b}{n\_v} \tag{29}$$

According to Equation (29), the ratio of *m<sup>b</sup>* and *m<sup>i</sup>* is mainly influenced by the roughness coefficients of the channel bed and ice cover. This result is also confirmed by the measured data of Li et al. (2020) [29]. They demonstrate that the vertical distribution of the velocity remains constant with changing water depth and flow rate under the same *m<sup>b</sup>* and *m<sup>i</sup>* or the same ratio of *m<sup>b</sup>* and *m<sup>i</sup>* . In the asymmetric flow, the maximum velocity tended to be closer to the smooth boundary with smaller roughness coefficient. In Table 3, the channel bed had a smaller roughness coefficient than the ice cover, which corresponded to the close maximum velocity to the channel bed, which was verified by *h<sup>m</sup>* < *H* − *hm*. In general, both open channel flow and symmetry flow can be considered as the special cases of asymmetric flow. For the open channel flow, the exponent *m<sup>i</sup>* tends to infinity. For the symmetry flow, *m<sup>b</sup>* = *m<sup>i</sup>* .
