Experimental and Numerical Determination of the Head Loss of a Pressure Driven Flow through an Unlined Rock-Blasted Tunnel

The friction loss in a part of the rock-blasted unlined tunnel of the Litjfossen hydropower plant in Norway was determined from experimental and numerical studies. Remote sensing data from the prototype tunnel provided the input data for both the numerical model and the construction of a 1:15 scale model with an innovative milling approach. The numerical simulations were based on the solution of the Reynolds-averaged Navier–Stokes equations using the CFD program OpenFoam. Head loss measurements in the scale model were carried out by means of pressure measurements for a range of discharges and were compared against the results of the numerical model. The measured data were used to determine the Darcy–Weisbach and Manning friction factors of the investigated tunnel reach. The high-resolution remote sensing data were also used to test the applicability of existing approaches to determine the friction factor in unlined rock blasted tunnels. The results of the study show the usefulness of the chosen hybrid approach of experimental investigations and numerical simulations and that existing approaches for the determination of head losses in unlined tunnels need to be further refined.


Introduction
Hydropower is the principal source of Norway's electricity supply with Norway being the world's seventh-largest hydropower producer and the largest in Europe in 2019 [1]. There exist more than 1000 hydropower reservoirs in Norway [2] including 4000 km of hydropower tunnels of which 96-98% are unlined [3,4]. Moreover, more than 50 embankment dams in Norway have closed conduit (shaft/tunnel) spillways, typically constructed as rock-blasted tunnels, which are used for flood protection. Compared to lined tunnels, unlined tunnels have a larger roughness due to irregularities in cross-sectional area and wall surface e.g., [5,6] which originates from the excavation method such as boring, drilling, and blasting. Despite the reduced conveyance capacity and larger energy losses, unlined tunnels are still considered to be economically advantageous in hard rocks because of reduced construction costs. Moreover, rock blasting represents a more rapid method for tunnel excavation than utilizing tunnel boring machines [7]. For all these reasons, unlined rock-blasted tunnels are a common feature in the Norwegian hydropower system [4]. Thus, there is the need to accurately forecast head losses and conveyance capacities of these waterways for the planning, design, and construction of hydropower projects. The scope of this paper is to highlight possibilities to determine the friction losses in such conduits using state-of-the-art measurement, experimental, and numerical techniques.
Frictional drag causing head losses in hydraulic waterways is typically characterized by "roughness" or "friction" factors such as Manning's n, Chezy's C, or Darcy-Weisbach's f, which are interconnected and can be converted into each other. The focus of this paper is on the Darcy-Weisbach friction factor f that is commonly used for the determination of the friction head loss h f along a certain reach length L of a pressurized conduit according to e.g., [8]: where D = 4A/P denotes the hydraulic diameter of the conduit (with A = cross-sectional area and P = wetted perimeter), V m the reach averaged flow velocity, and g the acceleration due to gravity. For pipe-flow, the friction factor f can be determined from the Moody diagram as a function of the relative roughness k/D p (i.e., the ratio of the roughness height k to pipe diameter D p ) and the Reynolds number Re = V m D p /ν [8], where ν denotes the kinematic viscosity of the fluid. For turbulent flow in rough pipes, the friction factor can be determined from [6]: For conduits with a cross-sectional shape deviating from a circular cross section, the pipe diameter D p is replaced by the hydraulic diameter D in Equation (2). The roughness height k is typically represented in terms of the equivalent sand-roughness k s . The latter is not a physical measure of wall roughness but represents the height of the uniformly distributed sand grains in Nikuradse's [9] experiments producing the same friction as the surface of interest in the fully rough regime [10]. However, unlined rock tunnels are characterized by a large variation in size and shape from one cross section to another, which results in surface irregularities e.g., [6]. Surface irregularities are the reason why the determination of k s is associated with large uncertainties in rough conduits, and despite significant research efforts in the past, available approaches linking the friction factor to flow and topographic features of the surface are largely of empirical nature that may not be theoretically justified by the underlying theory [11]. In other words, "engineers do not currently have the ability to accurately predict the drag of a generic rough surface" [12]. This is especially the case for the scale-dependent roughness of unlined rock-blasted tunnels, which can be divided into crystal roughness, break roughness, and overbreak roughness [13]. A similar classification was given by [14], who distinguished between material, surface, and cross-sectional roughness, respectively. The micro-scale crystal or material roughness depends on geological properties of the rocks such as the size of the rock crystals. The geological properties are important for the excavation process and tunnel stability (which is not in the scope of the present paper) and may, hence, indirectly affect the surface and cross-sectional roughness. For the evaluation of the friction factor, the effect of micro-roughness can be assumed to be negligible compared to the surface and cross-sectional roughness, which are typically considered as the main parameters for the energy losses. The surface (or break) roughness depends on the unevenness of the fracture surface, and overbreak roughness is defined as the deviation of the blasted tunnel profile from the theoretical profile [13], i.e., it is associated with the variation in cross-sectional area defining the cross-sectional roughness. It is worth mentioning that the tunnel invert is generally smoother than the tunnel walls and that its roughness depends on the treatment of the invert at the time of tunnel completion [14,15]. Several empirical approaches have been developed to determine the friction factor in unlined rock-blasted tunnels using field data and focusing on the fully rough regime. One of the first approaches was derived by Rahm [15], who used data from Swedish hydropower tunnels to characterize the cross-sectional roughness along a tunnel reach based on distribution curves of the cross-sectional areas by defining the parameter δ as a measure for the relative variation in the cross-sectional area: In this equation, A 99 denotes the tunnel cross-sectional area corresponding to a frequency of 99% and A 1 to 1%, respectively. These frequencies were derived by plotting the cross-sectional areas, which were determined approximately every 15 m along the tunnels, in their order of magnitude to obtain distribution curves on standard distribution paper and fitting a straight line to the distribution between A 99 and A 1 . Rahm [15] related the parameter δ to the friction factor according to: Reinius [16] found that this formula worked well for tunnels excavated before 1946 and that f = 0.02 + 0.00168δ (5) provided a better fit to data of [15] for tunnels excavated after 1946 because of the further development of the tunneling techniques. Investigations by Priha [13] regarding the friction of small unlined rock tunnels suggested that a 5-m spacing between the measured cross sections is sufficient for good results, and that the friction factor may be determined from f = 0.0033δ A 1 A 1 + 9 (6) An approach combining both surface and cross-sectional roughness was developed by Rønn and Skog [14]. The so-called IBA method, where IBA denotes the abbreviation of the author's university institute, is based on the analysis of longitudinal and cross-sectional profiles. The surface (or wall) roughness is described by the standard deviation s of at least three wall-profiles parallel to the tunnel axis (right wall, left wall, and roof). The profiles should be sampled over a length of 20-25 m at constant intervals of 0.25-0.5 m, and the final surface roughness s w is determined from [14]: where m denotes the number of longitudinal profiles. Analogously, the cross-sectional roughness is parameterized in terms of the standard deviation of the square roots of the areas: where n denotes the number of measured cross-sectional areas, A the mean cross-sectional area, and the factor 0.53 was given by the authors without further explanations. For the application of this approach, the cross-sectional areas should be determined at intervals between 0.5 and 1.0 m, and n > 50 measurements should be taken over a length of at least 25 m. If there are data for more sections available, the mean value of the cross-sectional roughness may be determined using the same formulation as in Equation (7). The final roughness to be used in Equation (2) is k = s w + s A .
The approaches presented above were specifically developed for unlined rock-blasted tunnels used in hydropower systems. The approaches focusing solely on cross-sectional roughness are based on manual measurements of tunnel topography, and the IBA-method may be applied using both manually collected data or data collected using state-of-the-art remote sensing instruments and methods (e.g., [17,18]). Recently, Watson, and Marshall [19] introduced a bi-directional method (which they named Queen's method) to estimate the friction factor in underground mine tunnels used for ventilation. This approach has been specifically designed to exploit point clouds determined from modern remote sensing instruments and it is also based on the determination of surface and cross-sectional roughness. For this purpose, the cross-sectional variance σ n 2 is calculated based on n cross-sectional scans and m longitudinal points [19]: where h i,j denotes the distance between a zero line parallel to the tunnel axis and the tunnel wall. The longitudinal variance is calculated similarly: The final roughness predictor to be used in (2) is then Further methods and studies exist dealing with the determination of roughness and friction factors of tunnels, which are partly based on similar statistical parameters and/or spectral analyses. However, most of these methods have been developed for other tunnel types such as unlined and shotcreted bored tunnels e.g., [5,20,21] and references therein. These methods are not applied here because of the focus of the present study on the friction losses in unlined rock-blasted tunnels.
An important aspect for the assessment of approaches for the determination of friction factors is the availability of validation data. Field data were used for the derivation of the approaches of [13][14][15][16]22], and the friction factors f reported in these studies varied between 0.045 and 0.116. These studies were carried out before the availability of modern remote sensing methods, which becomes also apparent from the spatial resolution of the field measurements. More recent data are sparse as field measurements of hydraulic parameters in rock-blasted tunnels in prototype conditions are demanding and complicated. Hydraulic scale models, on the other hand, represent an efficient way to determine friction losses in a controlled environment e.g., [23,24]. However, the complex topography of unlined rock-blasted tunnels is difficult to replicate in scale models and this is why many studies have been carried out using simplified roughness patterns e.g., [13,16,25]. A common approach used in Norway to model the cross-sectional roughness in physical models of rock-blasted tunnels has been to mimic the tunnel shape by smooth acrylic glass and the variation in cross-sectional geometry by strip-roughness (rubber bands glued into the tunnels) with a spacing so that the model friction factor is similar to an assumed prototype friction factor [25]. The physical roughness characteristics of strip-roughness is different from the prototype wall roughness, and these differences will result in a different near bed flow pattern which in turn affects hydraulic resistance. The assumption that such a roughness pattern represents the roughness associated with the blasting process is, hence, a limitation but is required before the excavation of the tunnel due the unknown geometry. Thus, such scale models are prone to both laboratory and scale effects e.g., [24,26] which may be minimized by further developing existing technologies for scale model production to avoid simplifying assumptions regarding wall roughness characteristics. This becomes possible through the use of computer numerically controlled (CNC) production techniques e.g., [27,28] and a recent study of [29] showed the possibilities to gain further insights into roughness effects in non-uniform hydraulic waterways of hydropower systems using this technique.
A further possibility to investigate the flow features and hydraulic losses in conduits is by means of computational fluid dynamics (CFD) which has become more and more popular because of the increased computer power over time. Simulations carried out by [30] using a computational grid derived from terrestrial laser scanning (TLS) data of a hydropower-tunnel and a 3D-modelling approach based on the Reynolds-averaged Navier-Stokes (RANS) equations with the k-ε turbulence model showed the possibility to study local pressure variations due to a complex tunnel geometry. A similar conclusion was reported based on the numerical simulations of physical scale model experiments by [29]. However, it should be kept in mind that CFD models require a rigorous verification and validation process. This calls for a closer inspection of the applicability of high-resolution CFD models to determine the friction factor of unlined rock-blasted tunnels with complex topography for which the actual friction factor is unknown.
This study deals with the determination of friction losses in an unlined rock-blasted tunnel. It was carried out in the framework of a research project which has the objective to link physical wall roughness of unlined rock-blasted hydropower tunnels to hydraulic resistance making use of the latest technological advances. Using TLS data of a straight tunnel section of the Litjfossen hydropower plant in Norway, friction factors were determined from both existing relationships and a scale model study. Moreover, the results from the scale model study were used to investigate the applicability of a RANS model to replicate the head losses in the tunnel. Section 2 presents the TLS data as well as the experimental and numerical setup and methods. The results from the scale model experiments are presented in Section 3 along with further considerations on the roughness geometry and the calibration of the numerical model. The results are discussed in Section 4 and the final conclusions and an outlook are given in Section 5.

