the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Flow recession behavior of preferential subsurface flow patterns with minimum energy dissipation
Jannick Strüven
Stefan Hergarten
Understanding the properties of preferential flow patterns is a major challenge in subsurface hydrology. Most of the theoretical approaches in this field stem from research on karst aquifers, where two or three distinct flow components with different timescales are typically considered. This study is based on a different concept: a continuous spatial variation in transmissivity and storativity over several orders of magnitude is assumed. The distribution and spatial pattern of these properties are derived from the concept of minimum energy dissipation. While the numerical simulation of such systems is challenging, it is found that a restriction to a dendritic flow pattern, similar to rivers at the surface, works well. It is also shown that spectral theory is useful for investigating the fundamental properties of such aquifers. As a main result, the longterm recession of the spring draining the aquifer during periods of drought becomes slower for large catchments. However, the dependence of the respective recession coefficient on catchment size is much weaker than for homogeneous aquifers. Concerning the shortterm behavior after an instantaneous recharge event, strong deviations from the exponential recession of a linear reservoir are observed. In particular, it takes a considerable time span until the spring discharge reaches its peak. The order of magnitude of this rise time is oneseventh of the characteristic time of the aquifer. Despite the strong deviations from the linear reservoir at short time spans, the exponential component typically contributes more than 80 % to the total discharge. This fraction is much higher than expected for karst aquifers and even exceeds the fraction predicted for homogeneous aquifers.
 Article
(5550 KB)  Fulltext XML
 BibTeX
 EndNote
The recession of spring discharge after recharge events can be seen as the fingerprint of an aquifer. In contrast to pumping tests at wells, it is a passive method based on data that are often recorded routinely. Additionally, spring discharges depend on the overall properties of the catchment, whereas pumping tests reflect the properties in a region around the well.
More than a century ago, Maillet (1905) proposed the linear reservoir, where discharge Q is directly proportional to the stored volume V:
The linear reservoir is described by a single parameter α (s^{−1}), and the stored volume follows the ordinary differential equation:
where R(t) (m^{3} s^{−1}) is the recharge. During periods with zero recharge, both the stored volume and the discharge decay exponentially with the same decay constant α, also called the recession coefficient:
The inverse of the recession coefficient, $\mathit{\tau}=\frac{\mathrm{1}}{\mathit{\alpha}}$, defines the efolding time: the time interval over which the discharge decreases by a factor of e.
The linear reservoir is not only appealing because it can be described by a single parameter but also because its behavior during periods of drought depends only on the actual amount of water and is independent of the recharge history. It often provides a reasonable approximation for long periods of drought. Deviations from the exponential decay at shorter timescales have been investigated and used for characterizing aquifers since the 1960s; in particular, karst systems have been addressed in numerous studies, and several theoretical concepts have been proposed.
Forkasiewicz and Paloc (1967) suggested a superposition of three distinct linear reservoirs with different decay constants describing three major flow components: a network of highly conductive conduits, an intermediate system of well integrated fissures, and a network of pores or narrow fissures with low permeability. The behavior of this model is dominated by the slowest reservoir during long periods of drought.
Mangin (1975) introduced a different approach using two components. The slow component was described as a linear reservoir, and a fast component with a limited range was added. The parameters of the fast components are the basis of the widely used karst classification system suggested in the same study. Several other modeling approaches which are similar in spirit were developed (e.g., Drogue, 1972; Atkinson, 1977; Padilla et al., 1994; Kovács and Perrochet, 2014; Xu et al., 2018; Basha, 2020; Kovács, 2021). Beyond these approaches, a multitude of numerical models designed for simulating realworld scenarios are now available. For deeper insights, readers are referred to the review paper by Fiorillo (2014) and to the model comparison by Jeannin et al. (2021).
While deviations from the linear reservoir are particularly relevant for karst systems, it should be noted that even the simplest Darcytype aquifers are not linear reservoirs. Assuming a given transmissivity T (m^{2} s^{−1}) and a given storativity S (–), the simplest Darcytype aquifer is described by the water balance equation:
where h (m) is the hydraulic head,
(m^{2} s^{−1}) is the 2D flux density (volume per time and crosssection width), ∇ is the 2D gradient operator, r is the recharge per area (ms^{−1}), and div is the 2D divergence operator. Inserting Eq. (6) into Eq. (5) yields a partial differential equation of the diffusion type for the hydraulic head h:
This model has been investigated in several studies for constant T and S in square or rectangular domains (e.g., Rorabaugh, 1964; Nutbrown, 1975; Kovács et al., 2005; Kovács and Perrochet, 2008). Applying spectral theory (see Sect. 2.4), it has been shown that the total discharge Q (q integrated over the entire boundary) can be described by an infinite series of exponential terms with different decay constants during periods of drought. As the longterm recession is dominated by the term with the smallest decay constant, the behavior approaches that of the linear reservoir. However, the question of how well the linear reservoir approximates the properties of real aquifers in a practical sense or whether the linear reservoir is even some kind of preferred state from a theoretical point of view has prompted debate (e.g., Fenicia et al., 2006; de Rooij, 2014; Kleidon and Savenije, 2017; Savenije, 2018).
Including the early recession phase in the analysis yields more information about the aquifer but also increases the dependence of the results on the recharge history. The instantaneous unit hydrograph, which dates back to concepts of Sherman (1932), is widely used in this context. It describes the discharge arising from a unit amount of recharge that is applied instantaneously at t=0 over the entire domain.
The unit hydrograph of the simple Darcy aquifer differs strongly from that of the linear reservoir as well as from the empirical approaches proposed by Forkasiewicz and Paloc (1967) and by Mangin (1975). While these models predict a finite discharge at t=0, the unit hydrograph diverges according to a power law,
in the limit t→0 for the simple Darcy aquifer (e.g., Hergarten and Birk, 2007). Such a powerlaw decrease also occurs in models consisting of porous blocks connected by highly conductive conduits (Kovács et al., 2005; Kovács and Perrochet, 2008). However, a finite conductance of the conduits limits the powerlaw divergence at short timescales (Kovács and Perrochet, 2014). Hergarten and Birk (2007) extended this concept using a fractal distribution of block sizes. While this model was able to explain a powerlaw recession with exponents different from $\frac{\mathrm{1}}{\mathrm{2}}$, deriving aquifer properties from the powerlaw behavior of recession curves turned out to be challenging. Birk and Hergarten (2010) investigated synthetic hydrographs for recharge events of finite duration and found that the properties of the recharge event likely obscure the shortterm dynamics of the porous blocks.
The quantification of heterogeneity and its representation in numerical models are still major challenges in hydrology. Carbonatic aquifers are particularly interesting in this context due to the interplay between fluid flow and structure. The spatial structure of the aquifer has a strong influence on fluid flow, which in turn controls dissolution and precipitation and, thus, the longterm change in the largescale structure. This interaction has been addressed in modeling studies of binary karst systems consisting of a porous medium and discrete conduits (e.g., Kaufmann and Braun, 2000; Kaufmann et al., 2010; Birk et al., 2003) as well as in the context of continuous changes in porosity and conductivity (e.g., Edery et al., 2021). Similar processes can also take place in soils by subsurface erosion (e.g., BernatekJakiel and Poesen, 2018).
As an alternative to forward modeling, there is an increasing number of studies attempting to derive preferred states of the atmosphere and the hydrosphere from principles of optimality (e.g., Kleidon and Schymanski, 2008; Kleidon and Renner, 2013; Kleidon et al., 2014; Kleidon and Savenije, 2017; Kleidon et al., 2019; Zehe et al., 2010, 2013, 2021; Westhoff and Zehe, 2013; Westhoff et al., 2014, 2017; Zhao et al., 2016; Schroers et al., 2022). In the context of fluid flow, minimum energy dissipation seems to be a promising concept, and it has been successfully applied to river networks (Howard, 1990; RodriguezIturbe et al., 1992a, b; Rinaldo et al., 1992; Maritan et al., 1996) and to the cardiovascular system of mammals (West et al., 1997; Enquist et al., 1998, 1999; Banavar et al., 1999; West et al., 1999a, b).
Hergarten et al. (2014) developed a theory for an optimal spatial distribution of porosity and hydraulic conductivity in the sense that the total energy dissipation of the flow is minimized. A major difference from the studies focusing on karst systems mentioned above is that this theory does not predict a binary or ternary system of distinct flow components but rather a continuous variation in the hydraulic properties over several orders of magnitude in combination with a highly organized spatial pattern. The predicted spatial pattern is dendritic and similar to drainage networks at the surface. Therefore, it does not entirely capture the variety of preferential flow patterns found in nature.
The potential relation of this theoretical concept to subsurface hydrology is still unclear. Validation is limited to the statistical distribution of catchment sizes and leaves several open questions (Hergarten et al., 2016). Even the properties of the model have not been analyzed thoroughly, except for residence times in a steady state (Hergarten et al., 2014).
This study investigates the dynamic properties of such preferential flow patterns based on the instantaneous unit hydrograph. The main research question is how preferential flow patterns with minimum energy dissipation are related to the different concepts discussed above. However, we also pose the following more detailed research questions:

How much does the behavior of preferential flow patterns differ from that of a homogeneous aquifer or a linear reservoir?

