Transferring the concept of minimum energy dissipation from river networks to subsurface flow patterns

Principles of optimality provide an interesting alternative to modeling hydrological processes in detail on small scales and have received growing interest in the last years. Inspired by the more than 20 years old concept of minimum energy dissipation in river networks, we present a corresponding theory for subsurface flow in order to obtain a better understanding of preferential flow patterns in the subsurface. The concept describes flow patterns which are optimal in the sense of minimizing the total energy dissipation at a given recharge under the constraint of a given total porosity. Results are illustrated using two examples: two-dimensional flow towards a spring with a radial symmetric distribution of the porosity and dendritic flow patterns. The latter are found to be similar to river networks in their structure and, as a main result, the model predicts a power-law distribution of the spring discharges. In combination with two data sets from the Austrian Alps, this result is used for validating the model. Both data sets reveal power-law-distributed spring discharges with similar scaling exponents. These are, however, slightly larger than the exponent predicted by the model. As a further result, the distributions of the residence times strongly differ between homogeneous porous media and optimized flow patterns, while the mean residence times are similar in both cases.


Introduction
Preferential flow due to heterogeneity is one of the most important topics in subsurface hydrology.In many cases, large parts of the uncertainty in modeling subsurface flow and transport arise from preferential flow patterns which are either not known in detail or too small-scaled to be included explicitly in the model.
Preferential flow patterns may be the result of external pre-design, but may also be created by the flow itself.Fractured aquifers and reservoirs mainly fall into the first category.Fractures may be opened at high fluid pressures (hy-35 draulic fracturing, fluid-induced seismicity), but the basic structure of the preferential flow pattern is still governed by the geologic pre-design here.In contrast, the role of predesign is less clear for preferential flow in soils and in karstified aquifers.In particular, the formation of conduit pat-40 terns in karst has been addressed by several modeling studies (e.g.Groves and Howard, 1994;Howard and Groves, 1995;Siemers and Dreybrodt, 1998;Kaufmann andBraun, 1999, 2000;Liedl et al., 2003;Dreybrodt et al., 2005;Kaufmann et al., 2010;Gabrovšek and Dreybrodt, 2011;Hubinger and 45 Birk, 2011), pioneered by early work of Kiraly (1979).Although there may be some pre-design, too, it is generally believed that the solution of material by the flow causes a strong tendency towards self-organization of the flow pattern.
The idea to derive self-organizing heterogeneity from prin-50 ciples of optimality instead of small-scale, physically-based models has received growing interest in the last years.While McDonnell et al. (2007) considered optimization in hydrology as a visionary concept some years ago, the rapidly increasing number of publications related to several topics in 55 hydrology (e.g.Kleidon and Schymanski, 2008;Kleidon and Renner, 2013;Kleidon et al., 2013Kleidon et al., , 2014;;Zehe et al., 2010Zehe et al., , 2013;;Westhoff and Zehe, 2013) shows that it has left this level and even gained predictive power now.
The idea to explain the morphology of river networks from 60 the concept of minimum energy dissipation is the perhaps oldest application of optimization principles in hydrology receiving considerable interest.It even dates back to the early 1990s (Howard, 1990;Rodriguez-Iturbe et al., 1992a, b;Rinaldo et al., 1992) and turned out to be rather successful in 2 S. Hergarten, G. Winkler, and S. Birk: Minimum energy dissipation in subsurface flow reproducing scale-invariant properties of drainage networks at the Earth's surface.
The probably most successful application of minimum energy dissipation to flow processes so far (West et al., 1997) emerged from the field of physiology.It considers the cardiovascular system of mammals as a space-filling hierarchical network of tubes.The principle of minimum energy dissipation was shown to result in a scale-invariant structure of this network, which finally reproduces several allometric scaling relations found in nature.This seminal article was followed by several publications also in very prestigious journals (Enquist et al., 1998(Enquist et al., , 1999;;Banavar et al., 1999;West et al., 1999a, b) within a few years.
Beyond minimum energy dissipation, the minimum and maximum production of entropy are the most important vari-80 ational principles in nonequilibrium thermodynamics.As explained, e.g. by Martyushev (2013), both are not opposed to each other, but refer to different thermodynamic constraints where maximum entropy production has a wider range of applicability.Minimum energy dissipation can also be related to maximum entropy production (e.g.Županović et al., 2010).So it is the most general concept of nonequilibrium thermodynamics for application to coupled or spatially distributed systems in hydrology, and all the recent papers on optimization in hydrology mentioned above hinge on maxi-90 mum entropy production in some sense.However, maximum entropy production is a thermodynamic principle, and delineating systems or subsystems that indeed act at their thermodynamic limit and their thermodynamic constraints is nontrivial.As pointed out by Westhoff and Zehe (2013), finding 95 the properties that have to be maximized or minimized in a given system is therefore still a challenge.
In this paper, we transfer the concept of self-organization towards minimum energy dissipation from river networks to subsurface flow patterns and directly derive equations for the 100 optimal spatial distribution of porosity and permeability.In principle this is the realization of the vision on preferential subsurface flow patterns described by McDonnell et al. (2007) for the case of saturated flow.A brief review on the established theory for river networks is presented in the fol-105 lowing section, while the more complicated theory for subsurface flow is presented in Sect.3.

The Concept of Minimum Energy Expenditure in River Networks
Scale-invariant properties of river networks have been known 110 for a long time (Horton, 1945;Hack, 1957), and the idea to relate these properties to minimum energy dissipation (Howard, 1990;Rodriguez-Iturbe et al., 1992b, a;Rinaldo et al., 1992Rinaldo et al., , 1998) ) seems to be the first substantial application of optimality principles in hydrology.The basic idea is 115 that river networks organize in such a way that the total energy dissipation of the water on its way through the entire domain is minimized.So it refers to the surficial part of the water cycle starting where precipitation strikes the surface and ending when the water reaches the ocean (or any other 120 fixed reference point).
If changes in kinetic energy are neglected, energy dissipation can be computed without any knowledge of flow dynamics since the potential energy of the water is dissipated when it flows downslope in a channel.Then the mean energy dissi-125 pation (energy per time) of an individual channel segment is where ρ is the density, g the gravitational acceleration, q the mean discharge (volume per time), l the length, and is S the 130 slope of the segment.In order to find those channel networks with the lowest energy dissipation (called optimal channel networks or OCNs), the domain is subdivided into discrete cells exposed to a uniform precipitation.The centers of the cells are linked by channel segments in such a way that each 135 site (except for one or more outlet sites at the boundary) drains to exactly one neighbored site, but can be supplied by an arbitrary subset of its neighbors.This leads to dendritic river networks.The energy dissipation of the entire domain is readily obtained by summing up the contributions of all 140 cells (Eq.1): However, minimizing P without any constraints makes no sense as the trivial solution is a flat topography where all slopes vanish, so that there is no potential energy to be dissi-145 pated.Therefore, the theory of OCNs introduces a constraint to the minimization of P by assuming a relationship between slope and discharge, namely with θ ≈ 0.5.This relationship was originally found empiri-150 cally by Hack (1957) from the analysis of longitudinal river profiles, but substantiated theoretically later by the so-called stream-power approach (Howard, 1994).This approach predicts the erosion rate as a function of slope and discharge and is the basis of almost all large-scale models of fluvial erosion 155 in the detachment-limited regime.In this context Eq. ( 3) can be interpreted as an equilibrium of erosion and a spatially constant uplift rate in a homogeneous material.In most considerations (including Hack's original work and the work on OCNs), the catchment size is used instead of the discharge q, 160 which makes no difference in case of uniform precipitation (minus evaporation) in the entire domain.
The OCN concept constrains the minimization of energy dissipation to steady-state topographies where homogeneous uplift is balanced by erosion, i.e. following Eq.(3).Inserting 165 this relation into Eq.(2) leads to with γ = 1 − θ ≈ 0.5.This expression allows a minimization of the energy dissipation at a given uplift rate (which itself hides in the constant of proportionality) from the network topology alone without explicitly considering the topography, provided that the precipitation (minus evaporation) is given.In the studies on OCNs, drainage networks that minimize Eq. (4) (strictly speaking, the version without the term l i that takes into account the length of the river segments) were determined numerically.
It was found that these OCNs reproduce some well-known statistical properties of real river networks, such as the relationship between river length and catchment size and the statistical distribution of the catchment sizes, quite well.But on the other hand, later studies using simple erosion models (Hergarten and Neugebauer, 2001;Hergarten, 2002) revealed that some of the statistical properties can be reproduced even a little better for river networks exposed to permanently changing boundary conditions which are clearly above the minimum energy dissipation.In general, the considered statistical properties of drainage networks seem to be rather robust, and it seems to be impossible to tell how close real-world river networks are to the thermodynamic limit of minimum energy dissipation.
In the following section we transfer the concept of minimum energy dissipation to subsurface flow.Afterwards we present two applications to different scenarios and a validation using spring discharge distributions in order to demonstrate that the concept may be as powerful as the original minimum energy dissipation idea for river networks.

Theory
In the following, we assume steady-state Darcy flow of an incompressible fluid in an isotropic, but inhomogeneous porous medium with a hydraulic conductivity K(x).The volumetric flow rate is then given by q(x) = −K(x)∇h(x) (5) with the hydraulic potential h(x).As Darcy's law neglects effects of inertia, changes in kinetic energy shall be neglected.Therefore, the potential energy of the water is completely dissipated, and the total energy dissipation per time is where the integral extends over the considered domain, and q(x) = |q(x)|.
The key point of this paper is determining the distribution of K(x) which minimizes the total energy dissipation.
This can be either achieved by minimizing Eq. ( 7) with re-215 gard to K(x) if h(x) is given or by minimizing Eq. (8) if q(x) is given.In principle, even a simultaneous minimization of Eq. ( 7) with respect to both K(x) and h(x) can be performed.However, the part referring to h(x) only yields the mass balance which is shown in the following.If the spatial 220 distribution of the conductivity K(x) is given, the field h(x) which minimizes P is given by the Euler-Lagrange equation of Eq. ( 7): where the derivatives are functional derivatives.This imme- so that the distribution of h(x) which minimizes the energy dissipation automatically satisfies the mass balance for an incompressible fluid under saturated conditions as long as there 230 are neither sources nor sinks.Source and sink terms can be included by taking into account the potential energy of the water entering or leaving according to its hydraulic potential.This minimization is the basis of the variational approach in finite-element groundwater modeling, But here it just il-235 lustrates that it makes no difference whether we optimize K(x) and h(x) simultaneously or K(x) alone and consider Darcy's law with the mass balance separately.As the latter will be easier in the applications presented in Sect.4, we focus on the optimization of K(x) where either h(x) or q(x) 240 is given.
Similarly to the surface drainage network optimization, the optimization of K(x) has trivial solutions, namely P → 0 if K(x) → 0 everywhere (h(x) given, Eq. 7) and K(x) → ∞ everywhere (q(x) given, Eq. 8).We thus need a constraint 245 in analogy to the equilibrium of uplift and erosion for river networks.The idea that the total conductivity, i.e. the conductivity integrated over the entire domain, is given, may be straightforward at first sight.Alternatively, we can assume a given total pore space volume where φ(x) is the porosity.With respect to preferential flow patterns generated by the flow itself, e.g. in karst systems, the second idea is more suitable because the change in V is the total amount of solid volume that has to be removed (e.g.

dissolved).
In order to minimize the energy dissipation under the constraint defined in Eq. ( 11), the Euler-Lagrange formalism must be applied to the functional P −λV instead of P where the number λ is a Lagrange multiplicator.As neither P nor 260 V contain derivatives of φ(x), the result is formally the same as if φ was a parameter instead of a function: While the functional derivative at the right-hand size is almost trivial, computing the left-hand side requires a constitutive law for the conductivity as a function of the porosity, i.e. a function K(φ).
Relations between porosity and conductivity have been ex-270 tensively studied since the seminal work of Kozeny (1927) and Carman (1937).The original Kozeny-Carman equation predicts where the factor of proportionality includes a spatial length 275 scale and the denominator becomes important in case of very high porosities.As the length scale may change if conduits are widened, applying Eq. ( 14) without further considerations may be misleading.So let us first consider the simple case of parallel tubular conduits of given radii r i .Accord- ing to the Hagen-Poiseuille law, the total conductivity of this system is while 285 The relation between porosity and conductivity depends on the way how additional pore space is distributed among the conduits.Increasing all radii by the same factor β is the perhaps simplest concept.Then φ increases by a factor β 2 , while K increases by a factor β 4 , resulting in a relation K ∝ φ 2 .

290
However, spending all additional pore space volume for widening the largest conduit is most efficient with regard to the conductivity.In this case we obtain 295 which means that the slope in a double-logarithmic K-φ plot is two or larger.It approaches the value two obtained for widening all conduits by the same factor if either all conduits are equally sized or if the largest conduit is much larger than the rest.So the relationship between conductivity and porosity is not an overall power law, but behaves like a power law with n = 2 over some range, while it may increase faster in other regimes.Therefore, an overall power-law relationship with n ≥ 2 seems to be the best tradeoff between a good 305 approximation and keeping the following considerations as simple as possible.However, most of the results obtained in the following only require n > 1, which just means that the conductivity increases more rapidly than the porosity.For simplicity, we neglect the strong increase in conductivity due 310 to the coalescence of conduits for φ → 1 reflected in the denominator of Eq. ( 14).
Inserting the expression for the energy dissipation (Eq.7), the K-φ relationship (Eq.18), and Eq. ( 13) into Eq.( 12) yields The Lagrange multiplicator can be eliminated using the total pore space volume V (Eq.11), leading to and finally to So the optimal distribution of the porosity if h(x) is given is proportional to |∇h(x)| − 2 n−1 , while the rest of Eq. ( 21) is just a normalization to maintain the given total pore space volume V .

325
Minimizing the energy dissipation if q(x) is given is similar.Computing the functional derivative of Eq. ( 8) leads to the condition The Lagrange multiplicator can again be eliminated using the 330 total pore space volume according to λ = − ρgn a q(x) finally leading to This dependency is not only similar to that obtained for the 335 case that h(x) is given (Eq.21), but even basically the same.This becomes obvious if we transform Eq. ( 21) to conductivities, leading to Applying the same to Eq. ( 24) leads to Replacing q(x) with K(x)|∇h(x)| in Eq. ( 28), it is recognized that Eqs. ( 26) and ( 28) are equivalent.
Thus, the optimal distribution of porosities and conductivities with respect to minimum energy dissipation is the same for h(x) given and for q(x) given and can be summarized in the following relations: In principle, we can now insert this result into the mass balance of Darcy's law (Eq.10) and obtain the highly nonlinear differential equation However, we should keep in mind that this is a Darcy equation where the conductivity increases even more than linearly with the flow rate.This results in a strong tendency to focus flow, so that this equation will not have a unique, regular solution in general.For this reason we refrain from considering this equation in detail here.Beyond this, the resulting optimal distribution of K(x) can be used to represent the total energy dissipation as a function of either h(x) or q(x) alone.Plugging Eq. ( 25) into Eq.( 7) yields Thus, minimizing the total energy dissipation is equivalent to maximizing the functional It is easily recognized that Eq. ( 31) is the Euler-Lagrange equation of this functional.As a consequence, maximizing F (h) approximately among a set of functions obeying some regularity should also provide an approximate solution of the nonlinear Darcy flow problem posed in Eq. ( 31).
An equivalent functional with regard to q(x) is obtained by inserting Eq. ( 27) into Eq.( 8), resulting in Therefore, minimizing P is equivalent to minimizing the functional This functional will be used in the application presented in Sect.4.2.

Radial Flow towards a Spring
As a simple application we consider radial symmetric flow 385 towards a spring in two dimensions.Due to the twodimensional structure, K must be interpreted as a transmissivity, i.e. the product of conductivity and thickness, here.
An appropriate spatial distribution of the transmissivity has to be assumed when modeling the water fluxes in karstified 390 aquifers without having detailed information on the subsurface structure in the vicinity of springs.If the transmissivity was constant, the hydraulic gradients would strongly increase close to springs.While steep hydraulic gradients occur in the vicinity of artificial pumping wells, the hydraulic gradients 395 around karst springs are in general moderate.Therefore, the effective transmissivity (without explicitly considering fractures or karst pipes in the model) must increase towards the spring.The simplest assumption is an increase proportional to 1 r (where r is the distance from the point-like spring) that 400 just compensates the decrease of flow cross section towards the spring.This scenario was investigated by Birk and Hergarten (2012).
Let us now determine the optimal distribution of the transmissivity for the radial symmetric case.For a domain of ra-405 dius R with uniform recharge, the volumetric flow rate is where the dominator increases towards the spring due to increasing catchment size, and the denominator decreases due to decreasing cross section area.From Eqs. ( 29) and (30) we 410 immediately obtain so that 415 This leads to and thus close to the spring (r R) if h(0) = 0.For all values n > 420 1, the predicted increase in K(r) even overcompensates the decrease in cross section area towards the spring, so that the hydraulic gradient even tends to zero when approaching the spring.As discussed in Sect.3, n should even be at least two in natural aquifers, so that the optimal K(r) should increase

Dendritic Flow Patterns
We now come to the direct analog of the optimal channel network approach reviewed in Sect. 2. As stated at the end of Sect.3, minimizing the energy dissipation causes a tendency to focus flow, so that it makes sense to assume that each site in a discrete system only drains towards one of it neighbors.This topology is obvious for large-scale surface river networks except for situations where sediment deposition is important, e.g.braided rivers and river deltas where flow paths may indeed split up in downstream direction.The following argument shows that splitting up the flow of a site is energetically unfavorable for subsurface flow, too.Let us assume an arbitrary network of flow paths connecting discrete sites with any of their neighbors.For such a discrete network, the functional to be minimized (Eq.35) turns into that for river networks (Eq.4), Let us assume that the discharges q i describe the configuration with the minimum energy dissipation.If this configuration contains any site that drains towards more than one neighbor, there must be a set of configurations with discharges q i + δq i that is allowed within a range of around zero.The condition that the functional G(q + δq) has a minimum at = 0 requires but computing the second derivative immediately yields for 0 < γ < 1.So the original configuration with discharge towards more than one neighbor cannot be a minimum of the energy dissipation.As a consequence, non-dendritic maze cave structures are not optimal in the sense of energy dissipation.
We therefore consider only dendritic networks where each site has a unique flow direction towards one of it neighbors, in analogy to the theory of OCNs.Similarly to this approach, we focus on two-dimensional flow patterns in plan view.
We subdivide the domain into square cells where each cell is exposed to the same recharge and assume that each cell drains to one of its eight nearest and second-nearest neighbors.Darcy's law is only applied in one dimension along the dendritic network, while arbitrary flow directions not necessarily related to a hydraulic gradient are allowed.In analogy to the numerical determination of OCNs, the technique of simulated annealing based on the Metropolis algorithm (Metropolis et al., 1953) can be used for approaching the minimum energy dissipation iteratively.It starts from an 475 arbitrary initial configuration that only has to be consistent in the sense that a flow path from each site to the boundary exists.In each iteration step, the flow direction of a randomly selected site is changed to another (consistent) random direction.If the resulting change in energy dissipation is neg-480 ative, the change is accepted anyway.However, rejecting all changes where the energy dissipation increases temporarily results in getting stuck at local minima.Therefore, changes which seem to be bad at first sight must be allowed with some probability.During the iteration, the acceptance of changes 485 increasing the energy dissipation is slowly reduced, which can be interpreted as a slow thermodynamic cooling of the system.
Figure 1 shows four realizations on lattices of 4096×4096 sites where outflow is allowed across the entire boundary.As 490 expected, dendritic network structures occur for all values of n considered here.The tendency towards highly branched networks increases with n.Vice versa, a significant tendency towards the shortest way to the boundary is visible for n = 4 3 .However, this value is outside the range n ≥ 2 obtained for 495 reasonable relations between porosity and permeability, but will become relevant for the discussion on turbulent flow in Sect. 5.As mentioned above, n should not be much greater than two, while the typical value γ = 0.5 (Eq.4) for surface river networks corresponds to n = 3 according to Eq. ( 43).500 Therefore, our model predicts that the tendency to search the shortest way to the boundary should be slightly stronger for two-dimensional subsurface flow patterns than for rivers at the surface.However, it was already pointed out by Rinaldo et al. (1998) that such small variations in γ do not result in 505 significant differences in the statistical properties of the networks.
As little is known about the internal structure of real subsurface flow patterns, spring discharges provide the most easily accessible information.Figure 2 shows that the distribu-510 tion of the spring discharges Q is a power law described by the probability density with τ = 1.5.Transformed to a cumulative distribution, i.e. the number N (Q) of springs with a discharge greater than or 515 equal to Q, it says This result is insensitive to n at least in the range n ≥ 4 3 considered here, and thus for all reasonable relations between porosity and permeability as discussed above.

Validation
In the following we attempt to validate our model by two data sets of spring discharges from Austria.The first data set comprises the discharges of 1675 springs in the Niedere Tauern Range, and the second one the discharges of 1001 springs in the Semmering region.The methods of data acquisition range from simple estimates to direct capture and tracer dilution techniques.The frequency of measurement also varies from spring to spring.Despite the potential errors in the individual spring discharges, both data sets reveal clear powerlaw distributions of the discharges above a minimum value of about 0.1 liters per second (Fig. 3).This result suggests that the subsurface flow patterns in the two regions are indeed strongly organized.
Both power-law distributions are strikingly similar even quantitatively.However, the scaling exponent τ appears to be greater than the value τ = 1.5 predicted by the model (Eq.46).The estimate τ = 1.8 given in Fig. 3 was determined by visual correlation.Fitting a power-law distribution taking into account the statistical variation in the data reveals a rather high uncertainty in the scaling exponent as determining the minimum discharge where the power law begins is difficult.Starting a power-law fit at a minimum discharge of 0.1 l/s even yields τ < 1.5, but visually, the data still suggest τ > 1.5 in these two regions.
However, the occurrence of an organized flow pattern does not imply that it is indeed self-organized.And even if it is, it does not proof that this self-organization follows the principle suggested in this study.If the principle of minimizing energy dissipation holds here, the deviations in the scaling exponent may, e.g.arise from a three-dimensional flow organization or from the location of the springs in relation to the topography.On the other hand, the power law distribution might also be the result of a scale-invariant tectonic predesign and not be related to any self-organization in the flow pattern.Therefore, more datasets are needed to validate or refute our idea of subsurface flow organization.

Residence Time Distributions
Residence times and their statistical distributions are among the most important properties in subsurface hydrology.On one hand, they reveal valuable information about storage and flow pathways, and on the other hand they are essential for predicting the propagation of pollutants and recovery.
The spatial distribution of porosity and conductivity has a strong influence on the residence time distribution in a catchment.The residence time distribution is obtained by integrating 1 v(x) over the flow path from each point of the catchment to the spring.Here, denotes the flow velocity.Using Eq. 24 we obtain Thus, the flow velocity increases with increasing flow rate for all values n > 1.However, the increase is weaker than the linear increase occurring in a homogeneous aquifer.The 575 exponent in Eq. ( 50) amounts to one third for n = 2, which means that one third of an increase in flow rate is due to an increase in velocity, while two thirds arise from the increase of conductivity.At this point it should be mentioned that the residence time 580 distribution obtained this way describes only the contribution of the flow paths since it is assumed that the entire water in a representative volume moves exactly at the velocity defined by Eq. ( 48).Effects of dispersion and the properties of specific tracers (e.g.sorption) must be taken into account addi-585 tionally when making a prediction.The relationship between flow rate and flow velocity is the perhaps most important difference of our model of subsurface flow towards rivers at the surface.Although the theory reviewed in Sect. 2 does not account for flow velocities ex-590 plicitly, it is clear that the flow velocities of large rivers having a small slope are lower than those of small, but steeper rivers.Therefore, surface rivers show a negative correlation between discharge and velocity in general.
Figure 4 shows the residence time distributions of the ex-595 amples considered in Sects.4.1 and 4.2, obtained by evaluating the path integrals over 1 v(x) according to Eq. ( 49) numerically.Both geometries yield exponential residence time distributions for a homogeneous aquifer, but completely different distributions otherwise.In the example of radial flow, the 600 distribution changes from being strongly negatively skewed to positively skewed for increasing n.The respective distributions of the dendritic networks are always positively skewed and become more symmetric for increasing n.For n = 2 or not much larger which was found to be the most 605 reasonable range, the distributions are moderately positively skewed in both examples.
The residence time distribution for radial flow can be explained directly from Eq. ( 50) since the flow rate is entirely determined by the recharge and thus independent of n and 610 the same as for a homogeneous aquifer.For small values of n, the increase in velocity towards the spring is much weaker than for in the homogeneous case.As this affects the residence times in the entire catchment, it results in high residence times for large parts of the catchment.In return, the 615 velocities approach those of the homogeneous distribution for n → ∞, resulting in a positively skewed residence time distribution.
For the dendritic pattern, the residence time distribution does not only depend on the relationship between flow rate 620 and velocity (Eq.49), but also on the flow pattern.We there-fore have two competing effects of the dendritic structure compared to direct flow towards the boundary.First, the flow length becomes longer for (almost) all sites, but in return, flow in the preferential flow paths is faster.In our example we found that the mean residence times in the dendritic networks are by about 10 % higher than for a homogeneous conductivity.The longest mean residence time was found for n = 2, but the dependence on n is rather weak.Notably, the residence time distributions obtained for n ≥ 2 are similar to those arising from an advection-dispersion model, although the theory behind our model is completely different from this.

Extension towards Turbulent Flow
The theory presented in this paper hinges on laminar flow according to Darcy's law, while flow becomes turbulent at large flow rates in reality.This particularly applies to large conduits in karstic systems.In the following we give a sketch of an extension of our theory towards turbulent flow.
Let us first assume fully turbulent flow through parallel conduits of radii r i , so that the flow in each of them is described by the Darcy-Weisbach law with the flow velocity v i .Then the volumetric flow rate is which can also be written as If we increase all radii by the same factor β, the denominator increases by a factor β 5 , corresponding to φ 5 2 .In analogy to the laminar case (Eq.17) we can alternatively assume that the only largest conduit is widened and obtain We therefore assume with m ≥ 5 2 .Inserting this approach into Eq.( 6) yields Minimizing this expression under the constraint of a given total pore space volume (Eq.11) requires and thus 660 φ(x) ∝ q(x) This increase of the optimal porosity with the volumetric flow rate is even stronger than the increase like q 2 n+1 with n ≥ 2 found for laminar flow (Eq.24).Inserting the respective version of Eq. ( 24) with 3 m+1 instead of 2 n+1 into 665 Eq. ( 56) leads to so that the functional to be minimized is This functional is qualitatively the same as the one obtained 670 for laminar flow.The exponent 3 m+1 is still smaller than one for all values m ≥ 5 2 , but closer to one than the exponent 2 n+1 for laminar flow.The lowest value, m = 5 2 , is equivalent to the laminar case for n = 4 3 (which is unrealistic for laminar flow as we found n ≥ 2 there).Repeating the calculations 675 leading to Eq. ( 50) for turbulent flow reveals that the same result also holds for the residence time distributions.
Due to this equivalence, the upper left network shown in Fig. 1 (n = 4 3 , equivalent to m = 5 2 ) illustrates that fully turbulent flow conditions result in less dendritic flow patterns 680 than laminar flow conditions.But as shown in Fig. 2, the size distribution of the spring discharges is not affected.On the other hand, the residence time distribution (Fig. 4b) becomes strongly skewed in the turbulent regime.So some of the properties of the flow network are basically the same for laminar 685 and turbulent flow, while some other properties change significantly.
However, the most significant effect of turbulence should occur at the transition from laminar to turbulent flow.This transition is in general accompanied by a strong increase in 690 energy dissipation.For an individual pipe, the transition occurs when the Reynolds number exceeds a given threshold.The Reynolds number itself is proportional to the radius of the pipe and to the flow velocity.In our approach, the radii increase with the porosity and thus with the discharge.Ac-695 cording to Eq. ( 50), the flow velocity also increases with the discharge, so that the Reynolds number increases with the discharge, too.As a consequence, the transition to turbulent flow shall occur at a given critical flow rate.This also implies that the turbulent regime can only cover a part of the domain 700 in general.Therefore, the functional to be minimized has will become more complicated than the one given in Eq. (4): or for discrete networks where f (q) = 1 and γ(q) = 2 n+1 if q is below the turbulent threshold, while f (q) = const > 1 and γ(q) = 3 m+1 else.As this functional is even discontinuous, the minimization becomes nontrivial, and it is not even clear whether the conjecture that each site has a unique flow direction holds in the region close to the transition.

Limitations
Although our first results are promising, some limitations of the new theory are already visible.While the theory in itself should be consistent, these limitations concern the relevance of the thermodynamic limit of minimum energy dissipation for real subsurface flow patterns.
Even if we assume that the evolution of porosity goes in direction of minimizing energy dissipation, the question remains whether a self-organizing porosity distribution starting from a given state is really able to come close to the ground state of minimum energy dissipation.In karst systems where the evolution of porosity due to solution is understood quite well at small scales, this evolution is always directed towards higher porosity.Even if we assume that any increase in porosity is spent in such a way that the energy dissipation is reduced most efficiently, more or less parallel preferential flow paths may arise first even at rather small distances.Removing one of them in order to obtain a more dendritic and thus energetically better network may require a strong further increase in porosity which may be impossible in reality.So the spatial scale up to which the optimization can proceed may be limited in reality, and so far there is no quantitative relationship between this scale and the total (increase in) porosity.At first sight this seems to be a striking difference towards the evolution of river networks at the surface.Due to erosion and uplift permanently counteracting to each other, a river network may in principle have unlimited time to self-organize.However, as shown by Hergarten and Neugebauer (2001), the ability of river networks to erase existing patterns is also limited, so that the problem is basically the same as for subsurface flow.And as pointed out before, the considered statistical properties of both subsurface and surface flow networks are not very helpful in this context as they are rather insensitive to the distance from the ground state of minimum energy dissipation.
As a second major aspect, the examples presented in Sect. 4 are horizontal two-dimensional flow patterns and do not include a vertical component.For three-dimensional flow patterns, the recharge applied to the entire domain in two di-750 mensions must be replaced by an influx at the upper boundary.It is easily recognized that the optimal flow pattern then starts with horizontal flow along the upper boundary keeping the lower regions dry as this is the fastest way to focus flow.If the level where water leaves the domain at the sides 755 is lower, the preferential flow paths will approach this lower level only close to the sideward boundaries.In reality, there is a tendency towards vertical flow from the beginning, which means that the idea of minimum energy dissipation can only work in three dimensions if the anisotropy due to gravity is 760 incorporated in an appropriate way.Such an extension should include a reduced energy dissipation for vertical flow, but its quantitative representation is not yet clear.

Conclusions
We have presented a new theoretical concept to derive the 765 properties of subsurface flow patterns from a principle of optimization.In its spirit, the idea is similar to the established idea of minimizing energy dissipation for river networks at the Earth's surface.Due to different physical and empirical laws for surface and subsurface flow, our theory significantly 770 differs from that of surface river networks and is more complicated.Beyond the basic equations, we have to introduce different constraints for the minimization.While equilibrium between tectonic uplift and erosion is assumed for river networks, we assume a given total pore space volume and a re-775 lationship between porosity and permeability.
Despite the differences in the theories, both concepts finally arrive at basically the same functional to be minimized.As a consequence, the predicted preferential flow patterns in two dimensions are quite similar to surface river networks.

780
The predicted size distribution of the spring discharges follows a power law with a (non-cumulative) scaling exponent τ = 1.5.This result reproduces the spring size distributions of two regions in Austria quite well, although the scaling exponent seems to be slightly higher there, namely τ ≈ 1.8.

785
The other application shown in this paper, radial flow towards a spring, makes a suggestion how the transmissivities in the region around a spring might be chosen in a model where preferential flow patterns are not taken into account explicitly.

790
The residence time distributions seem to be the most striking difference between our model for subsurface flow and rivers at the surface.In contrast to rivers, our theory predicts an increase of flow velocity with flow rate, which is, however, not as strong as in an homogeneous aquifer.As a con-795 sequence of the weaker increase, mean residence times are in general higher for optimized flow patterns than for homogeneous aquifers.
Finally it should not be forgotten that, similarly to all other principles of optimization in hydrology, the question remains 800 whether real systems indeed organize towards their thermo-dynamic limit and, if so, how close they can come to the optimized state.The validation based on spring discharge distributions is promising, but not a proof that flow patterns indeed self-organize according to minimum energy dissipation.
So further studies are needed in order to investigate the applicability of our hypothesis.These studies will also have to include three-dimensional flow patterns and the influence of external pre-design.doi:10.1038/25977,1998.Enquist, B. J., West, G. B., Charnov, E. L., and Brown, J. H.: Allometric scaling of production and life-history variation in vascular plants, Nature, 401, 907-911, doi:10.1038/44819, 1999. Gabrovšek, F. andDreybrodt, W.: Spreading of tracer plumes 425 at least like r − 4 3 towards the spring, and h(r) should decrease at least like r 6 S. Hergarten, G. Winkler, and S. Birk: Minimum energy dissipation in subsurface flow