**Implementation of a Local Time Stepping Algorithm and Its Acceleration E**ff**ect on Two-Dimensional Hydrodynamic Models**

#### **Xiyan Yang, Wenjie An, Wenda Li and Shanghong Zhang \***

Renewable Energy School, North China Electric Power University, Beijing 102206, China; 120192211870@ncepu.edu.cn (X.Y.); 120192211845@ncepu.edu.cn (W.A.); 120192111184@ncepu.edu.cn (W.L.) **\*** Correspondence: zhangshanghong@ncepu.edu.cn

Received: 11 February 2020; Accepted: 7 April 2020; Published: 17 April 2020

**Abstract:** The engineering applications of two-dimensional (2D) hydrodynamic models are restricted by the enormous number of meshes needed and the overheads of simulation time. The aim of this study is to improve computational efficiency and optimize the balance between the quantity of grids used in and the simulation accuracy of 2D hydrodynamic models. Local mesh refinement and a local time stepping (LTS) strategy were used to address this aim. The implementation of the LTS algorithm on a 2D shallow-water dynamic model was investigated using the finite volume method on unstructured meshes. The model performance was evaluated using three canonical test cases, which discussed the influential factors and the adaptive conditions of the algorithm. The results of the numerical tests show that the LTS method improved the computational efficiency and fulfilled mass conservation and solution accuracy constraints. Speedup ratios of between 1.3 and 2.1 were obtained. The LTS scheme was used for navigable flow simulation of the river reach between the Three Gorges and Gezhouba Dams. This showed that the LTS scheme is effective for real complex applications and long simulations and can meet the required accuracy. An analysis of the influence of the mesh refinement on the speedup was conducted. Coarse and refined mesh proportions and mesh scales observably affected the acceleration effect of the LTS algorithm. Smaller proportions of refined mesh resulted in higher speedup ratios. Acceleration was the most obvious when mesh scale differences were large. These results provide technical guidelines for reducing computational time for 2D hydrodynamic models on non-uniform unstructured grids.

**Keywords:** two-dimensional hydrodynamic model; local time stepping; unstructured grids; numerical simulation; computational efficiency

#### **1. Introduction**

The numerical simulation of shallow water flow is frequently used in flood forecasting, river regulation, and flood-control planning [1]. Given the need for better accuracy and breadth of shallow-water simulations in engineering applications, a balance between the number of cells, simulation accuracy, and computational time needs to be found. Addressing this issue using optimization of the solving algorithm [2–4] has received much research attention, as has the use of computer hardware acceleration and parallel computing technology [5–7].

The local mesh refinement technique is especially efficient for balancing the number of grids and the accuracy of the simulation. It is flexible and allows grids to be arranged according to different engineering concerns, while retaining a high level of refinement in regions requiring detail, such as where flows change sharply or where there are key features. A coarse mesh is used in general or slow-flowing areas. Consequently, the simulation accuracy can be ensured with only a small increase in the number of grids. This technique has been broadly used in engineering applications involving shallow-water flow simulations [8–12]. However, the time step must strictly satisfy the model solution, which is obtained using the Courant–Friedrichs–Lewy (CFL) stability condition. When using a locally refined mesh, the global time stepping (GTS) size is limited by the time step of the refined mesh, which markedly increases the overall computation time.

In response to this issue, Osher and Sanders [13] proposed a local time stepping (LTS) algorithm for solving the one-dimensional (1D) scalar conservation equations. In each cell, a locally allowable maximum time step satisfying the CFL stability condition was adopted to minimize the computational time. This effectively deals with the complexity of the lag caused by the heterogeneous time step. Tan and Huang [14] developed an efficient adaptive mesh refinement (AMR) with an LTS algorithm for 1D nonlinear hyperbolic problems. Crossley et al. [15,16] successfully applied the LTS technique to the modeling of open channel flows using 1D Saint-Venant equations; their results for unsteady shallow-water equations with LTS met the required accuracy. Consequently, the use of LTS algorithms has enabled significant progress in 1D hydrodynamic simulations. Compared with 1D models for simulating long and wide rivers, two-dimensional (2D) hydrodynamic models are better suited to complex terrains and provide more precise simulation results. Dazzi et al. [17,18] used AMR to implement an LTS algorithm in the numerical simulation of 2D shallow-water flow for a steady-state test case, a circular dam break case, and a Tacker test case. Using the time step of each cell and the CFL condition, they evaluated whether the cell's current time step met the CFL condition. If the condition was not met, then the rule of degrading the time-step rank was implemented, which involved a dynamic inspection of the time step for each cell over the whole computational domain. Wu et al. [19] adopted the LTS technique to avoid excessive storage when modeling three-dimensional (3D) fluid dynamics using the discontinuous finite element method. The powerful function of 3D models determines their complex and elusive performance. At present, hydrodynamic modeling is dominated by 2D modeling, which can provide sufficiently reliable results, such as for river navigation.

Nevertheless, recent research has paid increasing attention to this aspect with a suitable treatment of unstructured grids and has yielded promising results. Although the numerical procedure has unique advantages in terms of mesh generation, storage, use of elements, capturing of specific physical phenomena, and computational efficiency, the boundaries of such regular grids are too generalized and are not geometrically adaptable. Furthermore, the computational accuracy of the boundaries is extremely rough and neglects terrain changes, and is thus limited in its scope. The use of unstructured grids provides highly adaptable solutions for complex terrains. Unstructured grids can incorporate a region's boundaries and constraints using mesh density control in continuous regions. Consequently, unstructured grids play a dominant role in grid generation technology. Using the Runge–Kutta discontinuous Galerkin finite element method, Trahan and Dawson [20] applied LTS to the shallow-water equations and incorporated rainfall, wetting, and drying into the model. Using LTS, Hu et al. [21] established a new coupled model, which had high computational efficiency, for sediment-laden flows. The computational domain was discretized using unstructured triangular meshes. The inter-cell numerical fluxes were estimated using the Harten–Lax–van Leer–Contact (HLLC) approximate Riemann solver. The model used by Hu et al. was applied to the Taipingkou waterway of the Yangtze River and achieved a 92% reduction in computational cost without loss of quantitative accuracy. To solve the problem of mass conservation linked to time splitting in LTS algorithms, Michael [22] improved the finite volume method (FVM) using a deformation correction to achieve local mass conservation, while a row-oriented sub-stepping was proposed to deal with Courant numbers larger than one.

Despite several successful studies based on the LTS algorithm, the computation of flux at the cells with different time-step levels has not been dealt with clearly in recent research. The calculation of the unit variables at the interface is difficult because the use of the local time-step invariably leads to time faults. The calculation of element variables at the interface is crucial but is difficult because of time splitting errors in the LTS algorithm. Implementation of the LTS algorithm on unstructured grids is also valuable for advancing the knowledge base. No specific comparison of the acceleration efficiency of the LTS algorithm has been demonstrated in the existing research. To address these issues, the LTS technique was applied in the establishment of a high-efficiency, 2D, shallow-water, dynamic model using the FVM and unstructured grids. Experiments were performed to verify the flux approximation at the interface, to assess the accuracy of the results, and to analyze the feasibility analysis of the LTS algorithm for complex applications. Model evaluation was performed using 2D simulations of two instantaneous dam break test cases and engineering applications that simulate the navigable flow of the river reach between the Three Gorges and Gezhouba Dams. These tests show the acceleration effect, influencing factors, and the adaptive conditions of the LTS algorithm.

#### **2. Methods**

#### *2.1. Global Time Stepping Scheme*

#### 2.1.1. Governing Equation

The governing equations used herein are the 2D shallow-water equations (SWEs), written in conservation form [23]:

$$\frac{\partial \mathcal{U}}{\partial t} + \frac{\partial F(\mathcal{U})}{\partial \mathbf{x}} + \frac{\partial G(\mathcal{U})}{\partial y} = S\_0(\mathcal{U}, \mathbf{x}, y) + S\_f(\mathcal{U}, \mathbf{x}, y) \ (\Omega \times [0, T\_s]). \tag{1}$$

,

In the present study, the modified form of the SWEs was selected to guarantee that the scheme was well-balanced. The variables follow from

$$\begin{aligned} \mathbf{U} &= \begin{bmatrix} h \\ hu \\ hv \end{bmatrix} F(\mathbf{U}) = \begin{bmatrix} hu \\ hu^2 + \frac{1}{2}gh^2 \\ huv \end{bmatrix} G(\mathbf{U}) = \begin{bmatrix} hv \\ huv \\ hv^2 + \frac{1}{2}gh^2 \end{bmatrix}, \\\ \mathbf{S}\_0 &= \begin{bmatrix} 0 \\ -gh\frac{\partial z\_b}{\partial x} \\ -gh\frac{\partial z\_b}{\partial y} \end{bmatrix}, \mathbf{S}\_f = \begin{bmatrix} 0 \\ -gh\frac{u^2u\sqrt{u^2 + v^2}}{h^4/3} \\ -gh\frac{u^2v\sqrt{u^2 + v^2}}{h^4/3} \end{bmatrix}. \end{aligned}$$

where *U* is the vector of conserved variables; *F* and *G* are the tensors of fluxes in the *x* and *y* directions, respectively; *S*<sup>0</sup> and *Sf* are the bed and friction slope source terms, respectively; *x* and *y* represent the spatial horizontal coordinates; *u* and *v* represent the average velocity components from integration of water depth in the *x* and *y* directions, respectively; *h* is the flow depth; *zb* is the bed elevation above the datum; *g* represents the acceleration related to gravity; *n* is Manning's roughness coefficient; *Ts* is the duration.

#### 2.1.2. Numerical Technique

Because of the nonlinear characteristics of the SWEs, the analytical solution cannot be obtained generally but can be approximated by numerical discretization. The FVM was used to discretize the governing equation. The computational region was discretized into several points. With these points as the center, the whole computational domain was divided into several interconnected but non-overlapping control volumes. By integrating the basic equation with each control volume, a set of algebraic equations with the unknown quantity on the calculation node was obtained. A triangular mesh was used to discretize the computational domain to adapt the boundary conditions for a complex terrain. The cell-centered mode was used in the calculations, which means that the physical variables are defined in the center of the control volume [24], as shown in Figure 1.

**Figure 1.** Schematic diagram of the two-dimensional model with finite volume discretization: physical variables are defined at the grid center, where Ω*<sup>k</sup>* is the *k*th control volume; *n* is the unit normal vector; *p*(*i*, *j*) is a grid point; and Δ*lkj* is the length of a side.

Integration of Equation (1) over the area of the governing volume Ω, gives

$$\int\_{\Omega\_i} \left( \frac{\partial \mathcal{U}}{\partial t} + \nabla \cdot \mathbf{E} \right) d\Omega = \int\_{\Omega\_i} \left( \mathbb{S}\_0 + \mathbb{S}\_f \right) d\Omega. \tag{2}$$

In Equation (2), **E** = [**F**, **G**] represents a 2D matrix.

The volume integral is transformed to a line integral along the edge of the control volume using Gauss' formula 

$$\frac{\Delta\mu\_k}{\Delta t}\Delta\mathbf{s}\_k = -\int\_{l\_k} \mathbf{E} \cdot \mathbf{n} dl + \int\_{\Omega\_k} \left(\mathbf{S}\_0 + \mathbf{S}\_f\right) d\Omega\_\prime \tag{3}$$

where *Uk* is the average value of the governing cell; *lk* is the edge of the *k*th control volume Ω*k*; Δ*sk* is the area of each calculation cell; **n** is the unit vector of the normal direction outside the edge.

After discretization of Equation (3), the following formula is obtained:

$$
\Delta Ul = -\frac{\Delta t}{\Delta \mathbf{s}\_k} \sum\_{j=1}^{3} \left( \mathbf{E}\_{kj}, \mathbf{n}\_{kj} \right) \Delta l\_{kj} + \frac{\Delta t}{\Delta \mathbf{s}\_k} \overline{\mathbf{S}}\_{\prime} \tag{4}
$$

where Δ*lkj* is the length of each side (using triangular grids, *j* = 1, 2, 3); and **E***kj* and **n***kj* represent the numerical fluxes and the external normal unit vector of edge *<sup>j</sup>*. *<sup>S</sup>*<sup>0</sup> <sup>=</sup> Ω*i <sup>S</sup>*0*d*<sup>Ω</sup> and *Sf* <sup>=</sup> Ω*i Sf d*Ω = Δ*skSf* represent the integral values of the bed and friction slope source terms in the mesh cell, respectively. The notation **<sup>E</sup>***kj*, **<sup>n</sup>***kj* denotes the inner product.

Commonly used approximate Riemann solvers for the discontinuous problem include the Osher, Roe, Harten–Lax–van Leer (HLL), and HLLC schemes [25]. In the present study, the fluxes at the interface of the governing cell were estimated using the Roe scheme. The physical variables *UL* and *UR* on the interface of the left and right grid were determined to solve the Riemann discontinuity problem. The flux expression is defined as follows [26–29]:

$$\mathbf{E} \cdot \mathbf{n} = \frac{1}{2} [ (\mathbf{F}, \mathbf{G})\_{R} \cdot \mathbf{n} + (\mathbf{F}, \mathbf{G})\_{L} \cdot \mathbf{n} - |\widetilde{f}| (\mathcal{U}\_{R} - \mathcal{U}\_{L}) ],\tag{5}$$

where *<sup>J</sup>* <sup>=</sup> <sup>∂</sup>(**E**·**n**) <sup>∂</sup>*<sup>U</sup>* is the Roe average of the Jacobian matrix. The three eigenvalues of the Jacobian matrix are λ*k*(*<sup>k</sup>* <sup>=</sup> 1, 2, 3), while the related eigenvectors are *ek*. The left and right interface conservative form variables are *UL* and *UR*, respectively.

To improve the adaptability of the model to complicated terrain conditions, the source term reflecting the bottom topography was managed using an eigen-decomposition method, while upwind treatment was applied to balance the interface flux [30,31]:

$$\overline{S\_0} = \sum\_{j=1}^{3} \sum\_{k=1}^{3} \left[ \frac{1}{2} (1 - \text{sign}(\widetilde{\lambda}\_k)) \beta\_k \widetilde{\overline{c\_k}} l\_j \right]^j,\tag{6}$$

where *sign*() is the sign function. In this case, <sup>β</sup><sup>1</sup> = <sup>−</sup><sup>1</sup> <sup>2</sup>*cZ b*, <sup>β</sup><sup>2</sup> <sup>=</sup> 0, and <sup>β</sup><sup>3</sup> <sup>=</sup> <sup>1</sup> <sup>2</sup>*cZ b*. Furthermore, *<sup>c</sup>* is the average velocity given by the Roe scheme; and *Z <sup>b</sup>* <sup>=</sup> <sup>1</sup> <sup>2</sup> [(*Zb*)*<sup>L</sup>* + (*Zb*)*R*]; while (*Zb*)*<sup>L</sup>* and (*Zb*)*<sup>R</sup>* are the bottom elevations of the left and right side, respectively.

An explicit scheme that discretizes the frictional slope source term was carried out to increase the stability of the computation. Eventually, the variation of the conservation after undergoing Δ*t* was obtained.

$$
\Delta lI = \mathbf{I} \frac{\Delta t}{A} \left[ (-\sum\_{j=1}^{3} \mathbf{E}\_{j} \mathbf{n}\_{j} l\_{j}) + \sum\_{j=1}^{3} \sum\_{k=1}^{3} \left( \frac{1}{2} (1 - \text{sign}(\widetilde{\lambda}\_{k})) \beta\_{k} \widetilde{\mathbf{e}}\_{k} l\_{j})^{j} + \Delta t S\_{f}^{\text{u}} \right] \tag{7}
$$

where **I** is the unit matrix. Superscripts *n* and *n* + 1 refer to the current and updated time levels, respectively. Δ*t* is the globally permissible time step, which must be computed according to the CFL stability condition [32]:

$$
\Delta t\_i^{\mathcal{S}} = \text{Cr}\_{j=1,2,3}(\frac{d\_{i,j}}{\sqrt{u\_i^2 + v\_i^2} + \sqrt{g h\_i}}),
\tag{8a}
$$

$$
\Delta t = \min \left( \Delta t\_i^{\mathcal{J}} \right) \tag{8b}
$$

where *i* is the computational cell index, with values of 1, 2, ... , *Ne*. The term *Ne* refer to the total number of computational grid elements; *di*,*<sup>j</sup>* is the distance from the center of the cell (triangle centroid) to the three sides; *ui* and *vi* are the velocities in the *x* and *y* directions, respectively, of cell *i*; *hi* is the flow depth of cell *i*. The Courant number *Cr* was set to 0.8.

In the traditional GTS algorithm, the time step Δ*t* is equal to the minimum value for the whole computational domain and is proportional to the size of the grid elements. When the scale of each mesh element varies markedly, the uniform time-step restriction will severely reduce the efficiency of the model operation. Therefore, it is essential to explore an efficient method in which the time-step changes with the mesh size and satisfies the CFL conditions.

#### *2.2. Local Time Stepping Scheme*

The key approach in the LTS algorithm is to advance each cell with a time step as close as possible to its maximum allowable value while satisfying the stability of the CFL condition. The algorithm reconstructs the time-step level of each cell with the aim of reducing the computational load of coarse grids, consequently saving computational time. Figure 2 compares the implementation of GTS versus LTS algorithms. The LTS algorithm differs in two aspects: reconstructing the local time step of each cell (Section 2.2.1) and updating values of the conserved variables Δ*U* (Section 2.2.2).

**Figure 2.** Flow diagram showing global time stepping (GTS) and local time stepping (LTS) algorithms. In the GTS scheme, each cell is updated with the same time step Δ*t*, while in the LTS scheme, a locally allowable maximum time step (Δ*t b i* ) is used for updating variables in each cell.

#### 2.2.1. Reconstructing the Local Time Step of Each Cell

This procedure is defined as follows:

(1) Compute the initial time step Δ*t* ∗ *i*

The Δ*t* ∗ *<sup>i</sup>* of each cell satisfying the CFL condition is computed via Equation (9a):

$$
\Delta t\_i^\* = \text{Cr} \min\_{j=1,2,3} (\frac{d\_{i,j}}{\sqrt{\mathbf{u}^2 + \mathbf{v}^2} + \sqrt{gh}}),
\tag{9a}
$$

where *u* and *v* are the maximum velocities in the *x* and *y* directions, respectively; and *h* is the maximum water depth.

Next, calculate the minimum time step Δ*tr* in the whole computational domain:

$$
\Delta t\_r = \min \left( \Delta t\_i^\* \right). \tag{69}
$$

(2) Set an initial temporal resolution level for each cell

In terms of the CFL condition, the initial time step Δ*t* ∗ *<sup>i</sup>* in cell *i* is proportional to the mesh scale *di*,*j*. Taking the minimum time step Δ*tr* over the whole computational domain as the reference standard, a binary logarithmic function is selected as the conversion tool. Each cell is classified based on its initial time step Δ*t* ∗ *i* , with the allocation of an initial time-step level *m*∗ *i* , reflecting the mesh scale and having integer values of 1, 2, 4, 8, 16, 32, 64, ... The initial *m*-level distribution map is shown in Figure 3a.

$$m\_i^\* = \text{int}(\log\_2(\frac{\Delta t\_i^\*}{\Delta t\_r})),\tag{10}$$

where the function of *int*() is to obtain the maximum integer that does not exceed the given number.

**Figure 3.** Example of the *m*-level assignment to each cell in the unstructured grid. (**a**) Distribution showing initial values assigned by Equation (10); (**b**) distribution after being modified by Equation (11a), with the reassignment of buffer zone values based on neighboring cells values.

(3) Modify the time-step level

The initial *m*-level distribution calculated by Equation (10) is usually scattered and broken (see Figure 3a), especially within the transition zone between coarse and refined grids. This distribution varies dramatically because of the large difference in mesh size. Appropriate measures need to be taken to smooth the transition and to reduce mass and momentum errors in these zones related to the calculation of conserved variables at interfaces between cells with different *m*-levels. The initial *m*-levels are modified via Equation (11a), which produces a buffer zone.

$$m\_i = \min(m\_{i\prime}^\*, m\_{i,1\prime}, m\_{i,2\prime}, m\_{i,3})\_\prime \tag{11a}$$

where *mi*,*j*(j = 1, 2, 3) is the time-step level of the three neighboring cells of cell *i*. The evaluation of the time-step level *mi* in cell *i* is based on values of the neighboring cells: *mi*,1, *mi*,2, *mi*,3. If the time-step level of the current cell is larger than that of the neighboring cells, then the *mi* value will be reduced by one; otherwise, it will be retained. Thus, the difference between the *m*-levels of the current and neighboring cells cannot exceed one. In this way, the time-step level *mi* (see Figure 3a) is modified via Equation (11a), and the new *m*-level distribution is as shown in Figure 3b.

The maximum of time-step level *mmax* is given by

$$m\_{\text{max}} = \max(m\_i),\tag{11b}$$

where *L* = *mmax* + 1 refers to the local time-step rank of the LTS algorithm in the computation (L = 1 represents the GTS solution). Building on the requirement of simulation accuracy, the maximum value of the different time-step level *mu* is selected according to the simulation accuracy requirements, and the restriction condition *mi* <sup>=</sup> min *mu*, *m*<sup>∗</sup> *i* , *mi*,1, *mi*,2, *mi*,3 is added to Equation (11a).

In this case, the time step adopted by each calculation cell Δ*t b <sup>i</sup>* is

$$
\Delta t\_i^b = 2^{m\_i} \Delta t\_r. \tag{11c}
$$

The maximum time step Δ*tmax* is

$$
\Delta t\_{\text{max}} = \mathcal{Z}^{\text{tr}\_{\text{max}}} \Delta t\_r \tag{11d}
$$

#### 2.2.2. Calculation of Element Variables at the Interface

Once the local time-step level of each computation cell is reconstructed, no further cells are updated with the same time step but are instead updated with the locally allowable time step at each

cell. The algorithm focuses on the convergence of the element variables at different time levels along the interface. In a maximum time step Δ*tmax*, the number of updates at each cell *Ni* is calculated using

$$N\_{\bar{l}} = 2^{m\_{\max} - m\_{\bar{i}}}.\tag{12}$$

The expression for Δ*U* in Equation (7) is then replaced by Equation (13), in which Δ*t b <sup>i</sup>* from Equation (11c) takes the place of Δ*tr* in each computing cell. The improved expression is as follows:

$$
\Delta l I = \mathbf{I} \frac{\Delta l\_i^b}{A} [(-\sum\_{j=1}^3 \mathbf{E}\_j \mathbf{n}\_j l\_j) + \sum\_{j=1}^3 \sum\_{k=1}^3 \left(\frac{1}{2} (1 - \text{sign}(\widetilde{\lambda}\_k)) \mathbf{\beta}\_k \widetilde{\mathbf{e}}\_k l\_j)^j + \Delta l\_i^b S\_f^a]. \tag{13}
$$

Then, Equation (13) is iterated *Ni* times, and the conserved variables of each cell are updated to the same time level *n*, from which the variables at the next time level *n* + 1 are calculated until the calculation time terminates.

The fluxes on the interface of the governing cell are estimated from the approximate Riemann solution of the Roe scheme. However, the variable values of all cells at the current time level need to be given to confirm the *UR* and *UL* values of the adjacent states on the left and right sides of the interface. In the initial GTS model, the unified global time step Δ*t* can meet these requirements. In the LTS algorithm, however, when the *m*-level of the current and neighboring cells differ within the maximum time step Δ*tmax*, the update times are different, and time splitting occurs at the interface. In the present study, the value of the missing element variable at the interface was obtained by linear interpolation of the value of the element variable in two adjacent time layers. This is shown in Figure 4, where the maximum time-step level *mmax* = 2. At the starting time of 0, we assumed that Δ*tr* = 0.01 s. First, the element variables of the coarse grid (*m* = 2) were computed only once, which consumes 0.04 s. Next, the values for the intermediate meshes (*m* = 1) were calculated twice, costing 0.02 s for each calculation. The local time step of the coarse mesh (*m* = 2) at the interface is, however, 0.04 s, which suggests that the variable values at 0.02 s are missing. In the present study, the variable values U(*h*, *u*, *v*) were estimated approximately using the average value of the former and the latter time step (1/2 of the variable value at 0 and 0.04 s). Finally, the refined grids (*m* = 0) were calculated at four intervals, with time steps of 0.01 s. At 0.01 and 0.03 s, the values of the variables are absent at the interface of the cells having *m* = 1. The above approximation method was again used. At 0.04 s, all cell variables were updated to the same time level in preparation for the next iteration.

**Figure 4.** Sketch showing the calculation of element variables at an interface: the variable U(*h*, *u*, *v*) refers to the water depth *h*, and the velocity components *u* and *v* in the *x* and *y* directions, respectively. Within a maximum time step of Δ*tmax* = 0.04 s, the coarse mesh areas (time-step level, *m* = 2) were updated only once. The value of U1(*h*1, *u*1, *v*1) at 0 s was updated to U<sup>∗</sup> 1 *h*∗ <sup>1</sup>, *u*<sup>∗</sup> <sup>1</sup>, *v*<sup>∗</sup> 1 after 0.04 s. The blank value at 0.02 s was estimated using the average of the values at the former and the latter time-step level, giving U- 1 *h*1+*h*<sup>∗</sup> 1 <sup>2</sup> , *u*1+*u*<sup>∗</sup> 1 <sup>2</sup> , *v*1+*v*<sup>∗</sup> 1 2 . Immediately following this, the intermediate grids (*m* = 1) were updated twice over this maximum time step. The values of variables at 0.01 and 0.03 s were similarly approximated. Finally, the refined mesh grids (*m* = 0) were updated four times over this interval. All variables were ultimately updated to 0.04 s, which allowed the next time step to be performed.

#### **3. Numerical Tests**

Numerical tests were carried out to verify the reliability and adaptability of the LTS strategy within a 2D hydrodynamic model, and to make a comparison between the effects of the novel LTS and the original GTS strategies on the simulation results. To quantify these differences, we used the mean-square error (MSE)

$$L\_s(f) = \sqrt{\frac{1}{N\_\varepsilon} \sum\_{i=1}^{N\_\varepsilon} \left(f\_{i, GTS} - f\_{i, L}\right)^2} \tag{14}$$

where *Ls*(*f*) refers to the MSE of variables *f* (such as water depth *h*, and velocity components *u* and *v* in the *x* and *y* directions, respectively), and the subscript L refers to the local time-step rank of the solutions.

The principle purpose of using the LTS algorithm instead of the GTS algorithm is to reduce the computing time. The relative time saving ratio *Tr* and speedup ratio *Sn* were used as evaluation indices of the computing efficiency for different L values.

The time saving ratio *Tr* is given by

$$T\_r = \frac{T\_{GTS} - T\_L}{T\_{GTS}} \times 100.\tag{15a}$$

The speedup ratio *Sn* is given by

$$S\_n = \frac{T\_{GTS}}{T\_L} \tag{15b}$$

where *TGTS* is the computation time of the GTS strategy, and *TL* is the computation time of the LTS strategy using different local time-step rank values (L).

#### *3.1. Anti-Symmetric Dam Break Case*

A classic verification experiment that considers the anti-symmetric square dam break case has been described by Fennema and Chaudhry [33]. The computational domain consists of a 200 × 200 m square region with a flat bottom (bed elevation of 0 m). In this problem, the dam is assumed to fail instantaneously at a certain moment, producing an anti-symmetric breach with a width of 75 m at the center of the calculation domain, which is 95 m from the *x*-axis boundary. The detailed geometry of this problem is shown in Figure 5.

**Figure 5.** Sketch showing the plane geometry of an anti-symmetric dam break, with the breach occurring in the center.

The computational domain was discretized on triangular meshes. Local mesh refinement technology was adopted near the break to accurately reflect the important flow regions while controlling the number of grids. The spatial resolution level varied from 6 to 1 m, with a refined area of 5%. A total of 13,338 computing cells were generated. The LTS algorithm was applied to shorten the computing time and improve the calculation efficiency. The *m*-level distribution of the LTS strategy is presented in Figure 6. The simulation duration was *Ts* = 160 s. At the initial moment, upstream water depth was set to 10.0 m, while downstream water depth was assumed to be 5.0 m. Manning's roughness coefficient was set to 0.03 [34].

Before making a quantitative analysis of the MSE generated by the LTS strategy, it is useful to verify the validity of the GTS scheme applied to the original model. The numerical solution results at t = 7.2 s were obtained using the GTS strategy; the resulting water surface profiles are shown in Figure 7a, contours of water depth in Figure 7b, and velocity vector distribution in Figure 7c. These results effectively simulated the dam break waves and are consistent with previous research results [34–36].

**Figure 6.** Details of the non-uniform unstructured meshes and time-step level (*m*-level) distribution for the anti-symmetric dam break case.

**Figure 7.** Results of the simulation using the global time stepping strategy at dam break (*t* = 7.2 s). (**a**) Three-dimensional visualization of the water surface; (**b**) contour plot of water depth; (**c**) velocity vector distribution.

Table 1 shows the MSE of water depth *h*, as well as velocity components *u* and *v* along the *x* and *y* directions, respectively, of the simulation results. As the value of the local time-step rank L increased, the MSE increased accordingly. Both precision and efficiency are important when performing hydrodynamic numerical simulations. By reducing the value of L, the calculation efficiency can be improved when levels of high accuracy are reached. This observation agrees well with the research results of Hu et al. [37]. In the present study, the error of the solution using the LTS algorithm was within the 10−<sup>3</sup> to 10−<sup>2</sup> range. At the beginning of the calculation, the total water volume was 3 <sup>×</sup> 105 m3. When the value of L was 2, 3, and 4, the relative error of the water volume was 0.0017%, 0.0019%, and 0.0045%, respectively, which indicates that the simulation meets the accuracy and conservation of water volume requirements.

Table 2 shows the computing time, time-saving ratio, and speedup ratio for different values of L. The simulation time of the dam break was 160 s. A time saving of 40% to 50% was achieved using the LTS algorithm in the numerical experiment of the anti-symmetric square dam break. For L = 4, a speedup ratio of 2.1 was achieved.

This clearly shows that the calculation efficiency of the numerical simulation of the anti-symmetric square dam break is improved using the LTS technique (Figure 8). As the value of L increased, the speedup ratio also increased, thus shortening the model operation.

**Figure 8.** Computational cost (s) and speedup ratios for different local time-step rank values (L) for the anti-symmetric dam break case. As the value of L increases, the computational cost decreases and the speedup ratio increases.



**Table 2.** Computing time T (s), time-saving ratio *Tr* (%), speedup ratio (*Sn*) of different local time-step rank values (L) for the anti-symmetric dam break case (GTS: global time stepping strategy).

*3.2. Non-Flat Bottom Dam Break Case*

In this test case, the length of the rectangular channel was 75 m, and its width was 30 m. Three obstacles were placed on the bottom of the rectangular channel. Two obstacles had a radius of 6.5 m and a height of 1 m, and one had a radius of 8.6 m and a height of 2 m. The elevation expression is defined as follows:

$$Z\_{\Phi}(\mathbf{x},y) = \begin{cases} 1 - [(\mathbf{x}-30)^2 + (y-8)^2]/42.25, \ (\mathbf{x}-30)^2 + (y-8)^2 \le 42.25\\ 1 - [(\mathbf{x}-30)^2 + (y-22)^2]/42.25, \ (\mathbf{x}-30)^2 + (y-22)^2 \le 42.25\\ 2 - [(\mathbf{x}-52)^2 + (y-15)^2]/36.89, \ (\mathbf{x}-52)^2 + (y-15)^2 \le 73.98\\ 0, \ \text{else} \end{cases}$$

.

Local mesh refinement was applied with a spatial resolution level varying from 0.5 to 3 m. The grid generation results are shown in Figure 9. The total number of grid cells for the region was 6848. In the calculation, a solid wall boundary condition was adopted. Manning's roughness coefficient was set to 0.018 [30]. The water depth for *x* ≤ 16 m in the computational domain was set to 2 m and was assumed to be 0 m elsewhere. The variation of the ground level is shown in Figure 10a. In the GTS scheme, the time step was 0.005 s, the initial total water volume was 959.841 m3, and the final water volume was 959.978 m3. The relative error was 0.01%, which is within the discrete error range. The model dealt well with wetting and drying fronts over the complex terrain, satisfied computational stability, and conservation of water volume requirements, and provided results that were consistent with previous experimental results [38,39].

**Figure 9.** Details of the non-uniform unstructured meshes, time-step level (*m*-level) distribution, and bathymetric point arrangement for the non-flat bottom dam break case.

**Figure 10.** Three-dimensional visualization of the inundation process at various time points following the dam break. (**a**) *t* = 0 s; (**b**) *t* = 10 s; (**c**) *t* = 30 s; (**d**) *t* = 100 s. (GTS: global time stepping strategy, while L is the time-step rank value in the local time stepping strategy).

The distribution of *m*-values obtained by applying the LTS strategy is shown in Figure 9. The maximum time-step level *mmax* is 2. The simulation of the dam break water flow was modeled for *Ts* = 100 s. The inundation related to the dam break is shown in Figure 10. To observe the influence of the LTS algorithm on the simulation results, four bathymetric measuring points were arranged (see Figure 9 for the location of the measuring points). Figure 11 provides a comparison of water depths generated by the LTS algorithm and the classical algorithm during the process of the dam break. A visual inspection indicates that the simulated water depths from the two algorithms are

essentially the same. Thus, the LTS technique improves the efficiency of calculation without any reduction in accuracy.

**Figure 11.** Comparison of water levels at different local time-step rank values (L) for the non-flat bottom dam break case (GTS: global time stepping strategy).

According to the statistical results for this simulation (shown in Table 3), adopting the LTS technology in the numerical experiment of the non-flat bottom dam break created a saving of up to 26% in computation time, and a speedup ratio of 1.3, without compromising the solution accuracy. The best result was obtained for L = 3, which reduced the computational cost and improved the calculation efficiency.


**Table 3.** Computing time T (s), time saving ratio *Tr* (%), speedup ratio (*Sn*) for different local time-step rank values (L) for the non-flat bottom dam break case (GTS: global time stepping strategy).

#### *3.3. Navigable Flow Simulation Case*

There are several noticeable perilous waterways between the Three Gorges Dam and the Gezhouba Dam. The river section is about 38 km long [40,41]. Bed elevations measured in October 2008 provided by the Three Gorges Navigation Authority were used as the initial bed topography. In the present study, the three segments with the highest risk to shipping (the Letianxi, Liantuo, and Shipai riverways), and the outlet of the Three Gorges Reservoir are adopted for local mesh refinement. The spatial

resolution level ranged from 10 to 100 m, and 36,297 cells and 19,011 nodes were generated. The results of the grid generation and refined areas are shown in Figure 12.

**Figure 12.** Computational cells and details of refined areas for the river reach between the Three Gorges and Gezhouba Dams.

The LTS algorithm was applied to the numerical simulation of the navigable flow of the waterway. The results are provided in Table 4. The maximum time-step level (*mmax* = 4) and the minimum time step (Δ*tr* = 0.03 s) were calculated for the whole computational domain. The simulation duration was set to *Ts* = 4 h. The absolute deviation in the water level between the calculated and measured values was 0.06 to 0.07 m, and the relative deviation in flow velocity was between −3.8% and 4.6%. There is no significant difference between the LTS and GTS algorithm simulation results for water level and flow velocity. The mass errors are non-negligible for a large computational area and a long simulation. This was addressed by adopting a smaller maximum time-step level *mmax*. The obvious acceleration effect that this provides in the real application ensures that accuracy requirements are met. It is clear that the LTS algorithm has provided obvious improvements in the execution time without compromising the solution accuracy. This is of practical value and importance in navigation.

**Table 4.**Simulation results for the navigable flow of the waterway between the Three Gorges and Gezhouba Dams.

GTS: global time stepping strategy; L: local time-step rank in the local time stepping strategy; T: execution time (*h*); *Tr*: time saving ratio (%); and *Sn*: speedup

 ratio.

#### *Water* **2020** , *12* ,

#### **4. Discussion**

#### *4.1. The Influence of the Proportion of Refined Mesh on the Acceleration E*ff*ect*

In the anti-symmetric dam break case, the 75-m breach was used as the core area in which to analyze the acceleration effect of the LTS algorithm. In this area, we set the refined area to different proportions of the total area (4 hm2). The spatial resolution level varied from 1 to 6 m, with a minimum resolution close to the breach of 1 m, surrounded by a smooth transition up to a maximum resolution of 6 m in the outer domain. The mesh generation results for the refined areas for various proportions in the anti-symmetric dam break case is shown in Figure 13. The numerical simulation was performed for refined areas of 5%, 10%, 25%, 50%, and 75%. The impact of increasing the local refinement area on the acceleration effect was analyzed for a time-step level of *m* = 0 in the refinement area, *m* = 2 in the coarse grid area, and *m* = 1 in the buffer zone.

**Figure 13.** Locally refined areas of various proportions (%) in the anti-symmetric dam break case: the refined area of (**a**) 5%; (**b**) 10%; (**c**) 25%; (**d**) 50%, (**e**) 75%.

The duration of the dam break *Ts* was kept at 160 s. Table 5 shows the number of generated grids, computational cost, speedup ratio, and time-saving ratio for different refined area proportions. As the refined area was increased, the total number of meshes generated by partitioning increased. The grid number increased from 13,338 at 0.2 hm<sup>2</sup> to 73,306 at 3 hm2, which represents an increase by a factor of 5.5. The operating time of the model rose from 608 to 3394 s, which represents an increase by a factor of 5.6. The number of grids is a key factor determining the overhead of computing. There is a direct positive correlation between the number of grids and computing overheads.


**Table 5.** Simulation results for different refined area proportions in the dam break simulation.

NB: GTS: global time stepping strategy; L: local time-step rank in the local time stepping strategy; T: execution time (s), *Tr*: time-saving ratio, *Sn*: speedup ratio.

The speedup ratio was strongly affected by the ratio of the locally refined area to the total area. As the proportion of the refined area increased, the speedup ratio obtained by the LTS algorithm diminished, reducing the acceleration effect. When L = 3, the speedup ratio decreased from 1.39 times to 1.01 times; when the refined area was increased from 5% to 75%, there was scarcely any improvement in computational efficiency at high proportions of refined area. These patterns occurred because a high proportion of refined area brings about an unreasonable number of refined grids (with a local time-step level, *m* = 0), which reduces the number of coarse grids (*m* = 1, 2) that help optimize the operation. This suggests that little improvement is gained by implementing the LTS algorithm in cases that have a high proportion of refined grids. Conversely, when refined grids account for a small proportion of the total number of grids, the LTS algorithm can markedly shorten the computing time and achieve an excellent acceleration effect. Furthermore, the results in Figure 14 suggest that L has an influence on the calculation efficiency. The acceleration ratio increased as L increased.

**Figure 14.** Speedup ratio for implementing the LTS algorithm for various refined area proportions in the dam break model for different local time-step rank values (L). As the refined area increased, the speedup ratio decreased. The acceleration effect increased with increasing L.

#### *4.2. The Impact of the Di*ff*erent Scale of the Mesh on Acceleration E*ff*ect*

In the case study of the anti-symmetric dam break, the test was carried out using a refinement of 5% of the breach area, as shown in Figure 13a. A grid discretization method was employed to analyze the influence of using different grid scales on the acceleration effect while using the LTS algorithm. The refined grid had a minimum spatial step of 1 m in the burst section, while the coarse grids had spatial steps of 2, 3, 4, 5, and 6 m. The computational results are shown in Table 6 and Figure 15. These results indicate that the speedup ratio improves with the increasing differences in grid resolution. The CFL condition of Equation (9a) shows that the time step of the mesh is positively correlated with the mesh scale difference. As the difference between the maximum and minimum scales of the grids increases, the maximum time-step level *mmax* becomes higher, and a larger L can be adopted. In this test, once the spatial step of the coarse grid was 4, 5, and 6 m, it was possible to set the L value to 3, while the other two conditions were only subdivided into two grades. This explains why greater grid diversity leads to a higher speedup ratio.


**Table 6.** Simulation results for different spatial resolutions in the dam break simulation.

(GTS: global time stepping strategy, LTS: local time stepping strategy, T: execution time (s), *Tr*: time-saving ratio, and *Sn*: speedup ratio).

**Figure 15.** Graph showing the relationship between speedup ratio and grid scale difference.

#### **5. Conclusions**

Improvements in the computational efficiency of the 2D hydrodynamic model have a significant bearing on its engineering applications. In the present study, an improved 2D shallow-water dynamic model applying an LTS algorithm was established using an FVM on a triangular mesh. Results from simulation and analysis of three canonical test cases show that the LTS scheme has a satisfactory ability for adapting real complex applications and long simulations while meeting the required accuracy. The following conclusions were drawn:


**Author Contributions:** Conceptualization, S.Z. and X.Y.; Methodology, S.Z. and X.Y.; Software, X.Y., W.A. and W.L.; Validation, S.Z., X.Y., W.A. and W.L.; Formal Analysis, S.Z.; Investigation, X.Y. and W.A.; Resources, S.Z.; Data Curation, X.Y. and W.A.; Writing—Original Draft Preparation, X.Y.; Writing—Review & Editing, S.Z.; Supervision, S.Z. All authors have read and agreed to the published version of the manuscript.

**Funding:** This study was supported by the National Natural Science Foundation of China (51979105, 51722901).

**Acknowledgments:** The authors thank Trudi Semeniuk and Paul Seward from Liwen Bianji, Edanz Editing China (www.liwenbianji.cn/ac), for editing the English text of drafts of this manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


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

## *Article* **E**ff**ect of Open Boundary Conditions and Bottom Roughness on Tidal Modeling around the West Coast of Korea**

#### **Van Thinh Nguyen \* and Minjae Lee**

Department of Civil and Environmental Engineering, Seoul National University, 1 Gwanak-ro, Gwanak-gu, Seoul 151-744, Korea; mjlee14@snu.ac.kr

**\*** Correspondence: vnguyen@snu.ac.kr; Tel.: +82-2-880-7355

Received: 9 May 2020; Accepted: 11 June 2020; Published: 15 June 2020

**Abstract:** The aim of this study was to investigate the effect of open boundary conditions and bottom roughness on the tidal elevations around the West Coast of Korea (WCK) using an open-source computational fluids dynamics tool, the TELEMAC model. To obtain a detailed tidal forcing at open boundaries, three well-known assimilated tidal models—the Finite Element Solution (FES2014), the Oregon State University TOPEX/Poseidon Global Inverse Solution Tidal Model (TPXO9.1) and the National Astronomical Observatory of Japan (NAO.99Jb)—have been applied to interpolate the offshore tidal boundary conditions. A number of numerical simulations have been performed for different offshore open boundary conditions, as well as for various uniform and non-uniform bottom roughness coefficients. The numerical results were calibrated against observations to determine the best fit roughness values for different sub-regions within WCK. In order to find out the dependence of the tidal elevation around the WCK on the variations of open boundary forcing, a sensitivity analysis of coastal tide elevation was carried out. Consequently, it showed that the tidal elevation around the WCK was strongly affected by local characteristics, rather than by the offshore open boundary conditions. Eventually, the numerical results can provide better quantitative and qualitative tidal information around the WCK than the data obtained from assimilated tidal models.

**Keywords:** west coast of Korea (WCK); tidal elevation; open boundary condition; numerical modeling; TELEMAC-2D

#### **1. Introduction**

The Yellow Sea (YS) is located between the mainland of China and the Korean Peninsula, and its depth is about 44 m on average, with a maximum depth of 152 m. The sea bottom slope gradually rises toward the East China coast and more abruptly toward the Korean Peninsula, and the depth gradually increases from north to south. Southern YS is bounded by the East China Sea (ECS) along a line running northeastward from the mouth of the Changjiang River to the southwestern tip of Korea. ECS covers the area originating around the Taiwan Strait at about 25◦ N, extending northeastward to the Kyushu Island, Japan, and to the Korea Strait, and bounded by the Ryukyu Island chains on the southeastern open ocean (Figure 1). On the broad and shallow shelf under the Yellow and the East China Sea (YECS), the flow is dominated by strong semidiurnal M2 tidal currents with the superimposition of semidiurnal S2, diurnal K1 and O1 currents [1–4]. Approaching the WCK, the tidal amplitudes vary from 4 to 8 m, and the speed of the tidal current may increase to more than 1.56 m/s near the coasts.

**Figure 1.** Geophysical configuration of the Yellow Sea and the East China Sea.

Along with strong tidal currents, the WCK is surrounded by lots of islands including extremely irregular and indented coastlines, the tidal propagation and shoaling processes on the shore in this broad flat region are very complicated and sensitive to the variation of tidal elevation.

The study on the tidal circulation in this region has been carried out by various authors, such as Le Fevre et al. [5], Kim [6], Naimie et al. [7], and Suh [8] using 2D depth-averaged finite element models on unstructured grids, or Choi [2], Tang [9], Kantha et al. [10], Kang et al. [11] and Lee et al. [12] using 3-D finite difference models on structured grids. However, there are still very few authors using data obtained from different assimilated tidal models to set up the offshore boundary conditions. Most authors have applied a constant boundary forcing condition (BFC) or tidal charts, such as Ogura's chart [1,2,13,14]) and Nishida's chart [3,15,16] or applied a very coarse observation data combined with partly constant boundary forcing conditions [11,17]. Only Le Fevre et al. [5] and Suh [8] have applied tidal charts from FES 94.1 and FES 2004 models for OBC. Using FVCOM, we recently applied a tidal chart from NAO.99Jb. In principle, the tidal information around the WCK can be obtained from assimilated global and regional tidal models [18].

In principle, the tidal information around the WCK can be obtained from assimilated global and regional tidal models. The assimilation model, which is a combination of a dynamical model and observed data, was actively advanced in the research group of the modern global tide model, and pioneered by [19]. By the assimilation model, the problem of the empirical method with regard to the resolution can be compensated by a hydrodynamic model, and any inadequate potential of the hydrodynamic model can be recovered by observed data [20]. Recently, using the unprecedented sea surface height data of satellite altimeter TOPEX/POSEIDON (T/P), several research groups have developed highly accurate assimilation models especially in an open ocean ([5,20–24]). T/P was launched on 10 August 1992, which measures sea level relative to the ellipsoid with an overall accuracy of approximately 5 cm [25]. Due to the high accuracy of T/P altimeter data for those assimilation

models, it can be used to extract relatively accurate tidal boundary conditions around the marginal sea within YECS.

