Surface Water Flow Balance of a River Basin Using a Shallow Water Approach and GPU Parallel Computing—Pescara River (Italy) as Test Case

Abstract

The analysis and prevention of hydrogeological risks plays a very important role and, currently, much attention is paid to advanced numerical models that correspond more to physical reality and whose aim is to reproduce complex environmental phenomena even for long times and on large spatial scales. Within this context, the feasibility of performing an effective balance of surface water flow relating to several months was explored, based on accurate hydraulic and mathematical-numerical models applied to a system at the scale of a hydrographic basin. To pursue this target, a 2D Riemann–Godunov shallow-water approach, solved in parallel on a graphical processing unit (GPU), able to drastically reduce calculation time, and implemented into the RiverFlow2D code (2017 version), was selected. Infiltration and evapotranspiration were included but in a simplified way, in order to face the calibration and validation simulations and because, despite the parallel approach, it is very demanding even for the computer time requirement. As a test case the Pescara river basin, located in Abruzzo, Central Italy, covering an area of 813 km² and well representative of a typical medium-sized basin, was selected. The topography was described by a 10 × 10 m digital terrain model (DTM), covered by about 1,700,000 triangular elements, equipped with 11 rain gauges, distributed over the entire area, with some hydrometers and some fluviometric stations. Calibration, and validation were performed considering the flow data measured at a station located in close proximity to the mouth of the river. The comparison between the numerical and measured data, and also from a statistical point of view, was quite satisfactory. A further important outcome was the capability to highlight any differences between the numerical flow-rate balance carried out on the basis of the contributions of all known sources and the values actually measured. This characteristic of the applied modeling allows better calibration and verification not only of the effectiveness of much more simplified approaches, but also the entire network of measurement stations and could suggest the need for a more in-depth exploration of the territory in question. It would also enable the eventual identification of further hidden supplies of water inventory from underground sources and, accordingly, to enlarge the hydrographic and hydrogeological border of the basin under study. Moreover, the parallel computing platform would also allow the development of effective early warning systems, for example, of floods.

1. Introduction

The appropriate selection of an effective calculation tool constitutes the first action to be explored and pursued in order to obtain, as far as possible, a correct forecast of the evolution of possible impacting natural phenomena for the territory under observation. Several simplified methods have been proposed in the literature to evaluate the annual balances of water flows in river basins [1]. However, the approaches commonly used are based on the elaboration of experimental data expressed through empirical formulas. Obviously, the feasibility of using more advanced mathematical-numerical tools would be desirable. The most advanced approaches, based on the balance equations of physics are pursued by computational fluid dynamics (CFD). CFD relies on solving Navier–Stokes flow equations, including turbulence [2,3,4,5], the results of which are, consequently, very complex from a mathematical point of view. Therefore, to avoid at least the further complexity of the geometric domain discretization, the meshless approach has been proposed and adopted. In this context, the smoothed particle hydrodynamics (SPH) method appears to be quite well established [6,7,8].The mistaken was that Reduced complexity models (RCM), to which, for example, the cellular automata (CA) approach belongs [9,10], represent an important alternative to CFD, suitable to be applied on a large area and over relevant time scales (climate evolution as well), both at reach and at catchment scale. Therefore, for the study of the phenomena considered here, a suitable compromise can be found in approaches capable of simplifying the CFD technique, such as the widely used, simplified but effective ‘Shallow Water’ approach [11,12,13,14,15], as implemented in the RiverFlow2D calculation code [16] that was chosen to perform the simulations discussed in the following sections. The use of the GPU and, therefore, of parallel computing even on simple workstations and laptops, (by consequence not necessarily on very less accessible hardware such as supercomputers) was fundamental to fill the gap between two conflicting requirements. That is, the need for reasonable computation time and the need to use more realistic and advanced mathematical-numerical models than semiempirical or very simplified models such as CA, as well as for the study of phenomenology that evolve over a period of one year on a spatial scale of a river basin, similar to the selected test case. The results discussed in this article, which can be improved in the extension of this type of study and with more powerful but still accessible hardware, have demonstrated the concrete possibility of carrying out these types of simulations. The options available in RiverFlow2D allow the inclusion of various topic models that are useful for civil protection, from debris flow simulations to urban floods and, through the option of parallel computing accessible on desktop and laptop (GPU), also for a possible early-warning system [17,18,19].

The objective of the study discussed in this paper was to model the surface water flow dynamics of the Pescara River Basin, located in Abruzzo (Central Italy) and selected as a well representative test case of medium-sized fluvial basins. The territory is equipped with experimental measurement stations and has been studied over the years through several geological campaigns which have produced numerous available data. Detailed parametric sensitivity approaches, regarding some leading parameters that are difficult to be directly measured [20], were employed. The related assumed values, although reasonably justifiable from the physical point of view, were essentially used to calibrate the simulations. First, an extensive bibliographic analysis was conducted to have an appropriate geological knowledge of the study area [21] while also consulting the data relating to the hydrogeological balance [22,23,24,25]; something which is necessary to quantify the dynamics of surface waters within the basin. The available stations network consists of 11 pluviometers distributed over the entire basin area, whose data are then processed by using the Horton’s semiempirical model [26], whose parameter values in turn were selected from the literature, based on soil grain size, morphology and incidence of vegetation [27,28]. Evapotranspiration was also included, but in a simplified way, as detailed in Section 2.3.3., to reduce the number of parameters to be calibrated and, consequently, to more easily face the problem of computation time, which was very demanding despite the parallel approach. The inflow to the Pescara River basin was included considering the monitored data of the measurement stations, especially those at Pescara-Maraone and Tirino-Madonnina as main contributions. Minor water inventory sources, poorly defined not only during the reference periods for calibration and validation, but also currently, were not included. For example, the impact of snow occurrence and melting modalities (temperatures rising, rain events, etc.) [29] were not explicitly considered. However, in the calibration step, the not completely known contributions were included by adjusting some parameters, as detailed below. Accordingly, only reliable contributions were considered. In a subsequent phase, using the QGIS software, the morphological digital model of the basin under study, based on a 10 × 10 m DTM, was processed. ‘Shape files’ were created containing information relating to the area such as the perimeter of the basin area and the location of the 11 pluviometric stations defining their influence area by the Voronoy polygon technique. Subsequently these ‘shape files’ were inserted and processed on the River Flow 2D dedicated graphic user interface SMS, thus proceeding with the construction of the conceptual model, containing all the necessary information (perimeter of the study area, flow boundary conditions, Manning’s [30] roughness coefficient, rainfall distribution) with the related infiltration and evaporation rate as well. Once the conceptual model was acquired, the Delanuay [31,32] numerical tessellation of the selected DTM was performed through an unstructured mesh approach, based on linear triangular elements, applying the meshing tool accompanying the Riverflow2D computer code.

