Maximum entropy production: can it be used to constrain conceptual hydrological models?

In recent years, optimality principles have been proposed to constrain hydrological models. The principle of maximum entropy production (MEP) is one of the proposed principles and is subject of this study. It states that a steady state system is organized in such a way that entropy pro- duction is maximized. Although successful applications have been reported in literature, generally little guidance has been given on how to apply the principle. The aim of this paper is to use the maximum power principle - which is closely related to MEP - to constrain parameters of a simple con- ceptual (bucket) model. Although, we had to conclude that conceptual bucket models could not be constrained with re- spect to maximum power, this study sheds more light on how to use and how not to use the principle. Several of these is- sues have been correctly applied in other studies, but have not been explained or discussed as such. While other studies were based on resistance formulations, where the quantity to be optimized is a linear function of the resistance to be identified, our study shows that the approach also works for formulations that are only linear in the log- transformed space. Moreover, we showed that parameters de- scribing process thresholds or influencing boundary condi- tions cannot be constrained. We furthermore conclude that, in order to apply the principle correctly, the model should be (1) physically based; i.e. fluxes should be defined as a gra- dient divided by a resistance, (2) the optimized flux should have a feedback on the gradient; i.e. the influence of bound- ary conditions on gradients should be minimal, (3) the tem- poral scale of the model should be chosen in such a way that the parameter that is optimized is constant over the modelling period, (4) only when the correct feedbacks are implemented the fluxes can be correctly optimized and (5) there should be a trade-off between two or more fluxes. Although our ap- plication of the maximum power principle did not work, and although the principle is a hypothesis that should still be thor- oughly tested, we believe that the principle still has potential in advancing hydrological science.


Introduction
Traditionally, hydrology is concerned with predicting extremes or water balances for water resources management. But the simplified model structures induce errors, while model parameters are often calibrated on an observed integrated catchment response (often discharge). In recent years, optimality based principles have been suggested to be able to better estimate model parameters and thus model behaviour (e.g. McDonnell et al., 2007;Kleidon and Schymanski, 2008;Clark et al., 2011;Schaefli et al., 2011;Thompson et al., 2011).
The basic idea of optimality principles is that nature organizes itself in such a way that its functioning is optimal under given external forcing during steady state conditions. This can be simulated by taking into account competition between different plant species or trade-offs between fluxes that are driven by different gradients. This competition should then be translated into an objective function. For example, Rodriguez-Iturbe et al. (1999), Porporato et al. (2001) and Caylor et al. (2009) minimized water stress as the objective function for vegetation in (semi-)arid areas. Maximum transpiration and minimal water and oxygen stress have been used by Brolsma and Bierkens (2007), who simulated the competition between two vegetation species, while Published by Copernicus Publications on behalf of the European Geosciences Union. 3142 M. C. Westhoff and E. Zehe: MEP to constrain conceptual models Schymanski et al. (2009b) optimized net carbon profit under given environmental conditions. Another proposed organizing principle is the maximum entropy production principle (MEP) (McDonnell et al., 2007;Kleidon and Schymanski, 2008;Kleidon, 2009Kleidon, , 2010aZehe and Sivapalan, 2009;Schaefli et al., 2011). However, in the hydrological community it remained -so far -mainly on the visionary level.
The above mentioned studies are just a few examples of optimality based models (mainly) in ecohydrology. For a more extensive review, the reader is referred to e.g. Schymanski et al. (2009a).
Several different optimality principles have been used so far, and the question is which one to use. Paik and Kumar (2010) stated that many optimality principles are useful, but that MEP has at least a physical background. Although the physical background of MEP is still under debate, (e.g. Dewar, 2009, suggests that MEP is rather a statistical principle, where the state of MEP is just the most probable one), it seems a useful principle, and many of the proposed objective functions are related to MEP (Dewar, 2010). It is also the objective function used in this study (in which we use it as a physical principle).
The principle of MEP relies on the fact that a gradient drives a flux, while the same flux depletes the gradient. This has been clearly shown by Paltridge (1979) and Lorenz et al. (2001) who used the principle to explain the observed atmospheric movement from the Equator to the poles for planetary systems. Dewar (2010) showed how MEP could be used for plant optimization at different scales, but a comparison with observations was missing. A more hydrological application was formulated by Kleidon and Schymanski (2008), who used the principle to describe the partitioning between runoff and evaporation, but they also did not test their theory with observations. Schymanski et al. (2010) used a simple 2-box model and MEP to predict pattern formation of vegetation in semiarid regions, which gave similar results as the large-scale distributed model of Klausmeier (1999).
Comparison with observations has been done by Zehe et al. (2010), who used the principle of maximum energy dissipation (which is equivalent to MEP) to explain the observed larger density of worm burrows at the foot of the hillslope compared to the hill top. However, they were not able to explain the total number of observed worm burrows.
Another example was illustrated by Porada et al. (2011), who used MEP to constrain parameters for a physically based model based on multiple 1-D columns to simulate the water balance of the largest 35 catchments on Earth. Two parameters were constrained with respect to MEP; namely, resistance for root water uptake by maximizing entropy produced by rootwater uptake and hydraulic conductivity by maximizing entropy produced by baseflow. Subsequently they ranked the different hydrological processes with respect to the produced entropy. They found that transpiration has the strongest contribution to annual average entropy production. However, limited guidance was given on why only two parameters were constrained with the MEP principle, or why entropy production of only certain fluxes were optimized (and not the entropy produced by the entire system). Also, comparison between observed and simulated fluxes was only in terms of the long-term average water balance. Their work was a key motivation for our work, in which we use a conceptual bucket model that is easy to constrain to reproduce and predict rainfall runoff response, but appeared not so easy to constrain with respect to MEP.
Thus, so far, MEP has been suggested as a promising method in hydrology, but optimized model structures have not been extensively tested against observations: the added value of MEP for parameter estimation has not been explored yet. The aim of this paper is to use the MEP principle (although we used the equivalent maximum power principle) to constrain parameters of a bucket model. Guidance is given on which parameters can possibly be constrained with MEP and if these optimal parameter values indeed lead to realistic model performance. We start with constraining each parameter individually, and progressively try to constrain two or more parameters at a time.
This study reveals two main problems of using MEP to constrain conceptual models. The first is that most conceptual models are based on a set of connected (non-)linear reservoirs where outgoing fluxes only depend on the state of its reservoir of origin, which reduces the possibility for process interactions and trade-offs between fluxes. The second problem is that several fluxes in conceptual hydrological models are not driven by potential gradients, which is of prime importance for calculating entropy production or power produced by these fluxes.