Survey Data and Digital Elevation Model
This study is based on remote sensing data of a 230-m long straight and unlined rock-blasted hydropower tunnel reach belonging to the Litjfossen hydropower plant in Norway. A straight part of the tunnel system was chosen for the experimental and numerical investigations to minimize energy losses emerging from a curved tunnel planform. Point-clouds were measured with a Leica C10 Terrestrial Laser Scanner (TLS) in August 2016 in dry conditions during a maintenance period for the hydropower station. The surveyed tunnel reach had a small longitudinal slope, which helped to drain puddles of water on the invert so that corresponding disturbances were minimized. For the scan, the TLS was placed on a total of 19 positions along the 230 m long target section. Seven targets along the scanned tunnel section were used to establish a local coordinate system and for the merging of the overlapping point clouds. Point set registration and post-processing of the point-clouds were carried out using the Leica Cyclone software package (version 9.1.2) and the iterative closest point algorithm. The final point cloud consisted of 42,927,547 points and had an average resolution of approx. one point per cm 2 [27,28]. The surveyed section contained alternating alcoves on the tunnel sides, used for truck crossing during the tunnel excavation and maintenance. Therefore, the selected scanned section was reduced to a 120.8-m long singularity-free straight section, and the corresponding point cloud, visualized in Figure 1, contained 24,007,488 points. This point cloud was used to generate a digital elevation model (DEM) which in turn served as the basis for the scale model construction (cf. Section 2.2) and the geometry of the CFD-simulations (cf. Section 2.4). The triangular mesh of the DEM was generated by the 3D-Reshaper software package (library version 17.0.24386.0) and consisted of approx. 1 million triangles. A sub-section of the DEM is shown in Figure 1.

Scale Model Production
The geometric scale for the experimental investigations in the hydraulic laboratory of the Norwegian University for Science and Technology (NTNU) was set to 1:15 which renders the length of the tunnel model to 8 m (120 m in prototype scale). Using the scaled DEM as input, the tunnel model was milled at the CNC milling facility at SINTEF Ocean in Trondheim, Norway, using a 5-axis milling machine type Belotti MDL 12070. Divinycell H60, a low-density polyvinylchloride (PVC) foam material consisting of 2440-mm long, 1220 wide, and 85-mm high foam plates, served as the base material for the milling. The plates were glued together to form blocks of the required height for the milling (Figure 2a). The plate length defined the length of individual longitudinal sections and the DEM was cut accordingly using the Computer Aided Manufacturing (CAD) software. Each longitudinal tunnel section was divided into three segments ( Figure 2a) to ensure that all the edges and cavities could be milled out using 20 mm and 25 mm diameter ball-nosed mills for rough milling and a ball-nosed endmill with a diameter of 10 mm for fine milling. The effective milling time was approx. 70-80 h plus rigging and handling. A polyester coating was sprayed on the milled surfaces reinforcing the surface and improving the water tightness before gluing the segments together to reconstruct the 3D roughness. Figure 2b provides a view into the assembled tunnel model.
Two windows (150 mm × 250 mm) were manufactured and installed at the invert and at the tunnel side to gain optical access into the milled tunnel for fine-scale velocity measurements with particle laser velocimetry (PIV). The invert window was made of smooth perspex plates allowing for optical access with high-speed cameras. Its effect on the head loss measurements was assumed to be negligible as the milled tunnel invert was rather smooth. The surface of the wall-window was milled from a perspex block so that the wall roughness pattern was not changed. Thin horizontal grooves were milled out of the side window to enable the formation of a laser sheet in the tunnel without any impact of the irregular window roughness [28]. The presentation of the results from the fine-scale velocity measurements are beyond the scope of the present paper and will be reported in a followup study.
For the quantification of the accuracy of the milling process, the milled tunnel walls were scanned prior to surface treatment and assembly of the tunnel sections using the same total station and TLS as for the prototype-scan. The variable reflective properties of the foam material in the