To satisfy two conflicting requirements, namely the need for both numerical accuracy (high number of tessellation elements) and a less excessive execution time, the identification of an appropriate, balanced numerical discretization of the spatial domain was essential. Subsequently, after an accurate analysis of the regional annals, some complete time series, whose experimental data were consistent with each other, were identified. Therefore, once the mesh was generated, the calibration process was carried out on the basis of the experimental data of one month, measured in January 2003, while three semesters, two belonging to the full year 2000 and the third belonging to the first half of the year 2007 were selected for validation. In this work, in particular, we analyzed the impact that the infiltration-evapotranspiration rate and the Manning number have on the numerical results of the selected model and how the computation time can be reduced through the use of the GPU.

2. Materials and Methods

2.1. Geographical and Geological Setting of the Selected Study Area (Pescara River Basin, Abruzzo, Italy)

The Pescara River basin is in the foothills of the Abruzzo Region, extending from the eastern side of the mountain range (Massiccio della Majella, Monte Morrone and part of the Gran Sasso) to the coast of the Adriatic Sea, for a total length of about 45 km. This river is fed mainly by the Aterno and Sagittario rivers (whose contributions were included as flow-rate inputs measured at Tirino-Madonnina and Pescara-Maraone stations, Figure 1), as well as by various sources located on the right side of the river shaft. The investigated area of the hydrographic basin is about 813 km². The river enters the narrow gorges of Popoli considerably swollen, squeezed by the mountains of the Gran Sasso to the north and the Morrone massif to the south. Figure 1 shows the hydrographic boundary, selected for the study discussed in this paper. Maintaining the north-east direction, the river then crosses the alluvial plain of Pescara Valley, laps the lower part of the city of Chieti and enters to Pescara, dividing it in two and finally flowing into the Adriatic Sea. In this context, the mountain range has a high relief characterized by resistant preorogenic carbonates (Mesozoic—Lower Tertiary) and synorogenic siliciclastic valleys (Tertiary), as well as by intramountain basins of continental Quaternary deposits (post-orogenic) [23]. The foothill area is to the east of the mountain range and looks like a hilly environment, characterized by Miocene, Pliocene and Quaternary terrigenous deposits (related to sin- and post-orogenic phases of the Apennines). Here the marine environment persisted until the lower Pleistocene. Later, during the Middle Pleistocene, a phase of tectonic uplift began which, together with a series of glacial-interglacial cycles linked to the climatic transition of the period, favored a phase of fluvial incision and deposition. In general, however, we can distinguish the area into two fundamental units: the chain and the foredeep [24,25]. The first is characterized by carbonate masses, mostly Mesozoic, with variable facies, and by deposits of the Messinian terrigenous facies. The deposits of the carbonate platform, generally calcareous, are quite homogeneous. Instead, the Messinian facies (ranging from arenaceous, pelitic and evaporitic lithotypes) and pelagic (calcareous, flint, marl and clayey lithotypes) present an important diversification. The Avanfossa area, on the other hand, is characterized by clayey, marly and sandy sediments relating to the Plio-Pleistocene. The tectonic style of the area in question is characterized by imbricated tectonic scales in eastern vergence, with evident overlaps on the surface in correspondence with the chain, but also in depth where deposits dating back to the Middle Pliocene are involved. In this scenario, the Majella Massif is of particular importance as it is characterized by one of the water reservoirs of the Abruzzo region [23], carrying the potable water supply of the region, and has therefore been the subject of many studies for evaluation and protection. The Majella has a wide asymmetrical anticline, with an eastern slant that presents a variable inclination, which grows up to be sub-vertical. This anticline is characterized by a rather wide hinge, with an axis that progressively rotates towards the South, presenting NW-SE orientations in the northern part, to NS in the central Majella, and to NNE-SSW in the southern part. Observing the area to the west, from Caramanico to the south, the Majella anticline is cut off by a fault called “Faglia della Majella”, with a NS direction, with a rejection that varies from a few meters in the Caramanico area, up to 4 km in the area between Passo S. Leonardo and Campo di Giove. To the west of this fault, we can observe the depression of the Orta Valley, which is filled with terrigenous sediments, and which develops for about 20 km with a NW-SE direction.

2.2. Data Analysis

For all the aspects necessary to develop the simulations discussed in this paper, flow rate, rainfall, infiltration and evaporation were the fundamental parameters that were mandatory to be measured and/or evaluated. Since both infiltration and evapotranspiration cause a loss of the inventory of surface water that flows, to simplify the simulations, the two phenomenologies were fictitiously merged into a single model. In the calibration process, the model parameters were evaluated mainly on the basis of the infiltration values taken from the literature [1,16] and then modified appropriately to also include the evaporation process. The rainfall data (in mm/day) were acquired from the annals of ‘Regione Abruzzo’ [22], related to 11 pluviometric stations located within the selected area. To estimate the water balance and therefore the total flow that fed the entire basin during the selected time periods, only the reliable flow rates measured at Pescara-Maraone and Tirino-Madonnina, located at the borders of the hydrographic basin, were considered. The station at S. Teresa was selected as the “control point” for the comparison between the daily measured experimental data and the corresponding flow rates, numerically predicted through simulations. A hydroelectric power plant is located in the selected territory and various intake and release systems are placed along the river. Unfortunately, not much information is available regarding its hydrological management. However to simplify the modeling and since the global water balance of the activity of the plant is conserved so that the only possible impact could be just a minor variation on the flow measured downstream, this type of perturbation was not considered in the simulations. A careful preliminary study was conducted on the data reported in the annals [22]. A first analysis aimed to verify that in correspondence with time ranges affected by a peak in the rainfall data, a corresponding peak was, reasonably, recorded in the flow-rate data relating to the S. Teresa station. Accordingly, the time series from 2000 to 2007 was selected as the most suitable and reliable for the study, but nevertheless was affected by the occurrence of some discrepancies. Figure 2 shows the river flow rate (m³ s⁻¹) and the total rainfall (mm/d) measured at the experimental stations, ranging from 2000 to 2007. Further, inspection of the figures allow us to highlight when, in correspondence with a rainy event, there is the relative flood wave (or peak of flood).

