Evaluation on the Effect of Pressure Transients on Rock Joints in Unlined Hydropower Tunnels Using Numerical Simulation

Frequent pressure transients are identified as the cause of block failures in many unlined hydropower tunnels. The primary design objective of such tunnels is to prevent hydraulic jacking at design static pressure and mass oscillation but neglects the effect of short transients, i.e., water hammer. The issue has not been studied from the perspective of hydro-mechanical interactions due to frequent pore pressure changes in the rock mass. This article mainly focuses on the effect of pressure transients at different static heads, or different effective normal stresses across the joints and the effect of time period of pressure transient. Further, the change in such behaviour due to different mechanical properties of rock joints, such as stiffness, friction angle and dilation, is investigated. Numerical simulations of observed pore pressure response in the rock mass during a pressure transient are carried out using distinct element code 3DEC. The results show that relative joint deformation due to short pressure transients are the highest when joint normal stresses are 1.5–2.5 times higher than static water pressure in the tunnel and thus the vulnerability to weakening of such joints by hydraulic fatigue is higher. Further, results show that water hammers can travel up to 4 m into the rock mass even in stiff joint conditions and sufficiently high normal stresses. Results further indicate that the hydraulic impact due to water hammer is smaller as compared to mass oscillation. It is concluded that water hammers, wherever applicable along the waterway, can still contribute to hydraulic fatigue of rock joints in addition to the effect of mass oscillation and cannot be neglected when pressure transients occur frequently. Tunnel filling/dewatering and mass oscillations cause macroscopic joint displacements or block movements over long-term operation which is the major cause of block falls in unlined pressure tunnels.


Introduction
Hydromechanical coupling in fractured media is a major field of scientific research in rock engineering and has applications in engineering projects that involve fluid flow in rock mass, such as reservoir engineering, nuclear waste disposal, inflow during underground excavations, design of dam foundations, etc. A previously unexplored application is in the field of design of unlined pressure tunnels functioning as waterways for hydropower projects, where there is direct contact between flowing water and rock joints exposed to tunnel wall. Any change in water pressure in tunnel causes the rock mass pore pressure to change and thus involves the process of hydro-mechanical coupling. It should be noted that the term "Pore pressure" in this manuscript solely refers to water pressure in the rock joints and the pressure in the rock stratum is not considered.
The state-of-the-art design principle of unlined pressure tunnels (Palmstrøm and Broch 2017;Panthi and Basnet 2018) considers static tunnel pressure as the primary design load and the major design objective is to prevent hydraulic jacking. Most of the tunnel length is placed in the rock mass where in situ minor principle stress is higher than the static tunnel water pressure. Steel or concrete lining is carried out in places where the minor principal stress is smaller. This is a well-tested design criterion with applications all over the world (Panthi 2014;Broch 2010;Marwa 2004) and especially in Norway, where more than 95% of hydropower pressure tunnels are unlined, with a history of over 100 years.
However, the operation regime of hydropower power plants in Norway has changed drastically over the years, 1 3 especially since 1990s. The share of intermittent energy sources, such as solar and wind power, in the energy market is increasing every year (Charmasson et al. 2018). In EU member states, the share of renewable energy has almost doubled between 2004 and 2018 and wind and solar power share a large portion of this energy mix (Eurostat 2020). Hydropower plants must act as battery to supply energy when these intermittent sources are not producing. Thus, the operation of hydropower plants is becoming more dynamic, with load changes occurring on daily and even hourly basis. This induces large number of pressure transients, and frequent changes in rock mass pore pressure around the tunnel. It has been documented by Bråtveit et al. (2016) that block falls in tunnels have increased by a factor of 3-4 times as a result of dynamic operation. Such block falls are attributed to the cyclic fatigue of rock joints caused due to additional seepage forces acting on joint walls during pressure transients. Hence, in addition to hydraulic jacking, block falls during long-term operation of power plants are also becoming an important design issue which needs to be addressed.
The studies conducted so far regarding unlined pressure tunnel design and their failures (Brox 2019;Basnet and Panthi 2018;Rancourt 2010;Helwig 1987;Lang et al. 1976) point out the need for a detailed investigation in this topic. In Norway, there is now a renewed interest in understanding the consequences of hydromechanical effects regarding potential block falls in unlined pressure tunnels during long-term operation. The main reason being, majority of the 4300 km long waterway system in Norwegian hydropower plants consisting of unlined tunnels subjected to high static pressures maximum up to 1047 m.
A full-scale field instrumentation and monitoring is conducted in a headrace tunnel of a Norwegian hydropower plant to measure changes in rock mass pore pressure near tunnel walls during plant operation (Neupane et al. 2020). In this article, a brief summary of the instrumentation and some results relevant to the numerical simulation are presented. Numerical simulation of the observed pore pressure response in the rock mass during a pressure transient event is carried out using 3D distinct element code 3DEC. The model is then used to evaluate the impact at various static water heads and time periods of pressure oscillation in the tunnel, which are the key parameters for unlined pressure tunnel design. A parametric study has also been conducted to observe how the behaviour changes with some relevant rock mechanical parameters. Finally, implication of the observed results is linked to the actual block fall events.

Physical Behaviour During Tunnel Filling and Pressure Transients
An unlined water tunnel is filled at a slow pace so that there are no sudden changes in pore pressure and effective stresses in the rock mass. When a tunnel is being filled, the effective stresses are reduced, which causes the rock joints to deform due to normal opening and shear dilation. When it is ensured that the minor principal stress is higher than the maximum tunnel water pressure, the displacements are negligible, and joint apertures do not change significantly, thus limiting the leakage to an acceptable range. When the tunnel pressure is steady (during steady plant operation), the displacements remain constant. When turbine valves are closed to shut down or change the production, a pressure transient is induced in the tunnel. Pressure transient in a typical hydropower waterway consists of water hammer and mass oscillation. Water hammer is a pressure oscillation that occur in the conduit between the turbine and the surge shaft and have a short time period similar to sound waves in air or water. The water hammer is superimposed with the other, much slower but long-lasting oscillation (up to hours), i.e., mass oscillation, which gradually builds up when the water hammer has died down (Jaeger 1977). The time period of water hammer and mass oscillation for Roskrepp waterway are about 3.5 and 200 s, respectively. Mass oscillation can continue for hours after a load change event in the power plant, whereas water hammer dies down relatively faster.
During transient upsurge (increasing pressure), higher pressure is acting on the joint wall surfaces, further reducing the effective normal stress. A factor of safety is provided for static water head and for mass oscillation. Hence, the effective stress is still high enough to limit the joint deformations. On the other hand, during transient downsurge (decreasing pressure), higher pore pressure will be ¨trapped¨ in the rock joints due to pressure transmission delay between the tunnel and joint surfaces. This tends to push the rock blocks out of the tunnel wall. This phenomenon with higher rock mass pore pressure has been observed in the field experiment and is explained in Sect. 3 (also referenced to Neupane et al. 2020).
It is considered that pressure transients with a short time period, i.e. water hammer, does not require a factor of safety since the time of application is too short to cause hydraulic jacking (Benson 1989). Still, water hammer has a higher potential for pressure increase as compared to the mass oscillation, even though its zone of influence is lesser as compared to mass oscillation, due to smaller time period. Thus, it is believed that water hammer may also contribute to block falls, even though the magnitude is small as compared to mass oscillation. During first few years of operation, the joint stiffness in areas closer to the tunnel wall is usually reduced because the infilling material is washed away as a result of pore pressure changes. Yet, block falls in tunnels which are under operation for over 30 years, indicate that the joint properties continue to change as a result of prolonged cyclic loading and fatigue.