Short introduction to MEP
In this section we give a short introduction to MEP, with an emphasis to the application of this article. For a more detailed and fundamental description of MEP, the reader is referred to e.g. Kondepudi (2008), Kleidon (2010b) or Kleidon et al. (2013).
The second law of thermodynamics states that entropy cannot be consumed but only produced. This is the case during irreversible processes that cannot be reversed in time. For isolated systems this implies that the system evolves to a state of maximum entropy, which, for gases and fluids, implies perfect mixing. Perfect mixing means depletion of all gradients: the system is in a state of maximum disorder.
Open systems may, however, exchange energy and mass with their environment. Organized structures may form and persist when incoming fluxes provide the necessary free energy and entropy is exported to the environment. Persistent exchange of energy flows through the system requires a persistent macroscale driving gradient that spans across the entire system. The Bernard cell (Prigogine, 1989) but also planet Earth (Kleidon, 2010b)  fluxes along such macroscale gradients would deplete their driving gradient if there would be no external forcing. In the Bernard cell this is the heating at the bottom and cooling at the top of the cell. In planet Earth this is planetary energy exchange (the poles receive less radiation input than the Equator). The MEP principle now states that such an open system in steady state is structured in such a way that it maximizes entropy production. In case of no external forcing this implies that the system reaches thermodynamic equilibrium as fast as possible.
In the remainder of this paper we use the principle of maximum power, which is closely related to MEP. Following Kleidon et al. (2013), we start with expressing the internal energy dU [J] as where T is absolute temperature , M k is the mass of matter within the system [kg] and k its corresponding chemical potential (or energy level) [m 2 s −2 ]. If we now consider a steady state, so that dU = 0, dT = 0 (we neglected effects from sensible heating/cooling of the soil, as temperature gradients do not drive water fluxes in the subsurface), dp = 0, dV = 0, d k = 0 and incoming fluxes equal outgoing fluxes, then Eq. (1) can be simplified to Rearranging Eq.
(2) and dividing it by a change in time dt, yields Since k dM k /dt is a flux times a gradient (which is power), power is proportional to entropy production and maximizing one is the same as maximizing the other.

Model set-up
For the MEP principle to work the system should be in steady state. We therefore focused on the average yearly water balance, which we assume to reflect steady state partitioning of rainfall into hydrological fluxes (Kleidon et al., 2013). We used a simple bucket model to which we applied the maximum power principle. The model used is in accordance with several bucket models such as e.g. the HYMOD conceptual watershed model (Moore, 1985), the HBV model (Lindström et al., 1997), the GR4J model (Perrin et al., 2003) or the SU-PERFLEX model environment  which all have been successful in reproducing discharge in many hydrological settings. Testing the MEP hypothesis using a bucket model implied two main difficulties. The first was to design a soil moisture accounting scheme that is simple but nevertheless based on interpretable parameters (see Sect. 3.1). The second was to define proper descriptions of driving gradients since most bucket models do not account for the real gradients that drive these fluxes in an explicit manner. The used expressions for the driving gradients were thus not always derived from first principles, but they were constructed on plausible reasoning. An intensive discussion on the gradients is added in Sect. 6.2.1.
All fluxes in the model are described in m s −1 per unit area giving power the unit of W m −2 . The model was run on a hourly time step and solved with an implicit numerical scheme.

