Numerical study of Darcy's law of yield stress fluids on a deep tree-like network
Stéphane Munier Alberto Rosso
















arrows.meta,arrows

/tikzfeynman/warn luatex=false


    PP_satP_01]Stéphane Munier2]Alberto Rosso[1]CPHT, CNRS, École polytechnique, Institut Polytechnique de Paris, 91120 Palaiseau, France[2]LPTMS, CNRS, Université Paris-Saclay, 91405 Orsay, FranceNumerical study of Darcy's law of yield stress fluids
    
    on a deep tree-like network
        [
        September 9, 2024
    ==================================================================================
    
    § ABSTRACT 
    Understanding the flow dynamics of yield stress fluids in porous media presents a substantial challenge. Both experiments and extensive numerical simulations frequently show a non-linear relationship between the flow rate and the pressure gradient, deviating from the traditional Darcy law. In this article, we consider a tree-like porous structure and utilize an exact mapping with the directed polymer (DP) with disordered bond energies on the Cayley tree. Specifically, we adapt an algorithm recently introduced by Brunet et al. [Europhys. Lett.  131, 40002 (2020)] to simulate exactly the tip region of branching random walks with the help of a spinal decomposition, to accurately compute the flow on extensive trees with several thousand generations. Our results confirm the asymptotic predictions proposed by Schimmenti et al. [Phys. Rev. E  108, L023102 (2023)], tested therein only for moderate trees of about 20 generations.
    
    
    § INTRODUCTION
    
    
    The Darcy law describes the flow rate Q of a Newtonian fluid through a cylinder of radius R and length L, filled with a porous medium <cit.>:
    
        Q = κ R^2 P/η L .
    
    In this equation, P is the pressure drop between the two ends of the cylinder, and η is the viscosity of the fluid. The parameter κ, known as permeability, is dimensionally an area that characterizes the medium's ability to transmit fluid <cit.>. In the absence of a porous medium, the flow rate is described by the famous Poiseuille's law:
    
        Q_Pois.(R) = π R^4 P/8 η L .
    
    The permeability of a filled cylinder corresponds to an area much smaller than the  section, π  R^2. A simple model, proposed by H. Darcy, is based on the assumption that flow in a porous medium occurs through numerous thin, non-intersecting channels, each with a radius R_ ch≪ R and of number n^ ch per unit area. According to Poiseuille's law, the flow through the filled cylinder can be written as Q = π R^2 n^ ch× Q_Pois.(R_ ch). This leads to an expression for permeability as κ = (π/8) n^ ch R_ ch^4, which can be rewritten as the product of two factors: (i) the fraction of unit area occupied by the fluid, n^ ch R_ ch^2, a positive number smaller than one, and (ii) the area of the section of the microchannel π R_ ch^2. Although actual porous media feature a complex network of intersecting channels with varying shapes, the Darcy law holds true as long as the number of open channels n^ ch remains independent of pressure.
    
    This Darcy framework is significantly altered in the presence of yield stress fluids—fluids that behave like a solid below a given yield stress, σ_Y, and flow only when the stress exceeds this threshold. The consequence of yield stress is that the number of channels through which the fluid flows increases as the pressure drop applied at the ends of the throat increases. Three regimes are observed in experiments and numerical simulations:
    (i) no flow is observed below a critical pressure drop ;
    (ii)  a non linear flow is observed for P><cit.>, the non linearity being related to the growing number of open channels above threshold;
    (iii) only above a saturation pressure ≫, does the flow revert to linear growth,  the linear growth with the standard permeability κ <cit.>.
    
    Pore-network models are the simplest models that capture these three regimes. They feature a graph structure in which the nodes represent large open pores with well-defined pressure, and the edges are throats  connecting the pores. Establishing how the Poiseuille law is modified in the presence of a yield stress remains an open challenge. In this study, we examined the simplest modification to the Poiseuille law:
    
    Q= π R^4/8η L×(P-P'-τ)_+ ,
    
    where (x)_+≡max(0,x) and  the  threshold pressure drop is τ=L σ_Y /R<cit.>.
    
    
    
    We shall assume that all throats have the same geometry. Here τ represents the minimum pressure difference required to establish flow through the throat. We randomly draw these thresholds as independent and identically distributed variables.
    
    Following Refs. <cit.>, the hydraulic network is organized as a single inlet throat feeding a perfect binary tree with a total height of T-1 (see Fig. <ref>). The set of the 2^T-1 variables τ that specify each sample can be interpreted as “frozen disorder” in statistical physics terms.
    
    The tree-like network of pipes maps exactly onto the configuration space of directed polymers (DP) <cit.> in the infinite-dimensional limit. The elementary throats correspond to bonds between monomers, and the τ variables  represent their energies in that context. Therefore, as noted in Refs. <cit.>, many tools developed for the DP problem can be adapted to the Darcy problem. Understanding the former has benefited from the observation that configuration spaces are generated by branching random walks in the universality class of branching Brownian motion <cit.>. The main goal of this paper is to demonstrate how an algorithm developed in the latter context <cit.> can be applied to this fluid-mechanics problem to numerically investigate it on very large tree-like networks. Specifically, this algorithm will allow us to test a prediction proposed in Ref. <cit.>, which led to asymptotic analytical expressions for a set of relevant observables.
    
    
    
    
    
    The paper is organized as follows. Section <ref> provides an alternative and pedagogical derivation of the algorithm used in Ref. <cit.>. This algorithm computes the flow curve for moderate trees T∼ 20 because it requires the prior generation of the full network (2^T-1 branches). In Section <ref>, we introduce a new algorithm to compute the flow curve, inspired by the one described in Ref. <cit.>. Our new algorithm can compute the flow curve for trees with T∼ 10^3 because it draws only the sequence of the first open channels, without prior generation of the full network. Our results are presented and discussed in Section <ref>, and our conclusions are drawn in Section <ref>. Two appendices provide technical details.
    
    
    
    
    
    
    
    § BASIC ALGORITHM TO COMPUTE THE FLOW
    
    
    According to the discussion in the Introduction, we  model the pressure dependence of the flow in each elementary throat using  Eq. (<ref>). By an appropriate choice of units, we can set the overall constant to unity, to get
    
    Q=(P-P'-τ)_+ .
    
    Furthermore, the Kirchhoff's law <cit.> is assumed at each node of the tree that connects two outgoing throats to an incoming one. Specifically, the sum of the outgoing flows Q^(0) and Q^(1) is equal to the incoming flow Q:
    
    Q=Q^(0)+Q^(1).
    
    Equations (<ref>) and (<ref>) are sufficient to determine the properties of the flow for any disorder realization.
    
    The goal of this section is to design the simplest algorithm to compute the flow curve Q_T(P)  as a function of the pressure P applied at the inlet, for a given realization of the disorder in a network tree of height T. 
    
    
    
    
    
     §.§ Warm-up: the simple trees with T=1 and T=2
    
    
    We begin by presenting detailed analytical calculations of the flow in elementary cases with T=1 and T=2, formulated in a way that facilitates understanding the more general case.
    
    In the  case T=1, there is only one throat, characterized by the threshold τ_01; see Fig. <ref>. The flow as a function of the pressure is then just given by Eq. (<ref>), in which P' is the pressure at the free outlet of the network, namely P'=0:
    
    Q_1(P)=
    0    for P≤τ_01,
    
    P-τ_01   for P>τ_01.
      [baseline=0.4cm]
    [->] (0,0)–(0,1) node [anchor=south]Q_1(P);
    [->] (0,0)–(2,0) node [anchor=west]P;
    [thick] (0,0)–(1,0)–(2,1);
     (1,-0.1)node[anchor=north]P_0=τ_01–(1,0.1);
    
    Throughout, we denote P_0 as the disorder-dependent minimum pressure that needs to be applied at the inlet of a network to establish the flow. In the present case, obviously, P_0=τ_01.
    
    In the case T=2, the network has one node connecting three throats; see Fig. <ref>. 
    
    Let us call P' the pressure at the node, and call τ_01, τ_12^(0) and τ_12^(1) the thresholds that characterize the incoming and the two outgoing throats respectively. For the fluid to flow, the two following conditions must be satisfied: 
    
    P-P'>τ_01  and   P'>min(τ_12^(0),τ_12^(1)).
    
    This implies P>P_0, where
    
    P_0=τ_01+min(τ_12^(0),τ_12^(1)).
    
    In words, Eq. (<ref>) says that P_0 minimizes the sum of the thresholds τ along the directed paths connecting the root to the leaves of the tree. This definition of P_0 actually generalizes to trees of any height.
    
    Let us set the pressure P at the inlet just above the threshold P_0. Assuming for the time being that τ_12^(0)≠τ_12^(1) and hence that P_0 is not degenerate, the flow reads
    
    Q_2(P)=P-P'-τ_01=P'-min(τ_12^(0),τ_12^(1)).
    
    The first equality is just Eq. (<ref>) applied to the incoming throat. The second one follows from the conservation of the flow at the node Eq. (<ref>), in the case in which only one outgoing channel is open, and from Eq. (<ref>) applied to the open outgoing throat. We may solve for P', arriving at the following determination:
    
    P'=1/2[P-τ_01+min(τ_12^(0),τ_12^(1))].
    
    It is then straightforward to deduce Q_2(P):
    
    Q_2(P)=1/2[P-τ_01-min(τ_12^(0),τ_12^(1))]_+=1/2(P-P_0)_+.
    
    This equation holds for pressures below the threshold P_1 for the opening of the second channel. 
    
    The flow above P_1 is determined using again Eq. (<ref>) and the Kirchhoff law (<ref>) at the node, now assuming that both channels are open:
    
    Q_2(P)=P-P'-τ_01=2P'-τ_12^(0)-τ_12^(1).
    
    Solving for P', 
    
    P'=1/3(P+τ_12^(0)+τ_12^(1)-τ_01).
    
    We evaluate the threshold pressure P_1 based on the requirement that, for the flow to establish in both branches of the tree, the pressure P' at the node must exceed the threshold pressures of both outgoing throats:
    
    P'≥max(τ_12^(0),τ_12^(1)).
    
    The pressure P_1 coincides with the value of the pressure P at the inlet such that this inequality for P' is saturated. Equation (<ref>) then yields
    
    P_1=P_0+2[max(τ_12^(0),τ_12^(1))
    -min(τ_12^(0),τ_12^(1))].
    
    Now, inserting the expression for P' in Eq. (<ref>) into one of the equations (<ref>) for Q_2(P), we get the expression of the flow:
    
    Q_2(P)=2/3[P-τ_01-1/2(τ_12^(0)+τ_12^(1))].
    
    All in all, we have found that the flow depends on the applied pressure as
    
    Q_2(P)=
    0    for P≤ P_0
    1/2[P-τ_01-min(τ_12^(0),τ_12^(1))]    for P_0<P≤ P_1
    2/3[P-τ_01-1/2(τ_12^(0)+τ_12^(1))]
       for P>P_1,
      [baseline=0.6cm]
    [->] (0,0)–(0,1.5) node [anchor=south]Q_2(P);
    [->] (0,0)–(2.5,0) node [anchor=west]P;
    [thick] (0,0)–(0.5,0)–(1.5,0.5)–(2.5,0.5+1*2/3);
    [dashed] (1.5,0.5)–(2.5,1);
     (0.5,-0.1)node[anchor=north]P_0–(0.5,0.1);
     (1.5,-0.1)node[anchor=north]P_1–(1.5,0.1);
    [dotted](1.5,0.1)–(1.5,0.5);
    
    where the thresholds P_0 and P_1 are given in Eqs. (<ref>) and (<ref>) respectively. Let us comment that the degenerate case τ_12^(0)=τ_12^(1) does not require a special treatment: P_0=P_1=τ_01+τ_12^(0), and the formula we have just established applies, the second distinguished case being simply no longer necessary as the pressure interval in which it applies becomes trivial.
    
    We see that Q_2(P) is a piece-wise linear and continuous function. It is not difficult to figure out that this property generalizes for trees of arbitrary heights.
    
    
    
    
    
     §.§ Trees of arbitrary heights
    
    
    Given a network of size T with all the 2^T-1 elementary thresholds τ, we take the following generic Ansatz for the flow at the inlet of the complete tree, or of sub-trees of it, as that of an effective Darcy law:
    
    Q(P)=κ_eff(P)[P-P^*_eff(P)].
    
    We shall call the proportionality factor κ_eff(P) as the  effective permeability of the considered (sub-)network, and P^*_eff(P) as the  effective offset. These two flow parameters   are non-decreasing piece-wise constant functions of P.  Their  discontinuities occur at the pressures  at which new channels  open. Our goal is to compute recursively κ_eff(P) and P_eff^*(P) for a pressure P such that there are, say, n^ch open channels. 
    
    
    
    
    
      §.§.§ Flow parameters at a given node
    
    
    Let us start by assuming that the sub-tree of all open channels is known, and consider a generic node in the latter. We label the throats it connects as in Fig. <ref>, and write the Kirchhoff law (<ref>) at this node.
    
    We then replace the flows in the right-hand side of the latter by Eq. (<ref>), with P set to the pressure P' at the node and with appropriate superscripts for the parameters. As for Q in the left-hand side of Eq. (<ref>), we shall eventually replace it by Eq. (<ref>). 
    
    We first address the case in which only one of the outgoing channels is open, say the one labeled (0):
    
    Q=P-P'-τ=κ_eff^(0)(P'-P^*(0)_eff).
    Q in the left-hand side will eventually depend on P only. The effective permeabilities and offsets in the right-hand side depend on P' as piece-wise constant functions. Hence in the restricted range of P (and in the corresponding one of P') of interest, the second and third members form a linear equation for P'. Solving the latter, inserting the solution back into one of the alternative expressions for Q, we get for the flow
    
    Q=κ_eff^(0)/1+κ_eff^(0)(P-τ-P_eff^*(0)).
    
    Recalling that Q(P) itself obeys Darcy's law (<ref>), it is straightforward to identify the parameters that characterize the flow at the inlet:
    κ_eff^(0∧1)=κ_eff^(0)/1+κ_eff^(0)    
    P_eff^*(0∧1)=τ+P_eff^*(0) .
    
    We introduced self-explanatory superscripts. It is of course enough to exchange the superscripts (0) and (1) to get the case in which sub-channel (1) is open and sub-channel (0) is closed
    
    Next, let us assume that P is sufficiently large, such that the two outgoing channels are open. Instead of Eq. (<ref>), we now have the following equation:
    
    Q=P-P'-τ=κ_eff^(0)(P'-P^*(0)_eff)+κ_eff^(1)(P'-P^*(1)_eff).
    
    Solving for P' and inserting its solution back into the expression for Q, we arrive at
    
    Q=κ_eff^(0)+κ_eff^(1)/1+κ_eff^(0)+κ_eff^(1)(P-τ-κ_eff^(0)P_eff^*(0)+κ_eff^(1)P_eff^*(1)/κ_eff^(0)+κ_eff^(1)).
    
    Comparing to Eq. (<ref>), we easily obtain the following relations between the effective flow parameters:
    κ_eff^(0∧1)=κ_eff^(0)+κ_eff^(1)/1+κ_eff^(0)+κ_eff^(1)    
    P_eff^*(0∧1)=τ+κ_eff^(0)P_eff^*(0)+κ_eff^(1)P_eff^*(1)/κ_eff^(0)+κ_eff^(1).
    
    We take the convention that the effective permeability is null for closed channels. We easily check that in the case in which only channel (0) is open, Eq. (<ref>) boils down to Eq. (<ref>). Furthermore, the open throats located at the maximal depth, i.e. at the leaves of the sub-tree of open channels, have effective permeability κ_eff=1, and effective threshold P_eff^* equal to the elementary threshold of the considered throat. Then, iterating Eqs. (<ref>) and (<ref>) upward from the leaves of the open sub-tree to its root enables one to determine the parameters κ_eff and P_eff^* of all sub-channels, and, eventually, of the complete tree.
    
    
    
    
    
      §.§.§ Threshold pressures
     
    
    We recall that the expressions of the flow parameters just established are valid for pressures in the limited range in which the sub-tree we have singled out is that of all open channels. Let us now relate the threshold pressures of sub-trees rooted in the different throats joining at a given node. 
    
    We call P^(0)_thr and P^(1)_thr the minimum pressures to apply at the inlet of the outgoing throats, for the sub-trees they root to be independently open. We compute the corresponding minimum pressure P^(0∧1)_thr required at the inlet of the incoming throat to have a flow in both sub-trees. As above, we write the incoming flow for that pressure in two different ways:
    
    Q=P^(0∧1)_thr-max(P^(0)_thr,P^(1)_thr)-τ
    =(κ^(0)_eff+κ_eff^(1))max(P^(0)_thr,P^(1)_thr)-(κ^(0)_effP_eff^*(0)+κ^(1)_effP_eff^*(1)) .
    
    From the equality of the second and third members in this equation, we deduce the value of the threshold pressure:
    
    P^(0∧1)_thr=τ+(1+κ_eff^(0)+κ_eff^(1))max(P^(0)_thr,P^(1)_thr)-(κ_eff^(0)P_eff^*(0)+κ_eff^(1)P_eff^*(1)).
    
    
    The case in which only one outgoing channel, say (0), is open is obtained by substituting max(P^(0)_thr,P^(1)_thr) with P^(0)_thr, and setting κ_eff^(1) to zero in this equation. The corresponding threshold pressure, that we shall denote by P^(0∧1)_thr, reads
    
    P^(0∧1)_thr=τ+P^(0)_thr+
    κ_eff^(0)(P^(0)_thr-P_eff^*(0)).
    
    Obviously, the threshold pressure P^(0∧1)_thr, corresponding to the case in which (0) is closed and (1) is open, is obtained by interchanging the superscripts (0) and (1) in this formula.
    
    It will prove useful to have in mind the particular case of one-branch sub-trees, for which the threshold pressures P_thr coincide with the effective pressures P^*_eff. According to Eq. (<ref>), the latter just amount to the sums of the elementary thresholds τ along the branch of the considered sub-trees.
    
    We iterate Eqs. (<ref>) and (<ref>) from the leaves to the root to determine the minimum pressure that needs to be exerted at the inlet of the network to open all channels of that sub-tree.
    
    
    
    
    
      §.§.§ One- and two-branch trees
    
    
    Here, we shall work out the expression for the flow parameters and of the threshold pressures in the cases n^ch=1 and 2, but now for a tree of arbitrary height T. 
    
    Considering a generic directed path P that connects the root of the tree to some leaf, and its segment P_[t_1,t_2] between the depths t_1 and t_2, we shall introduce the convenient notation
    ε^( P)_t_1 t_2≡∑_I∈ P_[t_1,t_2]τ_I
    
    for the path sum of the elementary thresholds along this segment. The indices I over which the sum runs label the elementary segments the path is constructed of. We will also be led to use a notation such as 
    min_ P_[t,T]ε_0 T^( P) ,
    
    by which we mean the minimum of the quantity ε_0 T^( P), varying the segment of the directed path P between the depths t and T (i.e. the nodes between t+1 and T) in all possible ways on the given tree, while keeping fixed the segment and all its nodes between the depths 0 and t.
    
    Let us address the single-branch case n^ch=1, consisting of the directed path P_0. We may easily compute the flow parameters for sub-channels rooted at any node along the branch P_0. Iterating Eqs. (<ref>), we get
    κ_eff t=1/T-t,   P^*_eff t=ε^( P_0)_t T
    for the sub-channel P_0 [t,T] rooted at the node located at depth t on P_0, where t is such that 0≤ t<T. The corresponding threshold pressures P_thr t identify to P^*_eff t, as observed in the previous section. 
    
    If we want to require P_0 to be the first path (or one of the first paths, in case of degeneracies) that opens when the pressure at the network inlet is increased from zero, it must obviously be a path that minimizes P_eff t=0^*=ε^( P_0)_0 T. We denote the minimum effective threshold by P_0. If there is only one path that realizes the minimum, which is always the case if the distribution p(τ) of the τ's is continuous, then the flow as a function of the pressure P at the inlet follows from Darcy's law (<ref>), with the parameters read off the previous equation after having set t=0:
    
    Q_T(P)=1/T(P-P_0)_+.
    
    This relation holds only below the pressure above which the next channel opens.
    
    We now turn to the case n^ch=2. We shall assume that the second branch is attached at a node at depth t_1, and we denote by P_1t̂_1 the directed path from the root to the leaves of the tree that includes this branch. This path has obviously an overlap of size t_1 with P_0, which amounts to all throats between the inlet of the network and the node at depth t_1. 
    
    The flow parameters of the sub-trees rooted in nodes of depth t≥ t_1 are just those of the single-branch case for each of the channels. If t<t_1 instead, we apply Eq. (<ref>), and then Eq. (<ref>), possibly multiply. We find
    κ_eff t=2/T+t_1-2t,   P_eff t^*=ε^( P_0)_t t_1+1/2(ε^( P_0)_t_1 T+ε^( P_1t̂_1)_t_1 T).
    
    As for the threshold pressure to open the full two-channel sub-network rooted at the node at depth t<t_1, we check that the recursion built from Eqs. (<ref>),(<ref>) is solved by
    
    P_thr t=ε^( P_0)_t T+T-t/T-t_1(ε^( P_1t̂_1)_t_1 T-ε^( P_0)_t_1 T).
    
    
    The path that opens first when the pressure at the network inlet is increased from P_0 must obviously be a path that minimizes P_thr  t=0. Hence the expression of the second threshold pressure P_1 reads
    
    P_1≡ P_0+ min_t_1', P'_1t̂_1' [t_1'+1,T]T/T-t_1'(ε^( P'_1t̂_1')_t_1' T-ε^( P_0)_t_1' T).
    
    From now on, what we shall call P_1 will be a path that realizes this minimum of the threshold pressure, and t_1 will stand for its overlap with the path P_0. Note that in formula (<ref>), we may replace each sub-path sum ε^( P_0)_t_1 T,ε^( P_1)_t_1 T by the full path sum ε^( P_0)_0 T,ε^( P_1)_0 T: the difference would remain unchanged, since P_0 and P_1 overlap on [0,t_1].
    
    Applying to the network a pressure P chosen above P_1 but below the threshold pressure P_2 for the opening of the next channel (assuming P_2>P_1), the P-dependence of the flow reads
    
    Q_T(P)=2/T+t_1[P-ε_0 t_1^( P_0)-1/2(ε^( P_0)_t_1 T+ε^( P_1)_t_1 T)].
    
    We check that this formula boils down to Eq. (<ref>) for T=2 and when the overlap t_1 is set to its only possible value in this case, namely 1. 
    
    Alternatively, Q_T(P) may be expressed in terms of the threshold pressures and of the overlap t_1:
    
    Q_T(P)=2/T+t_1(P-T+t_1/2TP_0-T-t_1/2TP_1).
    
    This formula gives also the correct flow in the case so far disregarded in which the first two channels P_0, P_1 are degenerate and open simultaneously, namely if P_1=P_0.
    
    We may use the recursions we have established in this section and generalize the calculation we have just performed for the case of one and two branches to an arbitrary number of them. Based on the same method, it is now easy to design an algorithm that automates the search for open channels in a given tree-like network, and thus the calculation of the flow. This is what we will expose now.
    
    
    
    
    
     §.§ Algorithm to determine the open sub-tree and the corresponding flow
    
    
    We start with a fully-generated tree-like network of height T. Given this frozen disorder, we want to determine the sub-trees of open channels at any given pressure at the inlet, the thresholds P_0,P_1,⋯ at which the channels open successively, as well as the flow function Q(P). In this section, we shall follow closely in spirit the method exposed in Ref. <cit.>. 
    
    To begin, we search for a path P_0 that minimizes ε^( P_0)_0 T. Let us first assume, for simplicity, that the values of ε^( P)_0 T on all possible paths P in the tree are all distinct. Then, within this assumption, there is a single channel with lowest threshold pressure P_0. Once this first channel is determined, the flow parameters entering the Darcy law (<ref>) for this single-channel network are calculated iterating Eq. (<ref>) from the leaves to the root.
    
    To find the next channels to open and the next threshold pressure P_1, we loop over all nodes of the path P_0, namely over the depth t_1, with 0<t_1<T. For each t_1, we search for the path P_1 that minimizes ε_t_1+1 T^( P_1). Then, using Eqs. (<ref>) and (<ref>), we compute the minimum pressure P_thr 0 that has to be applied at the inlet of the full network to open the second channel at the outlet of this very node, assuming that all other closed channels of the tree remain closed. We eventually open the channel that minimizes this threshold pressure. We arrive at a two-branch subtree, such as the one singled out in Fig. <ref>. The flow parameters are subsequently computed by iterating Eqs. (<ref>) and (<ref>) going up the tree, from the leaves to the root.
    
    We repeat the steps just described to the two-branch sub-tree we have determined in order to search for the next channel to open. We test the opening pressure of all channels attaching to that sub-tree. These threshold pressures are determined iterating Eqs. (<ref>) and (<ref>). We open the channel endowed with the lowest threshold and compute its flow parameters.
    
    We iterate this overall procedure for the next branches, until some predefined stopping condition is satisfied, depending on the observable we aim at measuring. We see that in this way, we are eventually able to compute the flow Q_T(P) for any value of P.
    
    In the case in which there are n channels that have the same threshold pressure, the algorithm can be kept unchanged: we just open these degenerate channels successively, in random order to each other.
    
    
    
    
    
      
    Reduction to the directed polymer problem.
    
    This algorithm can be trivially modified to solve the problem of the determination of the successive configurations of the directed polymer problem, from the ground state to higher excited states. Instead of computing the threshold pressures P_0, P_1, P_2,⋯ and opening the channels in increasing order of these thresholds, we order the directed paths P_0, P_1, P_2⋯ on the tree in increasing sequence ε_0≡ε^( P_0)_0 T, ε_1≡ε^( P_1)_0 T, ε_2≡ε^( P_2)_0 T⋯ of the path sums of the τ's. The latter quantities correspond to the energies of the directed polymer, and the associated paths P_0, P_1, P_2,⋯ represent the corresponding polymer configurations.
    
    
    
    
    
    § OPTIMIZED ALGORITHM FOR SUCCESSIVE CHANNEL OPENING
    
    
    The algorithm just described requires the prior generation of the network. The complexity of the generation of a tree of height T grows like the number of bonds, namely exponentially in T. This definitely restricts the values of T that may be reached, in any practical implementation, to a few dozen units. With such an algorithm, whatever the computing power available, it will never be possible to increase T by more than a factor of order 1.
    
    However, there is a much more efficient way to proceed. Instead of generating all elementary thresholds τ a priori, we may grow the trees branch-by-branch. The new algorithm we shall introduce is inspired from that described in Ref. <cit.>, which was designed to generate exactly the tip region of a branching random walk at large times. The initial motivation came from particle physics: the problem there was to be able to efficiently generate scattering configurations relevant to electron-nucleus collisions, in some suitable asymptotic high-energy regime. (The interested reader is referred to <cit.> for a recent review on the connections between scattering in particle physics and general branching processes). But the algorithm may also be used to generate the lowest-lying energy levels of a directed polymer in a random medium in the mean-field approximation, which turns out to be a problem in the same class. Notably, the underlying method can be interpreted as a spinal decomposition of branching random walks, a technique previously recognized in mathematical literature (see e.g. Refs. <cit.>) but not previously developed into an algorithm.
    
    Let us start by describing the algorithm to solve the directed polymer problem, before explaining how to adapt it to address the Darcy problem.
    
    
    
    
    
     §.§ Energy-ordered tree generation
    
    
    Here, we explain how to sequentially construct the configurations P of increasing energy ε^( P)_0 T, where ε^( P)_t_1 t_2 was defined in Eq. (<ref>), starting from the ground state.
    
    The main idea is that the distribution of the ground state energy, min_ Pε_0 T^( P)≡ε_0, is determined by a non-linear evolution equation that is straightforward to solve numerically (at least when all energies are discrete). This allows us to determine the ground state energy without generating the entire random tree. Once ε_0 is known, a polymer configuration P_0 with that minimum energy can be exactly constructed using a branching random walk defined by two elementary processes with non-trivial but definite probabilities. Furthermore, we consider the first excited state as the configuration that minimizes the ground state energy among all polymer configurations branching off P_0. Its energy can also be determined without needing to generate all configurations a priori. This process can be iterated until all relevant configurations have been generated.
    
    We will successively discuss the probability distribution of the ground state energy, the construction of the ground state configuration(s), and then the excited states.
    
    
    
    
    
      
    Ground state energy.
    
    We denote by _T(ε_0>ε) the probability that the energy of the ground state for a tree of height T is greater than ε. We call p(τ) the distribution of the random energies τ on the bonds. Then
    _T(ε_0>ε)=∑_τ_01p(τ_01)_T-1(ε^bin_0>ε-τ_01) ,
    
    where _s(ε_0^bin>ε') stands for the probability that the energy of the ground state of a perfect binary tree of height s is larger than ε'. The discrete sum becomes an integral if p is a probability density. But for practical reasons, we will assume that τ takes only integer values. Other cases may be recovered through an appropriate scaling, and possible limiting procedures.
    
    The distribution of ε_0^bin obeys a recursion relation in the height of the binary tree:
    _s+1(ε_0^bin>ε')=(∑_τ p(τ) _s(ε_0^bin>ε'-τ))^2 .
    
    The ground state energy of a polymer made of a single monomer is of course null, and thus the initial condition simply reads
    _s=0(ε_0^bin>ε')= 1_{ε'<0} .
    
    The derivation of recursions such as Eq. (<ref>) is standard (see e.g. <cit.> for the discussion of a very similar model).
    
    It proves convenient to introduce the probability that the ground state energy ε_0^bin is not larger than some ε'. It is just the complementary of the previous probability:
    
    u_s(ε')≡_s(ε_0^bin≤ε')=1-_s(ε_0^bin>ε').
    
    From Eq. (<ref>), we see that it obeys the non-linear finite-difference equation
    
    u_s+1(ε')=2∑_τ p(τ)u_s(ε'-τ)-(∑_τ p(τ)u_s(ε'-τ))^2,  with u_s=0(ε')= 1_{ε'≥ 0}.
    
    Such equations can be shown to admit “pulled-front solutions”, and in this sense, belong to the universality class of the Fisher-Kolmogorov-Petrovsky-Piscounov (FKPP) equation <cit.> (see Ref. <cit.> for an extensive review, as well as Appendix <ref>).
    
    From Eqs. (<ref>) and (<ref>), the distribution of ε_0 is found to read
    _T(ε_0=ε)=∑_τ_01p(τ_01)
    [u_T-1(ε-τ_01)-u_T-1(ε-τ_01-1)].
    
    
      
    Building one configuration of minimal energy.
    
    Once the energy τ_01 of the first bond and the minimum energy ε_0 are known, the paths P_0≡ P_0 [0,T] in the tree corresponding to the minimal-energy configurations can be constructed iteratively, using a biased branching random walk that we shall now fully specify.
    
    We shall assume that the sub-path P_0 [0,t] of a configuration of minimum energy (which is not necessarily unique) is known. Given this, we can determine the law governing the next bond energy τ_t t+1. The sub-path P_0 [t,T], which we seek to construct, corresponds to a minimal-energy configuration of a polymer with T-t+1 monomers. This sub-path exists within the configuration space represented by the binary tree rooted in the node of P_0 [0,t] at depth t. The energy of this sub-path is fixed to ε^( P_0)_t T= ε_0-ε^( P_0)_0 t.
    
    We now proceed with the construction. The sub-polymer P_0 [0,t] has two possible continuation bonds, with energies denoted by τ^(0)_t t+1 and τ^(1)_t t+1, respectively (see Fig. <ref>). To determine their joint distribution, we need the probabilities R_t^ε_0(ε) and B_t^ε_0(ε) which describe the likelihood that a polymer configuration with a segment energy ε over [0,t] will have a minimum total energy either equal to ε_0 or greater than ε_0, respectively. These probabilities correspond to the likelihood that the ground state energy of sub-polymers rooted at (ε,t), and made of T-t+1 monomers, is equal to or greater than ε_0-ε. It is then enough to recall that u_s(ε') is the probability that the ground state energy of a binary tree of height s is not larger than ε'. The sought probabilities follow immediately:
    
    R_t^ε_0(ε)=u_T-t(ε_0-ε)-u_T-t(ε_0-ε-1)
      and  
    B_t^ε_0(ε)=1-u_T-t(ε_0-ε) .
    
    
    Let us now introduce the paths P_0^(0) and P_0^(1), defined such that they both overlap with P_0 on the segment [0,t] and then follow the two possible outgoing bonds from the node of P_0 located at depth t:
     P_0 [0,t]^(0)= P_0 [0,t]^(1)= P_0 [0,t],
      and  ε^( P_0^(0))_t t+1=τ_t t+1^(0),
     ε^( P_0^(1))_t t+1=τ_t t+1^(1).
    
    The labels are chosen so that P_0^(0) will eventually coincide with P_0. We must also impose the following conditions on the paths P_0^(0) and P_0^(1):
    min_ P_0 [t+1,T]^(0)ε^( P_0^(0))_0 T=ε_0,
      and  min_ P_0 [t+1,T]^(1)ε^( P_0^(1))_0 T≥ε_0.
    
    We are now ready to express the probability distribution for the pair (τ^(0)_t t+1,τ^(1)_t t+1).
    
    There are two cases to distinguish, depending on whether the inequality in Eq. (<ref>) holds as an equality, or as a strict inequality:
    
      *  Equality: the possible configurations P_0^(0) and P_0^(1) have the same minimum energy ε_0.  Then, the joint probability of the latter event and that the bond energies have the definite values τ, τ', given that the minimal energy of the overall polymer configuration P_0 is ε_0, reads
      
        (τ^(0)_t t+1=τ,
          τ^(1)_t t+1=τ';min_ P_0 [t+1,T]^(0)ε^( P_0^(0))_0 T=min_ P_0 [t+1,T]^(1)ε^( P_0^(1))_0 T
          =ε_0|min_ P_0[t,T]ε^( P_0)_0 T=ε_0)
        
          =p(τ)p(τ')R_t+1^ε_0(ε^( P_0)_0 t+τ)
            R_t+1^ε_0(ε^( P_0)_0 t+τ')/R_t^ε_0(ε^( P_0)_0 t).
    
      *  Strict inequality: one of the possible configurations, which we chose to label as P_0^(1), has ground state energy larger than ε_0. In this case,
      
        (τ^(0)_t t+1=τ,
            τ^(1)_t t+1=τ';min_ P_0 [t+1,T]^(0)ε^( P_0^(0))_0 T=ε_0,
            min_ P_0 [t+1,T]^(1)ε^( P_0^(1))_0 T
            >ε_0|min_ P_0[t,T]ε^( P_0)_0 T=ε_0)
        
            =2p(τ)p(τ')R_t+1^ε_0(ε^( P_0)_0 t+τ)
              B_t+1^ε_0(ε^( P_0)_0 t+τ')/R_t^ε_0(ε^( P_0)_0 t) .
    
    We check that the sum of these probabilities marginalized with respect to τ and τ' is unitary, as it should: it is enough to replace R and B appearing in the above expressions by their definitions (<ref>) in terms of u's and perform the sums, using the evolution equation (<ref>).
    
    We first determine if we are in case <ref> or <ref> using the probabilities (<ref>),(<ref>) summed over τ and τ'. In a second step, we draw τ and τ' from either the probability in Eq. (<ref>) or that in Eq. (<ref>), with appropriate normalization following from the conditioning to be in cases <ref> or <ref> respectively.
    
    In both cases, the outgoing bond labeled (0) extends the path P_0. If one ends up in the non-degenerate case <ref>, then once τ_t t+1^(1) is determined, one draws the minimal energy that can have the other polymer configuration: the next configuration to be constructed will be one among those realizing this minimum. The distribution of min_ P_0^(1)ε^( P_0^(1))_t+1 T reads
    (.min_ P_0 [t+1,T]^(1)ε^( P_0^(1))_t+1 T=ε|min_ P_0 [t+1,T]^(1)ε^( P_0^(1))_0 T>ε_0)
    =u_T-t-1(ε)-u_T-t-1(ε-1)/B^ε_0_t+1(ε^( P_0)_0 t+τ_t t+1^(1)) 1_{ε>ε_0-ε^( P_0)_0 t-τ_t t+1^(1)} .
    
    We shall call P_1t̂ a polymer configuration that realizes this minimum. Note that at this point, we don't need to have fully constructed this configuration: all we need for the next step is the value of its energy, and we have just shown that it can be drawn independently of knowledge of the detailed configuration.
    
    One iterates this construction of the segments of increasing depth until the first path P_0 is complete. Whether there are degeneracies or not, we shall always build the configuration labeled (0) first, leaving the construction of the configuration (1) for a next step that we shall now describe. We store τ^(1)_t t+1, and if there is no degeneracy, namely in case <ref>, we also store the values of the minimum energy of the polymer configurations branching off this node of P_0. 
    
    
    
    
    
    
      
    Next configurations.
    
    Once a ground state configuration has been constructed, we search for the next lowest energy level. There may be several configurations having the same ground state energy. In this case, we just pick randomly one of the nodes off which two sub-polymer configurations with degenerate energies branch, say at depth t, and we construct the second polymer configuration. The latter results from a branching random walk defined by the probabilities (<ref>) and (<ref>). We eventually get a path P_0', with overlap t with P_0. In the same way as we did for the path P_0, we draw and store the bond energies of the polymer configurations which branch off the newly-constructed nodes of P_0', as well as the minimum energies of these configurations.
    
    If there is no degeneracy of the ground-state energy ε_0, or after having opened all degenerate paths, we determine the energy of the first excited state. It must be the lowest energy of all configurations P_1t̂', where t' runs over all nodes of the path P_0 (or of the set of paths P_0, P'_0,⋯ of energy ε_0). In the non-degenerate case, in which the sub-tree of configurations of energy ε_0 consists of a single branch,
    ε_1=min_t'(ε_0 t'^( P_0)+ε_t' T^( P_1t̂')).
    
    The polymer configuration P_1 that branches off P_0 at the node labeled, say, t, and minimizes this energy is then constructed, again from a branching random walk. The latter is defined by the probabilities (<ref>) and (<ref>), up to the replacements P_0→ P_1 (as well as for the super-scripted paths with (0) and (1)) and ε_0→ε_1. A sketch of two configurations corresponding to the two lowest values of the energy in the non-degenerate case is displayed in Fig. <ref>.
    
    In the degenerate case, the minimum in Eq. (<ref>) is taken over all nodes of the sub-tree of the configurations of energy ε_0, and in the term ε_0 t'^( P_0), the path P_0 is replaced by a path of that tree that includes the considered node.
    
    The whole procedure is thus iterated until some stopping condition is reached, e.g. on the maximum energy or on the number of configurations.
    
    
    
    
    
     §.§ Threshold pressure-ordered tree generation
    
    
    The algorithm to open the successive channels as the pressure at the inlet is increased can almost be taken over from the one just described for the directed polymer problem. The only variation will turn out to be in the ordering of the channels to open.
    
    Let us start by recalling the dictionary between the directed polymer problem of finding the lowest-lying energy configurations, and the Darcy problem of opening channels when the pressure at the inlet of the network is dialed up. The bonds between monomers correspond to throats which join at the trivalent nodes of the tree-like network. The latter are identified to pores in the Darcy problem. The rank t of the monomers in the polymer is the depth of the throat joints, its size T is the height of the network. The bond energies τ are the elementary threshold pressure difference of the considered throats. 
    
    However, in the DP problem, the successive configurations that we generate are ordered in energy ε_0,ε_1,⋯. In the Darcy problem instead, we open the channels in increasing order of the applied pressure thresholds P_0, P_1,⋯. In general, the latter do not correspond to the successive energies ε_0 T^( P_0), ε_0 T^( P_1),⋯.
    
    To construct the sub-tree of open channels, we start with the first channel, which matches exactly a configuration that possesses the ground state energy in the DP case. We then search for the channels that open successively as the pressure is increased. The procedure to find the next channel, after P_0, to open is exactly the same as the one implemented in the naive algorithm in Sec. <ref>. Given the open sub-tree, instead of looking for the next minimum energy, we compute the threshold pressure at the inlet for each of the closed channels attaching to it, assuming that all other channels remain closed. Instead of the energy ε_1 given by Eq. (<ref>), the relevant variable P_1 is this minimum threshold pressure (<ref>). We construct and open one of the channels that possess the threshold pressure P_1. 
    
    We then iterate, using the more general equations (<ref>),(<ref>) in order to compute P_2,P_3⋯, until some stopping condition is satisfied.
    
    
    
    
    
      
    Performance of the new algorithm.
    
    The major improvement brought by our new method is obviously that we do not need to generate the full tree of throats (i.e. the random medium) a priori. We are able to generate the branches that correspond exactly to the channels that open, in successive order when the exerted pressure is increased. This is all the more useful as we only need a number of open channels that is much smaller than the total number 2^T-1 of possible channels.
    
    Consequently, the complexity of our new algorithm is linear in the height T of the tree-like network, instead of exponential as for the more straightforward algorithm described in Sec. <ref>. At each threshold, the calculation time is dominated by the time needed to search for the next channels to open. Hence for a total number of channels n^ch, the complexity is bounded by O(T×(n^ch)^2). This allows us to increase the height accessible to numerical investigations by several orders of magnitude.
    
    Note however that a solution u_t(ε') of Eq. (<ref>) needs to be pre-computed and stored for all t in [0,T], and for all relevant values of ε', the number of which becomes proportional to √(T) at large T. But in practice, this does not impose any major restriction on the calculations that are useful for the insight we wish to gain from a numerical study.
    
    All numerical data presented and discussed below were generated running our new algorithm on an off-the-shelf laptop computer, for not more than a few hours for each parameter setting. Note that we also implemented the naive algorithm, and performed accurate comparisons of the calculations of observables using both algorithms for low values of T, in order to validate our implementation.
    
    
    
    
    
    § ANALYSIS OF THE MODEL
    
    
    In this section, we apply our optimized algorithm to generate new relevant numerical data. Before presenting our results, we will briefly revisit a few conjectures proposed in Ref. <cit.> (with detailed derivations in Ref. <cit.>). These conjectures lead to asymptotic formulas that our data will allow us to test, thereby providing an opportunity to confirm or refute the underlying theoretical framework.
    
    
    
     §.§ Analytical formulas in two asymptotic limits
    
    
    According to Ref. <cit.>, there are essentially two limits in which the flow can be worked out: (i) when the pressure P exerted at the inlet is large enough so that all channels are open, (ii) for fixed P, but very large tree height T. Let us consider these limits in turn.
    
    
    
    
    
      
    High inlet pressure.
    
    We first assume that P is set above the threshold for the opening of the last of the 2^T-1 channels. Then the tree of height T is complete, as well as all its sub-trees. Looking at equation (<ref>) for κ_eff, it is easy to convince oneself that the latter depends only on the depth of the node one considers. It depends neither on the τ's, nor on the particular node at a given depth.
    
    Let us call κ_t the value of κ_eff at the top of a (sub-)tree of generic height t, and P_t^* the value of the effective threshold pressure at the inlet of the same (sub-)tree. [Note that at variance with κ_t, P_t^* is a random variable also in this case in which all channels are open, meaning that it depends on the particular realization of the (sub-)tree]. We then apply Eq. (<ref>) to the node of lowest depth, namely at the top of the tree: κ_T=2κ_T-1/1+2κ_T-1  and  
    P_T^*=τ_01+1/2(P_T-1^*(0)+P_T-1^*(1)) .
    
    The equation for κ_T is a recursion. With the initial condition κ_1=1, it can be easily solved:
    κ_T=2^T-1/2^T-1 .
    
    The equation for P_T^* can also be turned into a recursion by taking the average over realizations of the disorder. The obtained equation for P^*_T can be trivial solved, given the initial condition P^*_0=P^*_0=0:
    P^*_T=T×τ .
    
    Hence these simple calculations yield the following useful asymptotical expressions:
    κ_eff→1/2  and  P^*_eff/T→τ  when P→∞.
    
    Note that P^*_eff is a self-averaging quantity, so we expect that at large pressure P and large tree height T, all the flow curves converge to 
    
        Q_T(P)=1/2(P-τ T).
    
    
      
    Great tree height.
    
    When the pressure at the inlet is finite, a finite number of channels are open. Of course, the channels that are most likely to open when the pressure at the inlet increases are expected to sit at low ε_0 T. But it is also advantageous that they have a small overlap with each other, namely that the paths get quickly separated. This is easily seen from the minimization of P_thr 0 in Eq. (<ref>), realized in Eq. (<ref>) by the paths P_0 and P_1:
    
    P_1=P_0+T/T-t_1(ε_0 T^( P_1)-ε_0 T^( P_0)).
    
    To get the threshold pressure P_1 as close as possible to P_0, it is apparent that the path sum ε_0 T^( P_1) of the newly-opened channel needs to be close to that of the first channel ε_0 T^( P_0), and at the same time, the overlap t_1 of the two channels needs to be as small as possible. Consequently, we may think that the ε's of the open channels are essentially independent identically-distributed random variables, like the energy levels in the Random Energy Model (REM) <cit.> (see Appendix <ref>). 
    
    From this observation, one predicts a definite form for the distribution of the number of open channels in a given pressure range above P_0, see Eq. (<ref>), that can be checked with the help of our numerical code. The mean flow for a pressure P such that P-P_0=x, where x is kept fixed, can also be evaluated in this simplified picture <cit.>: 
    Q_T(P_0+x)=e^β_c x-1/β_c T .
    
    
    Let us give an alternative proof of this formula to that in Ref. <cit.>. For pressures different from the opening pressures of the various channels, the flow Q is the affine function of the pressure as described in Eq. (<ref>). Hence its derivative is simply the piece-wise constant function κ_eff:
    dQ_T(P_0+x)/dx=κ_eff(P_0+x).
    
    On the other hand, it is easy to see, iterating Eq. (<ref>), that for a deep tree made of n^ch branches which have negligible overlaps:
    κ_eff=n^ch/T.
    
    Let us now replace the right-hand side of Eq. (<ref>) by Eq. (<ref>), and take the expectation value:
    dQ_T(P_0+x)/dx=n^ch(x)/T.
    
    Consistently with our picture, we identify n^ch with n_REM provided in Eq. (<ref>). Integrating Eq. (<ref>) over x, taking account of the boundary condition Q_T(P_0)=0, one easily arrives at Eq. (<ref>). The latter is a prediction that we will be able to check numerically.
    
    
    
    
    
    
     §.§ Calculation of selected observables in a numerical implementation of the improved algorithm
    
    
    The numerical model is fully defined once the distribution p(τ) of the elementary thresholds τ is specified. We expect the asymptotic behavior of the observables we consider to be universal across a broad class of models. The qualitative behavior should not significantly depend on the specific choice of p(τ), aside from a small number of easily calculable constants, as long as the distribution is sufficiently regular and decreases “rapidly enough” for large values of |τ|.
    
    For simplicity, we take the uniform distribution on a set of the contiguous non-negative integers {1,2,⋯,N}. The parameter N was set to 20 for the flow calculations presented in Sec. <ref>, as these calculations are relatively insensitive to N. We increased N to 200 for the observables discussed in Sec. <ref> to minimize discretization artefacts, which these observables are more sensitive to. There would be no problem in principle in taking a continuous distribution p(τ), but the implementation would be significantly more complicated, technically speaking. The tricky part would be solving Eq. (<ref>) and storing the result for all values of the depth between 0 and T.
    
    The few model-dependent parameters needed as inputs to the analytical formulas are computed in Appendix <ref>.
    
    
    
    
    
      §.§.§ Flow
    
    
    As argued in Sec. <ref>, the flow for a given realization of the network is a piece-wise affine function of the pressure P applied at the inlet, described by the effective Darcy law (<ref>). In order to visualize Q_T(P), we generated 10 realizations of it for different values of T; see Fig. <ref>. We see that the realizations all exhibit the same global shape, essentially up to a T-dependent random shift. The main T-dependence must come from the random threshold pressure P_0, the expectation value of which roughly grows proportionally to T. We can also see clearly that at large P, namely for large number of open channels, the slope κ_eff of Q_T(P) for each of the realizations becomes close to the expected asymptotics given in Eq. (<ref>). It turns out that P_eff^* instead stays lower than the expected value T×τ, where τ=10.5 in our choice for the distribution p(τ). We shall discuss this fact further below. In Fig. <ref> is displayed the average of the flow Q_T(P_0+x) over many realizations, as a function of x. The curves that correspond to different values of T are very similar at least for large enough x, and almost superimpose, now that the realizations have been shifted by the threshold pressure -P_0. The inset shows that for large tree heights, the scaled mean flow T×Q_T(P_0+x) may converge to some scaling curve. We also see that the cluster of curves all tend to the function x↦ x at the origin x→ 0, consistently with the expression (<ref>) of Q_T just above the threshold for the opening of the network.
    
    
    
    Let us focus on this mean scaled flow for moderate values of x. We perform the calculation with the values of T used in Ref. <cit.> (i.e. T≤ 21), as well as for much larger T: for this observable, we were able to generate data for values of T that are about three orders of magnitude larger than in the latter reference. The result of our calculation is shown in Fig. <ref>. With the new data, we see a clear convergence to the expected asymptotic expression (<ref>), which was not possible to see with lower values of T: indeed, this convergence turns out to be very slow with T.
    
    
    
    
    
      
    Effective parameters.
    
    In order to investigate the properties of the flow in more details, we compute the permeability κ_eff and the threshold pressure P_eff^* entering the effective Darcy law (<ref>) as a function of the number of open channels. We then average these quantities over the realizations. 
    
    
    
    Our numerical data for κ_eff is displayed in Fig. <ref>. We see that the asymptotics (<ref>) are clearly approached when T assumes moderate values. For larger values of T, the expected convergence is slower. For the largest values of T we set (T=1000), the limit is hardly approached even when as many as n^ch=1000 channels are open.
    
    In order to quantify the convergence, we compute the expectation value of the minimum number of channels n^ch_SAT needed to reach an effective permeability κ_eff≥ 0.4; see Fig. <ref>. We find that the T-dependence of this quantity is linear to a good approximation when T becomes large. This clearly confirms the observation in Ref. <cit.>. The latter was inferred from data generated with the help of the equivalent of the basic algorithm described in Sec. <ref>, although we now see that these data were very far from any asymptotic regime.
    
    
    
    It is also useful to plot κ_eff as a function of the scaled variable n^ch/T; see Fig. <ref>. It seems clear that at large T, one tends to a limiting curve. The fact that the cluster of curves exhibits some significant dispersion is related to the slow convergence with T. The shape of the limiting curve for low and high numbers of open channels is consistent with the asymptotic form we may expect: we approach Eq. (<ref>) in the regime n^ch≪ T, and Eq. (<ref>) when n^ch≫ T. The transition between these two regimes is confirmed to happen when the number of open channels is on the order of T <cit.>.
    
    The mean threshold pressure P_eff^* is displayed in Fig. <ref>, as a function of the number of open channels n^ch and for different values of T. When n^ch is large, we expect that P_eff^*/T tends to τ. While the numerical results seem in agreement with this expectation for the lowest values of T, we see that when T gets large while keeping n^ch fixed, this quantity seems to converge instead to -v(β_c). Actually, this is the value that would be expected when only one channel is open, namely the ground state configuration in the DP problem, and in the limit T→∞. Indeed, according to Eq. (<ref>) and Eq. (<ref>),
    P_eff^*(n^ch=1)/T=1/Tε^( P_0)_0 T=-v(β_c)+3/2β_cln T/T+ O(1/T) .
    
    
      §.§.§ Statistics of open channels
    
    
    So far, we have focused on the flow, the most natural observable for this problem. Now we focus on n^ch(x), namely the number of open channels at the pressure P= P_0+x. Note that here we fix the value of x for all disorder realization while the threshold P_0 and thus the pressure P are realization dependent. In particular we are interested in its average n^ch(x) and its distribution. 
    
    In Fig. <ref> is displayed the histogram of n^ch(x) for  x=20. We see that the larger T, the closer one gets from what would be expected for the Random Energy Model (REM), namely
    
    P_n^ch(x)=e^-β_c x(1-e^-β_c x)^n^ch-1 ,
    
    consistently with the simplified picture (see Appendix <ref>). We have set N=200, instead of N=20 used for the flow calculation: in this way, the degeneracies  of the energy levels due to finite N are strongly suppressed and we are able to recover the asymptotic large T  behavior given in  Eqs. (<ref>),(<ref>).
    
    Figure <ref> shows the mean, n^ch(x), divided by what the expected  leading behavior, namely e^β_c x. We see that the larger T,  the curves get constant  and closer to 1. Again, the convergence is quite slow. We deemed it useful to compare the mean number of open channels to the mean number of configurations of the DP in the energy interval x above the ground state. We see that the latter also converges to the expected asymptotics ∝ x  e^β_c x.
    
    Note that the small  oscillations seen in the curves are an effect of finite  N and thus of the  fact that the energy levels are discrete. Their  amplitude increases for a smaller N. However, we have verified that also for N=20 used for the flow the trend of the curves is the same.
    
    
    
    § CONCLUSION
    
    
    We have designed and implemented a novel algorithm that exactly generates the open channels of a tree-like network of throats at a specified pressure. This method is effectively a spinal decomposition of the tree, with the primary spine corresponding to the first open channel. Leveraging this algorithm, we were able to analyze networks of unprecedented depth, achieving a two to three orders of magnitude increase in tree height T compared to previous studies, all while using a standard laptop. Given the algorithm's linear complexity with respect to T, even greater heights are within reach.
    
    
    The results obtained from our algorithm confirm with remarkable precision the conjectures proposed in Ref. <cit.>. Specifically, we have successfully verified the analytical expression for the mean flow [Eq.(<ref>)], demonstrating the feasibility of deriving an explicit expression for the non-linear Darcy law of a yield stress fluid within this network. Moreover, our findings strongly support the physical picture argued in Ref. <cit.>:
    
      *   The first channels that open are those with low opening thresholds and minimal overlap. As a consequence, the number of open channels, n^ch(x), identifies to the number of low-energy levels of the Random Energy Model (REM).
        
        
      *  The effective permeability, κ_eff, initially increases linearly with the number of open channels. Upon the activation of approximately t channels, κ_eff saturates to the Newtonian permeability of the medium. It is noteworthy that t channels represent a very small fraction compared to the total possible channels in the network, 2^t-1.
        
        
      *   The effective offset, P_eff^*, gradually evolves from the ground state energy of the associated directed polymer to the mean energy of a non-optimized polymer.
    
    
    These data and our algorithm open up the prospect of studying sub-asymptotic corrections, and to elaborate a more refined picture of the corresponding physical effects, beyond the simplified REM model for the tree-like network. These results provide a robust validation of the theoretical framework and offer significant insights into the behavior of yield stress fluids in complex networks.
    
    A longer-term and more ambitious project would be to understand if the behavior of the effective permeability and effective  offset is peculiar to tree-like networks or if it generalizes to networks of finite dimension <cit.>.
    
    
    
    
    
    
    § ACKNOWLEDGMENTS
    
    
    S.M. thanks Pascal Maillard, Bastien Mallein, Michel Pain for their interest and for related discussions. A.R. thanks Chen Liu, A. De Luca, V. Schimmenti and L. Talon for the useful discussions.  This project was initiated in the framework of the “GDR Branchement”. It has received financial support from the CNRS through the MITI interdisciplinary programs, and from the Agence Nationale de la Recherche (ANR), grant ANR-23-CE30-0031-04 (DISCREEP).
    
    
    
    
    
    § TRAVELING WAVE SOLUTIONS TO THE NON-LINEAR EQUATION AND ENERGY LEVELS OF THE DIRECTED POLYMER NEAR THE GROUND STATE
    
    
    From Eq. (<ref>) one derives the probability distribution of the ground state of the directed polymer. Let us discuss the solutions to this equation.
    
    It is well-known that for a wide class of initial conditions, including those relevant to this paper, the large-time solution of the non-linear equation (<ref>) forms a universal traveling wave (see e.g. Ref. <cit.>, and Ref. <cit.> for the mathematical proof in the case of the FKPP equation). The latter refers to a uniformly-moving front that monotonously connects u_s=0 (for ε'→-∞) to u_s=1 (for ε'→+∞). The standard way of determining the properties of the traveling wave is to solve the equation obtained from the linearization of Eq. (<ref>) in the region where u_s is small, where only the first term in the right-hand side is relevant. A suitable Ansatz is
    
    u_s(ε')∝ e^β[ε'+v(β)s],
    
    where β is a positive parameter. Requiring that this form be a solution of the linearized equation (<ref>) determines v(β):
    
    v(β)=1/βln(2∑_τ p(τ)e^-βτ).
    
    (See e.g. <cit.> for a discussion of a very similar model). For the initial condition in Eq. (<ref>), the shape of the asymptotic traveling wave in the region where u_s(ε')≪ 1 is known to be of the form (<ref>) with β=β_c, where β_c minimizes v(β). The velocity of the wave is then -v(β_c).
    
    Corrections for large but finite network heights s are also known. The front position reads <cit.>
    m_s=-v(β_c)s+3/2β_cln s+const+o(1) .
    
    This quantity also represents the expectation value (or, up to a constant, the median value) of the ground-state energy ε_0^bin of a directed polymer on a binary tree of height s. For ε'<m_s and in the limit m_s-ε'≫ 1, the traveling wave has the following shape:
    
    u_s(ε')≃const×(m_s-ε'+const') e^β_c(ε'-m_s)× e^-(ε'-m_s)^2/[2β_c v”(β_c)s] .
    
    In the region m_s-ε'≪√(β_c v”(β_c)s), the dominant ε'-dependence is indeed the exponential form (<ref>), with β=β_c, corrected by a linear factor. For m_s-ε'≳√(β_c v”(β_c)s), the Gaussian factor drives u_s(P') rapidly to very small values.
    
    A derivation of these results can be found e.g. in Ref. <cit.>. It is not difficult to recover them heuristically: it is sufficient to recognize that, at this level of accuracy and for the particular quantities considered, the non-linear evolution equation can be replaced by its linearization, supplemented with a moving absorptive boundary. We refer the reader to the established literature on this topic, e.g. Ref. <cit.> and references therein.
    
    
    
    
    
      
    Mean density of energy levels.
    
    Analyzing further the relationship between the FKPP equation and branching Brownian motion, or between the corresponding finite-difference equation and its associated branching random walk respectively, one may prove that the full information about the distribution of the energy levels near the ground state can be deduced from the solution to Eq. (<ref>), if one takes appropriate initial conditions <cit.>. Of particular interest for the present work is the expectation value n_BRW(x) of the number of energy levels at a distance x above the ground-state energy ε_0^bin <cit.>:
    n_BRW(x)≃const× x  e^β_c x .
    
    This formula holds for (formally) infinite-s, and x≫ 1/β_c.
    
    
    
    
    
      
    Choice of p(τ) in our numerical calculations.
    
    In our particular implementation, we chose p(τ) uniform on the first N non-zero integers, namely
    
    p(τ)=1/N   for τ∈{1,2⋯ N},
    
    0    else.
    
    Concretely, for the two values of N we have used for our numerical simulations, the relevant parameters are given in Tab. <ref>.
    
    
    
    § ENERGY LEVELS IN A RANDOM ENERGY MODEL
    
    
    Let us assume that our configuration space consists of 2^T-1 independent branches, and consider a branch labeled P. It is characterized by the random energy ε_0 T^( P). The generating function G(γ) of the probability (ε_0 T^( P)=ε) that this energy equals some given ε reads
    
    G(γ)≡∑_ε e^-γε (ε_0 T^( P)=ε)
    =∑_ε e^-γε[∑_τ'_0 1,⋯,τ'_T-1 T∈ P(∏_j=0^T-1 p(τ'_j j+1))δ_ε,∑_j=0^T-1τ'_j j+1],
    
    and is obviously the same for all branches. An easy calculation yields
    
    G(γ)=(∑_τ'p(τ')e^-γτ')^T.
    
    The sought probability is then obtained from an appropriate integration:
    (ε_0 T^( P)=ε)=∫_ Cdγ/2iπ e^γεG(γ)=∫_ Cdγ/2iπ e^γε(∑_τ'p(τ')e^-γτ')^T ,
    
    where the integration domain can be, for example, the segment C≡[γ_0-iπ,γ_0+iπ] in the complex plane, for any real number γ_0. In the limit of very large T, the integral may be evaluated by a steepest-descent method. We introduce
    χ(γ)≡ln(2∑_τ'p(τ')e^-γτ')
    
    in terms of which the probabiity (<ref>) can be conveniently re-expressed:
    (ε_0 T^( P)=ε)=∫_ Cdγ/2iπ e^γε+[χ(γ)-ln 2]T.
    
    We denote by γ_c the solution to the saddle-point equation
    ε+χ'(γ_c)T=0 .
    
    We then set γ_0=γ_c, and observe that the integration domain C may be extended to an infinite line parallel to the imaginary axis in the complex plane: the contribution to the integral of the added pieces of contour would be small for large T. The steepest-descent method evaluation of (ε_0 T^( P)=ε) yields
    (ε_0 T^( P)=ε)T≫ 1∼ e^γ_cε+[χ(γ_c)-ln 2]T .
    
    From Eq. (<ref>), we see that (ε_0 T^( P)=ε) decreases exponentially as ε decreases. Thus the typical ground state energy ε_0 is found by requiring that the probability that there is no state of that energy among the 2^T-1 independent branches is of order unity (for example 1/2; the precise value is actually irrelevant). This probability reads
    [1-(ε_0 T^( P)=ε_0)]^2^T-1 .
    
    Asking for this quantity to be of order unity for T large is equivalent to requiring that 2^T (ε_0 T^( P)=ε_0)= O(1). Taking the logarithm and using the saddle-point evaluation of (ε_0 T^( P)=ε_0) in Eq. (<ref>), neglecting terms that are small compared to T, we find
    γ_cε_0+χ(γ_c)T=0 .
    
    Combining Eq. (<ref>) and Eq. (<ref>), up to the replacement ε→ε_0 in the latter, we see that necessarily 
    γ_cχ'(γ_c)=χ(γ_c).
     
    Comparing χ defined in Eq. (<ref>) and v introduced in Appendix <ref> [Eq. (<ref>)], we observe that χ(γ)≡γ v(γ). Hence, Eq. (<ref>) is equivalent to v'(γ_c)=0, which implies that γ_c=β_c. The mean density of levels of energy ε near the typical ground state energy ε_0≃ -v(β_c)T then reads
    
    2^T-1(ε_0 T^( P)=ε)≃ e^β_c(ε-ε_0).
    
    
    If T is large enough, we can consider that the distribution of the energy levels is Poissonian <cit.>, the rate parameter of the Poisson law being the exponential in the right-hand side of Eq. (<ref>).
    
    One may compute the probability P_n(x) of having n energy levels in the interval of size x above the ground state energy, namely in [ε_0,ε_0+x]. When β_c is small, which is verified for large N, then x can be considered a continuous variable, since the rate parameter of the Poisson process varies slowly with x. The latter becomes the intensity of a Poisson point-like process on the line. In this continuous limit, we get a simple expression for P_n(x):
    
    P_n(x)=e^-β_c x(1-e^-β_c x)^n-1 .
    
    The mean number of levels which follows from this probability reads
    n_REM(x)=e^β_c x.
    
    Comparing this formula to the equivalent one for branching random walks (<ref>), we see that they differ essentially by a linear factor x.
    
    
    
    
    
    
    10Darcy:1856
    H. Darcy,  Les fontaines publiques de la ville de Dijon: Exposition et
      application des principes à suivre et des formules à employer dans les
      questions de distribution d'eau.
     Victor Dalmont, Libraire des Corps imperiaux des ponts et
      chaussées et des mines, 1856.
    
    Bear:1988
    J. Bear,  Dynamics of Fluids in Porous Media.
     Dover Civil and Mechanical Engineering Series, Dover, 1988.
    
    Sahimi:2011
    M. Sahimi,  Flow and transport in porous media and fractured rock: from
      classical methods to modern approaches.
     John Wiley & Sons, 2011.
    
    Blunt:2017
    M. J. Blunt,  Multiphase Flow in Permeable Media: A Pore-Scale
      Perspective.
     Cambridge University Press, 2017.
    
    Feder:2022
    J. Feder, E. G. Flekkøy, and A. Hansen,  Physics of Flow in Porous Media.
     Cambridge University Press, 2022.
    
    Roux:1987
    S. Roux and H. J. Herrmann, “Disorder-induced nonlinear conductivity,”
      Europhysics Letters, vol. 4, p. 1227, dec 1987.
    
    Liu:2019
    C. Liu, A. De Luca, A. Rosso, and L. Talon, “Darcy's law for yield stress
      fluids,” Phys. Rev. Lett., vol. 122, p. 245502, Jun 2019.
    
    Barnes:2000
    H. Barnes,  A Handbook of Elementary Rheology.
     Raymond F. Boyer Library Collection, University of Wales, Institute
      of Non-Newtonian Fluid Mechanics, 2000.
    
    Poiseuille:1840
    J.-L.-M. P. Poiseuille, “Recherches expérimentales sur le mouvement des
      liquides dans les tubes de très petits diamètre,” Comptes rendus des
      séances de l’Académie des Sciences, vol. 11, pp. 1041–1048, 1840.
    
    Schimmenti:2023
    V. M. Schimmenti, F. Lanza, A. Hansen, S. Franz, A. Rosso, L. Talon, and
      A. De Luca, “Darcy's law of yield stress fluids on a treelike network,” Phys. Rev. E, vol. 108, p. L023102, Aug 2023.
    
    Schimmenti:2023_supplemental
    V. M. Schimmenti, F. Lanza, A. Hansen, S. Franz, A. Rosso, L. Talon, and
      A. De Luca, “The darcy law of yield stress fluids on a tree-like network.” Supplemental material to Ref. <cit.>, available at
      http://link.aps.org/supplemental/10.1103/PhysRevE.108.L023102.
    
    Mezard:1986
    M. Mézard, G. Parisi, and M. Virasoro,  Spin Glass Theory and Beyond.
     WORLD SCIENTIFIC, 1986.
    
    Derrida:1988
    B. Derrida and H. Spohn, “Polymers on disordered trees, spin glasses, and
      traveling waves,” Journal of Statistical Physics, vol. 51,
      pp. 817–840, 1988.
    
    Brunet:2020yiq
    E. Brunet, A. D. Le, A. H. Mueller, and S. Munier, “How to generate the tip
      of branching random walks evolved to large times,” EPL, vol. 131,
      no. 4, p. 40002, 2020.
    
    Kirchhoff:1868
    G. Kirchhoff, “Ueber den Einfluss der Wärmeleitung in einem Gase auf die
      Schallbewegung,” Annalen der Physik, vol. 210, no. 6, pp. 177–193,
      1868.
    
    Angelopoulou:2023qdm
    A.-K. Angelopoulou, A. D. Le, and S. Munier, “Scattering from an external
      field in quantum chromodynamics at high energies: from foundations to
      interdisciplinary connections.” arXiv:2311.14796.
    
    Hardy:2009
    R. Hardy and S. C. Harris, “A SpineApproach to BranchingDiffusions
      with Applications to L^p-Convergence of Martingales,” in 
      Séminaire de Probabilités XLII (C. Donati-Martin, M. Émery,
      A. Rouault, and C. Stricker, eds.), pp. 281–330, Berlin, Heidelberg:
      Springer Berlin Heidelberg, 2009.
    
    Harris:2017
    S. C. Harris and M. I. Roberts, “The many-to-few lemma and multiple
      spines,” Annales de l'Institut Henri Poincaré, Probabilités et
      Statistiques, vol. 53, no. 1, pp. 226 – 242, 2017.
    
    Derrida:2016
    B. Derrida and P. Mottishaw, “On the genealogy of branching random walks and
      of directed polymers,” Europhysics Letters, vol. 115, p. 40005, sep
      2016.
    
    Fisher:1937
    R. A. Fisher, “The wave of advance of advantageous genes,” Annals of
      Eugenics, vol. 7, pp. 355–369, Jun 1937.
    
    KPP:1937
    A. Kolmogorov, I. Petrovsky, and N. Piscounov, “Étude de l'équation de la
      diffusion avec croissance de la quantité de matière et son application
      à un problème biologique,” Bull. Univ. État Moscou, A, vol. 1,
      no. 6, pp. 1–25, 1937.
    
    VanSaarloos:2003
    W. Van Saarloos, “Front propagation into unstable states,” Physics
      reports, vol. 386, no. 2-6, pp. 29–222, 2003.
    
    Derrida:1980
    B. Derrida, “Random-energy model: An exactly solvable model of disordered
      systems,” Physical Review B, vol. 24, no. 5, pp. 2613–2626, 1981.
    
    Bramson:1983
    M. D. Bramson, “Convergence of solutions of the Kolmogorov equation to
      travelling waves,” Memoirs of the American Mathematical Society,
      vol. 44, Jul 1983.
    
    Brunet:1997
    E. Brunet and B. Derrida, “Shift in the velocity of a front due to a cutoff,” Phys. Rev. E, vol. 56, pp. 2597–2604, Sep 1997.
    
    Brunet:2011
    E. Brunet and B. Derrida, “A branching random walk seen from the tip,”
      Journal of Statistical Physics, vol. 143, p. 420–446, Apr. 2011.
    
    Aidekon:2013
    E. Aïdekon, J. Berestycki, E. Brunet, and Z. Shi, “Branching brownian
      motion seen from its tip,” Probability Theory and Related Fields,
      vol. 157, no. 1-2, pp. 405–451, 2013.
    
    Mueller:2019ror
    A. H. Mueller and S. Munier, “Particle-number distribution in large
      fluctuations at the tip of branching random walks,” Phys. Rev. E,
      vol. 102, no. 2, p. 022104, 2020.
    
    Derrida:2015
    B. Derrida and P. Mottishaw, “Finite size corrections in the random energy
      model and the replica approach,” Journal of Statistical Mechanics:
      Theory and Experiment, vol. 2015, p. P01021, jan 2015.