Scale Model Production
The geometric scale for the experimental investigations in the hydraulic laboratory of the Norwegian University for Science and Technology (NTNU) was set to 1:15 which renders the length of the tunnel model to 8 m (120 m in prototype scale). Using the scaled DEM as input, the tunnel model was milled at the CNC milling facility at SINTEF Ocean in Trondheim, Norway, using a 5-axis milling machine type Belotti MDL 12070. Divinycell H60, a low-density polyvinylchloride (PVC) foam material consisting of 2440-mm long, 1220 wide, and 85-mm high foam plates, served as the base material for the milling. The plates were glued together to form blocks of the required height for the milling (Figure 2a). The plate length defined the length of individual longitudinal sections and the DEM was cut accordingly using the Computer Aided Manufacturing (CAD) software. Each longitudinal tunnel section was divided into three segments ( Figure 2a) to ensure that all the edges and cavities could be milled out using 20 mm and 25 mm diameter ball-nosed mills for rough milling and a ball-nosed endmill with a diameter of 10 mm for fine milling. The effective milling time was approx. 70-80 h plus rigging and handling. A polyester coating was sprayed on the milled surfaces reinforcing the surface and improving the water tightness before gluing the segments together to reconstruct the 3D roughness. Figure 2b provides a view into the assembled tunnel model.
Two windows (150 mm × 250 mm) were manufactured and installed at the invert and at the tunnel side to gain optical access into the milled tunnel for fine-scale velocity measurements with particle laser velocimetry (PIV). The invert window was made of smooth perspex plates allowing for optical access with high-speed cameras. Its effect on the head loss measurements was assumed to be negligible as the milled tunnel invert was rather smooth. The surface of the wall-window was milled from a perspex block so that the wall roughness pattern was not changed. Thin horizontal grooves were milled out of the side window to enable the formation of a laser sheet in the tunnel without any impact of the irregular window roughness [28]. The presentation of the results from the fine-scale velocity measurements are beyond the scope of the present paper and will be reported in a follow-up study.
For the quantification of the accuracy of the milling process, the milled tunnel walls were scanned prior to surface treatment and assembly of the tunnel sections using the same total station and TLS as for the prototype-scan. The variable reflective properties of the foam material in the infrared spectrum resulted in biased scans of the right-hand side tunnel segments (cf. Figure 2a), so that these scans were not further analyzed. The left-hand side tunnel segments, on the other hand, could be accurately scanned. The obtained point-clouds had a density of about 70 points per cm 2 (model scale) and could be registered with an accuracy close to the millimeter [28].
Water 2020, 12, x FOR PEER REVIEW 7 of 23 infrared spectrum resulted in biased scans of the right-hand side tunnel segments (cf. Figure 2a), so that these scans were not further analyzed. The left-hand side tunnel segments, on the other hand, could be accurately scanned. The obtained point-clouds had a density of about 70 points per cm 2 (model scale) and could be registered with an accuracy close to the millimeter [28].

Experimental Setup and Methods
The tunnel parts were assembled at the experimental stand on a special test-rig ( Figure 3). Water was recirculated through the tunnel by a centrifugal pump, which was installed at the downstream end of the tunnel section. The discharge was measured by an inductive discharge meter in the pipesection of the experimental stand (Siemens MAG 5000W, accuracy 0.2%). A honeycomb panel was installed in the pipe leading to the tunnel inlet and outlet to condition the flow. The transitioning inlet and outlet sections connecting the pipes to the milled tunnel were custom-made to guarantee a smooth transition from the circular pipe cross section to the horseshoe-shaped tunnel cross section. Pressure measurements were carried out at eight cross sections along the tunnel to determine the pressure gradient and, hence, the friction factors. The distance between the cross sections was 1 m each. A total of 17 pressure taps were installed around the peripheral of each cross section to

Experimental Setup and Methods
The tunnel parts were assembled at the experimental stand on a special test-rig ( Figure 3). Water was recirculated through the tunnel by a centrifugal pump, which was installed at the downstream end of the tunnel section. The discharge was measured by an inductive discharge meter in the pipe-section of the experimental stand (Siemens MAG 5000W, accuracy 0.2%). A honeycomb panel was installed in the pipe leading to the tunnel inlet and outlet to condition the flow. The transitioning inlet and outlet sections connecting the pipes to the milled tunnel were custom-made to guarantee a smooth transition from the circular pipe cross section to the horseshoe-shaped tunnel cross section.
Water 2020, 12, x FOR PEER REVIEW 7 of 23 infrared spectrum resulted in biased scans of the right-hand side tunnel segments (cf. Figure 2a), so that these scans were not further analyzed. The left-hand side tunnel segments, on the other hand, could be accurately scanned. The obtained point-clouds had a density of about 70 points per cm 2 (model scale) and could be registered with an accuracy close to the millimeter [28].

Experimental Setup and Methods
The tunnel parts were assembled at the experimental stand on a special test-rig ( Figure 3). Water was recirculated through the tunnel by a centrifugal pump, which was installed at the downstream end of the tunnel section. The discharge was measured by an inductive discharge meter in the pipesection of the experimental stand (Siemens MAG 5000W, accuracy 0.2%). A honeycomb panel was installed in the pipe leading to the tunnel inlet and outlet to condition the flow. The transitioning inlet and outlet sections connecting the pipes to the milled tunnel were custom-made to guarantee a smooth transition from the circular pipe cross section to the horseshoe-shaped tunnel cross section. Pressure measurements were carried out at eight cross sections along the tunnel to determine the pressure gradient and, hence, the friction factors. The distance between the cross sections was 1 m each. A total of 17 pressure taps were installed around the peripheral of each cross section to Pressure measurements were carried out at eight cross sections along the tunnel to determine the pressure gradient and, hence, the friction factors. The distance between the cross sections was 1 m each. A total of 17 pressure taps were installed around the peripheral of each cross section to measure the average static pressure akin to the principle of a piezometer-ring. Piezometer-rings are typically used in experiments at ducts characterized by pressure variations around the perimeter of the flow cross section e.g., [31], as can be expected in unlined tunnels because of both the irregular cross-sectional geometry and the surface-roughness. Using such a setup, it was assumed that the representative cross-sectionally averaged pressure can be determined without significant bias due to the surface irregularities. This assumption will be evaluated using the results of the numerical simulations in the discussion (Section 4). Differential pressure sensors (Aplisens APR-2000ALW, accuracy 0.075%) were used for the measurements. The pressure at each cross section was measured against a "reference pressure" formed by a constant water level in a small tank which was connected to the tunnel model. This constant water level was monitored using an acoustic sensor (Microsonic, accuracy 1%). The pressure at each section was measured for 500 s with a sampling frequency of 100 Hz.
The pressure measurements were carried out for nine steady discharges (35 l/s-113 l/s). An air valve installed in the model was used to evacuate the trapped air before data collection. Using the geometrical scale of 1:15 and similarity in the Euler number, which provides the same scale ratios as the Froude similarity for the problem at hand [32], these discharges correspond to 30.5 m 3 /s-98.5 m 3 /s in prototype conditions. We note that a discharge of 30 m 3 /s is the regular discharge of the hydropower station and that larger discharges were used to increase the head losses along the tunnel to determine the friction factor (cf. Section 3). The pressure measurements were repeated four times at different days. The observed variation in discharge was small (average standard deviation of 0.0002 m 3 /s) and the standard deviation of the mean pressure for the tested discharges at the measurement cross sections ranged between 0.0002 m and 0.0015 m (mean value 0.0007 m). These values indicate a high degree of reproducibility and the head loss considerations in Section 3 are hence based on the average of the four measurements.
Head losses were calculated based on the Bernoulli equation and the equation of continuity according to: where h f,i−j denotes the considered reach between cross-sections i and j, p the average pressure at the cross section, ρ the fluid density, Q the discharge, and α the kinematic correction factor. Because of the unknown velocity distribution at the measurement cross sections α = 1 was assumed. For the determination of the friction factor Equation (1) was rearranged to yield in which D i−j denotes the average hydraulic diameter over the reach i−j (determined from the DEM-analysis), x i,j the position of the cross sections i and j along the tunnel axis, respectively, and V m,i−j = Q/A i−j (with A i−j = reach averaged cross-sectional area) the reach averaged flow velocity.
The term h f,i−j /(x i − x j ) represents the slope of the energy line which can be determined by a linear regression analysis of the available head data instead of using only two local values.