Theoretical Aspect
The mechanical and hydraulic processes that govern fluid flow in rock joints are interdependent with each other, are called hydromechanical couplings and are divided into two types, i.e., direct and indirect couplings. As described by Wang (2000), direct coupling includes two basic phenomena: I. Solid-to-fluid coupling It occurs when a change in applied stress produces a change in fluid pressure or fluid mass. The applied stress produces displacement in the rock joints. This deformation generates surface stress on the fluid domain boundary, which deforms accordingly. A reduction in channel volume induces fluid outflow. II. Fluid-to-solid coupling It occurs when a change in fluid pressure or fluid mass produces a change in volume of the porous medium. A fluid inflow induces fluid pressure along the flow channels, which act on the channel boundaries and deforms the surrounding rock material. As a result of deformation, the rock counteracts the fluid pressure with surface stress at the fluid-rock boundary, which affects the fluid pressure and volume of fluid domain.
Mechanical and hydraulic processes can also affect each other through change in material properties, which are considered as indirect coupling. For example, the reduction in channel volume may increase contact area between the joint surfaces resulting a stiffer material. Indirect HM couplings tend to be most important in fractured rock mass or intact rock with flat inter-grain micropores, where changes in permeability caused by fracture or pore dilation can be dramatic (Rutqvist and Stephansson 2003). Indirect coupling composes of two basic phenomena: a solid-to-fluid coupling that occurs when an applied stress produces change in hydraulic properties; and a fluid-to-solid coupling that occurs when a change in fluid pressure produces a change in mechanical properties.

Instrumentation Site
The field instrumentation is carried out in the unlined headrace tunnel of Roskrepp power plant in southern Norway, which has been in operation since 1980. The instrumentation setup is presented in Fig. 1a; readers are referred to Neupane et al. (2020) for details. The rock mass at this location consists of weakly schistose granitic gneiss with three major joint sets. The fracture network at this location consists of two sets of low permeability joints, i.e., foliation joints and joint set J1. A relatively flat dipping joint set J2 is exposed in the tunnel crown and has high permeability. One single joint (Jf conductive ) of the same orientation as the foliation joints also has relatively larger permeability (Fig. 1b). Among these joint sets, Jf and J1 have a spacing 1-2 m whereas Jf conductive has spacing of about 10 m and has persistence of more than 10 m. The intact rock is of crystalline nature and thus can be considered impermeable such that the pore pressure in the rock mass is only governed by fluid flow through rock joints.

Instrumentation Setup
Five boreholes are drilled in the tunnel walls such that they intersect a particular joint set almost perpendicularly (Fig. 1b). Stainless steel pipes are fixed in the boreholes using packers at different lengths in the borehole and the pipes are taken out of the tunnel to a dry area in the access tunnel, where they are fitted with pressure sensors and data logger (Fig. 1a). The length of boreholes ahead of the packers (marked in blue in Fig. 1b) collects water from the rock joints and is connected to the pressure sensors through the steel pipes. Simultaneous readings of the tunnel water pressure are also taken from a pipe installed at the junction between the headrace tunnel and access tunnel.

Test Results
The pressure readings from boreholes BH1, BH2 and BH4 during a typical pressure transient event are shown in Fig. 2 along with the water pressure in the tunnel. It is seen that the boreholes which intersect the conductive joints, i.e., BH1 and BH4 strongly respond to pressure transients whereas BH2 is non-responsive, since it does not have direct hydraulic contact with the conductive joints. Even though BH3 intersects Jf conductive , it does not have a direct hydraulic connection because the intersection point is behind the packer. (Fig. 1b).
The delayed pore pressure response or timelag between pressure peaks in the two responsive boreholes is of relevance regarding potential effect on rock block stability. The shaded area in Fig. 2b shows the condition when the rock mass pore pressure is higher than the tunnel pore pressure during a negative pressure wave. It is observed that such situation occurs for the first few cycles of the pressure transient and then gradually decreases as the transient attenuates and both pressures become almost equal or the tunnel pressure starts to become higher (Neupane et al. 2020). Therefore, the first few cycles are of interest regarding block displacement due to pressure transient and will be simulated in the numerical models.
It is observed that BH1 registers a stronger response during water hammer ( Fig. 2b) but there is very little timelag during mass oscillation. For BH4, clear timelag is registered during both mass oscillation (Fig. 2a) and water hammer (Fig. 2b). But the amplitude of pore pressure in BH4 in response to the water hammer is smaller as compared to BH1. The difference in response is due to different resistance to flow through joints, which is a function of void geometry of joints and the length of flow path, i.e., joint length between tunnel wall and its intersection point with individual borehole.