The letters A, B, C, D, E in the graphs indicate exceptional river flow-rate peaks, which are reported out of the scale to allow a better visualization of the experimental time-series.

The flow rate measured at the S. Teresa station (red line Figure 2), located near the Pescara River mouth, collects the two contributions that are currently considered most significant and reliable for the water balance of the selected hydrographic basin. The first is provided by the inputs (black line Figure 2) measured at Pescara-Maraone and Pescara-Tirino which are the stations located at the border of the basin (Figure 1), while the other one is the contribution of all the precipitation inside of the basin during the observed transient. However, from a simple inspection of Figure 2, a systematic discrepancy emerges between the selected inflow data entering the basin (black line) and the out-flow measured at Santa Teresa (red line). It should be noted that the discrepancy is observed in all the time series of the analyzed annals [22] from 2000 to 2007 and persists even in the periods during which no rainy events are recorded. The causes can be multiple. It may be that the identified hydrographic basin may not coincide with a larger hydrogeological basin that can contain it; there is a differentiated aquifer–river exchange: upstream the river feeds the aquifer and downstream the aquifer feeds the river, and this could be a starting point for future research to be investigated [23].

At this point, we had to identify a calibration period whose values were used for validation in different years, preferably prior and subsequent to the calibration period. Accordingly, the January 2003-year time series was selected for calibration, while the first and the last time series resulting from the reliability analysis of the annals [22], respectively years 2000 and 2007, were selected for validation.

The systematic lack of water inventory was compensated by increasing the flow of water entering the basin. Consequently, as a parametric attempt, by a simple inspection of Figure 2, the following increases in inlet flow to the Pescara River were added: 28 m³ s⁻¹, 18 m³ s⁻¹ and 14 m³ s⁻¹, respectively for the time series of January 2003, the two semesters of the 2000 and the first half of 2007. These increases were sufficient to compensate, at least in part, the discrepancies.

2.3. Selected Mathematical–Numerical Approach

Among many available numerical tools, RiverFlow2D, 2017-year version, was selected. This code comes with a general-purpose graphical user interface (GUI), as pre and post processor. The numerical approach is based on the finite volume method (FVM) aimed at solving the depth-averaged 2D shallow water equations [11,12,13,14,15,16], assuming that the vertical depth of the flow is negligible compared with the other flow’s dimensions. This is an acceptable compromise between the need to adopt a numerical tool as realistic as possible and the need to avoid long computation times such as those required by 3D approaches. To also consider bed level jumps (considered as source terms in the related differential equations), due to the geomorphology of real complex landform, the code includes the ‘augmented approximate Riemann solvers’ approach [13]. In addition, the selected code allows numerical parallel simulations based on the graphic processing unit (GPU) card, which can run simulations many times faster than in single processor computers. In the version of the code selected to develop this work, important parameters such as the Manning coefficient, and other important physical-geo-mechanical parameters are considered constant throughout the transient. In the following, a summarized description of the utilized mathematical and numerical approaches is briefly discussed, freely adapted from the RiverFlow2D’ manual [16] in which many more details are reported.

2.3.1. Hydrodynamic Unsteady Flow Models

The resulting system of the coupled partial differential equations, based on the water volume conservation, water momentum conservation and shallow water assumption, are successively shown:

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

(1)

$$U={(h,q_x,q_y)}^T;F={\left(q_x,\frac{q_y^2}{h}+\frac{1}{2}g{h}_{}^2,\frac{q_x{q}_y}{h}\right)}^T\hspace{1em}G={\left(q_y,\frac{q_x{q}_y}{h},\frac{q_y^2}{h}+\frac{1}{2}g{h}_{}^2\right)}^T$$

(2)

$$S={\left(0,\frac{p_{bx}}{\rho }-\frac{{\tau }_{bx}}{\rho },\frac{p_{by}}{\rho }-\frac{{\tau }_{by}}{\rho }\right)}^T$$

(3)

where: $h$ is the water depth, $q_x=hu$ and $q_y=h\upsilon$ the unit discharges resulting from the shallow water approach; $(u,\upsilon )$ are the depth averaged components of the velocity vector $u$ along the $x$ and $y$ coordinates respectively; $g$ is the gravity acceleration; $g{h}^2/2$ the flux obtained after assuming a hydrostatic pressure distribution in every water column, as common practice in shallow water models. The source term $S(U,x,y)$ incorporates the effect of pressure forces $p_{bx}$ and $p_{by}$ over the bed and the tangential forces ${\tau }_{bx}$ and ${\tau }_{by}$. Moreover, $p_{bx}/\rho =gh{S}_{0x}$ and $p_{by}/\rho =gh{S}_{0y}$, while $S_{0x}=-\partial z/\partial x$ and $S_{0y}=-\partial z/\partial y$ are the bed slopes of the bottom level $z$. The tangential forces are characterized by selected rheological laws discussed briefly in the following sub-sections.

2.3.2. Manning’s Coefficient

The Manning’s “M” coefficient usually accounts for the effects of bed roughness, internal friction and variations in shape and size of the channel cross section, obstructions, or river meandering. For these reasons, the M coefficient was partially used as a tune parameter in the calibration process in order to adjust the numerical results to measured data. The initial guessed values, employed in this paper, were based on the tables reported in related studies [20,27]. Simulations aimed at validating the selected approaches were performed considering the value of 0.02, considering medium vegetation density, steep banks, a few trees, and brush. Many calibration tests were conducted considering other values, as shown below.

2.3.3. Infiltration and Evapotranspiration