Numerical Model
The numerical simulations were carried out using the open-source CFD model OpenFoam (v.1906) [33]. The program solved the RANS equations to find the water velocities in the tunnel. The equations were solved in three dimensions using a finite volume method and an unstructured grid with dominantly hexahedral cells. The use of hex-cells allowed the alignment of the grid with the flow direction and reduction of false diffusion. Tetrahedral cells were used close to the tunnel walls to describe the complex geometry. The pressure was computed with the Semi-Implicit Method for Pressure-Linked Equations (SIMPLE algorithm) [34]. The method is based on a guessed pressure field where a pressure-correction equation is used. The source term in the equation is the water continuity deficit in each cell. The turbulence was computed from the k-ε model, which is a two-equation model that calculates the turbulent kinetic energy, k, and its dissipation, ε, from two partial differential equations. The turbulent eddy viscosity is then computed from these two parameters. The OpenFOAM kqRWallFunction [35] was used as boundary condition at the wall for the k equation and the epsilonWallFunction [35] was used for ε.
This type of numerical model has been used in a large number of CFD projects over the last 50 years and is well-known. The main difference between the standard model and what was used in the current project was the treatment of the tunnel walls. The wall friction for rough walls was based on the nutkRoughWallFunction in OpenFOAM [35] where a reduced eddy viscosity ν Tw close to the wall was used: where E is a constant equal to 9.0 and y + = u * y/ν is the dimensionless distance to the wall (u * is the shear velocity and y is the wall distance). The parameter f n is given by: with k s + = u * k s /ν. C s is an empirical parameter chosen to be 0.5 in the current study. The roughness height k s is given in meters. Equation (14) is used for k s + > 2.25, and for k s + ≤ 2.25 f n is set to 1.0.
Equation (15) is only used for k s + < 90. For k s + > 90, the following equation is used instead: The eddy-viscosity in Equation (14) was further limited to be between 0.5 ν T and 2.0 ν T , where ν T is the maximum of ν and ν Tw .
The functions were used to account for the small-scale roughness on the surfaces of the computational cells. The largest roughness elements were resolved in the grid, where their effect on the head loss was computed directly by the solution of the RANS equations and the k-ε model.
The grid was generated using the blockMesh/snappyHexMesh method in OpenFOAM. This method allowed dominantly orthogonal hexahedral cells to be used in most of the computational domain, increasing stability and accuracy. The snappyHexMesh program reduced the hex-cells close to the boundary to tetrahedrals in order to capture the complex geometry (Figure 4b). The input for snappyHexMesh was the same DEM that was used for milling the physical model. The grid domain included the inlet and outlet sections as shown in Figure 4a. Following a grid refinement study (cf. Section 4), it was decided to use a basic grid size of 5 mm (with a grid size close to the wall of 1.25 mm), so that the final grid consisted of approx. 96 million cells.
The Gauss linearUpwind scheme, a second-order upwind-biased scheme, was used for discretization of the convective terms in the RANS equations. The limitedLinear scheme was used for the turbulence equations. This is a second-order central scheme with a limiter that prevents interpolated values to exceed the bounds of the surrounding cells.

Results
This section presents the results from the analysis of the accuracy of the milled model followed by the presentation and analysis of the geometrical data, the pressure measurements and the numerical simulations. The results are discussed in Section 4.

Milling Accuracy
The accuracy of the milled model was assessed in two steps. The first step consisted in the comparison of the original point-cloud with the generated DEM. This comparison revealed that 99.5% of the cloud points were less than 0.02 m away from the mesh in prototype dimensions (corresponding 0.0013 m in model dimensions) [27,28]. The DEM served as input for both the milling and the numerical model, and the data from the scan of the milled left-hand side tunnel walls were compared with the corresponding part of the DEM in the second step. This comparison showed that 94.6% of the scanned points were less than 3 mm away from the scaled DEM or, more precisely, 71% of these points were in the range ±1.5 mm from the scaled DEM [27,28]. The accuracy of the milled model was considered to be adequate for the subsequent hydraulic tests as discussed in Section 4.1.

DEM Analyses
The dimensions of the milled tunnel were determined from the DEM. The cross-sectional areas Ai of the eight measurement cross sections ranged between 0.1223 m 2 and 0.1447 m 2 (model dimensions), and the corresponding hydraulic diameters Di ranged between 0.346 m and 0.389 m ( Table 1). These values served as input for the determination of the head losses according to Equation (12). The reach averaged values of the cross-sectional area A and hydraulic diameter D were required for the determination of the friction factor according to Equation (13). The corresponding values are reported in Table 1 for two different reach lengths between the cross sections 1-8 and 2-7, respectively. The differences in the reach averaged cross-sectional areas

Results
This section presents the results from the analysis of the accuracy of the milled model followed by the presentation and analysis of the geometrical data, the pressure measurements and the numerical simulations. The results are discussed in Section 4.

Milling Accuracy
The accuracy of the milled model was assessed in two steps. The first step consisted in the comparison of the original point-cloud with the generated DEM. This comparison revealed that 99.5% of the cloud points were less than 0.02 m away from the mesh in prototype dimensions (corresponding 0.0013 m in model dimensions) [27,28]. The DEM served as input for both the milling and the numerical model, and the data from the scan of the milled left-hand side tunnel walls were compared with the corresponding part of the DEM in the second step. This comparison showed that 94.6% of the scanned points were less than 3 mm away from the scaled DEM or, more precisely, 71% of these points were in the range ±1.5 mm from the scaled DEM [27,28]. The accuracy of the milled model was considered to be adequate for the subsequent hydraulic tests as discussed in Section 4.1.