Unsaturated zone
The model is similar to the HBV model with the difference that the fast runoff reservoir is now fed directly from the unsaturated zone. (Fig. 1). For simplicity we assumed no snowfall and interception. The water balance of the bucket is given by   Table 1).
where E pot is the potential evaporation [m s −1 ], determined with the Penman formula (Monteith, 1981). Note that the albedo [−] is treated as a calibration parameter: a large albedo means more reflection of solar radiation leading to less potential evaporation. S max is the maximum storage of the bucket [m], and F C is the field capacity [−]. Although this is not a soil physical field capacity, this parameter can be interpreted in a soil physical sense. Soils with high clay content store a large amount of water against gravity. Hence, F C should be close to 0.8. Sandy soils have a low field capacity, thus F C should be around 0. Overland flow and its driving gradient are given by and soil porosity [−], respectively. The driving gradient is thus derived as the potential energy multiplied by the scaling factor (S M /S max ) β . Reasoning for this scaling factor is that we interpreted Q d as fast runoff from the hillslope, where saturation starts at the bottom of the slope. The value of β determines which fraction of the hillslope is saturated for a certain value of S M . The driving gradient is thus expressed as a sort of weighted average of these potential energies. Groundwater percolation and its driving gradient are given by where ∇ G is the gradient driving G [N m −2 ], and K s is the hydraulic conductivity of the soil [m s −1 ]. Note that the driving gradient is expressed as a potential gradient above field capacity, and that, though being simple, this conceptualization of percolation can be parameterized in a manner that is consistent with soil physics. For instance a sandy soil has a high K s , around values of 10 −4 m s −1 , in combination with a small value of F C which leads to a physically consistent behaviour of percolation and direct runoff production. The opposite is true for finer grained soils.

Fast runoff reservoir
The water balance for the fast runoff reservoir (middle bucket in Fig. 1) is given by The fluxes and driving gradients of Q 0 and Q 1 [m s −1 ] are given by where S 1 is the storage height [m], L is the storage threshold level [m] above which Q 0 becomes active and k 0 and k 1 are the reservoir constants [s −1 ]. The driving gradients are both given as a potential energy.

Slow runoff reservoir
The water balance for the slow runoff reservoir (lowest bucket in Fig. 1) is given by where Q 2 [m s −1 ] is given by where S 2 is the storage height [m] and k 2 is the reservoir constant [s −1 ].

Data and calibration procedure
We applied the model to two different watersheds: the first is watershed 02 of the HJ Andrews experimental forest (OR, USA, latitude 44 • 12 43.254 N, longitude 122 • 14 41.554 W). For a detailed description of the watershed the reader is referred to e.g. Rothacher (1965); Rothacher et al. (1967) or Tague and Band (2001). Here we give only a short description. The watershed has a surface area of 60 ha with elevations ranging from 545 to 1070 m, and a mean slope of 30 %. The catchment is completely forested, primarily by mature Douglas fir. The underlying bedrock is volcanic material from the Oligocene to lower Miocene which is overlain by a thick layer of weathered, unconsolidated material of generally high porosity (∼ 0.6) and hydraulic conductivity (∼ 90 m day −1 ) (Tague and Band, 2001).
Winters are generally cool and wet, while summers are warm and dry. Precipitation over the study period ranges from 1300 to 3000 mm yr −1 with a mean of 2222 mm yr −1 . Mean runoff is 1380 mm yr −1 . In general, there is no seasonal snow pack in winter, but daily to weekly snow accumulation does occur.
We used meteorological data from the PRIMET (primary meteorological) station located 900 m from the site (latitude 44 • 12 42.8148 N, longitude 122 • 15 21.3876 W; Daly and McKee, 2011). The data set contains air temperature, solar radiation, relative humidity, precipitation and wind speed on a 1 h time step. Discharge was taken from Johnson and Rothacher (2009), also at an hourly time step. For our model we used the time series from 14 September 1994 to 20 August 2006, of which the first 11 months were used as warming up period.
The second watershed is the upper half of the Weiherbach catchment in southwestern Germany. The watershed drains an area of 3.5 km 2 with gentle slopes characterized by a typical loess soil catena formed by erosion with soil depths ranging between 15 and 30 m. At the foot of the hillslopes the soil consists of moist but drained Colluvisols while in the upper hillslope sections drier Calcaric Regosols are found. Porosity of the soils is in the order of 0.3-0.4 with hydraulic conductivities ranging between 0.26 and 3.4 m day −1 (Zehe et al., 2006). The total catchment area used for cultivation of agricultural crops or pasture is 95 %, 4 % is forested and 1 % is paved area. Runoff is dominated by a relative constant baseflow combined with short peaks of Hortonian overland flow.
The climate is semi humid with an average annual precipitation of 800 mm yr −1 , average annual runoff of 150 mm yr −1 , and annual potential evapotranspiration of 775 mm yr −1 . For this study, we used discharge data from gauge Menzingen (latitude 49 • 08 09.63 N, longitude 8 • 44 44.51 E) and air temperature, solar radiation, relative humidity, precipitation and wind speed were measured within the catchment. All measuring intervals were ≤ 1 h. For this study we used the time series from 1 November 1990 to 1 November 1996, of which the first 12 months were used as warming up period.
In a first step, we varied only one parameter (while all other parameters were arbitrarily chosen) and determined the power of each flux. Because no feedbacks were implemented between the different reservoirs, we could analyse each reservoir separately. For the unsaturated zone, we compared the simulated cumulative yearly fluxes with the observed one, If the landscape is shaped according to maximum power, and the model appropriately represents for the energy conversions associated with partitioning of rainfall into evaporation, and fast and slow runoff components, a maximum in power should correspond with Q err ≈ 0.
To compare simulations of the fast and slow runoff reservoirs with observations, we nevertheless combined the three fluxes Q 0 , Q 1 and Q 2 to be able to compare it with the complete runoff time series. These runoff reservoirs are fed by fluxes from the unsaturated zone, which we parameterized in such a way that the long-term water balance is closed (Q err ≈ 0).
In a second step we only took those parameters that could be constrained in the first step and varied two of these parameters at the same time, provided that the two parameters parameterize the same reservoir. Again, the power of each flux was determined and the model output was compared with observations in the same way as in the first step.

