A thermodynamic approach to link self-organization, preferential flow and rainfall–runoff behaviour

This study investigates whether a thermodynamically optimal hillslope structure can, if existent, serve as a first guess for uncalibrated predictions of rainfall–runoff. To this end we propose a thermodynamic framework to link rainfall–runoff processes and dynamics of potential energy, kinetic energy and capillary binding energy in catchments and hillslopes. The starting point is that hydraulic equilibrium in soil corresponds to local thermodynamic equilibrium (LTE), characterized by a local maximum entropy/minimum of free energy of soil water. Deviations from LTE occur either due to evaporative losses, which increase absolute values of negative capillary binding energy of soil water and reduce its potential energy, or due to infiltration of rainfall, which increases potential energy of soil water and reduces the strength of capillary binding energy. The amplitude and relaxation time of these deviations depend on climate, vegetation, soil hydraulic functions, topography and density of macropores. Based on this framework we analysed the free energy balance of hillslopes within numerical experiments that perturbed model structures with respect to the surface density of macropores. These model structures have been previously shown to allow successful long-term simulations of the water balances of the Weiherbach and the Malalcahuello catchments, which are located in distinctly different pedological and climatic settings. Our findings offer a new perspective on different functions of preferential flow paths depending on the pedological setting. Free energy dynamics of soil water in the cohesive soils of the Weiherbach is dominated by dynamics of capillary binding energy. Macropores act as dissipative wetting structures by enlarging water flows against steep gradients in soil water potential after long dry spells. This implies accelerated depletion of these gradients and faster relaxation back towards LTE. We found two local optima in macropore density that maximize reduction rates of free energy of soil water during rainfall-driven conditions. These two optima exist because reduction rates of free energy are, in this case, a second-order polynomial of the wetting rate, which implicitly depends on macroporosity. An uncalibrated long-term simulation of the water balance of the Weiherbach catchment based on the first optimum macroporosity performed almost as well as the best fit when macroporosity was calibrated to match rainfall–runoff. In the Malalcahuello catchment we did not find an apparent optimum density of macropores, because free energy dynamics of soil water during rainfall-driven conditions is dominated by increases of potential energy. Macropores act as dissipative drainage structures by enhancing export of potential energy. No optimum macropore density exists in this case because potential energy change rates scale linearly with the wetting rate. We found, however, a distinguished macroporosity that assures steady-state conditions of the potential energy balance of the soil, in the sense that average storage of potential energy is compensated by average potential energy export. This distinguished macroporosity was close to the value that yielded the best fit of rainfall–runoff behaviour during a calibration exercise and allowed a robust estimate of the annual runoff coefficient. Our findings are promising for predictions in ungauged catchments (PUB) as the optimal/distinguished Published by Copernicus Publications on behalf of the European Geosciences Union. 4298 E. Zehe et al.: A thermodynamic approach to link self-organization model structures can serve as a first guess for uncalibrated predictions of rainfall–runoff. They also offer an alternative for classifying catchments according to their similarity of the free energy balance components.

Abstract. This study investigates whether a thermodynamically optimal hillslope structure can, if existent, serve as a first guess for uncalibrated predictions of rainfall-runoff. To this end we propose a thermodynamic framework to link rainfall-runoff processes and dynamics of potential energy, kinetic energy and capillary binding energy in catchments and hillslopes. The starting point is that hydraulic equilibrium in soil corresponds to local thermodynamic equilibrium (LTE), characterized by a local maximum entropy/minimum of free energy of soil water. Deviations from LTE occur either due to evaporative losses, which increase absolute values of negative capillary binding energy of soil water and reduce its potential energy, or due to infiltration of rainfall, which increases potential energy of soil water and reduces the strength of capillary binding energy. The amplitude and relaxation time of these deviations depend on climate, vegetation, soil hydraulic functions, topography and density of macropores. Based on this framework we analysed the free energy balance of hillslopes within numerical experiments that perturbed model structures with respect to the surface density of macropores. These model structures have been previously shown to allow successful long-term simulations of the water balances of the Weiherbach and the Malalcahuello catchments, which are located in distinctly different pedological and climatic settings. Our findings offer a new perspective on different functions of preferential flow paths depending on the pedological setting. Free energy dynamics of soil water in the cohesive soils of the Weiherbach is dominated by dynamics of capillary binding energy. Macropores act as dissipative wetting structures by enlarging water flows against steep gradients in soil water potential after long dry spells. This implies accelerated depletion of these gradients and faster relaxation back towards LTE. We found two local optima in macropore density that maximize reduction rates of free energy of soil water during rainfall-driven conditions. These two optima exist because reduction rates of free energy are, in this case, a second-order polynomial of the wetting rate, which implicitly depends on macroporosity. An uncalibrated long-term simulation of the water balance of the Weiherbach catchment based on the first optimum macroporosity performed almost as well as the best fit when macroporosity was calibrated to match rainfall-runoff. In the Malalcahuello catchment we did not find an apparent optimum density of macropores, because free energy dynamics of soil water during rainfall-driven conditions is dominated by increases of potential energy. Macropores act as dissipative drainage structures by enhancing export of potential energy. No optimum macropore density exists in this case because potential energy change rates scale linearly with the wetting rate. We found, however, a distinguished macroporosity that assures steady-state conditions of the potential energy balance of the soil, in the sense that average storage of potential energy is compensated by average potential energy export. This distinguished macroporosity was close to the value that yielded the best fit of rainfall-runoff behaviour during a calibration exercise and allowed a robust estimate of the annual runoff coefficient. Our findings are promising for predictions in ungauged catchments (PUB) as the optimal/distinguished 1 Introduction

Spatial organization and the search for organizing principles
A closer look at catchments often reveals a spatially organized "architecture" (Dooge, 1986;Schulz et al., 2006;Zehe and Sivapalan, 2009). Spatial organization has different fingerprints and affects different processes. An apparent spatial covariance of soil hydraulic properties (Zimmermann et al., 2008) is for instance a fingerprint of local spatial organization (Kondepudi and Prigogine, 1998), which translates into spatially organized storage of water (Brocca et al., 2007;Western et al., 2004;Zehe et al., 2010b). Hillslope scale spatial organization of soil types is often reflected in a typical soil catena , which can be of key importance for overland flow response and sediment yields. The most striking evidence for spatial organization is, however, the omnipresence of networks of preferential flow paths, which vein soils, unconsolidated rock and even solid bedrock. Independently from their genesis, whether they are bio-pores (worm or rodent burrows, root channels) or fingerprints of past dissipative processes (shrinkage cracks, pipes, erosion rills), they exhibit similar topological characteristics and similar functioning. Preferential flow networks allow for high mass flows at small driving gradients, because the flow resistance is continuously small in the entire network. They thus "organize" export and redistribution of water and matter across scales: locally as vertical macropores (Beven and Germann, 1982), at the hillslope scales as surface rills or subsurface pipe networks (Bull and Kirkby, 1997;Parkner et al., 2007;Weiler and McDonnell, 2007;Wienhöfer et al., 2009) or at the catchment scale as river networks (Kleidon et al., 2013;Howard, 1990). The rationale of this study is to better understand why these organized patterns of soils and preferential flow paths in catchments have formed and persisted up to now. In particular we explore whether the surface and subsurface architecture of hillslopes in distinctly different hydroclimatic settings is in accordance with thermodynamic optimality principles. This implies that these optimal hillslope architectures should allow for acceptable uncalibrated simulations of rainfall-runoff behaviour, as will be explained below.
The fact that persistent spatial organization in catchments exists has inspired many scientists to speculate whether this is a manifestation of self-reinforcing co-development (Dietrich and Perron, 2006) due to an underlying organizing principle. Rodriguez-Iturbe and Rinaldo (2001) explain river networks as "least energy structures" minimizing local energy dissipation. Maximum entropy production (MEP) suggests that steady-state system architectures are organized such that fluxes of mass and energy maximize production of entropy (Kleidon and Schymanski, 2008). The MEP hypothesis has been corroborated within studies that allowed (a) successful predictions of states of planetary atmospheres (Paltridge, 1979;Lorenz et al., 2001), (b) identification of parameters of general circulation models (Kleidon et al., 2006) or (c) identification of hydrological model parameters to estimate the annual water balances of the 35 largest catchments on Earth (Porada et al., 2011). Kleidon et al. (2013) recently explored whether the formation of connected river networks is in accordance with the maximum power principle (Lotka, 1922a, b) and thus whether "free" energy transfer to sediment flows is maximized, which implies fastest depletion of morphological gradients (compare Sect. 2.1 for a definition of free energy).

Objectives, motivation and hypotheses driving this study
The overall objective of this study is to test whether a thermodynamically optimal hillslope structure can, if existent, serve as a first guess for uncalibrated predictions of rainfall-runoff. To this end we propose a novel thermodynamic framework for analysing flows and conversions of kinetic energy, potential energy and capillary binding energy which occur during the rainfall-runoff processes. A thermodynamic treatment of these energy conversions is necessary because they reflect conservation of energy and the second law of thermodynamics, as they imply small amounts of dissipation of free energy into heat and thus irreversibility. This in fact explains why water naturally flows downhill and soil water flows naturally deplete gradients in soil water potential. The kinetic energy input from incoming rainfall is partly dissipated into heat and partly used to break up soil aggregates (Kleidon et al., 2013;Scherer et al., 2012). Depending on soil hydraulic properties and apparent macropores, incoming rainfall is split into surface water and soil water. Potential energy of surface water is partly transformed into kinetic energy of overland flow and exported by means of surface runoff/discharge and partly dissipated by frictional losses at the surface (Fig. 1a). The other fraction of rainfall infiltrates, thereby increasing potential energy of soil water by adding mass at the topsoil. This wetting process does at the same time reduce absolute values of negative capillary binding energy of soil water. These outlined energy conversions are two orders of magnitude smaller than those associated with the surface energy balance, which might explain why they have mostly been ignored in the past. These energy conversions are nevertheless of key importance because they are related to the partitioning of incoming rainfall mass into runoff components and storage dynamics. This splitting as well as the subsequent subsurface dynamics is strongly controlled by preferential flow paths. We thus focus on how macropores affect these E. Zehe et al.: A thermodynamic approach to link self-organization 4299 Fig. 1. Free energy dynamics of a hillslope during rainfall-/massflow-driven conditions; (a) potential energy input from rainfall (Eq. 5), kinetic energy export by means of overland flow (Eq. 6), changes in free energy of soil water arising from changes in soil water content (capillary binding energy and potential energy; compare Eqs. 2 and 7) as well as export of water mass at different potential energy levels (please refer to Eq. 7). Free energy dynamics of a hillslope during radiation-driven conditions; (b) components of the surface energy balance (i.e. net radiation, sensible and latent heat fluxes) as well as soil heat dynamics triggered by the soil heat flux. Note that the soil heat budget is linked to soil water dynamics as the volumetric heat capacity C v (θ ) and thermal conductivity α(θ) depend on the actual soil water content θ . conversions and whether this role is associated with the existence of an optimal macropore density which maximizes average entropy production or, equivalently, maximizes average dissipation of free energy of soil water at a temporal extent where quasi-steady state can be assumed. Zehe et al. (2010a) explored the role of macropores in these energy conversions. They showed that soil hydraulic equilibrium (HE) corresponds to local thermodynamic equilibrium: as soil water potential is uniformly zero in the entire profile, this is a state of maximum entropy or, equivalently, a state of minimum free energy of soil water. Note that free energy/chemical energy of soil water is equal to the sum of potential and capillary binding energy (Kleidon and Schymanski, 2008). Zehe et al. (2010a) provided evidence that vertical macropores accelerate depletion of matric potential gradients during single rainfall events after long dry spells, which implies accelerated reduction of free energy of soil water. In their study the reduction rate of free energy increased almost linearly with increasing surface density of macropores; an optimum density that maximized the reduction rate was however not found at this timescale.
The reported findings and the work of Kleidon et al. (2013) motivate our first hypothesis (H1): macropores (preferential flow paths in general) act as dissipative structures by enhancing reduction (dissipation plus export) of free energy of the hydrological system during the rainfall-runoff process. Preferential flow paths accelerate mass flows against driving gradients in soil water potentials or in surface water levels, which implies a faster depletion of these gradients. Accelerated reduction of free energy means that the system state comes closer (not necessary back) to local thermodynamic equilibrium. This is favourable for mechanic stability of the system because mass flows in the pipe or macropore networks exert less "hydro-mechanical stress" on the system itself and mechanic/hydraulic loads are quickly reduced (Lindenmaier et al., 2005;Ehlers et al., 2011). In the present study we test H1 within numerical experiments that focus on daily and annual timescales and perturb the macropore density of "behavioural" model structures. These have been previously shown to allow successful long-term simulations of the water balances of our test catchments Weiherbach and Malalcahuello. These model structures are thus regarded as best structural and functional representations of the hydrological systems, which are located in distinctly different pedological and climatic settings. The annual timescale is deemed to be appropriate to assess dynamics of free energy during rainfall-runoff processes in a representative manner and long enough that a slow negative feedback could lead to an optimum macropore density that maximizes entropy production or reduction of free energy of soil water (see Sect. 2.3.1 for further explanation).
In case such an optimum exists, will this be of any value for predictions? Here our second hypothesis (H2) comes into play: we postulate that a hydro-geo-ecosystem configuration is closer to a functional optimum than other configurations if it accelerates reduction of free energy of soil and surface waters during rainfall-driven conditions by accelerated redistribution of water against internal gradients in the soil and accelerated export of "excess" water from the system. Hypothesis H2 does not postulate that functionally optimal configurations necessarily exist, but in case they do, H2 implies maximization of export and dissipation of free energy during rainfall-driven conditions. In case an optimum macropore density exists and if the catchment developed indeed in accordance with H2, the optimum should allow for an acceptable uncalibrated simulation of the catchment behaviour to split rainfall into surface and subsurface runoff components. Validation of the uncalibrated prediction is thus a test of H2. In the remaining study we explain the thermodynamics of catchments with emphasis on the role of macropores during the rainfall-runoff processes, and develop the framework to link mass flows and energy conversions (Sect. 2). Section 3 introduces the study areas and the model. After presenting the results of the numerical experiments in Sect. 4, we close with discussion and conclusions in Sect. 5.

Classical thermodynamics in a nutshell
Thermodynamics is a fundamental theory of physics dealing with the general rules and limits for converting energy of different types, classically applied to conversions that involve heat such as heat engines. The first law of thermodynamics states that the change of total energy of a thermodynamic system is equal to the exchange in heat and the amount of work that is performed by the system. Due to the second law of thermodynamics, the entropy S (J K −1 ) in a thermodynamic system cannot be consumed, remains constant during reversible processes and is produced during irreversible mixing processes. The second law explains why heat flows naturally from warm to cold and thereby depletes temperature gradients: the reversal would consume entropy (Kondepudi and Prigogine, 1998). Isolated systems, which exchange neither energy nor mass nor entropy with their environment, evolve due to the second law to a "dead state" called thermodynamic equilibrium (TE), where entropy is maximized due to perfect mixing.
It is possible for heat to flow from cold to warm. However, this requires performance of work to steepen the temperature gradient as for instance realized in a refrigerator or a heat pump. Free energy is in the Oxford Dictionary defined as a thermodynamic quantity equivalent to the capacity of a system to perform work: i.e. to accelerate a (water) mass (as overland flow), to lift a (water) mass against gravity (as capillary rise) or to enlarge a potential gradient. Forms of free energy are kinetic energy, gravitational potential energy, electrical energy, chemical energy, and nuclear energy. Forms of free energy can ideally be entirely converted into another form of free energy. In reality these conversions are always associated with small amounts of dissipation. Heat is energy but not free energy. Heat can only be converted into work/free energy if an apparent temperature gradient can be exploited, but the converted amount is limited by the Carnot efficiency (T warm − T cold )/T warm (Kondepudi and Prigogine, 1998). Free energy can of course also entirely be dissipated into heat. Free energy is thus not a conservative property as dissipation is a "one-way avenue". Free energy is minimal in TE, as entropy and thus dissipation have been maximized.

Linking mass, energy and entropy balance of catchments
Catchments or hillslopes are open thermodynamic systems because they exchange mass, energy and entropy with their environment. Macropores, those organized structures discussed here, form and persist in the critical zone, which consists of the land surface including vegetation and the unsaturated zone down to the groundwater surface. Formation of organized structures and stabilization of such a system configuration far away from the TE require (a) that incoming energy and mass fluxes provide the necessary free energy to form and maintain structures, and (b) outgoing energy and mass fluxes export the entropy which is produced during irreversible processes to the environment (Kondepudi and Prigogine, 1998). Energy and mass inputs that sustain disequilibrium in the critical zone occur in two alternating regimes as radiation and rainfall (Figs. 1 and 2), thereby creating frequent oscillations of soil water potential around the local thermodynamic equilibrium (LTE)/hydraulic equilibrium. The amplitude and relaxation time of these oscillations depend on climate, vegetation, water retention curves of the soils (determining the "reset forces") and their flow resistance. The latter is controlled by unsaturated soil hydraulic conductivity and macropores. Depletion of re-established gradients in soil water potential implies dissipation of free energy and entropy production, because soil hydraulic conductivity specifies dissipative losses due to drag forces in the porous medium. Entropy production is, however, dominated by evapotranspiration, because energy conversions during the surface energy balance are much larger than those associated with rainfall-runoff processes. Evaporation is strongly dissipative due to the large specific heat of vaporization (Kleidon, 2012), thereby depleting steep surface gradients in air temperature and humidity (Fig. 1b). Entropy export from the critical zone is sustained by outgoing long-wave radiation and turbulent heat fluxes (Kleidon, 2012).

Free energy dynamics of soil water and the role of macropores
Free energy of soil water E soil chem (chemical energy according to Kleidon and Schymanski, 2008) is the sum of potential energy E soil pot (J) and the (negative) capillary binding energy E soil cap (J). Capillary pressure, the product of matric potential ψ (m) and the specific weight of water ρg, (kg (m 2 s 2 ) −1 ), determines how capillary binding energy changes with changing volume of the soil water phase V θ (m 3 ). Similarly we define hydrostatic pressure of soil water as the product of gravity potential z (m) and the specific

Fig. 2.
Local thermodynamic equilibrium in soil as state of zero soil water potential. Drying causes deviations to negative values by increasing absolute values of negative capillary binding energy; wetting causes deviations to positive values by increasing potential energy of soil water. Note that the change in capillary binding energy per unit soil moisture change scales non-linearly with the slope of the soil water retention curve, which is a fingerprint of the pore size distribution in soil (compare Fig. 3).
weight of water as change of potential energy per volume change of the soil water phase: As the groundwater surface defines the next relevant minimum of potential energy and matric potential is constantly zero there too, it is the natural choice for defining the reference level for potential energy in the critical zone (note that this selection is essential for defining soil hydraulic equilibrium). Free energy of soil water is thus indeed equal to the integral of soil water potential over the entire volume of the water phase, and its absolute value is equal to zero in soil hydraulic equilibrium/LTE (Fig. 2). As the volume of the water phase is the product of soil moisture θ (m 3 m −3 ) and a small control volume dV where soil moisture is spatially homogeneous, we obtain Infiltrating water from rainfall shifts near-surface soil water potential from LTE to positive values, which implies downward water fluxes. These reduce the excess potential energy in the system and shift the system back towards LTE (Figs. 1a and 2). Open-ended vertical macropores that connect either to the groundwater surface or via lateral preferential flow paths to the river body accelerate this reduction of potential energy in soil. They act as drainage structures by facilitating export of excess water. Note that this export stops at soil hydraulic equilibrium/LTE. Reduction rates of potential energy scale linearly with the drainage amount and topographic elevation above the groundwater level (Fig. 2).
Evapotranspiration decreases near-surface soil water potentials to strongly negative values, thereby enlarging absolute values of negative capillary binding energy (Figs. 1b and 2). This implies capillary rise, where upward water flows perform work against gravity (Fig. 2). Near-surface gradients in soil matric potential are, in cohesive soils during dry conditions, much steeper than the gradient in geo-potential. This is because the increase of the absolute value of capillary binding energy associated with a decrease of soil water content θ (m 3 m −3 ) increases non-linearly with decreasing soil moisture, reflecting the strongly increasing slope of the soil water retention curve C(θ ) = ∂ψ/∂θ (compare Fig. 3). Soil hydraulic conductivity in the matrix of cohesive soils decreases even stronger with decreasing soil moisture: net infiltration or capillary rise into such a dry "bottleneck" is thus very small (Zehe and Blöschl, 2004). Relaxation back to LTE by matrix flow is thus slow in cohesive soils (Zehe et al., 2010a). Dead-ended vertical macropores, worm burrows or cracks, allow infiltrating surface/rainfall water to bypass this bottleneck, thereby establishing saturated conditions at the macropore end ( Fig. 2) and successive wetting of the surrounding dry soil matrix. This process enhances depletion of steep gradients in matric potential, thereby shifting the system faster back towards LTE by reducing the strength of capillary binding energy.

Thermodynamic optimality and process trade-offs
in the critical zone

Macropores as wetting or drainage structures and related trade-offs
From a thermodynamic perspective plants act as "drainage structures". During long dry spells they "push" the root zone far away from LTE to states of strongly negative capillary binding energies, because their roots may extract water against large suction heads. Enhanced wetting and water storage is favourable for reducing the absolute value of capillary binding energies, especially in fine porous cohesive soils where capillarity dominates. This implies on average faster relaxation back to LTE. In contrast, enhanced reduction of potential energy due to enhanced drainage might be favourable to reduce free energy in non-cohesive soils -especially in steep terrain, where potential energy is the dominating term. With respect to these two "functions" we distinguish dead-ended macropores acting as wetting structures (similar to arteries) from open-ended macropores acting as drainage structures (similar to veins). Both types of macropores enlarge soil water flux against a driving gradient by reducing the flow resistance in the control volume. This reduction stems from the organized arrangement of soil material assuring connectedness of the flow paths, rather than from a change in the material properties.

E. Zehe et al.: A thermodynamic approach to link self-organization
The fact that depletion of gradients in capillary binding energy in cohesive soils is non-linear and facilitated by wetting structures implies a trade-off. Changes of the absolute values of capillary binding energy and potential energy during wetting/drying events occur always in opposite directions. An elevated density of macropores enlarges the wetting rate in the top soil and thus also the long-term average soil water content. This implies a long-term decreasing reduction rate of capillary binding energy associated with the same wetting amount. The potential energy change rate at a given depth depends in contrast exclusively and linearly on the wetting rate. This trade-off thus causes a long-term negative feedback that could result in an optimum density of dead-ended macropores that maximizes average reduction of free energy of soil water during rainfall-runoff processes. The break-even point is reached when the reduction rates of capillary binding energy are compensated by the increase rates of potential energy, because the slope of the soil water retention curve C(θ) = ∂ψ/∂θ in a soil layer becomes equal to its elevation above the groundwater surface. The former determines how strongly capillary binding energy is reduced, the latter how strongly potential energy increase with a certain wetting amount. This break-even point separates two different thermodynamic regimes in soil water dynamics: dominance of capillary binding energy, which implies that total free energy during rainfall-driven conditions is decreasing; from potential energy dominance, when total free energy of soil water is increasing during rainfall-driven conditions.
Closely related to this trade-off is the trade-off between soil water fluxes and gradients in soil water potentials. Elevated steady-state water fluxes due to an elevated density of macropores increase average near-surface soil moisture, thereby reducing the steepness of the gradient in soil water potential. Again, this implies a slow negative feedback, which could result in an optimum density maximizing the product of soil water flow times the potential gradient. Note that this product is due to the definition of power in fluid flow (as flow times pressure), in fact power in soil water flow. Also steady-state entropy production is in the literature often defined as being equal to the product of mass flow and the driving gradient, i.e. dissipation, divided by the absolute temperature (Porada et al., 2011;Kleidon and Schymanski, 2008).
Other studies have explored the second trade-off to search for an optimum soil hydraulic conductivity that maximizes entropy production or power (Kleidon and Schymanski, 2008;Porada et al., 2011;Westhoff and Zehe, 2013). Here, we propose a different avenue to search for an optimal macropore density, based on analysing the free energy balance of hillslopes during the rainfall-runoff process. The annual timescale is deemed as sufficiently long to analyse steady-state behaviour and to neglect changes in soil hydraulic properties (due to soil development) and nonstationarity in the density of biotic macropores.

Free energy balance of the critical zone during the rainfall-runoff process
Temporal changes in free energy density dE free /dt (W m −2 ) of the critical zone is determined by the rate of free energy input from the environment π in (W m −2 ) minus the rate of free energy that is dissipated D (W m −2 ) by irreversible processes and the export rate of free energy to the environment J out (W m −2 ): In steady state the storage term on the left-hand side of Eq.
(3) is zero; the entire free energy input is either dissipated or exported. A positive free energy balance implies that the surplus can be used to perform additional work for instance to create structures in a catchment/hillslope (Kleidon et al., 2013). The "environment" of the critical zone consists of the atmospheric boundary layer, the groundwater body and the river (at the hillslope scale) or the catchment outlet (at the catchment scale). Free storage in the critical zone is mainly storage of chemical energy (E Soil chem ) while free energy input is due to rainfall or radiation input. Free energy export (and entropy export) to the atmosphere is sustained by the surface energy balance components, free energy export in the form of kinetic energy occurs in the form of surface runoff or discharge. Our focus is on the part of the free energy balance associated with changes of soil moisture and the rainfall-runoff processes ( Fig. 1a): input of potential energy by rainfall, export of kinetic energy by surface runoff and export of potential energy by runoff components at distinct geo-potential levels. Thus Eq. (3) can be written as (ignoring D) The storage element for mass and energy in the critical zone is the soil, and free energy of soil water changes with changing soil moisture. The left-hand side of Eq. (3b) is hence the temporal derivative of Eq. (2) divided by the surface area of the hillslope/catchment A s (note we applied the chain rule to the term ψ(θ )): The first term in Eq. (4) implies that dE soil chem /dt depends nonlinearly on ∂θ/∂t. This is of key importance for the existence of possible optima in macropore density as explained in the next subsection. Free energy input from incoming rainfall P i (m s −1 ), the first term on right-hand side of Eq. (3b), depends on the incoming mass flux ρP i (kg (m 2 s) −1 ) and the geopotential height of a surface element gz s : This expression needs to be averaged across the entire surface area as neither gz s nor rainfall is homogeneous. Export of kinetic energy depends on the outgoing mass flux ρq OFL (kg (m 2 s) −1 ) of surface runoff and the squared flow velocity (m s −1 ). This is equal to q 2 OFL , because the specific runoff q R is the ratio of the volumetric flow divided by the wetted cross section. We obtain as the second term in Eq. (3b) Note that we neglect kinetic energy export by subsurface runoff components due to their small flow velocities. Specific mass losses and potential energy export from the system are finally determined by specific surface runoff q OFL , specific subsurface storm flow q SSF and groundwater recharge q GW , leaving the critical zone either at the geo-potential elevation of the river bank gz R , or the layer(s) gz L where interflow seeps into the river, or at the lower boundary gz 0 through a bottom element dA B (second line in Eq. 6). By inserting Eqs. (4)-(6) into Eq. (3b) we obtain the part of the normalized free energy balance of the critical zone that is determined by the rainfall-runoff processes: Please note that Eq. (7) mainly links specific rainfall-runoff components with appropriate coefficients to turn them into free energy fluxes and stocks. Equation (7) is not the complete free energy balance and thus not closed because the radiation/surface energy balance is missing. The left-hand side only accounts for changes in free energy due to changes in water mass and ignores those arising from temperature changes. Nevertheless the left-hand side is, as explained, sensitive to drying by evapotranspiration, which is a key component of the surface energy balance. Closing of the free energy balance equation requires furthermore explicit treatment of dissipation, because free energy is not a conserved property. Dissipation of free energy by surface runoff can be estimated by the geo-potential difference of the upslope source areas and the river bank, which determines the upper limit of the flow velocity. The squares of the upper limit compared with the actual flow velocities yield an estimate of dissipation. Steady-state entropy production is, in case kinetic energy export can be neglected, defined as the product of flow and the driving gradient divided by the absolute temperature (Porada et al., 2011). As kinetic energy of soil water flow can indeed be neglected, steady-state dissipation of free energy of soil water is equal to flow times the gradient in soil water potential, because dissipation is equal to entropy production times the absolute temperature at which this happens.

Maximum reduction of free energy of soil water
Steady-state conditions in Eq. (7) are reached at a timescale where changes in soil moisture storage can be neglected. We suggest this is the annual timescale as annual changes in soil water content are often close to zero, because of the repetitive annual cycles in short-wave radiation input, rainfall input and plant phenology. We aim to find an optimal macropore density (n * macro ) by focusing on free energy dynamics during rainfall-driven conditions, since this is when macropores accelerate infiltration and soil water flows, thereby reducing strengths of capillary binding energies and increasing potential energy of soil water. These two concurring processes determine in sum the change rates of total free energy of soil water E soil chem : (8) where dE soil chem /dt is determined by Eq. (4). Note that the change rate in E soil chem is in fact a second-order polynomial in θ, which depends implicitly on macropore density n macro . The quadratic term depends on the slope of the soil water retention curve, reflecting soil water retention properties and average soil moisture dynamics, which in turn are affected by macropores and climate. The first linear term is determined by matric potential, reflecting soil moisture (affected by macropores and climate) and the soil water retention curve. The last linear term depends on elevation above the groundwater level, reflecting the topographic setting and soil drainage properties (via depth to groundwater). In a given landscape this term gains weight when moving upslope.
In case free energy dynamics of soil water is dominated by dynamics of capillary binding energy, absolute values of E soil chem are, as explained in Sect. 2.3.1, decreasing during rainfall-runoff events. This implies relaxation back towards LTE. We expect this to happen in landscapes with cohesive soils, especially during summer conditions. The quadratic term in Eq. (8) dominates, which implies that two optima n * macro (1 m −2 ) could exist because a second-order polynomial may have two roots. These optima maximize the average reduction rate of free energy in soil water E soil chem during rainfall-driven conditions. This coincides with maximum energy dissipation by soil water flows, because the reduction rate in free energy reflects how fast the flux depletes the gradient. Please note that it is not unusual that several optima in environmental fluid mechanics exist (Kleidon et al., 2013); the most popular one is the hydraulic jump.
In the contrary case, when free energy dynamics is dominated by potential energy dynamics, absolute E soil chem is increasing during rainfall-runoff events. This implies that the system is pushed farther away from LET. We expect this to happen in landscapes with highly permeable, non-cohesive soils and steep topography, which implies that the third term dominates Eq. (8). Average increasing rates in total free energy are expected to change linearly with increasing macroporosity: there might be at best an optimum macroporosity at the upper margin of the parameter space. However, still a distinguished macroporosity can be defined as explained in the following. Due to the absence of overland flow production in such an environment the free energy input from rainfall contributes exclusively to the increase of potential energy of soil water (we neglect changes in capillary binding energy). Potential energy export is determined by subsurface storm flow and groundwater recharge, which are both affected by macroporosity. We can define a special/distinguished macroporosity as the one that leads on average to a steady state in the free energy/potential energy balance of the soil during the rainfall-runoff process. This is the case when the average export rates of potential energy are equal to the average increase rates in free energy storage in soil: Note that a steady-state free energy balance does not imply that average runoff equals average rainfall, because free energy export happens at a non-zero geo-potential. Next, we approximate change rates in total free energy by change rates in potential energy (second line of Eq. 9). We now use the fact that the integral of this first term is linear, the derivative and the integral can be exchanged and the mean is a linear operator as well. This allows us to derive an equation that sets the ratio of the mean change of areal averaged specific storage height (the vertical integral the of soil water content) over the mean runoff components equal to the ratio of the geo-potential levels of groundwater recharge and subsurface storm flow and the integral averaged geo-potential of upslope storage: where z 0.5 is the average depth over which soil moisture changes occur. Eq. (10), which approximates the mean of a product by the product of the means, estimates the steadystate storage-to-runoff ratio during rainfall-runoff conditions based on the ratio of the geo-potential levels where potential energy is exported and the upslope average geo-potential where storage happens. One minus this fraction yields an estimate of the annual runoff coefficient. Note that the proposed estimate relies essentially on the assumption that dynamics in capillary binding energies in soil can be neglected, and that the mean of the product in the first line of Eq. (10) can be approximated by the product of the means, which is valid as geo-potential changes at a much longer timescale than storage and discharge.
We can thus explore the existence of an optimum/distinguished density of macropores by analysing matric potentials, soil moisture and runoff components simulated within numerical experiments for varying macropore density, by means of Eqs. (8) and (10). We expect Eq. (10) to apply at our second test area (Malalcahuello), which is located in Chile, and the non-linear term in Eq. (8) to dominate in the Weiherbach catchment.

Scope and outline of numerical experiments
Based on the proposed framework we explore the existence of an optimum/distinguished density of macropores by analysing soil water flows, soil moisture and matric potentials from long-term numerical experiments with the physically based hydrological model CATFLOW. To this end we perturb the macropore density of "behavioural" model structures, which have been previously shown to allow successful long-term simulations of the water balances of the Weiherbach and Malalcahuello catchments. Both catchments are located in distinctly different hydro-climatic settings. Apparent soils differ strongly with respect to the degree to which capillary binding energy (Weiherbach Fig. 3a1 and a2) or potential energy (Malalcahuello Fig. 3b1 and b2) dominates free energy of soil water. Also the mechanisms causing vertical preferential flow are entirely different. In the Weiherbach catchment we expect that several optimal macropore densities can be defined based on Eq. (8). A generic test for our key hypothesis H2 is whether at least one of these optima allows acceptable uncalibrated predictions of rainfall-runoff behaviour. In the Malalcahuello catchment we search for a distinguished macroporosity as defined in Eq. (10). Also here we will test the predictive power of this macroporosity within an uncalibrated simulation of rainfall-runoff behaviour. The performance of this uncalibrated simulation and the goodness of the estimated annual runoff coefficient can be regarded as a test of H2 in the Malalcahuello catchment. Additionally we revisit numerical experiments carried out by Zehe et al. (2010a) for a small hillslope stretch close the Weiherbach brook. This is to shed light on how reduction of free energy in soil water during a single rainfall event increases with increasing surface density of macropores and whether this increase becomes, due to the above-explained trade-off, indeed smaller with increasing wetness of the soil.

The Weiherbach catchment
Geologically, the Weiherbach valley (southwestern Germany, 3.5 km 2 ) consists of Keuper and loess layers of up to 15 m thickness. The climate is semi-humid with an annual precipitation of 800 mm yr −1 , annual runoff of 150 mm yr −1 , and annual potential evapotranspiration of 775 mm yr −1 . More than 95 % of the total catchment area is used for cultivation 4306 E. Zehe et al.: A thermodynamic approach to link self-organization of agricultural crops, 4 % is forested and 1 % is paved area. The geomorphology is characterized by gentle slopes with a typical catena formed by erosion: moist but drained Colluvisols at the foot of the hillslopes and valley floors as well as drier Calcaric Regosols in the upper hillslope sectors. Capillary binding energy dominates free energy of soil water in these soils, which have a large amount of medium size pores. The slope of the soil water retention curve of, for instance, the Colluvisol ranges from values of 10 m close to saturation, over 50 m at 0.2 m 3 m −3 soil water content up to values larger than 1000 m when soil moisture decreases below a value of 0.15 m 3 m −3 (Fig. 3a2).
The Weiherbach catchment was intensely monitored in the period 1989-1996. This included (a) a soil and crop survey, derivation of typical soil hydraulic properties based on more than 200 soil cores (Table 1), and growth curves for arable crops; (b) monitoring of rainfall and runoff at six rain and two river gauges; and (c) observation of standard meteorological variables and ET from eddy correlation data. A soil map was compiled from texture information that was available on a regular grid of 50 m spacing. The macropore system consists of anecic earthworm burrows (Lumbricus terrestris and Aporrectodea longa) that build mostly deadended macropores. A survey of the macroporosity at 15 sites revealed a higher spatial density of larger and deeper worm burrows within the downslope Colluvisols (Zehe and Flühler, 2001;. One likely reason for this spatial organization is that earthworms find better habitat conditions due to the smoother moisture regime and the higher amount of soil organic material within the Colluvisols. These typical organized patterns of soils and macroporosity have caused a spatially organized pattern of infiltration at the slope scale (Fig. 3a, Zehe and Flühler, 2001), which in turn exerts a key influence on the catchment scale runoff response to extreme rainfall events . Due to the absence of lateral subsurface strata and the huge loess layers, Hortonian overland flow dominates event runoff generation. Because the apparent macropores elevate infiltrability of soils, rainfall-runoff events are rare. Event runoff coefficients range from 12 % during extreme thunder storms down to only 2 %. Base flow is almost constant throughout the year.

The Malalcahuello catchment
The Malalcahuello catchment is located on the southern slope of Volcano Lonquimay within the Reserva Forestal Malalcahuello (Fig. 3b) in the Chilean Andes. The catchment has a size of 6.26 km 2 and ranges in elevation from 1120 to 1856 m with mean slopes of 51 % (27 • ). Eighty (80) % of the catchment is covered with native forest of the types Araucaria and Nothofagus with a dense understorey of bamboo (Chusquea culeou). The climate is temperate humid with yearly rainfall ranging from 2000 to 3000 mm yr −1 . The soils are young, little developed volcanic ash soils with high porosities (60-80 %) and saturated hydraulic conductivities  van Genuchten, 1980, andMualem, 1976): saturated hydraulic conductivity k s , porosity θ s , residual water content θ r , air entry value α, shape parameter n.
Calc. (10 −5 -10 −3 m s −1 ). Note that capillary binding energy in these soils is rather small; matric potentials at 50 % saturation range between 0.1 and 1 m (Fig. 3b1), which is one to two orders of magnitude smaller than in the Weiherbach soils (Fig. 3a1). We expect thus that potential energy dominates free energy of soil water here.
Within an extensive field campaign, Blume et al. (2007Blume et al. ( , 2008a instrumented a transect in a focus area with soil moisture probes, shallow piezometers, and throughfall collectors and installed several rain and stream gauges in the catchment. These instruments yielded continuous data for the period 2004 to 2006. These observations were complemented by auger and electric resistivity tomography (ERT) profiles to explore strata and depth to bedrock, derivation of soil hydraulic parameters from soil cores, observations of stable isotopes, river chemistry, distributed dye tracer profiles and much more. Existing data prior to this study included data from two nearby climate stations where data have been measured (a) on a daily basis since 1989 and (b) on an hourly basis since 1999.
Note that runoff production within the Malalcahuello catchment functions distinctly different from the Weiherbach. Annual runoff coefficients are between 60 and 70 % and thus three times larger than in the Weiherbach, while event runoff coefficients are within 4 to 19 %: similar to the Weiherbach. Overland flow is of negligible importance; subsurface storm flow and mobilization of pre-event water from the riparian zone or the lower hillslope sector are the dominant runoff processes (Blume et al., 2008b). The vertical trigger is fast finger flow caused by the fairly strong hydrophobicity of these soils in summer (Blume et al., 2009). In this setting we expect that potential energy is largely dominating free energy of soil water. We thus search for a distinguished macropore density as defined in Eq. (10), which implies that the annual runoff coefficient estimated using Eq. (10) should compare to the observed values of 60-70 %.

Process representations in the CATFLOW model
The numerical experiments were conducted with the physically based model CATFLOW . The model represents a hillslope as a 2-dimensional cross section along the steepest descent line that is discretized by 2dimensional curvilinear coordinates. Hillslopes are assumed to be uniform perpendicular to the slope line. Soil water dynamics are described by the Richards equation in the potential form that is numerically solved by an implicit mass conservative Picard iteration (Celia and Bouloutas, 1990). The model can thus also account for soil water flows under saturated conditions . Soil hydraulic functions are described after van Genuchten (1980) and Mualem (1976). Evapotranspiration is represented by an advanced SVAT (Soil Atmosphere Vegetation Transfer) approach based on the Penman-Monteith equation, which accounts for growth of agricultural crops by annual cycles of LAI (Leaf Area Index) and ground cover of perennial plants.
Surface runoff routing down the hillslopes and flood routing in the river network is simulated based on the diffusion wave approximation of the Saint-Venant equations, which are numerically solved by an explicit upstream finite difference scheme. CATLFOW allows representation of vertical and lateral preferential pathways along two avenues. One is an explicit representation as connected flow paths of low flow resistance and low retention properties (Zehe et al., 2010a). This approach was shown to be suitable to predict preferential flow and tracer transport in the Weiherbach catchment Zehe, 2010, 2011) and a catchment in the Vorarlberg Alps (Wienhoefer and Zehe, 2013). The advantage of this approach is that it may be parameterized using observed data on density of worm burrows/macropores, their lengths and maximum water flow rates. Due to its computational expense, it is restricted to small spatial scales and short timescales. We thus use this approach within the small-scale numerical experiments focusing on single events.
For the experiments on the annual timescale we will thus use an effective approach that enlarges hydraulic conductivity by the macroporosity factor f mac when soil water saturation S exceeds a threshold S 0. Bulk hydraulic conductivity k B is then increased linearly with relative saturation: The macroporosity factor is the water flow in all macropores in a model element where the saturation is larger than the threshold S 0 , divided by the water flow rate in the soil matrix. Generally one can chose whether macroporosity affects only vertical soil water flows or whether it works in both directions. At small scales f mac can also be determined based on data on density of worm burrows, their lengths and maximum water flow rates (Zehe and Blöschl, 2004). In largescale studies we use f mac however as calibration parameter. This has been shown to work well during continuous simulation of the water balances of the Weiherbach catchment  and the Malalcahuello catchment (Blume, 2008). The range of f mac as a calibration parameter depends on whether the macropores are dead ended or open ended.
The former case implies that matrix hydraulic conductivity limits flow at the end of the macropores and across the macropore walls (Beven and Clarke, 1986). Hence, for the Weiherbach soils f mac ranges from 1 to 4. When vertical and lateral preferential pathways are interconnected and extend to the stream, as in Malalcahuello, matrix hydraulic conductivity is not limiting. In this case f mac values range up to values of 54 to reproduce observed flow reactions (Blume, 2008).

Numerical experiments: setup and results
4.1 Experiment 1: free energy reduction at the event scale as function of macropore density and initial saturation

Basic model setup and perturbations
In this experiment we explore whether the dynamics of free energy of soil water is accelerated by an increasing macropore density, and to which extent this sensitivity is reduced at higher initial saturations. This is to corroborate the existence of the trade-off discussed in Sect. 2.3.1, which in fact separates the two different thermodynamic regimes in soil water dynamics. To this end we employ a setup that has been proposed by Zehe et al. (2010a) to explore how dynamics of free energy of soil water is accelerated by an increasing macropore density during block rainfall events at the daily timescale. The simulation domain represents a field site in the meadows close to the Weiherbach. The stretch is 30 m long and 2 m deep, and total elevation difference is 2 m. Soil hydraulic functions of the apparent Colluvisol are listed in Table 1. The lateral and vertical grid width was 2 cm. Vertical worm burrows were represented as explicit connected structures, generated with a Poisson process and a stochastic agent which mimics the digging behaviour of the earthworms (Fig. 4a). The database for this consisted of observations of areal density and depth distribution of worm burrows at this field site. We generated five different setups that differed exclusively with respect to the areal density of worm burrows (5, 10, 15, 20 and 25 m 2 ); observed average macropore depth at this site was 0.8 m with a standard deviation of 0.2 m. Zehe and Blöschl (2004) determined the maximum water flux in worm burrows based on lab experiments as 1 × 10 −3 m s −1 . Additional details on the representation of macropores in the model can be found in Zehe et al. (2010a). We simulated soil water dynamics and rainfall-runoff of these hillslopes using a block rain of 10 mm in 2 h at different initial saturations, S, of 0.4, 0.6, and 0.8 and analysed the free energy balance using Eq. (8). . On the right one can see the sequence of soil layers; the red arrow symbolizes the large macroporosity factor that had to be assigned to the top 80 cm to fit observed specific discharge.

Results: event scale free energy balance and process trade-offs
Results of experiment 1 are presented in Figs. 5 and 6. Changes in free energy in the dryer experiment are strongly dominated by the reduction of the absolute value of capillary binding energy due to the accelerated wetting in macropores which peak at −1 W m −2 (Fig. 5f). The associated increase of potential energy is one order of magnitude smaller. Free energy of soil water is strongly negative and dominated by capillary binding energy; the absolute values are 5 to 6 times larger than potential energy of soil water (Fig. 5h). An increase of the initial saturation by 0.2 to 0.8 results in much flatter reduction rates of capillary binding energy, with absolute values only slightly larger than the increase rates in potential energy (Fig. 5e). Free energy of soil water is already positive at the beginning and only slightly increases during the course of the experiment. Note that the increase rates of potential energy are for both initial saturations nearly the same, because the net wetting amounts are pretty similar. The flattened reduction rate in capillary binding energy is thus indeed due to the fact that the slope of the soil water retention curve C(θ ) decreases with increasing wetness (compare also Fig. 3a2). This trade-off is even better illustrated in Fig. 6, which presents averaged reduction rates of capillary binding energy (blue dots, panels a, c, e) and increases of potential energy (black dots, panels b, d, f) which sum up to reduction rates of free energy of soil water (red dots, panels a, c, e). Reduction rates in capillary binding energy decrease by three orders of magnitude when increasing initial saturation from 0.4 to 0.8 (compare top, middle and bottom panels), while increases in potential energy remain much more stable.
We may thus state that the sensitivity of the reduction rates in free energy of soil water on macropore density does indeed decrease with increasing wetness, because the decreasing depletion rates in capillary binding energy are more and more compensated by the increases in potential energy of soil water. As proposed in Sects. 2.3.1 and 2.3.3 total free energy of soil water in cohesive soils (as the Weiherbach Colluvisol) is decreasing during rainfall-driven conditions, because reduction rates of capillary binding energies (in fact of their absolute values because they are negative) dominate against increases of potential energy. This dominance is strong, because the stretch under investigation is rather flat. As capillary binding energy is the dominating factor, we explore in the next subsection whether optima in macroporosity that maximize averaged reduction rates of free energy exist.

Basic model setup and perturbations
For the simulation in the Weiherbach catchment we perturbed a setup introduced in Zehe et al. (2001). The catchment was subdivided into 169 hillslopes and an associated drainage channel network (Fig. 4b). The hillslope model elements are 5-20 m wide and 10 m long, and vertical grid size varies from 10 cm close to the surface to 60 cm close to the lower boundary. Total soil depth represented by the model was 2 m. Manning roughness coefficients for the hillslopes and the channels were taken from irrigation experiments performed in the catchment as well as from the literature (Scherer et al., 2012). We used free drainage as the lower boundary condition for the hillslopes, a seepage boundary condition at the interface to the stream, an atmospheric boundary condition at the top and no flux at the wa- , temporally averaged change rates in potential energy, the absolute value of capillary binding energy and free energy of soil water (panels e and f, compare Eq. 7) and the time series of free energy components of soil water (panels g and f, after Eq. 2). The left panels a, c, e, and g belong to the wet initial saturation of S = 0.8, the right panels to the dryer initial saturation of 0.6. Please note that kinetic energy export by surface runoff deceases with macropore density, while changes in capillary binding energy increase with macropore density. , and (g) belong to the wet initial saturation of S = 0.8, the right panels to the dryer initial saturation of 0.6. Please note that kinetic energy export by surface runoff deceases with macropore density, while changes in capillary binding energy increase with macropore density. Time series of potential energy input by rainfall (a and b), kinetic energy export by surface runoff (c and d), temporally averaged change rates in potential energy, the absolute value of capillary binding energy and free energy of soil water (e and f; compare Eq. 7) and the time series of free energy components of soil water (g and f, after Eq. 2).
upper 80 % and Colluvisol in the lower 20 % of the hill. The corresponding van Genuchten-Mualem parameters are listed in Table 1. Vegetation was parameterized based on available crop patterns in the simulation period (April 1994 to September 1995) and tabulated data with typical annual cycles for LAI, plant height, root depth and values for stomata resistances.
Observations of macroporosity suggested a downslope increasing density of worm burrows with normalized macropore volumes of 1.5 × 10 −3 m 3 m −2 at the hillfoot and typically 0.6 × 10 −3 m 3 m −2 at the hilltops (Zehe, 1999;his Fig. 4.1). We represented this organized pattern by assigning fixed scaled values of macroporosity to each hillslope: 0.6 × f mac at the upper 70 %, 1.1 ×f mac at the range from 70 to 85 % of the hillslope, and 1.5 × f mac at the lowest 85 to 100 % of the slope length. The depth of the macroporous layer was assumed to be constant throughout the entire catchment and was set to 1 m to simulate dead-ended macropores. Threshold saturation was set to S 0 = 0.8. As worm burrows extend mainly in the vertical direction, we chose the anisotropic version that works only in vertical flow direction. The macroporosity f mac is the remaining free parameter, which was calibrated to match the hydrograph of the largest observed flood event in June 1994: f mac = 1.5 allowed the best of (a) observed discharge during a 1.5 yr-long simulation with a Nash-Sutcliffe coefficient (NSE) of 0.8, (b) soil moisture dynamics observed at 61 TDR ((Time Domain Reflectometry)) stations (R 2 between 0.5 and 0.7), and of evaporation observed at the central meteorological station (R 2 = 0.9, . During the numerical experiment we simulated the water balance of a single hillslope, which was randomly selected from the catchment model. Simulation started at 21 April 1994 and lasted until 15 September 1995, because this is the only period where rainfall-runoff events with runoff coefficients larger than 2 % occurred. In total Figure 6: Averaged reduction in absolute values of capillary binding energy (blue dots, panels a, c, e), enlargement of potential energy (black dots, panels b, d, f) and reduction of free energy of soil water (red dots, panels a, c, e Eq. 8) for numerical experiment 1. Average initial saturation decreases from top to the bottom panels as indicated in panels b, d, and f.  dots, a, c, e), enlargement of potential energy (black dots, b, d, f) and reduction of free energy of soil water (red dots, a, c, e, Eq. 8) for numerical experiment 1. Average initial saturation decreases from top to the bottom panels as indicated in panels from 0.8 over 0.6 to 0.4. we performed 30 simulations, which differed with respect to f mac , which was incrementally enlarged from a value of 1 to a value of 2.5. This gradual increase corresponds to a higher maximum infiltrability and implies a higher density of activated macropores/worm burrows. Each model structure was dynamically initialized by using the final state of a twoyear simulation period (21 April 1992-20 April 1994), which started with hydraulic equilibrium in soil and was driven using observed boundary conditions.
As the Weiherbach hillslope is composed of two different cohesive soils and the macroporosity increases downslope with fixed scaled ranges, separate optimum macroporosity values after Eq. (8) could exist for each soil. Please note that the macroporosity factor scales hydraulic conductivities, which is in the Calcaric Regosol, with 2×10 −6 m s −1 , twenty times smaller than in the downslope Colluvisol. Additionally the weight of potential energy increases in most upslope areas, due to the elevation difference of 40 m compared to the hillfoot sector. Free energy dynamics of soil water averaged over the entire hillslope is thus a volume-weighted average of what happens in upslope Calcaric Regosol and in downslope Colluvisol. The resulting optima of averaged reduction rates of free energy in the entire hillslope can thus be expected to be less sharp than in a system with uniform soil and uniform macroporosity. Figure 7b presents the accumulated fluxes of the hillslope water balance for the 30 model runs. The water balance in the Weiherbach is dominated by evapotranspiration, followed by groundwater recharge and occasional direct runoff events by Hortonian overland flow. Higher f mac values have thus only a small effect on the accumulated water balance components, because they mainly reduce event scale Hortonian overland flow production during rare events and slightly increase storage during the overland flow events (note that the macropores are dead ended). The difference in accumulated groundwater recharge among the different simulations is thus as small as 15 mm. Figure 7d presents the specific catchment scale runoff and the best fit obtained with f mac = 1.5 for the largest flood event. Please note that the corresponding runoff coefficient was only 12 %. Figure 8 presents the forcing data and temporal change in free energy balance components for the Weiherbach calculated from Eq. (7) as a function of time. Comparison of the upper three to the lower three panels corroborates that free energy fluxes and changes in free energy stocks related to the (large) mass fluxes during rainfall-driven conditions are two orders of magnitude smaller than during radiation-driven conditions (when mass fluxes are small). Syst. Sci., 17, 4297-4322 Fig. 10 for the maxima). Please note that base flow in the Weiherbach catchment is so small that it cannot be detected in (d) and that difference in accumulated surface runoff for the different model runs are not easy to detect at this scale (simply because of the fact that the number of overland flow events is 12 within this period).

Hydrol. Earth
As expected free energy dynamics of soil water is in these loess soils clearly dominated by capillary binding energies (Fig. 8d). Free energy of soil water is strongly negative in the dry summer 1994 (day 0-100), then oscillates around local thermodynamic equilibrium/hydraulic equilibrium during the wetting period from fall 1994 to spring 1995 (day 110 to day 200) and becomes positive from day 200 to day 420, as spring and summer 1995 were quite wet. Strongly negative values of total free energy occur again in late summer 1995 around day 450. The temporal changes in free energy of soil water are also dominated by reduction rates in capillary binding energy, which peak at values of −4 W m −2 for heavy infiltration events during summer conditions. Associated increase rates of potential energy are 50 % smaller; they gain, however, importance during wet conditions (day 100-400). Export of potential energy by means of groundwater recharge is, with values smaller than 10 −4 W m −2 , unimportant. Export of kinetic energy by surface runoff is mostly small except for the largest overland flow event, on 27 June 1994 (day 62). Figure 9 provides the free energy balance components after Eq. (7) for this largest event. Kinetic energy export peaks at 1 W m −2 and decreases clearly within increasing macroporosity. Complementary to this, the reduction rates of the absolute value of capillary binding energy increase with increasing macroporosity (Fig. 9c), while potential energy changes are positive and increase with increasing macroporosity as well.

Results: optima in macroporosity and their
performance during rainfall-runoff simulation Figure 10 presents the free energy balance components of the 30 model runs as function of the macroporosity factor. Note that the values were averaged across all time steps where rainfall occurred. Average changes of total free energy of soil water are indeed negative (panel a), which corroborates that reduction rates of capillary binding energy (panel b, in fact of their absolute values) dominate against the increases in potential energy (panel c). Kinetic energy export by means of overland flow decreases with increasing f mac because overland flow production decreases with enhanced infiltrability (Fig. 10c). Increase rates in potential energy in soil water (Fig. 10d) show at first a steep increase when infiltrability in upslope areas is increasing (which needs larger f mac values), followed by a slower increase with increasing macroporosity. Reduction rates in the strength of capillary binding energy are on average 3 times larger than changes in potential energy and they peak at two local maxima at f mac values of 1.45 and 2.0 respectively. The location of the first optimum is very close to the macroporosity that yielded the best fit of observed rainfall-runoff behaviour. The local maxima are, as expected, not very sharp. This may be explained by the fact that reduction rates in the two different soils peak at slightly different f mac values. The reduction rates of total free energy show a very similar dependence on macroporosity to reduction rates of capillary binding energy; the local maxima are located at the same positions. We may thus state that free energy dynamics of soil water during the rainfall-runoff process is in the Weiherbach soils on average dominated by capillary binding energy. As  and simulated latent heat flux λET (c) as function of time, trajectories of capillary binding energy, potential energy and total free energy of soil water (d) (derived from Eq. 2, the curves in (g) are the numerical derivatives of the curves in (d) with respect to time). Potential energy input by rainfall (e), kinetic energy export by surface runoff and potential energy loss by means of groundwater recharge (f), as well as temporal changes in free energy, capillary binding energy and potential energy of soil water (g) (compare Eqs. 7 or 4). the quadratic nature of Eq. (8) dominates we find two local optima that maximize average reduction rates of free energy of soil water during rainfall-driven conditions, with values of 0.028 and 0.0285 W m −2 . To assess the predictive value of these optima we assigned the corresponding f mac values of 1.45 and 2.0 to all 169 hillslopes of the catchment model and simulated the water balance for 1.5 yr. Figure 7d compares simulated catchment scale discharge based on the optimum f mac values to the results of the best-fit setup for the largest flood event, observed on 27 June 1994. The simulation based on the first optimum performs almost as well (NSE = 0.79 for the entire period) as the best fit (NSE = 0.8). The simulation based on the second optimum produces as expected much less direct runoff due to overland flow (Fig. 7d). Please note that observed discharge is not part of the presented analysis of the energy balance components; therefore this simulation is indeed uncalibrated.

Basic model setup and perturbations
For the simulation of rainfall-runoff in the Malalcahuello catchment, the monitored hillslope near the catchment outlet was considered typical for the entire catchment. Due to the very fast travel times within the stream network (event  response times as low as 30 min have been observed), it seemed feasible to neglect transport times within the stream. The model therefore consisted of only one hillslope of 80 m length and a total elevation difference of 60 m. Specific hillslope outflow was compared to catchment-specific discharge for validation purposes. The hillslope domain had a depth of 10 m (Fig. 4b). The soil layers were parameterized according to field observations and consisted of a 5 cm humus layer followed by several layers of volcanic ash soil. At 2 m depth a layer of gravel-size volcanic material is likely to function as a poorly permeable capillary barrier. This was represented as layer of lower permeability (compare Table 2). The lower boundary condition is no flow at the lowest 15 % of the hillslope; the rest was free drainage. The right boundary condition is no flow in the lower part (up to 2.5 m below ground), producing a saturated wedge in the lower part of the hillslope, above which a seepage boundary allows for hillslope outflow. Vertical and lateral hydraulic conductivity in the upper 1 m was enlarged by a macroporosity factor. Best fit with a NSE of 0.7 was achieved with a value of 54 and a threshold saturation of S 0 = 0.4, which corresponds to field capacity in these non-cohesive soils. Note that macroporosity acts isotropic and thus also enhances flow into the lateral downslope direction. The modelled vegetation was deciduous forest with parameter values either observed in the field or taken from the literature. We performed 20 runs with the basic model setup that differed with respect to the macroporosity factor, which was gradually increased from 1 to a maximum of 70 within 20 steps. This corresponds to maximum hydraulic conductivity of 2.1×10 −2 m s −1 in the case of activated macropores. Each model structure was dynamically initialized by using the final state of a one-year spin-up period, which started with hydraulic equilibrium in soil and was driven using observed boundary conditions.
As already explained in Sect. 3.1 we expect no optimum macroporosity but a distinguished macroporosity to exist, which brings potential energy in the system on average closest to steady-state. This idea allows, based on Eq. (10), an  van Genuchten, 1980, andMualem, 1976). Saturated hydraulic conductivity k s , porosity θ s , residual water content θ r , air entry value α, shape parameter n. show the minus one times the change rate to better visualize the maxima located at fmac values of 1.45 and 2.05. Fig. 10. Temporally averaged change rates in free energy during rainfall-driven conditions for numerical experiment 2 in the Weiherbach plotted against the macroporosity factor: total free energy of soil water (a) (as defined in Eq. 8), capillary binding energy (b), kinetic energy export by surface runoff (c) and potential energy (d).
Note that panels (a) and (b) show the minus one times the change rate to better visualize the maxima located at f mac values of 1.45 and 2.05. estimate of the annual runoff coefficient as the ratio of the sum of the geo-potential levels where mass is released from the system to the upslope-volume-averaged geo-potential where storage occurs. To estimate the annual runoff coefficient based on Eq. (10) we define z 05 as 50 % of the depth where relative soil moisture changes during rainfall-driven conditions were less than 1 %. The depth of the interflow layer z SSF corresponds to the depth of the impermeable layer. Figure 7a presents the accumulated fluxes of the hillslope water balance for the 20 model runs, which evidently strongly differ from the water balance in the Weiherbach: accumulated rainfall is twice as large (compare with Fig. 7b), though the simulation period is just 6 months compared to 18 months for the Weiherbach. Evaporation plays a minor role compared to groundwater recharge and direct runoff, which is fed by subsurface storm flow. For small macroporosity factors the simulated water balance is dominated by groundwater  Trajectories of capillary binding energy, potential energy and free energy of soil water (c) (derived from Eq. 2). Potential energy input by rainfall (d), potential energy export by means of subsurface storm flow groundwater recharge (e), as well as change rates in free energy, capillary binding energy and potential energy of soil water (f) (compare Eq. 7). Note that the curves in (f) are the numerical temporal derivatives of the curves in (c).

Results: rainfall-runoff process and related free energy dynamics
runoff, followed by subsurface storm flow. At larger macroporosity values subsurface storm flow takes over the dominant role. The best fit of observed specific discharge was achieved with a macroporosity factor of 54 (Fig. 5c). Figure 11 presents the free energy balance components calculated from Eq. (7) as a function of time. As expected, free energy of soil water is clearly dominated by potential energy. Capillary binding energy plays a minor role (Fig. 11c). The regime is characterized by an excess of potential energy through the entire simulation period. The enormous drainage capabilities of these soils and evapotranspiration during radiation-driven conditions are not sufficient to deplete this excess during the simulation period. Changes in free energy of soil water during rainfall events are clearly positive, because increases in potential energy dominate against associated much smaller decreases in capillary binding energy. Increase rates in potential energy of soil water peak at 2 W m −2 . Export of potential energy by means of subsurface storm flow events dominates reduction of free energy during events, while groundwater recharge adds a smaller but continuous component to potential energy export.

Results: distinguished macroporosity and its
performance during rainfall-runoff simulation Figure 12 present the energy balance components after Eq. (9) plotted against macroporosity. Note these components were averaged over all time steps when rainfall was non-zero. The change rates in total free energy in soil water are on average positive during rainfall-driven conditions (panel a), and decrease with increasing macroporosity. The latter appears at first sight counterintuitive, but can be understood by the fact that lateral soil water flow dominates against vertical flows with increasing macroporosity. This shift towards lateral soil water flows is also visible from the fact that (negative) export of potential energy by subsurface  Fig. 12. Temporally averaged change rates of free energy during rainfall-driven conditions for numerical experiment 3 in Malalcahuello plotted against the macroporosity factor: change rates of total free energy in soil water and potential energy of soil water (a), potential energy export by groundwater recharge (b) and subsurface storm flow (c) as well as the sum of the changes rates in free energy storage in soil and free energy export by runoff (c) defined after Eq. (9). storm flow increases with increasing f mac (panel c), which is fed by lateral soil water flows at the lower boundary of the macroporous soil layer. Consistently with this, the export by groundwater flow is decreasing with increasing macroporosity (panel b). Panel d presents the sum of average increase rates in potential energy (storage increase) and the export rates by subsurface storm flow and groundwater flow after Eq. (10). The distinguished macroporosity for which this sum becomes zero is f mac = 50. This macroporosity factor is close to the value of 54 which yielded the best fit. To assess the predictive value of the distinguished macroporosity we assigned the corresponding f mac value of 50 to the hillslope model and simulated the water balance for the snow-free period in 2006. The corresponding Nash-Sutcliffe efficiency is 0.65. Please note that this f mac balances on average the increases in potential energy during rainfall-driven conditions with the potential energy losses, and brings the system closest to a steady state in energetic terms. As observed discharge is not part of the analysis, the corresponding simulation is again uncalibrated. In addition we estimated the steady-state runoff coefficient based on Eq. (10) to 0.74. This value is comparable to the range of 0.6 to 0.7 that Blume et al. (2008b) derived from observations. Figure 13 presents box plots of the spatially averaged soil water potentials in the root zone (free energy of soil water per volume divided by the specific weight of water ρg) for the two optima in the Weiherbach (WB f mac = 1.45 and 2.0 respectively) as well as for the best fit model setup and the dis-

Basic model setup and perturbation
In the last numerical experiment we simulated the water balance for the entire Weiherbach for the period from 21 April 1994 to 15 September 1995 for two different cases: -Case 1: the basic model setup based on the typical catena and the observed spatially organized macroporosity pattern with small values at the top and largest values in the valley. As an average we selected f mac = 1.45 of the first thermodynamic optimum (close to the best fit).
-Case 2: we assume that the earthworms have the opposite habitat preference, building more worm burrows upslope and less downslope. We flipped the spatial pattern of macroporosity thus from hilltop to hillfoot and kept everything else as in case 1 (the fixed scaled ranges, the average macroporosity and depth as well as the soil catena).
Each model structure was dynamically initialized by using the final state of a two-year simulation period (21 April 1992-20 April 1994, which started with hydraulic equilibrium in soil and was driven using observed boundary conditions. Figure 14 presents the simulated catchment scale discharge for both cases (Fig. 14c) for the largest rainfall-runoff event observed in June 1994 (Fig. 14a presents the corresponding rainfall event), the trajectories of free energy of soil water and its components (Fig. 14b) as well as the averaged reduction rates of free energy in soil water (panel d after Eq. 8). Apparently case 2 produces more surface runoff. Due to the reduced macroporosity in downslope reaches, less of the upslope-generated overland flow infiltrates at the hillfoot and can thus runoff from the system. At the same time case 2 has a clearly higher free energy of soil water during almost the entire simulation period, because of a clearly enlarged potential energy (Fig. 14b). Because of the enlarged macroporosity at the hilltops, more water infiltrates in these reaches which are on average 20 to 40 m above the hillfoot sectors. Potential energy is thus clearly enlarged. Consistently, reduction rates of free energy of soil water are larger in case 1 than in case 2. This suggests that a worm species with the opposite habitat preference is from a thermodynamic perspective less favourable, as the system is on average farther away from local thermodynamic equilibrium and relaxation times are longer.

Optimum density and pattern of wetting structures in the Weiherbach soils
The results from the numerical experiments 2 and 4 suggest, as far as we know for the first time, that the spatially organized patterns of soils and dead-ended macropores that have been observed in a real-world landscape are in close accordance with one of two apparent local thermodynamic optima. These optima are expressed by locally minimized relaxation times towards local thermodynamic equilibrium, which is equivalent to maximize reduction rates of total free energy of soil water during rainfall-driven condition. This holds for both the f mac factor, which represents the enhancement of infiltration by activated macropores per unit area, as well as the spatial arrangement of the macropores along the lateral potential gradient driving overland flow at the hillslope scale. Note that potential energy in case 2 is during the entire period clearly larger as for case 1, which explain the difference in free energy of soil water Fig. 14. Precipitation (a) and simulated catchment scale discharge (c) for case 1 (optimum macroporosity, downslope increasing macroporosity; solid blue) and case 2 (optimum, but upslope increasing macroporosity; solid red). The graphs are zoomed to the largest runoff event, observed on 27 June 1994. Observed discharge is in solid black. Time series of total free energy of soil water, capillary binding energy and potential energy for the entire period, and averaged reduction rates of free energy in soil water for both cases after Eq. (8) (d). Note that potential energy in case 2 is during the entire period clearly larger than for case 1, which explains the difference in free energy of soil water.

Optimum density of wetting structures in cohesive loess soils
In numerical experiment 2 we found two apparent maxima in free energy reduction of soil water for macroporosity factors of 1.45 and 2.0. The existence of these optima is explained by the fact that accelerated reduction of capillary binding energy during rainfall-driven conditions dominate against accelerated increases of potential energy of soil water, when increasing the macroporosity factor. Free energy of soil water is thus on average reduced during rainfall-driven conditions and the reduction rates of free energy are a second-order polynomial of the wetting rate. This implies that two roots may exist. Capillary binding energy dominates against potential energy if the slope of the soil water retention curve C(θ ), which non-linearly decreases with increasing soil moisture, is on average larger than the depth to groundwater in a certain layer. This is because the former determines how strongly matric potential decreases while the latter determines how strongly potential energy increases with a certain amount of wetting during rainfall events. The break-even point is reached when the reduction in capillary binding energy is compensated by the increase in potential energy of soil water, which is the case when the slope of soil water retention curve C(θ ) in a soil layer is equal to its elevation above groundwater. This trade-off limits the degree to which a higher amount of active dead-ended macropores, which in the Weiherbach are created by anecic earthworms, increases reduction rates of free energy of soil water during rainfall-driven conditions. In fact reduction of free energy rates start to slowly decrease again when increasing f mac beyond 2.0 (Fig. 10a).
It is striking that the first optimum is, with a value of 1.45. close to the "best-fit macroporosity" of 1.5 and that simulated catchment scale discharge based on this thermodynamic optimal macroporosity is in pretty good accordance with observed discharge response. Let us recall that re-establishment of disequilibrium in soil is controlled by the climate regime in complement with vegetation, topography and soil water retention properties (Fig. 2). Transpiration and evaporation play a key role in drying of cohesive soils that would otherwise store most of their water against gravity and stay close to local thermodynamic equilibrium. Roots act as a short cut through dry top soil layers, thereby increasing dryness of deeper soil reaches during summer conditions. Plant roots are in fact those drainage structures that may extract water against very steep gradients of capillary binding energy and push the system far away from local equilibrium conditions (Fig. 14). The distance of the long-term average of soil water potential in the root zone from local thermodynamic equilibrium (0.48 m in the Weiherbach) as well as the deviation to stronger negative values farther away from LTE is thus mainly controlled by transpiration. This is in turn driven by global radiation input and plant physiology. Relaxation back to equilibrium is controlled by the amount and depth of wetting structures, worm burrows here. Hence, both the deviation from and relaxation back to equilibrium are controlled by biotic structures of different origin that dominate during different prevailing boundary conditions.
We thus conclude that in systems where capillary energy dominates, soil wetting during rainfall-driven conditions shifts the system back towards local thermodynamic equilibrium, because total free energy of soil water is reduced. The results from numerical experiment 1 fully support hypothesis H1, which states that vertical macropores act as dissipative structures by enhancing the reduction of free energy of soil water. Experiment 1 corroborates that this reduction rate depends strongly on the system's wetness even if the wetting rate is constant, because reduction in capillary binding energy depends on the wetness state. Also, hypothesis H2 cannot be rejected. The behavioural model structure is the best functional and structural representation of the Weiherbach catchment in a model, as it represents the main structural and textural features of the Weiherbach soils and reproduces different water balance components. This behavioural model structure is with respect to the calibrated f mac value (1.5) close to the thermodynamic optimum. Partitioning of rainfall input into Hortonian overland flow and infiltration is hence close to the optimum that assures fastest reduction rates of free energy of soil water and minimum average deviation from LTE.
We thus conclude that the Weiherbach itself has "coevolved" (Schaefli et al., 2011) close to a configuration that is thermodynamically optimal given the climate conditions, land use, and hydro-pedological setting and topographic setting. The second optimum is even more efficient in reducing free energy. Rainfall-runoff simulations produce less direct runoff due to overland flow and thus less soil erosion, which is at present a pressing problem in the Weiherbach (Scherer et al., 2012). The existence of these two optima makes sense when recalling that the landscape is not static but slowly evolving (Schaefli et al., 2011) and steady states can thus at best be metastable. The first optimum can be interpreted as such a metastable landscape configuration from a thermodynamic perspective, which can obviously be realized within the factors that limit population dynamics of the ecosystem engineers creating the dead-ended macropores. Note that the macroporosity factor is related to the areal density of worm burrows and thus to the abundance of anecic earthworms (van Schaik et al., 2013). In this sense the second optimum can be regarded as a probable metastable state that could be reached in a later stage of landscape evolution, which would be associated with less surface runoff and less erosion. Testing of such a hypothesis requires, however, development of coupled models for hydrology, population dynamics of earthworms and even geomorphic processes, or in general terms, of models that account for the relevant interactions and feedbacks between biotic and abiotic processes and limitations for biotic structure formation and their controls on water fluxes and related energy conversions.

The recent spatial pattern of macropores is efficient in reducing free energy of soil water
Numerical experiment 4 corroborates furthermore that the spatial pattern of worm burrow densities which reflect the habitat preference of the main ecosystem engineers is more efficient in reducing free energy of soil water than the artificial organized pattern that has been flipped from hilltop to hillfoot. The typical soil catena with accumulated fine materials at the hillfoot sectors is a fingerprint of past dissipative processes. Higher worm burrow densities downslope reflect their habitat preference for moist but drained Colluvisols and thus ecological rather than thermodynamic optimality. Note that case 2 in experiment 4 assumes that worms/ecosystem engineers with a (different) habitat preference for drier soils dominate macropore formation. This hypothetical form of spatial organization has been shown to produce more overland flow, and more wetting in hilltop areas, which increases potential energy and total free energy of soil water within the 1.5 yr-long simulation. This implies that the "reorganized" system is on average further away from LTE and needs more time for relaxing back towards LTE. We thus conclude that the recent spatial pattern of macroporosity, which evolved in response to past mass and energy flows and related dissipative processes, is indeed superior in the sense of H2 when compared to the other spatially organized pattern. This finding clearly supports hypothesis H2. We conclude that the combination of soil types and vertical dead-ended macropores and their spatial arrangement along a superordinate potential gradient driving overland flow is not a random product but might be the result of "co-evolution"  towards a configuration that increase reduction rates of free energy of soil water in this environment.
But what is the benefit for the ecosystem engineers? As there is no direct feedback between the energy conversions due during the rainfall-runoff procedure and population dynamics of earthworms, worms have no direct advantage from accelerated power in soil water flows. Nevertheless, these macropores apparently accelerate reduction rates of free energy in soil water, and two apparent optima exist, resulting from the mathematical fact that the free energy reduction is a second-order polynomial. One can speculate that coevolution selects species whose ecological optimum (pattern and density) coincides with the thermodynamic optimum partitioning of rainfall water into overland flow into infiltration in such a landscape. In this sense and the sense these optimum configurations can be deemed as most probable states in landscape evolution (Dewar, 2009). Possible advantages for the worms are as follows: -The optimum is a metastable state, establishes stable environmental conditions and reduces energetic extremes and thus disturbance as fast as possible.
-It assures moist but not too-moist conditions and avoids dry stress. This probably prolongs the phases in the year were worms can become active. They are at rest during dry summer conditions (van Schaik et al., 2013).
-As erosion events do still rarely occur in the first optimum, it assures that fine soil material from upslope regions accumulates at the hillfoot during erosion events. This increases fertility in downslope areas, as these guys live from litter.

Distinguished macroporosity establishing in young volcanic ash soils
In the case of the Malalcahuello catchment we did not find a macroporosity factor that maximized averaged free energy reduction. This is explained by the fact that accelerated increases of potential energy of soil water strongly dominate over reduction of capillary binding energy during rainfalldriven conditions, when increasing the macroporosity factor. Free energy of soil water is thus on average increasing during rainfall-driven conditions, and the exchange rates of free energy scale linearly with the wetting rate. This implies at best the existence of a trivial optimum at the upper margin of the parameter space, which would then maximize drainage. Soil water potentials in the upper 1 m of the soil are due to the large rainfall amount on average +22 m away from equilibrium (Fig. 13). Even at a macroporosity factor of 70 (run 20), which corresponds to maximum flow velocities in the macroporous layer of almost 2.1 × 10 −2 m s −1 , drainage is not sufficient to reduce this excess and to shift the soil closer to LTE. Export of potential energy by means of subsurface storm flow events dominates reduction of free energy during rainfall-driven conditions, and increases with increasing macroporosity, which in fact boosts lateral soil water flows that feed subsurface storm flow. This finding thus does fully corroborate hypotheses H1. Although no optimum exists we defined a distinguished macroporosity, which brings the system closest to steady state in the sense that average gains in potential energy of soil water during rainfall-driven conditions are equal to average export rates of potential energy by subsurface storm flow and groundwater recharge. This distinguished macroporosity factor is, with a value of f mac = 50, close to the value of 54 that yielded the best fit and performed reasonably well during an uncalibrated simulation.
We thus conclude that in systems where potential energy dominates soil wetting during rainfall-driven conditions, infiltration events push the system further away from local thermodynamic equilibrium, because free energy of soil water is increasing. In this case at best a trivial optimum macroporosity at the upper margin of the parameter space exists. Nevertheless, an energy-centred analysis of the rainfall-runoff process can be very helpful to explain the role of macropores in the rainfall-runoff process. The macroporosity in the Malalcahuello catchment seems to assure steady-state conditions with respect to the potential energy balance during the rainfall-runoff process, which implies that the rainfall input on average does not push the system further away from LTE than it already is: potential energy gains are compensated by the potential energy losses. Note that a steady-state potential energy balance is not equal to a steady-state mass balance. As water leaves the system with finite potential energy, steady state is determined by the geo-potential of upslope storage regions and the amount of stored water, which is to be balanced with the geo-potential levels where water leaves the system and the amount of runoff. Based on this idea one can estimate the storage to runoff ratio based on the ratio of geo-potential levels where water leaves the system to the geo-potential of upslope storage regions. This estimate yielded a storage to runoff ratio of 0.26. When assuming that steady-state soil moisture storage equals steady-state evaporation, this yields an annual runoff coefficient of 0.74, which is in reasonable accordance with the value of 0.65 that Blume et al. (2008b) derived from observations.

Concluding summary and outlook
The presented results corroborate that an energy-centred analysis of the rainfall-runoff process provides valuable information for understanding the role of preferential flow processes in different hydro-climatic settings. The proposed framework establishes a useful link between the water balance of the critical zone and conversions, export and dissipation of free energy (potential energy, kinetic energy, capillary binding energy). Based on this framework we may shed light on the dynamics of free energy in the critical zone by analysing water balance components simulated with physically based hydrological models.
Free energy of water stored in the soil plays a key role in this concert, as it determines local thermodynamic equilibrium in soil, when soil water potential is zero along the entire profile. Dynamics is associated with deviations from the local thermodynamic equilibrium either to states of excess negative capillary binding energy in case of evaporative losses during radiation-driven conditions or to states of excess positive potential energy in the case of infiltrating surface water. We found that soil water dynamics can -depending on the interplay of the soil water retention properties, the topographic and climate setting -be classified into different thermodynamic regimes. The end members of these regimes are characterized by dominance of capillary binding energy, which requires fine porous cohesive soils and apparent longer dry spells, or by dominance of potential energy, which is favoured by non-cohesive coarse porous soils, a steep topography and a wet climate. In the former case free energy of soil water is decreasing during rainfall input, which implies relaxation towards LTE. In the latter case free energy of soil water is increasing, meaning that the system is pushed further away from LTE.
When capillary binding energy is dominating free energy dynamics, dead-ended macropores act as wetting structures during rainfall-runoff, thereby reducing relaxation times back to LTE, especially after long dry spells (Fig. 2). Several optima in macropore density might exist in this case, because reduction rates of free energy are a second-order polynomial of the wetting rate, which implicitly depends on macroporosity. In contrast to this, we do not expect an optimum density of macropores when potential energy dominates free energy of soil water, because change rates of potential energy of soil water scale linearly with the density of drainage structures. However, we defined a distinguished macroporosity which assures that potential energy of the system is on average in steady state during rainfall-driven conditions.
Based on the proposed framework we analysed numerical experiments that perturbed behavioural model structures, which have been shown to closely portray system behaviour and its architecture in the Weiherbach and the Malalcahuello catchments. Both headwaters are located in strongly different hydro-climactic and hydro-pedological settings and can be deemed to represent end members of the possible thermodynamic regimes of soil water dynamics (Fig. 13). In the Weiherbach catchment we found indeed two local optima in macroporosity which maximized average reduction rates in free energy during rainfall-driven conditions, when analysing the free energy balance for varying density of wetting structures. The macroporosity of the first optimum comes close to the one that yielded the best fit of the water balance within a 1.5 yr-long simulation. A simulation of the catchment scale water balance for the same period based on the optimal macroporosity factor performed only a little worse than the best fit. Note that this is an uncalibrated simulation, as observed discharge did not enter the analysis of the free energy balance of the hillslope. In the Malalcahuello catchment, as expected, we did not find an apparent optimum of drainage structures. However, the distinguished macroporosity was close to the value that yielded the best simulation of the water balance in the snow-free period, and a related estimate of the annual runoff coefficient was in reasonable accordance with the observed one.
We suggest that this new thermodynamic perspective allows for better discrimination of the different functions of preferential flow paths as either of the following: -Wetting structures when capillarity dominates. They enhance redistribution and storage and thus dissipation of free energy by depletion of gradients in soil water potential on the negative branch in Fig. 2. This might imply the existence of local optima in their density.
-Drainage structures when potential energy dominates. They enhance export of excess water and thus export of free energy and dissipation of free energy of soil water by depletion of soil water potential gradients on the positive branch in Fig. 2. This does not imply the existence of local optima, as reduction of potential energy in soil water scales linearly with the soil moisture change.
The presented findings are also promising for predictions in ungauged catchments (PUB; Sivapalan et al., 2003) as the thermodynamic optimal/distinguished model structure can serve as a first guess for an uncalibrated simulation of rainfall-runoff. The proposed thermodynamic treatment does require physically based models that account for dynamics in soil water potentials, as this determines free energy of soil water.
The proposed approach might also be useful for discriminating behavioural combinations of soils, topography, and climate as suggested by Schaefli et al. (2011) in their behavioural model study, which ultimately aims to link the Darwinian and Newtonian paradigms (Harman and Troch, 2013). Finally we suggest that Fig. 13 could be used to classify catchments with respect to thermodynamic regimes of soil water dynamics reflecting soil properties, the hydroclimatic forcing and the ecological setting: -Does capillarity or potential energy dominate (left or right sector of Figs. 2 and 13)? This depends on the soil water retention properties, which are a fingerprint of parent rock and past weathering processes, i.e. biotic and climatic processes, climatic setting and topography.
-The distance of the long-term average of soil water potential to zero should reflect the constraining factors that determine long-term average distance to local thermodynamic equilibrium. This could be sensitive to a dryness index and prevailing vegetation but also to soil hydraulic properties and soil structure.
-Deviations to negative values from equilibrium are controlled by transpiration, reflecting climate and ecological setting but also depth to groundwater.
-Deviations to positive values from equilibrium are controlled by rainfall regimes and the partitioning of rainfall into storage and runoff and thus the density and topology of wetting and drainage structures in soil.
In this respect it seems interesting to explore the full concert of processes that affect the free energy balance of the critical zone, including the full land surface energy balance and soil heat budget. This requires, however, physically based hydrological models that close both the mass and the energy balance in an explicit manner. This implies that one could add a thermal energy density axis to Fig. 13 to account for the stocks of thermal energy, which depends not only on the energy balance of the land surface but also on soil water regimes which strongly affect soil thermal properties. We suggest that the proposed view opens a new alternative energy-centred approach to classify catchments according to the similarity of their free energy balance components, because this integrates both the mass balance and land surface energy balance components, which jointly control catchment functioning, and judges their contributions in a single "currency" named Joule.