Infiltration represents another component of the hydrological balance, and it can be defined as the process by which surface water enters the soil and, therefore, constitutes a loss for the surface balance, even if it could then be returned by internal sources or outside the selected catchment area. On the other hand, the phenomenon of evapotranspiration also involves a net loss of the water balance of the aquifer. However, in the calibration phase, to simplify the parametric analyses that are very demanding of computer time, the coefficients of the selected Horton model (Equation (4)) were considered as tuning parameters, including all the phenomena of surface water loss and therefore also the evapotranspiration. This expedient allows us to reduce the number of parameters to be calibrated by appropriately adjusting the parameters of the selected semiempirical infiltration model.

The Horton model suggests an exponential equation for modelling the soil infiltration capacity. This commonly used type of model is particularly suitable for the parametric analyses performed in this work. In the selected Horton’s infiltration model, [26], the soil infiltration capacity $f_p$ (m/s) is based on the following exponential law:

$$f_p=f_c+(f_0-f_c)e^{-kt}$$

(4)

where $f_0$ and $f_c$ are, respectively, the initial and final infiltration capacities, both measured in m/s and k (1/s) is the rate of decrease in the capacity. Actually, the development of the Equation (4) is based on an experimental approach, even if a reasonable justification can be acquired by considering the temporal saturation rate of the soil compared to the intensity of the rainfall. Accordingly, the resulting runoff inventory would be determined by the balance between rainfall and infiltrated water. At the same given rainfall intensity, the minimum value of the water that could runoff is obtained when the initial conditions of the soil are dry, while the maximum value is reached when the soil is in saturation conditions. However, the related coefficients must be determined from experimental data, found in the RiverFlow2D’ manual [16], which were considered as the first attempt values as they include the evapotranspiration, but just only by a numerical point of view. Considering its specific characteristics, the Pescara River basin was divided in two areas, according to soil and slope criteria, with an upstream area characterized by a lower infiltration rate value, and a downstream area with a higher one. Other values were explored during the calibration tests, as shown below.

2.3.4. Numerical Solver

The numerical solvers of the “conservation partial differential equations (PDE)”, resulting from the application of shallow water models [14], can be identified, in the context of ‘initial and boundary value’ problems, with ‘source terms’. However, the search for solutions must face the challenge posed by the presence and/or the occurrence of discontinuities, due to both the wet–dry state transitions of the control volume and the discontinuities of the topography. Therefore, the numerical context is that of the “Riemann Problems (RP)”, which naturally appears in finite volume methods, Godunov (FVM) type approaches. For its solution, the ‘Riemann solver’ method [13] was implemented into RiverFlow2D. Within this framework, quantities such as rarefaction and shock waves, generated on discontinuities, appear as ‘characteristics’, identified by the ‘eigenvalues’ approach, along which the PDEs are transformed into ordinary differential equations (ODE), which are easier solve than PDEs. Roe’s approach is also part of the RP solution [13], whereby the integral of the approximate solution of the linearized RP, on a suitable control volume, is equal to the integral of the exact solution on the same control volume. To pursue this requirement in Roe’s method, the approximate solution can be defined using an approximate Jacobian matrix of the non-linear normal flow and two approximate matrices, constructed using the Jacobian eigenvectors. To obtain a completely conservative method, RiverFlow2D considers the complete system, including the hydrodynamic and transport equations. Mathematically, the complete system retains the property of hyperbolicity, implying the existence of a 4 × 4 Jacobian matrix for the 2D model. Accordingly, the general underlaying numerical approach is based on the finite volume scheme, by the integration in a volume or grid-cell Ω using Gauss theorem:

$$\frac{\partial }{\partial t}{\displaystyle {\int }_{\mathsf{\Omega }}Ud\mathsf{\Omega }+{\displaystyle {\oint }_{\partial \mathsf{\Omega }}E\cdot ndl=}}{\displaystyle {\int }_{\mathsf{\Omega }}Sd\mathsf{\Omega }}$$

(5)

where $U,F,G$ and $S$ are the vectors already described in (1), ∂Ω is the volume domain boundary, $n(n_x,n_y)$ is the outward unit normal vector to the volume Ω. Then, a piecewise representation of the conserved variables and an upwind and unified formulation of fluxes and source terms is applied [14]:

$$\frac{\partial }{\partial t}{\displaystyle {\int }_{{\mathsf{\Omega }}_i}Ud\mathsf{\Omega }+{\displaystyle \sum _{k=1}^{NE}{(E\cdot n-\overline{S})}_k{A}_k=0}}$$

(6)

where Ω_i is the volume of the i-th computational cell, whose area of the k-th edge face is $A_k$. Then, the approximate solution is always constructed as a sum of jumps or shocks, also involving rarefactions [14].

2.3.5. Optimal Time-Step Computation

The problem of identifying the optimal time-step is part of the search for the requirements to be met to obtain the stability of the selected algorithms. Riverflow2D uses an automatic implementing procedure, based on the solution of the RP. In 1D simulations the time step is taken small enough to avoid interaction of “numerical waves” emerging from the application of this particular type of approach; then the “equivalent distance” Δx/2 between neighboring meshes is introduced. In the 2D framework, considering unstructured meshes, the equivalent distance, calculated for each i-th cell, must consider the “volume” of the cell and the length of the shared edges:

$$L_i=\frac{A_i}{\underset{k=1,NE}{\max }{l}_k}$$

(7)

where $A_i$ is the ‘volume’ (area in 2D) of the i-th cell, $l_k$ is the length of the k-th edge of the i-th cell. Considering that each k RP (along the direction perpendicular to the k-th edge) is used to deliver information to a pair of neighboring cells of different size, the associated distance between i-th and j-th cell, $\min (l_i,l_j)$ should affect the time-step selection. Accordingly, at this point, the procedure based on RP, in the case of triangular elements, requires that the time-step size is limited by:

$$\mathsf{\Delta }t\le CFL\cdot \frac{\min (l_i,l_k) }{\underset{m=1,2,3}{\max }\left|{\tilde{\lambda }}_m\right|}$$

(8)

