## Appendix B Numerical Method  
The equations describing the water hammer phenomenon, which are used in this work, have the following form [19]:  
* Continuity equation: \[\frac{\partial H}{\partial t}+\frac{a^{2}}{g}\frac{\partial V}{\partial x}=0\] (18)
* Momentum equation: \[\frac{\partial V}{\partial t}+g\frac{\partial H}{\partial x}+g\cdot J=0\] (19) where: \(V\)--flow velocity in the pipe [m/s];\(H\)--pressure head [m wc];  
\(J\)--unit pressure losses [-];  
\(a\)--pressure wave speed [m/s];  
\(x\)--coordinate along the axis of the pipe [m];  
\(t\)--time coordinate [s];  
\(g\)--acceleration of gravity [m/s\({}^{2}\)]. The component of unit pressure losses can be represented using the Darcy-Weisbach formula:  
\[J=f\frac{1}{D}\cdot\frac{V|V|}{2g} \tag{11}\]  
where the friction coefficient, \(f\), is defined depending on the friction loss model used in the calculation. The following models were used in this work:  
* Quasi-steady model, in which the value of the coefficient \(f\) is calculated using: \(\odot\) Hagen-Poiseuille law for laminar flows Re \(\leq\) 2320 [58]: \[f=\frac{64}{Re}\] (12)
* Colebrook-White formula ([PERSON] 1939) for transient and turbulent flow conditions \(Re>\) 2320: \[\frac{1}{\sqrt{f}}=-2\text{log}\left(\frac{2.51}{Re\sqrt{f}}+\frac{k_{r}}{3.71D}\right)\] (13) where: \(Re\)--Reynolds number: \(Re=VD\rho/\mu\) [-]; \(\rho\)--liquid density [kg\(\cdot\)m\({}^{-3}\)];  
\(\mu\)--dynamic viscosity coefficient [Pa\(\cdot\)s];  
\(k_{r}\)--absolute roughness of the inner surface of the pipe [m];  
* Unsteady friction models: \(\odot\) The Vardy and Brown model [52], in which the loss term is calculated based on the history of flow changes. The range of applications for this model includes smooth pipes and turbulent flows. The model has been positively verified by the authors in [3]. The efficient version of this model, prepared by [PERSON] et al. [53], which allows for a significant reduction in the requirements for computing, was considered in the calculations. According to this model, the unit friction loss term in equation (11) has the following form: \[J(t)=f\frac{1}{D}\cdot\frac{V(t)|V(t)|}{2g}+\frac{16\mu}{\rho gD^{2}}\int_{0}^ {t}\frac{\partial V}{\partial t^{*}}W(t-t^{*})\text{d}t^{*}\] (14) where: \(W\)--weighting function defined in the domain of dimensionless time \(\tau=4\mu t/\rho D^{2}\) in the form: \[W(\tau)=\frac{A^{*}e^{-B^{*}\tau}}{\sqrt{\tau}}\] (15) for: \(A^{*}=\frac{1}{2}\sqrt{\frac{\tau_{w}}{\pi\
u_{lam}}}\), \(B^{*}=\frac{Re^{*}}{12.86}\), \(\kappa=\text{lg}\left(15.29\cdot Re^{-0.0567}\right)\), \(\gamma_{w}\)--kinematic viscosity at the wall of the pipeline [m\({}^{2}\)/s],\(\
u_{lam}\)--laminar kinematic viscosity [m\({}^{2}\)/s]. In the proposal of [PERSON] et al., the weighting function \(W\) is approximated using the following expression:  
\[W(\tau)\approx W_{app}(\tau)=\sum_{k=1}^{N}A^{*}m_{k}^{*}e^{-(n_{k}^{*}+B^{*})\tau} \tag{10}\]  
where the coefficients \(m_{k}^{*}\) and \(n_{k}^{*}\) for \(k=(1,2,\ldots,10)\) take the values from Table 1:  
\(\circ\)  
In the model of [PERSON] et al. (2018, 2019), the loss term is calculated on the basis of instantaneous flow velocity and flow acceleration values. According to this model, the unit friction loss term in Equation (11) has the following form:  
\[J(t)=f\frac{1}{D}.\frac{V(t)|V(t)|}{2g}+\frac{k_{3}}{g}\left(\frac{\partial V} {\partial t}-a\frac{\partial V}{\partial x}\right) \tag{11}\]  
Two versions of this model were considered in the calculations:  
In the first (universal) version, the \(k_{3}\) coefficient, which is the basic element of this model, is calculated from the formula proposed by [PERSON] and [PERSON] (2018), which makes its value dependent on the flow rate value:  
\[\frac{\partial V}{\partial t}+g\frac{\partial H}{\partial x}+g\cdot J=0\]0  
In the second (the most popular) version, the value of the \(k_{3}\) coefficient is selected in a way that guarantees the best possible adjustment of the calculation course to the measured course.