Intact Rock Properties
Laboratory tests were carried out on rock samples of granitic gneiss collected during field mapping showed mean values of uniaxial compressive strength and modulus of elasticity of 148 MPa and 65 GPa, respectively. A summary of lab test results on intact rock samples and other estimated properties are presented in Table 1. The intact rock is modelled as an isotropic, homogeneous, linearly elastic material.
Normal and shear stiffness of joints are among the key parameters needed for numerical simulation. However, there is a great deal of uncertainty regarding these parameters since these are not easily measurable in the field and very limited amount of data are available. Zangerl et al. (2008) concludes that laboratory and in situ normal closure experiments in granitic rock showed a very Fig. 2 a Pore pressure response of BH1, BH2 and BH4 during a typical pressure transient event in the headrace tunnel and b response during water hammer peak large range of stiffness characteristic values, even for welldefined laboratory tests within the same rock type. This is because, fracture normal stiffness is highly affected by several extremely complex interacting factors, such as fracture surface geometry, asperity deformability, fracture interlocking, testing condition, etc., which are difficult to determine, if not entirely impossible. In addition, joint properties measured in the laboratory typically do not represent field joint conditions because of scale effects. For this analysis, the joint stiffness is calculated assuming that the jointed rock mass has same deformational response as an equivalent elastic continuum, for uniaxial loading of rock containing a single set of uniformly spaced joints oriented normal to the direction of loading (Barton 1972). To account for the uncertainty, a range of rock mass deformation modulus has been estimated using three different methods, i.e., Barton (2002), Hoek and Diederichs (2006) and Panthi (2006) and are presented in Table 2. Minimum and maximum values of rock mass deformation modulus have been calculated using each of these methods.

Joint Constitutive Model and Properties
An elasto-perfectly plastic Mohr-Coulomb model is chosen for simulating the fracture behavior. During simulation with this joint model, the fracture stiffness is kept constant, assuming it to be independent of normal stress, and friction and dilation angle are also kept constant. This is a reasonable assumption since the joint deformation during one pressure transient event can be very small. In this model, the fracture normal stiffness controls the normal deformation due to change in normal effective stress and shear stiffness controls the elastic shear behaviour of the fracture. The plastic shear behaviour is governed by fracture shear strength, which is a function of friction angle, cohesion and effective normal stress. The aperture changes as a result of joint displacement/dilation, which governs the flow rate and pore pressure along the joint. The maximum and minimum values of joint stiffness have been calculated based on rock mass modulus (Table 2), and presented in Table 3, along with other parameters used in the numerical simulation.

Joint Hydraulic Properties
Water inflow into the tunnel was noted until a month after tunnel dewatering through the joints Jf conductive and J2, which shows that these joints are highly permeable and act as direct hydraulic connection between natural ground water table and the tunnel. Other joints are tight and do not indicate any seepage through them, indicating low permeability. This is also established by the fact that BH1 (intersecting Jf conductive ) responds quickly to the pressure transients. Further, outflow was measured at the pipe outlet of all boreholes to have an idea of the relative difference in permeability between the joint sets. The de-aeration valve connected at the pipe outlets was opened and flow was measured for about 5 min when the tunnel registered a static water pressure of 9.45 bars. The steady flow rate observed during such event gives the flow occurring through the length of joints between the tunnel and the boreholes at a known tunnel water pressure. The outflow from BH1 was measured to be 0.0012 L/s. The outflow from BH2, BH3 and BH5 was negligible and was difficult to be measured with relative accuracy. This shows that the differences in permeability between the conductive Jf conductive and other joints are in several orders of magnitude.

In Situ Rock Stresses
An estimate of the in-situ rock stress at the instrumentation location is made based on 3D stress measurement carried out at Holen power plant in 1980 (Myrvang 2019), located about 30 km from Roskrepp power plant. Both these locations are in the same geological and geo-tectonic setting consisting of massive Precambrian gneisses. Hence, the tectonic stresses are assumed to be similar in these locations. At Holen, the maximum principal stress is approximately horizontal and oriented towards east-west direction and magnitude vary between 12 and 23 MPa. The minimum principal stress is approximately horizontal also and oriented towards north-south direction and has magnitude of 5 MPa. Intermediate principal stress is approximately vertical and has magnitude of 7.5 MPa in accordance with the gravity induced vertical stress ( v ) . Since major (σ 1 ) and minor (σ 3 ) principal stresses are horizontal, they are equal to maximum ( hmax ) and minimum ( hmin ) horizontal stresses at this location. Hence, the maximum ( tmax ) and minimum ( tmin ) tectonic stresses are calculated using the following equation (Panthi 2012): These tectonic stress values are then used to calculate maximum and minimum total horizontal stresses at the instrumentation location using Eq. (1) and are presented in Table 4.

Numerical Simulation
The numerical simulation is carried out using three-dimensional distinct element code 3DEC (Itasca 2018). 3DEC can calculate fluid flow and effect of fluid pressures on rock/soil, based on specified material properties and fluid/mechanical boundary conditions using coupled hydro-mechanical calculations through a network of fractures between deformable rock blocks. The flow rate through contacts is calculated by assuming interface as two parallel plates with a defined aperture width, where flow is laminar and is governed by the modified form of cubic law (Witherspoon et al. 1980) given by the following expression: where q is the flow rate, b hi is the initial hydraulic aperture at initial effective stress, f is a factor reflecting the influence of roughness on the flow tortuosity, ΔU n is the change in fracture normal displacement, w is the fracture width, is the fluid density, is the fluid dynamic viscosity, Δp is the pressure difference; and l is the length of joint. Fracture deformation is calculated as a function of effective normal stress assuming a constant fracture stiffness. The domain pressure is then updated, considering the net flow into the domain and possible changes in domain volume due to incremental motion of surrounding blocks. The new domain pressure is computed using following equation: where p 0 is the initial pressure, q is the flow rate defined from Eq. (2), k w is the bulk modulus of the fluid, V is the initial fracture volume and ΔV is the change of fracture volume due to deformation. As a result of new domain fluid pressure, the forces exerted by the fluid on the edges of surrounding blocks are obtained. New mechanical calculations are carried out to update the geometry and the process is continued until equilibrium is reached.

Model Geometry and Boundary Conditions
The model geometry consists of a rock mass block of 60 × 60 × 60 m with an inverted D-shaped tunnel section with height and width of 8 m. The rock mass block consists of a single fracture resembling Jf conductive . This model size was found to be sufficient to prevent boundary effects since the radius of influence around tunnel periphery regarding pore pressure changes during tunnel infilling is smaller than the joint area formed by this model size. It is also seen that the radius of influence during a pressure transient is even smaller as compared to the radius affected during tunnel infilling. Figure 3 shows the model block and the joint plane showing pore pressure buildup and joint displacements during tunnel filling. The model consists of a graded mesh with the finest discretization of 1 m zone size around an area equal to one tunnel diameter in all directions from the tunnel boundary and direction perpendicular to the rock joint. Rest of the model is discretized with 2 and 4 m zone size with the coarsest discretization along outer boundaries.
Fixed boundaries are applied on all faces of the rock block. The initial stresses are set as depth-dependent with the minimum in situ stress values from Table 4 acting at the centre of the model. The boundary pore pressure is set in all faces according to the hydrostatic pressure gradient corresponding to water pressure of 78 m (0.76 MPa) at the tunnel centre, with a linearly varying pressure gradient. The model allows outflow from the boundary. The same magnitude of initial pore pressure is applied on the joint with a linearly varying pressure gradient. This value is estimated based on the elevation of water surface in small lakes above the instrumentation location, which are hydraulically connected to the tunnel through conductive joints.

