Qom—A New Hydrologic Prediction Model Enhanced with Multi-Objective Optimization

applied sciences Article Qom—A New Hydrologic Prediction Model Enhanced with Multi-Objective Optimization Gustavo R. Zavala 1,†, José García-Nieto 2,† and Antonio J. Nebro 2,*,† 1Subsecretaría de Obras Hídricas y Proyectos Especiales, W3400BCN Corrientes, Argentina; [email protected] 2Departamento de Lenguajes y Ciencias de la Computación, ITIS Software, University of Málaga, E.T.S. de Ingeniería Informática, Campus de Teatinos, 29071 Málaga, Spain; [email protected] *Correspondence: [email protected] † These authors contributed equally to this work. Received: 6 November 2019; Accepted: 24 December 2019; Published: 28 December 2019   Abstract: The efficient calibration of hydrologic models allows experts to evaluate past events in river basins, as well as to describe new scenarios and predict possible future floodings. A difficulty in this context is the need to adjust a large number of parameters in the model to reduce prediction errors. In this work, we address this issue with two complementary contributions. First, we propose a new lumped rainfall-runoff hydrologic model—called Qom—which is featured by a limited set of continuous decision variables associated with soil moisture and direct runoff. Qom allows to separate and quantify the volume of losses and excesses of the rainwater falling in a hydrographic basin, while a Clark’s model is used to determine output hydrograms. Second, we apply a multi-objective optimization approach to find accurate calibrations of the model in a systematic and automatic way. The idea is to formulate the process as a bi-objective optimization problem where the Nash-Sutcliffe Efficiency coefficient and percent bias have to be minimized, and to combine the results found by a set of metaheuristics used to solve it. For validation purposes, we apply our proposal in six hydrographic scenarios, comprising river basins located in Spain, USA, Brazil and Argentina. The proposed approach is shown to minimize prediction errors of simulated streamflows with regards to those observed in these real-world basins. Keywords: hydrologic model; prediction; multi-objective optimization; metaheuristics 1. Introduction The modeling over time of the volumes of rainfall collected in a watershed in terms of losses and excesses, quantifying the accumulation on the surface and the movement by gravity superficially to topographically lower areas is currently an indispensable task in hydrology, since it allows experts to efficiently project water management infrastructures and design ecological strategies. Examples of these applications include the generation of hydroelectric power facilities, water reservoirs for human consumption, agricultural and livestock exploitation, water quality management, and so forth. It can also be used as a base and support for other studies, such as the dissolution of pollutants in a river, the erosion and sedimentation in riverbeds and basins, the impact of urbanization on the increase of the waterproof surface, as well as the formulation of contingency plans for emergencies due to deficiencies (droughts) and excesses of water (floods). However, realistic hydrologic models often require many parameters to be properly tuned to reduce prediction errors with regards to observed scenarios. Finding the most appropriate settings of these parameters can be formulated as a complex optimization problem [ 1 ]. This was outlined in early studies [ 2 , 3 ] and more recently in Reference [ 4 ], where the importance of using bio-inspired Appl. Sci. 2020,10, 251; doi:10.3390/app10010251 www.mdpi.com/journal/applsci Appl. Sci. 2020,10, 251 2 of 28 multi-objective approaches was emphasized due to their ability to explore the search space generated for each problem river basin scenario and for several objectives at the same time. The main contribution of this work is the combination of a new hydrologic model—called Qom—with the application of a multi-objective approach to find efficient calibrations of the model when applied to real-world scenarios. The motivation driving us is to provide an efficient tool, able to provide accurate predictions when applied to heterogeneous river basin real scenarios, considering different time periods and climatological conditions. Compared to other models, Qom is featured by requiring a small number of parameter to be tuned—maximum storage water volume on surface, maximum storage volume in the soil and volumetric conductivity coefficient. Qom allows the efficient separation and quantification of the volume of losses and excesses of the rainwater falling in a hydrographic basin. Qom considers antecedent moisture and it is applicable to urban and rural watersheds, for permeable and impermeable soils. In addition, it determines the runoff in three layers (two superficial ones and a groundwater one) and it considers simple and continuous rainfall events for long periods of time of several years, including dry and wet periods with very short time records from several minutes to one day. A Clark’s [ 5 ] model is then used to obtain predicted output hydrograms that include the transit of excesses with the effects of delays and storage. The calibration of Qom is carried out by formulating it as a bi-objective optimization problem, which is solved with the combination of a number of bio-inspired algorithms belonging to metaheuristics, a family of non-exact optimization techniques [ 6 ]. We have taken seven algorithms that are representative of the state-of-the-art, namely NSGAII [ 7 ], OMOPSO [ 8 ], SMSEMOA [ 9 ], MOEAD [ 10 ], AbYSS [ 11 ], SMPSO [ 12 ] and MOCell [ 13 ]. These algorithms are provided by the jMetal framework for multi-objective optimization with metaheuristics [ 14 ], an open-source project, developed in the Java programming language, and released under an MIT (Massachusetts Institute of Technology) licence. Qom has also been implemented in Java to be integrated into jMetal as a black-box optimization problem, with vectors of decision variables as inputs and objective function values as outputs. The version of jMetal including Qom as well as all the utilities and scripts (most of them implemented in R) we have used in this work, are freely available to the community of interested users. For validation purposes, six realistic hydrographic scenarios comprising river basins, located in different regions of Spain, USA, Brazil and Argentina, have been modeled and evaluated with different size conditions, climates, topographies, heterogeneous soils and with a series of very short time periods of 10, 15, 30 and 60 min and 24 h. The results are mathematically validated by three gauges parameters and exponential functions that relate evapotranspiration, surface stored water, water stored in the soil with the infiltration and percolation processes. The conducted experimentation shows that the combination of Qom with the multi-objective based calibration procedure can provide experts with successful trade-off solutions, which minimize the prediction errors of simulated streamflows with regards to those observed in real-world basins. The remainder of this paper is organized as follows. Section 2provides the algorithmic background. Section 3presents a series of related works in the state-of-the-art. In Section 4, the proposed Qom hydrologic model is detailed and formulated in the form of multi-objective optimization problem. Experimental results and discussions are presented in Section 5. In Section 6, a series of discussions have been included. Finally, conclusions and future work are commented in Section 7. 2. Background on Multi-Objective Optimization with Metaheuristics This section is devoted to providing a basic background of multi-objective optimization with metaheuristics, which is needed to understand the calibration scheme we propose. Many real-world optimization problems can be formulated as the maximization and/or minimization of two or more conflicting functions or objectives at the same time. These problems are known as multi-objective optimization problems, and they can be found in many disciplines such as civil engineering [ 15 ], economics [ 16 ], telecommunications [ 17 ], bioinformatics [ 18 , 19 ], Appl. Sci. 2020,10, 251 3 of 28 agriculture [ 20 ], and so forth. A multi-objective optimization problem can be defined as follows (we assume minimization of all the functions, without loss of generality): Definition 1 (Multi-objective optimization problem) . A multi-objective optimization problem is a tuple (S , f , g , h) , where S6=∅ is called the solution space (or search space), f= [f1 , f2 , . . . , fk] , with k>= 2, is a vector of functions, where fi:S→R , are the objective functions and g= [g1 , g2 , . . . , gm] and h= [h1 , h2 , . . . , hp] are also vectors of functions, where gi:S→R and hi:S→R are the constraint functions. Thus, solving an optimization problem consists of finding a set of solutions X∗⊆S such that, for all x∗∈X∗: fj(x∗)≤fj(x),∀x∈S, (1) for some 1≤j≤k, subject to: gi(x∗)≤0i=1, 2, . . . , m(2) hi(x∗) = 0i=1, 2, . . . , p, (3) where gi,hj:S→R, i =1, ..., m, j =1, ..., p are the constraint functions of the problem. The set of all values satisfying all the equality and inequality constraints defines the feasible region Ω , so any point ~ x∈Ω is a feasible solution. For simplicity, in the following definitions we will consider Ω=S. A key difference between mono- and multi-objective optimization is that in the latter, vectors of values (instead of a single value) must be compared to determine which is the best of them. In this sense, a basic concept in multi-objective optimization is Pareto dominance, which is defined as follows: Definition 2 (Pareto dominance) . Given two vectors ~ x , ~ y∈Rk , we say that ~ x≤~ y if xi≤yi for i= 1, ..., k , and that ~ x dominates~ y (denoted by ~ x≺~ y) if ~ x≤~ y and ~ x6=~ y. Given two solutions in objective space, Pareto dominance indicates whether one dominates the other one or, on the contrary, that the compared solutions are non-dominated, that is, none of them is strictly better than the other in all the objective values. This way, solving the multi-objective problem can be viewed as the process of finding the set of solutions that dominates every other point in the solution space and all the solutions of that set is said to be Pareto optimal for that problem. Formally: Definition 3 (Pareto Optimality) . A solution x∗∈ S is Pareto optimal if it is non-dominated by any other solution x0∈ S. The set of solutions composed of all the Pareto optimal solutions is known as Pareto optimal set (or Pareto set): Definition 4 (Pareto Optimal Set).The Pareto Optimal Set P∗is defined as: P∗={x∈ S|xis Pareto optimal.} The solutions in the Pareto set belong to the variable space ( S ), and their correspondence in the objective space (Rk) is a set known as Pareto front: Definition 5 (Pareto Front).The Pareto Front PF∗is defined by: PF∗={f(x)∈Rk|x∈ P∗.} Appl. Sci. 2020,10, 251 4 of 28 The solutions in the Pareto front are usually referred to as non inferior,acceptable or efficient. The Pareto front is also known in some contexts as efficient frontier. From these definitions, we can conclude that the main goal of multi-objective optimization is to find the Pareto set of a given problem so that an expert in the problem domain can choose one or more solutions from the corresponding Pareto front according to some preferences. However, when dealing with real-world problems finding the Pareto front can be unpractical for a number of reasons [ 21 ] affecting the objective functions, such as non-linear linearity, NP-hard complexity, epistasis, and so forth. Furthermore, a Pareto front can contain a very large (or even an infinite) number of solutions. So, in practice, the goal usually is to find a high quality approximation of the true Pareto front containing a limited number of solutions (e.g., 50, 100, 500). The quality of these approximations is measured in terms of two properties—convergence, which implies that they must be as close as possible to the Pareto front, and diversity, which means they are uniformly spread. The concept of Pareto front approximation is illustrated in Figure 1, which includes the solutions found by an optimization algorithm when solving a multi-objective problem having two objective functions (named f1 and f2 in the figure). They are trade-off solutions in the sense that neither is better than the other in the two objectives. The true Pareto front in the example, depicted as a continuous line, has a convex shape and it can be observed that some of the obtained solutions do not completely overlap with the line (i.e., they have not fully converged to the Pareto front, so they are not optimal solutions) and the set of solutions covers all the Pareto front although it does not have a perfect spread out. 0 0.2 0.4 0.6 0.8 1 1.2 0 0.2 0.4 0.6 0.8 1 f2 f1 Pareto front approximation Pareto front Figure 1. Example of Pareto front approximation of a multi-objective optimization problem. The most popular techniques to deal with multi-objective problems are metaheuristics [ 22 ], a family of non-exact, stochastic optimization algorithms which do not guarantee to find the optimum of a problem, but do usually provide high quality solutions (near optimal and sometimes optimal) within a reasonable amount of time and resources. There are different ways of classifying metaheuristics, one of which distinguishes between bioinspired (or nature-inspired) and non-bioinspired algorithms [ 6 ]. The first category is by far the most well known and used, including techniques such as evolutionary algorithms, particle swarm optimization, or ant colony optimization. The reference algorithm in multi-objective optimization, NSGA-II [7], is an evolutionary algorithm. A metaheuristic typically follows an iterative process in which a set of tentative solutions P is manipulated somehow by a number of variation operators, aiming at progressively generating better solutions, as shown in the pseudo-code in Algorithm 1. In the case of evolutionary algorithms, Appl. Sci. 2020,10, 251 5 of 28 the solutions and the set P are named, respectively, individuals and population, and the variation operators are crossover and mutation. The new produced solutions are evaluated to compute their corresponding objective function values and they are used to update the population. The main loop of the algorithm is executed until a stopping condition is met. Algorithm 1 Template of a metaheuristic. 1: Set control parameters 2: P(0)←GenerateInitialSolutions() 3: t←0 4: Evaluate(P(0)) 5: while not StoppingCriterion( ) do 6: Q(t)←Variation(P(t)) 7: Evaluate(Q(t)) 8: P(t+1)←Update(P(t),S(t)) 9: t←t+1 10: end while 11: Return P(N) A reason for the popularity of metaheuristics is the development of software frameworks that provide algorithms that are representative of the state of the art plus benchmark problems and utilities for performance assessment [ 23 ]. As commented in the introduction, in this work we use jMetal, a Java-based framework for multi-objective optimization with metaheuristics. 3. Related Work Since the last decade, the optimal parameter calibration of hydrologic models has been tackled with mono- and multi-objective metaheuristic algorithms, usually providing successful results for different approaches of decision variables and objective functions. In particular, two early surveys [2,3] outlined the importance of using bio-inspired multi-objective approaches, as they are specially adapted to work as black-box optimizers with capabilities to explore the search space generated for each problem river basin scenario and for several objectives at the same time. In this regard, an interesting proposal was presented in Reference [ 24 ], where a sensitivity analysis and automatic calibration of a rainfall–runoff NedborAfstromnings model (NAM) [ 25 ] was carried out with nine decision variables to measure the quantitative and qualitative variation of results. In this work, an Elitist Non-dominated Sorting Differential Evolution (NSDE) [ 26 ] algorithm was used to optimize two objective functions consisting on Average Root Mean Squared-Error (RMSE) of peak and low flow Events. In Reference [ 27 ], a strategy that combined the principles of multi-objective optimization and depth-based sampling was presented with the title of Multi-Objective Robust Parameter Estimation (MOROPE). This algorithm applied multi-objective optimization to identify non-dominated robust model parameter vectors for the calibration of a distributed hydrologic model with a focus on flood events in a small pre-alpine and fast responding catchment in Switzerland. A thorough comparison was conducted in Reference [ 28 ] with a series of popular multi-objective algorithms (Borg, AMALGAM, GDE3, NSGA-II, OMOPSO, and SPEA2) for the calibration of the HBV (Hydrologiska Byråns Vattenbalansavdelning) conceptual rainfall-runoff model [ 29 ] with 14 real-valued decision variables and considering four objective functions. In this study, Borg and OMOPSO obtained, in overall, the best performance for the tackled instances. Also in 2013, the work presented in Reference [ 30 ] was centered on a river basin according to the distributed hydrologic model HEC-HSM (Hydrologic Engineering Center - Hydrologic Modeling System, https://www.hec.usace.army.mil/), for which a number of 10 decision variables and 4 objectives were approached with a multi-objective particle swarm optimization algorithm. Appl. Sci. 2020,10, 251 6 of 28 More recently, several variable balancing approaches for the exploration and exploitation of the NSGA-II, in the automatic parameter calibration of a HYdrological MODel (HYMOD) were evaluated in Reference [ 31 ]. These balancing approaches were compared with traditional static balancing methods (the two values are fixed during optimization) in a benchmark hydrological calibration problem for the Leaf River (1950 km2) near Collins, Mississippi. Another related study is found in Reference [ 32 ], where authors analyzed how adding daily Total Water Storage (dTWS) derived from the Gravity Recovery And Climate Experiment (GRACE) as extra observations, besides the traditionally used runoff data, improved the calibration of a conceptual hydrological model within the Danube River Basin. As calibration approach, four popular evolutionary optimization techniques (NSGA-II, MOPSO, PESA-II, SPEA2) were tested to calibrate the model. Results indicated that NSGA-II performed better than the other techniques for calibrating GR4J using GRACE dTWS and in situ runoff data. A common aspect in all these works is the use of already existing multi-objective optimizers (NSGA-II, SPEA2, MOPSO, etc.) with well-known hydrologic models usually requiring large number of parameters (HBV, HEC-HSM, GRACE, etc.), trying to find the best performing algorithm. Our proposal differs not only in the definition of a novel hydrologic model (Qom), but also in how the calibration with metaheuristics is conducted. Instead of looking for a particular technique, we take the practical approach of combining the results of seven algorithms, configured with commonly used settings, with the idea of combining the best solutions provided by all of them. 4. Hydrologic Calibration Strategy In this section, the complete calibration strategy of the hydrologic model is detailed. First, the Qom model is explained with the main parameters and formulations. Then, the model calibration is formulated as a multi-objective optimization problem, so its objective functions and constraints are defined. Finally, the optimization approach is described from the perspective of software and computational issues. 4.1. Qom Model The Qom model considers moisture, evapotranspiration, surface stored water and infiltration, and percolation processes, that is, those actions that cause water storage in the soil. It is applicable to permeable and impermeable soils and it determines the runoff in three layers, two of them superficial (permeable and impervious soils) and one another of underground flow. In order to determine the outlet flows in a given place in the basin in response to rainfall, we separate the total volumes of rainfall into volumes of losses and excesses discreetly over time, simulating for excess water the way in which runoff occurs, considering the transit and delay in the basin. Therefore, besides the rainfall data and morphometric characteristics of the basin, the model needs three initial tuning parameters, namely the maximum storage water volume on surface Rmax , the maximum storage volume in the soil Smax and the volumetric conductivity coefficient VC , which represents the volume of water that would move in three-dimensions in the soil towards the exit per unit of surface in a certain time, expressed in volume of water sheet per unit of time. These variables are emphasized in red color in the flow chart of Figure 2, which represents the mathematical model that corresponds to the equations commented in lines (1) to (12) of the this figure, involving the hydrologic natural processes of Qom. The computation cycle starts obtaining those required data from the problem specification given from the scenario instance, as well as the parameters (decision variables), to calculate the soil moisture states. A part of the rainfall defined with the variable P in the unity of the volume of the water sheet. Retentions in the basin are produced by surface depressions of the soil R , by the infiltration of water in the porosity of the soil S and percolates in permeable soils D . Another part of the volume is lost in the atmosphere through evaporation and plant transpiration EVT . The rest of the volumes would leave the basin superficially by a portion of impermeable areas Vi and permeable ones Vp . In addition, another part Vg of the volume can be Appl. Sci. 2020,10, 251 7 of 28 percolated into deep layers remaining in the basin as soil moisture or as a mass of water that will be added to the underground outflow of the river. Clark's method (runoff Vi, Vp, Vg) CT,Shape,Eq F low Translation (Synthetic Histogram) Flow Attenuation (Simple Linear Reservoir) PEVT[i] - = R[i] EVT [i] = R[i] R[i] = 0.0 R[i] - = PEVT[i] EVT [i] = PEVT [i] PEVT[i] = 0.0 R[i] += P[i] P[i] = 0.0 i=0->N, Δt, A, %P, P, PEVT, So, Ro, Qo, Smax, Rmax, VC Data and moisture initial PEVT[i]<= R[i] Potential evapotranspiration cover with excess water on the surface WITHOUT RUNOFF P[i] - = PEVT[i] EVT [i] += PEVT [i] PEVT[i] = 0.0 PEVT[i] - = P[i] EVT [i] += P[i] P[i] = 0.0 PEVT[i] < P[i] Total evapotranspiration of surface water Surface water accumulation with runoff Surface water accumulation without runoff P[i] = P[i]-(Rmax - R[i]) R[i] = Rmax P[i] >= (Rmax - R[i]) True False True aux = - R[i]/(Rmax-R[i]) z = 1-2.0 aux Sp = z * Smax Saux = max(S[i], Sp) ΔS = Saux -S[i] R[i] - ΔS > 0.0 R[i] - = ΔS S [i] = S[i] + ΔS ΔS = 0 S[i] += R[i] ΔS = ΔS - R[i] R[i] = 0.0 S[i] += P[i] P[i] = 0.0 P[i] - = ΔS S[i] += ΔS Total evaporation of precipitation with remainder of potential evapotranspiration False False True Covered evaporation quota with the precipitation and excedent on the surface P[i] - ΔS > 0.0 False True False True Potential soil moisture storage capacity and variation in volume of water stored Volume of water that moves in the soil pores Total infiltration of stored surface water Partial infiltration of stored surface water Direct infiltration of the rest of rainwater without runoff Direct infiltration of the rest of rainwater runoff z = 1-e -aux PEVT [i] = PEVT[i] - z EVT [i] = EVT[i] + z S[i] = S[i] - z Vi_[i] = P[i] * (1.0 - %Permeable) Vp_ [i] = P[i] * %Permeable Vg_[i] = D[i] False True aux = -S[i] / (Smax - S[i]) z = 1-e aux Dp = z * VC * Δt D[i] = Dp S[i] - = Dp S[i] > Dp Effective percolation. Increase of humidity in the soil or remnant of humidity in the soil Potential percolation i < N False True i = i + 1 Finding moisture in the soil and potential evapotranspiration Moisture reduction and increased EVT losses. Drainage volumes: Vi on impermeable soil (1st layer) Vp in permeable soil (2nd layer). Vg of groundwater (3rd layer) or according to consideration storage of humidity in soil (2nd layer). If there are more records, the process is repeated. PEVT[i] > S[i] aux= S[i] / (PEVT[i] - S[i]) False True aux= PEVT[i] / (S[i] - PEVT[i]) An additional fraction in evapotranspiration is considered WITH RUNOFF 1 2 3 4 5 6 7 8 9 10 11 1 12 Figure 2. Flowchart of the Qom model. Appl. Sci. 2020,10, 251 8 of 28 When the Qom process is finished, the Vi,Vp and Vg volumes are passed to model the hydrograms. Clark’s model is then used for modeling direct runoff of watershed rainfall response. It is based on the concept that instantaneous unit hydrograph can be derived by routing the unit excess rainfall in the form of a time-area histogram (translation) through a single linear reservoir (attenuation). As a consequence of the application of this model, a new set of decision variables is taken into account for the calibration of the model, comprising: Ki , Kp , Kg (emphasized in blue in Figure 3) and additionally the variables CT , SC and Eq can be enabled and considered (emphasized in grey); which leads to obtaining the solution vector Qe that includes the estimated water flow as formalized in Equation (4). Qei=Qii+Qpi+Qgi, 1 ≤i≤N. (4) A description of all the terms and factors used in equations of Qom model can be found in Table 1. The coding of the solution consists of a real value vector that includes each decision variable as a parameter for calibration, three decision variables namely: Rmax , Smax and VC , which represent the soil moisture estimated by Qom. Finally, in order to estimate the output flow Qe through Clark, we need 3 other decision variables: Ki , Kp , Kg , alternatively CT,Shape,Eq. Actual data values Qo and those estimated by the computation of Qe are used to calculate prediction errors in the two objective functions that will be optimized, as explained next. Table 1. Description of the terms in flowchart of Figure 2and in Equations. Term Description Abasin area (km2) Shape shape coefficient (dimensionless) + CT concentration time (h) + Eq other alternative exponetial base to e=2.711828 + %ppercent of pervious area %ipercent of impervious area Ro initial surface storage (mm) So initial soil storage (mm) Ndata register amount Rmax maximum surface storage (mm) * Smax maximum soil storage (mm) * VC volumetric conductivity coefficient (mm/h) * ivector register index Pvector of precipitation, volume expressed in height (mm) ∆ttime step inter registers (s) PEVT vector potential evapotranspiration (mm) Rvector of surface storage (mm) EVT vector of evapotranspiration (mm) Svector of soil storage (mm) Dvector of percolated volume (mm) zpercolation ratio between Rmax and Ri(dimensionless) z2soil moisture ratio between Smax and Si(dimensionless) Sp potential soil storage (mm) ∆Ssoil moisture storage between the current state and the previous (mm) Dp percolated potential volume (mm) CR soil moisture by capillary rise (mm) Vi vector of runoff volume of impervious area (mm) Vp vector of runoff volume of pervious area (mm) Vg vector of percolated volume for losses or deep runoff (mm) Ki time delay for emptying reservoir Vi (h) ] Kp time delay for emptying reservoir Vp (h) ] Kg time delay for emptying reservoir Vg (h) ] Qo vector of observed output streamflow (m3/s) Qe vector of estimated output streamflow (m3/s) Qs initial superficial streamflow (m3/s) Appl. Sci. 2020,10, 251 9 of 28 Table 1. Cont. Term Description Qd initial river base streamflow (m3/s) Qm average output streamflow observed (m3/s) Ovector of observed data Evector of estimated solution sign: * optimizing variables of Qom, ]optimizing variables of Clark, +optional decision variables. 4.2. Multi-Objective Formulation To measure the quality of solutions obtained thorough the optimization process, a widely used statistic metric in hydrology is the Nash-Sutcliffe Efficiency coefficient ( NSE ) [ 33 ], which is computed by means of Equation (5). It is a normalized metric that determines the relative magnitude of the residual variance compared to the measured data variance, evaluating the model fit primarily during high-flow periods. NSE ranges between −∞ and 1.0, with NSE = 1.0 being the optimal value. Values between 0.0 and 1.0 are generally viewed as acceptable levels of performance, whereas values < 0.0 indicate that the mean observed value is a better predictor than the simulated one, hence leading to unacceptable performance. Based on this standard measure, our multi-objective approach is focused on minimizing two objective functions f1 and f2 , by taking into account the time-series data from both, observed and estimated water flows in a given hydrologic scenario. The first objective function is based on NSE , although using the formulation f1= 1 −NSE , to promote minimization (optimum in 0.0, jMetal assumes that all the objective functions are to be minimized.) The second objective function is based on the percent bias ( PBIAS ) [ 34 ], which is formulated as f2=Abs(PBIAS) and measures the average tendency of the estimated data to be larger or smaller than their observed counterparts. PBIAS is calculated as shown in Equation (6), with optimum value in 0.0. Low-magnitude, positive, and negative values indicate, respectively, accurate model estimation, model underestimation bias and model overestimation bias. NSE =1−∑N i=1(Qei−Qoi)2 ∑N i=1(Qoi−Qm)2(5) PBIAS =∑N i=1(Qei−Qoi) ∑N i=1Qoi.(6) The problem formulation includes constraints based on the relative error in terms of the percentages RE , with AE being the absolute error between observed and estimated values according to Equation (7). A constraint considers the mass balance function MB (water balance or soil moisture accounting) that measures the river basin drainage with Equation (8) [UNESCO, 1971] for bounding the relative error MBRE . Other two constraints are defined for bounding the fitting errors in hyetograms as well as in losses and excesses volumes. In this work, solutions are constrained to f1<= 0.3, f2<= 0.3 and MBRE <=0.3 (30%). RE =AE Obeing AE =O−E(7) MB =(P−EVT −R−S−D) ∆t+ (Qs +Qd)−Qe ±AE =0. (8) 4.3. Optimization Approach Once we have formulated the calibration of the Qom model as a multi-objective optimization problem, the next step is to solve it to find an accurate approximation of the Pareto front and then take from it that solution providing the best overall result. Determining the best metaheuristic for a given problem instance is far from being a trivial process; thus, the No Free Lunch Theorem [ 35 ] states Appl. Sci. 2020,10, 251 16 of 28 5.3. Results After Choosing Solutions Once the resulting solutions are selected after calibration and validation, their corresponding decision variables are registered and analyzed in terms of the hyetographs and hydrograms, where the volumes of losses and excesses, as well as the fitting error of observed versus estimated data series are plotted. These decision variables are shown in Table 2for each problem instance, which have been used to set the Qom model. Table 2. Selected solutions from resulting RF approximations after experiments and validations. Problem Instance Rmax (mm) Smax (mm) VC (mm/h) Ki(h) Kp(h) Kg(h) 1. Asua River 16.0 131.0 0.09 8,613.0 6.6 20,446.0 2. Artibai River 12.7 45.0 0.02 653.0 7.6 75,596.0 3. San Antonio River, Loop 410 13.7 46.0 0.29 132.0 1.1 3.0 ×106 4. San Antonio River, Goliad 218.0 435.0 0.29 107.0 50.0 16,150.0 5. Juquiá River, 4F-018R 7.0 10.0 0.01 50.0 4,500.0 3,000.0 6. Bermejo River, Balapuca 17.0 349.0 0.27 1829 60.0 50,398.0 The corresponding NSE and PBIAS obtained by using these variables on the different calibration and validation data series are shown in Table 3, where additional information concerning the time period of observation, time slot of registration (10 min, 15 min, 1 h, and 24 h) and the number of registries is also provided. In general, values of NSE ≥ 0.75 and PBIAS ≤ 10% are obtained for calibration series, while NSE ≥ 0.6 and PBIAS ≤ 25% can be observed for validation series of data. In accordance with reference studies in the literature [ 39 ], where registration periods are conducted monthly, values of NSE ≥ 0.50 and PBIAS ≤ 25% can be considered as satisfactory. In the light of these results, we can claim that the prediction errors obtained by our proposed Qom model based on multi-objective calibration are more than satisfactory, since it is able to deal with complex problem instances with observation time periods significantly lower than 1 month, as done in Reference [39]. Table 3. Objective values NSE and PBIAS of selected solutions according to each data series, time period, and registry. Problem Instance Data Series Time Period # Registries NSE PBIAS 1 (Asua river) 1.1CC ∆t= 24 h 05-02-2007:10-08-2007 187 0.82 8.7% 1.2CV ∆t= 24 h 06-06-2005:08-01-2007 582 0.56 9.4% 1.3CV ∆t= 24 h 09-09-2007:11-4-2009 581 0.75 9.4% 1.4CV ∆t= 24 h 08-06-2009:30-09-2009 115 0.82 18.0% 2 (Artibai river) 2.1CV ∆t= 10 min 17-06-2002:31-12-2002 29,224 0.63 23.9% 2.2CV ∆t= 10 min 01-02-2003:04-03-2003 4588 0.68 0.5% 2.3CV ∆t= 10 min 01-04-2003:08-01-2004 40,747 0.53 12.8% 2.4CC ∆t= 10 min 17-06-2002:21-06-2002 637 0.85 22.3% 3 (San Antonio, Loop 410) 3.1CC ∆t= 15 m 15-11-2010:21-12-2010 300 0.85 24.6% 3.2CV ∆t= 15 min 17-11-2009:12-06-2012 94,315 0.43 0.5% 3.3CV ∆t= 15 min 04-06-2013:24-03-2016 98,328 0.64 53.3% 4 (San Antonio, Goliad) 4.1SC ∆t= 24 h 15-10-1998:02-11-1998 19 0.85 4.3% 4.2SV ∆t= 24 h 23-06-2002:31-08-2002 70 0.91 4.5% 4.3SV ∆t= 24 h 09-03-2007:09-03-2007 18 0.84 6.2% 5 (Juquiá, 4F-018R) 5.1CC ∆t= 1 h 09-07-2005:09-07-2005 2,351 0.79 5.0% 5.2CV ∆t= 1 h 17-10-2005:04-01-2006 1634 0.40 3.7% 5.3CV ∆t= 15 min 01-05-2013:13-06-2013 3987 0.10 6.6% 5.4CV ∆t= 15 min 11-08-2013:22-12-2013 17,551 0.60 1.05% 6 (Bermejo river, Balapuca) 6.1CC ∆t= 24 h 01-12-1971:02-07-1972 214 0.89 1.3% 6.2CV ∆t= 24 h 02-07-1972:31-08-2015 15,766 0.26 1.9% Appl. Sci. 2020,10, 251 17 of 28 More in depth, when focusing on a specific solution for a river basin, a series of hydrologic results are extracted, which leads the expert to define a hydrology engineering project. An example of these values are shown in Tables 4and 5, which include the specific mass balance obtained by Qom-Clark and estimated peak values with relative errors for problem instance 6.2CV (Rio Bermejo, Argentina), and data period 1971 09:00 to 3 April 2014 09:00 step 24.0 h. Table 4. Hydrologic results for instance 6.2CV (Cuenca del Alta del Rio Bermejo, Argentina), Qom-Clark model for period 2 February 1971 09:00 to 3 April 2014 09:00 step 24.0 h. Concept Parameter Volume (mm) (Hm3) Initial condition Rain 80,544.35 354,300.07 Superficial Storage 0.00 0.00 Soil Storage 0.00 0.00 Qom solution EVT 8855.45 38,953.54 Superficial Storage 0.00 0.00 Soil Storage 20.13 88.55 Groundwater Runoff 47,329.02 208,191.81 Superficial Runoff 24,339.74 107,066.15 Total volume estimated 80,544.35 354,300.07 Balance AE EVT 0.00 0.00 Relative Error (RE) 0.00 - Total Volume Runoff In pervious superficial - 96,359.53 In impervious superficial - 10,706.61 Off growndwater - 208,191.81 Corff. Runofff - 0.30 Comparison of water volumes Observed hydrograph 28,287.22 124,430.37 Estimated for Clark (RE = 2.05) 28,866.98 126,980.63 Estimated for Qom (RE = −13.95) 24,339.74 107,066.15 Table 5. Specific hydrograph results for instance 6.2CV (Rio Bermejo, Argentina), Clark model for period 2 February 1971 09:00 to 3 April 2014 09:00 step 24.0 h. Concept Parameter Streamflow (m3/s) Hydrograph Estimated 2726.94 Superficial impervious 41.26 Superficial pervious 2677.15 Groundwater 17.58 Flow peak of Hydrograph Observed maximum 3278.74 Relative error max peak −16.83 Observed second 2464.00 Relative error second max peak 10.67 Minimum flowrates Observed 6.58 Estimated 7.00 From a visual perspective, it is possible to check the quality of the hydrographic models obtained by plotting the hyetographs, which represent the humidity losses (HL) and excesses (HE) values from input data, together with the hydrographs, with observed (QO) and estimated (QE) outputs. In this way, Figure 8shows the resulting rainfall and streamflow by means of the hyetograph and hydrograph of the 1.1CC (Asua river basin) calibration data series, where we can observe the high correspondence between water losses/excesses and the predicted series (estimated data), which indeed fits the curve of observed data. For validation data series 1.3CV, a similar plot is made in Figure 9to better show the precision in estimated streamflow (QE), so even for external data (not used in calibration) a high prediction can be reached with the proposed approach. Appl. Sci. 2020,10, 251 18 of 28 1.1CC−Spain Period for 2007/02/06 09:00 to 2007/08/11 09:00 step 24.0h Number of days / Date 0 13 25 38 50 62 75 87 100 112 124 137 149 162 174 186 0 13 25 38 50 62 75 87 100 112 124 137 149 162 174 186 44 40 35 31 27 22 18 14 9 5 0 Height (mm) 1.00 2.19 4.34 6.48 8.63 10.77 12.92 15.06 17.21 19.35 21.50 Discharge (m3s) 2007/02/06 2007/02/18 2007/03/02 2007/03/15 2007/03/27 2007/04/09 2007/04/21 2007/05/03 2007/05/16 2007/05/28 2007/06/10 2007/06/22 2007/07/04 2007/07/17 2007/07/29 2007/08/11 Rainfall and Streamflow Reference HL HE QO QE QG Figure 8. Hyetograph and Hydrograph from calibration data series 1.1CC of Asua basin, Bizkaia, Spain. 1.3CV−Spain Period for 2007/09/10 09:00 to 2009/04/12 09:00 step 24.0h Number of days / Date 0 39 78 116 155 194 232 271 310 348 387 426 464 503 542 580 0 39 78 116 155 194 232 271 310 348 387 426 464 503 542 580 79 72 64 56 48 40 32 24 16 8 0 Height (mm) 1.00 6.78 13.51 20.24 26.97 33.70 40.42 47.15 53.88 60.61 67.34 Discharge (m3s) 2007/09/10 2007/10/18 2007/11/26 2008/01/04 2008/02/11 2008/03/21 2008/04/28 2008/06/06 2008/07/15 2008/08/23 2008/09/30 2008/11/08 2008/12/17 2009/01/24 2009/03/04 2009/04/12 Rainfall and Streamflow Reference HL HE QO QE QG Figure 9. Hyetograph and Hydrograph from validation data series 1.3CV of Asua basin, Bizkaia, Spain. A similar observation can be extracted from problem instance 2 (Artibai reiver basin), for which Figures 10 and 11 show the corresponding hyetographs and hydrographs of calibration and validation. In this case, although the amount of data used for calibration is lower than for validation, the precision power of the model is satisfactory, as changes in streamflow are usually detected and modeled for several years. With respect to the San Antonio River, two hydrometric stations were considered, in Loop 410 and Goliad. The calibration and validation hietograms and hydrograms for the two stations are represented in Figures 12–15, respectively. It is worth noting that for San Antonio and Juquiá problem instances, similar hydrographs are obtained (see Figures 16 and 17), which lead us to suggest that Qom is useful and accurate when applied to multiple and heterogeneous river basins and time periods. The hyetograph and hydrograph of the basin that belongs to Balapuca (Bermejo river) are shown in Figure 18. In this regard, an interesting comparison can be also extracted from the duration curve of the observed and estimated flows in Figure 19. It is another way to evaluate the behavior of the simulation model and the parameters, allowing to appreciate for which groups of flows better adjustments were obtained. Appl. Sci. 2020,10, 251 19 of 28 2.4CC−Spain Period for 2002/06/07 00:00 to 2002/06/11 10:00 step 0.16h Time (h) / Date 0 8 15 22 29 36 43 50 57 64 71 78 85 92 99 106 0 8 15 22 29 36 43 50 57 64 71 78 85 92 99 106 333 2222 111 0 Height (mm) 1.00 6.70 13.25 19.80 26.36 32.91 39.46 46.01 52.56 59.12 65.67 Discharge (m3s) 2002/06/07 2002/06/07 2002/06/07 2002/06/07 2002/06/08 2002/06/08 2002/06/08 2002/06/09 2002/06/09 2002/06/09 2002/06/09 2002/06/10 2002/06/10 2002/06/10 2002/06/11 2002/06/11 Rainfall and Streamflow Reference HL HE QO QE QG Figure 10. Hyetograph and Hydrograph from calibration data series 2.4CC of Artibai barin, Bizkaia, Spain. 2.3CV−Spain Period for 2003/04/01 00:00 to 2004/01/08 23:00 step 0.16h Time (h) / Date 0 453 906 1359 1811 2264 2717 3170 3622 4075 4528 4981 5433 5886 6339 6791 0 453 906 1359 1811 2264 2717 3170 3622 4075 4528 4981 5433 5886 6339 6791 88 7 6 5 44 3 2 1 0 Height (mm) 1.00 10.35 20.67 30.98 41.30 51.61 61.93 72.25 82.56 92.88 103.19 Discharge (m3s) 2003/04/01 2003/04/19 2003/05/08 2003/05/27 2003/06/15 2003/07/04 2003/07/23 2003/08/11 2003/08/29 2003/09/17 2003/10/06 2003/10/25 2003/11/13 2003/12/02 2003/12/21 2004/01/08 Rainfall and Streamflow Reference HL HE QO QE QG Figure 11. Hyetograph and Hydrograph from validation data series 2.3CV of Artibai barin, Bizkaia, Spain. 3.1CC−USA Period for 2010/02/07 02:15 to 2010/02/09 09:15 step 0.25h Time (h) / Date 0 4 8 11 15 19 22 26 30 33 37 41 44 48 52 55 0 4 8 11 15 19 22 26 30 33 37 41 44 48 52 55 24 22 19 17 15 12 10 8 5 3 0 Height (mm) 1.00 77.63 154.87 232.10 309.33 386.57 463.80 541.03 618.26 695.50 772.73 Discharge (m3s) 2010/02/07 2010/02/07 2010/02/07 2010/02/07 2010/02/07 2010/02/07 2010/02/08 2010/02/08 2010/02/08 2010/02/08 2010/02/08 2010/02/08 2010/02/08 2010/02/09 2010/02/09 2010/02/09 Rainfall and Streamflow Reference HL HE QO QE QG Figure 12. Hyetograph and Hydrograph from validation data series 3.1CC of San Antonio river, at Loop 410, USA. Appl. Sci. 2020,10, 251 20 of 28 3.3CV−USA Period for 2013/06/04 00:00 to 2016/03/24 06:30 step 0.25h Time (h) / Date 0 1639 3278 4917 6556 8195 9833 11472 13111 14750 16389 18028 19666 21305 22944 24583 0 1639 3278 4917 6556 8195 9833 11472 13111 14750 16389 18028 19666 21305 22944 24583 28 25 22 20 17 14 11 9 6 3 0 Height (mm) 1.00 144.99 289.94 434.89 579.84 724.79 869.73 1014.68 1159.63 1304.58 1449.53 Discharge (m3s) 2013/06/04 2013/08/11 2013/10/18 2013/12/25 2014/03/04 2014/05/11 2014/07/18 2014/09/24 2014/12/02 2015/02/08 2015/04/17 2015/06/25 2015/09/01 2015/11/08 2016/01/15 2016/03/24 Rainfall and Streamflow Reference HL HE QO QE QG Figure 13. Hyetograph and Hydrograph from validation data series 3.3CV of San Antonio river, at Loop 410, USA. 3.3CV−USA Period for 2013/06/04 00:00 to 2016/03/24 06:30 step 0.25h Time (h) / Date 16962 17028 17094 17159 17225 17290 17356 17421 17487 17552 17618 17684 17749 17815 17880 17946 16962 17028 17094 17159 17225 17290 17356 17421 17487 17552 17618 17684 17749 17815 17880 17946 26 23 21 18 16 13 11 8 6 3 0 Height (mm) 1.00 144.99 289.94 434.89 579.84 724.79 869.73 1014.68 1159.63 1304.58 1449.53 Discharge (m3s) 2015/05/11 2015/05/14 2015/05/17 2015/05/19 2015/05/22 2015/05/25 2015/05/28 2015/05/30 2015/06/02 2015/06/05 2015/06/08 2015/06/10 2015/06/13 2015/06/16 2015/06/18 2015/06/21 Rainfall and Streamflow Reference HL HE QO QE QG Figure 14. Hyetograph and Hydrograph from validation data series 3.3CV-Zoom of San Antonio river, at Loop 410, USA. 4.2SC−USA Period for 2002/06/04 09:00 to 2002/08/12 09:00 step 24.0h Number of days / Date 0 5 10 14 19 23 28 33 37 42 46 51 56 60 65 69 0 5 10 14 19 23 28 33 37 42 46 51 56 60 65 69 137 123 110 96 82 69 55 41 28 14 0 Height (mm) 3.00 177.70 353.03 528.36 703.69 879.01 1054.34 1229.67 1405.00 1580.32 1755.65 Discharge (m3s) 2002/06/04 2002/06/08 2002/06/13 2002/06/17 2002/06/22 2002/06/27 2002/07/01 2002/07/06 2002/07/10 2002/07/15 2002/07/20 2002/07/24 2002/07/29 2002/08/02 2002/08/07 2002/08/12 Rainfall and Streamflow Reference HL HE QO QE QG (a)4-2SC-HH Figure 15. Cont. Appl. Sci. 2020,10, 251 21 of 28 4.1SV−USA Period for 1997/12/06 09:00 to 1997/12/24 09:00 step 24.0h Number of days / Date 0 2 3 4 5 6 8 9 10 11 12 14 15 16 17 18 0 2 3 4 5 6 8 9 10 11 12 14 15 16 17 18 170 153 136 119 102 85 68 51 34 17 0 Height (mm) 5.00 161.81 319.40 476.98 634.57 792.16 949.74 1107.33 1264.92 1422.50 1580.09 Discharge (m3s) 1997/12/061997/12/071997/12/081997/12/091997/12/10 1997/12/121997/12/131997/12/141997/12/151997/12/161997/12/17 1997/12/191997/12/201997/12/211997/12/221997/12/23 Rainfall and Streamflow Reference HL HE QO QE QG (b)4-1SV-HH 4.3SV−USA Period for 2007/03/10 09:00 to 2007/03/27 09:00 step 24.0h Number of days / Date 0 2 3 4 5 6 7 8 10 11 12 13 14 15 16 17 0 2 3 4 5 6 7 8 10 11 12 13 14 15 16 17 80 72 64 56 48 40 32 24 16 8 0 Height (mm) 6.00 53.05 100.15 147.25 194.35 241.45 288.55 335.65 382.75 429.85 476.95 Discharge (m3s) 2007/03/10 2007/03/11 2007/03/12 2007/03/13 2007/03/14 2007/03/15 2007/03/16 2007/03/17 2007/03/19 2007/03/20 2007/03/21 2007/03/22 2007/03/23 2007/03/24 2007/03/25 2007/03/26 Rainfall and Streamflow Reference HL HE QO QE QG (c)4-3SV-HH Figure 15. Hyetograph and Hydrograph from calibration 4-2SC-HH and idem for validation 4-1SV-HH and 4-3SV-HH data series. 5.1CC−Brazil Period for 2005/07/09 00:00 to 2005/10/14 22:00 step 1.0h Time (h) / Date 0 157 314 470 627 784 940 1097 1254 1410 1567 1724 1880 2037 2194 2350 0 157 314 470 627 784 940 1097 1254 1410 1567 1724 1880 2037 2194 2350 16 15 13 12 10 8 7 5 4 2 0 Height (mm) 63.00 78.20 93.74 109.27 124.80 140.33 155.87 171.40 186.93 202.47 218.00 Discharge (m3s) 2005/07/09 2005/07/15 2005/07/22 2005/07/28 2005/08/04 2005/08/10 2005/08/17 2005/08/23 2005/08/30 2005/09/05 2005/09/12 2005/09/18 2005/09/25 2005/10/01 2005/10/08 2005/10/14 Rainfall and Streamflow Reference HL HE QO QE QG Figure 16. Hyetograph and Hydrograph from validation data series 5.1CC of Juquiá river at 4F-018R. Appl. Sci. 2020,10, 251 22 of 28 5.4CV−Brazil Period for 2013/08/01 17:00 to 2014/01/31 12:30 step 0.25h Time (h) / Date 0 293 585 878 1170 1463 1755 2048 2340 2633 2925 3218 3510 3803 4095 4388 0 293 585 878 1170 1463 1755 2048 2340 2633 2925 3218 3510 3803 4095 4388 24 21 19 17 14 12 10 7 5 3 0 Height (mm) 59.00 83.19 108.13 133.06 158.00 182.93 207.86 232.80 257.73 282.67 307.60 Discharge (m3s) 2013/08/01 2013/08/13 2013/08/26 2013/09/07 2013/09/19 2013/10/01 2013/10/13 2013/10/26 2013/11/07 2013/11/19 2013/12/01 2013/12/13 2013/12/25 2014/01/07 2014/01/19 2014/01/31 Rainfall and Streamflow Reference HL HE QO QE QG Figure 17. Hyetograph and Hydrograph from validation data series 5.4CV of Juquiá river at 4F-018R. 6.2CV−Argentina−Bermejo−Balapuca Period for 1971/02/02 09:00 to 2014/04/03 09:00 step 24.0h Number of days / Date 0 1052 2103 3154 4205 5256 6307 7358 8409 9460 10511 11562 12613 13664 14715 15766 0 1052 2103 3154 4205 5256 6307 7358 8409 9460 10511 11562 12613 13664 14715 15766 250 225 200 175 150 125 100 75 50 25 0 Height (mm) 7.00 333.79 661.01 988.23 1315.44 1642.66 1969.87 2297.09 2624.31 2951.52 3278.74 Discharge (m3s) 1971/02/02 1973/12/19 1976/11/04 1979/09/21 1982/08/07 1985/06/23 1988/05/09 1991/03/26 1994/02/09 1996/12/26 1999/11/12 2002/09/28 2005/08/14 2008/06/30 2011/05/17 2014/04/02 Rainfall and Streamflow Reference HL HE QO QE QG Figure 18. Hyetograph and Hydrograph from validation data series 6.2CV of Balapuca, Bermejo river, Bolivia-Argentina. % Time Discharge (m3s) 0.0 5.0 10.0 20.0 30.0 40.0 50.0 60.0 70.0 80.0 90.0 100.0 3278.7 1532 1102 711 540 346 227 135 85 53 34 24 18 14 12 6.6 Duration Curve 3278.7 1532 1102 711 540 346 227 135 85 53 34 24 18 14 12 6.6 Reference Observed Estimated Figure 19. Duration Curve from validation data series 6.2CV of Balapuca, Bermejo river, Bolivia-Argentina. Appl. Sci. 2020,10, 251 23 of 28 Resulting hydrographs in Figures 8and 19 reveal that observed streamflows are properly fitted by estimated ones in different problem instances and time periods. 6. Discussion The main goal driving us in this work has been to provide a practical tool for hydrologic prediction based on the combination of two components—the Qom model and its automatic calibration using multi-objective optimization algorithms. The results reported in the previous section indicate that the Qom model is stable for a set of heterogeneous scenarios with different conditions and it does not produce large changes to small variations of the parameters. Therefore, it can be considered an adequate model for its use in the basin hydrological risk estimation, in basin climate change assessments, and in changes by use or soil practices. A question that can arise is why have we proposed a new hydrologic model instead of using an existing one? The reason is related to the expertise of the paper authors, comprising hydraulic and computer science engineers. We are the creators of the jMetal software framework for multi-objective optimization metaheuristics. jMetal includes a large amount of algorithms representative of the state of the art but, as pointed out in Section 4.3, choosing which is the most promising one to solve a given optimization problem is not a easy task, in particular for experts in the application domain, which are usually not experts in the optimization techniques. In this sense, a solution selection strategy is also suggested for guiding the expert in hydrology in the decision making process, which comprises the use of independent data series for calibration, as well as for validation, each of them involving different time periods of a given problem instance. This leads the proposal to select successful solutions that could be generated by calibration algorithms with low contribution in the Pareto reference front approximation, but with very stable results in terms of practical validation. Our starting point was to use jMetal for the calibration of a hydrologic model; although this algorithmic library is implemented in the Java programming language, it hinders using existing models; for example, the Sacramental Soil Moisture Accounting (SAC-SMA) model is available in FORTRAN, Matlab, and R, but to the best of our knowledge there is no a Java version. Consequently, from our experiences in water resources we decided to define a new model—Qom—based on a reduced number of parameters to be adjusted, to implement it in Java, and integrating it into jMetal to use a multi-objective approach for automatic calibration, all aimed at providing a generic and accurate software tool that is freely available to the community. Our objective in this paper has not been to compare Qom with existing models, nor to compare the automatic calibration scheme with existing strategies. These issues remain open and are included as lines of future work in the next section. 7. Conclusions In this work, a novel hydrologic model called Qom is proposed to efficiently separate and quantify the volumes of losses and excesses of the rainwater falling in a hydrographic basin. In combination with it, an evolutionary multi-objective approach is used for parameter calibration focusing on indicators to qualify the adjustments between the observed and estimated hydrograms. A number of multi-objective algorithms in the scope of the jMetal framework are used to deal with the hydrologic model as a black-box optimization problem, with a vector of decision variables as inputs and objective function values as outputs. The Qom model, featuring a reduced set of parameters, together with the possibility of using existing multi-objective optimizers without the need to tune each of them, just combining the solutions they produce, make this proposal of special interest for both academy and industry communities. For testing purposes, six realistic hydrographic scenarios, comprising river basins located in Spain, USA, Brazil and Argentina, have been modeled and evaluated with different size conditions, climates, topographies, heterogeneous soils, and with series of very short time periods of 10, 15, 30 Appl. Sci. 2020,10, 251 24 of 28 and 60 min and 24 h. The results indicate that the calibration of Qom of each problem allow to predict future storm flows as it is possible of accurately reproduce the phenomenon rain-runoff (cause-effect), thus fulfilling the goal of our work. In general, the NSE and PBIAS obtained by Qom involved highly successful values of roughly 75% and 10% (respectively) in calibration, as well as in validation. This could be explained in terms of low errors of the predicted series with regards to observed ones, which can be classified as very satisfactory according to Reference [ 39 ]. Moreover, these results could even be improved if larger time periods of registration are considered, as we measure changing conditions from 10 min to 24 h, whereas references in the literature worked within monthly periods [39]. As lines of future work, an open issue to explore is the analysis of the Qom model when compared with other hydrological models, such as SWAT, VIC or SAC-SMA. Our proposal to combine different multi-objective metaheuristics can be time consuming, as a large number of independent runs of configurations of the algorithm/problem must be carried out. It would be more efficient to explore the use of automatic parameter tuning tools to determine if a single algorithm could be good enough. In this sense, a further step would be to self-generate an ad-hoc metaheuristic suited for the calibration of Qom. In both cases, the idea would always be to use automatic tools to avoid the expert to deal with algorithmic issues. The performance of the resulting method would be validated by comparing it with other parameter estimation techniques. Author Contributions: G.R.Z. has contributed to the conceptualization, methodology, software development and validation; J.G.-N. has contributed to methodology, formal analysis, writing—original draft preparation, writing—review, and editing and visualization; A.J.N. has contributed to supervision, software development, writing—review, and editing. All authors have read and agreed to the published version of the manuscript. Funding: This work has been partially funded by Grants TIN2017-86049-R (Spanish Ministry of Education and Science). José García Nieto is the recipient of a Post-Doctoral fellowship of “Captación de Talento para la Investigación” Plan Propio at University of Málaga. Conflicts of Interest: The authors declare no conflict of interest. Appendix A. River Basins Description A description of the main features of the studied river basins is included next: • Asua river basin, Bizkaia, Spain. It is located in the north of the Iberian Peninsula. The climate is temperate, oceanic, with frequent and abundant rainfall, especially in autumn and winter, with an annual average of 1200 mm. It is a rural basin comprising a drainage area equal to 50.8 km 2 , corresponding to the Derio (meteorological) and San Groniz (capacity) stations, with slopes not exceeding 1%. It is one of the most open valleys in Bizkaia. It is bordered by mountains of low altitude without exceeding 360 m in height. The lands through which the courses of this basin run are mainly made up of marl and limestone, crossing the main river alluvial lands from the middle stretch to the mouth.The shape can resemble a rectangular figure of proportions 1.5 to 1 with respect to length and width. • Artibai, Bizkaia, Spain. The Artibai river is about 20 km long, and extends in a S-NE direction. Its origin are two groups of streams from mountains between 1029 m 793 m altitude. The geological substratum of the basin is predominantly limestone at the head and later sandstone and clay. The fluvial bed is stony with a predominance of blocks or boulders both in the streams and in the main channel. The shape can resemble a rectangular figure of proportions 2 to 1 with respect to length and width. • San Antonio river sub-basin, station Loop 410, Texas, USA. It is a sub-urban basin, having an average slope of 0.2%, and a climate that alternates between dry and humid, summers are hot and winters vary from mild to cold; the average annual rainfall is 738 mm. It can resemble a rectangular shape of proportions 2 to 1 in relation to length and width. • San Antonio river basin, station Goliad, Texas, USA. The morphometry of the San Antonio River basin includes the contribution to the Goliad station. The basin has been considered regional as a Appl. Sci. 2020,10, 251 25 of 28 single unit of drainage area, most of its surface being permeable. It has slopes from 0.2% to 3%, has an elongated shape from northwest to southeast, the climate is varied, alternating between humid and predominantly dry. Being an extensive drainage area of the regional type, the area of incidence of storms has been distributed spatially. The shape is very elongated with a ratio of 5 to 1 with respect to the width. • Juquiá river basin, station 4F-018R, Ribeira Do Iguapé, São Pablo, Brasil: We have considered the sub-basin of the Juquiá River in the hydrometeorological station 4F-018R, which includes a drainage area equal to 4360.0 km 2 , code 81679000, basin of the Ribeira Do Iguapé River, São Pablo, Southeast Atlantic Hydrographic Region, Brazil. The topography is mountainous with the characteristics of a somewhat flowing and turbulent river. The region has a hot tropical climate, with high temperatures in summer. The high rainy season occurs in winter with mild temperatures. Formed by two sub-basins, its shape is elongated at the headwaters and widens at the exit. • Bermejo river upper basin, station Balapuca, Argentina-Bolivia. The basin is located in the north of Argentina, province of Salta, sharing a part with Bolivia; the sub-basin belongs to the Bermejo river basin. It has the mountainous ecosystem of the Andes Mountains with an average slope of 35%. The hydrological regime of the rivers is purely pluvial with a well-defined seasonal variety, characterized by a period of important flows in the rainy season. The shape of the basin can be assimilated to a square shape of relation 1 to 1 with respect to the length and width. The basins are depicted in Figure A1. Figure A1. Cont.