Constraining one free parameter
In our first step, we tried to constrain one parameter at a time while keeping all other parameters constant.

Unsaturated zone
For the unsaturated zone, results for two different parameter sets are shown (Figs. 2, 3, solid and dashed lines; the parameter values for each run are listed in Table 1). For the HJ Andrews (Fig. 2), the shown simulations are done with parameter sets that resulted in a closed water balance that coincides with maximum power along a range of β (panel a): for the solid lines this is maximum power in G (P G ), and for the dashed line this corresponds to maximum power in Q d (P Q d ) and in total power (P tot ). Note, however, that in many other combinations of the fixed parameters, a closed water balance did not coincide with maximum power in any of the fluxes and since parameters for these kind of models can often not be directly measured in the field, the values of each fixed parameter are in principle arbitrary.
Most power was produced by either G or Q d . Power produced by E a (P E a ) was relatively small and insensitive over the parameter spaces. The reason for this is that none of the parameters influenced the gradient driving E a and no feedback between the flux (E a ) and its driving gradient (∇ E a ) was implemented in the used model set-up. In contrary to the model of Porada et al. (2011), conceptual models do not account for the capillary rise and the associated power in this flux.
For the parameters F C , S max and albedo, power is maximum at the boundaries of the parameter space ( Fig. 2b, d, e), which means they could not be constrained with respect to power. Nevertheless, two parameters do result in maxima in power, namely β and K s (Fig. 2a, c).
β results in both maximum power by G (for large β) and maximum power by Q d (for small β): both fluxes that are directly influenced by β. For β larger than 50, the water balance becomes insensitive, which is due to the fact that for such a large β the water level (S M ) should be so close to the maximum water level (S max ) to produce direct runoff that it practically functions as a binary threshold where β becomes insensitive. Such a large β represents a catchment with a uniform soil depth.
Along K s , a maximum in power by G appeared. But with the chosen fixed parameters, the water balance is rather insensitive over the full range of K s , which is caused by the relatively large value of F C . Although we tested many different combinations of fixed parameters for which a maximum in power along K s appeared (results not shown here), we could not find a parameter set at which a closed water balance coincided with a maximum in power at the observed value (by Tague and Band, 2001) For the Weiherbach catchment we could only find a β where a closed water balance corresponded with a maximum in P G (Fig. 3a), while a maximum in P Q d always corresponded with a large error in the water balance: in the shown case this error is −362 mm yr −1 , which is almost half of the yearly rainfall. Nevertheless we used this optimum value for β as fixed parameter when testing the other parameters.
Similar to the results of the HJ Andrews, maxima in power could only be found for β (Fig. 3a) and for K s (Fig. 3c). However, for the latter this is not always the case (see solid line in  Table 1).
in such a way that the water level will never be larger than the threshold F C S max that should be exceeded to get groundwater percolation G. In these cases, K s is insensitive. What is also striking is that when G = 0, the error in the water balance is rather insensitive for S max (Fig. 3d). This means that for larger values of S max , Q d remains almost constant, while power only increases due to the larger gradient driving Q d . Contrary to the model results of HJ Andrews is that in many cases most power is produced by evaporation (P E a ), while power produced by G is always small and power produced by Q d is only in some cases large. But the latter cases correspond with a large error in the water balance. The main reason for these results is that the runoff coefficient is rather low (< 20 %) meaning that the evaporative flux is relatively large, leading to large power by E a .