Calculation Sequence
The calculation sequence for each simulation consists of the following steps: Fig. 3 a Model geometry (60 × 60 × 60 m) and b-d fracture plane with pore pressure, normal displacement and discharge rate contours during tunnel filling 1. Solving for initial state to achieve mechanical equilibrium of the rock block. 2. Mechanical equilibrium after tunnel excavation. 3. Flow only calculations to achieve steady state fluid pressures. 4. Application of water pressure in the tunnel and run to steady state (hydromechanical coupling) to simulate tunnel filling. 5. Time-dependent pressure variation in the tunnel boundary (hydromechanical coupling) to simulate pressure transient (Fig. 4).
The comparison between results of various simulations is done using three output parameters, i.e., maximum normal and shear displacements, and the area bounded between the pressure curves of the tunnel and the rock mass, when rock mass pore pressure is higher than tunnel pressure (shaded area in Fig. 2b). This area is the result of transmission delay between two pressure peaks, during the transient event which is the additional seepage force acting on the joint surfaces during pressure transients. This term hereafter named as ¨hydraulic impact¨ has a unit of pa s and is similar to dynamic viscosity or the force acting on the joint surface per unit area over a certain time.

Fluid Bulk Modulus
The bulk modulus of water at 20 °C is 2.2 GPa, which if used in the calculations, results in very small time steps, making the calculations time-consuming. To mitigate this problem, Itasca (2018) suggests that fluid bulk modulus be reduced such that the apparent stiffness of fluid-filled joint (fluid bulk modulus divided by hydraulic aperture) is approximately equal to the equivalent stiffness of the adjacent zone in the block representing the rock. The average length of zone in the model is 1 m, which gives equivalent stiffness of the zone adjacent to the joint to be in the range of 8.4-16.7 GPa. The apparent stiffness of joint, on the other hand, is at least three orders of magnitude higher than these values. Comparative simulations were run using bulk modulus values of 0.02 and 0.2 GPa, which resulted time steps of 9.5e−4 and 9.5e−5 s, respectively. The difference in calculated timelag, maximum pore pressure magnitude, maximum normal and shear displacements between these zone sizes was less than 10% whereas the runtime was almost 10 times higher for 0.2 GPa. Ivars (2006) mentions that there was no noticeable effect on the equilibrium results when conducting simulations with reduced fluid bulk modulus and that it only exhibits transient differences. Hence, fluid bulk modulus of 0.02 GPa was deemed to be sufficient and has been used for all calculations.

Validation of BH1 Pore Pressure
The conductive joint Jf conductive , was modelled for the largest water hammer peak (Fig. 2b) that occurred after complete closure of the valve. The water hammer peak has a cycle time of about 3.5 s, which travels through joints into the borehole, which registers a pore pressure response with a timelag of 0.75 s between the pressure peaks.
In the simulation, the total normal stress acting across the joint before tunnel filling is 6.25 MPa and effective normal stress after tunnel filling is 5.25 MPa. The result shows that exact simulation of both pore pressure magnitude and timelag was not possible with a single set of parameters. Hence, two different values of joint normal stiffness were used to simulate these two parameters separately (Fig. 5). It Fig. 4 Joint pore pressure contours a when water hammer is at peak in the tunnel, b when water hammer is at negative peak was found that low normal stiffness resulted a better simulation of the timelag but the pressure magnitude was overestimated. Higher normal stiffness resulted a better match for pore pressure magnitude but increased the timelag. Cappa et al. (2008) reports similar pore pressure behavior as a result of change in aperture during HM coupling when modelling pressure pulse in rock fractures. This is because low normal stiffness results in larger hydraulic aperture during tunnel filling and thus the response is faster during a pressure transient and vice versa. At steady state during filling, the maximum hydraulic aperture around the tunnel contour is 99 µm for high normal stiffness and 170 µm for low normal stiffness cases, respectively. The joint displacements can also be seen in Fig. 5.
Since timelag and pore pressure magnitudes were separately simulated, the hydraulic impacts for these two simulations are also different. At joint length of 2.3 m, the measured hydraulic impact is 8.1e4 pa s while the simulated hydraulic impact for low and high normal stiffness cases are 8.9e4 and 7.8e4 pa s, respectively. This shows that the simulation with high stiffness is a better simulation in terms of total hydraulic impact per cycle.

Hydraulic Impact Along Joint Length
The effect of pressure pulse at different lengths along the joint was also analyzed with the same model (Fig. 6). During the pressure pulse, pressure magnitude is higher at locations near the tunnel wall with the smallest timelag. As we go deeper, the pressure magnitude decreases as a result of frictional loss due to increased resistance to flow. The timelag increases as the pressure pulse travels deeper into the rock mass. Since the hydraulic impact is a combination of these two parameters, it is not the highest closest to the tunnel wall since the timelag is too small. On the other hand, deeper into the rock mass, it is again small because the pressure Fig. 5 Simulation results for BH1 during pressure pulse for a high normal/shear stiffness and b low normal/shear stiffness Fig. 6 a Effect of pressure pulse at different lengths of Jf conductive, and b hydraulic impact and maximum joint displacements during the pulse along joint length magnitude is too small even though the timelag is larger. In this case, the highest hydraulic impact is between 2 to 4 m from the tunnel wall, along the joint length.
It is also seen that, the maximum hydraulic impact for both stiffness values is almost equal but the curve for high stiffness joint is much steeper and the zone of influence is much smaller as compared to the less stiff joint. This means that the pressure pulse will impact a larger area of the rock mass when the joint has less stiffness. This seems logical since for a smaller stiffness, larger hydraulic aperture is formed at the end of tunnel filling. Hence, the flow is larger during a transient, and the pressure can travel deeper into the rock mass and thus the total impact will be distributed over a deeper area as compared to the stiffer joint.
A detailed analysis of such impact is carried out in the following sections. The area of interest with respect to block stability is the joint length of less than one tunnel diameter from the tunnel wall, which is the most vulnerable area in terms of block stability and also the starting point for a larger cave-in as a result of consecutive block failures over a period of power plant operation. To gain a better understanding of the behaviour under varying conditions of power plant operation and joint properties, several similar simulations are carried out by varying the parameters which are of specific interest, particularly the tunnel static head and time period of pressure wave oscillation.