where CFL is Courant–Friedrichs–Lewy number (for 1D case CFL = (u∆t)/∆x), u is the actual flow velocity, ∆t the numerical time-step and ∆x is the mesh-size) to be defined by the user (for all simulations discussed in this paper the selected value was: CFL = 1), ${\tilde{\lambda }}_m$ is a velocity of the flow perpendicular to the m-th edge, emerging from the application of the 1D RP solver along the three directions in the case of triangular elements. Accordingly, the selected software automatically computes the time step that results to be variable during the transitory, proportional to the local grid-cell width, but also inversely proportional to velocity and depth.

2.4. Conceptual Model and Selected Inputs

The topographic data set was based on a 10 × 10 m grid resolution DTM obtained by merging all the smaller DTM files regarding the area, available on the Geo-site of Regione Abruzzo (Italy), UTM-WGS84 33N [24]. The boundary of the geometric domain was traced considering the extension of the water catchment area of the Pescara River, embedding all the measurement station (pluviometric, thermometric and hydrometric in each of the 11 points identified in Figure 1) placed within the area. Discharge and Free-Outflow conditions were imposed, respectively, at upstream and downstream border of the water catchment area, as commonly needed in computational codes that solve problems with partial differential equations with boundary conditions. No lateral boundary flow was assumed. The Discharge condition placed upstream accounts for the Pescara River inflow point into the study area (Figure 1). No water sources were assumed, in particular along the lateral borders where the Closed Boundary condition was set. This assumption is possibly the most important cause of the discrepancy between the incoming water flow and the measured water flow at the S. Teresa station. The beginning of the numeric transient, for the simulations related to years 2000, 2003 and 2007, was set respectively at 1 January 2000, 00:00 a.m., and 1 January 2007, 00:00 a.m. Moreover, the terrain under study was assumed in initial dry condition.

2.4.1. Numerical Meshing

The numerical tessellation of the selected DTM was performed through an unstructured mesh approach, based on linear triangular elements, applying the meshing tool that is supplied with the Riverflow2D code. Given the resolution of the DEM file (10 × 10 m), the maximum level of detail that could be imposed was of 10 m for the mesh as well. A first attempt was to assign this value to the entire study area, obtaining a number of cells of about 7 million. This number of cells, as expected, made the computational charge too onerous, expanding the calculation times excessively. A too-high level of detail would require a very high number of cells, which would result in excessive computing time; however, a too-low level would not produce appreciable results. By a consequence, a differentiation of the level of the numerical meshing approximation was carried out, so that it was kept detailed close to the Pescara River course and the surrounding area, setting the resolution of the mesh at 15 m, with a minimum of 10 m in the location of the Pescara measurement station at S. Teresa. In the rest of the area, a first attempt was to assign a value of 200 m in average, obtaining a number of cells of about 700,000. Moreover, it was necessary to create transition bands between the two areas, as a too-abrupt variation would have led to errors during the mesh generation process. After various tests, it was decided to place two transition bands, which are set along both the sides of the 15-m area, the first with a resolution set at 30 m, and right after the second set at 60 m. In this way, the computation time was faster, without affecting the river flow values at Pescara-S. Teresa (10 m average wide meshes); however, some points of water accumulation do appear along the drainage courses because such a mesh resolution value does not allow for a detailed enough rendering of the altimetry of the area. So, in some points, where the water was supposed to drain, nonrealistic blocks emerged as if there were a wall. This solution was not acceptable. Therefore, it was necessary to manually modify the automatic discretization of the spatial domain, processed by the RiverFlow2D code, by means of an appropriate position shift of the mesh nodes (selected by mouse) in order to obtain, in some critical area by a morphological point of view, more regular tessellation triangles. Furthermore, in addition to S. Teresa’s area, it was decided to create one more area set with the maximum level of detail allowed by the files resolution (10 × 10 m), at the canyon along the Orfento Valley (Figure 1), to improve the accuracy of numerical simulations.

Several numerical tests were carried out, not reported in order to avoid any further overload of the paper, aimed at obtaining the right compromise between the average aspect ratio of the elements of geometric tessellation (that defines quantitatively the distortion of a triangular element with respect to an equilateral triangle, optimal by a numerical point of view [12,13,14,15]), the level of detail and the computational cost. The optimal value of the total amount of cells was identified at 1.5 million (1,448,264 cells). Moreover, during the calibration of the automatic tessellation of the spatial domain, virtual, non-real points of water accumulation emerged which were eliminated by the final refinement

Figure 3 displays the numerical discretization of the spatial domain, employed for all the performed simulations, in agreement with the purpose of the paper, to meet both the requirement for numerical accuracy (high number of tessellation elements) and the request of an acceptable running time.

From Figure 3, it is possible to observe how the mesh-cell size decreases towards the river. The same figure, by way of example, also displays the plots of two final simulations of surface water flows within the basin, due to two rainy events relating to the year 2000, as, moreover, will be discussed in greater detail in the following.

2.4.2. Rainfall

Rainfall data was obtained by consulting the annals of the ‘Regione Abruzzo’ [22,23], which report the measurements of the 11 pluviometers within the selected area (Figure 1 and Figure 2). The point values measured in each of the 11 rainfall stations were assumed as the precipitation values within each of the 11 areas identified by the well-known technique of Voronoy tessellation.

2.4.3. Discharge Inputs

For the discharge of the Pescara River, the data measured by the hydrometric stations of Pescara-Maraone and Tirino-Madonnina were evaluated (Figure 1). Another station used as “control point” was Pescara S. Teresa. In fact, the experimental data of the Pescara River flow reported by the latter, are used as a reference in the comparison with the numerical results obtained from the simulation. These measurement stations (Figure 1) record the daily data of the flow rate of the watercourse in cubic meters per second (m³/s).

2.5. The Calibration Process

In preliminary tests it was observed that, with the same time series, there are two fundamental parameters that significantly influenced the results of the simulations based on the selected approach: infiltration’s parameters and the Manning’s number M.

The infiltration-evapotranspiration significantly affects the height of the river flow peaks, while the Manning coefficient affects the arrival times of the peaks and therefore the synchrony in the comparison between the results obtained numerically and the corresponding measured data. Accordingly, for calibration, parametric tests were performed. In the first step the Manning coefficient was kept constant to obtain the best values for infiltration-evapotranspiration, while in the successive step it was varied.