Among the global assimilation models, FES2014 (https://www.aviso.altimetry.fr/en/data/products/ auxiliary-products/global-tide-fes.html) and TPXO9.1 (http://volkov.oce.orst.edu/tides/tpxo9\_atlas. html) atlases are open to the public and can be freely downloaded. The FES2014 atlas was produced by Legos and CLS Space Oceanography Division, and distributed by Aviso, with support from CNES (https://datastore.cls.fr/catalogues/fes2014-tide-model/). It is distributed as a tidal chart on grid resolution of 1/16◦ based on the dynamical finite element tide model CEFMO, and data assimilation code CADOR. The finite element grid of the model CEFMO has about 2,900,000 computational nodes and covered 34 tidal constituents. In addition to above two global tide models, the regional tide model NAO.99Jb (https://www.miz.nao.ac.jp/staffs/nao99/index\_En.html) was developed by Matsumoto et al. (2000) for the region around Japan with a resolution of 1/12◦, covering from 110 to 165◦ E and from 20–65◦ N so that YECS was entirely included in this region. This model was nested with the global model NAO.99b to provide open boundary conditions to the regional model and the assimilation was performed with both T/P data and 219 coastal tidal gauges obtained around Japan and Korea. Most T/P-derived altimetric tides can be inaccurate in shallow coastal oceans because the resolution enforced on data analysis by the spatial resolution of T/P is simply inadequate near ocean margins.

Therefore, we first validated the tidal information of tidal constituents *K*1, *O*1, *M*<sup>2</sup> and *S*<sup>2</sup> provided by three assimilated tidal models; FES2014, NAO99Jb and TPXO9.1 against the observation at 21 tidal gauges along WCK. Figure 2 shows the comparison of amplitude (left) and phase (right) obtained from the assimilated tidal models with observation data obtained from gauging stations along the WCK. It shows that while the phases are better agreement with observation, the amplitudes obtained from the assimilated tidal models still poorly agreed with the observed data. More detailed quantitative comparisons between the amplitude and phase obtained from the assimilated tidal models and observation data are shown in Table A1 (in Appendix A). Wherein it shows NAO99Jb provides a better agreement of the amplitude for the constituents *M*<sup>2</sup> and *S*2, while FES2014 provides a little better fit of the amplitude for *K*<sup>1</sup> and *O*1. It also shows the constituent M2 has the largest amplitude in comparison with other tidal constituents; about 5–8 times, 2.5–2.8 times, and 6.6–11.8 times larger than K1, S2 and O1 amplitudes, respectively.

In addition, the resolutions (1/6◦, 1/12◦ and 1/16◦) of three tidal models are not fine enough along the coastlines, so they cannot provide enough detailed quantitative and qualitative tidal information along the WCK. Therefore, in order to obtain more detail on tidal information along the WCK, an accuracy numerical simulation with higher grid resolutions is needed.

This study is to introduce several numerical investigations using higher grid resolutions and applying different open boundary conditions extracted from three well-known assimilated tidal models FES2014, TPXO9.1 and NAO.99Jb then calibrated with different uniform and non-uniform bottom roughness values for different sub-regions within WCK to investigate the effect of open boundary condition and bottom roughness on the tidal elevations from major constituents (M2, K1, S2 and O1) around the West Coast of Korea.

Based on the form number (*F*), the ratio of the amplitudes of *K*1, *O*1, *M*2, *and S*<sup>2</sup> suggested by [26] is shown in Equation (1).

$$F = \frac{K\_1 + O\_1}{M\_2 + S\_2} \tag{1}$$

whereby the mixed tide is classified following the value of the form number *F*; if the value of *F* less than 0.25 or larger than 1.25, then it is classified as semidiurnal or diurnal, respectively. Otherwise, if the value of *F* is located between 0.25 and 1.25, it is considered mixed.

**Figure 2.** Validation of amplitudes and phases obtained from assimilated tidal models (FES2014, NAO99Jb, NAO99 and TPXO9) against observed data at tidal gauges; (**a**): M2's amplitude, (**b**): M2's phase, (**c**): K1's amplitude, (**d**): K1's phase, (**e**): S2's amplitude S2, (**f**): S2's phase, (**g**): O1's amplitude, (**h**): O1's phase.

Figure 3 shows the values of the form number along the WCK are below 0.25 according to the data obtained from observation at 21 tidal gauges, as well as from three assimilated tidal models, FES2014, TPXO9.1 and NAO.99Jb; therefore, the mixed tide prevails semidiurnal in the WCK.

**Figure 3.** The value of the form number (*F*) in the WCK. Grey color in the background shows the WCK, from left to right corresponding to north to south.

#### **2. Numerical Modeling**

As mentioned above, since the study region was quite shallow, the numerical simulation was based on an open-source software, the TELEMAC-2D (http://www.opentelemac.org) model. It is well-known software with an abundant history of development and application over many years in fluvial and maritime hydraulics [27,28]. The TELEMAC-2D is capable of capturing the complicated bathymetry and coastlines of the WCK surrounded by lots of islands including extremely irregular and indented coastlines by using unstructured grids.

#### *2.1. Basic Equations*

The TELEMAC-2D has solved the depth-averaged Navier-Stokes equations with two-equation turbulence closure models (k-ε) using the finite element method, as follows:

The continuity equation

$$\frac{\partial h}{\partial t} + \frac{\partial (h \, \mathrm{U})}{\partial \mathbf{x}} + \frac{\partial (h \, \mathrm{V})}{\partial y} = 0 \tag{2}$$

and momentum equations

$$\frac{\partial(h\boldsymbol{l}\boldsymbol{l}\boldsymbol{l})}{\partial\boldsymbol{t}} + \frac{\partial(h\boldsymbol{l}\boldsymbol{l}\boldsymbol{l}\boldsymbol{l}\boldsymbol{l})}{\partial\boldsymbol{x}} + \frac{\partial(h\boldsymbol{l}\boldsymbol{l}\boldsymbol{l}\boldsymbol{V})}{\partial\boldsymbol{y}} = -\boldsymbol{h}\cdot\boldsymbol{g}\frac{\partial\boldsymbol{Z}}{\partial\boldsymbol{x}} + \boldsymbol{h}\cdot\boldsymbol{F}\_{\boldsymbol{x}} + \mathrm{div}\Big{(}\boldsymbol{h}\cdot\boldsymbol{\nu}\_{\boldsymbol{t}}\cdot\overrightarrow{\boldsymbol{\nabla}}(\boldsymbol{l}\boldsymbol{l})\Big{)}\tag{3}$$

$$\frac{\partial(hV)}{\partial t} + \frac{\partial(hLV)}{\partial \mathbf{x}} + \frac{\partial(hVV)}{\partial y} = -h \cdot \mathbf{g} \frac{\partial \mathbf{Z}}{\partial y} + h \cdot \mathbf{F}\_y + \text{div} \left( h \cdot \nu\_{\boldsymbol{\varepsilon}^\*} \overrightarrow{\nabla}(V) \right) \tag{4}$$

$$\frac{\partial \mathbf{k}}{\partial t} + \mathcal{U} \frac{\partial \mathbf{k}}{\partial \mathbf{x}} + \mathcal{V} \frac{\partial \mathbf{k}}{\partial y} = \frac{1}{h} \operatorname{div} \Big( h \cdot \frac{\nu\_t}{\sigma\_k} \, \overrightarrow{\nabla}(\mathbf{k}) \Big) + P - \varepsilon + P\_{kv} \tag{5}$$

$$\frac{\partial \varepsilon}{\partial t} + \mathcal{U} \frac{\partial \varepsilon}{\partial \mathbf{x}} + \mathcal{V} \frac{\partial \varepsilon}{\partial y} = \frac{1}{h} \operatorname{div} \Big( h \cdot \frac{\nu\_t}{\sigma\_k} \, \overrightarrow{\nabla} (\varepsilon) \Big) + \frac{\varepsilon}{k} (c\_{1\varepsilon} P - c\_{2\varepsilon} \varepsilon) + P\_{\varepsilon \upsilon} \tag{6}$$

where *h* is the water depth; *U* and *V* are the depth averaged velocity components;

*U* = <sup>1</sup> *h Z Zf udz* and *V* = <sup>1</sup> *h Z Zf vdz*

*Z* and *Zf* are the free surface and bottom elevations, respectively;

Fx and Fy are body forces in x and y directions including Coriolis (*F<sup>c</sup> <sup>x</sup>*, *F<sup>c</sup> <sup>y</sup>*), bottom friction *<sup>F</sup><sup>f</sup> <sup>x</sup>* , *<sup>F</sup><sup>f</sup> y*), and surface wind *Fw <sup>x</sup>* , *Fw <sup>y</sup>* ) forces determined as *Fx* = *F<sup>c</sup> <sup>x</sup>* <sup>+</sup> *<sup>F</sup><sup>f</sup> <sup>x</sup>* + *F<sup>w</sup> <sup>x</sup>* and *Fy* = *Fc <sup>y</sup>* <sup>+</sup> *<sup>F</sup><sup>f</sup> <sup>y</sup>* + *F<sup>w</sup> y*

Coriolis force; *F<sup>c</sup> <sup>x</sup>* = 2ω sin(λ)*v* and *F<sup>c</sup> <sup>y</sup>* = <sup>−</sup>2<sup>ω</sup> sin(λ)*u*; <sup>ω</sup> = 7.292 <sup>×</sup> 10−<sup>5</sup> is angular and <sup>λ</sup> is latitude;

Bottom friction force; *F<sup>f</sup> <sup>x</sup>* = <sup>−</sup> <sup>1</sup> <sup>2</sup>*hCfU* <sup>√</sup> *U*<sup>2</sup> + *V*<sup>2</sup> and *F<sup>f</sup> <sup>y</sup>* = <sup>−</sup> <sup>1</sup> <sup>2</sup>*hCf <sup>V</sup>* <sup>√</sup> *U*<sup>2</sup> + *V*2, where *Cf* is the bottom friction coefficient, *Cf* <sup>=</sup> *gn*<sup>2</sup> *<sup>h</sup>*1/3 , n is Manning's coefficient;

Wind forcing on the free surface is neglected; i.e., (*Fw <sup>x</sup>* , *Fw <sup>y</sup>* ) = 0.

And ν*<sup>e</sup>* is the effective viscosity, the summation of the molecular viscosity ν and the turbulent viscosity ν*<sup>t</sup>* ν*<sup>e</sup>* = ν + ν*t*; where

$$\nu\_t = c\_\mu \frac{k^2}{\varepsilon}. \tag{7}$$

The depth averaged kinetic energy *k* and its dissipation rate ε are:

$$k = \frac{1}{h} \int\_{Z\_f}^{Z} \frac{1}{2} \overline{u\_i' u\_i'} dz \text{ and } \varepsilon = \frac{1}{h} \int\_{Z\_f}^{Z} \nu \frac{\overline{\partial u\_i'} \, \overline{\partial u\_j'}}{\partial \mathbf{x}\_j} dz$$

where *u*- *<sup>i</sup>* is the fluctuating velocity and the overbar represents an average over time; i.e., *ui* = *Ui* + *u*- *i* ; *Ui* (=*U* or *V*) is an average over time of velocity components; and *P* is the production term:

$$P = \nu\_l \left(\frac{\partial \mathcal{U}\_i}{\partial \mathbf{x}\_j} + \frac{\partial \mathcal{U}\_j}{\partial \mathbf{x}\_i}\right) \frac{\partial \mathcal{U}\_i}{\partial \mathbf{x}\_j} \tag{8}$$

Pkv and Pε<sup>v</sup> are vertical shear terms: *Pkv* = *ck u*3 ∗ *<sup>h</sup>* ; *P*ε*<sup>v</sup>* = *c*<sup>ε</sup> *u*4 ∗ *<sup>h</sup>*<sup>2</sup> where *ck* <sup>=</sup> <sup>√</sup><sup>1</sup> *c f* and *<sup>c</sup>*<sup>ε</sup> <sup>=</sup> 3.6*c*2<sup>ε</sup> √*c*<sup>μ</sup> *c f* 3/4 *Cf*

*u*<sup>∗</sup> is the friction velocity: *u*<sup>∗</sup> = <sup>2</sup> (*U*<sup>2</sup> + *V*2)

The empirical constants in Equations (5) and (6) *C*<sup>μ</sup> = 0.09, *C*1<sup>ε</sup> = 1.44, *C*2<sup>ε</sup> = 1.92, σ*<sup>k</sup>* = 1.0 *and* σε = 1.3 of the *k-*ε model are obtained from classical test cases.

#### *2.2. Boundary Conditions*

#### 2.2.1. At the Solid Boundary

On the solid boundary, such as coastline and bottom, the no-slip condition is applied; i.e., velocity → *U* = (*U*, *V*) = (0, 0).

In this region, particularly near the shallow coastline, the effect of bottom friction becomes prominent, the bottom friction is calculated by:

$$
\tau\_b = \begin{pmatrix} \tau\_{bx'} & \tau\_{by} \end{pmatrix} = -\frac{1}{2} \rho \mathbf{C}\_f \sqrt{\mathcal{U}^2 + V^2} (\mathcal{U}, V) \tag{9}
$$

The value of the coefficient *Cf* is obtained from the calibration procedure, which is shown in Section 3.1 below.

The boundary condition for turbulence kinetic energy and its dissipation are defined as

$$k\_{\delta} = \frac{u\_{\ast}^{2}}{\sqrt{\mathbb{C}\_{\mu}}} + \frac{\mathbb{C}\_{2\varepsilon}}{\mathbb{C}\_{\varepsilon}\mathbb{C}\_{f}}u\_{\ast}^{2} \text{ and } \varepsilon\_{\delta} = \frac{u\_{\ast}^{3}}{\kappa\delta} + \frac{1}{\sqrt{\mathbb{C}\_{f}}}\frac{u\_{\ast}^{3}}{h} \tag{10}$$

where δ is a distance from the wall, δ = 0.1 of a local mesh size; *u*<sup>∗</sup> is friction velocity of the wall.

#### 2.2.2. At the Free Surface

We assume the wind force on the free surface is neglected as mentioned above, i.e., *Fw <sup>x</sup>* = *F<sup>w</sup> <sup>y</sup>* = 0.

#### 2.2.3. At the Open Boundary

Tidal harmonic parameters obtained from each assimilated tidal models FES2014, NAO99Jb and TPXO9.1 were combined into the open boundary conditions, as follows: *H* = *<sup>i</sup> Hi*

$$H\_i(M, \mathbf{t}) = A\_{F\_i}(M) \cos \left(2\pi \frac{\mathbf{t}}{T} - q\nu\_{F\_i}(M)\right) \tag{11}$$

where *Hi* is the water elevation, *AFi* is amplitude; *T* is period and ϕ*Fi* is phase of each tidal constituent.

In this study, four major tidal constituents M2, K1, S2 and O1, were considered. The harmonic constants of tidal constituents for open boundary conditions were obtained from three assimilated tidal models; FES2014, TPXO9.1 and NAO.99Jb.

Turbulence kinetic energy and its dissipation at the open boundary were set as:

$$k = \text{C2}\_{\text{2}} \frac{P\_{kv}^2}{P\_{uv}} \text{ and } \varepsilon = P\_{kv}\dots$$

#### *2.3. Study Region*

The study area covered the entire YS and most of the ECS in spherical coordinate. The bathymetry data of 30 arc-seconds resolution was downloaded from the General Bathymetric Chart of the Ocean (GEBCO2014: https://www.ngdc.noaa.gov/mgg/ibcm/ibcmdvc.html) as shown in Figure 4. The shorelines with an accuracy of about 250 m were downloaded from the NOAA shoreline website (https://shoreline. noaa.gov). This region included the flat and broad continental shelf of YECS, bounded by Taiwan Strait (A), a line of shelf-break determined by 200 m isobaths (B), and Korea Strait (C), as shown in Figure 5. The open boundary was expanded to the shelf edge of YECS which was advantageous to use offshore co-oscillating tidal conditions as a boundary forcing in order to minimize disturbances by the geographical configuration and nonlinear tidal response near-shore. Open boundaries adjacent to the offshore also enabled the numerical results around WCK to minimize the influence of unexpected numerical instabilities at open boundaries.

**Figure 4.** Bathymetry of the YECS Region.

**Figure 5.** Horizontal unstructured grid of the YECS (**a**) and an enlargement of WCK (**b**).

In this study, the numerical simulation applied a hybrid horizontal resolution of an unstructured grid using a free pre-processing Blue Kenue software tool (https://nrc.canada.ca/en/research-development/ products-services/software-applications/blue-kenuetm-software-tool-hydraulic-modellers). Herein the unstructured triangular grid had a horizontal resolution varying from 0.2 to 1 km in the WCK, 3–5 km in the East coast of China, and 9–12 km in the interior of YS and ECS. After running a number of simulations to check for grid independence, the final grid of 203,927 grid cells was used for the simulation in this study. Bathymetry was interpolated to each nodal point of the horizontal grid using Blue Kenue software as well (Figure 5).

A total of 21 tidal gauges were collected around the WCK to provide the observed tidal elevation data of semidiurnal M2 and S2 and diurnal K1 and O1 constituents in YECS. As shown in Figure 6, the WCK can be divided into three sub-regions based on the common characteristics of their bathymetry and coastlines—the first region was Kyunggi Bay near Incheon (I), the second region was around the estuary of the Geum River near Gunsan (G), and the third region was around the south-western tip near Mokpo (M). The numbers of tidal gauges located in three regions Mokpo, Gunsan and Incheon were 2, 8 and 11, respectively. The most recent observation data in 2017 downloaded from Korea Hydrographic and Oceanographic Agency (http://www.khoa.go.kr/koofs/kor/observation/obs\_real.do) have been selected for the calibration and validation of numerical models.

**Figure 6.** Tidal gauge locations in the WCK selected for the evaluation of the modeling; three sub-regions: Mokpo (M), Gunsan (G) and Incheon (I) Regions.

#### **3. Numerical Results and Discussion**

The numerical simulations were carried out on our high-performance computing cluster (http: //cfdlabsnu.com) for the real-time simulation of 60 days with time step of 20 s. The numerical results obtained from the last 30 days simulation were used for the following analyses.

#### *3.1. Response of the Tide around WCK to the Bottom Roughness*

#### 3.1.1. Applying Uniform Roughness Coefficient

Previous tidal calculations in the YECS region were mostly using the quadratic bottom friction law, in which the typical value of roughness coefficient *Cf* was chosen as a uniform value of 0.0025, as shown in [2,16,29,30]. We again have verified the numerical model with various uniform roughness coefficients that ranged from 0.0015 to 0.03 for different offshore open boundary conditions (OBCs) interpolated from three assimilated tidal models (FES2014, TPXO9.1 and NAO.99Jb).

The WCK was divided into three different sub-regions—Mokpo (M), Gunsan (G), and Incheon (I)—as shown in Figure 7, for applying different bottom friction coefficients. The borderline of this sub-regions was set based on a water depth contour of around 50 m. Different values of bottom friction coefficient were applied with regards to the geographic characteristics.

**Figure 7.** Sub-dividing regions in WCK based on roughness coefficients.

As shown in Tables A2–A4 (in Appendix A), based on the evaluation of root mean square deviation (RMSE) of tidal amplitudes between the observation and the numerical results obtained from applying various uniform roughness coefficients, the results show that the numerical results were very distinct from each other. For the tidal constituent M2, the bottom roughness coefficient *Cf* = 0.0023 provides the best result obtained from FES2014 and TPXO9.1, while the numerical results obtained from NAO99Jb were not significantly different between *Cf* = 0.0023 and 0.0025. For the tidal constituent K1, the bottom roughness coefficient *Cf* = 0.002 provides the best result obtained from NAO99Jb, while the numerical results obtained from FES2014 were not significantly different between *Cf* = 0.0015 and 0.002; the numerical results obtained from TPXO9.1 were not significantly different between *Cf* = 0.002, 0.0023 and 0.0025. For the tidal constituent S2, the bottom roughness coefficient *Cf* = 0.0015 provides the best result obtained from all three assimilated tidal models. For the tidal constituent O1, the bottom roughness coefficient *Cf* = 0.0015 provides the best result obtained from FES2014 and TPXO9.1, while the numerical results obtained from NAO99Jb were not significantly different between *Cf* = 0.0015 and 0.002. In addition, based on the RMSE values in Tables A2–A4, the results show that for M2 the most appropriate roughness values for Mokpo, Gunsan and Incheon regions were 0.0025, 0.003 and 0.0023, respectively. Whereas, for S2, the most appropriate roughness values for Mokpo, Gunsan and Incheon regions were 0.002, 0.0023 and 0.0015, respectively. For O1, the most appropriate roughness value for Mokpo, Gunsan and Incheon regions was 0.0015 and for K1 the most appropriate roughness values with FES2014 and NAO99Jb for Mokpo, Gunsan and Incheon regions were 0.0015, 0.002 and 0.002, respectively. However, with TPXO9.1 the values for Mokpo, Gunsan and Incheon regions were 0.002, 0.0023 and 0.002, respectively.

#### 3.1.2. Applying Non-Uniform Roughness Coefficient

As shown above, applying a uniform bottom roughness coefficient for the entire WCK was inappropriate. According to the bathymetry data, the bottom slope was relatively milder in the Incheon region than in Mokpo and Gunsan regions. In addition, there were lots of islands around the Mokpo region rather than the Incheon region. We can expect that more energy dissipation occurred in the Mokpo and Gunsan regions than in the Incheon region. The larger values of the bottom friction coefficient were applied to Mokpo and Gunsan regions than the Incheon region. In addition, the roughness coefficient values for the Gunsan region were set to be the largest because its bed form was more variable than the Mokpo region.

Table A5 shows the list of simulation cases applying non-uniform bottom roughness coefficients to different sub-regions Mokpo, Gunsan, Incheon and other regions within WCK. As shown in Tables A6–A9, once the non-uniform bottom roughness coefficients were applied, the mean RMSE values for M2 from the varied ranged of 14–36 cm down to a value of about 11 cm with OBCs obtained from FES2014 and TPXO9.1. These values varied from 16–41 cm down to 12–13 cm with OBCs obtained from NAO99Jb while, there was no significant improvement for K1, S2 and O1. In the WCK, M2 was the dominant constituent, followed by S2 and K1, as mentioned by Teague et al. (1998). In addition, the semidiurnal tides M2 and S2, as well as the diurnal tides K1 and O1, had similar tendencies responding to the local effects and OBCs. Therefore, M2 was presented for the semidiurnal tide, and K1 presented for diurnal tides in most following numerical results. The comparisons of M2 and K1 amplitudes between the observations and simulation results with uniform and non-uniform bottom roughness coefficients are shown in Figures 8 and 9. It clearly shows the simulation results obtained from non-uniform bottom roughness coefficients have significantly improved the results for M2 and K1.

**Figure 8.** The numerical results of tidal amplitudes M2 and K1 obtained from different uniform bottom roughness; (**a**), (**c**) and (**e**): M2's amplitude; (**b**), (**d**) and (**f**): K1's amplitude.

**Figure 9.** The numerical results of tidal amplitudes M2 and K1 obtained from non-uniform bottom roughness; (**a**), (**c**) and (**e**): M2's amplitude; (**b**), (**d**) and (**f**): K1's amplitude.

The bottom friction coefficient is a physical parameter characteristic of bottom materials and bathymetries. Therefore, it cannot vary following each tidal constituent. From Tables A5–A9, an acceptable selection from all test cases, which can deliver reasonable results for all tidal constituents (M2, K1, S2 and O1), was the case name "FES/NAO/TPX 2-4" whereby the bottom friction coefficient *Cf* took a value of 0.0025 for the Mokpo region, 0.0035 for the Gunsan region, and 0.002 for Incheon and other regions. These values are in the range suggested by [31,32] applying also for four tidal constituents (M2, K1, S2 and O1).

#### *3.2. Response of the Tide around the WCK to the Open Boundary*

As mentioned in Section 2.2 above, the modeling region has been bounded by three different open boundaries; A, B and C, as shown in Figure 3. Averaged tidal amplitudes at each boundary A, B and C obtained from three tidal models are shown in Table A10, in which the nodal boundary amplitudes are averaged with weighting by the length of each segment of boundary cells. The tidal amplitudes along the boundary A were significantly larger than the values along the boundaries B and C. Particularly, the tidal amplitude of M2 along the boundary A were larger than the values along the boundaries B and C; about 3.3 times and 9.3 times, respectively. The tidal amplitudes at the boundary C were the

smallest because when the tides propagated from the deep region of the East Sea (its depth was about 2000 m) through boundary C via the shallow region at the Korea Strait (its depth was less than 80 m), the tidal amplitudes were significantly damped due to shoaling process on this shallow shelf region.

Table A11 shows the difference of the interpolated tidal amplitudes at three open boundaries (A, B and C) obtained from two global tidal models (FES2014 and TPXO9.1) in comparison with those obtained from the regional tidal model (NAO.99Jb). It shows that the amplitudes of M2, K1, S2 and O1 obtained from three tidal models were not substantially different at the open boundaries A, B and C.

Nevertheless, it poses a question whether the open boundary condition at A, B or C will play a substantial role on the tide around the WCK. To identify such influence, we have evaluated the response of coastal tide to open boundary forcing in the WCK, using the sensitivity analysis in the following sections.

#### 3.2.1. Sensitivity Analysis

Response of the Tide around the WCK to Individual Boundary Forcing

The first numerical test was performed to find out how the co-oscillation of each boundary effected the tidal elevation along the WCK using uniform forcing at open boundaries A, B and C.

Assuming uniform boundary tide constituent for M2 or K1, an amplitude of 100 cm was forced on each boundary A, B, and C, individually. The resultant amplitudes for all gauge locations obtained from each open boundary are shown in Figure 10.

**Figure 10.** Computed amplitudes at gauge locations with respect to the assumed forcing with an amplitude of 100 cm at each open boundaries A, B and C; (**a**): M2's amplitude; (**b**): M2's amplitude.

It shows that only the tidal wave from the open boundary B produced much higher coastal amplitude than other incident boundary waves from open boundaries A and C. Whereby the tidal amplitudes of M2 and K1 approached 350 and 160 cm in the Kyunggi Bay within area I, respectively; whereas the force propagated from the boundaries A and C was damped on the way to the WCK before it approached the coast. Defant [33] also noticed that the tidal phenomenon in the entire ECS was almost exclusively conditioned by those water-masses which penetrated through the canals between the Ryukyu Islands. Figure 10a, shows the trend of M2 along the WCK obtained from the forcing at only open boundary B was similar to the results shown in Figure 9a,c,e, in which the real forcing at all three open boundaries A, B and C was taken into account. Therefore, the accuracy of boundary condition on B played a major role in the accuracy of the whole modeling results even when the real tidal condition was fully specified on three boundaries. As a result, in the following section, we have estimated the sensitivity of open boundary forcing at the boundary B on the tide around the WCK, instead of taking all three open boundaries into consideration.

Response of the Tide around the WCK to the Tidal Amplitude at Open Boundary

For the quantitative assessment of the influence of tidal amplitudes at the open boundary on the coastal tide elevation, a mean increase of amplitude at each gauge location with respect to the variation of open boundary amplitude was calculated by

$$D\_{a,i} = \frac{r\_i(a\_2) - r\_i(a\_1)}{a\_2 - a\_1} \tag{12}$$

where *a*<sup>1</sup> and *a*<sup>2</sup> are different boundary tidal amplitudes in comparison, and *ri aj <sup>j</sup>*=1,2 is a response of coastal amplitude at the *i-th* gauge location.

An increase of tidal amplitudes of constituents M2 and K1 was averaged over gauge locations within a specified region (I, G or M). As shown in Table A12, once the boundary amplitude of M2 constituent increased from 20 to 40 cm, 40 to 60 cm, 60 to 80 cm, and 80 to 100 cm, the mean coastal amplitude increased 3.25, 2.62, 2.01 and 1.91 cm per 1 cm increase in the open boundary amplitude, respectively. It shows that the higher the increment of the boundary amplitude, the lower the mean increase of the coastal amplitude per 1 cm. Meanwhile, the mean increase of the coastal amplitude of K1 was about 1.0 cm (ranged from 0.91 to 1.10), once the boundary amplitude of K1 constituent increased from 20 to 40 cm, 40 to 60 cm, 60 to 80 cm, and 80 to 100 cm.

Figure 11 shows the coastal M2 and K1 amplitudes at gauge locations corresponding to forcing with increasing amplitudes of 20, 40, 60, 80 and 100 cm at open boundary B. It shows that as the amplitudes increased at open boundary B, the M2 amplitude significantly increased once it entered the Mokpo region (M), and held this tendency until it reached to the estuary of Geum River, then it dropped until leaving the Gunsan region (G). Thereafter, it started to increase again when entering the Incheon regions (I) reaching the maximum value within Kyunggi Bay, and then significantly decreased. It means that the M2 amplitude was very sensitive within Kyunggi Bay. While K1 was linearly increased a little when it entered the WCK.

**Figure 11.** M2 and K1 amplitudes at gauge locations with respect to the forcing with the amplitudes of 20, 40, 60, 80 and 100 cm at the open boundary B; (**a**): M2's amplitude; (**b**): M2's amplitude.

(**a**) M2 (**b**) K1

In addition, Figure 11 also shows that once the tidal amplitude of M2 at boundary B was about 20 cm, the resultant amplitudes propagated within Kyunggi Bay had little different tendencies in comparison with those obtained from the amplitude larger than 20 cm, since it was not remarkably damped.

Table A12 also shows that regions I and G were most sensitive to the boundary amplitude for M2 simulation, once its amplitude was less than 60 cm. Once the amplitude of M2 increased more than 60 cm at open boundary B, it caused less impact on the increase of M2 amplitude along the WCK, since the value of *Da* did not significantly increase. Meanwhile, an increase of K1 amplitude at the boundary B with an increment of 20 cm caused less significant change on mean *Da* within all three regions of WCK.

Due to the nonlinearity in the advection terms of the governing equation, and the tidal dissipation by the nonlinear bottom shear stress on the shallow water, the increasing coastal amplitude was found to not be proportional to the increase of offshore tidal amplitude. For a more detailed analysis of the response of WCK to the tidal amplitude at open boundary B, the nonlinear response of tidal amplitude within WCK was evaluated by the following sensitivity estimation:

$$R\_{a,i} = \frac{r\_i(a\_2) / r\_i(a\_1)}{a\_2 / a\_1} \tag{13}$$

As shown in Figure 12, until the tide approached the tip of the region M, in which *Ra*,*<sup>i</sup>* was close to 1, the M2 amplitude increased proportionally to the boundary amplitude. However, when it propagated along the coastline, the increase of amplitude started to depress. As it was pointed out, a substantial depression of amplitudes occurred at some locations in the region M. The sensitivity *Ra*,*<sup>i</sup>* dropped mildly during tide propagation through the region G, then it significantly dropped in the region I. Being different from the M2 constituent, the coastal K1 amplitude was less sensitive in the entire WCK after the tip of region M.

**Figure 12.** Value *Ra* at tide station locations with respect to the increase of the amplitude of M2 (**a**) and K1 (**b**) at the open boundary B; (**a**): *Ra* of M2's amplitude; (**b**): *Ra* of K1's amplitude.

From the sensitivity analysis in this section, we can estimate the response of coastal amplitude in the WCK to the difference in open boundary amplitudes of the three established boundary conditions. As shown in Table A10, the averaged amplitude obtained from three different offshore models on the boundary B was in the range from 55.5 to 56.9 cm for M2, and from 21.4 to 21.6 cm for K1. As mentioned in the previous section, the variation of mean coastal amplitude on total gauges was approximated from 1.91 to 3.25 cm for M2, and from 0.91 to 1.10 cm (about 1 cm) for K1 per 1 cm of boundary amplitude variation shown in Table A12, respectively. In Table A11, the RMSE difference in tidal amplitude between the interpolated boundary conditions by two global models (FES2014 and TPXO9.1) and the regional model (NAO99Jb) was found to be 1.85 and 2.78 cm for M2; 0.41 and 0.42 cm for K1; respectively at the open boundary B.

#### *3.3. Comparison between the Numerical Results and Observations*

As mentioned above, many authors such as [1–3] have carried out the tidal simulations in the Yellow Sea and East China Sea (YECS) region using only four dominant tidal constituents; M2, K1, S2 and O1. Particularly, [4] analyzed 13 significant tidal constituents (M2, K1, S2, O1, MM, P1, MU2, N2, MKS2, L2, K2, M4 and MS4) to conclude that M2 was the dominant constituent, followed by S2 and K1 in this region. In this study, we continued to conduct the tidal modeling for this region, taking into account four dominant tidal constituents.

Figure 13 shows a comparison of amplitudes between the observations and numerical results of M2 and K1 at gauge stations. It shows a great improvement in comparison to Figure 2, and the numerical results implying different OBCs obtained from three assimilated tidal models (FES2014, TPXO9.1 and NAO.99Jb) were similar.

**Figure 13.** Comparison of the amplitudes between the observation and numerical simulation at tidal gauge stations (refer to Figure 6); (**a**,**c**,**e**): M2's amplitude; (**b**,**d**,**f**): K1's amplitude.

A question has arisen regarding the difference between the real water level observed around the WCK in comparison with the water level obtained from the numerical results contributed to by only four dominant tidal constituents. Figure 14 below shows a comparison of water level between numerical results and observations at typical gauge stations in three different regions (Incheon, Gunsan

and Mokpo). In detail, Table A14 shows a comparison of water levels between the observations and numerical results at 21 gauge stations. Overall, it shows that the water levels obtained from the simulation were overestimated for the neap tide, and underestimated for spring tide. Most stations show the difference in water level was below 20%; at some specific stations located in the Incheon Region, the difference was up to 30%. Nevertheless, as mentioned by the study of Teague et al. [4], after analysis on the contribution of 13 significant tidal constituents on the water level in this region, they noticed that the tides were found to account for at least 85% of the sea surface height variance in the Yellow Sea. Therefore, the water level obtained from the numerical results was reasonable, since its differences compared to the observations can be partly compensated with the amount of 15% mentioned in [4].

**Figure 14.** A comparison of water level between numerical results and observations at typical gauge stations; (**a**,**b**): water level at gauge station I3 and I5 (located in Incheon Region); (**c**,**d**): water level at gauge station G3 and G5 (located in Gunsan Region); (**e**,**f**): water level at gauge station M1 and M2 (located in Mokpo Region).

#### **4. Conclusions**

The main goal of this study was to use various numerical investigations to figure out the response of coastal tides in the WCK to the open boundary conditions and bottom roughness. After the application of three different open boundary conditions interpolated from three different assimilated tidal models (FES2014, NAO.99Jb and TPXO9.1), it has been shown that there were no significant differences between the responses of tidal amplitudes in the WCK induced by three open boundary conditions obtained from three assimilated tidal models. In addition, the numerical simulation of the tidal flow in the WCK should not use a uniform bottom roughness coefficient. Due to the complicated bathymetry, indented coastlines and bed variability of the WCK, it caused strong local effects on the tides in this region. Therefore, a non-uniform bottom roughness should be applied to the modeling whereby the smallest value can be applied for Incheon, a larger value for Mokpo, and the largest value for Gunsan. The largest value of the bottom roughness coefficient was applied to the Gunsan region because its bed form was more variable than other regions. The numerical results show that the accuracy of the modeling of the tidal elevation around the WCK was strongly dependent on the bottom roughness rather than the offshore tidal boundary conditions. Moreover, the numerical results can provide not only a better fit to the observations but also higher spatial resolutions in comparison to the results obtained from assimilated tidal models around the WCK.

However, it should be noted that the numerical results obtained from this study were still limited due to the coarse resolution (30 arcs/second) of bathymetry obtained from GEBCO2014, which was not sufficient to capture the real geometries, whose sizes were less than such resolution. Therefore, a further study with a higher resolution is necessary in order to obtain a more precise prediction of the tidal current and its elevation around the WCK. Furthermore, the wind forcing on the sea surface and the tidal energy dissipation should be taken into account.

**Author Contributions:** Conceptualization, V.T.N.; Data curation, V.T.N. and M.L.; Formal analysis, V.T.N.; Funding acquisition, V.T.N.; Investigation, V.T.N. and M.L.; Methodology, V.T.N.; Software, V.T.N. and M.L.; Supervision, V.T.N. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was funded by the National Research Foundation Grant of Korea (NRF-2018R1D1A1A09083747).

**Acknowledgments:** The authors also would like to thank the anonymous reviewers for their valuable and constructive comments to improve our manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A**


**Table A1.** Validation of amplitude and phase obtained from three tidal models against observation data.


**Table A2.** RMSE results obtained from uniform roughness bottom for FES2014.

**Table A3.** RMSE results obtained from uniform roughness bottom for NAO99Jb.



**Table A4.** RMSE results obtained from uniform roughness bottom for TPXO9.1.



**Table A6.** RMSE results obtained from non-uniform roughness bottom for M2.




**Table A8.** RMSE results obtained from non-uniform roughness bottom for S2.


**Table A9.** RMSE results obtained from non-uniform roughness bottom for O1.



**Table A10.** Averaged interpolated amplitudes of tidal constituents at open boundaries.

**Table A11.** Averaged difference of interpolated amplitude on open boundaries (cm).


**Table A12.** *Da* per 1cm increase of offshore tidal amplitude on boundary B.


**Table A13.** *Da* per 1 cm increase of the tidal amplitude of K1 on the open boundary B.



**Table A14.** A comparison of water surface level at the gauge stations.

X(\*): observation data were missing at this station.

#### **References**


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

## *Article* **Numerical and Experimental Study of Classical Hydraulic Jump**

#### **Eugene Retsinis \* and Panayiotis Papanicolaou**

School of Civil Engineering, National Technical University of Athens, 9 Iroon Polytexneiou St., 15780 Athens, Greece; panospap@mail.ntua.gr

**\*** Correspondence: retsinis76@hotmail.com; Tel.: +30-6944-916-187

Received: 20 May 2020; Accepted: 19 June 2020; Published: 21 June 2020

**Abstract:** The present work is an effort to simulate numerically a classical hydraulic jump in a horizontal open channel with a rectangular cross-section, as far as the jump location and free surface elevation is concerned, and compare the results to experiments with Froude numbers in the range 2.44 to 5.38. The governing equations describing the unsteady one-dimensional rapidly varied flow have been solved with the assumption of non-hydrostatic pressure distribution. Two finite difference schemes were used for the discretization of the mass and momentum conservation equations, along with the appropriate initial and boundary conditions. The method of specified intervals has been employed for the calculation of the velocity at the downstream boundary node. Artificial viscosity was required for damping the oscillations near the steep gradients of the jump. An iterative algorithm was used to minimize the difference of flow depth between two successive iterations that must be less than a threshold value, for achieving steady state solution. The time interval varied in each iteration as a function of the Courant number for stability reasons. Comparison of the numerical results with experiments showed the validity of the computations. The numerical codes have been implemented in house using a Matlab® environment.

**Keywords:** hydraulic jump; Boussinesq equations; MacCormack; Dissipative two-four; method of specified intervals

#### **1. Introduction**

The flow in prismatic open channels as a civil engineering specialty has been studied for over 100 years via laboratory experiments, theoretical one-dimensional analysis, as well as with numerical modeling. The flow may be characterized as a gradually varied (GVF) or a rapidly varied (RVF) flow if changes of flow depth are gradual over a significant length or abrupt over a short length of the channel respectively. A characteristic dimensionless parameter describing the type of flow is the Froude number, Fr = u/ gh that expresses the ratio of inertia over the gravity forces, u being the average, over a cross section velocity; g is the gravitational acceleration and D is the hydraulic depth, i.e., the depth of the equivalent orthogonal channel with bottom width equal to the width of free surface of the prismatic channel. When inertia forces are dominant, Fr > 1 and the flow is called supercritical, otherwise when Fr < 1, the flow is named subcritical. Critical flow (Fr = 1) may occur, but it is quite unstable, oscillating between a subcritical and a supercritical state. A special case of a RVF is the hydraulic jump, which is formed in an open channel when supercritical flow turns into subcritical, resulting in an abrupt rise of the free surface, formation of a surface roller and high energy dissipation locally. The hydraulic jump in civil engineering applications can be used as an energy dissipation mechanism, for mixing purposes as well as for the aeration of the flow.

Hydraulic jumps appear usually in stilling basins for energy dissipation of spillways when required. The sequent (subcritical) depth, the length (distance between subcritical and supercritical flow) and the location where the jump occurs are important parameters for the safe design of such structures. Among the experiments performed to determine these parameters are the pioneering works [1,2], where the one-dimensional momentum equation was derived relating the sequent depth ratio to the upstream Froude number. Moreover, in more recent experiments on hydraulic jumps [3–7], basic flow variables have been determined such as the sequent depth ratio, the length of the jump and the geometry of the surface roller.

The one-dimensional depth-averaged equations of mass conservation (continuity) and momentum (Navier–Stokes), describing the unsteady, rapidly varied open channel flow with non-hydrostatic pressure distribution, form a two-equation system named the Boussinesq equations. Assuming parallel flow and hydrostatics over the depth pressure distribution, the simplified system is called the Saint Venant (shallow water) equation. In both systems of equations the variables are the depth of flow and the average over depth velocity. A hydraulic jump can be modeled using the Boussinesq equations, which is the main objective of this work. So far, the shallow water Saint Venant equations have been used in several investigations for the numerical simulation of the classical hydraulic jump. For example, the finite element method was used [8] to solve the Saint Venant equations up to the point that steady state was reached. The assumption of hydrostatic pressure distribution employed is not valid in the region of the hydraulic jump, because the substantial curvature of the streamlines results in vertical accelerations that cannot be neglected. Moreover, the energy losses were neglected (weak hydraulic jump), making questionable the validity of simulation results for practical purposes. An attempt was made [9], to combine the knowledge regarding vorticity generation in a hydraulic jump, using the Navier–Stokes and shallow water equations. The Saint Venant equations were also modified [10], including terms related to turbulent shear stress and non-uniform velocity distribution, applying a depth-averaged approach. The equations were solved by the finite element method for the location and the profile of a hydraulic jump with upstream Froude numbers in the range 2 to 7. An effort to simulate the classical hydraulic jump using the shallow water equations and solve them using the MacCormack numerical scheme [11], produced results that compare well with the theoretical ones. The finite volume method [12] has been used to solve the governing shallow water equations. This model was applied to a hydraulic jump, and the numerical results are in agreement with experiments.

The numerical solution of Navier–Stokes equations has been widely used in modeling of the hydraulic jump. The use of the Reynolds-Averaged Navier–Stokes (RANS) approach and the k-ε turbulence closure model [13], was applied to study the turbulence characteristics of the hydraulic jump, which had to be adapted to flow cases of a moving free surface by means of the partially explicit approach. The study includes both free and forced hydraulic jumps with Froude numbers in the range from 2.1 to 7.6. Comparison of the numerical predictions with velocity measurements using one-dimensional Laser Doppler Anemometry, showed that the results qualitatively compare well, regarding the mean velocity and turbulence intensity in the flow direction. RANS equations were solved using a turbulence closure k-ε model [14], to study the hydraulic jump in a straight horizontal channel with Froude numbers 2.0 and 4.0. The mixed Eulerian-Lagrangian method for the calculation of the free surface of the flow is utilized to overcome the problem posed by a moving free boundary. A detailed study of the internal and external characteristics of a hydraulic jump was made and compared to experiments where possible. The surface roller and recirculation zone are found to play a dominant role in turbulence generation and dissipation. The study of free and submerged hydraulic jumps in a rectangular open channel downstream a sluice gate [15] has also been considered. Air-water two-phase flows are considered in the numerical simulations using RANS and turbulence closure k-ε modeling [15], and compared to velocity measurements made using Acoustic Doppler and Particle Image Velocimetry. The numerical results concerning water depths, the hydraulic jump length and the velocity profiles compare well with laboratory measurements. A three-dimensional model was implemented [16], using OpenFOAM software to analyze a hydraulic jump with upstream Froude number 6.1 in a horizontal, smooth, rectangular open channel. RANS equations were solved using three turbulence closure models, namely the Standard k-ε, the Re-Normalization Group (RNG) k-ε,

and the Shear Stress Transport (SST) k-ω model. Direct Numerical Simulation (DNS) was used [17], in a hydraulic jump with Froude number equal 2.0, using the volume of fluid method for tracking the free surface. Results have been produced regarding the mean velocity field, Reynolds stresses, turbulence production and dissipation, velocity spectra and air entrainment concentration. Finally the two recent works [18,19] present a very thorough investigation on datasets for numerical modeling assessment of the hydraulic jump.

In summary, the shallow water one-dimensional equations are capable of modeling only some characteristics of a transient hydraulic jump, while RANS combined with turbulence closure models or DNS can capture the turbulent structure of a steady hydraulic jump more accurately. Implementation of such algorithms though is limited by the Reynolds number of the flow that usually exceeds 106. Moreover, the computational cost of such models is very high, because the computational time if a decent computing machine is very long. For practical, mainly civil engineering applications, shallow water modeling is much simpler to use, while the results can be acceptable. In the present work, we will show that one-dimensional modeling of a hydraulic jump is improved using the more accurate, unsteady, one-dimensional Boussinesq equations. They are solved numerically using the MacCormack and the Dissipative two-four finite difference schemes. The practicality of such a model used in applications like the design of stilling basins, has led us to choose this model for analysis. The maximum error of the flow depth is used as an iterative parameter, for achieving the steady state solution for the free surface profile and the location of a hydraulic jump, in a horizontal rectangular open channel for a wide range of Froude numbers. To apply the Boussinesq model at the downstream boundary node, the method of specified intervals was used for the evaluation of the unknown flow variables. The results are compared to experiments performed at the Laboratory of Applied Hydraulics, School of Civil Engineering, National Technical University of Athens, Greece, showing that such numerical schemes can accurately simulate the classical hydraulic jump.

#### **2. Governing Equations**

The equations modeling unsteady one-dimensional rapidly varied flows are the Boussinesq equations [20]. These include additional terms for the non-hydrostatic pressure distribution due to curved streamlines. The basic assumptions made are: (1) the vertical velocity varies from zero at the channel bottom to its maximum value at the free surface, (2) the velocity in the main flow direction is uniformly distributed over the depth, (3) the velocity in the lateral direction is zero, (4) the fluid is incompressible, (5) the channel is prismatic of a rectangular cross section with rigid bottom and sides, (6) the longitudinal bottom slope is small, and (7) the formulas regarding energy friction slope for steady flow can be used for the unsteady flow as well. The one-dimensional Boussinesq equations for mass and momentum conservation are written in vector form as:

$$\frac{\partial \mathbf{G}}{\partial \mathbf{t}} + \frac{\partial \mathbf{F}}{\partial \mathbf{x}} = \mathbf{S}, \tag{1}$$

where:

$$\mathbf{G} = \begin{bmatrix} \mathbf{h} \\ \mathbf{u}\mathbf{h} \end{bmatrix} \mathbf{F} = \begin{bmatrix} \mathbf{u}\mathbf{h} \\ \mathbf{u}^2 \mathbf{h} + \begin{pmatrix} \frac{1}{2} \end{pmatrix} \mathbf{g} \mathbf{h}^2 - \begin{pmatrix} \frac{1}{2} \end{pmatrix} \mathbf{h}^3 \mathbf{E} \end{bmatrix} \mathbf{S} = \begin{bmatrix} \mathbf{0} \\ \mathbf{g} \mathbf{h} (\mathbf{S}\_0 - \mathbf{S}\_\mathbf{f}) \end{bmatrix} \tag{2}$$

$$\mathbf{E} = \frac{\partial^2 \mathbf{u}}{\partial \mathbf{x} \partial \mathbf{t}} + \mathbf{u} \frac{\partial^2 \mathbf{u}}{\partial \mathbf{x}^2} - \left(\frac{\partial \mathbf{u}}{\partial \mathbf{x}}\right)^2,\tag{3}$$

x is the longitudinal distance along the channel bottom as shown in Figure 1, t is the time, h = h(x,t) and u = u(x,t) are the unknown flow depth and mean velocity components in the main flow direction respectively, Sf is the energy grade line slope, So is the longitudinal bottom slope and g is the gravitational acceleration. E = E(x,t) is the Boussinesq term, which makes the difference if compared to the Saint Venant equations where the pressure distribution is assumed to be hydrostatic. The energy grade slope can be computed using the Manning formula in SI, u = (R2/3S1/2)/n, where u is the mean

over the wetted cross-section velocity, R is the hydraulic radius (cross sectional area over the wetted perimeter ratio), S = Sf the energy grade slope and n the Manning friction coefficient as Sf = n2 <sup>f</sup> <sup>u</sup>2/R4/3. Alternatively, the Darcy–Weisbach formula could be used for computing Sf, but we used the Manning formula in order to simplify calculations since the friction coefficient is kept constant for every step and the results are not affected.

**Figure 1.** Physical and computational domain of the classical hydraulic jump.

#### **3. Numerical Schemes**

The classical hydraulic jump Boussinesq equations cannot be solved analytically, therefore a numerical method has to be implemented. In this work, the MacCormack [21] and the Dissipative two-four [22] finite difference schemes have been employed for the solution of the flow equations with the appropriate initial and boundary conditions. The first scheme is second order accurate both in space and time, and the second one is fourth order accurate in space and second order in time, allowing for proper simulation of the Boussinesq terms. The developed algorithm iterates until the change of the flow depth between two successive iterations is less than a fixed convergence value. Then the jump forms a part of the steady state solution.

A sketch of the jump with the computational grid is shown in Figure 1. The hydraulic jump is formed combining the use of a sluice gate upstream with a weir downstream. The origin of the spatial coordinate (x = 0) is set at the sluice gate location, while x is the longitudinal distance along the channel bottom measured from the origin. The distance between the sluice gate and the weir is L (equal to 5.20 m, the region of interest where the computational solution is sought), and it is discretized by a total number of n nodes including the boundary nodes, thus creating a uniform grid of points at distance Δx = L/(n−1) in between. The index i denotes a spatial node, while the upstream and downstream boundaries correspond to nodes i = 1 and i = n with flow depths denoted as hup and hdo respectively. A brief presentation of the numerical schemes for the discretization of Equation (1) follows.

#### *3.1. MacCormack Scheme*

This scheme is a two-step algorithm. For the spatial derivatives of Equation (1) forward finite differences are used in the predictor step and backward finite differences in the corrector step including in the computational stencil two spatial nodes in the predictor step, the nodes i + 1, i and in the corrector step the nodes i and i − 1 as follows:

Predictor step:

$$\mathbf{G}\_{\rm i}^{\star} = \mathbf{G}\_{\rm i}^{\rm k} - \lambda \mathbf{\dot{\lambda}} \mathbf{\dot{F}}\_{\rm i+1}^{\rm k} - \mathbf{F}\_{\rm i}^{\rm k} \tag{4} \\ \text{s.t.} \tag{4}$$

Corrector step:

$$\mathbf{G}\_{\mathbf{i}}^{\*\*} = \mathbf{G}\_{\mathbf{i}}^{\*} - \lambda \mathbf{\bar{\epsilon}} \mathbf{F}\_{\mathbf{i}}^{\*} - \mathbf{F}\_{\mathbf{i}-1}^{\*} \mathbf{\bar{\epsilon}} + \Delta t \mathbf{S}\_{\mathbf{i}'}^{\*} \tag{5}$$

The flow variables at the next iteration level k+1, and grid point i, are given by:

$$\mathbf{G}\_{\mathbf{i}}^{k+1} = \frac{1}{2} (\mathbf{G}\_{\mathbf{i}}^{k} + \mathbf{G}\_{\mathbf{i}}^{\*\*}),\tag{6}$$

where λ = Δt/Δx, Δt being the time interval and the superscripts k and k + 1 refer to the two successive levels of iteration. All variables with a single asterisk (\*) refer to those computed at the predictor step where all variables with double asterisk (\*\*) refer to those computed at the corrector step. The Boussinesq terms are neglected in this scheme.

#### *3.2. Dissipative Two-Four Scheme*

The Dissipative two-four scheme consists of a predictor and a corrector step, again for the spatial derivatives of Equation (1) forward finite differences are used in the predictor step and backward finite differences having in the corrector step, as including in the computational stencil, three spatial nodes in the predictor step, the nodes i + 2, i + 1 and i and in the corrector step the nodes i, i − 1 and i − 2:

Predictor step:

$$\mathbf{G}\_{\mathbf{i}}^{\*} = \mathbf{G}\_{\mathbf{i}}^{\mathbf{k}} + \frac{\lambda}{6} \Big( \mathbf{F}\_{\mathbf{i}+2}^{\mathbf{k}} - 8\mathbf{F}\_{\mathbf{i}+1}^{\mathbf{k}} + 7\mathbf{F}\_{\mathbf{i}}^{\mathbf{k}} \Big) + \Delta t \mathbf{S}\_{\mathbf{i}}^{\mathbf{k}}.\tag{7}$$

Corrector step:

$$\mathbf{G\_i^{\*\*}} = \frac{1}{2} (\mathbf{G\_i^k} + \mathbf{G\_i^\*}) + \frac{\lambda}{12} \Big( -\mathbf{7F\_i^\*} + 8\mathbf{F\_{i-1}^\*} - \mathbf{F\_{i-2}^\*} \Big) + \frac{1}{2} \Delta t \mathbf{S\_{i'}^\*} \tag{8}$$

The vector Gk+<sup>1</sup> <sup>i</sup> at the next iteration level k + 1 and grid point i, is given by:

$$\mathbf{G\_i^{k+1}} = \frac{1}{2} (\mathbf{G\_i^k} + \mathbf{G\_i^{\*\*}}) \,\tag{9}$$

Regarding Boussinesq terms, the second order derivative ∂2u/∂x2 is approximated using a three point central finite difference in both steps, while for the first order derivative ∂u/∂x, forward finite difference is used in the predictor step, and a backward finite difference in the corrector step. The partial derivative ∂2u/∂x∂t in the Boussinesq terms is ignored, since it is zero at steady state. Regarding the code implementation, Equations (4)–(9), have to be written in non-vector form.

Due to computational reasons in the case of the Dissipative two-four scheme, the flow equations are discretized using Equations (7) and (8) at the interior nodes, i.e., for i = 4, 5, ... , n−2, the Boussinesq terms are either considered or neglected, while at nodes i = 2, 3, n−1, the MacCormack scheme is used for discretizing Equations (4) and (5) and omitting the Boussinesq terms. Apparently, this is allowed since the jump appears far from the boundary nodes, therefore the Boussinesq terms are not significant since the flow is parallel.

#### *3.3. Initial and Boundary Conditions*

The appropriate initial and boundary conditions must be set in a well-posed problem. The number of characteristics that enter the computational domain determine the auxiliary conditions to be specified [23]. The initial condition includes steady, supercritical, gradually varied flow for the entire length of the channel, resulting in two characteristic curves entering the computational domain, thus two flow parameters must be specified at each grid point, the flow depth and the velocity. The depth at each grid point can be computed by numerical integration of the first order ordinary differential equation:

$$\frac{\text{dh}}{\text{dx}} = \frac{\text{S}\_{\text{o}} - \text{S}\_{\text{f}}}{1 - \text{Fr}^{2}},\tag{10}$$

of the steady state gradually varied flow, Fr being the Froude number, Fr = u/ gh. The upstream depth hup is known (obtained from the experiments downstream of the sluice gate) and the computation is performed with a Kutta-Merson method, beginning the calculations with the flow depth hup.

The flow conditions are fixed at the two boundaries of the channel. For a given discharge and settings of the sluice gate upstream and the thin crested weir downstream, the flow depths are known and take constant values hup and hdo at nodes i = 1 and i = n respectively. The mean velocity is known at i = 1 but it has to be computed at end node i = n. Apparently, the flow is assumed to be supercritical at node i = 1 and subcritical at node i = n.

The velocity at the end node i = n will be computed using the method of specified intervals and the positive characteristic equation discretized by finite differences. From Figure 1, points A and B correspond to nodes n−1 and n respectively at the known time level k, while the positive characteristic through point P with the unknown velocity at the downstream boundary at time level k+1 is indicated. The method of specified intervals is used to compute the velocity, the celerity and the flow depth at point R, which is the intersection of the positive characteristic through point P with the grid line of the known time level k using the following relationships [24]:

$$\mathbf{u}\_{\mathbb{R}} = \frac{\mathbf{u}\_{\mathbb{B}} + \lambda \left(\mathbf{c}\_{\mathbb{B}} \mathbf{u}\_{\mathbb{A}} - \mathbf{c}\_{\mathbb{A}} \mathbf{u}\_{\mathbb{B}}\right)}{1 + \lambda \left(\mathbf{u}\_{\mathbb{B}} - \mathbf{u}\_{\mathbb{A}} + \mathbf{c}\_{\mathbb{B}} - \mathbf{c}\_{\mathbb{A}}\right)},\tag{11}$$

$$\mathbf{c}\_{\rm R} = \frac{\mathbf{c}\_{\rm B} + \lambda \mathbf{u}\_{\rm R} (\mathbf{c}\_{\rm B} - \mathbf{c}\_{\rm A})}{1 + \lambda (\mathbf{c}\_{\rm B} - \mathbf{c}\_{\rm A})},\tag{12}$$

$$\mathbf{h}\_{\rm R} = \mathbf{h}\_{\rm B} - \lambda (\mathbf{u}\_{\rm R} + \mathbf{c}\_{\rm R}) (\mathbf{h}\_{\rm B} - \mathbf{h}\_{\rm A}),\tag{13}$$

where c = gh is the celerity of a small amplitude wave propagating inside a rectangular open channel in shallow water. The energy grade line slope at point R is estimated from the following Equation:

$$\mathbf{S}\_{\mathbf{f}\_{\mathbb{R}}} = \mathbf{n}^2 \mathbf{u}\_{\mathbb{R}}^2 / \mathbf{R}\_{\mathbb{R}}^{4/3},\tag{14}$$

Then the velocity uk+<sup>1</sup> <sup>n</sup> at point P can be computed from the following relationship:

$$\mathbf{u}\_{\rm P} = \mathbf{u}\_{\rm n}^{k+1} = \mathbf{u}\_{\rm R}^{k} - \frac{\mathbf{g}}{\mathbf{c}\_{\rm R}^{k}} (\mathbf{h}\_{\rm n}^{k+1} - \mathbf{h}\_{\rm R}^{k}) + \mathbf{g} \Delta \mathbf{t} (\mathbf{S}\_{\rm o} - \mathbf{S}\_{\rm f\_{\rm R}}^{k}) \tag{15}$$

#### *3.4. Courant Condition Artificial Viscosity*

The variable time step Δt in each iteration is restricted according to the Courant criterion for stability purposes, and it is computed from:

$$
\Delta \mathbf{t} = \frac{\mathbf{c}\_{\mathbf{n}} \Delta \mathbf{x}}{\max \left( |\mathbf{u}| + \sqrt{\mathbf{g} \mathbf{h}} \right)} \text{ \tag{16}
$$

where cn is the Courant number that must be less than or equal to 0.65 [24], and Δx is the constant spatial step as shown in Figure 1.

To dampen the oscillations that occur around the jump region, artificial viscosity has to be added to the numerical schemes. The formulas used are those proposed by [25].

First the parameter ξ<sup>i</sup> at the computational node i is calculated:

$$\mathfrak{E}\_{\mathbf{i}} = \frac{\Delta \mathbf{x}}{\Delta \mathbf{t}} \frac{\left| \mathbf{h}\_{\mathbf{i}+1} - 2\mathbf{h}\_{\mathbf{i}} + \mathbf{h}\_{\mathbf{i}-1} \right|}{\left| \mathbf{h}\_{\mathbf{i}+1} \right| + 2|\mathbf{h}\_{\mathbf{i}}| + |\mathbf{h}\_{\mathbf{i}-1}|}, \text{ for the nodes, } \mathbf{i} = \mathbf{2}, \dots, \mathbf{n} - 1 \tag{17}$$

$$\mathbf{E}\_{\mathbf{i}} = \frac{\Delta \mathbf{x}}{\Delta \mathbf{t}} \frac{\left| \mathbf{h}\_{\mathbf{i}+1} - \mathbf{h}\_{\mathbf{i}} \right|}{\left| \mathbf{h}\_{\mathbf{i}+1} \right| + \left| \mathbf{h}\_{\mathbf{i}} \right|}, \text{ for the node, } \mathbf{i} = \mathbf{1} \tag{18}$$

$$\left|\xi\_{\rm i}\right| = \frac{\Delta \mathbf{x}}{\Delta \mathbf{t}} \frac{|\mathbf{h}\_{\rm i} - \mathbf{h}\_{\rm i-1}|}{|\mathbf{h}\_{\rm i}| + |\mathbf{h}\_{\rm i-1}|}, \text{ for the node, } \mathbf{i} = \mathbf{n} \tag{19}$$

Then at the center of segment i + 1/2, between nodes i and i + 1:

$$
\xi\_{i+(1/2)} = \mathbf{k}\_{\text{art}} \text{max}(\xi\_{i\prime}\xi\_{i+1}),
\tag{20}
$$

Similarity between nodes i − 1 and i:

$$
\xi\_{i-(1/2)} = \text{\tiny{\bf{k}}\_{\text{art}} \text{max}(\xi\_{i-1}, \xi\_i)},
\tag{21}
$$

where kart is a coefficient adjusting the amount of dissipation.

Finally the flow variables h and u, at iteration level k+1 are modified to the new ones according to the following relationships:

$$\mathbf{f\_{new\_i}}^{\mathbf{k}+1} = \mathbf{f\_{old\_i}}^{\mathbf{k}+1} + \boldsymbol{\xi\_{i+(1/2)}} \mathbf{(f\_{old\_{i+1}}}^{\mathbf{k}+1} - \mathbf{f\_{old\_i}}^{\mathbf{k}+1}) - \boldsymbol{\xi\_{i-(1/2)}} \mathbf{(f\_{old\_i}}^{\mathbf{k}+1} - \mathbf{f\_{old\_{i-1}}}^{\mathbf{k}+1}),\tag{22}$$

where the variable f stands either for the flow depth or for the velocity. The numerical codes have been implemented in house using Matlab® computational environment, while the input data for the numerical simulations are the channel geometry, the flow depths hup, hdo and the flow rate Q.

#### **4. Results**

The experiments were performed at the Laboratory of Applied Hydraulics of the School of Civil Engineering at the National Technical University of Athens, Greece, in a 10.50 m long flume with a rectangular cross section 0.248 m wide and 0.50 m deep. The channel has a steel bottom and vertical side walls made of glass. The water supply system consists of a constant head tank, a pipeline and a regulating valve to adjust the flow rate that was measured with a Venturi meter and did not exceed 50 L/s. The flow was controlled by a sluice gate setting the supercritical flow conditions upstream, and a thin crested weir at the end of the channel located 5.20 m downstream from the sluice gate for controlling the tailwater, subcritical flow. The hydraulic jump formed at some location between the sluice gate and the weir, depending upon the upstream, supercritical and tailwater, subcritical depths. The flow depths and free surface profile of the jumps were measured by a point gauge with accuracy ±0.1 mm.

Six experiments were implemented, the flow conditions of which are shown in Table 1. In the table appear the flow rate q, the supercritical upstream and subcritical downstream depths hup and hdo respectively, and the Froude number of the flow, Fr, at the toe of the jump from the experimental measurements. The six different jump cases measured, were computed using same conditions (upstream and tailwater depths and the flow rate) and the numerical results are compared to experiments in the following paragraphs.


**Table 1.** Experimental flow variables of the hydraulic jumps.

For the numerical modeling of the experiments shown in Table 1 three discretizing schemes of the Boussinesq equations were used: (i) the Dissipative two-four scheme with the influence of Boussinesq terms, (ii) the Dissipative two-four scheme without the influence of non-hydrostatic pressure distribution terms, and (iii) the MacCormack scheme without Boussinesq terms. The time interval computed from Equation (16) varied in every iteration. Moreover, Equation (22) was applied to dampen the oscillations near the jump. In all test cases 100 spatial nodes were used for discretization (n = 100), resulting in Δx = 0.0525 m with acceptable spatial resolution. All the results regarding the computed dimensionless flow depth h/hdo, alongside the dimensionless longitudinal distance of the channel, x/L, for each test case are shown in Figures 2a, 3a, 4a, 5a, 6a and 7a. The measured and computed flow profiles by the three numerical algorithms are plotted together along the channel, and the significance of the Boussinesq terms (due to non-parallel streamlines) inside the region of the hydraulic jump are depicted sideways in Figures 2b, 3b, 4b, 5b, 6b and 7b. Outside the jump these terms are almost zero (small enough), thus confirming the hydrostatic pressure distribution. From the numerical results, it must be noted that the flow depth increases from the vena contracta downstream of the sluice gate up to the jump, due to energy losses encountered in the governing equations. One may note that in test case 1, the dissipative behavior of the numerical results is probably due to the artificial dissipation added in the numerical scheme. In addition, from the Boussinesq term distribution along the channel, one may note that this term is not zero at the upstream end, therefore we would not expect a smooth transition in free surface from supercritical to subcritical flow. Furthermore, note that the Froude number is quite low and the hydraulic jump is characterized as weak [26], since 1.7 < Fr < 2.5, which means that it is a non-fully developed jump because of the weak surface roller with low energy loss. From these figures we conclude that the agreement between experiments and computations is satisfactory.

**Figure 2.** Experimental and Numerical Results for test case 1: (**a**) Free surface profile; (**b**) Boussinesq term versus distance.

**Figure 3.** Experimental and Numerical Results for test case 2: (**a**) Free surface profile; (**b**) Boussinesq term versus distance.

**Figure 4.** Experimental and Numerical Results for test case 3: (**a**) Free surface profile; (**b**) Boussinesq term versus distance.

**Figure 5.** Experimental and Numerical Results for test case 4: (**a**) Free surface profile; (**b**) Boussinesq term versus distance.

**Figure 6.** Experimental and Numerical Results for test case 5: (**a**) Free surface profile; (**b**) Boussinesq term versus distance.

**Figure 7.** Experimental and Numerical Results for test case 6: (**a**) Free surface profile; (**b**) Boussinesq term versus distance.

Figure 8a,b exhibit the variation of the dimensionless spatially (over a cross-section) averaged flow velocity in the main flow direction x, u/uup, where uup is the average velocity at x = 0, with the dimensionless longitudinal distance, x/L, for all three algorithms for test cases 1 and 6 respectively. It is evident that the MacCormack scheme overestimates the flow velocity at the upstream end of the channel, while the results from the other two methods are almost identical. Figure 9a,b depict the evolution of convergence criteria of the maximum absolute change in flow depth between two successive iterations, until a steady state is obtained for test case 2 using the Dissipative two-four scheme with Boussinesq terms, and test case 5 using the MacCormack scheme without Boussinesq terms respectively. Figures 10 and 11 present the temporal evolution in "computational-pseudo time" of the free surface elevation for the Dissipative two-four scheme including the Boussinesq terms for test cases 3 and 4 respectively, both in dimensionless form. Note that the jump moves upstream until it is stabilized. Similar results have been produced for all other test cases.

**Figure 8.** Variation of the mean velocity alongside the channel: (**a**) for test case 1; (**b**) for test case 6.

**Figure 9.** Evolution of the convergence criteria: (**a**) for test case 2 with Dissipative scheme with Boussinesq terms; (**b**) for test case 5 with MacCormack scheme without Boussinesq terms.

**Figure 10.** Temporal evolution of jump formation for test case 3 using the Dissipative two-four scheme including the Boussinesq terms.

**Figure 11.** Temporal evolution of jump formation for test case 4 using the Dissipative two-four scheme including the Boussinesq terms.

The specific force of the flow (momentum flux per unit weight) in a channel with rectangular cross section written as <sup>M</sup> = 0.5bh2 + Q2/gbh , where b is the width of the rectangular cross section, is conserved across a hydraulic jump, or alternatively, it takes the same value upstream and downstream of the jump. From the numerical results it turns out that the relative error of the specific force between the two end cross-sections of the jump for all test cases and the three numerical schemes applied ranged from 0.62 % to 12.92 %, results that are acceptable for the accurate location of the jump. The number of iterations required to obtain a steady state with accuracy of order 10−<sup>4</sup> m or 5 <sup>×</sup> <sup>10</sup>−<sup>5</sup> m, along with the closure of continuity equation for each test case and all numerical schemes are shown in Table 2, where the flow rate was computed from the depth and the average velocity over the cross section. It is evident that for all test cases the MacCormack algorithm gave the highest error in mass conservation, if compared to the other algorithms, due to the lower order of spatial accuracy.


**Table 2.** Required iterations for convergence and mass balance error.

An iterative algorithm was implemented for both numerical schemes which converged to the steady state solution, when the maximum value of the difference of the flow depth between two successive iterations was smaller than a threshold value. The sequent depths as well as the location where the jump forms for given flow and boundary conditions, resulted from the steady state solution of the governing equations. The numerical results showed that the magnitude of the Boussinesq terms is greater at the jump regime, whereas outside this it is almost zero, as expected. It was also shown that the Boussinesq terms do not affect practically the location where the jump forms for test cases 1–6, so that their exclusion eases the discretization without significant error. Note that to obtain the results presented in this paper artificial viscosity has been used. A coefficient adjusting the amount of dissipation to reproduce the experimental measurements by a trial and error procedure was set to the value 0.011.

#### **5. Case Studies**

For further validation of the numerical results obtained from the six test cases examined, let us present a case study of practical interest. We consider a 100 m long and 2 m wide (narrow channel, depth to width ratio around one) rectangular, horizontal open channel, with Manning roughness coefficient 0.013 and flow rate 7.2 m3/s. In Figure 12a–c we present the computed dimensionless flow depth, h/hdo, versus the dimensionless distance x/L, for three Froude numbers at the vena contracta cross-section namely 2, 4 and 6, including or excluding the Boussinesq terms, using a dissipation parameter equal to 0.011. In Figure 12d we demonstrate the temporal evolution of the jump until it is stabilized, for Froude number equal to 4. In all three cases one may notice the formation of H3 and H2 types of free surface curves upstream and downstream of the jump respectively, due to energy loss. From Figure 12b,c it is evident that the hydraulic jump is much longer when computed using Boussinesq terms.

**Figure 12.** Dimensionless free surface elevation for 2 m wide channel: (**a**) for Fr = 2; (**b**) for Fr = 4; (**c**) for Fr = 6; (**d**) temporal evolution of the jump for Fr = 4.

Let us now examine the influence of the channel width on the length of the hydraulic jump, Lj. We consider a 100 m long channel and 5 m wide (wide section, depth to width ratio smaller than one) with orthogonal cross section conveying 10 m3/s, with same dissipation factor equal to 0.011. In Figure 13a–c the computed dimensionless flow depth, h/hdo, is plotted versus the dimensionless distance x/L, for three Froude numbers at the vena contracta namely 2, 4 and 6, with or without the Boussinesq terms incorporated. The flow depth at x/L = 0 is the boundary condition, and it is the same with and without Boussinesq terms, as indicated in Figures 12 and 13. Near the supercritical boundary the depth deviates from that at x/L = 0, with deviation being greater when Boussinesq terms are not included. This is probably due to numerical instability from the boundary condition, also shown clearly in Figures 2b, 3b, 4b, 5b, 6b and 7b, where Boussinesq terms are plotted vs x/L. The non-Boussinesq

terms showed greater instability because the depth near the supercritical boundary deviates more than that with Boussinesq terms.

The length of the hydraulic jump computed with Boussinesq terms, is still greater than that computed without incorporating Boussinesq terms, but the difference is not as sound as it was in the previous example. In Tables 3 and 4 we present the length of the computed hydraulic jump and compare it to that from experiments [27] (p. 14). The results in the 5 m wide channel with or without the Boussinesq terms are in satisfactory agreement with measurements in orthogonal channels, if compared to the results of the 2 m wide channel. The noticeable differences between the flow profiles in Figures 12 and 13 are probably due to the non-hydrostatic pressure distribution encountered in Boussinesq terms. The difference is sound when the hydraulic jump is "weak" for low Froude numbers, where the transition from supercritical to subcritical flow is not as vigorous as for greater Froude numbers. Therefore there is a smoother free surface elevation when Boussinesq terms are included, resulting in longer jump lengths.

**Figure 13.** Dimensionless free surface elevation for 5 m wide channel: (**a**) for Fr = 2; (**b**) for Fr = 4; (**c**) for Fr = 6.


**Table 3.** Length of jump for 2 m wide open channel.

From Table 3, it is obvious that comparing the numerical results with experiments, exclusion of Boussinesq terms leads to more accurate results concerning the length of the jump. Alternatively, from Table 4 the opposite is true. Apparently, the flow in a "wide" orthogonal open channel is not affected from the side walls while in a "narrow" channel, it is. Therefore, the one-dimensional shallow water equations that are improved with Boussinesq terms can predict the characteristics of a hydraulic jump accurately, if the channel is wide if compared to the flow depth.


**Table 4.** Length of jump for 5 m wide open channel.

#### **6. Conclusions**

The one-dimensional Boussinesq equations have been solved numerically using the MacCormack as well as the Dissipative two-four finite difference schemes, for the simulation of hydraulic jump formation in a horizontal rectangular open channel and for upstream Froude numbers Fr in the range 2.44 to 5.38. The governing equations have been enriched with additional terms if compared to the Saint Venant equations, to account for the non-hydrostatic pressure distribution in the regime of rapidly varied flow. Terms related to the energy loss and the gravity forces have been also included. The initial condition was a steady supercritical gradually varied flow along the whole length of the channel modeled. The upstream and downstream boundary conditions regarding the flow depth remained constant during the iteration process, and equal to the values measured in experiments. The method of specified intervals was used for the calculation of the velocity at the downstream end, assuming that the positive characteristic through a point does not intersect with already established grid points. Variable time step was used in every iteration according to the CFL (Courant-Friedrichs-Lewy) stability criterion, along with artificial viscosity for smoothing of the oscillations occurring in the jump. The computational results compare well with experiments since the specific force was computed from the depth and mean velocity at both ends of the hydraulic jump with acceptable tolerance, and the mass conservation equation was verified for all numerical schemes and all test cases.

From such a model one can determine the sequent depth ratio as well as the length of the jump, results that are useful in the design of stilling basins (geometrical properties). Given a stilling basin with a known inflow Froude number and flow depth, the engineer must decide the end sill dimensions and the basin length, so that the hydraulic jump is contained in the stilling basin. Finally, from comparison of the numerical results and experiments, it can be concluded that the aforementioned numerical modeling schemes can predict the basic features of the classical hydraulic jump with acceptable accuracy.

**Author Contributions:** E.R. has performed the experimental measurements and implemented the numerical algorithms. E.R. wrote the paper with the contribution of P.P. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** The authors would like to thank J. Demetriou for his contribution in the implementation of the experiments.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


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

## *Article* **An Optimized and Scalable Algorithm for the Fast Convergence of Steady 1-D Open-Channel Flows**

#### **Louis Go**ffi**n \*, Benjamin Dewals, Sebastien Erpicum, Michel Pirotton and Pierre Archambeau**

Hydraulics in Environmental and Civil Engineering (HECE) Research Unit, Urban and Environmental Engineering (UEE) Department, Faculty of Applied Sciences, University of Liege (ULiege), Allée de la Découverte, 9-4000 Liège, Belgium; b.dewals@uliege.be (B.D.); s.erpicum@uliege.be (S.E.); michel.pirotton@uliege.be (M.P.); pierre.archambeau@uliege.be (P.A.)

**\*** Correspondence: l.goffin@uliege.be

Received: 11 October 2020; Accepted: 11 November 2020; Published: 17 November 2020

**Abstract:** Calculating an open-channel steady flow is of main interest in many situations; this includes defining the initial conditions for the unsteady simulation or the computation of the water level for a given discharge. There are several applications that require a very short computation time in order to envisage a large number of runs, for example, uncertainty analysis or optimization. Here, an optimized algorithm was implemented for the fast and efficient computation of a 1-D steady flow. It merges several techniques: a pseudo-time version of the Saint-Venant equations, an evolutionary domain and the use of a non-linear Krylov accelerator. After validation of this new algorithm, we also showed that it performs well in scalability tests. The computation cost evolves linearly with the number of nodes. This was also corroborated when the execution time was compared to that obtained by the non-linear solver, CasADi. A real-world example using a 9.5 km stretch of river confirmed that the computation times were very short compared to a standard time-dependent computation.

**Keywords:** shallow water; CasADi; fast computation; 1-D

#### **1. Introduction**

Channelized flows can be simulated using 1-D, 2-D or 3-D models depending on the level of flow detail that is required [1]. Results from hydrodynamic numerical models are used in multiple domains including flood risk analysis [2] or real-time control of river facilities [3], for instance.

One-dimensional models are used when a dominant direction can be assumed in the velocity field. This may be the case when the flow is restricted to the main riverbed. A 1-D model can still be used in the case of out-of-bank flooding, although they are unable to represent complex 2-D flow patterns in the floodplain [4]. Several practical cases have shown that flood mapping can be performed using 1-D models [5,6] and they can also be used in other fields, such as flood routing for hydropower plant operations [7] and mixed flows in pipes [8].

In fact, 1-D models are still used extensively even though 2-D and 3-D models are currently widely available. There are various reasons for their use, for example, digital elevation models (DEM) and bathymetry data are not available in some regions of the world. When only cross-section profiles are available to represent the geometry of a riverbed, 1-D models can make direct use of such profiles whereas they have to be interpolated to be used in 2-D or 3-D models. Besides, many applications do not require a detailed description of the flow features in the floodplain, and thus, a 1-D modeling approach is sufficient.

Large-scale hydraulic modeling of river networks [9–13] makes heavy use of 1-D models because of their ability to compute long stretches of the river at a reasonable cost. To initiate a computation, one needs boundary conditions as well as the initial condition. This initial condition is often a steady water profile, which can be obtained by performing a time-dependent simulation with steady boundary conditions over a period of time that is long enough to reach a steady solution [14]. This step may consume a considerable amount of time before the main problem can be addressed. The initial condition should be computed with the same numerical scheme as the one used in the unsteady model in order to ensure the steadiness of the first step of the unsteady model. The ability of these models to quickly obtain a steady initial solution is also of great importance.

Optimization is another field where obtaining a steady result as quickly as possible is important. Indeed, most optimization techniques require a large number of runs in order to figure out the optimal solution. In order to ensure that the overall computation time is as short as possible, techniques that are quick should be utilized including parallelization [15] or the use of fast computing models.

Although 1-D simulations are known to provide results in a short period of time, accelerating the computation of the 1-D steady solution is of great interest in the fields mentioned above. This can be achieved by using two main strategies. The first consists of exploiting the resources of modern computers more efficiently. Such techniques are more frequently applied to 2-D cases, which naturally require more computing resources. Common hardware acceleration strategies include parallelizing codes on several CPU cores [16,17] or on a GPU [18–20]. The second method consists of designing algorithms in order to converge with less effort toward the solution. To the authors' knowledge, there is no such work available in the literature. Both strategies can be combined in order to obtain the best performance.

The purpose of this paper is to propose algorithmic strategies for the fast computation of 1-D steady solutions. First, the equations are presented. Then, we introduce two original strategies in order to reduce the overall computation time. These strategies can easily be implemented in other hydraulic codes. An alternative non-linear solver is introduced for comparison purposes. Finally, the validation results, the optimal parameters for the algorithm, the performance of the model and its application to a real-world problem are presented.

#### **2. Materials and Methods**

#### *2.1. Equations*

One-dimensional water flow is described by the St-Venant equations, which are as follows [21]:

$$\begin{aligned} \frac{\partial \underline{A}}{\partial t} + \frac{\partial \underline{Q}}{\partial \mathbf{x}} &= q\_l\\ \frac{\partial \underline{Q}}{\partial t} + \frac{\partial}{\partial \mathbf{x}} \left(\frac{\underline{Q}^2}{A}\right) + gA \frac{\partial z\_s}{\partial \mathbf{x}} &= -gAS\_f + u\_x q\_l \end{aligned} \tag{1}$$

where *A* is the cross-section area (m2), *Q* is the discharge (m3/s), *ql* is a lateral discharge per unit length (m2/s), *g* is the gravity acceleration (m/s2), *Sf* is the friction slope (-), *ux* is the velocity along the *x* direction of *ql* (m/s) and *zs* is the free surface elevation (m). The numerical resolution of this set of non-conservative equations has been shown to provide accurate results in many practical cases, including in the presence of discontinuities [22,23].

Assuming a steady flow, temporal derivatives vanish in Equation (1). Since solving the first equation is straightforward, the discharge distribution is known on the entire domain for a given distribution of *ql*. We assume that the sign of *Q* is independent from *x*. It means that a single equation remains with a single unknown *A*. In order to keep the same numerical scheme as the one used for the unsteady system (which is important when a steady solution is used as the initial condition of an unsteady problem and to be able to use a similar algorithmic strategy), Kerger et al. [14] added a pseudo-temporal term (pseudo-time is τ (s)) to the steady form of Equation (1):

$$
\beta \frac{\partial A}{\partial \tau} + \frac{\partial}{\partial \mathbf{x}} \left( \frac{Q^2}{A} \right) + gA \frac{\partial z\_s}{\partial \mathbf{x}} = -gAS\_f + \mu\_\mathbf{x} q\_l \tag{2}
$$

where β = −sign(*Q*). Kerger et al. [14] justify this choice for β by analyzing the characteristic velocity of Equation (2):

$$
\lambda = \frac{c^2 - \mu^2}{\beta} = \frac{c^2 \left(1 - \text{Fr}^2\right)}{\beta} \tag{3}
$$

where *c* (m/s) is the wave celerity. In subcritical flows (Fr < 1, Fr = *Q*/ *gA*3/*b* 0.5 is the Froude number (-), *b* is the width of the free surface (m)), and the sign of λ is the sign of β. If Fr > 1, then sign(λ) = −sign(β). For critical flows (Fr = 1), the characteristic velocity is 0, independently from β. In order to keep some form of consistency with the general model, if we choose β = −sign(*Q*), an upstream boundary condition is required when Fr > 1 and a downstream boundary condition is required when Fr < 1. This is equivalent to the position of the water depth boundary condition for the numerical resolution of the full 1-D set of equations.

With a single boundary condition and a discharge distributed in the channel, solving Equation (2) determines the cross-section area (and subsequently, the water depth) all along the stretch. Equation (2) is discretized according to the finite volume method. It is solved according to the same numerical scheme as the one used for the full unsteady model [14,24].

For a node *i*, Equation (2) is discretized in finite volumes as:

$$\frac{Q\_{i+1/2}^2/A\_{i+1/2} - Q\_{i-1/2}^2/A\_{i-1/2}}{\Delta x} + gA\_i \Big(\frac{z\_{s,i+1/2} - z\_{s,i-1/2}}{\Delta x} + S\_{f,i}\Big) - u\_{x,i}q\_{l,i} \approx -\beta \frac{\partial A}{\partial \tau} \tag{4}$$

where Δ*x* (m) is the spatial discretization step and subscripts refer to the position of variables values.

For the sake of clarity, we consider *Q* > 0 and a constant reconstruction of the flux at finite volume boundaries to explicate the numerical scheme. When applying the considered upwinding directions [14] for a node *i* not located next to a boundary, Equation (4) is equivalent to:

$$\frac{Q\_i^2 / A\_i - Q\_{i-1}^2 / A\_{i-1}}{\Delta x} + gA\_i \Big( \frac{z\_{s,i+1} - z\_{s,i}}{\Delta x} + S\_{f,i} \Big) - \mu\_{x,i} \eta\_{l,i} = \frac{\partial A}{\partial \tau} \tag{5}$$

This flux vector splitting method has been shown to be unconditionally stable [14]. For the node located at the downstream boundary (*i* = *N* − 1), if a weak water level boundary condition *zs,BC* is imposed at the border, Equation (4) becomes:

$$\frac{Q\_{N-1}^2/A\_{N-1} - Q\_{N-2}^2/A\_{N-2}}{\Delta x} + gA\_{N-1} \left(\frac{z\_{s,BC} - z\_{s,N-1}}{\Delta x} + S\_{f,N-1}\right) - \mu\_{x,N-1}q\_{l,N-1} = \frac{\partial A}{\partial \tau} \tag{6}$$

Without a boundary condition imposed on the value of *zs* at the external border, a nil *zs* gradient is imposed and Equation (4) becomes:

$$\frac{Q\_{N-1}^2/A\_{N-1} - Q\_{N-2}^2/A\_{N-2}}{\Delta x} + gA\_{N-1}S\_{f,N-1} - u\_{x,N-1}q\_{l,N-1} = \frac{\partial A}{\partial \tau} \tag{7}$$

At the upstream node (*i* = 0), if a weak boundary condition of the water level is imposed at the border, Equation (4) becomes:

$$
\log A\_0 \left( \frac{z\_{s,1} - z\_{s,BC}}{\Delta x} + S\_{f,0} \right) - \mu\_{x,0} q\_{l,0} = \frac{\partial A}{\partial \tau} \tag{8}
$$

Without a boundary condition at the upstream border, Equation (4) becomes:

$$
\delta\_0 g A\_0 \Big(\frac{z\_{\mathbf{s},1} - z\_{\mathbf{s},0}}{\Delta x} + S\_{f,0}\Big) - u\_{\mathbf{x},0} q\_{l,0} = \frac{\partial A}{\partial \pi} \tag{9}
$$

#### *2.2. Original Solving Strategy*

Solving Equation (2) instead of Equation (1) decreases the computation time since the number of equations and unknowns is reduced. In order to save even more time, two additional strategies were implemented: (a) a non-linear Krylov accelerator was used to promote fast convergence and (b) the computation was only performed on a sliding part of the full domain.

#### 2.2.1. Non-Linear Krylov Acceleration

Numerically solving a non-linear system can be performed by different means, including Newton's method and Broyden's method [25]. More sophisticated methods exist in order to solve non-linear systems faster, such as the Jacobian-free Newton–Krylov method [26] and Anderson acceleration [27]. The Anderson acceleration method uses the results from successive iterations in order to adapt the new approximation. Walker and Ni [28] showed that this method can be considered equivalent to the well-known GMRES method [29] when applied to linear systems. The nonlinear Krylov acceleration (NKA) [30,31], which is similar to Anderson acceleration, has been found to be more efficient in some applications than more recent methods such as the Jacobian-free Newton–Krylov method [32]. NKA was used for a faster convergence of our pseudo-time model.

Since NKA only relies on the results directly produced by the hydraulic model, it can be easily applied to other algorithms or other domains. Indeed, no gradient evaluation (nor other prerequisite) is required before calling on the NKA algorithm.

The NKA algorithm records *<sup>N</sup>* (*<sup>N</sup>* <sup>∈</sup> <sup>N</sup>>0) previous moves of the root finding process. Based on these previous moves, NKA adapts its guess for the new root. One of the main assumptions is that the Jacobian matrix remains constant within the scope of *N* moves. We briefly explain the method here and extended details can be found in [30–33].

A Newton–Raphson iteration process computes the *n* + 1st guess of the root based on the *n*th guess **x***n*, on the value of the function *f* at **x***n* and on the invert of the Jacobian matrix J:

$$\mathbf{x}\_{n+1} = \mathbf{x}\_n - \mathbf{J}^{-1} f(\mathbf{x}\_n) \tag{10}$$

Instead of evaluating J−<sup>1</sup> at each iteration, NKA evaluates it from *N* previous moves or sets it to the identity matrix, in which case the method degenerates to the fixed-point method.

NKA takes advantage of a history of *N* corrections of **x** (denoted **v**) and *N* evolutions of *f*(**x**) (denoted **w**) at iterate *n*:

$$\begin{aligned} \mathbf{v}\_i &= \mathbf{x}\_i - \mathbf{x}\_{i-1} \\ \mathbf{w}\_i &= f(\mathbf{x}\_i) - f(\mathbf{x}\_{i-1}) \end{aligned} \quad i = n - N + 1, \dots, n \tag{11}$$

The method assumes that J is constant and invertible within the scope of the *N* previous iterations, which is written as follows:

$$\begin{array}{l} \mathbf{J} \mathbf{v}\_{i} = \mathbf{w}\_{i} \\ \mathbf{v}\_{i} = \mathbf{J}^{-1} \mathbf{w}\_{i} \end{array} \tag{12}$$

Mathematical developments described in the references cited above lead to the expression:

$$\begin{cases} \mathbf{v}\_{n+1} = \sum\_{i=n-N+1}^{n} z\_i \mathbf{v}\_i + \left( f(\mathbf{x}\_n) - \sum\_{i=n-N+1}^{n} z\_i \mathbf{w}\_i \right) \\ \mathbf{x}\_{n+1} = \mathbf{x}\_n + \mathbf{v}\_{n+1} \end{cases} \tag{13}$$

where the coefficients *zi* are the solution of the projection:

$$\mathbf{z} = \operatorname\*{argmin}\_{a \in \mathbb{R}\_0} \left\| f(\mathbf{x}\_n) - \sum\_{i=n-N+1}^n a\_i \mathbf{w}\_i \right\| \tag{14}$$

Equation (13) shows that the correction of the variable **x** is decomposed into two components:


Note that the Jacobian matrix is not used in this iterative process.

#### 2.2.2. Evolutionary Domain

The second strategy was designed to reduce computational cost. It consists in reducing the size of the computation domain and sliding it along the river stretch in order to evaluate the cross-section areas from downstream to upstream.

The development of this strategy has arisen from the long history of numerical hydrodynamics in various flow regimes. Indeed, many flows that are solved for rivers are subcritical (Fr < 1) at downstream and upstream boundary conditions. In such a situation, the cross-section area information is propagated from downstream to upstream. For a fixed flow direction, the upwinding direction of the scheme takes all unknowns and properties (except those linked to the discharge) downstream. This means that when a new node is computed, it depends only on the downstream nodes (Figure 1). It should be noticed that when a node *i* is added, the upstream border of node *i* + 1 produces a change in the upwinding direction of property *Q* and the unknown *Avel* (cross-section area used to compute the velocity). The node *i* + 1, which was supposed to be converged, undergoes a new convergence process that indirectly affects nodes *i* + 2, *i* + 3, ... Theoretically, all the nodes located downstream of node *i* should be kept in the computational domain.

**Figure 1.** Upwinding directions of the unknown and flow properties at the upstream limit of the computation domain.

The boundary condition on the border between the node *i* and node *i* − 1 (see Figure 1) is called a "mirror border". On this border, all the unknowns and properties are reconstructed from the computed inner node. This method is equivalent to reproducing the node *i* in *i* − 1, like a mirror. For node *i* in Figure 1, the discretization of Equation (2) is:

$$\lg A\_i \left( \frac{z\_{s,i+1} - z\_{s,i}}{\Delta x} + S\_{f,i} \right) - \mu\_{x,i} q\_{l,i} = 0 \tag{15}$$

The evolution of the domain relies on three strategies. First, all nodes in the computation domain should have reached a partial residual threshold ξ*<sup>p</sup>* (m2/s) before considering the extension of the computation domain: ∂*A*/∂τ ≤ ξ*p*. Then, the most downstream node can be removed from the partial domain once it drops under the final residual threshold <sup>ξ</sup>*<sup>f</sup>* (m2/s): <sup>∂</sup>*A*/∂τ <sup>≤</sup> <sup>ξ</sup>*<sup>f</sup>* . Finally, we have to make sure that the nodes that were removed from the domain are not impacted by the computation domain in such a way that their residuals ∂*A*/∂τ become higher than the final precision threshold ξ*<sup>f</sup>* . Dimensionless parameters are discussed further in Section 3.2.

In order to assess whether the removed nodes ∂*A*/∂τ remain lower than ξ*<sup>f</sup>* , an analytical analysis of the discretized model is performed. Let us consider three nodes and their borders (Figure 2). The discretized formulation of Equation (2) at node *i* is:

$$\begin{aligned} \left[ \frac{\partial A}{\partial \tau} \right]\_i &= \left[ \frac{\partial}{\partial x} \left( \frac{Q^2}{A} \right) + \mathfrak{g}A \frac{\partial z\_r}{\partial x} + \mathfrak{g}AS\_f - \mathfrak{u}\_x q\_l \right]\_i \\ &\approx \frac{\left( \frac{Q^2\_i}{A\_i} - \frac{Q^2\_{i-1}}{A\_{i-1}} \right)}{\Delta x} + \mathfrak{g} \frac{A^2\_{i+1} - A^2\_i}{A\_{i+1} - A\_i} \left( \frac{z\_{s,i+1} - z\_{s,i}}{\Delta x} \right) + \left[ \mathfrak{g}AS\_f - \mathfrak{u}\_x q\_l \right]\_i \end{aligned} \tag{16}$$

**Figure 2.** Upwinding directions of the unknown and flow properties in the computation domain.

The influence of a change in the cross-section area in *i* − 1 on the value of ∂*A*/∂τ in *i* is given by the derivative:

$$\frac{\partial}{\partial A\_{i-1}} \left( \left[ \frac{\partial A}{\partial \tau} \right]\_i \right) = \frac{Q\_{i-1}^2}{\Delta x} \frac{1}{A\_{i-1}^2} \tag{17}$$

From the result in (17), one can predict the evolution of the residual ∂*A*/∂τ of a node that is not in the computation domain anymore. If node *i* is planned to be removed from the computation domain, one can check that its ∂*A*/∂τ remains below the threshold ξ*<sup>f</sup>* even with an evolution of the cross-section area in *i* − 1:

$$
\left[\frac{\partial A}{\partial \tau}\right]\_{i, \text{uvct}} \approx \left[\frac{\partial A}{\partial \tau}\right]\_{i, \text{pvcv}} + \frac{Q\_{i-1}^2}{\Delta x} \frac{1}{A\_{i-1}^2} \Delta A\_{i-1} \tag{18}
$$

where Δ*Ai*−<sup>1</sup> is the evolution of the cross-section *A* in *i* − 1 from the last evaluation of [∂*A*/∂τ]*i*,*prev*.

If the evolution of the cross-section area in the computation domain (*i* − 1) implies that [∂*A*/∂τ]*<sup>i</sup>* exceeds the threshold ξ*<sup>f</sup>* , then node *i* should be added again in the computation domain in order to make sure it remains below the threshold until the end of the entire computation. It should be noted that node *i* + 1, *i* + 2, ... might also be impacted. This technique guarantees that once the sliding domain has moved along the entire domain, all nodes have reached at least the final precision required.

As stated earlier, the method described here was designed within the framework of subcritical flows. In order to be able to deal with a larger range of flow regimes, several adaptations were made. When a supercritical node is detected downstream (let say at position *m*), it is not computed and the computation domain is extended until a subcritical node is found upstream (say at position *n*). Then, the domain starting from *m* to *n* is computed and converged. This technique avoids boundary condition problems. Indeed, a supercritical flow requires an upstream BC since the characteristic velocity is directed towards the downstream.

#### 2.2.3. Combination of Solving Strategies

The combined use of the sliding domain and NKA involves several specific considerations. The use of NKA is implemented in the code with a safety coefficient that deactivates this optimizing technique in some cases. Indeed, it was experienced that NKA could lead to some instabilities when there was a sudden change in the cross-section area. This behavior is due to the assumption in NKA that the Jacobian matrix is constant and invertible locally [32,33]. In order to avoid such a situation, the accelerator is deactivated when *Ai*+<sup>1</sup> <sup>−</sup> *Ai* > η*Ai*+1, where *Ai* and *Ai* <sup>+</sup> <sup>1</sup> refer to the cross-section area of the nodes *i* and *i* + 1 as depicted in Figure 1, and 0 < η ≤ 1 is the safety coefficient (-). After several tests, we found that η = 0.5 provides stability with a limited impact on the computation time.

The method can be summarized with the following pseudo-code (Algorithm 1):


#### Remove this node from computation list

**If** upstream nodes remain to be added to computation list:

Add 1 upstream node to the computation list

#### Finalize

#### *2.3. Alternative Non-Linear Solver*

Various techniques can be used to solve nonlinear Equation (2). Up to this point, we have chosen to discretize the equation using a finite volume scheme and solve it with an explicit time scheme, which is consistent with the unified strategy of WOLF [24]. Other techniques can be used, including finite difference schemes and/or implicit time schemes. Another possibility is to use an optimization algorithm for nonlinear systems. One of these is the recently developed CasADi software [34,35].

CasADi first started as an algorithmic differentiation tool. During its evolution, developers chose to shift the focus toward optimization. From non-linear expressions, CasADi is able to generate all the information needed by a nonlinear solver in order to return a solution to the problem. CasADi provides interfaces to MATLAB or Python for easy use.

The purpose of using CasADi is to show how our algorithm performs compared to a state-of-the-art solver.

The implementation in CasADi was done through Opti stack, a collection of helper functions used for nonlinear programming problems. It is possible to define variables to optimize, parameters, an objective function and constraints. The solving of a 1-D steady open channel flow can be done thanks to this framework.

The constraints of the problem are discretized in Equation (2) for each node and a water depth above 0 everywhere. The downstream boundary condition is imposed through Equation (6). If no boundary condition is set, then a flow condition can be imposed through a constraint on the Froude number for the downstream node and the flow head is minimized at the upstream node. If the flow presents a critical section, minimizing the head upstream is equivalent to finding the section with the highest critical head.

Another way to solve a flow with a critical section is to prescribe a Froude number transition fromFr < 1 to Fr > 1 at that critical section. This is done by setting a constraint on the Froude number on the nodes upstream and downstream of the critical section. The identification of the critical section should be done prior to the computation on the basis of a critical head analysis.

#### **3. Results and Discussion**

This section presents the validation of the results and focuses on the optimal parameters and performance. The geometries of the tests were different in order to examine as many cases as possible.

#### *3.1. Validation*

The validation of the models was performed on a bump placed in a straight horizontal channel, considering three different flow conditions. The bump and channel geometry have been described previously in [36]. The whole domain ranges from 0 to 20 m with the following bed elevation:

$$z\_{\mathbb{B}}(\mathbf{x}) = \begin{cases} \ 0.8 \Big( 1 - \frac{\left(x - 10\right)^2}{4} \Big) 8 \text{ m} \le x \le 12 \text{ m} \\ 0 \end{cases} \tag{19}$$

The channel is considered to be rectangular. The discretization step was chosen as 0.1 m.

The different flow conditions are described in Table 1. All tests were performed with ξ*<sup>p</sup>* = 10−<sup>6</sup> m2/s and ξ*<sup>f</sup>* = 10−<sup>10</sup> m2/s.


**Table 1.** Boundary conditions for three test cases.

The objective was to show that the model is able to deal accurately with various flow regimes and transitions. Test A simulates a fully subcritical flow with no transition. Test B creates a subcritical flow upstream, a subcritical flow downstream and a hydraulic jump in between, downstream of the bump. Finally, the goal of test C is to show the robustness of the method for a downstream supercritical flow and an upstream subcritical flow.

The analytical solutions for tests A, B and C were computed from the Bernoulli principle (head conservation) [22] and the conjugate water depth formula was used for test B. A graphical comparison of the analytical values and numerical results obtained by our algorithm and CasADi is given in Figures 3–5. It appears that the new algorithm and CasADi provide results that fit well with the analytical solution. However, a small difference in energy can be noticed between the analytical solution and the numerical results and also between both numerical methods (see Table 2). This was quantified and explained by Bruwier et al. [22]. Moreover, a more noticeable localized difference appears between CasADi and the new algorithm results at the hydraulic jump. Even if the numerical scheme is the same for the new algorithm and CasADi, each method has its own convergence threshold. Altogether, the analysis validates both models.

**Figure 3.** Comparison between the analytical solution and the numerical results produced by the new algorithm and CasADi for validation test A.

**Figure 4.** Comparison between the analytical solution and the numerical results produced by the new algorithm and CasADi for validation test B.

**Figure 5.** Comparison between the analytical solution and the numerical results produced by the new algorithm and CasADi for validation test C.

**Table 2.** Upstream head values for tests A, B and C and differential to the analytical value.


#### *3.2. Optimal Setting of the New Algorithm*

Five test cases were defined in order to specify the optimal values for the parameters for the new algorithm. These five tests were designed in order to induce changes in the flow characteristics due to topographic or cross-section variations. The first three tests (1 to 3) concern a channel with a rectangular cross-section and a bed slope that follows a sine function:

$$z(\mathbf{x}) = a \sin(\beta \pi \mathbf{x}) + \mathbf{y} \tag{20}$$

with *x* ∈ [0; 20], α = 0.05, γ = 0.05, β = 1/2 (for tests 1 and 2) and β = 2 for test 3. Two hundred nodes were used to discretize the 20 m long channel, resulting in a 10 cm spatial step. For test 1, the downstream boundary condition is a free surface elevation imposed at 1.2 m. For tests 2 and 3, the same type of boundary condition was imposed with a smaller value of 0.55 m, which results in a higher Froude number downstream. The specific discharge was imposed upstream at 1 m2/s for tests 1 to 3. The Manning equation was used for friction in tests 1 to 3, with the Manning coefficient *n* = 0.04 s/m1/3.

The topography and the hydraulic solutions were found thanks to the principle of head conservation (Bernoulli) and are shown in Figures 6–8 for tests 1 to 3. Table 3 summarizes the characteristics of each test. The objective of tests 1 to 3 was to analyze the influence of a variation of the bed topography on the behavior of the sliding domain. For test 1, the irregularity of the bed has only a slight influence on the water level. For the other tests, the higher Froude number and less spaced bed elevation variations were meant to investigate the possible influence of oscillations in the water level on the sliding domain performance.

**Figure 6.** Hydraulic solution for test 1.

**Figure 7.** Hydraulic solution for test 2.

**Figure 8.** Hydraulic solution for test 3.

**Table 3.** Summary of the 5 tests used to investigate the best parameters for the new algorithm.


Tests 4 and 5 were performed on regular bottoms. The discontinuities that we wanted to explore here are linked to a change in the flow regime due to the presence of a weir (test 4) or a severe change in the cross-section (test 5). For test 4, the topography was set at *z* = 0 m for all nodes except for three of them: *z* = 0.5 m at *x* = 9.85 m and *x* = 10.05 m, and *z* = 1 m at *x* = 9.95 m. Friction was computed with the Manning formula and a coefficient *n* = 0.04 s/m1/3. The cross-section is uniform along the channel and is rectangular with a width of 1 m. A discharge of 1 m2/s was injected upstream and a water depth equal to 0.7 m was imposed downstream.

Test 5 deals with a severe change in the cross-section on an inclined bottom. The channel extends 100 m, discretized with 200 nodes. The slope is 0.2%. The cross-section is trapezoidal upstream (*x* ≤ 47.5 m), then suddenly becomes rectangular in the middle of the channel (47.5 m < *x* < 52.5 m), and finally returns to a trapezoidal shape in the downstream part (*x* ≥ 52.5 m). The trapezoidal sections have a width at the bottom of 2 m and the banks are inclined with an angle of 45◦. The rectangular cross-sections are 1 m wide. The friction and boundary conditions are the same as in test 4.

The topography and final water levels for tests 4 and 5 are depicted in Figures 9 and 10. A summary of these tests is given in Table 3.

**Figure 9.** Hydraulic solution for test 4 on the entire domain (**a**) and zoomed on the weir (**b**).

**Figure 10.** Hydraulic solution (**a**) and cross-section change (**b**) for test 5.

Our algorithm includes many parameters that need to be specified. These parameters are the partial residual threshold ξ*p*, the final residual threshold ξ*<sup>f</sup>* , the temporal scheme to solve Equation (4) and the coefficient for the deactivation of the nonlinear Krylov accelerator. The final and partial residuals are two closely linked parameters. They also have a direct impact on the computation time. For a given final residual, which has to be parametrized by the user, the partial residual influences the number of iterates required to converge the partial domain and the size of these domains. After investigation, the other two parameters were shown to have almost no influence on the computation time. The

following results focus on the best value to use for the partial residual threshold ξ*<sup>p</sup>* for a fixed value of ξ*<sup>f</sup>* = 10−<sup>8</sup> m2/s.

CPU times were measured on a desktop computer (Intel i7 3.5 GHz CPU) for the five tests and various partial residual thresholds. The results are reported graphically in Figure 11. It appears that the overall computation time decreases with an increase in the partial residual threshold. Some stagnation appears around 10−<sup>2</sup> m2/s for tests 3 and 5. This can be explained by the fact that the residual naturally decreases at each iteration. Keeping some nodes in the computation domain results in a decrease in the residual for each node included in the computation domain.

**Figure 11.** CPU time according to the partial residual threshold for tests 1 to 5.

In order to assess the efficiency of the new algorithm, two scalability tests were performed. The first one consisted of extending the domain upstream, with a constant spatial discretization. The second test consisted of keeping the same channel but refining the discretization and then increasing the number of computation nodes.

The first test took place on a frictionless sine bed elevation described by:

$$z(\mathbf{x}) = 0.05 \sin(2\pi \mathbf{x}) + 0.05 \tag{21}$$

The downstream boundary condition is a water level imposed at 0.65 m. The discharge is constant along the channel stretch and is equal to 1 m2/s. Cross-sections are rectangular and 1 m wide. Four domain lengths were tested (20 m, 200 m, 2 km and 20 km) with a spatial step of 0.1 m, meaning that these domains include 200, 2000, 20,000 and 200,000 nodes.

The computation results (Figure 12a) showed that the computation time per node decreases when the number of nodes increase. This can be explained by the fact that when the domain gets longer, the flow conditions upstream are smoother than downstream. Longer domains undergo fewer changes than shorter domains, leading to shorter computation time per node. This example shows that higher partial residual values provide the best computation times.

In order to complete this scalability study, we looked at the behavior of the algorithm when the spatial step decreases for a given domain length. In classical explicit schemes, this case leads to a quadratic increase in the computation time. Indeed, when the spatial step decreases, the number of nodes increases and the time step decreases.

A 100 km-long channel with a constant 0.025% slope was chosen to illustrate the behavior of the new algorithm. The cross-sections are trapezoidal and are described using tabular values (1 m wide at the bottom of the section and 5 m wide at 1 m above the bottom). Friction was generated using a Manning law with *n* = 0.03 s/m1/3. The downstream boundary condition is a water level set at 1 m, and 1 m3/s is injected at the upper node and the injection of 4 m3/s is shared amongst the other nodes (through the *ql* term in Equation (1), *ux* = 0).

**Figure 12.** CPU time evolution with the number of nodes when (**a**) the spatial step is kept constant and when (**b**) the domain span is fixed.

The computation times are showed in Figure 12b. The behavior is slightly different from that observed in the previous test. Indeed, the evolution of the computation time is linear in regard to the number of nodes (the CPU time per node is globally constant) when the partial residual is set at 10−<sup>6</sup> and 10−<sup>8</sup> m2/s. For partial residual values of 10−<sup>2</sup> and 10−<sup>4</sup> m2/s, the evolution is linear up to 50,000 nodes; however, for the finest discretization, the computation time increases in a nonlinear way.

This point was investigated further. It appears that at some moment in the computation, the algorithm needs to increase its computation domain size without being able to reduce it quickly (i.e., upstream nodes are added to the computation list while downstream nodes cannot be removed for residual values reasons). This increase in the number of nodes in the computation list was nonlinear compared to the situation with coarser discretization.

We looked at dimensionless values for parameters <sup>ξ</sup>*<sup>p</sup>* and <sup>ξ</sup>*<sup>f</sup>* by dividing them by <sup>√</sup>*gA*3/4 <sup>0</sup> , *A*0, which is a characteristic cross-section area. The results obtained showed rather constant values, suggesting that a coherent choice for ξ*<sup>f</sup>* / <sup>√</sup>*gA*3/4 0 should be around 10−10. We also found that an efficient ratio ξ*<sup>f</sup>* /ξ*<sup>p</sup>* is around 10<sup>−</sup>6.

#### *3.3. Performance of the Models*

The test case chosen for the comparison between CasADi and our algorithm is a rectangular channel (1 m wide) with a 100 m long sine wave bottom as shown in the following equation:

$$z(\mathbf{x}) = 0.05 \sin(\frac{\pi}{4}\mathbf{x}) + 0.05\tag{22}$$

A Manning friction formula with *n* = 0.04 s/m1/<sup>3</sup> was used to estimate friction losses. The downstream boundary condition is a water depth equal to 0.6 m. The unit discharge is uniform and is equal to 1 m2/s. The precision parameters are as follows: ξ*<sup>p</sup>* = 10−<sup>6</sup> m2/s and ξ*<sup>f</sup>* = 10−<sup>8</sup> m2/s.

We compared the CPU time spent in the solving stages of our new algorithm and CasADi. The goal was to compare the evolution of the computation time with the number of nodes, rather than comparing absolute values. The results are given in Figure 13. They confirm that the computation time evolves almost linearly with the number of nodes when the new algorithm is used. This is not the case with CasADi: the computation time increases following a power law *N*α, α > 1. This is due to the matrix operations that CasADi has to perform. Increasing the number of nodes leads to a non-proportional increase in computation time. The speed up factor, which is the ratio between the CPU time spent when using our new algorithm and the CPU time spent with CasADi, ranges from in order of 10<sup>0</sup> to in order of 102 according to the number of nodes considered.

**Figure 13.** Comparison of the CPU time evolution with the number of nodes between our algorithm and CasADi and the associated speed up factor.

#### *3.4. Real-World Application*

The Romanche River in the French Alps is currently facing a number of significant changes. A new hydropower facility is being built in order to replace older power plants. In this context, the dam operator needs a fast computation routine in order to operate its facilities in an optimized and safe way. To do so, an unsteady 1-D model was implemented in Fortran and integrated in a Simulink S-Function in order to be compatible with the operator model [7]. Our new algorithm was used for the fast computation of an initial condition.

The studied part of the Romanche River stretches over 9.5 km (Figure 14), from the new Livet Dam to Gavet. The geometry and the calibration of the friction coefficients are detailed in [7]. The river was discretized with 191 nodes (50 m-long meshes). The downstream boundary condition was set at an elevation of 436.8 m. Two discharges were tested: 10 and 40 m3/s. Since the details of the hydraulic results are of limited interest for this paper, we focused on the execution time and showed that the Froude number remains under 1 for both discharges (Figure 14). For the same machine as the one used earlier, the execution times (CPU time) are 0.018 s and 0.022 s, respectively (ξ*<sup>f</sup>* = 10−<sup>8</sup> m2/s and ξ*<sup>p</sup>* = 10−<sup>2</sup> m2/s). In comparison, the computation time of a five-minute full simulation was about 0.2 s. Given the length of the stretch and the celerity of a wave, the time for a wave to travel from downstream to upstream is of the order of 80 min. This means that approximately 3.2 s of computation time are required is order to simulate the propagation of a wave downstream to upstream. It is clear that the gain of time offered by our original technique is significant.

**Figure 14.** Bottom elevation and Froude number distribution for both simulated discharges.

#### **4. Conclusions**

Several innovations are introduced in this paper, including the use of the non-linear Krylov accelerator in open-channel flows, an evolutionary domain algorithm and the use of CasADi to solve steady 1-D flows. These improvements lead to an algorithm that is able to quickly solve steady open-channel flows. Therefore, optimization problems and uncertainty analyses that require many evaluations, become more tractable.

An original algorithm was implemented in order to significantly improve the computation time of a steady 1-D open-channel flow problem. It includes two main optimizing strategies: a non-linear Krylov accelerator and an evolutionary domain algorithm. This new algorithm was validated against the academic benchmarks of flows over a bump. The results showed a good agreement between the numerical and analytical values.

The performance of the suggested algorithm was evaluated against the non-linear optimization software CasADi. It showed good scalability properties. Indeed, the execution time of the proposed algorithm evolves linearly with the number of nodes. This is not the case with other techniques when the mesh is refined and/or when the number of nodes increase.

Finally, we demonstrated the capabilities of our algorithm in a real-world case. We used the optimized algorithm in order to compute quickly the initial condition required by the operational model for the Romanche River in France. Our technique was able to provide a steady state solution to the unsteady model in a very short period of time.

**Author Contributions:** Conceptualization, M.P. and P.A.; methodology, L.G. and P.A.; software, L.G. and P.A.; validation, L.G.; writing—original draft preparation, L.G.; writing—review and editing, P.A., M.P., B.D. and S.E.; supervision, M.P. and P.A. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflict of interest.

**Computer Code and Software:** The following software and codes were used for this paper: (a) WOLF was developed by the HECE research group (http://www.hece.ulg.ac.be/cms/) at the University of Liège since 2000 and is not freely available, (b) CasADi is freely available at https://web.casadi.org/, and (c) the routines used to test CasADi are freely available at https://gitlab.uliege.be/HECE/HydroCasADi.

#### **References**


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

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

### *Article* **Efficient Reservoir Modelling for Flood Regulation in the Ebro River (Spain)**

**Isabel Echeverribar 1,2,\*, Pablo Vallés 1, Juan Mairal <sup>1</sup> and Pilar García-Navarro <sup>1</sup>**


**Abstract:** The vast majority of reservoirs, although built for irrigation and water supply purposes, are also used as regulation tools during floods in river basins. Thus, the selection of the most suitable model when facing the simulation of a flood wave in a combination of river reach and reservoir is not direct and frequently some analysis of the proper system of equations and the number of solved flow velocity components is needed. In this work, a stretch of the Ebro River (Spain), which is the biggest river in Spain, is simulated solving the Shallow Water Equations (SWE). The simulation model covers the area of river between the city of Zaragoza and the Mequinenza dam. The domain encompasses 721.92 km<sup>2</sup> with 221 km of river bed, of which the last 75 km belong to the Mequinenza reservoir. The results obtained from a one-dimensional (1D) model are validated comparing with those provided by a two-dimensional (2D) model based on the same numerical scheme and with measurements. The 1D modelling loses the detail of the floodplain, but nevertheless the computational consumption is much lower compared to the 2D model with a permissible loss of accuracy. Additionally, the particular nature of this reservoir might turn the 1D model into a more suitable option. An alternative technique is applied in order to model the reservoir globally by means of a volume balance (0D) model, coupled to the 1D model of the river (1D-0D model). The results obtained are similar to those provided by the full 1D model with an improvement on computational time. Finally, an automatic regulation is implemented by means of a Proportional-Integral-Derivative (PID) algorithm and tested in both the full 1D model and the 1D-0D model. The results show that the coupled model behaves correctly even when controlled by the automatic algorithm.

**Keywords:** reservoir model; numerical simulation; shallow water equations; PID regulation

#### **1. Introduction**

As extreme phenomena, flood events raise concern among governments, institutions and general society. The European Union has been developing plans and directives during the last decades focusing on the control of their impact [1]. River overflows cause the flooding of adjacent lands, urbanised areas and other infrastructures. Additionally, floods can also take human lives, as reported by the UN [2], specially in areas with poor prevention plans and a lack of predictive tools. Frequently, dams and reservoirs are present in river basins as hydraulic elements with different functions. Not only to ensure enough water supply for agricultural activities or energy production, but also as hydraulic structures for discharge adjustment and control during flood events. Basin authorities manage their operation focusing on available space in the reservoir, maximum acceptable downstream discharges, and peak arrival times.

In this context, the development of predictive tools that provide information about the temporal and spatial evolution of water level and discharge along a river during flood events can help to quantify the damage caused and has been widely addressed in last decades [3]. Some works are focused on urban areas coupling their overland models with sewer systems [4,5]. Some others are more oriented to large scale floods quantifying

**Citation:** Echeverribar, I.; Vallés, P.; Mairal, J.; García-Navarro, P. Efficient Reservoir Modelling for Flood Regulation in the Ebro River (Spain). *Water* **2021**, *13*, 3160. https:// doi.org/10.3390/w13223160

Academic Editors: Anargiros I. Delis and Ioannis K. Nikolos

Received: 22 September 2021 Accepted: 12 October 2021 Published: 9 November 2021

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

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

inundated areas in crops and surrounding fields [6,7], including components such as weirs, gates, dams and reservoirs [8,9]. Nowadays, there are even operational tools developed to simulate extremely large domains using massive parallel algorithms [10]. In either case, all models are then used to generate information and data for other models and tools that analyse and classify flood-prone areas [11,12]. Fast and efficient numerical models for the resolution of the equations that govern free surface flows in rivers have been developed and improved in recent years [6,7,10,13–17]. However, including a reservoir in the numerical model of the river reach is an additional difficulty. Even today, it is widely accepted that three dimensional models (3D) are computationally expensive and the phenomenon of a flood in a river at a large scale can be addressed by averaging the equations vertically (2D approximation) [6,7]. However, flood event simulations often involve large domains and long time scales, and practical applications require a compromise between spatial resolution and computational efficiency [18,19]. To achieve the necessary spatial resolution, in many cases quite fine computational grids are needed, so more data storage is required, proportionally increasing the number of operations and reducing the size of the time step allowed for explicit calculations. Therefore, flood risk evaluations are often performed considering the average in the cross section to reduce the phenomenon to a 1D approximation [13]. Finally, depending on terrain morphology, some particular river reaches might transport hydrographs almost immediately, as reservoirs.

On the other hand, reservoirs can be assessed with different approximations depending on the variables of interest. Reservoirs can be solved either as part of the river reach; this is, discretised. Alternatively, they can be considered as a storage volume with a constant level [20,21]. When the detailed phenomena that might occur within a reservoir are of interest, complex numerical models are developed. In [22], a 3D model based on the Reynolds-Averaged Navier-Stokes equations is used to study secondary currents and three dimensional behaviour of the velocity field. When the details of the 3D velocity field are not required but the longitudinal profile of the water surface is of interest, they can be incorporated as part of 1D discretised models [23]. Additionally, when some other additional phenomena must be simulated, such as eutrophization [24] or sediment transport [25], also a spatial discretisation of the reservoir is needed, whether in one or three dimensions. An alternative is an aggregated reservoir routing where only a volume balance is considered [26–28]. This may include runoffs, evaporation and some other mass exchanges. In any case, each reservoir must be specifically analysed. Aggregated models may be suitable, due to the representation of the reservoir as a unique volume, providing CPU times in the order of seconds [28], and a discretised model must compute as many operations as grid elements, leading to higher computational times [22]. However, the main disadvantage of these simplified approaches is that they do not represent in detail the flow behaviour of the river and floodplains [29–31].

Concerning optimization of the reservoir as a hydraulic structure, several works have focused on their hydropower potential. In [32], for instance, a linear optimization model with three different objective functions was implemented to automatically manage the reservoir in order to maximize total energy.

The aim of the present work is to couple recent research tools based on shallow water numerical models for flood forecasting with an aggregated model for the reservoir, developing a complete efficient simulation tool during flood events.

First, a comparison of the results obtained with a 2D and a 1D model in the middle reach of the Ebro river is carried out for validation purposes. Both, the 2D and the 1D model, are based on a finite volume scheme, which uses terrain data for the 2D mesh and 1D bathymetry creation. The 1D modelling is likely to lose the detail of the floodplain, but nevertheless the computational cost is expected much lower compared to the 2D model. Additionally, the particular nature of this reservoir, which is highly channelised, might turn the 1D discretisation into a more suitable option than the 2D discretisation. Therefore, the quality of the 1D results should be checked and the calculation times of both models should be compared.

Due to the particularity of the simulated section and the Mequinenza reservoir, which transports flood waves almost immediately, an aggregated alternative technique is applied. This approach formulates the reservoir flow globally by means of a volume balance model. Our focus is on checking if the results obtained are similar to those provided by the fully 1D model and comparing the computation times of both simulations. This later option is completed with a PID algorithm for regulation purposes.

#### **2. Study Area**

The Ebro River basin is one of the largest drainage areas in the Iberian Peninsula, as seen in Figures 1a,b, where the Ebro River represents the biggest river in Spain. In this work, a stretch of this Ebro River is simulated solving the Shallow Water Equations (SWE). The simulation model, delimited in dashed line in Figure 1c, covers the area between the city of Zaragoza and the Mequinenza dam, encompassing 721.92 km2. Between the inlet and outlet locations, there are 221 km of river bed of which the last 75 km belong to the Mequinenza reservoir, where the dynamics of the river changes to be nearly at rest. During the entire stretch, the river descends from 208 m.a.s.l. of the elevation in Zaragoza up to approximately 60 m.a.s.l. at the bottom of the riverbed in the Mequinenza dam, leaving an average slope of 6 per 10,000.

**Figure 1.** Location of Spain in Europe (**a**); location of the Ebro River basin in Spain (**b**) and location of the computational domain of the study in the basin (**c**).

The Ebro river is managed by the Ebro River Authority (CHE, www.chebro.es), which controls, rules and prepares the reports of the basin (http://www.chebro.es/contenido. visualizar.do?idContenido=14093&idMenu=3048; accessed on 20 October 2021). CHE monitors the evolution of the flow discharge and water levels at a few control stations along the river course storing data every 15 min. Figure 2 represents the Ebro River reach simulated in this work. In the figure, the most important cities and the CHE gauging stations available for data comparison are marked. The represented domain coincides with the 2D domain used for simulations. The gauging stations in this reach are located in Pina, Villafranca and Gelsa. Additionally, a point of estimation exist near the Mequinenza dam. Each of these stations has an official label that can be seen in the same figure. This region is of special interest due to its agricultural activity, and frequently suffers flooding with important damages. It is a river reach where two different parts can be identified: the first part of the region is characterised by marked meanders and large flooding areas; the second part, around 75 km of the reach, is dominated by the large Mequinenza reservoir provided with vertical walls. The dynamics of the river changes in the reservoir: its velocity is reduced until water is practically at rest and flood waves are transported almost instantly.

The Mequinenza reservoir, the largest in the entire region, covers a surface area of about 7540 hectares, with a maximum capacity of 1530 hm3 at a maximum normal surface level of 121 m.a.s.l. The reservoir is exploited for hydroelectric production and irrigation to nearby agricultural areas. At the same time, with its 124 m crest above sea level and its 6 gates, the dam is used to regulate the water storage in order to dump peak discharge during floods and to guarantee hydroelectric generation. At the reservoir, there is only a measurement point for water level that is transformed by CHE into a discharge value through volume estimations. CHE uses two different approaches for discharge estimation, so that the results in the Mequinenza reservoir are compared not with observed but with 2 different estimated data.

**Figure 2.** Representation of the 2D simulation domain of the Ebro River with the most important cities and gauging stations of CHE. The labels correspond to the official names of the gauging stations.

Two historical events of the Ebro River, the 2015 and the 2018 floods, have been identified as relevant. Information concerning discharge hydrographs as well as time evolution of the water surface level are available at the gauging stations. Additionally, the European Emergency System (EMS) provides data of flooded area extensions, as seen in Figure 3 (https://emergency.copernicus.eu/; accessed on 20 October 2021). This helped to choose the domain extension setting the boundaries far enough not to interfere the flow.

**Figure 3.** Zoom view of an ortophoto with measured extension of the flooded area (blue) in 2018 flood event.

#### **3. Methodology**

Derived from the Navier-Stokes equations by depth averaging and assuming hydrostatic pressure, the Shallow Water Equations (SWE) can be considered to govern the free surface flow of a river.

#### *3.1. Two Dimensional (2D) Model*

The 2D model can be compactly formulated as

$$\frac{\partial \mathbf{U}}{\partial t} + \frac{\partial \mathbf{F}(\mathbf{U})}{\partial x} + \frac{\partial \mathbf{G}(\mathbf{U})}{\partial y} = \mathbf{S}(\mathbf{U}) \tag{1}$$

with:

$$\mathbf{U} = \begin{pmatrix} h \\ hu \\ hv \end{pmatrix} \quad \mathbf{F}(\mathbf{U}) = \begin{pmatrix} hu \\ hu^2 + gh^2/2 \\ huv \end{pmatrix} \tag{2}$$

$$\mathbf{G(U)} = \begin{pmatrix} hv \\ huv \\ hv^2 + gh^2/2 \end{pmatrix} \quad \mathbf{S(U)} = \begin{pmatrix} 0 \\ gh(S\_{ox} - S\_{fx}) \\ gh(S\_{oy} - S\_{fy}) \end{pmatrix} \tag{3}$$

in terms of the water depth, *h*, the depth averaged unit discharges *hu* and *hv* in the *x* and *y* directions respectively. The slopes *Sox* and *Soy* are the two components of the bottom surface gradient *zb*(*x*, *y*):

$$S\_{ox} = -\frac{\partial z\_b}{\partial x} \quad S\_{oy} = -\frac{\partial z\_b}{\partial y};\tag{4}$$

and *Sf x* and *Sf y* represent friction slopes, that are here formulated as:

$$S\_{fx} = \frac{n^2 \mu \sqrt{\mu^2 + \upsilon^2}}{h^{4/3}} \quad S\_{fy} = \frac{n^2 \upsilon \sqrt{\mu^2 + \upsilon^2}}{h^{4/3}} \tag{5}$$

where *n* stands for the semiempirical Manning friction coefficient ([33]).

#### *3.2. One Dimensional (1D) Model*

When the equations are averaged over the cross sectional area of the flow, a 1D model is obtained, representing changes along the longitudinal direction of the river. The obtained system is analogous to the 2D system, with a mass conservation equation and a linear momentum equation along the river channel:

$$\frac{\partial \mathbf{U}}{\partial t} + \frac{\partial \mathbf{F}(\mathbf{U})}{\partial \mathbf{x}} = \mathbf{S}(\mathbf{U}) \tag{6}$$

with:

$$\mathbf{U} = \begin{pmatrix} A \\ Q \end{pmatrix} \quad \mathbf{F} = \begin{pmatrix} Q \\ Q^2/A + gI\_1 \end{pmatrix} \quad \mathbf{S} = \begin{pmatrix} 0 \\ g\left[I\_2 + A(\mathbb{S}\_0 - \mathbb{S}\_f)\right] \end{pmatrix} \tag{7}$$

where *Q* stands for transversal discharge, *A* is the cross section wetted area and *I*1, *I*<sup>2</sup> are hydrostatic pressure integrals. *S*<sup>0</sup> is the bottom slope along the longitudinal coordinate of the channel:

$$S\_0 = -\frac{\partial z\_b}{\partial \mathbf{x}}\tag{8}$$

and *Sf* is the friction slope, that is also formulated through the Manning law as:

$$S\_f = \frac{Q^2 n^2}{A^2 R^{4/3}}\tag{9}$$

where *R* is the hydraulic radius, defined as *R* = *A*/*P*, being *A* the wetted area and *P* the wetted perimeter. Finally, the Manning friction coefficient is obtained empirically ([33]).

#### *3.3. Finite Volume Model for the 1D Flow Equations*

In this work, an explicit upwind first order finite volume method is used for both systems of equations ([34–36]). The systems of Equations (1) and (6), can be generally expressed as:

$$\frac{\partial \mathbf{U}}{\partial t} + \overrightarrow{\nabla} \mathbf{E} = \mathbf{S} \tag{10}$$

where **E** = **F G** in the 2D model, **F** and **G** are given by (2) and (3). Moreover, **E** = **F** in the 1D model, with **F** given by (7). When integrating (10) into a control volume or cell, Ω and applying the divergence theorem, the following expression is obtained:

$$\frac{d}{dt} \int\_{\Omega} \mathbf{U} \, d\Omega + \oint\_{\partial \Omega} \mathbf{E}(\mathbf{U}) \cdot \hat{n} \, dl = \int\_{\Omega} \mathbf{S}(\mathbf{U}) \, d\Omega \tag{11}$$

where *n*ˆ is the outward unit vector in the normal direction to the volume Ω. From this, the 1D approach is next developed and details for the 2D approach can be found in [7,37]

Due to the hyperbolic character of the 1D equations, the numerical scheme used to solve them is based on the Jacobian matrix of the fluxes:

$$\mathbf{J} = \frac{\partial(\mathbf{E} \cdot \hat{\mathbf{n}})}{\partial \mathbf{U}} \xrightarrow[\mathbf{I} \cdot \mathbf{D}]{} \mathbf{J} = \frac{\partial \mathbf{F}}{\partial \mathbf{U}} = \begin{pmatrix} 0 & 1\\ c^2 - u^2 & 2u \end{pmatrix} \tag{12}$$

whose eigenvalues are:

$$
\lambda\_1 = \mathfrak{u} - \mathfrak{c} \quad \lambda\_2 = \mathfrak{u} + \mathfrak{c} \tag{13}
$$

with *u* and *c* given by:

$$
\mu = \frac{\mathcal{Q}}{A} \quad \mathcal{c} = \sqrt{\mathcal{g} \, A/B} \tag{14}
$$

being *A* the cross section and *B* the free surface width. The celerity *c* characterizes the speed of the infinitesimal surface deformation waves defining the dimensionless Froude number *Fr* = *<sup>u</sup> c* .

Following [16,35,38], the final updating scheme for a cell *i* of the domain in the time *t <sup>n</sup>*+<sup>1</sup> takes into account the contributions of neighbour cells containing fluxes and source terms as:

$$\mathbf{U}\_{i}^{n+1} = \mathbf{U}\_{i}^{n} - \frac{\Delta t\_{1D}}{\Delta x} \left[ \sum\_{m=1}^{2} \left( \tilde{\lambda}^{+} \tilde{\gamma} \tilde{\mathbf{e}} \right)\_{i-1/2}^{m} + \sum\_{m=1}^{2} \left( \tilde{\lambda}^{-} \tilde{\gamma} \tilde{\mathbf{e}} \right)\_{i+1/2}^{m} \right]^{n} \tag{15}$$

being *λ*˜ and **e**˜, respectively, the eigenvalues and eigenvectors of the Jacobian matrix of the flux, **J**, linearised on the cell edge. Additionally, Δ*x* stands for the cell size. The upwind scheme sends the information to the wave propagation direction through the eigenvalues and their sign:

$$
\tilde{\lambda}\_{i+1/2}^{\pm m} = \frac{1}{2} (\tilde{\lambda} \pm |\tilde{\lambda}|)\_{i+1/2}^m \tag{16}
$$

This scheme for the ordinary computational cells must be complemented with proper initial and boundary conditions. The numerical scheme is stabilised by dynamically limiting the time step size, Δ*t*1*D*, with the CFL condition:

$$
\Delta t\_{1D} = \text{CFL } \min\_{m,k} \left( \frac{\Delta x}{|\vec{\lambda}\_k^m|} \right) \tag{17}
$$

where 0 < CFL ≤ 1 [39].

#### *3.4. Reservoir Model*

In the context of the 1D shallow water model, the reservoir can be modelled with two different approaches. In both cases, sketched in Figure 4, the upstream river reach is discretised with a 1D finite volume method. However, the two approaches differ in reservoir representation:


As depicted in Figure 4, the aim of approach (b) is to remove the computational cells needed for the reservoir. Using the sketch as an example, while the fully 1D distributed model encompasses from x*<sup>L</sup>* = 0 to x*<sup>L</sup>* = *L*, the coupled 1D-0D model has a discretised domain only from x*<sup>L</sup>* = 0 to x*<sup>L</sup>* = *L* , so that the computational cost is reduced.

**Figure 4.** Representation of the two different approaches for reservoir representation: 1D distributed discretisation as the rest of the domain (**a**); or coupling the 1D model of the river with an aggregated 0D model of the reservoir (**b**).

In near-rest flows, such as those in a reservoir, the velocity field is negligible so that it is likely that the flow behaviour is properly solved only with volume balance. This represents a 0D approximation. The *Modified Puls Method* [21] is based on the hypotheses that the flow surface is always horizontal, the stored volume in the reservoir can be formulated as a function of water level (*V* = *V*(*H*)) and the outlet discharge can also be expressed as a water level function (*Qout* = *Qout*(*H*)). Thus, the volume variation is given by the difference between the reservoir inlet, *Qin*, and outlet, *Qout*, as:

$$\frac{dV}{dt} = Q\_{in} - Q\_{out} \tag{18}$$

Discretising Equation (18) in time by assuming <sup>Δ</sup>*<sup>V</sup>* <sup>=</sup> *<sup>S</sup>*(*Hn*)(*Hn*+<sup>1</sup> <sup>−</sup> *<sup>H</sup>n*) with *<sup>S</sup>* the free surface reservoir area (*S* = *f*(*H*)) leads to:

$$H^{n+1} = H^n + \frac{\Delta t}{S(H^n)} \left( \frac{Q\_{in}^{n+1} + Q\_{in}^n}{2} - \frac{Q\_{out}^{n+1} + Q\_{out}^n}{2} \right) \tag{19}$$

When combining this formulation with the 1D model to lead to the 1D-0D model, the water level calculated with expression (19) is set at *L* (Figure 4b). It is important to note, that the resolution of (18) requires knowing *Qin*(*t*) and the relation between *Qout* and volume *V* at the reservoir. The inflow discharge to the reservoir is directly given by the computation at the last cell of the 1D model. On the other hand, the reservoir outflow discharge depends on the geometry and characteristics of the dam. In the present work the downstream boundary condition, either for pure 1D or for 1D-0D model, is based on a weir/dam law of the form [40]:

$$Q^{n+1} = \frac{2}{3} \sqrt{2g} b \mathcal{C} H\_w^{4/3} + \frac{8}{15} \sqrt{2g} \tan\left(\frac{\theta}{2}\right) \mathcal{C} H\_w^{5/2} \tag{20}$$

where *Hw* = *H* − *hCrest* is the water depth above the weir crest. Assuming a trapezoid shape, *b* is the width of the minor length of the horizontal sides of the weir and *θ* the opening angle of the trapezoid. Finally, *C* is an energy loss coefficient here taken as *C* = 0.611 [9,26].

It is important to note that for both approaches, the full 1D model and the 1D-0D model, the 1D river reach upstream the reservoir must be identically discretised, this is, using the same Δ*x*, so that the analysis can show the differences provoked by the reservoir model.

#### *3.5. PID Regulation*

The Mequinenza dam gates can be manually operated at present according to energy, agricultural or safety criteria. In this work, a control Proportional-Integral-Differential (PID) algorithm is implemented to show the possibility to dynamically include in the simulation model the control of the gate opening during a flood. In particular, a specific maximum reservoir surface level is set as target in the automatic algorithm so that the gate opening must change under discharge variations during the flood event.

The PID controller computes the error between the predicted value of the variable water surface level and the stated reference value, and uses it to compute a change on the free parameter, gate opening, using an algorithm based on:


The equation that describes those PID terms is:

$$h\_{Crest}(t) = K \left\{ e(t) + \frac{1}{T\_i} \int\_0^{T\_i} e(t) \, dt + T\_d \, \frac{d\varepsilon(t)}{dt} \right\} = P + I + D \tag{21}$$

where *e*(*t*) is the control error (*e*(*t*) = *Href* − *H*(*t*)), *Href* is the reference value (or *setpoint*) of water level and *H*(*t*) is the current water level at time *t*. *K* is the proportional coefficient, *Ti* and *Td* are integration and derivative times, respectively.

Equation (21) is discretised as:

$$h\_{\rm Crest}(t^n) = a\_1 \, K \left( 1 + \frac{T\_s}{T\_i} + \frac{T\_d}{T\_s} \right) \varepsilon(t^n) - a\_2 \, K \left( 1 + \frac{2T\_d}{T\_s} \right) \varepsilon(t^{n-1}) + a\_3 \, K \, \frac{T\_d}{T\_s} \, \varepsilon(t^{n-2}) \tag{22}$$

where *e*(*t <sup>n</sup>*) = *Href*(*t <sup>n</sup>*) − *<sup>H</sup>*(*<sup>t</sup> <sup>n</sup>*), *e*(*t <sup>n</sup>*−1) = *Href*(*t <sup>n</sup>*−1) − *<sup>H</sup>*(*<sup>t</sup> <sup>n</sup>*−1) and *e*(*t <sup>n</sup>*−2) = *Href*(*t <sup>n</sup>*−2) − *H*(*t <sup>n</sup>*−2), being *n* the current time step. Parameters *α*1, *α*<sup>2</sup> and *α*<sup>3</sup> are the weights given to each of the time steps that are included on the controller operation. Parameter T*<sup>s</sup>* stands for the sampling period for the input data to the algorithm.

The values for *K*, *Ti*, *Td* and *Ts* directly affect the *hCrest* evolution and, thus, have an effect on the speed at which the controlled variable (i.e. water level) reaches the *setpoint*. Therefore, proper determination of those parameters is essential to optimize and stabilize the algorithm. Otherwise, the controller could lead to extreme gate opening values hence destabilizing the system.

#### **4. Model Application**

#### *4.1. Discretisation of the Domain*

The Ebro River reach has been first simulated to validate the 1D model comparing with results obtained with the 2D model [7] applied to the same stretch. In both cases the discretisation of the reservoir is included within the computational grid. The aim is to evaluate if the 1D approach provides reliable results improving the efficiency of the 2D model. Once the 1D model is demonstrated to be reliable enough, the coupled 1D-0D model and the control algorithm are evaluated.

The domain discretisation is different depending on the numerical scheme. In the 2D model, the mesh is an unstructured triangulation of the (x, y) domain with piecewise uniform values of terrain elevation and roughness. The triangles can be of variable size and adapt to the terrain topography. The mesh used was generated from a DTM in RASTER format and included 949,445 triangular elements. In the 1D model, the domain is discretised into a set of cells separated by cross sections along the riverbed evenly spaced at a distance Δ*x*.

The 1D mesh was generated from topographic information of the field compound by 433 cross sections. Among them, 100 sections are within the reservoir region (as in example in Figure 5). The lateral span of the sections must capture the shape of the river bed to avoid losing information relevant to the evolution of the variables but avoiding overlapping. It is of vital importance that the sections are always normal to the river in curved regions, as seen in Figure 5. Finally, a 1D mesh of 2000 cells is obtained. Additionally, a uniform roughness coefficient n = 0.032 is chosen ([33]). A steady flow with a discharge value that matches the value at the initial time of the inlet hydrograph is set as initial condition. The upstream boundary condition is a hydrograph, while the downstream boundary condition is a spillway condition, which represents the presence of the Mequinenza dam.

**Figure 5.** Top view of several sections over the raster with different river meanders.

#### *4.2. Performance Analysis of the 1D and 2D Models*

Two historical events of the Ebro River, the 2015 and the 2018 floods, have been simulated with both the 1D and the 2D models. The 2015 inlet hydrograph, obtained from the Zaragoza gauging station (see Figure 2), is set in Zaragoza as inlet boundary condition and can be seen in Figure 6. The comparison between results obtained with both simulation models and real observation data for this event can be seen in Figures 7 and 8. They show, respectively, the time evolution of discharge and water level at Gelsa (A263) station and Mequinenza dam (E003). It is important to note that the data provided by the CHE gauging station are for water depth (*h*) and not for water level (*H* = *h* + *zb*), and the actual bed elevation of the station is unknown.

It can be seen in Figures 7 and 8 that the 1D and 2D simulations produce remarkably similar data which follow the tendency of the actual data. However, neither of the two models is able to reproduce the detail of the curves at *t* = 380 h. This is possibly due to a dynamic change in the terrain such as a levee breach not included in the static terrain representation of the models owing to the lack of available information. It should be noted that, in the first 250 h, it is the 1D model which offers data closer to the reality. This suggests that the flow was more channeled in this period of time and the overflow is estimated by the 2D model too early. From *t* = 300 h it is the 2D model which provides a behaviour closer to the real one, possibly due to the fact that, near the peak flow, the floodplains are inundated.

**Figure 6.** Inlet hydrograph for the Ebro River flood event in 2015 in Zaragoza (A011).

**Figure 7.** *Cont*.

**Figure 7.** Discharge temporal evolution comparison between 1D model, 2D model and observation at Gelsa (A263) gauging station (**upper**) and comparison between models and estimations at Mequinenza dam (E003) (**lower**) for 2015 event.

**Figure 8.** Water level temporal evolution comparison between 1D model, 2D model and observation data at Gelsa (A263) gauging station for the 2015 event.

Figure 9 shows the 2018 inlet hydrograph used as upstream boundary condition. Figure 10 displays the discharge time evolution as computed with the 1D and the 2D models together with the observation at Gelsa (A263) gauging station (upper) and Mequinenza dam (E003) (lower) for this event. The time evolution of the surface water level at Villafranca gauging station (the only measured variable) is displayed in Figure 11. The evolution of the water levels predicted by the two models is again quite similar to that of the real data until approximately *t* = 100 h. Around this time, the measured data suffer a slight decrease not predicted by the 2D model.

**Figure 9.** Inlet hydrographs for the Ebro River flood event in 2018 in Zaragoza (A011).

**Figure 10.** *Cont*.

**Figure 10.** Discharge temporal evolution comparison between the 1D model, the 2D model and the observation at Gelsa (A263) gauging station (**upper**) and comparison between models and estimations in Mequinenza dam (E003) (**lower**) for the 2018 event.

**Figure 11.** Water level temporal evolution comparison between the 1D model, the 2D model and the observation in Gelsa (A263) gauging station for the 2018 event.

Regarding computational times, for the 2015 case the 1D model required 511 s to compute the full event, while the 2D model took 47 h, as seen in Table 1. These computations were performed with High Performance Computing (HPC) techniques for the 2D model. In particular, a NVIDIA GeForce GTX 1080 Ti GPU was used to compute the 2D cases, while a simple paralellization into 8 Intel Xeon X5650 CPU's was necessary for the 1D computations. For the 2018 event, which is a bit shorter, the computational time required by the 1D model was 364 s, while that of the 2D model was 23.8 h.


**Table 1.** Comparison of computational times for the two flood events simulated in the Ebro River.

#### *4.3. Performance Analysis of the 1D and 1D-0D Models*

The comparison between the full 1D model and the 1D-0D model is next carried out for the 2015 flood event. The discretisation of the full 1D model is the same used in the former subsection, while the 1D-0D model is characterised by a partial discretisation of the domain embedding the reservoir zone within the outlet boundary condition. In this last case, only the first 352 sections and 1511 mesh cells are necessary for the discretised part. This is, following Figure 4, for the river reach from x*<sup>L</sup>* = 0 to x*<sup>L</sup>* = L'.

Figure 12 represents the river profile for different times of both the full 1D model (fine and dark lines) and the 1D-0D model (wider and light lines). Bottom elevation, *zb*, as well as water level, *h* + *zb*, profiles are represented for both models. The initial condition can be seen at the upper picture of the figure, with a low water depth profile upstream the reservoir area, where the level remains uniform and the water depth increases. The middle picture corresponds to *t* = 300 h, when the discharge is increasing and a higher water depth can be seen. The lower picture coincides with the discharge peak of the inlet hydrograph (see Figure 6), reaching the highest value of water level. In the three cases, the level reached by the last cell of the 1D-0D model is almost the same as the value of the full 1D model, which discretises the whole reservoir. Therefore, it can be said that embedding the reservoir in the outlet boundary condition provides very similar results to those of a full model, without the necessity of such amount of computational cells.

As illustrated in Figure 12, the section located at x*<sup>L</sup>* = *L* is not exactly the beginning of the reservoir, as the length of the reservoir varies throughout the simulation depending on the level of the water surface. For that reason, this section is chosen displaced forward ensuring that it always belongs to the reservoir. Thus, for high level values, there is part of the reservoir being discretised and also simulated by the 1D numerical scheme.

The temporal evolution of the water level at x*<sup>L</sup>* = *L* can be seen in Figure 13. Although there is a very good agreement, the figure shows that the level of the 1D-0D model is slightly below the value of the full model at that location (x*<sup>L</sup>* = *L* ). This is because there is not a uniform level in the entire reservoir and the 1D model represents this behaviour. However, the 0D model assumes a constant level in the reservoir that matches quite exactly the level at the end of the reservoir in the 1D model (x*<sup>L</sup>* = *L*). In addition, the lag previously obtained by the model 1D-0D is no longer present. The reservoir surface area function, *S*(*H*), corresponds to the entire reservoir. However, as part of the reservoir is being discretised by the 1D model, the boundary condition of the 1D-0D model causes a slower evolution of the level, as it is considering that there is a larger modeled reservoir than there should be and, therefore, it overestimates the value of *S*(*H*).

The computational times of these simulations can be seen in Table 2. It can be seen how the 1D-0D model allows for considerable time reduction due to the due to the absence of the cells representing the reservoir. These results show that the coupled model results accurate enough to predict water levels along the river providing a performance improvement.

**Table 2.** Computational times for 2015 flood event simulated in the Ebro River with the 1D-0D model and the pure 1D model.


**Figure 12.** Longitudinal profile of bottom level, *z*, and water surface elevation (WSE) at different times computed with the 1D model and the coupled 1D-0D model for the 2015 case.

**Figure 13.** Comparison of computed water surface elevation (WSE) at x*L*=L' using 1D model and 1D-0D model and computed WSE at x*L*=L using 1D model.

#### *4.4. Performance Analysis of the 1D and 1D-0D Models Including DAM Regulation*

Once the coupled 1D-0D model has been proved to perform properly, a dam regulation algorithm is implemented in both the coupled and the full 1D model. This is, both models include in their outlet boundary condition a PID algorithm that updates the dam crest, *hcrest*, to approach or to maintain a reference water level or *setpoint*, regardless of the inlet discharge.

For that purpose, both models have been discretised exactly as in the preceding section. The parameters used in the PID controller must be first calculated and validated. In this work, the chosen values are *K* = 11576, *Ti* = 12 s, *Td* = 3 s, *Ts* = 1000 s, *α*<sup>1</sup> = 0.5, *α*<sup>2</sup> = 0.3 and *α*<sup>3</sup> = 0.2. Those values were obtained following [41]. The dam movement is limited by *vmax* = 0.01 m per time interval and the water surface level is limited by a maximum and minimum value of 115 m and 105 m respectively. The target water surface level is *Href* = 112 m. The comparison is done at x*<sup>L</sup>* = *L* (see Figure 4).

It is worth mentioning that variations of dam crest, *hcrest*, are a practical representation of variations in cross section area of spillways. In reality, dams can not change their crest, but the gate opening of their structure. However, the discharge law implemented would be the same and the dam crest results in a very representative parameter of the dam opening.

The Figure 14 shows the temporal variation of water level at x*<sup>L</sup>* = *L* for both models. Besides that, the figure also represents the time evolution of the dam crest throughout the simulation. At the same time, Figure 15 depicts, also for both models, the temporal evolution of outlet discharge. It worth mentioning that this discharge is the flow rate passing through the dam.

Figure 14 shows that, at the beginning of the simulation, the level of the reservoir is much lower than the *setpoint* (*Href* = 112 m), so the dam crest is at its maximum height preventing the water leaking (see Figure 15) and provoking an increase of water level. Once the reference is reached, the dam crest varies to maintain a constant water level while the inlet flow rate changes. At that time, the outflow hydrograph tends to resemble the inflow rate.

The time evolution of *hCrest* displayed in Figure 14 is rather similar for both models, reaching the target value in a short time. It is worth noting that the 1D model reaches the objective earlier than the 1D-0D model. This is provoked by the lower level obtained with

the 1D-0D model at x*L*=L' (see Figure 4) in comparison with that computed with the 1D model. This causes a delay in the reservoir filling up to the *setpoint*.

**Figure 14.** Temporal evolution of the water level and dam crest computed with the fully 1D model and the coupled 1D-0D model for the 2015 case with the control of a PID algorithm.

**Figure 15.** Temporal evolution of the outlet discharge at Mequinenza dam (E003) computed with the 1D model and the coupled 1D-0D model for the 2015 case with the control of a PID algorithm.

The computational times for the model with the PID algorithm are shown in Table 3. The same trend found for the comparison between the different approaches for the reservoir modellization results into an improvement of the optimization when using the 1D-0D model. Besides that, the PID algorithm does not penalize the computational time, but it makes the model more efficient.

**Table 3.** Computational times for 2015 flood event simulated in the Ebro River with the 1D-0D model and the pure 1D model, both with the dam regulated with a PID algorithm.


#### **5. Conclusions**

In this work, the performance of several modelling approaches has been compared in order to evaluate their results and computational requirements in a transient river flow event in a reach of the Ebro river (Spain) that includes a reservoir covering a large area. A 2D distributed shallow water model solved over a triangular grid and a 1D shallow water model have been used to discretise the full domain. Additionally, an aggregated volume balance model has been implemented to model the reservoir region in order to allow computational saving. This has led to a coupled 1D-0D model. Finally, a PID control algorithm has been implemented as a regulation technique at the dam location and has been combined with both the 1D model and the 1D-0D model.

From the comparison of the performance of the 2D and 1D models, it can be be concluded that the results of the 1D model for the recent flooding events at the considered Ebro River reach are very similar to those provided by the 2D model. The water level and discharge data predicted by both models follow the same trend. The cross sections used to build the 1D model computational mesh were carefully located to reproduce the river curvature in detail, which is important to obtain a realistic evolution of the hydraulic variables. This effort is justified by the immense computational saving that the use of the 1D model offers, as long as there is no interest in representing the floodplain flow, that the 1D model does not take into account.

The coupling of the 1D model for the river flow at the upstream reach and the 0D model for the reservoir (1D-0D model) offers results very similar to those from the full 1D model. There is some lag due to the instantaneous propagation of the hydrograph in the reservoir assumed by the 0D model but this is acceptable considering the computational savings that the use of this model implies compared to the full 1D model. The computational times observed with the 1D-0D model justifies the use of this combined approach. Therefore, the coupling of a 0D model for the reservoir with the 2D model for the upstream river reach is envisaged as future work since this will lead to high computational savings, something very positive for simulations with 2D models as well as the possibility to simulate the floodplain flow behaviour.

The PID control algorithm has been implemented with the objective to ensure a fixed surface water level at the dam. The results show that this target level value is never reached despite the time variable discharge, which means that the implementation of the control algorithm is a correct security measure to avoid exceeding certain levels in the reservoir. It will be convenient in the future to implement an algorithm that takes into account more realistic and complex objectives.

**Author Contributions:** Conceptualization, P.G.-N. and I.E.; methodology, P.G.-N. and I.E.; software, I.E.; 2D simulations J.M.; 1D model validation, P.V., I.E. and P.G.-N.; formal analysis, I.E., P.G.-N.; writing original draft preparation, P.V.; writing review and editing, I.E. P.G.-N.; visualization, P.V.; supervision, P.G.-N.; project administration, P.G.-N.; funding acquisition, P.G.-N. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was partially funded by the PGC2018-094341-B-I00 research project of the Ministry of Science and Innovation/FEDER. The authors would like to thank also the Confederación Hidrográfica del Ebro staff for their availability and for the supply of the data. Additionally, Isabel Echeverribar was wants to thank to the MINECO for his Research Grant DIN2018-010036.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** River measurements are available at http://www.saihebro.com/ saihebro/index.php?url=/datos/mapas/tipoestacion:A; accessed on 20 October 2021.

**Acknowledgments:** The authors acknowledge the CHE for the data availability and their support. The authors also would like to thank all collaborators for their help performing the 2D simulations: Pilar Brufau and Mario Morales.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **References**


MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Water* Editorial Office E-mail: water@mdpi.com www.mdpi.com/journal/water

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34 Fax: +41 61 302 89 18

www.mdpi.com

ISBN 978-3-0365-3318-6