Effect of Varying Static Tunnel Water Pressure
More than 200 unlined pressure shafts and tunnels with a combined length of over 4300 km are currently in operation in Norwegian hydropower plants with a maximum static head of up to 1047 m. Palmstrøm and Broch (2017) presents a graph which shows that more than 95% of these tunnels/ shafts in Norway have static heads below 600 m. In this analysis, simulations are carried out to see the effect of increasing static head up to 600 m. The joint properties are the same as for high stiffness case for the model used for validation in Sect. 6.1, and the normal and shear joint stiffness values of 30 and 11 GPa/m, respectively, are referred to as stiffness ST1 hereafter. The time period of pressure pulse is kept constant at 3.5 s and the peak magnitude is increased linearly, relative to the static head increase (14.5% of the static head) in all simulations. The normal stress acting across the joint before tunnel filling is 6.25 MPa and effective normal stress reduces during tunnel filling according to the static head for each case. In terms of conventional design approach, the Factor of Safety (FoS) of design is the ratio between the minimum principal in situ stress and the static water pressure in the tunnel. Hence, the FoS gradually decreases as the design static water pressure increases.
It is important to study the changes in joint aperture during tunnel filling to assess the impact of pressure transients. This is because, the final hydraulic aperture after tunnel filling governs the flow, pore pressure and the transmission delay (timelag) during a transient. Hydro-mechanical behaviour during tunnel filling and pressure transients affects rock mass around the tunnel in different ways. Tunnel filling event impacts a larger area of rock mass around the tunnel and even affects the local groundwater level in the vicinity (Vigl and Gerstner 2010). A very small but steady outflow is maintained through the joint network at the end of filling up to a pre-designed static head. However, transients have short period of application and thus the effect travels only within shallow depths as compared to filling and causes joint deformations in a localized area. Hence, the results are explained by differentiating the effect of these two phases of tunnel operation. Figure 7a, b and Table 5 present normal and shear displacements occurred during tunnel filling for various static heads and maximum normal and shear displacements during respective transient events. The solid lines in Fig. 7a, b show joint normal and shear displacements at different static heads along the length of the joint. During filling, it is observed that, both normal and shear displacements are relatively small for up to 400 m head but then increase significantly for higher heads. As per Benson (1989), the required FoS for static (steady) condition is 1.3 and no FoS is required for water hammer. The FoS for steady state at 400 m head is 1.6 and then reduces to 1.27 and 1.06 for 500 and 600 m, respectively. Hence, according to current design practice, 500 m case is at the boundary of allowable limit whereas 600 m case is an unsafe design. Therefore, the observed joint deformations are as expected since there is large joint shear and dilation when the FoS is equal to higher than allowable limit (Table 5).
During respective transient events, the joint deformations seem to be gradually increasing with static head. However, it is seen that the percentage increase in joint deformations at different static heads Fig. 7c, d is non-linear. It is seen that, they are the highest at 300 m static head, when the FoS is close to 2. The change in normal deformation is the lowest at 500 m, (FoS 1.3) and then again begins to increase at 600 m static head, which indicates hydraulic jacking of the joint since the factor of safety is less than 1 for transient in this case. This jacking, however, occurs for a very short time since the time of application of this pressure transient is small and hence does not travel deeper into the rock mass.
The reason for this behavior where the impact of pressure transients is higher in the mid ranges of static head is intuitively linked to the joint aperture at the end of tunnel filling. At smaller static heads, the aperture after tunnel filling is still small, causing high flow resistance during a transient event, which then causes joint deformation as a result of hydraulic-mechanical interaction. For larger static heads, the joint aperture after tunnel filling is larger, which allows a relatively higher flow out through the joint and thus causing less seepage forces acting on the joint surface and hence smaller joint deformation during transient. Further, it is observed that the effect of pressure transients reduces as we go deeper into the rock mass but there were no joint deformations at areas deeper than 6 m.
Hydraulic impact along the joint during pressure transients at varying static heads is presented in Fig. 8. It is seen that the hydraulic impact increases almost linearly until 400 m along the whole joint length. The highest impact is between 2 to 4 m and the curve is flat at locations deeper than 5 m because the transient does not travel beyond this point. At 500 m, the increase in hydraulic impact is smaller at locations close to the tunnel wall. This is because, large joint aperture after filling causes less flow resistance, leading to smaller timelag as compared to  filling phase itself. Similar behaviour is seen in 600 m case but is more pronounced due to even higher joint aperture. Joint deformation magnitudes during pressure transients are very small when the FoS is at a safe limit (Table 5) and are thought to be insignificant in relation to tunnel stability as per current design practices. However, these observations may have some implications in relation to design consideration regarding long-term operation scenario, when considering the effect of pressure transients, which will be presented in the discussion section.