2.5.1. Infiltration-Evapotranspiration Calibration; Merged Model (InF_Ev)

In this first step the Manning’s coefficient was not yet calibrated and was chosen equal to 0.1, as an attempt based on the generic RiverFlow2D’s guideline. As expected, a discrepancy was observed between the arrival times of the numerically calculated peaks and those measured experimentally, equal to 48 h.

In order to also graphically facilitate the calibration of the numerical model with respect to the experimental peaks of the water flow (variation of the InF_Ev parameters only) a fictitious time translation of the original numerical curves was carried out. Figure 4 shows the dashed green curve obtained without translation and its translated solid blue curve, synchronized with the experimental data (red curve).

A first attempt (S1) was made assigning the following values: decay coefficient $k=0.00115$, final InF_Ev rate $f_c=1.94\cdot {10}^{-6}$; initial InF_Ev rate $f_0=1.41\cdot {10}^{-5}$. The results show how the set values caused a too-high InF_Ev process, preventing the drainage phenomena to be correctly reproduced and obtaining low amplitude river flow peaks (Figure 4a, S1). Under these considerations, another simulation, S2, was performed with the following, some decreased, InF_Ev rate values: $k=0.00115$, $f_c=2.78\cdot {10}^{-7}$ $f_0=2.78\cdot {10}^{-8}$. This time the river numerical flow peak amplitude was higher and closer to the measured values, but still not acceptable (Figure 4b, S2). Accordingly, we refined the model, considering that different areas could reasonably have different infiltration coefficients. Given this consideration, the Pescara River basin was divided into two areas, each defined according to a slope criteria. Reasonably, the steeper a territory, the higher the occurrence of more runoff than the infiltration process. Therefore, an upstream area, including both the territory shared with part of the Apennine chain and its piedmont area, characterized by a lower InF_Ev rate value and a downstream area including the coastal plain with a higher values, were defined. Hence, after some attempts the following numerical values were finally selected: for the downstream area, decay coefficient $k=0.00115$, final InF_Ev rate $f_c=2.78\cdot {10}^{-8}$, initial InF_Ev rate $f_0=2.78\cdot {10}^{-9}$, while for the upstream area, $k=0.00115$, $f_c=1.816\cdot {10}^{-8}$, $f_0=1.816\cdot {10}^{-9}$. The results of the last attempt allowed us to obtain a better fitting between the numerical values and the measured values.

Therefore, the values of the InF_Ev parameters identified through the last simulations were selected to carry out the further calibration and validation (Figure 4c, S3).

2.5.2. Manning’s Number Parametric Calibration

After having acquired the Inf_Ev calibration, the next fundamental step was the identification of an adequate value of the Manning number for the roughness of the soil. Initial attempts showed that this parameter significantly influences the temporal correspondence of the numerical flux peaks with the measured data, as explained below.

The first attempt was to assign $\mathrm{M}1=0.1$ on the entire area. However, the time of arrival at the S.Teresa station of the peak flow, predicted numerically, was 48 h late compared with the actual measured value (Figure 5). Another attempt was explored dividing the territory under study into two zones: to the main drainage area surrounding the riverbed path was assigned $\mathrm{M}2=0.035$, while to the rest of the basin $\mathrm{M}2=0.03$. The delay of the flow peak arrival time decreased significantly, equal to 16 h. It was a better result but not acceptable yet. Moreover, the separation of the basin into two areas, for Manning’s number calibration, did not appear as necessary. Therefore, in the next attempt, a single homogeneous value $\mathrm{M}3=0.02$ was decided to assign to the whole area, making the model more streamlined and therefore decreasing the computational charge. Figure 5 shows that a good match was obtained between numerical simulations and actual data. The parameters acquired through the calibration tests were then used for all subsequent validation simulations.

2.5.3. Useful Features of the Selected Numerical Tool

In order to give some operational features, useful for the efficient use of the selected code, some further details are provided. RiverFlow2D provides some tools that have proved fundamental to perform the simulations developed and discussed in this paper.

The HOTSTART file may be used to restart a simulation from previously computed results. By default, the file contains the name of the last report. This command was particularly useful when it was necessary to restart the simulation from a very specific moment (after a system crash, for example).

Once the computation phase was over, the software generated the HDF5 output file containing a code producing an animated file and therefore a different and more interactive form of visualization of the results. Accordingly, it was possible to set visualization filters to make it easier to focus on the aspects of the simulation that were considered most relevant. In particular, the “Dry-Wet” filter was activated since we were interested in a convenient display of the behavior of surface waters over the entire study area (Figure 3). This visualization mode allows to observe the water drainage distribution.

2.5.4. Parallel Platform through GPU

Riverflow2D allows to take advantage of the parallel calculation by making use of the numerous processors present in the GPU, to speed up the computational time. Of course, the processing speed depends on the type of hardware available and therefore also on the graphics card used. Accordingly, some details of the computer used to perform the simulations discussed in this paper are given:

• Central Processing Unit (CPU) Intel i7-7700 3.6 GHz;
• #Cores: 4;
• #Threads: 8;
• GPU Nvidia GeForce GTX 1060 6 GB;
• #GPU Cores: 1280;
• Clock: 1506 MHz;
• Ram: 32 Gb.

To compare the differences in computation times, a test was conducted by launching two simulations of the developed model, whose details are reported in the following Figure 6: simulation A (CPU only): actual transitory time: 216 h from 1 January to 9 January 2003, computation time 8226 min (~137 h); simulation B (CPU + GPU): actual transitory time 216 h from 1 January to 9 January 2003, computation time: 272 min (~4.5 h).

Parallel computation (CPU + GPU), applied to the analyzed phenomena, required only $2.1\%$ of the actual transient time ($4.5/216=1/48$, accordingly 48 times faster), while, for the same simulation, only CPU approach required over $63.4\%$ ($137/216\cong 0.63425$). Accordingly, the CPU + GPU technique, was about thirty times faster than the CPU technology ($137/4.5\cong 30.4$). Under these considerations, it can be stated that, regarding the utilized hardware, the computation time for simulating a period of 1 full year can be estimated to be equal to just 7.7 days including the parallel calculation (CPU + GPU) instead of 231.5 days (~7.7 months) using only the CPU.

