A constraint-based search algorithm for parameter identiﬁcation of environmental models

. Many environmental systems models, such as conceptual rainfall-runoff models, rely on model calibration for parameter identiﬁcation. For this, an observed output time series (such as runoff) is needed, but frequently not available (e.g., when making predictions in ungauged basins). In this study, we provide an alternative approach for parameter identiﬁcation using constraints based on two types of re-strictions derived from prior (or expert) knowledge. The ﬁrst, called parameter constraints , restricts the solution space based on realistic relationships that must hold between the different model parameters while the second, called process constraints requires that additional realism relationships between the


Introduction
Environmental systems models, such as conceptual rainfallrunoff (CRR) models, are abstract simplifications of real system behavior. Often, the parameters in such models cannot be specified through direct measurements of physical properties of the system. Further, even when a parameter is related to measurable quantities, its value in the model typically represents an integrated value over a much larger scale than the measurement scale. For this reason, such models typically rely upon calibration (tuning to match system input-output behavior for a given historical data period) to ensure satisfactory performance when applied to specific hydrological systems of interest (Wheater et al., 1993;Beven, 2001).
In the case of CRR, parameter values are typically specified through a process of calibration that seeks to match the model runoff simulations to observed hydrographs by the use of an objective function (e.g., root mean square error, RMSE). Expert knowledge is brought to bear implicitly, by the prior selection of the model and the specification of parameter ranges that define the prior parameter space. Recently, several studies have tested strategies that relate the parameters of CRR models to catchment physical characteristics (Koren et al., 2000(Koren et al., , 2003Anderson et al., 2006;Yadav et al., 2007;Pokhrel et al., 2008Pokhrel et al., , 2012Kling and Gupta, 2009;Hrachowitz et al., 2013). The general picture that emerges from these studies is that exploiting expert knowledge (by somehow imposing more rigorous constraints 4862 S. Gharari et al.: Constraint-based parameter identification on the parameters) has the potential to result in more realistic models (Martinez and Gupta, 2011) that in several cases has practical benefits.
For example, Pokhrel et al. (2008Pokhrel et al. ( , 2012 linked the parameters of a spatially distributed model to catchment physical characteristics via a set of regularization (or regionalization) relationships, and thereby converting the original high-dimensional parameter space to a much reduced set of "super-parameters", which then leads to a dramatic simplification for a calibration problem. Similarly, Merz and Blöschl (2004), Kling and Gupta (2009) and Yadav et al. (2007), amongst others, investigated explicit links between catchment characteristics and the parameters of a simple lumped conceptual model; they concluded, however, that such relationships are difficult to establish and may not be often possible given the available data. Other studies have used a comparison of catchment characteristics based on similarities between catchment responses to constrain parameter values. For example, Zhang et al. (2008), imposed a set of three constraints to infer the runoff characteristics of catchments, following the concept of regionalization of hydrological signatures (initially developed by Yadav et al., 2007). Recently, Kapangaziwiri et al. (2012) constrained the Pitman monthly rainfall runoff model (Hughes et al., 2006) based on a regionalization of runoff signatures. Perrin et al. (2008) proposed a method called discrete parameterization based on the use of parameter sets compiled a priori via calibration to other catchments. Their approach "abandon(s) the idea of searching for an optimum parameter set in the continuous, n-dimensional parameter space" and instead "limit(s) the calibration process to a search within a finite collection (a library) of predefined parameter sets". More recently, Samaniego et al. (2010) and Kumar et al. (2010Kumar et al. ( , 2013 demonstrated that a multi-scale approach to parameter regionalization can provide consistent model performance for both gauged and ungauged catchments. In a complementary direction, the use of multiple objective functions or multiple system responses for calibration (Gupta et al., 1998) has been shown to result in more realistic parameter sets that achieve improved simulations of system dynamics. The multi-objective approach seeks to identify parameter sets that simultaneously provide optimal performance for different aspects of system response (Gupta et al., 1998;Boyle et al., 2000Boyle et al., , 2001. This can include constraining the model to reproduce multiple system fluxes and state variables such as runoff, evaporation, groundwater levels or tracer concentrations (e.g., Gupta et al., 1999;Bastidas et al., 1999;Freer et al., 2002;McDonnell, 2002, 2013;Khu and Madsen, 2005;Fenicia et al., 2008;Winsemius et al., 2008;Birkel et al., 2011;Hrachowitz et al., 2013).
While the aforementioned studies have demonstrated that incorporation of expert and a priori knowledge can help improve the realism of models, to our knowledge, no systematic strategy has been presented in the literature for constraining the model parameters to be consistent with the patchy understanding of a modeler regarding how the real system might work. Part of the difficulty in doing this is that expert knowledge may not directly translate to quantifiable relationships (e.g., between catchment physical characteristics and model parameters) rather, it may consist of conceptual understanding about consistency relationships that must exist among modeled state variables and/or fluxes, as well as, among various model parameters. For example, the geology of a given catchment may suggest that the catchment response during intense rainfall events is characterized by a slow responding groundwater component accompanied by fast responding Hortonian overland flow. In this case, any model result that implies peak flows are composed of a strong groundwater response should be discarded or should be given low importance.
An example of such an approach toward modeling was mentioned by Götzinger and Bárdossy (2007), where they impose the Lipschitz and monotonic conditions to avoid the abrupt jump in soil moisture values for the neighboring cells of a distributed model based on the physical premises that such jumps are numerical artifacts and hydrologically unrealistic. Such kinds of information, which are not explicitly provided during model calibration, act as constraints to limit the feasible extent of the model parameter space, thus resulting in physically more meaningful model simulations. As pointed out by Efstratiadis and Koutsoyiannis (2010) "It also offers a means to partially handle the huge uncertainty resulting from the complexity of model parameterizations in contrast to data scarcity, which is a global engineering problem that is getting increasingly severe. Actual research should provide more guidance on the effective combination of statistical and expert-based evaluation procedures." In this study, we present a constraint-based algorithm to limit the feasible parameter space of a conceptual hydrologic model, based on relational constraints inferred from expert knowledge regarding plausible catchment behavior. The approach is applicable to both lumped/semi-distributed and spatially distributed catchment models.

Constraints in environmental models
In the companion paper, Gharari et al. (2014), we test the importance of constraining the feasible parameter space. The case study focuses on the meso-scale catchment of the Wark in Luxembourg covering an area of approximately 82 km 2 . Three different models with varying (spatial) complexity are developed to capture hydrological processes across different delineated landscapes, wetland, hillslope and plateau, following the work by Gharari et al. (2011). The first model, FLEX A , is a lumped model which considers the whole catchment as a single unit. A more complex model, FLEX B , takes into account wetlands (i.e., riparian zones), and considers the rest of the catchment as another unit. The most complex model of all, FLEX C , distinguishes between the aforementioned landscape classes (wetland, hillslope and plateau).
A set of parameter and process constraints is defined for each of the three landscape entities, based on the expert knowledge. Parameter constraints are considered to be a priori because they can be imposed without actually running the model, while process constraints can only be imposed after a model is run with selected parameter sets. The number of constraints may vary from model to model, and here depends on the model complexity. For example all of the predefined constraints (see the following section for more detail) can be applied to the most complex model, FLEX C , while only a subset of constraints can be applied to FLEX B and just few constraints to FLEX A . Using the proposed search algorithm, explained later in this study, we define the behavioral set of parameters that satisfy all the parameter and process constraints. These resulting parameter sets are termed as constrained but un-calibrated. The performance of the constrained but un-calibrated parameter sets for the three models could then be evaluated based on Nash-Sutcliffe efficiency, E NS (Nash and Sutcliffe, 1970). Figure 1 illustrates the 95 % uncertainty interval of the model simulations generated using the constrained but un-calibrated parameter sets, and reports the median Nash-Sutcliffe efficiency (E NS,median ) for three model structures.
The results show that the constrained and calibrated FLEX C , the most complex model, has lower predictive uncertainty compared to a (loosely) constrained and calibratedlumped model, FLEX A . This result contradicts to what appears to be a general belief about the higher predictive uncertainty of complex models, and shows that if a complex model is properly constrained then it can perform well during the evaluation period while providing even lower predictive uncertainty than a simpler model. Result of the recent study by Gao et al. (2014) had shown that a more complex, but simultaneously constrained and calibrated model, can provide better predictions of extreme events outside of the calibration period as well as capturing the flow duration curve in nested catchments for the large scale Heihe basin in China. Meanwhile, Hrachowitz et al. (2014) reported similar results for a small-scale catchment in France, using 11 different model structures and 20 different runoff signatures. In these studies, the introduction of constraints was found to be crucial for ensuring better model performance, particularly outside of a calibration period.
As mentioned earlier, two types of model constraints can be distinguished: a priori constraints applicable to model parameters (i.e., parameter constraints) and a posteriori constraints on model states and fluxes (i.e., process constraints; e.g., Bulygina and Gupta, 2009. These two types of constraints are elaborated upon in the following section.

Parameter constraints
Parameter constraints provide information regarding the relationships between parameters of the same process that correspond to different spatial components (or units or grid cells) of a (semi-) distributed model. Such constraints can be expressed by inequality (or equality) constraints; for example where G 1 is a dimensionless constant where G 2 has the same units as A 1 and A 2 .
As a simple illustration of this concept, the maximum interception capacity of a forested area (I max,forest ) can typically be assumed to be larger than the maximum interception capacity of a grassland area (I max,grass ).

Process constraints
Process constraints provide comparative information regarding the fluxes and/or states of a model at each time step, or integrated over some specific time period. Examples of such constraints include the following: where G 3 and G 4 are dimensionless constants.
As an illustration, one can compare the transpiration fluxes from different spatial entities of a (semi-) distributed model. For example for two regions having similar soil type and aspect, the region with smaller normalized difference vegetation index (NDVI) is expected to transpire at a lower rate.
It is worth noting that in either case, parameter sets that satisfy the constraints are not conditional on the information provided by observations (or measurements) of the output response of the system (e.g., the runoff hydrograph), and these can therefore be determined without resorting to model calibration. Moreover, parameter sets that satisfy all of the constraints (within some acceptable range) can provide insights into how the real system can be expected to behave, assuming that it corresponds to the expert's perception of realistic (behavioral) system properties and dynamics.
Unfortunately, the use of available evolutionary algorithms to search for parameter sets that satisfy such constraints is complicated by the non-convex and potentially non-continuous parameter search space that results. In the following section, we propose a stepwise search strategy that can be used to identify behavioral parameter sets fulfilling the constraints imposed by expert knowledge.

Methodology and algorithm -constraint based search (CBS)
The constraint-based search (CBS) algorithm is based on a simple parameter sampling approach to identify the parts of the feasible parameter space that satisfy the set of constraints as discussed in the previous section. At each step, the algorithm tries to generate new parameter sets that satisfy the parameter constraints, while only violating the process constraints to an acceptable level. This level could be set-up based on the desired model performance. The process of search continues until all of the parameter sets generated properly satisfy the set of imposed process constraints. In the following description M refers to the total number of process constraints and m, m ∈ {1, . . ., M}, is an index indicating how many of the process constraints are satisfied by a given parameter set; for example, if a parameter set satisfies two process constraints then m will be equal to 2. The algorithm ultimately generates a set P of n number of behavioral parameter sets that satisfy all of the parameter and process constraints (i.e., m = M for all of members of P ).
The search algorithm is as follows; let C be a counter which identifies the minimum number of process constraints to be satisfied (at a given stage of the search) by the parameter sets considered to be behavioral, N is the initial sample size, and S is the sample size at each stage. The user-specified parameter f represents the memory in CBS corresponding to which of the existing parameter sets should be used to generate the new parameter sets at each stage; specifically, it indicates the extent to which the algorithm is allowed to relax the acceptable minimum number of process constraints (C) when considering parameter sets for use in generating the new samples 1 .
-Step 1: generate N random samples (parameter sets) across the entire feasible parameter space using uniform prior distributions. - Step 2: evaluate the parameter constraints and identify the K 1 samples that satisfy them. - Step 3: run the model K 1 times, one time for each of the samples identified in step 2, evaluate the M process constraints for each sample and assign a value of m to each parameter set corresponding to the number of process constraints satisfied.

-
Step 4: place the K 2 samples that satisfy C or more process constraints in set P , and the K 3 samples that satisfy exactly C − f − 1 to C − 1 process constraints in set P . Discard samples that satisfy C − f − 2 or fewer process constraints. - Step 5: use the members of sets P and P to generate S new samples by applying each of the three Monte Carlo based rules below to generate S/3 of the samples, where θ new is the newly generated sample. θ P and θ P are samples selected randomly from sets P and P respectively and α is a random value between 0 and 1 which differs for each newly generated sample. Figure 2 shows a graphical illustration of these rules.
Step 6: discard all existing members of set P which satisfy exactly C − f − 1 process constraints 2 .

-
Step 7: increase C by one and return to step 2. Repeat this process until C becomes equal to the total number of process constraints (i.e., C = M), or the maximum number of model evaluation is exceeded, or the size of P is satisfactory.
In step 5, the selection of parameter sets from P and P can also be based on the number of satisfied process constraints or any other distribution different from a uniform distribution.
Note that any member of set P is within the space marked by members of set P . Using members of P to generate new parameter sets (step 5) helps to identify the boundary between the parameter space that satisfies exactly C − f − 1 to C − 1 (set P ) and C or more (set P ) constraints. The intention is to obtain a diverse parameter representation for set P by including the set P (Fig. 2).
The final set P contains parameters that simultaneously satisfy all of the parameter and process constraints. These parameter sets can be referred to as constrained but uncalibrated, as they are not calibrated to match observed data sets or target variables (e.g., discharge time-series). Figure 3 presents a graphical illustration of the steps of CBS to find set constrained but un-calibrated parameter sets.
Note that the set P can also be used to constrain a search for optimal parameter sets within this space of constrained but un-calibrated parameter sets. This can be easily achieved by evaluating them based on model performance in regard to a target variable (e.g., observed runoff). The set P can be used as an initial sample for any evolutionary algorithm. In 2 In the case that f is set to 0 there will not be any existing parameter sets in P (i.e., P = ).
Eq. 6 Eq. 7 Eq. 7 Eq. 8 P Pʹ P P θ 2 θ θ 1 Figure 2. A conceptual illustration of possible positions of newly generated parameter sets based on parameter sets randomly drawn from P and P for a two dimensional parameter space. The area indicated in yellow represents the set P that satisfies exactly C − f − 1 to C − 1 process based constraints. The area indicated in green represents the set of P that satisfies C or more process constraints. The circles indicate randomly selected parameter sets drawn from the sets P or P . Different line styles indicate different parameter generation rules (corresponding to the different equations). Solid lines represent the first rule, where the parameter sets are randomly selected from set P (Eq. 6); dashed lines illustrate the second rule, where one parameter set is randomly selected from P and one from P (Eq. 7); the dash-dot line represents the third rule, where both randomly selected parameter sets are selected from set P (Eq. 8). Note that due to possible non-convexity of sets of P and P the newly generated parameter sets based on the three rules can be outside of sets P and P . this case, any new parameter set generated by an optimization algorithm would need to be checked for both parameter and process constraints and would be retained only if they satisfy the entire set of constraints.

Case study
A synthetic case study was designed to illustrate the efficiency of the proposed constraint-based search algorithm. First, the lumped model FLEX A (Fig. 4; see companion paper Gharari et al., 2014) was calibrated to runoff data from the Wark catchment for the year 2002-2005, using year 2001 as warm up period. A multi-criteria calibration was conducted without imposing any parameter or process constraints, based on three objective functions: (i) the Nash-Sutcliffe efficiency of the flows (E NS ), (ii) the Nash-Sutcliffe efficiency of the logarithm of the flows (E NS,log ) and (iii) the Nash-Sutcliffe efficiency of the flow duration curve (E NS,FDC ). From the parameter sets belonging to the Pareto front, the parameter set having minimum distance to the origin was taken as the best behavioral parameter set (θ best ).   Figure 4. The models structure of FLEX A . Precipitation is indicated by P ; I and T represent components of evaporation which are interception and transpiration respectively. Q f and Q s indicate the fast and slow components that make the total model output (Q m ). The water balance equations and constitutive functions are presented in the companion paper .
A set of (parameter and process) constraints was then designed based on the model simulations (fluxes and states) provided by this best parameter set (θ best ). Since the constraints were constructed so that the best performing parameter set (θ best ) is behavioral, we are guaranteed the existence of feasible parameter sets that satisfy all of those constraints. Table 1 summarizes the parameter and process constraints designed in this manner; note that the value for each constraint corresponding to θ best is mentioned, along with the acceptable range limits for each constraint to be used during the search.
The proposed algorithm was then applied to search for parameter sets that satisfy the constraints mentioned in Table 1. The initial sample size (N) was taken as 50 000 parameter sets. S, f and C were set to be 5000, 1 and 2, respectively. The search was terminated when the number of generated parameter sets satisfying all of the parameter and process constraints exceeded 8000.
The parameter sets identified by the proposed search algorithm were then compared with the best parameter set (θ best ). Figure 5 shows the results using a normalized parameter plot (Gupta et al., 1998), where the red trajectory indicates θ best , and the hue of the grey-to-black lines indicates the number of process constraints satisfied by a given parameter set (the darker the line, the higher number of process constraints satisfied). The result clearly shows that with the increasing number of process constraints satisfied, the feasible solution space converges towards the best parameter set (θ best ).
A similar comparison (Fig. 6) was conducted for the modeled hydrographs associated with these parameter sets. The hydrograph for the best performing parameter set (θ best ) is indicated in red, and the simulated hydrographs for different parameter sets generated by the proposed search algorithm are depicted in different shades ranging from grey-toblack (darker colors indicate that larger numbers of process constraints are satisfied). Clearly, with increasing number of Table 1. The "true" values for each of the design constraints (corresponding to θ best .) along with the acceptable limits assigned to each constraint for use during the search. Dry and wet periods are defined as May to October and November to April respectively. Events in which the discharge increases with a rate of more than 0.2 mm (3 h) −2 are defined as peak flows.

No
Quantification of designed constraints Acceptable range for each designed based on θ best constraints Parameter constraint 1 Process constraints < 0.6 * t indicates the number of time steps over the entire period of simulation. * * f (S u ) = 1, 0.5 S u,max ≤ S u 0, S u < 0.5 S u,max satisfied process constraints, the hydrographs progressively approach the one simulated using the θ best .
Overall the search algorithm generated and evaluated 102 106 parameter sets to find 8000 feasible solutions that satisfy all of the nine constraints imposed, which corresponds to approximately 8 % efficiency. In comparison, when a conventional Monte-Carlo sampling approach was applied using 102 106 samples, none of the samples were found to be able to satisfy all of the constraints, while only two of the samples were able to satisfy at least 7 of the 9 process constraintsimposed. Clearly, the proposed search algorithm is capable of relatively rapid convergence towards the region of the parameter space where all of the constraints are satisfied. Figure 7 illustrates how quickly the search algorithm is able to locate the behavioral parameter sets. It depicts the percentages of generated samples satisfying a given number of process constraints at each step of the search. Darker colors are used to indicate the proportions of parameter sets that satisfy progressively more of the process constraints (white indicates none of the parameter constraints being satisfied). The initial sample of 50 000 parameter sets, see region in Stepwise increase of C, thus requiring an increasing number of process constraints to be satisfied at each stage, and (c) the final stage of the search when C is set to the maximum number of design constraints (here 9), at which point the search continues until 8000 behavioral parameter sets have been located. Fig. 7a, consists of samples drawn uniformly from the entire parameter space, of which less than 10 % satisfy any of the imposed constraints, and only a very few satisfy 1, 2 and 3 constraints (progressively darker shades of grey). Each progressive step (see region in Fig. 7b) then consists of S = 5000 samples generated according to the strategy discussed above, in which the value of C is increased by one. Clearly, at each step, the proportion of newly generated samples satisfying larger numbers of constraints increases rapidly. Note that when C was set equal to 9 (i.e., the maximum number of design constraints), the search was continued until the prespecified 8000 behavioral parameter sets had been found, and then terminated. At this point the fraction of number of newly generated samples (behavioral parameter sets) that satisfy all of the constraints has reached approximately 40 % (region in Fig. 7c).
Of course, in this illustrative case study, the constraints were specifically designed to guarantee that the observed hydrograph corresponding to the "best" performing parameter set lie within a predetermined feasible space. In principle the constraints can be specified without recourse to the information contained in the discharge time series (as discussed earlier). The main purpose of this synthetic case study was to illustrate the capability of the proposed CBS algorithm to efficiently locate behavioral parameter sets that satisfy user-specified a priori parameter and a posteriori process constraints.

Conclusions
One of the most challenging tasks in the development of complex conceptual hydrological models for simulation of catchment responses to inputs is the realistic specification of parameter values. We have presented a constraint-based search strategy that facilitates the incorporation of expert knowledge (i.e., the modeler's perception of catchment behavior and characteristics) into the parameter specification process. Because the CBS algorithm does not require observational data regarding the target system output (e.g., runoff) it can provide a way forward for better prediction in ungauged basins in absence of streamflow (or other system output) data for model calibration. As constraints are much easier for understand, rather than parameters, when discussing system behavior, they can potentially be used as an efficient tools to bridge the gap in the dialogue between modelers and experimentalists. Further, the approach can help to provide behaviorally more conceptually realistic parameter sets when used in conjunction with model calibration. Future study may apply the proposed CBS algorithm to identify behavioral parameter sets for different kind of hydrologic models in different regions, and compare the results with other existing parameter specification methods and algorithms. A Matlab TM code of the CBS algorithm can be obtained through personal communication with the lead author.