Fast runoff reservoir
We followed the same procedure for the fast runoff reservoir. Since the fast runoff reservoir is fed by Q d from the unsaturated zone, the unsaturated zone had to be parameterized as well. Here, we focus only on results from HJ Andrews because for the Weiherbach catchment the maximum in P Q d only existed for large water balance errors, making results for the fast runoff reservoir by definition unreliable, while the maximum in P G only appeared for Q d = 0, meaning that  0.71 0.7 1 × 10 −7 1200 0.05 100 0.01 HJ-Andrews * The slow runoff reservoir was added to be able to compare simulated runoff with the observed one. * * All parameters were fixed, except the one varied along the x axis. no water flows into the fast runoff reservoir. The latter case is, however, consistent with the way the catchment functions.
For the HJ Andrews, the chosen parameter sets for the unsaturated zone came from the solid and dashed line in Fig. 2a, where a closed water balance coincides with maximum power for G (solid line in Fig. 2a and panels a 1 -a 3 in Fig. 4) or Q d (dashed line in Fig. 2a and panels b 1 -b 3 in Fig. 4). For both cases we showed again results of two sets of fixed parameters (see Table 1).
For parameter L, which represents the threshold above which Q 0 becomes active, no maximum in power occurs (Fig. 4a 1 , b 1 ). For large values of L, power only tops off and becomes flat. This is due to the fact that, given certain forcing and reservoir constants (k 0 and k 1 ), the water level in the reservoir will have a certain maximum. If L is larger than this water level, it has no influence on the fluxes anymore and thus no influence on the power either.
For k 0 and k 1 , a maximum in power is observed in the corresponding fluxes: thus a maximum in P Q 0 for k 0 (Fig. 4a 2 , b 2 ) and in P Q 1 for k 1 (Fig. 4a 3 , b 3 ). Note that a maximum in power in panels a 1 -a 3 are an order of magnitude smaller than in panel b 1 -b 3 : this is due to the fact that the input flux Q 0 is much smaller for panels a 1 -a 3 .
Besides the maxima in power by the Q 0 and Q 1 , a maximum in total power (P tot ) only appears (but not always) for k 1 (Fig. 4a 3 , b 3 ) and the location of this maximum is always close to the location of maximum power by Q 1 . This is not strange, since P Q 1 is for most cases larger than P Q 0 . Only when the water level becomes sufficiently higher than threshold L (which only occurs for small values of k 1 ), P Q 0 may become larger than P Q 1 .
The link to observations is made by looking at the Nash-Sutcliffe model efficiency (NSE) of simulated vs. observed discharge. To do this properly, total runoff should be compared which means that runoff by the slow runoff reservoir should also be added. Therefore, we gave k 2 the arbitrary value of 0.02 s −1 . Note that this value has little influence for the simulations in panels b 1 -b 3 , since G (the flux flowing to the slow runoff reservoir) is only about 2 % of Q d : thus 98 % of the total runoff is produced by the fast runoff reservoir.
For most cases, the maximum in power is far away from the maximum in NSE. Only for one (shown) case the two correspond ( Fig. 4b 3 , dashed line). Although the variation in NSE is rather small over the full parameter space, this parameter set may be an optimal one. However, this may also be a coincidence (see Sect. 6.1).  Table 1).

Slow runoff reservoir
The same procedure was also followed for the slow runoff reservoir, albeit that this reservoir has only one parameter and one flux leaving the reservoir. This reservoir is fed by G coming from the unsaturated zone and the presented results were obtained with the same two unsaturated zone scenarios as for the fast runoff reservoir (see also Table 1). Although we found a maximum in P G that corresponded to a closed water balance for the Weiherbach catchment, we only show results from the HJ Andrews because results from both catchments are similar for this reservoir. The obtained power by Q 2 is at its maximum at the lower end of k 2 (Fig. 5). When thinking this over, this is not surprising, since the maximum power principle assumes a steady state system, meaning that on average G = Q 2 . This means that the flux is fixed and that power is only influenced by the driving gradient. Thus the smaller k 2 is, the larger the gradient should be to obtain a large enough flux. Power is thus maximum at the asymptote k 2 → 0. This asymptote vanishes  Table 1 for the used parameters). The fast runoff reservoir was also parameterized to be able to determine the NSE of the total outflow. when there are more fluxes leaving the reservoir, since Q 2 is then not fixed anymore.
When comparing these results with observed discharge, maximum NSE occurs at relatively large values of k 2 . Although the NSE is also influenced by the other reservoirs, the low NSE for small values of k 2 comes from the fact that smaller k 2 values lead to more constant outflow, meaning that at maximum power (lim k 2 →0 ) the outgoing flux Q 2 is completely constant.

Constraining two free parameters
For both the unsaturated zone and the fast runoff reservoir, there were two parameters that showed a maximum in power. For the unsaturated zone, we found several parameter sets where a closed water balance coincided with a maximum in power along one variable, and also for the fast runoff reservoir we found a parameter set where maximum NSE coincided with maximum power. However, this approach means that only one parameter per submodule can be constrained with the maximum power principle while the other parameters should be chosen in another manner. Porada et al. (2011) proposed a method to be able to constrain two parameters with respect to entropy production. They created surface plots with on the x and y axis the two parameters and on the z axis the entropy produced by the two fluxes that were controlled by the two parameters. They concluded that the optimal parameters were the location where the "ridges" of the surface plots of the two fluxes intersected.
We did the same here, and combined the two parameters that showed an optimum; i.e. β and K s for the unsaturated zone (Fig. 6) and k 0 and k 1 for the fast runoff reservoir (Fig. 7). If the approach of Porada et al. (2011) works here, the intersect of the two ridges should coincide with a closed water balance for the unsaturated zone and maximum NSE for the fast runoff reservoir.
Also here, we focus only on results from the HJ Andrews, because no reliable optimum for P Q d could be found for the Weiherbach catchment.

Unsaturated zone
For the unsaturated zone, this resulted in both local and global maxima in total power (P tot ) and power produced by Q d (P Q d ), a global maximum in power produced by G (P G ) and an almost flat surface in power produced by evaporation (P E a ). P Q d is characterized by a ridge at β ≈ 0.5 where K s is insensitive, while at the largest values of β power is maximum for the smallest values of K s . This edge of the parameter space coincides with the global maximum in power.
For the first parameter set ( Fig. 6a; see Table 1 for the used parameters), the water balance is closed for β larger than 50, which corresponds with the global maximum in P tot , P Q d and P G . However, the global maxima of P Q d and P G are located at opposite ends of the parameter range of K s (panels a 2 and a 3 ). Also, the "ridges", as used by Porada et al. (2011), do not intersect: for P Q d , the ridge only exists for K s > 1×10 −8 m s −1 (for smaller values a ridge exists along the full range of β at K s = 1 × 10 −10 m s −1 , which is the edge of the parameter space), while for β < 1, no maximum in P G exists anymore (P G = 0). The point where the ridges almost overlap is around K s = 1 × 10 −8 m s −1 and β = 0.5, but here the water balance is far from being closed (Q err = −283 mm yr −1 ).
For different values of S max , F C and albedo, the surfaces of the plots are similar to those created with the first parameter set ( Fig. 6b; see Table 1 for the used parameters). The main difference is that the water balance is now closed around the ridge in P Q d . Thus the point where the ridges almost overlap coincides with a closed water balance. Following the approach of Porada et al. (2011), this should be the optimal parameter set.

Fast runoff reservoir
As input to the fast runoff reservoir we used the same two time series of Q d as in Sect. 5.1, of which the second is simulated with parameters close to the optimum defined after the approach of Porada et al. (2011). For each of the two time series of Q d , results are shown for two different values of L: Fig. 7a Table 1 for the used parameters). The colour coding shows the average error in the yearly water balance (Q err ).
Although the "ridges" in panels 2 and 3 are sometimes difficult to see they do exist. The intersect of the ridges in P Q 0 and P Q 1 are always at values of k 0 and k 1 that are smaller than 1 × 10 −2 s −1 , where the NSE is relatively low. This also applies for the case where Q d is (close to) optimal (Fig. 7c,  d). We can thus state that even if the unsaturated zone can be correctly constrained by using the maximum power principle, it does not seem to work for the fast runoff reservoir.

Discussion
The question of whether MEP and maximum power are indeed fundamental principles -such as the first or second law of thermodynamics -is still under debate. Nevertheless, several authors already showed their potential (e.g. Kleidon and Schymanski, 2008;Schymanski et al., 2010;Zehe et al., 2010;Porada et al., 2011;Kleidon et al., 2013), but generally little guidance was given on how to use the principle. For example, Kleidon and Schymanski (2008) set up a simple model to illustrate how MEP could be used to estimate the resistance for root water uptake. Although they did  Table 1 for the used parameters). The colour coding shows the NSE of the simulated runoff.
apply the principle in a correct way, they mainly stressed that the principle works because there is a trade-off between a flux and its gradient, but they did not explain explicitly that there should also be a trade-off between two or more fluxes -which is clearly demonstrated with our model results of the slow runoff reservoir, where power is maximum at the asymptote of k 2 → 0 s −1 . Also the study of Porada et al. (2011), which provides a major motivation for the present study, applied and tested MEP in a sophisticated manner to large catchments; however, some of their findings were reported in a confusing or imprecise manner: they reported produced entropy by all fluxes, while they only constrained two parameters of two different fluxes using MEP. Thus only the produced entropy of these two fluxes were relevant for their aim to constrain parameters of their hydrological model. This left room for this study in which we tried to use the maximum power principle to constrain a conceptual bucket model. During this exercise we were confronted with the questions of what exactly we should optimize, and if the model used was appropriate to apply the principle.

What should and what can we optimize?
All of the above mentioned studies constrained a resistance (or conductance) of a flux by maximizing entropy production (or power) of that flux. But they neither explained clearly why they only focused on resistances, nor did they explain why they maximized entropy production of single fluxes (as opposed to MEP of the whole system).
The answer of why only resistances were constrained is at first glance rather simple: because resistances can mathematically be optimized with MEP (the word mathematically is essential here, because mathematical optimum is not necessarily the same as the optimal structure or state of the real world system). In this study we also showed that resistances (or in our case conductances) could be constrained with respect to power, but we also showed in this study that β could (at least mathematically) be constrained as well. Where resistances influence fluxes in a linear manner, the influence of β on the flux is only linear in the log-transformed space. We consider this as an important finding, since this may broaden the applicability of MEP; for example, the soil moisture retention curve is also non-linear, and could thus be a candidate to be (partly) parameterized using the MEP principle.
Parameters describing a process threshold (e.g. S max , F C and L) or parameters that influence boundary conditions (e.g. albedo) could, on the other hand, not be constrained. The reason that process thresholds cannot be constrained is that a threshold only influences the driving gradient: a lower threshold leads to a larger driving gradient of the flux that is controlled by the threshold (see e.g. Fig. 4a 1 , b 1 , light cyan lines). As a consequence, these parameters that cannot be constrained with respect to power should, if they cannot be directly measured in the field, be estimated in another way. In this research we estimated them by constraining them in such a way that the water balance was closed.
For the unsaturated zone this reduced the amount of behavioural parameter sets significantly, but did not lead to a unique solution. However, when we added the runoff reservoirs, we did find a unique solution where a maximum in power by Q 1 corresponded with a maximum in NSE (Fig. 4b 3 ). Although this looks promising at first sight, we believe it is rather coincidence than an inevitable result of our approach. Because why does this optimum only appear when optimizing β and k 1 ? And how should one know on the beforehand that one should optimize these two parameters, as opposed to optimizing K s and k 0 ?
The reason why only entropy production of single fluxes are maximized can be answered in a similar way: because only then a maximum exists. But this question goes a little deeper: we believe that if natures strives to a state of maximum power, nature strives to a state of maximum in total power. Now let us consider the approach of Porada et al. (2011). They constrained two parameters (resistance of two different fluxes) at the same time and considered the intersect of the "ridges" of produced entropy by the two different fluxes as the optimal one. In their model, produced entropy by one flux was two orders of magnitude higher than produced entropy by the other flux. Thus the question arises, why does nature not maximize entropy production by the first flux? This may lead to lower entropy production by the second flux, but since this was two orders of magnitude smaller, total entropy production will still be larger. However, if this reasoning would be applied, entropy production by the first flux would be maximum at the border of the surface plot, meaning that no maximum exists. The results of our own model show the same: the (global) maximum in total power is often at the edge of the parameter space and can thus not be optimized. This is especially clear in Fig. 7, panels 1, where maximum total power is produced by the smallest values of k 0 and k 1 . This is due to the fact that there is no feedback implemented in the model that limits the power (or entropy production) when there are two free parameters. Thus, we cannot mathematically maximize total power and we come back to the simple first answer: it is a mathematical issue. This is in our opinion not strange, because if we want to maximize total power, we should take all processes into account; at all spatial and temporal scales. This includes erosion, continental uplift, power created by biota, dissipation of energy into heat, etc. But our models are simplifications and many processes are considered constant at the temporal or spatial scale of our model (e.g. continental uplift) or function as boundary conditions (e.g. potential evaporation).
The question is thus: does nature acts in a similar way, where a fast process considers a slow process as constant and where processes that are not (or to a minimum) influenced by other processes see the other processes as boundary conditions? If this is indeed the case, maximizing power of single fluxes is indeed correct and the approach of Porada et al. (2011) is a clever way to constrain two parameters simultaneously. However, although Porada et al. (2011) reported seemingly good results with their approach, they still reported differences in mean long-term observed and simulated runoff of up to 300 % for temperate climates (estimated from their Fig. 5b) which is for hydrological models too large (although this difference may also be caused by other aspects).
Thus, to conclude, maximizing total power (or entropy production) is practically impossible since all processes and feedbacks should be taken into account, while maximizing power of single fluxes show at least mathematically an optimum. But it still has to be foreseen if nature indeed behaves like this.

How to apply the candidate principle?
In this study we applied the maximum power principle to a conceptual bucket model and we found maxima in power for several different parameter sets, of which one parameter set also coincided with a maximum in NSE. Also when constraining two parameters simultaneously (with the approach of Porada et al., 2011) we found an optimum for the unsaturated zone that maximized power and had a closed water balance. However, applying the same method for the fast runoff reservoir, the parameters that were optimal with respect to power were far from optimal with respect to reproducing observed discharge (Fig. 7). This means that the mathematical optima we found, did not reflect the catchment functioning.
Therefore we conclude that either the maximum power hypothesis is not valid or the model is not suitable for testing the hypothesis. If we assume for now that the principle is correct, our model is thus not suitable. This means that conversions of energies associated with rainfall runoff processes are not represented in a feasible manner: conceptual models have not been developed for this issue, but with the aim to reproduce runoff.
Nevertheless, by applying the principle to our model we gained more insight in how the principle should be applied and why other authors applied the principle in the way they did. We will discuss the following issues: (1) thermodynamic consistency; i.e. fluxes should be described as gradients divided by resistances, and (2) system boundaries and feedbacks.

Flux = gradient divided by resistance
To be consistent with thermodynamic principles, each flux in the model should be defined as a gradient divided by a resistance (or multiplied by a conductance), and power is then the flux times the gradient; thus power is the gradient squared divided by the resistance. In this way the flux has a feedback to its gradient (as stressed by e.g. Kleidon and Schymanski, 2008); i.e. a larger flux depletes the gradient faster, and, consequently, the flux will be smaller. The resistance is then treated as the unknown that can be constrained with the MEP principle.
Our model concept is often not consistent with thermodynamics. An obvious one is Q d . Of course, S M is part of the driving gradient, but as long as it does not rain, Q d and S M have no connection at all, let alone a feedback. The reason that this flux could be optimized is because the equation was non-linear.
Although K s could also be constrained, the formulation of G was not consistent either: where K s is the conductance, the force that drives the flux ((S M /S max ) β ) is not the same as the defined driving gradient (∇ G = ρg/η max (0, S M − S max F C )). Power is thus not described as a conductance times a gradient squared.
The fluxes leaving the fast runoff reservoir were described in a consistent manner, where the gradient is only scaled by ρg/η to get the correct units for power. However, the formulation of the flux and gradient is a description of a conceptual linear reservoir. The real world gradient is different as well as the depletion of the gradient: the runoff from the fast runoff reservoir can be interpreted as overland flow or rapid subsurface flow. For these flows, the real-world driving gradient is thus the slope of the surface or of the subsurface lateral flow network (influenced by e.g. bedrock topography). The depletion of this real-world gradient is thus not so much the depletion of the water level (as in the linear reservoir), but mainly the depletion of the slope by erosion. This is a much slower process and can often be considered constant over the timescale of the modelling period. The optimum parameters for our linear reservoir model will thus be totally different than the real-world resistances. We may thus conclude that the maximum power principle (or MEP) cannot be applied to conceptual models, due to their lack of realism.
Taking a better look at the study of Porada et al. (2011), we learned that although they calculate entropy production of surface runoff and river discharge, they did not optimize these fluxes. The fluxes they did optimize were rootwater uptake and baseflow which was actually the partitioning of soil water into transpiration and runoff, while they did not partition total runoff into different runoff components (probably because of the above mentioned reason that depletion of the driving gradient is a too slow process).
Another flux that we aimed to parameterize with the maximum power principle is evaporation E a . However, to do this properly, this flux should also be described as a gradient divided by a resistance, which was not the case in our model set-up. The main parameter that influenced evaporation was the albedo, which only influences the potential evaporation. As shown in Figs. 2e and 3e this parameter could not be constrained and potential evaporation should thus be treated as a boundary condition.

System boundaries and feedbacks
Since we cannot model all processes and feedbacks at all scales, the boundaries of the model are important. Boundary conditions should be chosen carefully as well as the feedbacks in the model. Both depend on the spatial and temporal scale of the model and on the notion that the system should be at (or at least close to) steady state. The maximum time span of the model is influenced by the process described by the optimized parameter. This is because the optimized parameter should be constant for the whole modelling period.
In this study we took the average yearly water balance as timescale. We thus implicitly assumed that the soil structure and vegetation pattern is constant over the 11 yr of simulation and that it has enough memory not to be significantly influenced during rain storms or dry spells. Unless information is available about periodicity of certain parameters, this is the best one can do and we believe that the yearly timescale is defensible as appropriate timescale for this study. However, other systems require other timescales; for example, for arable land that is ploughed seasonally, a better timescale may be a single rainstorm (e.g. Zehe et al., 2010) and to estimate long-term erosion rates much longer timescales are needed (e.g. Kleidon et al., 2013). 5. there should be a trade-off between two or more fluxes.
This follows directly from the fact that the principle only holds for systems in steady state, meaning that incoming fluxes equal outgoing fluxes.
Furthermore, we would like to stress that with this study we do not claim that the MEP principle does not work, but we showed that applying the principle is not as straightforward as it may seem. Only if it is applied it in a correct way, while having the correct physics, one will be able to test the MEP principle. If that point is reached, and if the MEP principle still holds, it would mean a big step forward in the hydrological science.