DEM Analyses
The dimensions of the milled tunnel were determined from the DEM. The cross-sectional areas A i of the eight measurement cross sections ranged between 0.1223 m 2 and 0.1447 m 2 (model dimensions), and the corresponding hydraulic diameters D i ranged between 0.346 m and 0.389 m (Table 1). These values served as input for the determination of the head losses according to Equation (12). The reach averaged values of the cross-sectional area A and hydraulic diameter D were required for the determination of the friction factor according to Equation (13). The corresponding values are reported in Table 1 for two different reach lengths between the cross sections 1-8 and 2-7, respectively. The differences in the reach averaged cross-sectional areas A 1−8 = 0.1357 m 2 and A 2−7 = 0.1366 m 2 are small and the corresponding reach averaged hydraulic diameters D 1−8 = 0.359 m and D 2−7 = 0.360 m are almost identical. Table 1. Distance x to the beginning of the tunnel, cross-sectional area, wetted perimeter, and hydraulic diameter of the measurement cross sections (see Figure 3) and reach averaged values for the total and reduced model length. The variability of the cross-sectional area along the model tunnel is presented in Figure 5a based on the analysis of 1000 cross-sectional areas (i.e., the spacing between two subsequent cross sections was ∆ = 8 mm). The cross-sectional areas varied between 0.12 m 2 and 0.15 m 2 (27 m 2 and 34 m 2 in prototype dimensions) and several near-periodic peaks became apparent from Figure 5a which were investigated further by calculating the power spectrum shown in Figure 5b. The latter is characterized by a peak at a wavelength of about 0.23 m (ca. 3.5 m in prototype dimensions), which we associate with the blasting frequency of the tunnel.    Figure 5c shows the cumulative distribution functions for the cross-sectional areas between cross sections 1-8 and 2-7, respectively, which were used to determine A 1 and A 99 according to [15]. The corresponding values are presented in Table 2 together with the δ-values according to Equation (3). In addition, the average standard deviation s, skewness Sk, and kurtosis Ku were determined from the analysis of 360 longitudinal profiles (see Table 2). Table 2. Characteristic values derived from the distribution curves of the cross-sectional areas, standard deviation, skewness, and kurtosis for reaches 1-8 and 2-7 (in model scale).  Figure 6 presents the measured pressure head p i /(ρg) for the different discharges (see Table 3 for discharges). There is a clear separation in the pressure lines for the different discharges, and the gradient of the pressure lines increases with increasing discharge. The local pressure heads for a given discharge deviate from the straight best-fit line because of the variation in cross-sectional area (see Figure 5). This can be explained by the Bernoulli principle, i.e., for a steady flow, the pressure increases with increasing cross-sectional area and the cross-sectionally averaged velocity decreases. The largest increases in pressure head were typically observed at cross sections which were characterized by rather large increases in cross-sectional area (e.g., from 2-3, 3-4, and 6-7; see Table 1).

Head Loss Measurements
The pressure values were used to determine the total head at the eight cross sections. The results are presented in Figure 7 which includes the best fit line from the linear regression for the reach 1-8. The figure shows that, especially for the low discharges, small gradients of the energy line are observed. It needs to be mentioned that small negative head losses were determined from the experimental data between sections 3-4 (x 3 = 2.5 m and x 4 = 3.5 m). This non-physical result can be associated with the resolution and uncertainties of the measurements. However, these negative head losses were small (<0.2 mm) which is within the margin of error of the measurements.
Water 2020, 12, x FOR PEER REVIEW 13 of 23 The figure shows that, especially for the low discharges, small gradients of the energy line are observed. It needs to be mentioned that small negative head losses were determined from the experimental data between sections 3-4 (x3 = 2.5 m and x4 = 3.5 m). This non-physical result can be associated with the resolution and uncertainties of the measurements. However, these negative head losses were small (<0.2 mm) which is within the margin of error of the measurements. Friction factors were derived for the reaches between cross sections 1-8 and 2-7, respectively, using Equation (13). The 5 m long reach between cross sections 2 and 7 was considered as a compromise between reach length for the determination of friction losses and the flow development at the entrance (and accounting for outflow conditions). This issue will be discussed further in Section 4. Using reach averaged geometrical values and the mean gradient of the energy line in these sections, friction factors were determined which ranged between 0.047 and 0.056 for reach 1-8 and between 0.046 and 0.058 for reach 2-7. The Reynolds-numbers = ν m Re V D / ranged between 9.4 × 10 4 and 3.0 × 10 5 , and the largest friction factor was obtained for the lowest discharge (Table 3). Table 3. Friction factors f and head losses determined from the pressure measurements for the reaches 1-8 and 2-7, and hf,n from the five numerical simulations (between cross sections 1-8).   Table 3. Friction factors f and head losses determined from the pressure measurements for the reaches 1-8 and 2-7, and h f,n from the five numerical simulations (between cross sections 1-8).

Numerical Model
The scope of the numerical simulations covering the tunnel scale model including the transitional in-and outlet zone was to investigate the capability of the model to reproduce the head loss between cross sections 1-8. No information about the kinematic correction factor α Equation (12) was available for the scale model, so that it was also set to 1 for the calibration of the numerical model. However, the results of the simulations allow for the determination of the kinematic correction factor for each cross section, and the determined values for α ranged between 1.11 and 1.24. These values were used for additional considerations further below to investigate the importance of α for the determination of the head losses in unlined rock-blasted tunnels.
The first step of the numerical investigations consisted in a grid refinement study to identify the required grid size by quantifying the numerical uncertainty. Using a discharge of 0.0353 m 3 /s (model scale), cell sizes from 100 mm to 5 mm were tested. The best result was obtained for the grid size of 5 mm and simulations with smaller grid sizes became unstable. It was therefore decided to use the grid size of 5 mm for the further simulations. Because of the high resolution of the computational grid (approx. 96 million cells), a small roughness value ks had to be used for the simulations. A value of ks = 0.5 mm was identified during the calibration of the model, for which the difference in head loss between the physical and numerical model was less than 0.12 mm (corresponding to a difference of 3.1%). This difference was considered to be accurate enough given the achievable accuracy of the scale model tests.
The calibrated model was used to evaluate the head loss along the tunnel for four more discharges after the grid refinement study (Table 3). Figure 8 shows the comparison of the CFDresults with the scale-model data and illustrates that it was possible to reproduce the head losses between cross sections 1 and 8 for all discharges with a high degree of accuracy. For the discharges of 0.054, 0.074, and 0.093 m 3 /s, the differences between the CFD and scale model head losses were less Friction factors were derived for the reaches between cross sections 1-8 and 2-7, respectively, using Equation (13). The 5 m long reach between cross sections 2 and 7 was considered as a compromise between reach length for the determination of friction losses and the flow development at the entrance (and accounting for outflow conditions). This issue will be discussed further in Section 4. Using reach averaged geometrical values and the mean gradient of the energy line in these sections, friction factors were determined which ranged between 0.047 and 0.056 for reach 1-8 and between 0.046 and 0.058 for reach 2-7. The Reynolds-numbers Re = V m D/ν ranged between 9.4 × 10 4 and 3.0 × 10 5 , and the largest friction factor was obtained for the lowest discharge (Table 3).