Effect of Joint Stiffness
Similar simulations as in the preceding sections were carried out using the lower range of stiffness values presented in Table 3 and the normal and shear joint stiffness values of 4 and 1.5 GPa/m, respectively, are referred to as stiffness ST2 hereafter. The results are presented in Fig. 9 and Table 6 and show similar behaviour but with different magnitudes. For static heads up to 400 m, the normal deformations during filling are higher than the high stiffness case, which is logical (Tables 5, 6). However, for 500 and 600 m, it is seen that the deformations are smaller than the high stiffness case. This seems like an anomaly but can be explained by the fact that in low stiffness case, the flow can travel deeper into the rock mass, causing higher flow out of the system. It is seen for 500 m, the maximum flow out of the tunnel for low stiffness case is almost ten times higher than the high stiffness case (0.005 and 0.04 L per second for ST1 and ST2, respectively). Such higher flow out of the system causes pressure release during filling, thus causing smaller deformation as compared to high stiffness case.
Similar to high stiffness case, the percentage change in deformation due to transient is higher in the mid ranges within safe range of FoS but the magnitude of such change is higher because of low joint stiffness. Figure 10 shows hydraulic impact along the joint at different static heads for low stiffness case. It is seen that the pressure transient travels deeper into the rock mass than the high stiffness case and thus the hydraulic impact is more evenly distributed over a larger area. The peak hydraulic impact is at a deeper location as compared to the high stiffness case. Up to 400 m static head, the maximum value of hydraulic impact is almost equal for both cases. But for 500 and 600 m cases, this value is slightly smaller. This is because of the high joint aperture after filling which allows, more outflow in these static heads and lesser pressure buildup during transients.
Three major conclusions can be drawn from the simulations conducted at varying static heads with different stiffness values.
1. The change in joint normal displacement during transient is higher in the mid ranges of static head within acceptable FoS. Hence, pressure transients affect the rock joints more when the FoS within 1.5-2.5. 2. For low stiffness cases, the transient travels deeper into the rock mass and thus affects a larger area of the rock mass as compared to a high stiffness case. The magnitude of peak hydraulic impact is almost the same as high stiffness case but is distributed over a large area. For high stiffness case, the effect of pressure transient is confined to a smaller area. 3. Results indicate that hydraulic impact caused due to short pressure transients are between 2 to 6 m depth. This area seems of greater significance in relation to contribution of fast transients to potential block falls. However, the depth of impact is greatly affected by the time of application of pressure transient. A larger area of rock mass will be affected by a slower pressure pulse; i.e. mass oscillation. The next section presents results of simulations with various time periods of oscillation.

Effect of Time Period of Pressure Transient
Simulations were run for four different time periods; i.e. 3.5, 60, 120 and 200 s at a static head of 100 m for both stiffness cases, to analyze the effect of time period of oscillation. The peak amplitude of the pressure pulse is kept constant at 114.5% of the static head for all time periods. The main objective of these simulations is to observe the changes in zone of influence, joint displacements and hydraulic impact with different time periods. The results are presented in Figs. 11 and 12. In Fig. 11, it is seen that increase in time period from seconds to minutes causes a significant increase in both normal and shear displacements. The effect of time period increase from 60 to 200 s, however, is relatively lesser as compared to the same between 3.5 to 60 s. For high stiffness case, the impact is concentrated within 8 m depth and maximum displacement is seen between 2 to 4 m depth. There is no significant effect in rock mass deeper than 8 m. For low stiffness case, similar behavior is seen regarding the effect of time period increase from seconds to minutes. However, the impact is more uniformly distributed though out the joint length and the effect also travels deeper as compared to high stiffness case. In Fig. 12, it is seen that the zone of influence increases significantly with increasing time period since the transient can travel deep into the rock mass. It is also seen that the peak hydraulic impacts due to time period of 60 s and higher are significantly higher than the peak hydraulic impact due to time period of 3.5 s. The peak hydraulic impact is comparatively smaller in case of low stiffness case, because of higher flow out of the system (and pressure release) due to larger hydraulic aperture and deeper zone of influence.
It should be noted that in these simulations, the pressure amplitudes for all simulations are kept constant by varying the time period only. However, water hammer has a higher potential for pressure increase as compared to the mass oscillations. This can be seen in Fig. 2a where the pressure rise due to mass oscillation is 7% of the static head before shutdown and the pressure rise due to water hammer is 14.5%. Hence, simulation is carried out with a pressure amplitude increase of 7% during transient and time period 200 s to compare the actual difference in hydraulic impact between a water hammer and mass oscillation for this Fig. 11 a Normal displacement and b shear displacement for high stiffness case and c normal displacement and d shear displacement for low stiffness case during transients with different time periods Fig. 12 Hydraulic impact along the joint due to pressure transients with different time periods for a high stiffness (ST1) and b low stiffness (ST2) particular waterway system. The results are indicated by dotted lines in Figs. 11 and 12. It is seen that the displacements are still much larger than 3.5 s case, but the hydraulic impact is reduced by almost half in case of both stiffness cases. Hence, the results show that the effect of mass oscillations is much higher than that of water hammer, when compared in terms of a single pressure pulse. However, the cumulative impact is more relevant in terms of block stability over longterm operation, which will be presented in the discussion section. Three major conclusions can be drawn from the above results, which are as follows: 1. The effect in terms of both joint displacements and hydraulic impact increases significantly when the time period increases from seconds to minutes (between a typical water hammer and a mass oscillation event). Joint displacement does not increase significantly when the time period is increased further. The hydraulic impact, however, keeps on increasing when the time period is increased. This is because, the pressure is ¨trapped¨ in the joint for a longer period, even though the aperture remains almost same.
2. In high stiffness case, the joint displacement is smaller but concentrated within shallow depths from the tunnel walls. For low stiffness, the displacements are much larger, and the effect is dispersed more uniformly over a deeper area in the rock mass. 3. The hydraulic impact is smaller in low stiffness case, but the zone of influence is larger.

Effect of Joint Friction Angle
Results of simulations with joint friction angles of 40° and 25° are presented in Figs. 13 and 14. It is seen that joint friction angle has little or no impact in both joint displacement and hydraulic impact when the normal stress across joint is high. With reduced normal stress, the joint displacements during transients increase significantly. Similar results as in previous simulations are seen along the length of joint because of change in joint stiffness. However, at 300 m static head, the joint deformations are significantly high when the friction angle is reduced, and the effect is seen deep into the rock mass for both stiffness cases.

Fig. 13
Normal displacement along the joint due to pressure transients with different joint friction angle for a high stiffness (ST1) and b low stiffness (ST2) Fig. 14 Hydraulic impact along the joint due to pressure transients with different joint friction angle for a high stiffness (ST1) and b low stiffness (ST2) Figure 14 shows that hydraulic impact decreases with decrease in friction angle for higher static heads. This is the result of high apertures formed during the end of tunnel filling. During tunnel filling, smaller friction angle causes larger shear displacement of the joint. Joint dilation during shear displacement causes higher joint aperture at the end of tunnel filling. This phenomenon is more pronounced when the static heads are higher (lower effective normal stresses across joints).
When normal stresses are high, for example at 100 m static heads, the hydraulic impacts are same for both friction angles since hydraulic apertures are also the same.