3. Results and Discussion

After the calibration step, some simulations were performed to validate the feasibility of the discussed approach. The parameters identified through calibration may vary depending on the year chosen and the different seasons of the year. However, in this paper we considered all parameters constant over time. The effectiveness of the procedure will be demonstrated in the validation step. Of course, further analysis could refine the methodology.

The validation of the process certainly had to cover a wide time interval, therefore the two semesters of the year 2000 and the first semester of the year 2007 were selected, i.e., respectively before and after the calibration year 2003. As previously observed, from a first remark of Figure 2, a systematic discrepancy emerges between the inflow data at the entrance to the basin (flow rate measured at Pescara-Maraone summed to the measurements at Tirino-Madonnina; black line) and those detected at Santa Teresa. The origin of these kind of anomalies are not yet completely known. However, to balance the water flow correctly, the inlet flow had to be synthetically increased. As a result, 18 m³/s were added to the base-flow of the 2000 time series, while 14 m³/s was added to the base-flow of the 2007 time series. Figure 7 shows the good match between the measured values and the numerical values predicted after calibrating the model and adjusting adequately the input, as described above.

Regarding the year 2000, the simulations for the period of November to mid-December overestimate the measured values, while the peak of mid-December is well simulated, even if the predicted value persists to be slightly higher. The blue arrows show some peaks resulting from the numerical model, which, however, do not match adequately the measured actual values, although in correspondence some not negligible rainfall events occurred. The reasons, in these cases, could be attributed to the measuring instruments, management of the hydropower plant and/or to an inadequate input distribution of rainfall, whose measurements were localized and not diffused.

To quantify the deviations between the numerically predicted values and the measured values, some statistical analyses were carried out. Figure 8 shows the linear correlations for the simulations regarding the calibration (Figure 8a) and the validation of the first half of 2007(Figure 8a,d), while both linear and cubic regression the validation processes relating to year 2000 (Figure 8b,c). By the calibration step a rather high value of ${\mathrm{R}}^2$, equal to $0.9378$, was gained (Figure 8a).

The validation step, performed by adopting the physical parameters already calibrated, indicated that the simulations of the first half of the two selected years showed a better correspondence with the measured values (${\mathrm{R}}^2=0.8386$ year 2000, Figure 8b; ${\mathrm{R}}^2=0.9112$ year 2007, Figure 8d) compared to the only second half considered (${\mathrm{R}}^2=0.6765$ year 2000, Figure 8c). On the other hand, this circumstance may not be surprising since the calibration step considered only the month of January, therefore belonging to the first half of the year. Evidently, for the second half of the year, the need for a different calibration could emerge. By means of cubic correlation, slightly better values were obtained (${\mathrm{R}}^2=0.8491$ first semester year 2000, Figure 8b; ${\mathrm{R}}^2=0.6979$ second semester year 2007, Figure 8c). In this case, ${\mathrm{R}}^2$ indicates the degree of reliability (from 0 to 1) of the hypotheses, the assumptions, the applied modeling and the numerical description of the territory for obtaining a correct correspondence between the values of the flow water predicted numerically and measured values.

Moreover, the validation carried out shows that the predictions of the flow rate were particularly satisfactory for the first two semesters of the year 2000 (${\mathrm{R}}^2=0.8386$) and 2007 (${\mathrm{R}}^2=0.9112$), while the matching was quite satisfactory for the second half of the year 2000 (${\mathrm{R}}^2=0.6765$). This observation will also provide further indications to direct future investigations.

The regression provides global indications. However, it does not indicate nor describe which are the distributions of the deviations. In this regard, the following index has been introduced:

$$\mathrm{RDMN}(\%)=\frac{{\dot{Q}}_{meas}-{\dot{Q}}_{pred}}{{\dot{Q}}_{meas}}$$

(9)

where $\mathrm{RDMN}(\%)$ is the acronymous of ‘Relative Difference between Measured and Numerical values’, ${\dot{Q}}_{meas}$ measured flow rate, ${\dot{Q}}_{pred}$ numerically predicted flow rate.

Accordingly, Figure 9a shows the deviations related to calibration, while Figure 9b,c indicate the deviations resulting from the validation step.

A first inspection reveals that the most important deviations are observed at higher rainfall peaks for both calibration and validation, highlighted in the figures by red and blue arrows. Furthermore, it appears that the predicted values overestimated the measured values (negative RDMN), specifically when more intense rain peaks are observed, apart from some exceptions (three events, during 2000, blue arrows in Figure 9b). This characteristic of the applied model certainly satisfies the safety requirements, but could be too conservative, leading to dangerous announcement of false flood alarms.

As already highlighted above, the reasons of these circumstances should be explored both by a mathematical-phenomenological level (important phenomena could have been excluded or the numerical discretization of the territory may not have been sufficiently adequate) and by rain gauges network and measuring instruments levels.

To better visualize the details of the distributions of the deviations, Figure 10 shows the box plots that highlight the quartiles, the medians, the average values, the minimum and maximum values. The distribution spread is highlighted through the range between the first quartile to the third quartile ($(-25\%÷+75\%)$, maximum for calibration $(-5.87\%÷+18.85\%)$, minimum $(-3.0\%÷+0.9\%)$) for the validation of the first half of 2007. The resulting average was positive for calibration, while all negative for validations ($+8\%$ calibration, January 2003; $-1.28\%$ validation, first semester 2000; $-5.52\%$ validation, second semester 2000; $-1.19\%$ validation first semester 2007). It is also noted that the smallest average deviation surveyed was obtained for the validation of the first half of year 2007.

In order to get a better visualization of the RDMN distribution, by means of histograms, the Freedman–Diaconis rule [33], which not only considers the sample size, but also considers the spread of the sample, was selected:

$$binwidth=\frac{2(q_3-q_1)}{\sqrt[3]{n}};bins=ceil\left(\frac{max(x)-min(x)}{binwidth}\right)$$

(10)

where (referred to Figure 10): $q_3$ stands for third quartile $(75\%)$, $q_1$ stands for the first quartile $(25\%)$, $n$ stands for sample size and $\left(max(x)-min(x)\right)$, respectively the maximum and minimum values of the elements belonging to the sample. Furthermore, the representation by means of histogram could identify analytic functions that can possibly describe the statistical distribution under consideration. Accordingly, from Figure 10 and Equation (10), the