Numerical Model
The scope of the numerical simulations covering the tunnel scale model including the transitional inand outlet zone was to investigate the capability of the model to reproduce the head loss between cross sections 1-8. No information about the kinematic correction factor α Equation (12) was available for the scale model, so that it was also set to 1 for the calibration of the numerical model. However, the results of the simulations allow for the determination of the kinematic correction factor for each cross section, and the determined values for α ranged between 1.11 and 1.24. These values were used for additional considerations further below to investigate the importance of α for the determination of the head losses in unlined rock-blasted tunnels.
The first step of the numerical investigations consisted in a grid refinement study to identify the required grid size by quantifying the numerical uncertainty. Using a discharge of 0.0353 m 3 /s (model scale), cell sizes from 100 mm to 5 mm were tested. The best result was obtained for the grid size of 5 mm and simulations with smaller grid sizes became unstable. It was therefore decided to use the grid size of 5 mm for the further simulations. Because of the high resolution of the computational grid (approx. 96 million cells), a small roughness value k s had to be used for the simulations. A value of k s = 0.5 mm was identified during the calibration of the model, for which the difference in head loss between the physical and numerical model was less than 0.12 mm (corresponding to a difference of 3.1%). This difference was considered to be accurate enough given the achievable accuracy of the scale model tests.
The calibrated model was used to evaluate the head loss along the tunnel for four more discharges after the grid refinement study (Table 3). Figure 8 shows the comparison of the CFD-results with the scale-model data and illustrates that it was possible to reproduce the head losses between cross sections 1 and 8 for all discharges with a high degree of accuracy. For the discharges of 0.054, 0.074, and 0.093 m 3 /s, the differences between the CFD and scale model head losses were less than 2%. On the other hand, for the largest discharge of (Q = 0.103 m 3 /s) we observed a slightly larger deviation (7.3%). Note that the k s + -values ranged between 9 and 25 for the used discharges and that these values were lower than the corresponding y + -values for the smallest cell-size at the tunnel wall (23 ≤ y + ≤ 62).

Accuracy of the Milled Model
The quality of the DEM generated from the point cloud determines the accuracy of the numerical and physical model tests. The fact that 99.5% of the cloud points were less than 0.02 m (in prototype dimensions) away from the generated mesh is a strong indicator for the accuracy of the DEM. A full evaluation of the accuracy of the milled model was not possible because of technical problems which prevented the generation of a point cloud for the right-hand tunnel segments. The comparison of the point cloud of the left-hand tunnel segments with the DEM (94.6% of the scanned points were less than 3 mm away from the scaled DEM) showed larger discrepancies than the comparison of the DEM with the original point cloud. These differences can be associated with the scanning and scale-model production processes. The latter are more important in this study, as the size of the finest milling head is a main parameter for the smoothing of the surface during model production. Milling with a finer milling head increases the production time (and cost) considerably, so that the size of the finest milling head is usually restricted to a diameter of 10 mm for large model sizes [pers. comm. SINTEF

Accuracy of the Milled Model
The quality of the DEM generated from the point cloud determines the accuracy of the numerical and physical model tests. The fact that 99.5% of the cloud points were less than 0.02 m (in prototype dimensions) away from the generated mesh is a strong indicator for the accuracy of the DEM. A full evaluation of the accuracy of the milled model was not possible because of technical problems which prevented the generation of a point cloud for the right-hand tunnel segments. The comparison of the point cloud of the left-hand tunnel segments with the DEM (94.6% of the scanned points were less than 3 mm away from the scaled DEM) showed larger discrepancies than the comparison of the DEM with the original point cloud. These differences can be associated with the scanning and scale-model production processes. The latter are more important in this study, as the size of the finest milling head is a main parameter for the smoothing of the surface during model production. Milling with a finer milling head increases the production time (and cost) considerably, so that the size of the finest milling head is usually restricted to a diameter of 10 mm for large model sizes [pers. comm. SINTEF Ocean]. The model scale of 1:15 results in difficulties in reproducing roughness elements in the order of 15 cm and less in prototype scale. For the present model, larger deviations were observed in roughness cavities reflecting the limitations of the used milling head to reproduce small cavities accurately [27,28]. On the other hand, 71% of the model cloud points were still in the range ±1.5 mm from the scaled DEM, and this is a strong indicator that most of the surface features have been reproduced well. Since the main inaccuracy associated with the milling is related to deeper cavities where the flow velocity is small, it can be hypothesized that the impact of these inaccuracies is negligible for the outcomes of the hydraulic tests.

Geometrical Values Derived from DEM
The spatially heterogenous topography of the tunnels determines the roughness of the tunnels, and the associated cross-sectional area variations, defining the cross-sectional roughness, become visible from the longitudinal profile of the cross-sectional areas shown in Figure 5a. The power spectrum of the longitudinal area-variation in Figure 5b revealed a small peak at a wavelength of about 0.23 m. For wavelengths smaller than 0.23 m, the spectrum may be characterized by a power law (slope-3.2). For wavelengths larger than 0.23 m, the spectrum reaches a plateau value which may be expected as the largest cross-sectional areas will most likely be associated with the blasting frequency. The detailed investigation on how a tunnel is shaped by blasting is beyond the scope of the present paper, as this also depends on geotechnical and technical (explosives etc.) boundary conditions, see e.g., [36,37] for more details.
The distribution of the cross-sectional areas shown in Figure 5c differ for the reaches 1-8 and 2-7. The reason for this difference becomes apparent from Figure 5a showing that the smallest cross-sectional areas are observed at the beginning of the tunnel (before cross section 2 at x = 1.5 m, Table 1). These smaller areas are not considered in the analysis of the area-distribution for reach 2-7 so that the corresponding A 1 value is larger than for reach 1-8. On the other hand, both the mean area A and A 99 are similar for the two reach lengths ( Table 2). The different A 1 -values are the reason for the difference in the δ-value (Table 2), which in turn is required for the determination of the friction factor by Equations (4)- (6). The determined friction factors for the reaches 1-8 and 2-7 are summarized in Table 4 together with the friction factors according to the IBA-and Queen's method. The table shows that the friction factors determined by Equations (4)-(6) are in the same order of magnitude for the same reach but depend on the reach length due to the aforementioned reasons. This in turn illustrates the dependency of the approaches on spatial scales for a given reach, which is characterized by large-scale variations in cross-sectional area. We further note that care needs to be taken when applying Equation (6), as it should only be used for data in prototype dimensions due to the addition of 9 m 2 to A 1 +9 in the denominator and that we used a higher spatial resolution for the application of the approaches than recommended by the corresponding authors. Table 4. Friction factors f according to different approaches.

Reach f [-]
Equation ( Compared to the approaches quantifying only the cross-sectional roughness, the methods based on the additional analysis of longitudinal profiles (IBA and Queen's) result in larger friction factors. Moreover, the friction factors determined by the Queen's method are larger than the ones for the IBA-method, but the differences for the two different reach lengths are small. The reason for the increased friction factors may be associated by the consideration of the surface roughness which is derived from the analysis of longitudinal profiles. This issue is further discussed below with regard to the measured friction factors. The statistical values standard deviation s, skewness Sk, and Kurtosis Ku for the reaches 1-8 and 2-7 from the profile analysis are characterized by only minor differences (Table 2), and the Sk-and Ku-values indicate that the distribution of the elevations along the longitudinal profiles area is almost bell-shaped and similar to a normal distribution.

Experimental Results of Pressure, Head Loss, and Friction Factor
The results from the pressure measurements revealed the effect of the irregular area distribution on the pressure head ( Figure 6). Cross sections with a larger area were characterized by a smaller velocity head and a larger pressure head than cross sections with a smaller area according to the Bernoulli-principle. Despite the head losses along the tunnel, we found increasing pressure heads between the cross sections 2 to 4 and 6 to 7, indicating the effect of the cross-sectional variations on the flow. The increase in pressure head was compensated by a reduction in the velocity head between the cross sections, as the cross-sectional areas increased. However, we still calculated small negative total head losses between cross sections 3 and 4 for some discharges, which cannot be physically based. These negative head losses were determined under the assumption that the kinematic correction factor α is equal to unity. Increasing this factor would mean that the velocity head loss between the cross sections would slightly increase so that positive head loss values may be obtained. However, we note that the determined negative head losses were small (<0.2 mm) and are thus within the margin of error (or uncertainty) for the measurements. In order to minimize the influence of measurement errors and uncertainties, we decided to use a linear regression to determine the reach averaged energy slope S E . This slope was then used in Equation (13) instead of h f,i−j /(x i − x j ) to determine the friction factors for the different discharges (Table 3). Two reach lengths (1-8: 7 m and 8 data points; 2-7: 5 m and 6 data points) were used for the regression analysis to account for the observed behavior of the total head at cross section 4 ( Figure 7). The results presented in Table 3 indicate only minor differences between the friction factors for the two reach lengths. We acknowledge that a certain conduit-length is typically required for the development of the flow, but the heterogeneous roughness of the tunnel hampers the assessment of the required length. Moreover, because of the close similarity of the obtained friction factors for the two reach lengths, it can be hypothesized that the smooth transition from the circular-pipe cross section to the horseshoe tunnel-shape, together with the flow conditioning, helped in reducing the required length for flow development. It is worth mentioning that a more pronounced effect would become visible if only the head loss data of cross sections 1 and 8 as well as 2 and 7 were used for the analysis. Using these data, which is not shown here, would result in a smaller friction factor for the shorter reach which would also be more prone to measurement errors as only two local values were used. Overall, the analyzed reach length ranged between 20 and 15 times the hydraulic diameter (reach 1-8 and 2-7, respectively) which we consider to be a fair length so that potential errors due to the in-and outlet sections are expected to be small.
The determined friction factors show a dependency on the Reynolds number for Re < 1.44 × 10 5 as visualized in Figure 9. For Re > 1.44 × 10 5 , on the other hand, the friction factor becomes independent of the Reynolds number and corresponds in average to f = 0.048 and f = 0.047 for the reaches 1-8 and 2-7, respectively (these values correspond to Manning's n values of 0.026 s/m 1/3 ). The order of magnitude is comparable to friction factors of prototype tunnels reported by [13,15,22] which ranged between 0.045 ≤ f ≤ 0.116. The fact that the present estimates are at the lower end of this range is not surprising because of the straight model reach. The experimentally determined friction factor of f = 0.047 for Re > 1.44 × 10 5 is almost identical to the friction factors estimated via the Equations (4)-(6) for the reach 2-7 (Table 4). For the reach 1-8, the approaches result in larger estimates which, as already mentioned, may be associated with the inhomogeneity of the cross-sectional area at the model entrance area (from 0 ≤ × ≤ 1 m). Both the IBA and Queen's method overestimate the determined friction factors indicating that the combination of cross-sectional roughness and surface roughness based on analyses of high-resolution data may lead to a biased estimate of the friction factor, as also observed by [19]. We hypothesize that the variation in cross-sectional area is correlated to the standard deviation of longitudinal profiles and that this correlation is the reason for the overestimation of the friction factor by the IBA method. The reason for the observed deviation of the Queen's method, however, remains unclear as this approach was also derived from the analysis of TLS data. We will investigate this issue further using data from further tunnels that were acquired in the framework of the present research project.

Numerical Simulations
The objective of the numerical modelling study was to replicate the head loss in the scale model between cross sections 1 and 8. The numerical model was calibrated for the lowest discharge (0.035 m 3 /s) by adjusting the roughness height ks and the results presented in Figure 8 showed that the model performed well for the other discharges without re-calibrating ks. Slightly larger deviations were observed for the largest discharge used, but the corresponding differences were still below 10%.
A closer inspection of the numerical data revealed that the total heads between cross sections 2 and 7 deviated slightly from the measured ones, especially at the cross sections 2, 4, and 6. This is shown in Figure 10, for which the total head of the numerical model at cross section 1 was set equal to the corresponding total head from the scale model.  The reason for the increase of the friction factor with decreasing Re cannot be fully resolved in the present study but may be associated with the influence of flow conditions similar to the flow in pipes (transition from smooth to rough conditions in the Moody diagram). The corresponding line delineating these two regimes in the Moody diagram is also included in Figure 9, although the cross-sectional shape of the tunnel is not circular. Assuming the validity of the line, the determined friction factors are in the fully turbulent hydraulically rough regime. We note that a dependency of the friction factor on Re has also been reported by Priha [13], although he identified mostly a decreasing friction factor with decreasing Reynolds number. Our results, on the other hand, show a similar trend as the lines in the Moody diagram for given relative roughness, see [8].
The experimentally determined friction factor of f = 0.047 for Re > 1.44 × 10 5 is almost identical to the friction factors estimated via the Equations (4)-(6) for the reach 2-7 (Table 4). For the reach 1-8, the approaches result in larger estimates which, as already mentioned, may be associated with the inhomogeneity of the cross-sectional area at the model entrance area (from 0 ≤ × ≤ 1 m). Both the IBA and Queen's method overestimate the determined friction factors indicating that the combination of cross-sectional roughness and surface roughness based on analyses of high-resolution data may lead to a biased estimate of the friction factor, as also observed by [19]. We hypothesize that the variation in cross-sectional area is correlated to the standard deviation of longitudinal profiles and that this correlation is the reason for the overestimation of the friction factor by the IBA method. The reason for the observed deviation of the Queen's method, however, remains unclear as this approach was also derived from the analysis of TLS data. We will investigate this issue further using data from further tunnels that were acquired in the framework of the present research project.

Numerical Simulations
The objective of the numerical modelling study was to replicate the head loss in the scale model between cross sections 1 and 8. The numerical model was calibrated for the lowest discharge (0.035 m 3 /s) by adjusting the roughness height k s and the results presented in Figure 8 showed that the model performed well for the other discharges without re-calibrating k s . Slightly larger deviations were observed for the largest discharge used, but the corresponding differences were still below 10%.
A closer inspection of the numerical data revealed that the total heads between cross sections 2 and 7 deviated slightly from the measured ones, especially at the cross sections 2, 4, and 6. This is shown in Figure 10, for which the total head of the numerical model at cross section 1 was set equal to the corresponding total head from the scale model. The head values for the subsequent cross sections were then calculated from the computed head losses along the tunnel. Compared to the scale model data, the total head values from the numerical simulations followed the measured gradient without larger deviations from the fitted straight line.  The numerical model determines friction losses by computing form drag and skin friction. The form drag exerted by the larger roughness elements, which are resolved by the grid, is computed directly with the help of the turbulence model. If the simulations were carried out with an extremely fine grid, the finest cells close to the wall would be in the laminar sub-layer and the roughness height would hence not be needed, as smooth wall laws could be used. However, such fine grids were not feasible because of the limited computing resources. The chosen basic grid size of 5 mm was determined in the grid refinement study for the calibration discharge, and gave good results for the head loss between cross sections 1 and 8 (see Figure 8 and Table 3). Finer basic grid sizes could not be used as the corresponding simulations became prohibitively expensive. Therefore, very fine geometrical roughness features need to be included in the numerical model in terms of skin friction which was computed via the wall function. In the calibration process, the roughness height of ks = 0.5 mm gave the best results.
We note that the calibration procedure is also affected by the ability of the model to compute the correct drag on the large roughness elements. If the grid is too coarse to resolve recirculation zones behind larger roughness elements, the form-induced drag cannot to be computed correctly. Moreover, the k-ε turbulence model will not give an exact estimate of the drag on the obstacles. Thus, in the current approach where the roughness height ks was calibrated, the procedure will not only account for the actual skin friction but also compensate for the uncertainties in drag estimates. This means that the calibration procedure can give roughness heights which are not physically based. However, the calibrated roughness height of ks = 0.5 mm is small and 10% of the basic cell size (or The numerical model determines friction losses by computing form drag and skin friction. The form drag exerted by the larger roughness elements, which are resolved by the grid, is computed directly with the help of the turbulence model. If the simulations were carried out with an extremely fine grid, the finest cells close to the wall would be in the laminar sub-layer and the roughness height would hence not be needed, as smooth wall laws could be used. However, such fine grids were not feasible because of the limited computing resources. The chosen basic grid size of 5 mm was determined in the grid refinement study for the calibration discharge, and gave good results for the head loss between cross sections 1 and 8 (see Figure 8 and Table 3). Finer basic grid sizes could not be used as the corresponding simulations became prohibitively expensive. Therefore, very fine geometrical roughness features need to be included in the numerical model in terms of skin friction which was computed via the wall function. In the calibration process, the roughness height of k s = 0.5 mm gave the best results.
We note that the calibration procedure is also affected by the ability of the model to compute the correct drag on the large roughness elements. If the grid is too coarse to resolve recirculation zones behind larger roughness elements, the form-induced drag cannot to be computed correctly.
Moreover, the k-ε turbulence model will not give an exact estimate of the drag on the obstacles. Thus, in the current approach where the roughness height k s was calibrated, the procedure will not only account for the actual skin friction but also compensate for the uncertainties in drag estimates. This means that the calibration procedure can give roughness heights which are not physically based. However, the calibrated roughness height of k s = 0.5 mm is small and 10% of the basic cell size (or 40% of the wall-near cell size), which is not unreasonable from a physical point of view, as most of the roughness was directly resolved by the grid, see also [38]. Note also that the k s + -values were always larger than the y + -values (computed using the smallest cell-size at the wall). This indicates that the computed drag on the roughness elements is reasonable, and the obtained results for the bulk head losses between cross sections 1-8 support this statement. We note that the used wall functions have many empirical coefficients that may be tuned in the calibration process. However, a detailed investigation of this issue was beyond the scope of the study and these empirical coefficients were not tuned in the calibration process. We further note that the study carried out by [30] revealed the applicability of the k-ε model for modelling the flow in an unlined tunnel. Localized pressure variations due to the variation in hydraulic diameter and cross-sectional areas could be predicted. On the other hand, Andersson et al. [30] commented critically on the performance of the model to capture the head loss at larger cross sections, which is in line with the results presented above. Finally, another possible reason for the observed deviations between the numerical simulations and the scale model may be associated with uncertainties from the scale model. However, the pattern of the head values along the tunnel was similar for all discharges (and identical for measurements carried out for the same discharge at different days), and this indicates that the measurements were reproducible and reliable.
The results from the numerical study were also used to cross-check the experimental setup for the pressure measurements. For this purpose, pressure values were extracted at the locations of the pressure taps at each cross section. These values were then averaged and compared to the modelled cross-sectionally averaged pressure. This comparison revealed almost identical results which is a strong indicator of the reliability of the pressure measurements in the scale model. We also note that the results from the CFD study follow the shape of the f-Re relationship from the experiments indicating the existence of a transitionally rough regime and hence confirming the experimental results.
The aforementioned differences in the total head between the numerical and experimental investigations resulted in slightly different estimates of the friction factors based on the slope of the energy line along the tunnel for the reach between cross sections 1 and 8. For the numerical simulations carried out with a kinematic correction factor of α = 1, the friction factors varied between 0.048 and 0.054 for the different discharges (see Figure 9) and similar friction factors were obtained for the reach between cross sections 2 and 7. It is interesting to note that smaller friction factors were obtained when the kinematic correction factor was considered in the calculations. This is indicated in Figure 9, which includes the friction factors that were determined using the numerical simulations to determine α (a reach averaged α value was used by averaging the eight available α i values at the cross sections). The corresponding differences ranged between 13% and 15%. On the other hand, the differences in head losses between cross sections 2 and 7 were only around 4% (all discharges) and around 1% for the discharges up to 0.093 m 3 /s and 9% for Q = 0.1032 m 3 /s between cross sections 1 and 8. This result indicates that the consideration of the kinematic correction factor may affect the slope of the energy line which in turn affects the value of the friction factor.
Despite the critical discussion of the results, it needs to be stated that the differences in friction factors between the numerical simulations with α = 1 (f = 0.049 in average for Re > 10 5 ) and the scale model (f = 0.048 for cross sections 1-8) is small for larger Reynolds numbers. The difference becomes larger for lower Reynolds numbers (and if α would be considered), but may still be interpreted as small.

Conclusions and Outlook
A combined scale model and numerical study was carried out to evaluate the friction losses in a pressure-driven flow through an unlined rock-blasted hydropower tunnel. Remote sensing data from the prototype tunnel captured by terrestrial laser scanning provided the input data for both the construction of a 1:15 scale model with an innovative milling approach and a numerical model. The accuracy of the milled model was high so that the large-and medium-scale roughness features of the prototype model were adequately reproduced in the 1:15 scale model. The high-resolution spatial data of the tunnel were used for the analysis of the performance of existing approaches to determine the friction factor in unlined tunnels. For this purpose, pressure measurements were carried out for different discharges at eight cross sections so that the total head and hence the friction factor could be evaluated. This was done for two different reach lengths of the tunnel which accounted for the variation in cross-sectional geometry, but we found that the corresponding friction factors for the two reach lengths were similar. We also found a dependency of the friction factor on the Reynolds number for lower discharges which was no longer visible for larger discharges (i.e., larger Reynolds numbers), which is in line with the results reported in the literature. The friction factor determined from the physical model study for the 5-m long reach corresponded to f = 0.047 which is in the range of friction factors reported in the literature for unlined rock-blasted tunnels.
Existing approaches for the determination of the friction factor based on the cross-sectionally roughness performed rather well, although the considered reach length was limited to 8 m (120 m in prototype scale), and it was pointed out that local reaches with considerably deviating cross-sectional areas may bias the friction factor estimates. On the other hand, we found that more sophisticated approaches accounting for both surface and area irregularities overestimated the friction factor. We therefore hypothesize that the surface and cross-sectionally roughness are to some extent correlated.
The numerical simulations were based on the solution of the Reynolds averaged Navier-Stokes Equations using the CFD program OpenFoam and the k-ε model. The numerical model was calibrated for the head loss between the first and last measurement cross section in the scale model, and the calibrated model performed very well for other discharges. However, slight differences between the experimental and numerical head losses were observed at other cross sections, and this was critically discussed regarding the performance of the k-ε model, wall-laws, and grid size. The performance of the k-ε model will be evaluated further by using velocity data that was acquired in the tunnel by means of PIV. These data were not presented in the paper as its focus was on the evaluation of the friction losses. One of the main objectives of the numerical study was to investigate the applicability of methods that can also be used by practitioners. Our results confirm the conclusions of [30] that it is possible to calibrate a computationally efficient model for a specific target value (here: head loss between two specific cross sections). We will further investigate the applicability of the numerical model with data from three other tunnels, which were acquired in the framework of the present research project. These investigations will also include further detailed analyses of the roughness geometry to link frictional losses to physical characteristics of the wall topography. Overall, the results of the study show the usefulness of the chosen hybrid approach and that existing approaches for the determination of head losses in unlined tunnels need to be further refined.