Effect of Joint Dilation
Simulations for 500 m static head for both stiffness cases were also run without allowing dilation. This static head is chosen to see maximum effect since it is the maximum allowable static head as per the conventional FoS principle. In Fig. 15a, it is seen that when dilation is restricted, the peak hydraulic impact increases significantly and is shifted closer to the tunnel wall, with a narrower zone of influence. This is logical because the hydraulic aperture after tunnel filling (steady state in Fig. 15b), without dilation will be smaller and thus creates larger flow resistance, which prevents the pressure pulse to travel deep into the rock. This causes high pressure buildup and thus higher hydraulic impact closer to the tunnel wall as compared to when dilation is allowed. In Fig. 15a, it is also seen that the effect of dilation is larger in high stiffness case. This is because the onset of dilation is more sudden and pronounced as compared to a low stiffness joint which allows gradual displacement. Higher dilation effect causes large change in hydraulic aperture and hence, larger difference in hydraulic impact as compared to low stiffness joint.

Discussion
In contrast to established understanding that the additional seepage forces during pressure transients (or hydraulic impact) are higher at locations closer to the tunnel, it is seen that the highest impact occurs where there is sufficient pore pressure buildup and also enough flow resistance to cause a significant timelag. Such locations are few meters into the tunnel wall, depending upon joint properties. Further deep into the rock mass, the hydraulic impact starts to decrease as the pore pressure buildup is reduced, even though the timelag is higher. Therefore, in addition to rock mass pore pressure, a more important parameter that needs to be added to this phenomenon is the timelag or transmission delay of the pressure pulse into the rock mass, which depends upon the joint aperture created at the end of tunnel filling.
In relation to tunnel static heads, it is obvious that during water filling, the magnitude of joint deformation will increase with increasing static head. However, it is seen that the percentage increase in joint deformation during pressure transient is found to be the highest when the FoS is between 1.5 and 2.5. This means that tunnels designed within this FoS range will be the one that are most impacted by pressure transients. Palmstrøm (1987) shows that critical locations along the unlined tunnels usually have FoS lying within this range. FoS larger than 2.5 are usually undesirable with respect to economic reasons, whereas FoS smaller than 1.5 are uncommon, since they would require specific stress measurements during construction.
Block failures occur at a certain amount of absolute joint displacement. Such absolute displacement is the cumulative effect of a number of pressure transients over the tunnel operation period, referred to as hydraulic fatigue. For such fatigue to occur due to transients, the joint has to deform further from its initial deformation state at static tunnel pressure or steady state. At the static condition, when the FoS Fig. 15 Effect of joint dilation at 500 m static head on a hydraulic impact for both stiffness cases and b max. normal displacement along the joint for high stiffness case during steady state and pressure transient is greater than 1, the block will remain stable. This applies for all static head cases in Fig. 7, even if absolute deformations are higher for higher static heads. During transients, the more a joint deforms from its initial displacement state, the more it contributes to fatigue. As seen in Fig. 7, the cases with FoS between 1.5 and 2.5 show larger relative displacement during transients and hence will have more "loosening effect" on the joint. Hence, percentage increase of deformation during pressure transients is used for comparison because it shows which case of FoS shows more relative joint displacement from their respective steady states after tunnel filling.
This shows that most of the unlined tunnels which are under operation may be impacted significantly in terms of fatigue due to cyclic loading of joints over long-term operation of tunnels, eventually leading to block falls.
This conclusion holds true for both hard crystalline rocks and schistose rocks as represented by two sets of joint stiffness values used in the simulations. Figure 16 shows a typical example of block falls in rock mass with varying joint conditions registered during headrace tunnel inspections in Norwegian power plant under operation over 30 years.
A commonly established concept about impact of time period of pressure transient is that only slower pressure oscillations, i.e., mass oscillations are of significance, when it comes to instabilities occurring during plant operation. Prevalent argument for this is that the time period of oscillation for water hammer is too short to travel deep into the rock mass. However, from the field data and simulations, it is seen that such transients can also travel up to 4 m deep into the tunnel wall, even in stiff joint conditions and sufficiently high normal stresses. Hence, the next question is whether the displacements and hydraulic impact are significant to contribute to block falls. It is seen that the maximum hydraulic impact along the joint due to each cycle of mass oscillation is about 25 times higher than the impact due to one water hammer pulse for high stiffness case and 15 times higher for low stiffness case. The displacements for both cases due to mass oscillation are only about 2 times higher than water hammer. At higher static heads, when the FoS is lesser but still within allowable range, it is likely that hydraulic impact ratio will be almost equal or slightly less since mass oscillation may cause increased leakage and pressure release because of longer time period. However, this still supports the previously established knowledge that mass oscillations cause significantly larger hydraulic impact as compared to water hammers.
The hydraulic impact ratio can vary depending upon different parameters, such as waterway lengths, area, valve closure time, etc. Longer tunnel length between turbine and surge tank (pressure tunnels/shafts) contributes to higher water hammer magnitude and longer tunnel length between surge tank and reservoir (headrace/ tailrace tunnels) causes larger mass oscillations. For Roskrepp waterway, these two lengths are 560 m and 3.5 km, respectively.
It is desirable to place the surge shaft as close as possible to the turbine to reduce the detrimental effect of large water hammer. This means that most waterways have relatively longer headrace/tailrace tunnels and shorter pressure tunnels/shafts. Further, the pressure increase due to water hammer is only limited within the length of pressure tunnels/ shafts, but the pressure rise due to mass oscillations applies to the full length of waterway. This means the hydraulic impact of mass oscillations are much large in magnitude and felt deeper into the rock mass. Also, they apply to longer stretch of tunnels as compared to water hammer. Faster valve closure time may contribute to large water hammer magnitude but still the hydraulic impact is smaller because of their small time period. However, the effect of water hammer may not be fully ignored. This is because, for every shutdown event, the number of pressure pulses due to water hammer with significant timelag is higher than mass oscillations. In the case of Roskrepp, for every shutdown event, there are 3-4 mass oscillation cycles with timelag whereas the Fig. 16 Block falls registered in a hard crystalline rock (Trondhjemite) and b schistose phyllite during inspections of unlined pressure tunnel of a power plant in Norway number of water hammer pulses is much higher (Fig. 17a). Further, during a pressure transient, it is very difficult to make a distinction to when water hammer ceases and mass oscillation start. Hence, the cumulative hydraulic impact of all such pulses is significant regarding long-term stability of waterways systems.
The effect of small but frequent pressure pulses, such as water hammer on the opening of joints, can be compared to the concept of fatigue hydraulic fracturing (FHF) introduced by Zang et al. (2013), to enhance the permeability of fractured rocks in petroleum reservoirs. Hofmann et al. (2018) describes the cyclic soft stimulation (CSS) technique which consists of injection protocol with three types of cycles with different time periods, i.e. Long-term cycles (LTCs, hours or more), medium-term cycles, (MTCs, minutes to hours) and short-term cycles (STCs, minutes and less).
In Fig. 17a, the duration between two valve closure events can be compared to the LTCs, mass oscillations are equivalent to the MTCs and water hammers can be compared to the STCs. In CSS, the short-term cycles are applied to amplify the fatigue and weakening of the rock by inducing additional small fissures before and besides the macroscopic fracture development. The time scale and magnitude of pressures in these different two processes can vary but the basic principle in both cases is hydraulic fatigue.
Hence, in the authors' opinion, water hammers, wherever applicable along the waterway, must be taken in consideration to assess long-term stability of water tunnels subjected to pressure transient. Therefore, it is important to consider the cumulative impact of both water hammer and mass oscillation because the impacts due to water hammers at one cycle, even though small, are measurable and more frequent, as shown by these experiment and simulations. A deeper understanding of hydraulic fatigue is necessary to quantify the effect of water hammers. Till recent time, very limited research has been conducted in this field. A much larger number of theories and experiments exist regarding ¨mechanical fatigue¨ of rock without fluid injection. But whether and how classical fatigue concepts can be applied to hydraulic fatigue remains a subject of further research.
It can still be stated that fatigue causes gradual deterioration of joint surfaces over long-term resulting in reduction of joint stiffness and friction angle. Simulation results with decreasing stiffness and friction angle can provide some insight on how the behaviour changes with change in joint properties. Table 7 presents some results for 300 m static head (FoS 2.1) for different stiffness and friction angles.
The values in Table 7 indicate that when the joint stiffness gradually deteriorates from ST1 to ST2, transients cause much larger joint deformation (7.5 times) than tunnel filling (3.5 times) for 40° friction angle case. Similar results are seen for 25° friction angle as well. But for the same stiffness, when friction angle decreases from 40° to 25°, larger joint deformation occurs during tunnel filling itself and thus the effect during transient is relatively small. This shows that transients tend to cause larger deformation as the joint Fig. 17 Comparison between a recorded pressure fluctuation in an unlined hydropower tunnel and b cyclic fluid injection protocol used in petroleum and geothermal reservoirs (Hofmann et al. (2018)  3.14 2.12 1.63 1.43 undergoes gradual reduction in stiffness. But for gradual reduction of friction angle, joints are likely to deform more during filling itself. In real scenario, a simultaneous reduction of these parameters is a more likely case. Hence, a possible combined effect would be such that transients with small time periods, cause ¨gradual loosening¨ of joints in addition to the effect of mass oscillation, thereby reducing the joint strength. Macroscopic movement of blocks occurs during filling/dewatering and mass oscillations. The orientation of the joint with respect to major principal stresses also plays an important role in the behaviour during tunnel filling and pressure transient because it affects the normal stresses across the joint. If the minimum principal stress is parallel to the joint, the normal stress across joint will be higher because the major or intermediate principal stress will be acting perpendicular to the joint. In such condition, the FoS will be higher and hence the joint deformation during tunnel filling and transient will be reduced. This may result in a non-conducting, tightly closed joint which shows very little or no response during pressure transient. On the other hand, if the minimum principal stress is perpendicular to the joint, this will create a larger joint displacement because of reduced normal stress. If such normal stress is critically low (lower than tunnel water pressure), then the joint will be hydraulically jacked during tunnel filling itself.

Conclusion
The results of numerical simulation show that relative joint deformations due to short pressure transients are highest when normal stresses acting across the joints are 1.5-2.5 times higher than the static tunnel water pressure. Since most critical locations in the tunnels are designed to be in this range of factor of safety, it is concluded that such tunnels are significantly affected by pressure transients, eventually leading to hydraulic fatigue. This conclusion holds true for varying joint stiffness values.
The results also confirm that mass oscillations can cause significantly large hydraulic impact on block stability. It has been demonstrated that for a typical hydropower waterway system, the hydraulic impact ratio between mass oscillation and water hammer per cycle can vary between 15 and 25. This ratio can vary, mainly depending on the length of headrace/tailrace tunnel and high-pressure tunnel/shaft. Simulation results show that joint deformation does not increase significantly, when the time period of oscillation is increased from one minute to several minutes. The higher impact with increasing time period is mainly because the pressure is ¨trapped¨ for a longer time.
More importantly, the analysis results show that water hammers, wherever applicable along the waterway can contribute to gradual ¨loosening¨ of the joints since they can travel up to 4 m into the rock mass even in stiff joint conditions and sufficiently high normal stresses. As indicated by results with various stiffness and friction angles, such ¨loosening¨ becomes more prominent when the stiffness of joint gradually reduces. This phenomenon, combined with reduced friction angle over long term, causes larger displacements or block movements during tunnel dewatering/filling and pressure oscillations with large time period, i.e., mass oscillations. Thus, it is concluded that water hammers and mass oscillations have a cumulative effect on the long-term stability of blocks. However, the effect of water hammer is only limited in the waterway length between the surge shaft and turbine but effect of mass oscillations applies to the full length of waterway.
The analyses and results presented above are based on numerical models assuming the rock joint as the interface between two parallel plates. Limitation of these models is that it does not consider the effect of flow tortuosity in rock joints caused by joint roughness, which is still an outstanding issue in numerical simulation of flow processes in rock. However, the achieved results are consistent with basic phenomenon of joint fluid flow and some well-established trends from previous researches. Hence, these findings can be of significance to further understand the issue of block failures caused by frequent pressure transients. Another issue of interest is the ¨hydraulic fatigue¨ due to which the joint properties change gradually over time, which is a large field of research.
Long-term stability of unlined hydropower water tunnels subjected to medium to high static water heads in changing energy market conditions is an immerging challenge, since frequent start stop sequences will cause intensified occurrence of pressure transients. In authors' knowledge, this is the first time such study has been conducted involving hydro-mechanical interactions of frequent pore pressure changes in the rock mass around hydropower tunnels. The effect of frequent pressure transients of both short-and long-time periods (water hammer and mass oscillation) are the main cause for frequent block falls in unlined pressure tunnels.
Funding Open access funding provided by NTNU Norwegian University of Science and Technology (incl St. Olavs Hospital -Trondheim University Hospital)..

Conflict of interest
The authors wish to confirm that there are no known conflicts of interest associated with this publication and there has been no financial support for this work that could have influenced its outcome.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.