Table 1. RDMN statistical parameters related to Freedman–Diaconis rule.

	$\mathit{n}$	${\mathit{q}}_1$	${\mathit{q}}_3$	$\mathit{m}\mathit{i}\mathit{n}(\mathit{x})$	$\mathit{m}\mathit{a}\mathit{x}(\mathit{x})$	Bins
January 2003	31	−5.87	18.85	−31.12	35.04	4
1 semester 2000	181	−3.77	2,41	−23.69	22.21	21
2 semester 2000	183	−9.13	−1.15	−53.69	28.65	29
1 semester 2007	181	−3	−1	−15	10	35

follows.

The following are the selected functions:

$$\begin{array}{l}Gaussian4\mathrm{parameters}=a\cdot e^{\left[-0.5\cdot {\left|\frac{(x-\mu )}{\sigma }\right|}^c\right]};\\ Gaussian{5\mathrm{parameters}=\mathrm{y}}_0+a\cdot e^{\left[-0.5\cdot {\left|\frac{(x-\mu )}{\sigma }\right|}^c\right]}\end{array}$$

(11)

where ${\mathrm{y}}_0$, $a$, $\mu$, $\sigma$, $c$ are parameters numerically identified by one-way ANOVA test [34].

It is worth underlining that the values of the statistical parameters obtained with the technique shown in Figure 11 can be compared with the results obtained by constructing the box plots shown in Figure 10.

4. Conclusions

Despite the considerable complexity of the studied phenomena, the identified approach based on the ‘shallow water modelling’, an advanced mathematical-numerical model able to reproduce reality in almost the best possible way, proved to be feasible to simulate well the behavior of the surface water flow and the hydrological balance of a catchment area similar to the Pescara basin, selected as a test case. The discussed modeling can also be a useful integrative tool in designing further hydrogeological experimental campaigns, over the different seasons, with different weather conditions, aimed to acquire a more robust water basin inventory. It is worth underlining that one of the important objectives of the work was to verify the possibility of using a valid and advanced modeling tool to explore the consistency of all available data and to identify which of them should be subjected to further investigations, also by an experimental point of view. The results, shown in the graphs and in the discussion, were remarkably encouraging.

On the other hand, from the inspection of Figure 7 and in particular from Figure 9, the emergence of statistical outliers in correspondence to intense rainfall rates has suggested the need for further investigations in order to improve not only the mathematical-numerical model selected, but also the analysis and methods of acquiring experimental data. This point is particularly important as the current climate change and its developments evolve, and lead to further occurrences of extreme events such as very violent rainfalls, which are limited in time and would need specific models.

Therefore, in short, based on the results and in the face of different needs (suggestions for further measurement campaigns, management of the territory, etc.) we could define the quality of the applied modeling between good and very good.

Of course, the simulation of a complex phenomenology affecting a large area, within a period of one year or more, combined with the necessity of an accurate level of numerical discretization of the topography, appears to be very demanding in terms of calculation time, when only CPU platforms are employed. The use of parallel computing, even if on desktop or laptop hardware, revealed mandatory to perform parametric simulations on such large spatial scales and over long periods of time in, not only satisfactory, but acceptable computation times. By taking advantage of the GPU processors, computation times were made about 30 times faster by the hardware available to us. A simulation re-quiring 30 days running only on CPU hardware would be performed in only 1 day by the inclusion of the parallel calculation (CPU + GPU). Of course, a more powerful hardware, like Tesla Nvidia GPU Card A100, in accordance with what is discussed in this study, would be more suitable allowing us to reduce the computation time even more.

A detailed description of the hydrology, rainfall, water sources, flow rate input, were also fundamental information. In our paper, the data were acquired from the documentation made available by the Regione Abruzzo and subsequently analyzed, processed, and introduced into the model.

The numerical outcomes were strongly affected by the correct calibration of the spatial discretization, the Manning’s roughness coefficient and the infiltration-evapotranspiration rate. The simplified inclusion of the evapotranspiration effect as an adjustment of the actual parameter values of the selected infiltration model to obtain the match between the numerical simulations and the experimental results, proved to be effective.

Many calibration tests were performed covering 744 h of January 2003, compared successfully with the experimental measurements obtained from the available annals. So, a further outcome of this work was achieved by understanding and identifying how the different parameters of the selected model affected and could affect quantitatively the simulated phenomena.

Then, the validation simulations involved three semesters, the entire year 2000 and the first half of 2007. A more detailed morphological digital model would have been desirable, to perform more accurate numerical simulations of the phenomena.

Moreover, the concentration time of the basin was numerically estimated in about 22 h, consistent with the actual value. The application of the CPU-GPU allows us to simulate a transient of this duration in about 30 min (22/48 = 0.46 h, on our hardware). Accordingly, another important outcome of this study is certainly the possible utilization of the developed model as a part of an early warning system of floods.

A further consideration emerged from this study. The selected hydrographic basin may not coincide with a possible larger hydrogeological basin that could contain it. There is an aquifer–river exchange, in which the river feeds the aquifer up-stream, while downstream the aquifer feeds the river, and this is certainly a topic to be studied and deepened in future research activities that can clarify the lack of water inventory between the flow rate at the mouth of the Pescara river and the main tributaries identified.

Another important opportunity, not discussed in depth in this paper, would also be the water management of the area, for example irrigation and the impact of hydroelectric plants.

The developed model so far, based on GPU, on an advanced numerical approach, and on hardware such as a simple workstation and laptop, could be useful not only for the validation of the available data, such as for the water balance of a river basin, but also as a basis for applications in sediment transport, influence (or instabilities) of hydric structures such as dams or bridles, debris flow phenomena [7,8,17], transport of pollutants, floods [18,19,35] and many other aspects involved in ‘Mitigation Projects of Territory Risk’ and could therefore, by consequence, be useful in the context of the public territory management (‘Civil Protection’ in Italy).

This will surely be a starting point for future works, including also further analysis over other years, on other fluvial basins, aimed at providing more exhaustive validation of the model.