Can we explain the typical behavior of karst springs without assuming two or three distinct flow components?

Which physical properties do short and longterm properties of such an aquifer depend on?
2.1 Linearized consideration
The instantaneous unit hydrograph describes the response of an aquifer to adding a unit amount of water instantaneously. It is particularly useful for linear systems because the response to any recharge curve r(t) can be obtained by superposing multiple recharge events at different times (formally, the convolution integral of r(t) and the instantaneous unit hydrograph).
For unconfined aquifers, however, the transmissivity is proportional to the height of the water column:
where K is the hydraulic conductivity (m s^{−1}) and b is the base of the aquifer (m). As this dependence introduces a nonlinearity in Eq. (7), we can only make use of the linear theory for small disturbances. For this, we assume a steady state at a recharge r_{0} and apply a small additional recharge δr, so that $r={r}_{\mathrm{0}}+\mathit{\delta}r$. We can then write the actual head in the form $h={h}_{\mathrm{0}}+\mathit{\delta}h$ with the steadystate head h_{0}. Equation (5) retains its shape,
but with
for small δh. The first term is the effect of a change in the hydraulic gradient at constant transmissivity ${T}_{\mathrm{0}}=K({h}_{\mathrm{0}}b)$, which has the same form as Eq. (6). The second term describes the change in flux arising from the change in the thickness of the water column at constant hydraulic gradient. This term is particularly relevant for sloping aquifers with thin water layers.
As the concept developed in the following only captures the first term, the second term is neglected in this study. This means that the steadystate water table should be almost horizontal and the respective water column should not be too thin. In this case (neglecting the second term in Eq. 12), δh is described by the same equation as h in Eq. (7) but with δr and T=T_{0}. For simplicity, δh, δq, and δr are labeled h, q, and r, respectively. Thus, we have to keep in mind that these properties are not absolute values but rather small deviations from a steady state. Furthermore, T is the steadystate transmissivity T_{0} in the following.
2.2 Basic setup
Following the considerations outlined in the previous section, we consider a 2D aquifer with a given transmissivity T (m^{2} s^{−1}) and a given storativity S (–), described by Eqs. (5) and (6). As the focus is on strongly organized preferential flow patterns, both T and S are not spatially uniform and may even vary over several orders of magnitude. As discussed in the previous section, the head values h, the flux density q, and the recharge r are not absolute values but rather deviations from a steady state. As the simplest boundary condition, it can be assumed that the hydraulic head at springs remains constant, which is equivalent to the boundary condition h=0.
On a regular grid with a uniform grid spacing d, Eqs. (5) and (6) can be discretized by a finitevolume approach according to
with the fluxes
The nodes of the grid were numbered by a single index i, where N(i) denotes the nearest neighborhood of the node i, consisting of four neighbors except for boundary nodes. The symbol q_{ij} (m^{2} s^{−1}) refers to the flux from the node i to the node j, while T_{ij} is the respective transmissivity. Note that q_{ij} already includes the length of the edge between the two nodes (first term d in Eq. 14), so that it is no longer a flux per unit width. Inserting Eq. (14) into Eq. (13) yields the respective discrete form of Eq. (7).
Solving this equation numerically on large grids is challenging if the transmissivity varies over several orders of magnitude. Using an explicit scheme for the time step requires very small time increments because small changes in hydraulic head cause large fluxes. Employing a fully implicit scheme overcomes this limitation, as fully implicit schemes are stable at arbitrary time increments for diffusiontype equations. However, most of the available algorithms for solving the resulting linear equation system do not perform well if T varies strongly. This also applies to multigrid schemes (e.g., Hackbusch, 1985), which are the only schemes with a numerical effort that increases only linearly with the number of nodes.
2.3 Dendritic flow patterns
The theory of minimum energy dissipation in Darcy flow proposed by Hergarten et al. (2014) approximates preferential flow patterns by dendritic structures. This means that each node i delivers its entire discharge to one of its neighbors, b, referred to as the flow target in the following (strictly speaking, it should be labeled b_{i}). The flow target is defined as the neighbor with the steepest descent in hydraulic head, which is the same as the neighbor with the lowest head value for a grid with uniform spacing in both directions. Then, the neighborhood consists of three groups of nodes: (i) one flow target; (ii) some neighbors that deliver their discharge to the considered node, called donors in the following; and (iii) some nodes that do not interact with the considered node. The last group of nodes results in the difference from the original model where all neighboring nodes interact.
As we consider only small deviations from a steady state, we can assume that the topology of the flow pattern does not change through time. The flow target of each node is determined from the steady state and persists. As a major advantage, this simplification inhibits the occurrence of nodes without a flow target.
The discrete version of the balance equation (Eq. 13) becomes
for a dendritic flow pattern, where the notation q_{i} (with a single index) describes the flux from the node i to its flow target (so q_{ib} in the general notation). The sum extends over all donors of the node i, denoted here as D(i).
Because the topology of the flow pattern is static, the head value of a node may become lower than that of its flow target, in which case Eq. (14) predicts a negative flux. In the context of small deviations from a steady state, a negative flux would not have an immediate meaning because it just says that the flux towards the flow target is lower than in the steady state. However, there is no problem, even if backward flux occurs on an absolute scale.
Dendritic networks are widely used in the context of surface flow patterns, in particular river networks at large scales. In order to reduce the effects of anisotropy, the socalled D8 scheme is typically used on a regular, 2D grid. Here, eight neighbors (four nearest neighbors and four diagonal neighbors) are considered, so that river segments are either parallel to one of the coordinate axes or diagonal. While the numerical construction of optimized drainage patterns by Hergarten et al. (2014) also used the D8 scheme, we consider only the four nearest neighbors (D4 scheme) in the following. The main reason for this limitation is the comparison to the original model in Sect. 3.1. The D4 scheme can be seen as a restriction of the original model, whereas the D8 version would allow for additional flow paths. In particular, a diagonal line of points with a high transmissivity would be a preferential flow path with regard to the D8 scheme, although not in the original model.
While the concept of dendritic flow patterns was used by Hergarten et al. (2014) to construct patterns of porosity and conductivity (see also Sect. 2.7), recent developments in numerics make this concept interesting for timedependent modeling. In the following sections, two numerical approaches that are robust against strong variations in transmissivity are presented.
2.4 Spectral theory
Large parts of this study are based on spectral theory. Spectral theory decomposes the solution into a set of functions with a simple behavior. The functions are simple in the sense that they can be written in the form
which means that their shape does not change through time. For r=0, Eq. (7) can then be separated into one equation for h(x,0) and another equation for the time dependence f(t). This result is recognized by inserting h(x,t) into Eq. (7) and dividing both sides by h(x,t) and by S, which yields
As the lefthand side is independent of x and the righthand side is independent of t, both sides must be constant. If we call the respective constant −α, the solution for f is
because Eq. (16) requires f(0)=1. Thus,
decreases exponentially for α>0 and increases for α<0. The righthand side of Eq. (17) yields the eigenvalue equation
for h(x,0). This means that the differential operator $\frac{\mathrm{1}}{S}\mathrm{div}T\mathrm{\nabla}$ applied to the function h just multiplies h by the factor α. Mathematically, h is an eigenfunction of the differential operator and α is the respective eigenvalue.
While the eigenfunctions and the respective eigenvalues can be computed analytically for some simple geometries (e.g., Rorabaugh, 1964; Nutbrown, 1975; Kovács et al., 2005; Kovács and Perrochet, 2008) and constant parameters, a heterogeneous distribution of S and T requires a numerical treatment. The numerical treatment requires a transition from the continuous function h to the values h_{i} on the discrete grid. Using the same finitevolume discretization as above (Eqs. 13 and 14), the term $\frac{\mathrm{1}}{S}\mathrm{div}\left(T\mathrm{\nabla}h\right)$ can be approximated at the ith node by
where unit grid spacing (d=1) was assumed for simplicity. Aligning all values h_{i} in a column vector h, this relation can be written in matrix form
with a square matrix A. It is easily recognized that the nondiagonal elements of A are
whereas the diagonal elements are
This matrix is the same as that needed for solving the steadystate version of Eq. (7) numerically. The respective form of the matrix for dendritic flow patterns is basically the same, where only the neighborhood has to be reduced to those points that are connected, i.e., where either j is the flow target of i or vice versa.
As all software packages for linear algebra provide functions for computing eigenvalues, the steps are technically not very challenging up to this point. However, each eigenvalue and the respective eigenfunction is just one solution of Eq. (7) with r=0 for a specific initial condition h(x,0). In turn, we need the solution for any given initial condition. Spectral theory expresses a given initial condition as a superposition of eigenfunctions and then uses the knowledge about the individual eigenfunctions for predicting the full solution.
Let us assume that we computed all eigenfunctions and that e_{ki} is the value of the kth eigenfunction at the node i. If h_{i}(0) is the initial value at this node, we need coefficients λ_{k} so that h_{i}(0) is the superposition of the respective values e_{ki},
Then, the values of h at time t are
where α_{k} is the kth eigenvalue of A. Thus, all head values can be written as a sum of exponentially decaying terms, provided that all eigenvalues are positive. The latter can be shown for the matrix defined by Eqs. (23) and (24) with the help of Gershgorin's circle theorem.
Because all head values can be decomposed into a sum of exponential functions, the fluxes and, thus, the discharge of a spring can be expressed the same way:
where Q_{k} is the discharge of the kth eigenfunction. If we assume that the eigenvalues are sorted in increasing order, the longterm recession coefficient is α=α_{1}.
The nontrivial question is whether each initial condition can be expressed as a superposition of eigenfunctions according to Eq. (25). Mathematically, this property relies on the symmetry of the differential operator or, for the discrete problem, on the symmetry of the matrix (A_{ij}=A_{ji}). If A is symmetric, a basis of eigenvectors (the representations of the eigenfunctions on the discrete grid) can be found, which ensures the applicability of Eq. (25) for any initial condition with unique coefficients λ_{k}.
The symmetry of A even guarantees the existence of an orthonormal basis. As a main advantage, an orthonormal basis allows for the computation of each coefficient λ_{k} from the inner product (scalar product) of the initial head distribution and the respective eigenvector. Therefore, we can obtain the longterm recession coefficient α=α_{1} and the contribution of the respective component to the total (timeintegrated) discharge from the lowest eigenvalue and the respective eigenvector alone. Otherwise, we would need all eigenvectors. As their number equals the number of nodes, this would be expensive on large grids and would cost the advantage over a direct forward simulation.
Because T_{ij} is the transmissivity between the nodes i and j, T_{ij}=T_{ji}. Without the term S_{i} in Eq. (23), this property would already guarantee the symmetry of A. With this term, however, A is only symmetric if S_{i}=S_{j}, i.e., if S is the same everywhere.
To overcome this limitation, we need to get a bit deeper into linear algebra and consider a wider definition of symmetry coming from the concept of linear maps. Here, symmetry is related to an inner product, and the condition for symmetry is
where ⋅ is the inner product. This condition has to be satisfied for all h and $\stackrel{\mathrm{\u0303}}{\mathit{h}}$. It is easily recognized that the criterion A_{ij}=A_{ji} describes the specific case for the standard scalar product
For deeper insights into the fundamentals, readers are referred to textbooks of linear algebra (e.g., Halmos, 1958).
To examine the symmetry of A in this sense, we need the custom inner product
which differs from the standard scalar product (Eq. 29) by its weighting with the storativity S. Using the definition of A (Eqs. 23 and 24), it is easily recognized that
and thus
which proves the symmetry of A with regard to the custom inner product.
Technically, the normalization of the eigenvectors obtained numerically is the only point that requires attention. The eigenvectors have to be normalized with respect to the custom inner product and not with respect to the standard scalar product, i.e., according to the condition
If h_{i} represents the head values at t=0, the coefficients λ_{k} in Eqs. (25) and (26) are then obtained from the inner products in the following form:
This relation becomes particularly simple for the instantaneous unit hydrograph. Applying a unit amount of water per area to all nodes results in S_{i}h_{i}=1 for all i, and thus
Thus, the coefficient λ_{k} is obtained by adding the head values of the respective normalized eigenfunction at all nodes.
2.5 Forward modeling
Spectral theory is very efficient for characterizing the longterm recession, as only the components with the lowest recession coefficients α_{k} are relevant here. In turn, computing the instantaneous unit hydrograph at short timescales requires a large number of eigenfunctions, so that the advantage over direct forward modeling is lost. Therefore, we also implemented a numerical scheme for forward modeling of the dendritic model. For this purpose, we adopted the very efficient, fully implicit scheme that was recently proposed in the context of fluvial erosion and sediment transport by Hergarten (2020). As a major advantage, the scheme is robust against strong spatial variations in transmissivity. For simplicity, the scheme is only described for the D4 neighborhood with unit grid spacing (d=1) in the following.
Let us consider a time step from t to t+δt. Replacing the time derivative on the lefthand side of Eq. (15) with a difference quotient and evaluating the fluxes on the righthand side at the time t+δt (for the fully implicit treatment) yields
The key idea behind the scheme is that all fluxes respond linearly to changes in h of the flow target b due to the linearity of the differential equation. In particular, the flux q_{i}(t+δt) depends linearly on h_{b}, so that
Here, ${q}_{i}^{\mathrm{0}}$ is the hypothetic flux at time t+δt under the assumption that h_{b} remains constant (${h}_{b}(t+\mathit{\delta}t)={h}_{b}\left(t\right)$). It is not the same as q_{i}(t) because q_{i} may also change if q_{b} remains constant. The second term in Eq. (41) describes the effect of changes in h_{b}(t+δt) on q_{i}(t+δt), where ${q}_{i}^{\prime}$ is the partial derivative of q_{i}(t+δt) with respect to h_{b}(t+δt).
As shown in Appendix A, ${q}_{i}^{\mathrm{0}}$ and ${q}_{i}^{\prime}$ can be computed from the respective properties of the donors and from the known values of h at time t by the following expressions:
Thus, ${q}_{i}^{\mathrm{0}}$ and ${q}_{i}^{\prime}$ can be computed successively in downstream order, starting from the nodes without donors.
However, q_{i}(t+δt) cannot be computed during the downstream sweep because the actual values h_{b}(t+δt) required in Eq. (41) are still unknown. Computing the values h_{i}(t+δt) and q_{i}(t+δt) requires a second sweep over all nodes, which has to be performed in the opposite direction (upstream, starting from the boundary). As soon as h_{b}(t+δt) is known (as the node b is treated prior to the node i), q_{i}(t+δt) can be computed from Eq. (41). Finally, the hydraulic head can be calculated from Eq. (14) according to
This scheme provides a direct solver (without the need for iterations) for the fully implicit discretization of the flow equation on a dendritic network. As mentioned at the end of Sect. 2.2, fully implicit schemes are stable for arbitrary time increments δt. As a main advantage over other efficient schemes (e.g., multigrid schemes), the performance of this scheme is not affected by spatial variations in transmissivity and storativity. The numerical complexity is O(n), which means that the computing effort increases only linearly with the total number of nodes. It is, therefore, perfectly suited for simulations on grids with several million nodes.
2.6 Nondimensionalization
The linearized problem (Eq. 7, where T does not depend on h) can easily be transformed to nondimensional properties, which makes the results independent of the absolute values of S and T and of the spatial scale. For this purpose, arbitrary reference values S_{0} and T_{0} and an arbitrary length scale l can be defined. Then, the respective nondimensional properties
are defined. If we further introduce
and define the nondimensional properties
the differential equation for the nondimensional properties is the same as the original equation (Eq. 7). This result can easily be verified by expressing the original properties in Eq. (7) in terms of the respective nondimensional properties. The scaling factor for the flux density is T_{0} then, so $\mathit{q}={T}_{\mathrm{0}}\widehat{\mathit{q}}$.
As a consequence of these scaling properties, the shape of the instantaneous unit hydrograph depends only on the shape of the aquifer and on the spatial pattern of S and T. It is independent of the absolute values of S and T and on the absolute size of the aquifer. This property will become useful in the following section.
2.7 Patterns of transmissivity and storativity
While systematic knowledge about the spatial structure of preferential flow paths is still limited, we need a model for generating spatial patterns of transmissivity and storativity. In this study, we adopt the theory based on minimum energy dissipation proposed by Hergarten et al. (2014). This theory addresses the best spatial distribution of a given total pore space volume in the sense that the total energy dissipation of steadystate flow is minimized. The central assumption behind this theory is a powerlaw relation between the hydraulic conductivity K on the porosity ϕ,
with a given exponent n. Some ideas regarding how dissolution could increase K were discussed in the aforementioned study. If all pores have the same diameter and grow uniformly, an exponent n=2 was obtained. The same result was obtained for an arbitrary distribution of diameters if all pores grow by the same factor with respect to diameter. In turn, a stronger increase in K with ϕ (n>2) occurs if dissolution mainly affects the largest pores. Overall, these simple models yield a weaker increase in conductivity with porosity than predicted by the fundamental relation of Kozeny (1927) and Carman (1937). This relation predicts n=3 at small porosities, which was also adopted in the numerical study of dissolution and precipitation by Edery et al. (2021). In turn, the theoretical study of fractal pore spaces by Costa (2006) points towards n≈2.
The exact value of the exponent n is, however, not important in the following. Theoretically, it is only important that the increase in K with ϕ is stronger than linear (n>1) because concentrating porosity along preferential flow structures is only energetically favorable under this condition.
Based on the powerlaw relation, Hergarten et al. (2014) derived the following relations:
They also computed optimized dendritic networks on a discrete grid. This result is, however, limited to confined aquifers with a constant thickness. Assuming S=ϕ for an unconfined aquifer is the simplest and somehow most straightforward approach. In turn, the transfer from K to T is not trivial and depends on the topography of the base b. Even if the bottom is flat, points with the same q will have different values of h and, thus, different $T=K(hb)$, depending on their position in the network. As this dependency would make the theory of minimum energy dissipation more complicated, it is assumed in the following that the dependence of T on q is the same as that of K (Eq. 50). As a modification, we might take into account that sites with high q and, thus, high K may have a lower thickness h−b. The increase in T with ϕ will then be weaker than the increase in K, which could be included by reducing the exponent n. In each case, however, this 2D approach only captures horizontal heterogeneity. A high value of T requires that K is high over the entire thickness, so that neighboring sites with a high transmissivity are connected well.
As the most important aspect, however, Eqs. (49) and (50) refer to a steady state. The aquifer's properties do not change on the timescale of individual events. Instead, we assume that S and T are adjusted to a steady state that reflects some longterm average. At this point, we can make use of the nondimensional formulation developed in Sect. 2.6. Singlepixel catchments (nodes without donors) have the lowest values of S and T. Let us use their storativity and transmissivity as the reference values S_{0} and T_{0}, respectively. If we assume a uniform longterm average recharge, the mean flux is proportional to the catchment size A_{i} of the respective node. If we further use the grid spacing d as the characteristic length scale l, the respective nondimensional generalizations of Eqs. (49) and (50) are as follows:
Owing to the choice of l=d, the catchment size A_{i} has to be measured in grid pixels (so in units of d^{2}) here. Then, singlepixel catchments are characterized by ${S}_{i}={T}_{i}=\mathrm{1}$, corresponding to S_{i}=S_{0} and T_{i}=T_{0} in physical units.
The largest catchments in the examples considered later consist of slightly more than 10^{6} pixels, which means that the catchment size varies by a bit more than 6 orders of magnitude. For n=2, this corresponds to a variation in T by 8 orders of magnitude. If we assume, e.g., ${T}_{\mathrm{0}}={\mathrm{10}}^{\mathrm{9}}$ m^{2} s^{−1}, the values of T would cover the range from fresh limestone to strongly karstified limestone at a thickness of some meters. The respective range of S is only slightly more than 4 orders of magnitude for n=2. However, assuming ${S}_{\mathrm{0}}={\mathrm{10}}^{\mathrm{5}}$ would cover the range up to S=0.1.
Practically, the characteristic timescale t_{0} (Eq. 46) is the central property to be taken into account when transferring nondimensional hydrographs to realworld coordinates. With the values S_{0} and T_{0} defined above, this timescale would be t_{0}=10^{4} s at a grid spacing d=1 m and t_{0}=10^{6} s at d=10 m. As d seems to be a numerical parameter rather than a physical parameter, the strong dependence on d may be confusing at first. However, d is a property of the spatial pattern. It defines the smallest cell that cannot be subdivided by a preferential flow pattern and, thus, some kind of representative elementary volume (in combination with the thickness).
The characteristic timescale t_{0} can be interpreted physically in the setup used here. According to Eqs. (14) and (15), the smallest spatial unit (a singlepixel catchment) is described by the equation
for zero recharge. Therefore, the smallest spatial units behave like a linear reservoir with a recession coefficient
which means that t_{0} is the characteristic time of a singlepixel catchment.
2.8 Considered scenarios
Based on the assumptions described in the previous section, a numerically obtained flow pattern for n=2 on a 4096×4096 grid is used in this study. The algorithm for generating such patterns was described by Hergarten et al. (2014). All points at the boundaries are considered to be springs where the discharge is measured. Figure 1 illustrates the catchments of the springs.
Several scenarios will be considered for the same geometry in the following. Beside the reference scenario with n=2, the influence of n will be investigated. These investigations also include uniform transmissivity and storativity. The suitability of the description as a dendritic network will be tested against full Darcy flow.
In order to enable a comparison of the recession curves of individual springs, the same catchments are used in all simulations. This means that the boundaries between the catchments shown in Fig. 1 are enforced by cutting the respective connections, i.e., by inhibiting flow across the watersheds.
3.1 Dendritic flow patterns vs. full Darcy flow
As a first step, we investigate under which conditions the reduction from full Darcy flow (taking all four neighbors into account) to a dendritic flow pattern (a single flow target for each node) provides a suitable approximation. We start from the optimized distribution of S and T for n=2, i.e., from a strongly preferential flow pattern. Figure 2 compares the longterm recession coefficients α (α_{1} in Eq. 27) obtained from the dendritic flow pattern to those of full Darcy flow. The recession coefficients of both approaches agree well, which means that the longterm recession behavior is captured well by the dendritic flow pattern. Some deviations occur at rather small catchment sizes from about 10 to 200. Here, allowing only one single flow direction leads to a slight underestimation of α, which means that the recession is slightly too slow. This underestimation is highest at catchment sizes A≈16 and reaches about 10 % there.
In order to investigate the effect of the restricted flow directions in more detail, we analyzed the flow pattern of the full Darcy scenario. For this purpose, we define the major flux component of each node as the flux towards the given flow target (according to the dendritic pattern) and the minor flux component as the sum of the fluxes towards all other neighbors with lower head values.
Figure 3 shows the relative contribution of the minor fluxes as a function of the size of the respective upstream catchment. The upstream catchment size refers to the individual pixels here and not to the embedding catchment. Thus, the data point for a A=1 describes the average over all pixels without donors, regardless of whether they are draining directly to the boundary (i.e., are indeed singlepixel catchments) or are part of a larger catchment. As points at drainage divides are already restricted concerning their flow direction in the full Darcy scenario, these points are analyzed separately from inner points.
The contribution of the minor fluxes decreases with increasing catchment size, i.e., downstream along the preferential flow paths. It becomes negligible for all considered scenarios at catchment sizes A⪆200. The contribution of the minor flux is particularly small immediately after a short, uniform recharge pulse. If all sites receive the same amount of water, the resulting change in hydraulic head is inversely proportional to the storativity. Because high transmissivity goes along with high storativity, the hydraulic gradient between sites with low transmissivity and sites with high transmissivity increases. As a consequence, the fluxes are concentrated towards the preferential flow paths, which reduces the minor fluxes.
The strongest contribution of the minor fluxes is found for the longterm recession, characterized by an exponential recession curve. The minor fluxes contribute even almost 60 % for A=1 here (only inner points, about 35 % for drainage divides). It may be surprising that the minor fluxes do not affect the recession coefficient strongly, as shown in Fig. 2. In an extreme scenario where all four neighbors are at the same head values, the fluxes would be 4 times higher in the full Darcy model than in the singleflowtarget realization, which would result in a 4 times faster recession than predicted by Eq. (54). Although the majority of the nodes have small upstream catchment sizes (e.g., A=1 for 43 % of all nodes and A≤10 for 81 % of all nodes), their behavior is obviously not crucial for the properties of the entire catchment.
These results suggest that preferential flow patterns can be represented well on a discrete grid by a dendritic structure where each node drains only towards one of its neighbors. In turn, the approximation using a dendritic flow pattern does not work well for a spatially uniform distribution of T and S, as shown in Fig. 4. Here, α is underestimated by more than 1 order of magnitude for large catchments, which means that the recession is more than a factor of 10 too slow. The coefficients agree well only for small catchments, where the minor fluxes are small or even absent due to the restricted flow directions at drainage divides.
3.2 Scaling properties of the recession coefficient
The fact that the recession coefficient α decreases with increasing catchment size is already visible in Figs. 2 and 4. The expected scaling behavior is $\mathit{\alpha}\propto {A}^{\mathrm{1}}$. This result is formally obtained from the characteristic timescale t_{0} (Eq. 46 with A∝l^{2}), but it can also be derived directly from the occurrence of a firstorder time derivative on the lefthand side of Eq. (7) and secondorder spatial derivatives on the righthand side of Eq. (7). If we rescale the entire catchment including the pattern of S and T by a factor β, the righthand side of Eq. (7) changes by a factor of β^{−2} for r=0. Then, the timescale must change by the same factor. As the catchment size A increases quadratically with β, the timescale is proportional to the catchment size and, thus, $\mathit{\alpha}\propto {A}^{\mathrm{1}}$.
As shown in Fig. 5, this simple scaling behavior does not hold for the patterns of S and T considered here. While the scaling behavior follows a power law
reasonably well for all considered values of n, it seems that an exponent γ=1 is only achieved for n=1. For all values n>1, we find γ<1, where γ decreases with increasing n. This means that the recession of large catchments is still slower than the recession of small catchments, but the effect is considerably weaker than for a simple Darcytype aquifer (γ=1).
Qualitatively, this weaker increase is the expected behavior for a preferential flow pattern. Preferential flow paths are able to transport water rapidly over large distances, so that an increasing spatial extension does not slow down the recession as strongly as in a homogeneous aquifer.
The opposite behavior is observed for n<1. Here, large catchments become extremely slow. Revisiting Eqs. (51) and (52), it is recognized that T still increases with catchment size, but it is weaker than S for n<1. Therefore, flow is still facilitated along the preferential flow paths, but the respective cells become slow due to their high storativity. This becomes visible if we apply Eq. (54) formally to such cells, which yields $\mathit{\alpha}\propto {A}^{\frac{\mathrm{2}(n\mathrm{1})}{n+\mathrm{1}}}$, so that α decreases with increasing A. However, n<1 would be unrealistic for porous media, where n≥2 should instead be expected. Beyond this, Hergarten et al. (2014) demonstrated that dendritic flow patterns are energetically favorable only for n>1. Thus, preferential flow patterns with n<1 can be constructed, but they do not make much sense. Nevertheless, the curves are almost indistinguishable not only for n≥1 but also for n<1. The results reveal that the numerical approximation using a dendritic flow pattern also works well for n<1.
As already recognized in Sect. 3.1, the approximation using a dendritic flow pattern does not work for constant transmissivity and storativity. Figure 5 reveals that the unusual scaling behavior with γ>1 also occurs for the description using full Darcy flow. This result is related to the subdivision of the domain into fixed catchments drained by distinct springs. In order to test this hypothesis, we computed the recession coefficients for radial flow towards a spring in polar coordinates. Keeping the radius of the spring constant and varying the total size of the catchment (the outer radius), we obtained roughly the same scaling behavior, γ≈1.1, found for the catchments considered in Fig. 5. This scaling behavior indicates that the region around the spring is some kind of bottleneck that becomes increasingly relevant for large catchments. This result aligns well with the result of Hergarten et al. (2014), who found that minimum energy dissipation for radial flow requires an increase in permeability towards the spring.
The nontrivial dependence of α on the catchment size not only affects spring hydrographs but also the response time of aquifers to changes in climate or hydraulic conditions at the boundaries. As a recent example, Cuthbert et al. (2019) provided worldwide estimates of these groundwater response times. These estimates are based on the 1D version of Eq. (7) for homogeneous aquifers, with the size of the aquifer assumed to be twice the distance L of perennial rivers. The scaling arguments given above (and formally in Sect. 2.6) then yield $\mathit{\alpha}\propto {L}^{\mathrm{2}}$ and, thus, a response time proportional to L^{2}. Because L increases with decreasing precipitation, this model predicts a strong negative correlation between precipitation and groundwater response time with typical response times of several thousand years in arid and semiarid regions. If the respective aquifers are not homogeneous but rather have a preferential flow structure as considered here, the dependence of the groundwater response time on L and, thus, on climate will be considerably weaker than predicted by Cuthbert et al. (2019).
3.3 Shortterm recession
The differences between the model considered here and simple Darcytype models are not limited to the scaling properties of the longterm recession coefficient. As exemplified in Fig. 6 for the biggest catchment, the unit hydrograph shows a clear rising limb at short timescales. In contrast, the 1D Darcytype aquifer shows a powerlaw decrease in discharge at short timescales, so that Q→∞ for t→0.
For a better comparison, size and parameter values of the 1D aquifer were adjusted in such a way that the longterm recession coefficient α and the total amount of water supplied by the recharge event are the same as for the considered catchment (see Appendix B). For such a simple aquifer, a lag between a short precipitation event and the peak discharge would typically be attributed to the infiltration process. We would then assume that the time lag is related to the transit time of the water from the surface to the aquifer. In a model consisting of individual porous blocks connected by highly conductive fractures (e.g., Kovács et al., 2005; Hergarten and Birk, 2007), the time lag may also be owing to the finite conductivity of the fracture system. In contrast, the finite rise time is an inherent property of the structure of the aquifer in the model considered here. It cannot be attributed uniquely to any individual component of the system.
Figure 7 provides an analysis of the rise time t_{rise} for the considered catchments and different values of the exponent n. For n=2, the data suggest a linear relationship between t_{rise} and the longterm efolding time $\mathit{\tau}=\frac{\mathrm{1}}{\mathit{\alpha}}$. The ratio of t_{rise} and τ varies from about 9 % to 19 % for the individual catchments. The rise times become slightly shorter in relation to the efolding time for lower n. For n=1.5, we found ratios in a range from about 4 % to 22 %. The data obtained for n>2 are nonunique, in particular for small catchments (small τ). Overall, there seems to be a weak increase in the ratio of t_{rise} and τ with increasing n. This trend aligns well with the absence of a finite rise time in the simple 1D and 2D aquifers without a preferential flow pattern. However, the trend is much weaker than the variation among individual catchments; therefore, it is not investigated further. In each case, however, t_{rise} is not very small compared with τ. Oneseventh of the efolding recession time seems to be a reasonable order of magnitude. Taking into account that longterm efolding times of karst aquifers are typically in an order of magnitude of several weeks, the rise time is of an order of magnitude of several days.
The occurrence of a rising limb in the unit hydrograph requires that some of the coefficients λ_{k}Q_{k} in Eq. (27) are negative. This can be seen formally by computing the derivative of Eq. (27). The derivative is similar to Eq. (27) itself, except that the coefficients are −α_{k}λ_{k}Q_{k}. As the sum cannot become zero if all coefficients have the same sign, at least one of the original coefficients λ_{k}Q_{k} must be negative. Figure 8 shows the coefficients for the largest catchment. While the coefficients of the slowest components (small α_{k}) are positive, there is no obvious preference for either sign in the faster components.
The spectrum of the considered catchment differs fundamentally from that of the 1D aquifer. While the spectrum of the 1D aquifer consists of distinct components with recession coefficients ${\mathit{\alpha}}_{k}\propto (\mathrm{2}k\mathrm{1}{)}^{\mathrm{2}}$ (see Appendix B), the spectrum of the catchment becomes more or less continuous at large k. However, the smallest recession coefficients are still distinct. For the catchment analyzed in Fig. 8, the ratio is $\frac{{\mathit{\alpha}}_{\mathrm{2}}}{{\mathit{\alpha}}_{\mathrm{1}}}\approx \mathrm{2}$. While the lowest ratios among all simulated catchments are about 1.5, the ratio is typically in a range from about 2.5 to 4 for smaller catchments. Thus, the ratio $\frac{{\mathit{\alpha}}_{\mathrm{2}}}{{\mathit{\alpha}}_{\mathrm{1}}}$ differs from catchment to catchment but is always clearly above unity. This also holds for n≠2. This property ensures that the recession curve approaches a single exponential function at a reasonable time and that the coefficient λ_{1}Q_{1} captures the longterm recession well. However, the ratio is much lower than the ratio $\frac{{\mathit{\alpha}}_{\mathrm{2}}}{{\mathit{\alpha}}_{\mathrm{1}}}=\mathrm{9}$ of the 1D aquifer. Although the deviation from an exponential recession also depends on the coefficients λ_{k}Q_{k}, the difference is already visible in Fig. 6. While the unit hydrograph of the 1D aquifer is almost indistinguishable from the exponential decay at t=0.5τ, the difference is more than 15 % at this time for the largest catchment. Therefore, the model with the preferential flow pattern approaches the longterm exponential recession much more slowly than the simple 1D aquifer.
Figure 9 illustrates the approximation of the unit hydrograph by a finite number of exponential components. In this example, all coefficients for k≤9 are positive. In the range $\mathrm{9}<k\le \mathrm{23}$, there are also negative coefficients. However, the positive coefficients still dominate here, so that the highest peak discharge is achieved when including the components with k≤23. The negative components become increasingly relevant for larger k. This results in a decreasing peak discharge, however, that is still at t=0 for k≤60. The component with k=61 is the first that shifts the peak discharge to a time t>0 and, thus, produces a rising limb in the unit hydrograph. However, approximating the behavior around the peak discharge reasonably well requires several hundred components.
The occurrence of negative coefficients in Eq. (27) is not a fundamental problem, but it impedes a simple interpretation of the decomposition into exponentially decreasing terms. If all components were positive, we could imagine the aquifer as a set of linear reservoirs drained in parallel. However, negative components would correspond to reservoirs with a negative amount of water. Formally, the slow component (first exponential function) could even contain more water than available in total, and the fast components (all higher exponentials) could be negative in total. In Fig. 6, this would mean that the area of the unit hydrograph below the exponential component at short timescales was larger than the area above the exponential component. In this case, the widely used characterization of karst aquifers according to the contribution of the slow (exponential) component to the total discharge (Mangin, 1975; Jeannin and Sauter, 1998) would be taken ad absurdum.
Figure 10 shows the contribution of the slow component (first exponential) to the total amount of water for all catchments considered here (n=2). Depending on the catchment size, this contribution is between 77 % and 104 %. Thus, there are indeed catchments where the effect of the rising limb is so strong that the slow component is formally larger than the total recharge and the fast components are negative in sum. However, this effect is found only for some rather small catchments. For the largest catchments, the contribution of the slow component is about 90 %.
The obtained contributions of the slow component are high compared with other models. For the 1D aquifer used here as a reference, this contribution is $\frac{\mathrm{8}}{{\mathit{\pi}}^{\mathrm{2}}}\approx \mathrm{81}$ %. For an aquifer consisting of square porous blocks connected by conduits with infinite conductance, it is $\frac{\mathrm{64}}{{\mathit{\pi}}^{\mathrm{4}}}\approx \mathrm{66}$ % (e.g., Birk and Hergarten, 2010). The contribution of the slow component was not explicitly investigated by Hergarten and Birk (2007) in their fractal model with powerlawdistributed block sizes. However, as the slowest flow component arises from the largest blocks, the total contribution of the slow component must even be considerably smaller than the 66 % obtained for an aquifer with a uniform block size. In the widely used classification scheme of karst aquifers proposed by Mangin (1975), even aquifers with a contribution of the slow component of less than 50 % are considered poorly karstified (see also Jeannin and Sauter, 1998). Therefore, our continuous model of preferential flow patterns predicts an even higher contribution of the slow component to the unit hydrograph than other models and is very different from what is typically assumed for karst systems.
At first sight, the observed high contribution of the slow component could be an effect of the finite rise time. As recognized in Figs. 6 and 9, the unit hydrograph is below the slow component for t⪅0.5t_{rise}. When investigating recession curves in reality, the analysis typically starts from the peak in discharge (t=t_{rise}). The orange markers in Fig. 10 show the respective contributions of the slow component. Starting the analysis from t=t_{rise} instead of t=0 indeed yields lower contributions of the slow components. However, the effect is only of the order of magnitude of a few percent. Thus, the result that continuous preferential flow patterns predict a high contribution of the slowest flow component is not an artifact of the analysis.
3.4 The effect of nonuniform recharge
While the unit hydrograph describes uniform recharge, the spatial distribution of the recharge may have a strong influence on recession curves of large catchments. Therefore, a question arises regarding whether the strong deviations from the exponential behavior found for karst springs could be related to the occurrence of local rainstorms that affect only a small part of the catchment.
As a simple example, we separated the domain into a proximal part and a distal part. Both parts are equally sized, and a distinction is made based on the distance from the boundary of the domain. Because the overall domain is the same, the recession coefficients α_{k} of all flow components are the same as for the entire domain. Only the coefficients a_{k} in Eq. (27) differ. This difference, however, has a strong effect on the rise time and on the contribution of the slow flow component.
Figure 11 compares the results obtained numerically for the largest catchment to those obtained by spectral decomposition for a 1D aquifer (for details, see Appendix B). While the hydrograph starts with a peak at t=0 after a spatially uniform recharge event, applying recharge only to the distal part of the domain introduces a finite rise time. This rise time is, however, shorter (in relation to the efolding recession time) than for the aquifer with the preferential flow pattern. The rise time of the preferential flow pattern also changes considerably if recharge is only applied to a part of the domain. The strong influence may be surprising at first because we could expect the signal to propagate rapidly through the preferential flow structure from the distal part of the domain to the spring. However, preferential flow paths not only have a high transmissivity here but also a high storativity. Thus, the propagation of signals is not as fast as we might expect.
The contribution of the slowest component also changes if recharge is only applied to a part of the domain. Similarly to the rise time, it increases if only the distal part of the domain is filled because the instantaneous recharge signal has already been smoothed when it arrives at the spring. The contribution of the first exponential component is formally higher than 100 % for the distal part. For the largest simulated catchment, it is 137 %, whereas it is 115 % for the 1D aquifer. In turn, the contribution of the first exponential component is lower for the proximal parts: 50 % for the largest catchment in the simulation and 47 % for the 1D aquifer. Therefore, the contribution of the first exponential component is always higher for the preferential flow pattern than for the 1D aquifer, and the effect of only applying recharge to a part of the domain is similar.
On average, however, the effect vanishes. For the largest catchment, the 50 % for the proximal region and the 137 % for the distal region yield a mean value of 94 %, in agreement with Fig. 10 (the rightmost blue circle). This is a general property of the linear model, which allows for the superposition of recharge events not only temporally but also spatially. It could even be generalized to arbitrary parts of the domain down to recharge events that are limited to a single node.
While there should theoretically be no effect on average, recharge events in the proximal region will be more prominent in the hydrograph than those in the distal region. Therefore, the analysis of realworld hydrographs might be biased towards the more prominent, proximal recharge events, which would yield a decreasing contribution of the slowest component for large catchments. However, explaining the often small contribution of the exponential component this way seems quite a stretch.
This study is a first attempt to describe the dynamics of aquifers with continuous preferential flow patterns. In contrast to approaches based on two or three distinct flow components widely used in the context of karst aquifers, the concept used here assumes a continuous spatial variation in hydraulic properties over several orders of magnitude.
A 2D aquifer with a spatially variable but timeindependent transmissivity was considered. This scenario corresponds to the application of small disturbances to a steady state of an aquifer with an almost horizontal water table. Synthetic spatial patterns of transmissivity and storativity were obtained from principles of minimum energy dissipation based on the theory proposed by Hergarten et al. (2014).
As a major technical result, it was found that such aquifers can be approximated well by dendritic flow patterns, in which the entire discharge of each cell is delivered to the neighbor with the steepest gradient in hydraulic head. This approximation has been widely used for channelized flow patterns at the surface. The dendritic structure enables an efficient, fully implicit numerical scheme with a numerical effort that increases only linearly with the number of cells, also known as O(n) complexity. This property allows for simulations on grids consisting of several million nodes and, thus, for a reasonable spatial resolution of the preferential flow pattern.
As a second, rather theoretical result, it was shown that spectral theory is not restricted to homogeneous aquifers and can also be applied to aquifers with any spatial distribution of transmissivity and storativity. Although the eigenvalues and the respective eigenvectors have to be computed numerically, this approach allows for a fast computation of the longterm recession coefficient without forward modeling over a long time span. In addition, the contribution of the slowest flow component to the instantaneous unit hydrograph (and also to any other initial state) can be computed easily. However, the efficient numerical scheme and spectral theory rely on the assumption of timeindependent transmissivity and cannot be extended easily towards sloping aquifers.
The longterm recession coefficient α depends on the catchment size. The dependency is, however, weaker than for homogeneous aquifers and follows a power law $\mathit{\alpha}\propto {A}^{\mathit{\gamma}}$ (Eq. 55). The exponent γ depends on the assumed relation between transmissivity T and storativity S. It approaches 1 for T∝S, which is also the limit where dendritic flow patterns are energetically favorable. In this case, the scaling is the same as for homogeneous aquifers (γ=1). For relations T∝S^{n}, γ decreases with increasing n. As a typical value, γ=0.4 was found for n=2. Therefore, the decrease in the recession coefficient with catchment size is typically less than half as strong as for homogeneous aquifers. This finding challenges previous results on very long groundwater response times of large aquifers (e.g., Cuthbert et al., 2019).
Because flow patterns obtained from minimum energy dissipation are typically scaleinvariant, a powerlaw decrease in the discharge after a short recharge pulse might be expected at first. However, the respective instantaneous unit hydrograph shows a completely different behavior. The discharge immediately after the recharge event is quite small, and it takes a considerable time until it reaches its peak. The order of magnitude of this rise time is oneseventh of the characteristic time of the aquifer ($\mathit{\tau}={\mathit{\alpha}}^{\mathrm{1}}$). It seems not to depend strongly on the catchment size nor on the relation between transmissivity and storativity.
The contribution of the slowest component to the unit hydrograph is of an order of magnitude of 90 % for large catchments and even larger for small catchments. This contribution increases further if recharge is only applied to a part of the domain at a distance from the spring and can even exceed 100 % in these cases. Formally, this result arises from the occurrence of negative coefficients in the decomposition of the unit hydrograph into exponentially decaying components. The occurrence of negative coefficients also inhibits the simple interpretation as a set of linear reservoirs draining in parallel. Measuring the contribution of the slowest component from the peak of the unit hydrograph instead of the time at which the instantaneous recharge occurs reduces the contribution of the slowest component only slightly. This contribution is higher than for homogeneous aquifers and much higher than typically assumed for karst aquifers (less than 50 %).
Thus, we have to conclude that preferential flow patterns arising from a strongly organized pattern of transmissivity and storativity differ fundamentally from karst aquifers in their properties. For future work, it would be interesting to find out whether the difference mainly concerns the contribution of the slowest flow component or also the scaling of the recession coefficient with catchment size.
With respect to further development, an extension of the numerical scheme towards unconfined sloping aquifers (e.g., Rupp and Selker, 2006; Pauritsch et al., 2015) would be particularly useful. Although there is ongoing development in this field (e.g., Alemie et al., 2019; Pathania et al., 2019), including preferential flow patterns at a reasonable spatial resolution is still a challenge. Extending the implicit scheme for dendritic flow patterns towards unconfined sloping aquifers would still be challenging, but it might considerably contribute to understanding the response of hillslopes to precipitation events and phenomena such as subsurface storm flow (e.g., Chifflard et al., 2019).
In this section, Eqs. (42) and (43), which are the basis of the implicit scheme discussed in Sect. 2.5, are proven. Inserting Eqs. (14) and (41) into Eq. (40) yields
and thus
Using Eq. (14), we can then compute the flux according to
Then, ${q}_{i}^{\mathrm{0}}$ (Eq. 42) is obtained by setting ${h}_{b}(t+\mathit{\delta}t)={h}_{b}\left(t\right)$ and ${q}_{i}^{\prime}$ (Eq. 43) is obtained by taking the derivative with respect to h_{b}(t+δt).
Let us consider a 1D aquifer with a length L, where the spring is located at x=0 and the drainage divide at x=L. The boundary conditions are then $h(x,t)=\mathrm{0}$ at x=0 and $\frac{\partial}{\partial x}h(x,t)=\mathrm{0}$ at x=L, and h(x,t) is periodic with a wavelength of 4L. Thus, h(x,0) can be written as a Fourier series:
where the respective terms with the cosine function are zero due to the boundary conditions and a_{k}=0 for even values of k. The coefficients a_{k} are given by the relation
If we assume that the distal region, $\mathit{\lambda}L\le x\le L$ with $\mathit{\lambda}\in [\mathrm{0},\mathrm{1}]$, is initially filled to a given head value h_{0}, we obtain
for uneven values of k. The timedependent solution h(x,t) must satisfy the 1D version of Eq. (7) with r=0,
It is easily recognized that the solution of this equation with the initial condition defined by Eq. (B1) is
where
The flux per unit width across the boundary is then
with the coefficients a_{k} from Eq. (B4). The respective expression for the proximal region, $\mathrm{0}\le x\le \mathit{\lambda}L$, is readily obtained by subtracting this expression from the same expression with λ=0.
All codes and simulated data are available in a Zenodo repository at https://doi.org/10.5281/zenodo.7050521 (Strüven and Hergarten, 2022). The authors are happy to assist interested readers with reproducing the results and performing subsequent research.
JS developed the numerical codes and analyzed the results. SH developed the theoretical framework. Both authors wrote the paper.
The contact author has declared that neither of the authors has any competing interests.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
The authors would like to thank Erwin Zehe and an anonymous reviewer for their constructive comments and Monica Riva for the editorial handling.
This research has been supported by the Deutsche Forschungsgemeinschaft (grant no. 414728598).
This openaccess publication was funded by the University of Freiburg.
This paper was edited by Monica Riva and reviewed by Erwin Zehe and one anonymous referee.
Alemie, T. C., Tilahun, S. A., OchoaTocachi, B. F., Schmitter, P., Buytaert, W., Parlange, J.Y., and Steenhuis, T. S.: Predicting shallow groundwater tables for sloping highland aquifers, Water Resour. Res., 55, 11088–11100, https://doi.org/10.1029/2019WR025050, 2019. a
Atkinson, T. C.: Diffuse flow and conduit flow in limestone terrain in the Mendip Hills, Sommerset (Great Britain), J. Hydrol., 35, 93–110, https://doi.org/10.1016/00221694(77)900798, 1977. a
Banavar, J. R., Maritan, A., and Rinaldo, A.: Size and form in efficient transportation networks, Nature, 399, 130–132, https://doi.org/10.1038/20144, 1999. a
Basha, H. A.: Flow recession equations for karst systems, Water Resour. Res., 56, e2020WR027384, https://doi.org/10.1029/2020WR027384, 2020. a
BernatekJakiel, A. and Poesen, J.: Subsurface erosion by soil piping: significance and research needs, Earth Sci. Rev., 185, 1107–1128, https://doi.org/10.1016/j.earscirev.2018.08.006, 2018. a
Birk, S. and Hergarten, S.: Early recession behaviour of spring hydrographs, J. Hydrol., 387, 24–32, https://doi.org/10.1016/j.jhydrol.2010.03.026, 2010. a, b
Birk, S., Liedl, R., Sauter, M., and Teutsch, G.: Hydraulic boundary conditions as a controlling factor in karst genesis: A numerical modeling study on artesian conduit development in gypsum, Water Resour. Res., 39, 1004, https://doi.org/10.1029/2002WR001308, 2003. a
Carman, P. C.: Fluid flow through granular beds, Trans. Inst. Chem. Eng. Lond., 15, 150–166, 1937. a
Chifflard, P., Blume, T., Maerker, K., Hopp, L., van Meerveld, I., Graef, T., Gronz, O., Hartmann, A., Kohl, B., Martini, E., ReinhardtImjela, C., Reiss, M., Rinderer, M., and Achleitner, S.: How can we model subsurface stormflow at the catchment scale if we cannot measure it?, Hydrol. Process., 33, 1378–1385, https://doi.org/10.1002/hyp.13407, 2019. a
Costa, A.: Permeabilityporosity relationship: A reexamination of the Kozeny–Carman equation based on a fractal porespace geometry assumption, Geophys. Res. Lett., 33, L02318, https://doi.org/10.1029/2005GL025134, 2006. a
Cuthbert, M. O., Gleeson, T., Moosdorf, N., Befus, K. M., Schneider, A., Hartmann, J., and Lehner, B.: Global patterns and dynamics of climate–groundwater interactions, Nat. Clim. Change, 9, 137–141, https://doi.org/10.1038/s4155801803864, 2019. a, b, c
de Rooij, G. H.: Is the groundwater reservoir linear? A mathematical analysis of two limiting cases, Hydrol. Earth Syst. Sci. Discuss., 11, 83–108, https://doi.org/10.5194/hessd11832014, 2014. a
Drogue, C.: Analyse statistique des hydrogrammes de decrues des sources karstiques, J. Hydrol., 15, 49–68, https://doi.org/10.1016/00221694(72)900753, 1972. a
Edery, Y., Stolar, M., Porta, G., and Guadagnini, A.: Feedback mechanisms between precipitation and dissolution reactions across randomly heterogeneous conductivity fields, Hydrol. Earth Syst. Sci., 25, 5905–5915, https://doi.org/10.5194/hess2559052021, 2021. a, b
Enquist, B. J., Brown, J. H., and West, G. B.: Allometric scaling of plant energetics and population density, Nature, 395, 163–165, https://doi.org/10.1038/25977, 1998. a
Enquist, B. J., West, G. B., Charnov, E. L., and Brown, J. H.: Allometric scaling of production and lifehistory variation in vascular plants, Nature, 401, 907–911, https://doi.org/10.1038/44819, 1999. a
Fenicia, F., Savenije, H. H. G., Matgen, P., and Pfister, L.: Is the groundwater reservoir linear? Learning from data in hydrological modelling, Hydrol. Earth Syst. Sci., 10, 139–150, https://doi.org/10.5194/hess101392006, 2006. a
Fiorillo, F.: The recession of spring hydrographs, focused on karst aquifers, Water Resour. Manage., 28, 1781–1805, https://doi.org/10.1007/s112690140597z, 2014. a
Forkasiewicz, J. and Paloc, H.: Le régime de tarissement de la FouxdelaVis, Etude préliminaire, Chron. Hydrogéol., 3, 61–73, 1967. a, b
Hackbusch, W.: MultiGrid Methods and Applications, Springer, Berlin, Heidelberg, New York, https://doi.org/10.1007/9783662024270, 1985. a
Halmos, P.: FiniteDimensional Vector Spaces, in: 2nd Edn., Springer, New York, https://doi.org/10.1007/9781461263876, 1958. a
Hergarten, S.: Transportlimited fluvial erosion – simple formulation and efficient numerical treatment, Earth Surf. Dynam., 8, 841–854, https://doi.org/10.5194/esurf88412020, 2020. a
Hergarten, S. and Birk, S.: A fractal approach to the recession of spring hydrographs, Geophys. Res. Lett., 34, 11401, https://doi.org/10.1029/2007GL030097, 2007. a, b, c, d
Hergarten, S., Winkler, G., and Birk, S.: Transferring the concept of minimum energy dissipation from river networks to subsurface flow patterns, Hydrol. Earth Syst. Sci., 18, 4277–4288, https://doi.org/10.5194/hess1842772014, 2014. a, b, c, d, e, f, g, h, i, j, k
Hergarten, S., Winkler, G., and Birk, S.: Scale invariance of subsurface flow patterns and its limitation, Water Resour. Res., 52, 3881–3887, https://doi.org/10.1002/2015WR017530, 2016. a
Howard, A. D.: Theoretical model of optimal drainage networks, Water Resour. Res., 26, 2107–2117, https://doi.org/10.1029/WR026i009p02107, 1990. a
Jeannin, P.Y. and Sauter, M.: Analysis of karst hydrodynamic behaviour using global approaches: a review, Bull. Hydrogeol., 16, 31–48, 1998. a, b
Jeannin, P.Y., Artigue, G., Butscher, C., Chang, Y., Charlier, J.B., Duran, L., Gill, L., Hartmann, A., Johannet, A., Jourde, H., Kavousi, A., Liesch, T., Liu, Y., Lüthi, M., Malard, A., Mazzilli, N., PardoIgúzquiza, E., Thiéry, D., Reimann, T., Schuler, P., Wöhling, T., and Wunsch, A.: Karst modelling challenge 1: Results of hydrological modelling, J. Hydrol., 600, 126508, https://doi.org/10.1016/j.jhydrol.2021.126508, 2021. a
Kaufmann, G. and Braun, J.: Karst aquifer evolution in fractured, porous rocks, Water Resour. Res., 36, 1381–1391, https://doi.org/10.1029/1999WR900356, 2000. a
Kaufmann, G., Romanov, D., and Hiller, T.: Modeling threedimensional karst aquifer evolution using different matrixflow contributions, J. Hydrol., 388, 241–250, https://doi.org/10.1016/j.jhydrol.2010.05.001, 2010. a
Kleidon, A. and Renner, M.: Thermodynamic limits of hydrologic cycling within the Earth system: concepts, estimates and implications, Hydrol. Earth Syst. Sci., 17, 2873–2892, https://doi.org/10.5194/hess1728732013, 2013. a
Kleidon, A. and Savenije, H. H. G.: Minimum dissipation of potential energy by groundwater outflow results in a simple linear catchment reservoir, Hydrol. Earth Syst. Sci. Discuss. [preprint], https://doi.org/10.5194/hess2017674, 2017. a, b
Kleidon, A. and Schymanski, S. J.: Thermodynamics and optimality of the water budget on land: a review, Geophys. Res. Lett., 35, L20404, https://doi.org/10.1029/2008GL035393, 2008. a
Kleidon, A., Renner, M., and Porada, P.: Estimates of the climatological land surface energy and water balance derived from maximum convective power, Hydrol. Earth Syst. Sci., 18, 2201–2218, https://doi.org/10.5194/hess1822012014, 2014. a
Kleidon, A., Zehe, E., and Loritz, R.: ESD Ideas: Structures dominate the functioning of Earth systems, but their dynamics are not well represented, Earth Syst. Dynam. Discuss. [preprint], https://doi.org/10.5194/esd201952, 2019. a
Kovács, A.: Quantitative classification of carbonate aquifers based on hydrodynamic behaviour, Hydrogeol. J., 29, 33–52, https://doi.org/10.1007/s1004002002285w, 2021. a
Kovács, A. and Perrochet, P.: A quantitative approach to spring hydrograph decomposition, J. Hydrol., 352, 16–29, https://doi.org/10.1016/j.jhydrol.2007.12.009, 2008. a, b, c
Kovács, A. and Perrochet, P.: Well hydrograph analysis for the estimation of hydraulic and geometric parameters of karst and connected water systems, in: H2Karst Research in Limestone Hydrogeology, edited by: Mudry, J., Zwahlen, F., Bertrand, C., and LaMoreaux, J., Springer, Cham, 97–114, https://doi.org/10.1007/9783319061399_7, 2014. a, b
Kovács, A., Perrochet, P., Király, L., and Jeannin, P.Y.: A quantitative method for the characterisation of karst aquifers based on spring hydrograph analysis, J. Hydrol., 303, 152–164, https://doi.org/10.1016/j.jhydrol.2004.08.023, 2005. a, b, c, d
Kozeny, J.: Über kapillare Leitung des Wassers im Boden, Sitzungsber. Akad. Wiss. Wien, 136, 271–306, 1927. a
Maillet, E.: Essais d'Hydraulique souterraine et fluviale, Librairie Scientifique, A. Hermann, Paris, ISBN 9780282516925, 1905. a
Mangin, A.: Contribution a l'etude hydrodynamique des aquifères karstiques, PhD thesis, Université de Dijon, Dijon, https://hal.science/tel01575806 (last access: 17 August 2023), 1975. a, b, c, d
Maritan, A., Colaiori, F., Flammini, A., Cieplak, M., and Banavar, J. R.: Universality classes of optimal channel networks, Science, 272, 984–986, https://doi.org/10.1126/science.272.5264.984, 1996. a
Nutbrown, D. A.: Identification of parameters in a linear equation of groundwater flow, Water Resour. Res., 11, 581–588, https://doi.org/10.1029/WR011i004p00581, 1975. a, b
Padilla, A., PulidoBosch, A., and Mangin, A.: Relative importance of baseflow and quickflow from hydrographs of karst spring, Groundwater, 32, 267–277, https://doi.org/10.1111/j.17456584.1994.tb00641.x, 1994. a
Pathania, T., BottacinBusolin, A., Rastogi, A. K., and Eldho, T. I.: Simulation of groundwater flow in an unconfined sloping aquifer using the elementfree Galerkin method, Water Resour. Manage., 33, 2827–2845, https://doi.org/10.1007/s11269019022614, 2019. a
Pauritsch, M., Birk, S., Wagner, T., Hergarten, S., and Winkler, G.: Analytical approximations of discharge recessions for steeply sloping aquifers in alpine catchments, Water Resour. Res., 51, 8729–8740, https://doi.org/10.1002/2015WR017749, 2015. a
Rinaldo, A., RodriguezIturbe, I., Bras, R. L., IjjaszVasquez, E., and Marani, A.: Minimum energy and fractal structures of drainage networks, Water Resour. Res., 28, 2181–2195, https://doi.org/10.1029/92WR00801, 1992. a
RodriguezIturbe, I., Rinaldo, A., Rigon, R., Bras, R. L., IjjaszVasquez, E., and Marani, A.: Fractal structures as least energy patterns: The case of river networks, Geophys. Res. Lett., 19, 889–892, https://doi.org/10.1029/92GL00938, 1992a. a
RodriguezIturbe, I., Rinaldo, A., Rigon, R., Bras, R. L., Marani, A., and IjjaszVasquez, E.: Energy dissipation, runoff production, and the threedimensional structure of river basins, Water Resour. Res., 28, 1095–1103, https://doi.org/10.1029/91WR03034, 1992b. a
Rorabaugh, M. I.: Estimating changes in bank storage and groundwater contribution to streamflow, IAHS Publ., 63, 432–441, 1964. a, b
Rupp, D. E. and Selker, J. S.: On the use of the Boussinesq equation for interpreting recession hydrographs from sloping aquifers, Water Resour. Res., 42, W12421, https://doi.org/10.1029/2006WR005080, 2006. a
Savenije, H. H. G.: HESS Opinions: Linking Darcy's equation to the linear reservoir, Hydrol. Earth Syst. Sci., 22, 1911–1916, https://doi.org/10.5194/hess2219112018, 2018. a
Schroers, S., Eiff, O., Kleidon, A., Scherer, U., Wienhöfer, J., and Zehe, E.: Morphological controls on surface runoff: an interpretation of steadystate energy patterns, maximum power states and dissipation regimes within a thermodynamic framework, Hydrol. Earth Syst. Sci., 26, 3125–3150, https://doi.org/10.5194/hess2631252022, 2022. a
Sherman, L. K.: Stream flow from rainfall by the unit graph method, Eng. NewsRec., 108, 501–505, 1932. a
Strüven, J. and Hergarten, S.: Flow recession behavior of preferential subsurface flow patterns, Zenodo [code and data set], https://doi.org/10.5281/zenodo.7050521, 2022. a
West, G. B., Brown, J. H., and Enquist, B. J.: A general model for the origin of allometric scaling laws in biology, Science, 276, 122–126, https://doi.org/10.1126/science.276.5309.122, 1997. a
West, G. B., Brown, J. H., and Enquist, B. J.: A general model for the structure and allometry of plant vascular systems, Nature, 400, 664–667, https://doi.org/10.1038/23251, 1999a. a
West, G. B., Brown, J. H., and Enquist, B. J.: The fourth dimension of life: fractal geometry and allometric scaling of organisms, Science, 284, 1677–1679, https://doi.org/10.1126/science.284.5420.1677, 1999b. a
Westhoff, M. C. and Zehe, E.: Maximum entropy production: can it be used to constrain conceptual hydrological models?, Hydrol. Earth Syst. Sci., 17, 3141–3157, https://doi.org/10.5194/hess1731412013, 2013. a
Westhoff, M. C., Zehe, E., and Schymanski, S. J.: Importance of temporal variability for hydrological predictions based on the maximum entropy production principle, Geophys. Res. Lett., 41, 67–73, https://doi.org/10.1002/2013GL058533, 2014. a
Westhoff, M. C., Erpicum, S., Archambeau, P., Pirotton, M., and Dewals, B.: Maximum energy dissipation to explain velocity fields in shallow reservoirs, J. Hydraul. Res., 56, 221–230, https://doi.org/10.1080/00221686.2017.1289268, 2017. a
Xu, B., Ye, M., Dong, S., Dai, Z., and Pei, Y.: A new model for simulating spring discharge recession and estimating effective porosity of karst aquifers, J. Hydrol., 562, 609–622, https://doi.org/10.1016/j.jhydrol.2018.05.039, 2018. a
Zehe, E., Blume, T., and Blöschl, G.: The principle of maximum energy dissipation: a novel thermodynamic perspective on rapid water flow in connected soil structures, Philos. T. Roy. Soc. B, 365, 1377–1386, https://doi.org/10.1098/rstb.2009.0308, 2010. a
Zehe, E., Ehret, U., Blume, T., Kleidon, A., Scherer, U., and Westhoff, M.: A thermodynamic approach to link selforganization, preferential flow and rainfall–runoff behaviour, Hydrol. Earth Syst. Sci., 17, 4297–4322, https://doi.org/10.5194/hess1742972013, 2013. a
Zehe, E., Loritz, R., Edery, Y., and Berkowitz, B.: Preferential pathways for fluid and solutes in heterogeneous groundwater systems: selforganization, entropy, work, Hydrol. Earth Syst. Sci., 25, 5337–5353, https://doi.org/10.5194/hess2553372021, 2021. a
Zhao, J., Wang, D., Yang, H., and Sivapalan, M.: Unifying catchment water balance models for different time scales through the maximum entropy production principle, Water Resour. Res., 52, 7503–7512, https://doi.org/10.1002/2016WR018977, 2016. a
 Abstract
 Introduction
 Approach
 Results and discussion
 Conclusions and perspectives
 Appendix A: The fully implicit scheme for a dendritic network
 Appendix B: The unit hydrograph of a homogeneous 1D aquifer
 Code and data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
 Abstract
 Introduction
 Approach
 Results and discussion
 Conclusions and perspectives
 Appendix A: The fully implicit scheme for a dendritic network
 Appendix B: The unit hydrograph of a homogeneous 1D aquifer
 Code and data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References