We sincerely Daniele Pedretti (DP) for his thorough and thoughtful assessment of our 

manuscript.

DP: The manuscript ("Preferential Pathways for Fluid and Solutes in Heterogeneous 

Groundwater Systems: Self-Organization, Entropy, Work") provides a new framework 

that combines the use of free energy and entropy to characterize and quantify the 

emergence of preferential flow and channelled transport in heterogeneous media. 

Although the specific components of this framework are not novel themselves ,as 

correctly acknowledged by the authors, their combined use makes it a novel way that 

help disentangling open questions regarding the mechanisms of transport in 

heterogenous porous media.

The paper is well written. The objective, methodologies and conclusions are clear. I 

have summarized a line-by-line set of comments for the authors, which makes me 

recommending me accepting this manuscript after major revisions. 

RESPONSE: We thank the reviewer for the positive comments.

DP: My major concern is that, while I consider this approach excellently explained, it 

should have been demonstrated on a 3D heterogeneous system. Percolation 

thresholds are different in 2D and 3D systems. As such, the results of this study could 

have been very different if drawn from 2D or 3D stochastic models. I ask the authors to 

at least comment on this issue critically in their manuscript.

RESPONSE: We agree that a transfer of the proposed assessment to a three 3D 

stochastic media would yield interesting insights, which will for sure differ somewhat 

from what we found in 2D. However, we respectfully disagree that this might imply a 

different qualitatively behavior, as long as we work in a confined system (no flow 

boundary conditions for the upper, lower, inlet and outlet boundaries). The local 

changes in power arise from the local feedback on the pressure head gradient in front 

of the low conductivity bottlenecks. Gradients steepen ahead of the bottlenecks, which 

implies a higher power, locally. This feedback will also occur in a 3D confined system, as 

it is a direct result of the boundary conditions. It would not likely occur in an aquifer 

with a free surface neither in 2D or 3D. In the revised manuscript we will state that an 

analysis in 3 D is the next step of a forthcoming study, and explain why think that we 

expect qualitatively similar results due to the above stated reasons.

We furthermore like to stress, that percolation considerations are not relevant here, as 

the domains are well connected and well above a percolation threshold. Moreover, the 

original study by Edery et al. (2014), which we cite, presented already a critical path 

analysis in reference to the percolation threshold, based on the common assumption 

that preferential flows are a manifestation of percolation, controlled by the lower cut-

off for the hydraulic conductivity from which a path is possible. The limitations of 

percolation theory in evaluating the preferential flow are presented therein, and as 

such, the equivalence or connection between percolation and entropy is not 

straightforward.

DP: Thanking the authors for considering our 2017 WRR publication, I also suggest 

having a look at our follow-up manuscript (Bianchi and Pedretti 2018 WRR 

https://doi.org/10.1029/2018WR022827) where we extend our previous theory by 

computing the geological entrogram on evolving sampling scales. I think that most 

conclusions we got in those studies there are very much in line with those obtained 

through this study. Indeed, in the 2018 paper we also address the question of 2D vs 3D 

models, and described at page 4444 how solute particles tend to sample specific K 

clusters when travelling in the heterogeneous media. 

RESPONSE: We will be pleased to read your 2018 paper and refer to the study and 

conclusions in our revised manuscript, if/where appropriate. 

Best Regards 

 

Daniele Pedretti 

 

University of Milan - UNIMI (Italy).

Line-specific comments

DP: L72-76 I wouldn’t be so strict. Someone has succeeded in this task (Zhang et al 

2013 JH for instance). What is really complicated is finding a universal way to predict 

solute transport based on the aquifers geological structures. In Bianchi and Pedretti’s 

works on geological entropy we found an explanation for that: the lower the structure’s 

Shannon entropy, the more organized the flow and transport patterns in the field. In 

that set of works our aim was to start from the geology and not from self-customed 

flow fields (e.g. power-law distributed seepage velocities). 

RESPONSE: We thank DP for pointing this out, and will revise this passage to be less 

strict. 

DP: L97 "the probability of solutes to pass through HIGH (not low) conductivity regions". 

Please , fix it.

RESPONSE: We will fix this; thanks for noting the typo.

DP: L98 Please consider also our follow-up study (Bianchi and Pedretti 2018 WRR), 

where we study evolving scales of geological entropy rather than studying fixed-size 

blocks (as we did in the 2017 paper). In the 2018 paper we developed the concept of 

entrogram scale, which is also nicely correlated with the emergence of preferential flow 

and solute channelling.

RESPONSE: We will be pleased to read your 2018 paper and refer to the study and 

conclusions in our revised manuscript, if/where appropriate. 

DP: L101 enigmatic OK, emergent not really I would say.

RESPONSE: Agreed. We will delete the term emergent.

DP: L104 Again, I wouldn’t be so strict ("virtually impossible").  I’d rather just say that 

such predictions remain challenging.  For instance, Bianchi and Pedretti works or Zhang 

et al 2013 showed that it can be done. There is also a set of works by Rizzo and de 

Barros showing that predictions can be made starting from the aquifer structures. 

RESPONSE: Agreed. We will revise this passage accordingly. We note, too, that we can 

also achieve useful predictions of preferential transport of solutes in the partially 

saturated zone, in cases where detailed information about the pdf of macropores in 

soil and their connectivity are available. 

DP: L158 please see comment at L98.

RESPONSE: As for the L98 comment, we will read your 2018 paper and refer to the 

study and conclusions in our revised manuscript, if/where appropriate.

DP: L165-on. Rather than Objectives, these are Results and Conclusions.

RESPONSE: We agree that we provide here some foreshadowing on the results. The 

idea is to interest and motivate the reader to study, in particular, the theory section. In 

the revised manuscript, from the end of line 164 ("Specifically"), we will define the 

text as a new paragraph.

DP: L185 I wonder if all these nice concepts can be exported directly to 3D models, 

considering the different percolation thresholds between 2D and 3D models. Could the 

authors discuss on this?

RESPONSE: We think that the above mentioned local feedback on the head gradient 

will also occur in a 3D confined system. As such we expect qualitatively similar behavior 

in three 3D, as noted in our response to the second comment above. In the revised 

manuscript, we will include a brief consideration of these comments/responses in the 

conclusion and outlook section. 

DP: L216 why no local dispersion? This is a physical mechanism, which can substantially 

modify the solute pathways by increasing mixing and coalescence among the so-called 

"lamellas". Why neglecting it? I think this should have been investigated, from low to 

high Peclet numbers.

RESPONSE: The modeling approach used here is well accepted and accurate/sufficient 

for our purposes. Hydrodynamic dispersion is a macroscale fingerprint of diffusive 

transversal mixing of solute between flow lines of different fluid velocities. Adding a 

local dispersion term would incorporate implicit assumptions about subscale 

heterogeneity in velocities. 

DP: L256 I totally agree, but again, I think that transversal dispersion could have a big 

impact here. 

RESPONSE: We agree that transversal dispersion would of course work against 

steepening of lateral gradients. But lateral dispersion is the result of variability in the 

transversal flow field. We do not see much option for that, because the system is 

confined and flow is at steady state. 

DP: L388 see comment at L98

RESPONSE: We will address this point, as noted in our response to the L98 comment.

DP: L395 how this new theory could be connected to previously developed ones, such 

as the "lamella" description of solute transport in heterogenous media (e.g. 

https://doi.org/10.1017/jfm.2015.117) which also strongly depends on concentration 

gradients transversal to the main flow directions?  or the concepts of "least-resistance 

paths" (e.g.  https://doi.org/10.1002/2017WR020418)?

RESPONSE: We think these concepts are well connected. At the end of the day, an 

increase in spatial organization implies to a steepen gradient. There are many ways to 

express this, but there is only one mechanism to explain why this implies ordering and 

this is the second law of thermodynamics. This is because production of physical 

entropy implies essentially to deplete gradients.  

DP: L403-405 this looks like a conclusion of this work, rather than a result. Consider 

moving it to the appropriate sections

RESPONSE: Indeed we do a little bit of foreshadowing here, but we think this is helpful 

to keep the reader on track. We state clearly that "In the following, we demonstrate".

DP: L442 this is very similar to a conclusion by Bianchi and Pedretti 2018 (page 4444), 

which reads " These observations are further confirmed by the CDFs of the subsamples 

of K values, which show a significantly higher variability for 2-D fields compared to 3-D 

(Figures 6c and 6d). These results strengthen the knowledge that solutes tend to travel 

in the upper 15–20% of the K distribution (K classes 32–40 in Figure 6), which tend to 

fully percolate in 3-D correlated random fields regardless of their overall structure 

(Fogg et al., 2000; Fogg & Zhang, 2016; Harter, 2005). " Consider commenting on that

RESPONSE: We will happily include this point in our discussion here.

DP: L457 I  find it a bit complicated to relate "watts" to something related to groundwater. I 

mean, the entire derivation of power is clearly described in the previous sections, but 

as a hydrogeologist I have some problem to understand, for instance, if this is a 

high/low power. For instance would 2 watts be a high or lower power in this context?  

 

RESPONSE: We absolutely agree that power and W/m is not very common 

groundwater and vadose zone hydrology. Zehe et al. (2013) analyzed energy 

conversion associated with infiltration and found that macropores increase power in 

the infiltrating water flux. The maximum values during a rainfall event were of order 2 

W/m2. The use of W/m2 is much more common when dealing with the landsurface 

energy balance, as evaporation as water flux can be expressed as energy flux as well. In 

this context it is interesting to recall that the climatological land-surface energy balance 

is of order 100 W/m2. By comparison, a difference of 2 W/m2 is hence quite significant.  

In the revised manuscript, we will add text to include this background information, and 

comment on the significance of a 2 W per unit width increase as noted on L457.

DP: L455-478: these are great findings. It is particularly interesting that at some point 

the system behaves effectively as a 1D system. I’m wondering if they also hold for a 3D 

system, which is in general more percolated than a 2D one. 

 

RESPONSE: Thanks. You raise a very good point. We expect that the 3D system will 

deviate even more from the 1D approximation, but the local feedback on the head 

gradient will remain in a confined aquifer.  As mentioned above, we will address this 

point in the conclusions and outlook section of the revised manuscript.

DP: L487-488 for the same token, then solute injection mode (i.e resident vs flux-

averaged) could be also important to control the "effective" system power, right? Even if 

the flow field is the same, you change the way particles are already injected into certain 

flow zones. 

 

RESPONSE: Good point, thanks. We will mention this in the revised manuscript.

DP: L505-507 this is well aligned with the result of Bianchi and Pedretti 2017,2018, who 

expressed geological entropy with the BTC moments. 

 

RESPONSE:  Agreed. We will mention this in the revised manuscript.

DP: L537 also similar to Bianchi and Pedretti 2017,2018.  

 

RESPONSE: Agreed, See reply to the previous comment.

DP: L558 notice an underlined word: just a writing typo? 

 

RESPONSE: Indeed, we will fix this.

DP: L564 and the 2018 work, which extends the previous one. 

 

RESPONSE: Yes, we will read this work and refer to it where appropriate. 

DP: L571 I agree that the CTWR beta could be a good indicator to connect to entropy. In 

general, any indicator that explains the departure from the BTC symmetry could be 

also good. This is why in Bianchi and Pedretti 2017, 2018 we did not fit a power-law 

curve to our BTCs, but limited ourselves to the third moment of the BTC, which is 

strongly correlated with geological entropy indicators. However, in Pedretti and Bianchi 

2018 ADWR (https://doi.org/10.1016/j.advwatres.2018.01.023) we found that for a 

system with very-low geological entropy (i.e. high spatial order) all BTC tailing tended to 

the same power-law value, close to one. Any comment on that?   

 

RESPONSE: We agree that various measures of non-symmetry in the BTC may act as 

good measures of non-Fickian transport. Zehe and Flühler (2001), for example, used 

also the skewness of vertical travel distance as a measure of preferential flow in 

infiltration patterns. In fact, one could also use, for example, the difference in the 

entropies of the actual BTC and the (homogeneous case) Fickian BTC. However, these 

additional comments are beyond the focus of our current study.

DP: L604 I think it could have been worth looking at the probability of transitioning of a 

particle from a low to a high flow zone, and relate it to the derived power P. 

 

RESPONSE: This is certainly an interesting point, but beyond the scope of the current 

study. We plan to explore this and related aspects in a future study.

DP: L606-608 to me the framework should be also demonstrated on a 3D simulated 

aquifer to properly claim that this approach "holds the keys" to disentangle a problem 

that have been a great challenge for many researchers until now. I may agree with that, 

and it would be amazing, but it has to be proved.

RESPONSE: We agree that further investigations in 3D are needed, but we also 

encourage study of systems in which K is evolving. We will modify the text here to note 

these points, and to express that this approach potentially holds the key.

On behalf of all co-authors I thank Daniele Pedretti again for his helpful comments.

Best regards, Erwin Zehe 

References:

Zehe, E., and Fluhler, H.: Slope scale variation of flow patterns in soil profiles, Journal of 

Hydrology, 247, 116-132, 2001.

Zehe, E., Ehret, U., Blume, T., Kleidon, A., Scherer, U., and Westhoff, M.: A 

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

behaviour, Hydrology And Earth System Sciences, 17, 4297-4322, 10.5194/hess-17-

4297-2013, 2013.

We sincerely Hubert Savenije (HS) for his thorough and encouraging assessment of our 

manuscript. We also thank him for the highlighted edits, which we will happily include in the 

revised manuscript.

We agree that the interplay of dissolution and precipitation of minerals such as silicate or 

carbonate rock, and the related local feedbacks on saturated hydraulic conductivity, will 

certainly affect and change the distribution of entropy and power in fluid flow (and in flow of 

chemicals). Aspects of dissolution and precipitation and their feedbacks on hydraulic 

conductivity, flow patterns and chemical transport are currently addressed in hess-2021-238 

by Edery et al. (in review 2021). A precondition for dissolution of e.g., carbonate, is to 

maintain a local "saturation deficit", i.e., the actual carbonic acid concentration in relation to 

the pH must be lower than the local equilibrium concentration established by the pH. Such 

local reaction-limited conditions arise if the low pH is funneled towards the preferential flows 

and is not distributed homogeneously over the inlet. This is expected to happen in high 

conductivity regions, where preferential flow occurs, and one could imagine that these regions 

grow preferentially due to the dissolution (maybe even backwards as rill systems) and 

establish a more connected drainage network. Precipitation requires exactly the opposite 

conditions: local oversaturation of carbonic acid relative to the pH, because the concentration 

is larger than the equilibrium concentration. This could happen upstream of "low conductivity 

bottle necks" or perpendicular to the preferential flows, which might imply that conductivity, 

locally, declines even more. 

Overall, this could imply that preferential pathways become more preferential, while local 

bottle necks become even more "narrow". Such a system could, overall, still evolve to a more 

organized dynamic behavior and we agree with HS that the key to assess this is to include 

molar entropy " but also free energy differences associated with the chemical reactions and 

chemical energy fluxes associated with chemical transport " into the entropy and energy 

balances. 

We will reflect on these aspects in the conclusion and outlook section of the revised 

manuscript; we are currently working on them as a continuation of our study and that of 

Edery et al. (in review 2021).

Reference:

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. Discuss. [preprint], https://doi.org/10.5194/hess-2021-238, in review, 2021.