2D Numerical Simulation of Floods in Ebro River and Analysis of Boundary Conditions to Model the Mequinenza Reservoir Dam

: The computational simulation of rivers is a useful tool that can be applied in a wide range of situations from providing real time alerts to the design of future mitigation plans. However, for all the applications, there are two important requirements when modeling river behavior: accuracy and reasonable computational times. This target has led to recent developments in numerical models based on the full two-dimensional (2D) shallow water equations (SWE). This work presents a GPU accelerated 2D SW model for the simulation of ﬂood events in real time. It is based on a well-balanced explicit ﬁrst-order ﬁnite volume scheme able to run over dry beds without the numerical instabilities that are likely to occur when used in complex topography. The model is applied to reproduce a real event in the reach of the Ebro River (Spain) with a downstream reservoir, in which a study of the most appropriate boundary condition (BC) for modeling of the dam is assessed (time-dependent level condition and weir condition). The whole creation of the model is detailed in terms of mesh optimization and validation. The simulation results are compared with ﬁeld data over the ﬂood duration (up to 20 days), allowing an analysis of the performance and time saved by different GPU devices and with the different BCs. The high values of ﬁt between observed and simulated results, as well as the computational times achieved, are encouraging to propose the use of the model as a forecasting system.


Introduction
For centuries, natural disasters have been a source of concern for human beings due to the damage and losses they cause. In addition, in recent years, these losses and their frequency have been on the rise [1], leading to increased concern from governments, institutions and society in general. Within natural disasters, floods are one of the most destructive extreme events [2] and are the second leading cause of natural-disaster-related deaths in Spain, with 209 deaths between 2000 and 2019 [3]. They also entail a high expense due to damage repair, with losses in Spain amounting to EUR 12,000 million in the period between 2016 and 2020, with Zaragoza-located at NE of Spain-being the province with the third highest economic damage in the agricultural sector in this period. Figure 1 shows examples of such losses, demonstrating images of the surroundings of Zaragoza during flooding of the Ebro River, which is the largest river in terms of discharge in Spain. In view of these numbers, governments and public institutions require tools [4][5][6][7][8][9] and plans [10][11][12] to foresee and mitigate the damage caused by these events. One of these tools is the use of predictive models based on numerical simulations which are able to provide an accurate description of the spatial and temporal evolution of flow [13][14][15][16][17][18][19]. Due to the physical complexity of these phenomena, it is common to consider approximations to simplify the equations describing the flow. Although a river flood event is naturally a 3D problem, it is common to study it by averaging the equations in the vertical coordinate to reduce the problem to two horizontal dimensions [13,20]. Practical applications require a trade-off between spatial accuracy and computational efficiency [21], so approximations that reduce the dimensions of the problem are frequently used. The shallow water model (SWE) is a widely used approach to simulate surface geophysical flows in situations that involve large domains and long time scales [22][23][24]. This approach is based on the assumption that horizontal scales are larger than vertical scales, which leads to the possibility of neglecting vertical accelerations and assumes a hydrostatic pressure distribution in the vertical direction. This is the basis of the dynamic non-linear 2D SWE formulation. In some cases, the full dynamic equations are reduced in complexity by neglecting inertial terms. Reductions of the SWE system can be used to model floods with zero inertia models, which maintain the 2D framework and neglect terms that do not govern the phenomena, although special attention must be paid to their limits [25]. Further dimensional simplifications consider the average of the equations in the cross-section to reduce the formulation to a 1D approximation [14,26,27]. In large and complex flow domains, as is the case of a river in a flood event, two-dimensional models are the most frequently used to obtain the temporal and spatial description of the flow [13,20,21,25,[28][29][30][31][32][33].
When modeling hydraulic structures present in the domain, the adopted different numerical strategies are of wide variety not only when these structures govern all the flow [34], but also if they only affect a part of the flow. For instance, when modeling complete reservoirs affecting a part of river, several strategies can be considered. On one hand, as the reservoir dynamics are close to quiescent equilibrium, some aggregated dimensionless models can be used to model their presence [35][36][37][38][39]. When more detail is required, the reservoir can be discretized [40] as a river extension by incorporating it to the computational domain and the dam presence is introduced via boundary conditions with a different hypothesis, as applied in this study.
To represent the behavior of a dam spillways, several boundary conditions can be used. If the main objective is to model the backwater generated by the dam, a constant water level with the main surface elevation of the reservoir can provide good results in terms of discharge and modeling the river. However, if the level of the reservoir is a variable to compute, which is a requirement of many applications focused on dam regulation, other, more sophisticated BCs must be applied. In particular, when the physical characteristics of the dam spillway are known, the outflow provided by the general discharge law of a dam spillway can be imposed as a boundary condition, so the backwater is still generated but the level is allowed to vary. Both strategies are analyzed in the this study.
Therefore, the main aim of this work is to study the region of the middle reach of the Ebro river between Zaragoza and Mequinenza, both located in Aragón (Spain), as seen in Figure 2. This region is not only of special interest due to the significant damage caused by large floods that occur in the meandering flood prone areas in the first half of the river in flood events, as seen in Figure 1a, but also because it contains a long reservoir, the Mequinenza reservoir, limited by the Mequinenza dam downstream. In this work, a 2D model is used for the discretization of both the river and the reservoir in order to compromise between computational efficiency and accuracy of the results. Therefore, the main objectives are to obtain an accurate computational model setup to study the Ebro River region and to optimize it, obtaining a predictive tool to foresee and mitigate potential damage caused by flooding events. Moreover, the presence of a reservoir allows the study of different boundary conditions that model dam spillways in order to obtain a realistic temporal evolution of the surface water level in the reservoir. For the simulation, the PEKA2D program [16], developed at the University of Zaragoza, is applied to the mentioned river reach. A version of this program is included as computational core of the brand name RiverFlow2D ® (Hydronia LLC, https://www.hydronia.com/, accessed on 1 March 2023).

Governing Equations and Numerical Model
The governing equations and the numerical aspects of the scheme used can be found in full detail in [41,42] and it has been extended to 2D unstructured meshes in [20]. Their most important details are outlined in the following subsections.

2D Shallow Water Equations
The mathematical model that describes the surface flow is given by the hyperbolic 2D shallow water system of equations based on mass and momentum conservation [43]: where the conserved variables: are the water depth, h, and the unit discharge in x and y direction, q x = hu and q y = hv, respectively, with (u, v) being the depth averaged components of the velocity. The fluxes of these conserved variables are The source term for the mass conservation equation is zero because neither precipitation, infiltration nor evaporation are included, assuming their contribution is practically negligible during a flooding event. Finally, the momentum source terms are related to the bed slopes and friction stresses: The bed slopes represent the variation in the x and y directions of the bottom level, z b : (5) and the friction stress components are given by: where n is Manning's roughness coefficient [44].

Numerical Scheme
The system in Equation (1) must be solved numerically due to the lack of an analytical solution. For this purpose, an explicit upwind finite volume scheme, based on the Roe-Riemann solver [45,46], is used in this case. From (1), in compact form, this can be expressed as: ∂U ∂t where E = F, G . Integrating (7) in a control volume or cell, Ω, and applying the divergence theorem to the second term, we obtain: where ∂Ω is the contour of the control volume andn is the outgoing unit vector normal to the Ω volume. By discretizing (8) in time and space, the basis of the numerical method in finite volumes is given by: where Ω i is the area of cell i, n is the current time level and the number of neighboring cells is 3 because triangular cells are used in this work (see, for example, Figure 3). In addition, the fluxes E are evaluated at cell boundaries: where E j is the flux value at cell Ω j , and shares a wall k of length l k with cell Ω i with a flux value of E i . Considering the hyperbolic character of system (7), the Jacobian matrix normal to the flow direction, E, can be defined as: The local value of the Jacobian matrix (11),J n k , at wall k is whereΛ k is the diagonal matrix whose non-zero elements are the eigenvalues of the system λ m , andP m is the matrix containing the eigenvectors of the systemẽ m , providing a matrix with three eigenvalues to the 2D model. The eigenvalues and eigenvectors of the Jacobian matrix are:λ 1 =ũ ·n −c ,λ 2 =ũ ·n ,λ 3 =ũ ·n +c (13) whereũ ·n =ũ n x +ṽ n y andc is the celerity of the infinitesimal surface deformation waves. The tilde variables represent an average state at each cell edge. Therefore, starting from the expression (9), and using the eigenvalues and eigenvectors of the Jacobian matrix, an updated expression at cell i is obtained: The time step ∆t is calculated dynamically throughout the simulation by: with 0 < CFL ≤ 1 to guarantee stability in the numerical scheme, where CFL is the Courant-Friedrichs-Lewy number [47] and where Equation (15) is only solved in those cells where water is present and the flooded boundary is not regulated but advances according to the flow dynamics. The method is well designed to deal with this type of calculation ensuring accurate and robust results.
The time step depends, on one hand, on the dynamics of the problem to be solved throughλ m k and, on the other hand, on the cell size chosen for the computational grid, given by l k . Therefore, from (16), (17) and (18), it becomes clear that increasing the refinement of the computational grid leads to smaller time steps, and therefore a higher computational cost. For this reason, an unstructured mesh is used, in order not to restrict either the accuracy of the results with a very coarse mesh over the whole domain or the time step ∆t with a very fine mesh [28].

Study Case
The studied reach of the Ebro River is located between the city of Zaragoza and the Mequinenza Dam. It is more than 200 km long and covers 722 km 2 of surface area. Along the river there are a few gauging stations managed by the Ebro River Authority (CHE, www.chebro.es, accessed on 1 March 2023), where the evolution of the flow and surface level are continuously (fortnightly) recorded. Their locations are represented in Figure 4. The area is continuously suffering from the damage provoked by these events, which have even flooded the A-1107 motorway and collapsed the regional highway (ARA-1) in Villafranca de Ebro (www.heraldo.es/noticias/aragon/2018/04/12/crecidas-del-ebro-lasultimas-riadas-aragon-1234800-300.html, accessed on 1 March 2023) during the 2015 event (see Figure 1). In addition, it is an area of great importance due to the presence of the Mequinenza reservoir, which runs along about 75 km of the river, and the regulation of its dam. The first 25 km of the studied reach is characterized by meanders and associated with important inundation areas. The final part, on the other hand, is characterized by the reservoir vessel where water is practically at rest.
The Mequinenza reservoir, whose satellite view can be seen in Figure 5, covers a surface area of approximately 7540 hectares, with a maximum capacity of 1530 hm 3 when the surface level is 121 m.
The Mequinenza dam (see Figure 6) has a crest height of 124 m and a single spillway with six gates located at an elevation of 106.5 m whose discharge limit is 11,000 m 3 /s. Data on the reservoir and dam were provided by the Ebro River Authority. The gates at the spillway are used to regulate the volume of water in the reservoir for various purposes, but mainly to control floods and guarantee hydroelectric generation.

Computational Model Setup
The configuration of the computational model is based on the following parts: a surface roughness map; 3.
The creation of a triangular mesh; 4.
Boundary conditions and initial conditions.

Topography: DTM
A raster digital terrain model with a resolution of 5 × 5 m provided by the IGN (http://centrodedescargas.cnig.es/CentroDescargas/locale?request_locale=en, accessed on 1 March 2023) and obtained by interpolating data from flights with LIDAR sensors in 2010 was used as a base. However, the used data do not contain a reliable representation of the riverbed. Thus, this DTM is used only for floodplains, where a proper representation of the terrain can be found, as seen Figure 7, while complementary information from measured cross-sections is used to reconstruct the river bed, since the LIDAR provides a uniform free surface at the river, as seen in Figure 7.

Reconstruction of the Riverbed from Sástago
To achieve a reliable riverbed, cross-sections of the river are taken at the end of the Osera-Sástago DTM. These cross-sections, which are groups of coordinates (x, y, z), are duplicated along the riverbed up to the Chiprana area. Afterwards, the elevation of the copied sections is corrected, since the terrain through which the river flows descends in altitude. To do this, two points are taken, one at the beginning of the part to be reconstructed and the other at the end. The difference between their elevations is calculated, and the distance that separates them is a straight line. The quotient of these two values is an average slope in the river segment, so that the elevation of each of the sections can be recalculated: where the index i denotes the section number and the index j denotes the point within each section. m denotes the average slope. This methodology is followed from Sástago to the Mequinenza reservoir, where, since there are no sections to start from, it cannot be used.

Reconstruction of the Reservoir
River cross-sections cannot be used to reconstruct the bottom of the reservoir. Therefore, to incorporate that information to the global DTM, geo-referenced historical maps prior to the construction of the dam, also available on the IGN website (http://centrodedescargas. cnig.es/CentroDescargas/locale?request_locale=en, accessed on 1 March 2023) were used. Using their contour lines, new cross-sections are obtained formed by groups of five points with coordinates (x, y, z), as shown in the example in Figure 8, which, after interpolation, produces a DTM with the appropriate reservoir bed elevations. Figure 9 shows an image of the historical map of part of the reservoir, comparing it with the current state. The historical topographic map clearly shows how the Ebro riverbed ran under what is now covered by the reservoir. The center points of each section, which were not placed on grade lines, were corrected using Equation (19). The lowest elevation of the reservoir bottom, z = 60 m, was taken in the first section next to the dam, and from there, the minimum elevation of each section was raised. An example of the result is shown in Figure 10. Using the presented strategy, a channel in the reservoir region similar to the one shown in Figure 11 is obtained, containing the information of the interpolated sections shown in Figure 8. By doing this, a new DTM fo the bottom of the reservoir is obtained and can be added to the global DTM.

Surface Roughness Map
Each computational grid cell is associated with a value of the Manning roughness coefficient, n, which is assigned from a terrain roughness raster map provided by CHE. This map covers the potential floodplain from Zaragoza to Escatrón, where the reservoir starts and the dynamics are not as affected by roughness. Thus, the region covered by the Mequinenza reservoir has a simpler roughness distribution reconstructed following the criteria in the literature [44,48,49]. This distribution is shown in Table 1, and the whole distribution map can be seen in Figure 12.

Meshing
The numerical simulation is based on the discretization of the terrain into cells that form the so-called computational grid. Relevant magnitudes are associated with each of its cells during the simulation, such as bed elevation, water depth or velocity. The computational mesh has a triangular geometry and is spatially adaptive, i.e., the cell size varies in space. The riverbed, the adjacent levees and, in general, all the areas to be represented in detail require a much finer mesh. Less relevant areas or areas with no abrupt changes in elevation can be discretized with larger cell sizes. Figure 13 shows an example of this, where it can be seen that the cells close to the channel have smaller sizes than those of the fields with more or less uniform elevation. In addition, a greater refinement of the mesh is also observed at the levees.

Initial Condition
In river problems, the initial and boundary conditions are of great importance. In order to perform a flooding event simulation, the initial condition must correspond to the steady state state of the river flow before the flood event occurs. For this reason, the initial condition comes from the convergence of the model to a steady flow condition, where the entire river flows with a constant flow of the same value as the discharge at the initial instant of the flood.

Boundary Conditions
The boundary conditions provide information about how the flow enters and leaves the domain. For the river inflow, the commonly used rating curve boundary condition is applied, i.e., a time variation in discharge at the inlet section of the reach, provided by the gauging station. On the other hand, for the downstream boundary conditions, river simulations typically use a gauging curve. However, the domain under study ends at the Mequinenza dam, so it is necessary to study the most appropriate outflow boundary conditions.

Upstream Boundary Conditions
In the present case, there are inlet sections: the main channel of the Ebro river and the mouths of two tributary rivers (the Gallego River and the Huerva River). A hydrographtype boundary condition is imposed for all of them. The temporal evolution of the flow imposed in each of the regions is obtained from data provided by CHE at gauging points.

Downstream Boundary Conditions
The boundary conditions considered to model the behavior of the dam are detailed below.

Time-Dependent Level
Considering that the flow is almost at rest when it reaches the reservoir, it can be considered that the level in this region is practically uniform and constant [50]. The way to use this boundary condition in our case will be to impose a constant downstream level throughout the simulation, the value of which will depend on the dimensions of the flood event under study.

Dam Spillway Condition
A weir boundary condition is implemented in this work as a function, Q out = f (H), to model the dam outflow [51]. Considering certain simplifications, the outflow of a trapezoidal weir follows the expression: where H w = H − h Crest is the thickness of the sheet of water above the weir crest h Crest ; H is the level of the free surface (H = h + z b ); α is twice the angle that the lateral sides of the trapezoid make with the vertical; b is the width of the minor base of the trapezoid (see Figure 14); and C d = 0.611 [52]. Taking into account the shape of the gates of the Mequinenza dam (see Figure 6), the latter case will be imposed on the output of our problem. In this work, the spillway law is applied to the dam boundary condition as a generic discharge law. Depending on the value of H w , different outflows through the spillway will be obtained. When the level in the cross-section in contact with the weir is below the crest height, there will be a zero outflow (H w ≤ 0), while if the level is above (H w > 0) the crest height, the outflow will follow the expression (20). Thus, the flow function Q out = f (H) through the weir is defined as follows:

Calibration and Optimization of the Computational Mesh
Although the numerical model used has been validated on numerous occasions and in different domains [20,53], each model configuration requires a calibration process to adjust parameters and correct data that may contain errors. Moreover, even if the model has been calibrated and validated, this does not mean that it is optimal, leading to high computational consumption if the discretization is excessively fine. Therefore, the configuration of the model must be optimized to reduce computational consumption as much as possible.

Calibration of Mesh Refinement
In this work, an unstructured triangular mesh is used with the advantage of being able to adapt to the topography. This mesh contains smaller cell sizes in those areas with a need for detail, while coarser cells are used in areas where there are few or no relevant terrain irregularities.
Each floodplain has its own topographical particularities that may be relevant not only because of their geometrical characteristics, but also because of the effect they usually have on flooding. There are always overtopped river banks in the first instants of a flood, the modeling of which may have a lower impact than other more distant levees that retain water for long periods of time until they are over passed by water. The detection of these structures is crucial for a correct modeling of the terrain through refinement.
In order to choose the most appropriate mesh, satellite observations of the flood extension are used. By comparing the flooded areas provided by the different computational meshes, sensible zones where refinement was needed were detected. Thus, the cell size distribution of the preliminary mesh, M1, was modified, leading to a refined mesh, M2. The differences in resolution between the meshes are shown in Figure 15, where it can be seen that in M1, some of the floodplain levees are not well represented. This results in different flood extensions that can be seen in Figure 16, where the water depth for a certain time is represented. Finally, in the same figure, those areas provided by the two different meshes are compared to the observed flood envelope points. The preliminary mesh does not retain the water as it should because of the coarse representation of the floodplain irregularities.
To carry out this analysis, the 2018 discharge evolution was used as the inlet boundary condition. As the reservoir has no influence in the flood prone area upstream, a simple constant level is set to the reservoir as the downstream boundary condition.  In Table 2, some data related to the meshes and their performance can be seen. In view of these results, it can be concluded that, although M2 provides more accurate results, its computational cost is unfeasible. For this reason, a second optimization step was carried out and is detailed in next subsection.

Optimization
When the mesh is calibrated to refine the results and provide more accurate results, optimization should be performed to reduce the computational consumption without reducing the accuracy of the results. In this work, the optimization of the model was carried out by using internal boundary conditions that allow modeling some hydraulic structures of the channel, such as dikes and motes, without the need for very fine computational meshes [34,54]. In this way, a mesh with a smaller number of cells, M3 mesh, was obtained, with a consequent reduction in computational consumption. The validation of the optimization of the model is carried out by simulating the historical flood event of 2018 and the comparison in the Gelsa gauging station (A263). The results of water level and discharge temporal evolution are shown in Figures 17 and 18, where it can be seen that M1 is much less accurate than M2 and M3, and these two are similar, although M3 is somewhat less accurate than M2. However, this loss of accuracy is justified by the reduction in computational consumption obtained, going from 23.80 h with mesh M2 to 12.48 h with mesh M3, as can be seen in Table 3. Moreover, Table 4 shows the Root Mean Square Error (RMSE) of the discharge values to demonstrate how accurate each mesh is. The Root Mean Square Error is calculated by the following equation [55]: where x f ,n is the simulated value at time n, x o,n is the measured value at time n and N is the number of measures at a point. Table 4 shows that the accuracy of M1 is lower than M2 and M3, and these two offer a very similar accuracy.  Finally, to provide information on the cell size distribution of the meshes, Figure 19 shows the cell area distribution for each of the analysed meshes (M1, M2 and M3). As the three meshes share the maximum cell size, the x-axis of the figure is bound to 1000 m 2 , with the purpose of focusing on the refinement distribution. The M1 mesh contains less refined elements than the other two, with only small cells in the river bed area. The most refined mesh, M2, has a higher concentration of small cells, as the refinement affects not only the river bed but also the levees in the floodplains. Finally, the figure shows how the optimized mesh, M3, is in the middle of the two previously mentioned distributions. The use of an internal boundary condition allows a coarser mesh in the floodplains without loosing accuracy.

Numerical Results
A relevant historical event of flooding of the Ebro River which occurred in 2018 was simulated with both analysed downstream boundary conditions: the constant level and the dam model.

Event
The 2018 inlet hydrograph (see Figure 20), obtained from the Zaragoza gauging station (A011), is set as the upstream boundary condition. As an initial condition, a steady flow with a discharge value that coincides with the initial discharge of the inflow hydrograph is established. A comparison between numerical results and real observation was performed using data from the Gelsa gauging station (A263) (Figures 21 and 22) and the Mequinenza gauging station (E003) (Figures 23 and 24). The results show that at the Gelsa gauging station, both boundary conditions provide the same results since the reservoir does not affect upstream areas such as Gelsa (see Figure 4), as the flow regime is very different from the reservoir in the first half of the studied reach, as can be seen in Figures 21 and 22. At the Mequinenza gauging station, the constant level boundary condition provides a discharge temporal evolution more similar to the actual data than that given by the weir boundary condition (see Figure 23). However, the water level temporal evolution provided by the weir boundary condition is much more realistic than the results from the constant level boundary condition, which cannot represent the temporal evolution of this value, as can be seen in Figure 24.

Conclusions
In this work, a finite volume numerical model, well designed for the resolution of unsteady 2D SWE on flexible and adpatative triangular meshes, is used for the simulation of flood events in a specific region of the Ebro River (Spain). This region, which contains a very unique reservoir that changes the river dynamics, presents a challenge for the model regarding dam modeling. For this reason, this work has been able not only to prove the good performance of the model in a particular region, but also to explore new strategies for dam simulations.
First, this study highlights the importance of the choice of computational model. This implies the calibration and optimization of the computational mesh that is first transformed from a coarse to a refined mesh, showing the importance of adaptive meshes and the proper refinement at relevant topography points. From this calibration, it can be concluded that the discretization in certain areas, such as the limits of riverbeds, is crucial to obtain accurate results. However, this produces a very fine mesh with a high number of cells at a high computational cost.
The application to the simulation of flood events in a reach of the Ebro river highlights that the use of real data introduces some uncertainties related with coarse discretization of the terrain measurements, non-detailed characterization of bed roughness, spurious points on the discharge time series and other problems that may provoke errors in the results regardless of the numerical scheme. Hence, a calibration processes must be carried out. An optimal computational mesh has been generated and calibrated for the 2018 flooding event in the Ebro River. Additionally, although real test cases introduce some errors due to the lack of available data, the model is still able to provide very good results. Finally, the benefits of an accurate and fast numerical method are not only desirable for flood prediction, but also the generation of an appropriate computational mesh involving internal boundary conditions. This strategy provides results with enough accuracy, allowing the use of a coarser mesh where the small loss of accuracy obtained is justified by the 12 h of computational savings. These strategies, together with an adequate representation of land use maps, are demonstrated to be necessary in order to carry out computations leading to accurate numerical results.
The historical event was simulated using two different downstream boundary conditions at the Mequinenza dam location. The predictions of discharge and water level at the gauging locations improve when the outlet boundary condition is expressed as a discharge rating curve. This work has highlighted the need for the modeling of hydraulic structures integrated within simulation models instead of simpler boundary conditions. It is worth noting that the downstream boundary condition has no influence on the first half of the river reach, where the cross-sections are shallow and well connected to the floodplain; however, in the second half of the reach, dominated by the long reservoir, a weir-type boundary condition should be used if the temporal evolution of the level is of interest, thus obtaining more realistic results. At this point, it should be noted that this work highlights this need, but opens it up to improvements in terms of dam modeling, where more complex discharge laws could be used and combined to represent bottom spillways, gates and other hydraulic structures.
The numerical results obtained have been presented and compared with measurements. It should be emphasized that the use of HPC technologies is vital when carrying out simulations with large domains and long event durations. In this work, the results for the Ebro River, containing a large domain and great number of computational cells, were obtained using a GPU-parallelized well-balanced upwind numerical scheme which simulates a hydrograph of 14.5 days in a little more than 12 h with a GPU NVIDIA GeForce RTX 3070. Thus, it has been made feasible to reproduce such events on a real-time basis.  Data Availability Statement: River measurements are available at http://www.saihebro.com/saihebro/ index.php?url=/datos/mapas/tipoestacion:A; accessed on 1 February 2023. 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.

Abbreviations
The following abbreviations are used in this manuscript: