Vincenzo Maria Schimmenti

Federico Lanza

Alex Hansen

Silvio Franz

Alberto Rosso

Laurent Talon

Andrea De Luca

Darcy's law of yield stress fluids on a treelike network

Understanding the flow of yield stress fluids in porous media is a major challenge. In particular, experiments and extensive numerical simulations report a non-linear Darcy's law as a function of the pressure gradient. In this letter, we consider a tree-like porous structure for which the problem of the flow can be resolved exactly thanks to a mapping with the directed polymer (DP) with disordered bond energies on the Cayley tree. Our results confirm the non-linear behavior of the flow and expresses its full pressure-dependence via the density of low-energy states of DP restricted to vanishing overlap. These universal predictions are confirmed by extensive numerical simulations.

In a series of experiments during the nineteenth century, Henry Darcy studied the flow of water in a cylinder filled with sand [1] and established the empirical law for the flow rate Q as a function of the pressure difference P between the two ends of the cylinder

Q = κR 2 P/(ηL) (1)
here R and L are respectively the radius and the length of the cylinder, η is the viscosity of the fluid and κ is the permeability of the medium. Darcy gave an interpretation to the permeability κ starting from the Poiseuille's law which holds for empty cylinders and predicts Q Pois = πR 4 P/(8ηL). He assumed that, in a medium, the flow is possible only along non-intersecting thin channels, each of radius R c R, and the permeability can be identified as κ ∼ n ch R 2 c , with n ch the number of open channels per unit surface [2].

This empirical law is valid not only for water flowing in sand but for all Newtonian fluids [3] embedded in a porous medium (namely a complex structure with strong heterogeneities such as soils, rocks or sand [4][5][6][7]). Indeed for such fluids n ch is pressure-independent, as in each channel flow occurs for arbitrary weak pressure. This is not the case for another class of fluids, such as suspensions [8], gels [9], heavy oil [10], slurries or cement [11] for which a minimum yield stress σ Y is needed to flow [12]. For these fluids n ch grows with the applied pressure P . Experiments [13,14] and numerical simulations [15][16][17] indicated that the Darcy law is indeed modified: below a threshold pressure P 0 no flow occurs while, above it, the flow grows non-linearly with P . Three flowing regimes are observed [18,19] : i) initially the flow grows linearly in P -P 0 with a very small permeability ∼ 1/R 2 ; ii) for larger pressure the flow grows non-linearly as (P -P 0 ) β (with β ≈ 2) [20,21]. iii) above a saturation pressure P sat P 0 , the flow shows again a linear growth but with the much larger Newtonian permeability. Despite these detailed observations, a theoretical explanation for the non-linearity is still lacking. In this letter we propose both an explanation and provide an explicit prediction for the modified Darcy law. We consider a porous structure with the geometry of a binary Cayley tree with t levels and 2 t-1 channels (see Fig. 1). This geometry describes the flow in several biological networks (e.g. alveoli system in lungs [22] or leaf veins) and it is the simplest geometry with overlapping channels. Here we build on a method proposed in [23] to determine analytically n ch as a function of pressure. Our model captures the three regimes of the flow, in particular the non-linear one is described by Eq. ( 13)) and the saturation pressure reads P sat ∼ P 0 + cst • ln t . At this pressure there are ∼ t independent open channels and increasing further the pressure does not affect sensibly the permeability. The key point of our method is the identification of the first ∼ t open channels with the zero-overlap low energy states of the directed polymer on the Cayley tree, a model displaying one-step replica symmetry breaking (1-RSB) [24,25]. Those low-overlap excitations are abundant in mean-field glassy disordered systems but their number is suppressed in finite dimension and their role in realistic set-ups is controversial. However, in the Darcy problem, excitations with high overlap are strongly penalized as they are inessential in increasing the flow independently of the spacial dimensionality. For this reason, low-overlap excitations play an important role also in finite dimension for the flow of yield stress fluids in porous media.

The model. -Our model is a Cayley tree pore network filled by a Bingham fluid described by the modified Poiseuille's law

Q Pois (P ) = πR 4 8ηL (P -τ ) + (2)
Here we denote (x) + = max(0, x) and τ = Lσ Y /R and consider P > 0. In this model large open pores are con-

qα1α2 = q q0α1 = q0α2 0 α1 α2
FIG. 1. Top: Binary Cayley tree with t = 8 levels. Each path α has an energy α given by the sum of the thresholds along to the path. We assume a fixed pressure difference P between the top (root) and the leaves. Bottom: energies of the paths sorted as 0 < 1 < 2, . . . . At the pressure P0 = 0 the channel indicated in blue opens. Increasing the pressure, the orange channel opens at pressure P1 > P0 (see ( 5)), the green one at pressure P2 > P1. The pairwise overlap are indicated in figure with q the maximum among them. The energies corresponding to open channels are indicated in red in the axis at the bottom. In the large-t limit we show that q t and P1 = α 1 , P2 = α 2 , . . . .

nected by tubes with random radius and unit length so that the threshold τ fluctuates from tube to tube. For simplicity we set to 1 the prefactor of (2) and we consider the Bingham rhelogy between two pores i ad j as Q ij = (P ij -τ ij ) + with P ij the pressure difference between the pores and the yield threshold τ ij is drawn independently from a distribution p(τ ). We are interested in determining the global flow function Q t (P ; τ ) where τ represents the collective set of thresholds on each bond. Thank to the hierarchical structure, it is easy to obtain a recursive relation increasing the number of branches:

Q t (P ; τ ) = (P -P -τ 0 ) + , (3)
(P -P -τ 0 ) + = Q t-1 (P ; τ (1) ) + Q t-1 (P ; τ (2) ) . (4)

Eq. ( 3) comes from the constitutive relation applied on the topmost bond. Here P is the pressure on the node connected to the root and τ 0 the corresponding threshold. The second relation descends from the conservation of the flow where τ (i) , i = 1, 2 are the sets of thresholds restricted to the two subtrees with t -1 levels. Thus, assuming Q t-1 (P ; τ ) is known, one can use Eq. (4) to determine P = P (P ) and inject it in Eq. (3) to compute the flow function of the larger Cayley tree. Equations (3) and (4) are exact and amenable for numerical calculations, but their analytic treatment is possible only for vanishing thresholds, leading to the Newtonian behavior for P > 0, Q t (P ) = κ Newton t P , with the Newtonian permeability κ Newton t = 2 t-1 /(2 t -1). To obtain analytical expression of the flow in the disordered case one has to exploit the mapping with the directed polymer in the limit of large t.

Mapping to the directed polymer. -For each realisation of the thresholds τ , one has a sequence of pressures {P 0 , P 1 , P 2 , . . . , } at which a new channel opens. As originally observed in [21], this sequence of pressures is also related to the energy levels of an associated directed polymer (DP) problem. Indeed, consider the directed paths on the Cayley tree connecting the root to the leaves. Interpreting the threshold τ ij as the energy associated to the bond (ij), we can associate an energy α to each path α given by the sum of the thresholds belonging to the path. Ordering such energies as 0 < 1 < . . ., it is easy to verify that P 0 = 0 (see figure 1). Moreover, since there is only one channel open, only one term in the right-hand side of Eq. ( 4) is non-vanishing and for P ∈ [P 0 , P 1 ], one gets simply Q t (P ; τ ) = 1 t (P -P 0 ) + . The computation of the flow for P > P 0 has a twofold difficulty: first, one needs to characterise the sequence of pressures P 1 , P 2 , . . . at which new channels open; second, for a given geometry of open channels, one needs to evaluate the effective permeability κ ≡ dQ t /dP . Indeed for the subsequent channels, the relation between P a involves not only the energies but also the overlaps (i.e. the number of common bonds) of the new channel with the previously open channels. For instance to open two channels one has to solve the following minimization

P 1 = 0 + min α =0 α -0 1 -q0α /t = 0 + α1 -0 1 -q0α1 /t (5)
where α 1 labels the path realising the minimum and q0α stands for the overlap between the α-th channel and the ground state. It is crucial to remark that in the minimization, among the low-energy excited states, the ones with high overlap with the ground state are strongly penalized by the factor 1/(1 -q0α /t). The same property holds for the higher pressures P 2 , P 3 , . . . which display similar but more cumbersome dependence (see section A of [26]). Large-t limit. -A major simplification occurs in the limit t → ∞ as it is known that the directed polymers on Cayley tree display one-step replica symmetry breaking (1-RSB). For these systems the scaled overlap q 0α = q0α /t converges, at large t, to a bimodal distribution, with either q → 0 or q → 1 [24] (finite t corrections are also known [27,28]) implying q0α ∼ 1 or q0α ∼ t. Hence, in this large-t limit, among the two kinds of low-energy excited states, the open channels are of low-overlap type and basically uncorrelated. From this major simplification two other important properties arise: First, the opening pressures P 1 , P 2 , . . . coincide with the low energy levels α1 , α2 , . . . . Second, the total flow reduces to the sum of the contribution from each single open channel. As a consequence the computation of the flow problem reduces to determine the growth of the number of open channels n ch as a function of x = P -P 0 . This number identifies with the number of low-energy excited states of the directed polymer, provided they have low overlap among them. We compute this number adapting the tools introduced in [23] and based on the mapping to the discrete Kolmogorov Petrovsky Piskunov (KPP) equation [25,[29][30][31][32][33].

KPP equation. -To begin, we introduce the number of states with energy smaller than 0 + x, namely m

(full) t (x) = α ϑ(x -( α -0 )) (ϑ(x)
is the Heaviside theta function). The statistics of m (full) t (x) was studied in [23] and its large-t behavior can be computed by making use of the mapping to the discrete KPP equation (see also section C of [26] for a summary). An important result is that the t → ∞ limit exists finite, m

(full) ∞ (x) ≡ lim t→∞ m (full) t (x).
The same result holds for the random energy model (REM) [24,34], where the 2 t states are independent Gaussian random numbers of variance σ 2 t. In this case the average number of states is:

m (REM) ∞ (x) = e βcx ( 6
)
with β c = √ 2 ln 2/σ (the overline stands for the average over the disorder). For the DP a closed-form expression for the average number of states is not known [35]. Numerically it was shown in [23] that their number is larger than the corresponding REM model with the following asymptotics

m (full) ∞ (x) x 1 = Axe βcx ( 7
)
where A is an O(1) non-universal constant. Here we are interested in a subset of states around the minimum, namely the number of open channels below the pressure

P 0 + x, n ch t (x) = a ϑ(x -(P a -P 0 )). Since one al- ways has the upper bound m (full) t (x) > n ch t (x) the limit n ch ∞ (x) = lim t→∞ n ch t (x) exists finite.
It is tempting to identify their mean number n ch ∞ (x) with the number of states of the corresponding REM ∼ e βcx [24], as the zero overlap states are basically uncorrelated paths.

A direct numerical investigation of this claim is possible only for moderate t and it is carried out by the exact flow problem using an algorithm derived from [21] and explained in section B of [26]. Our results are plotted in Fig. 3 left and display a slow convergence to our claim, but with strong finite size effects. To make progress we introduce the quantity m q,t (x) which counts, in a tree with t levels, the number of paths with energy α ≤ P 0 +x and maximum overlap among themselves q = 1, . . . , t. In practice, this can be obtained by a pruning procedure (see Fig. 2): at the level q of the full Cayley tree, there are 2 q subtrees labelled by a = 1, . . . , 2 q . We replace each of these subtrees with a singe bond with an effective energy that equals (a) 0 , the minimal energy of the a-th subtree. These 2 q remaining states correspond to the ones with maximal overlap q. Clearly, for q = t, we FIG. 2. Left: An example of Cayley tree with a total of t = 5 generations and 2 4 leaves. The minimal path of the topmost subbranch after q = 2 generations is shown in bold. Right: Pruning of the full Cayley tree, where within each subtree starting from the q-th generation, only the minimal path is retained.

recover m (full) t (x). If, on the contrary, we take first the limit t → ∞, leading to m q (x) = lim t→∞ m q,t (x), we are effectively considering energies with vanishing scaled overlap q0α /t → 0. This remains true even if the limit q → ∞ is taken thereafter and we claim the following equality

n ch ∞ (x) = lim q→∞ m q (x) . ( 8
)
For t → ∞ and finite q, the pruning procedure is equivalent to growing a tree with q levels where the leaves thresholds are drawn from the distribution of of the minimum of an infinite tree. As a consequence, m q (x) can be computed by adapting the equations employed in [23] to compute m (full) t (x) as an integral of a function r t (x ; x):

m (full) t (x) = dx r t (x ; x) (9a) r t+1 (x ; x) = 2 dτ p(τ )Ω t (x -τ )r t (x -τ ; z) (9b) with initial conditions r 1 (x ; x) = p(x + x ) and Ω 1 (x) = -p(x).
Here -Ω t (x) is the distribution of the minimal energy of the directed polymer of t levels and it satisfies the discrete KPP equation:

Ω t+1 (x) = dτ p(τ )Ω t (x -τ ) 2 (10)
The recursive equations Eqs. ( 9) corresponds to growing a t + 1-level tree starting from two t-level trees. Thus, to compute m q (x), one has to adjust the initial conditions (see section C.3 of [26]), starting from a single bond with a threshold drawn from the distribution of the minimum of an infinite tree -w min (x) [36]. Here w min (x) solves the fixed point KPP equation where c(β c ) is the minimal value for which (11) has a solution (see section C.1 of [26]). Thus the equations for m q (x) are still equations ( 9) and (10) with the modified initial conditions Ω q=1 (x) = w min (x) and r q=1 (x ; x) = -w min (x + x ). In Fig. 3 middle we compute numerically both m (full) t (x) for large t and m q (x) for various q. Their difference is clear, moreover we see that by increasing q the solution approach the analytical prediction lim q→∞ m q (x) = e βcx (12) We refer to section D of [26] for an analytical derivation.

w min (x + c(β c )) = dτ p(τ )w 2 min (x -τ ) (11)
We are now in the position of evaluating explicitly the average flow at large t. For a fixed x = P -P 0 , there will be a finite number of open channels, sharing low overlap one with the other and each supporting a flow (x -x )/t, being P 0 + x its opening pressure. In this regime the total flow reduces the sum of the contribution from each single open channel, leading to

Q t (P 0 + x) t→∞ = x 0 dx t n ch ∞ (x -x ) = e βcx -1 β c t . ( 13
)
Equation ( 13) is the main analytical result of the work and it is valid for moderate pressure above P 0 , describing both the linear single-channel as well as the non-linear behaviours observed in experiments. Saturation pressure. -We now give an estimate for P sat . Let us first observe that from Eq. ( 13), we obtain an exponentially growth in the permeability κ t = dQ t /dP ∼ e βc(P -P0) /t, however its value should saturate to κ Newton t = 2 t-1 /(2 t -1). Equating the two regimes, one obtains an saturation already at P sat -P 0 ∼ ln t/β c . Let us comment about this result: when the pressure is slightly above the minimal value P 0 , only a single channel os is open and the permeability is ∼ 1/t. However, due to the fast growth of the density of available channels (12), a logarithmic increase in the pressure is enough to open ∼ t channels and saturate the permeability to its Newtonian value κ Newton t . A manifestation of this effect is shown in Fig. 3 right, where one observes the fast growth of the permeability followed by a negligible increase. Finally, in the inset of Fig. 3 right, we show the average number of channels n ch SAT required for the permeability to surpass a certain fraction δ • κ Netwon t and the linear growth in t is apparent.

Conclusions. -In this work, we have presented the study of yield stress fluids on the Cayley tree. We show that the problem is closely related to the directed polymer in the same geometry. In particular, in the limit of large trees, a direct mapping emerges between n ch and low-energy excited states with zero-overlap of the directed polymer. Thanks to this identification, we derive a simple universal expression for the flow as a function of the applied pressure. Equation ( 13) is independent of most microscopic details and threshold distribution and the tree branching ratio set only the parameter β c . The next big challenge would be to solve the problem of the flow in finite dimension. In particular it would be interesting to see if the low-overlap excitations of the associated directed polymer play a role also in the rheology of a finite dimensional porous medium.

Supplementary Material

The Darcy's law of yield stress fluids on a tree-like network

Few-channel solution

In this section we report the derivation of the expression for P 1 used in the main text and P 2 . We also write down the expression of the flow for one and two channels. You can check that in the limit t → ∞, P 1 = α1 , P 2 = α2 and tQ(P ) = (P -P 0 ) + (P -P 1 ) for P 0 < P < P 1 .

Flow for a single channel

At the critical pressure P 0 = 0 , corresponding to the ground state of the directed polymer, the first channels opens and the flow simply writes:

Q 0,t (P ) = P -0 t (1)
The subscript 0 indicates that only one channel (i.e. the ground state) is open. Such a formula holds for P > P 0 = 0 but less than P 1 , the critical pressure at which a second channel opens. We will see in detail in the next section how it can be determined. Let us now consider the case of a tree with two open channels (see figure 4).One is the ground state with energy 0 depicted in blue, the other is a channel with energy α depicted in orange. These two channels have a common part of length q0α and we denote with 0α the sum of the thresholds along this common portion. Along it, the fluid flows along a single channel with a pressure difference = P -P , being P the pressure at the bottom of the common part. We can thus write the flow adapting (1), upon replacing the length t → q0α , the energy E 0 → E 0α and the pressure difference P → P -P . Given these ingredients, the flow reads:

Q 1,t (P ) = P -P -0α q0α (2)
The pressure P can still be determined using the conservation of the flow, i.e. Eq. ( 4). Indeed:

P -P -0α q0α = Q (0) 0,t-q0α (P ) + Q (1) 0,t-q0α (P ) (3)
where Q (0) 0,t-q0α (P ) is the flow along the subtree containing the ground state, Q

(1) 0,t-q0α (P ) containing the other channel (see the two branches of figure 4). Since each of these is a single channel of length t -q0α , we can use once again Eq. ( 1). One has

Q (0) 0,t-q01 (P ) = P -( 0 -0α ) t -q0α , Q (1) 0,t-q0α (P ) = P -( α -0α ) t -q0α (4)
From equation (3) we derive first the expression for P (P ) and then Q 1,t (P ):

P (P ) = q0α ( 0 + α ) + (t -q0α )P t + q0α -0α , Q 1,t (P ) = 2 t + q0α P - 0 + α 2 (5)
Finally we impose the continuity of the flow Q 0,t ( P1 ) = Q 1,t ( P1 ) to find the crossover pressure between one and two channels

P1 = 0 + t t -q0α ( α -0 ) (6)
The criterion to select the first excited channel that opens above P 0 , is that the pressure P 1 is the smallest among all the P1 computed for all possible two-channel geometries. This translates into

P 1 = min P1 P1 = 0 + min α =0 t t -q0α ( α -0 ) (7)
The channel satisfying the minimum condition is denoted by α 1 . As discussed in the main text, one sees from this expression that the minimization procedure involves a competition between the energy cost α and the overlap q0α .

In the next section, we show that a similar behavior holds even where three channels are open.

Expression of the flow for three channels

In this section, we report the behavior of the flow when three channels are open. There are three possible configurations for the position of the second excited channel with respect to the ground state and the first one. They each lead to a slightly different expression for the pressure P 2 , but all simplify to P 2 = α2 in the limit t → ∞.

Case I

The first case is the simplest: the second channel opens with a common overlap q0α2 = qα1α2 with the ground state and the first channel. See figure 5 left. By construction, q0α2 < q0α1 . The pressure P 2 reads:

P 2 = min α2 ={α1,0} α2 - q0α2 t + q0α1 -2q 0α2 ( α1 + 0 -2 α2 ) (8)
In the limit t → ∞ we saw that q0α1 = O(1); from this and q02 < q01 , it follows that q02 = O(1) and P 2 = α2 .

algorithm valid for any directed network and discussed in the supplementary material of [21]. In practice, the critical pressure P 0 is the ground state of the associated directed polymer model. The other pressures are found following an iterative procedure. At each step we denote by C s the subtree made of open channels exclusively for P ∈ [P s , P s+1 ] (where P s+1 is not yet determined). Then, 1. Using Kirchhoff's equations, we determine the pressure at each node of C s as a function of the applied pressure P (these functions are linear in P ).

2. We find all the subtrees connected to a node in C s which are still closed. In figure 6 they are highlighted by ovals.

3. For each subtree we find its ground state, denoted by

E (1)
s , E

s , . . . . They are indicated with dashed lines in figure 6.

4. We find the minimal applied pressure P (1) s , P (2) s , . . . at which the ground state of each subtree opens. This is done by setting the node pressure at the root of each subtree (determined in step (1) using Kirchhoff's equation)

equal to E (1)
s , E (2) s , . . . . 5. We pick as a new channel the one with minimal opening pressure, namely P s = min α P (α) s . In numerical simulations we stop at t = 23 due to the exponential growth in t of the number of configurations on the Cayley tree.

Directed polymer on the Cayley tree: the KPP approach

The problem of a directed polymer on the Cayley tree can be studied using the formalism of the discrete KPP equation (see [23,25] for the original literature). A binary Cayley tree of t levels has 2 t-1 distinct paths denoted by α. Each path α has an associated energy given by the sum of the thresholds along the path:

α = i∈α τ i ( 13
)
Since any two paths α, α have some bonds in common, the energies are correlated random variables. As a concrete example, for Gaussian thresholds τ i ∼ N (0, σ 2 ), the correlation of the energies is given by

α α = qσ 2 (14)
where q is the overlap between the two paths, namely the number of common bonds. The overbar indicates the average over the disorder. A directed polymer on the Cayley tree with t + 1 levels can be constructed by considering two independent trees with t levels and joining them together by adding a bond with threshold τ . This hierarchical structure can be exploited when we want to study functions of the energies with a multiplicative form

G t (x) = α g(x -α )) (15)
for some g(x). Indeed G t+1 (x) and G t (x) can be related by the following discrete KPP recursive equation:

G t+1 (x) = dτ p(τ )G 2 t (x -τ ) (16)
where p(τ ) is the probability distribution of a single threshold and the initial conditions reads G t=0 (x) = g(x).

A few relevant examples of functions g(x) are:

• g(x) = θ(-x) (with θ(x) the Heaviside step function) used for the study of the minimum energy of the directed polymer. In particular we have

Ω t (x) = α θ( α -x) = θ( 0 -x) , 0 = min α α ( 17
)
where 1 -Ω t (x) is the cumulative distribution of the minimum energy of a directed polymer of t levels.

• g(x) = λ θ(x) = 1 + θ(x)(λ -1) for the study of the generating function Ψ t (x; λ) of n t (x), namely the number of states with energy smaller than x (see [23]):

Ψ t (x; λ) = λ nt(x) , n t (x) = α θ(x -α ) (18)
In the following subsection we exploit these recursive KPP equations to describe various properties of the directed polymer on the Cayley tree which are relevant for the study of the Darcy law.

Discrete KPP equation

At large t, the discrete KPP equation ( 16) allows for travelling wave solutions of the form

G t (x) = w(x + ct) when t → ∞ ( 19
)
where c is the velocity of the traveling wave. When the initial conditions G 0 (x) = g(x) satisfy lim x→∞ g(x) = 0 and lim x→-∞ g(x) = 1, c is positive and the front moves backward (when the limits are reversed is c negative and the front propagates in the forward direction). The form of the front w(x) satisfies the following fixed point equation:

w(x + c) = dτ p(τ )w 2 (x -τ ) (20)
We consider initial conditions of the form g(x)

x→-∞ = 1 -e βx . By substituting this form into equation ( 20) and expanding at first order in e βx one obtains an equation for c as a function of β:

c(β) = 1 β log 2 dτ p(τ )e -βτ (21)
The velocity is shown, as a function of β, in figure () and presents a minimum at β c . It was proven [29,38,39] that the increasing branch of c(β) is unstable and that for β > β c both the velocity and the front shape freeze, namely c(β > β c ) = c(β c ) and the front shape takes the form w min (x) satisfying the fixed point equation:

w min (x + c(β c )) = dτ p(τ )w 2 min (x -τ ) (22)
This equation corresponds to equation (11) of the main text. In the Gaussian case, τ ∼ N (0, σ 2 ), we find from equation ( 21) c(β) = ln 2/β + βσ 2 /2 with a minimum at βc = 2 ln 2/σ. For β > βc the value of the velocity freezes at c(βc). We report in dashed lines the unstable branch of c(β) for β > βc. Right panel: numerical solution for Ωt(x) (17) for different values of t (black lines). At large t, Ωt(x) approaches the stationary solution wmin(x) (red line), up to a t-dependent translation. Note that wmin(x) moves exactly at velocity c(βc) while the location of the front Ωt(x) is -c(βc)t up to logarithmic corrections (23).

To study the minimum energy of the directed polymer the appropriate initial condition is g(x) = θ(-x) which corresponds to β → ∞. Hence the directed polymer with minimal energy on a Cayley tree with t levels for large t is

P 0 = 0 = -c(β c )t + 3 2β c log t + χ 0 (23)
The extensive part of the minimal energy, -c(β c )t, comes from the velocity selection criterion. The subextensive logarithmic correction cannot be obtain by our first order analysis of equation ( 20) but was proven by Bramson [39].

The term χ 0 is a random variable of order 1 whose cumulative distribution is 1 -w min (x).

Average number of states above the minimum

Using the KPP formalism we can also study the number of states with an energy less or equal than x:

n t (x) = α θ(x -α ) (24)
Indeed by choosing as initial condition g(x) = 1 + θ(x)(λ -1) we can obtain the generating function of n t (x), namely Ψ t (x) = λ nt(x) . However we are interested in the average number of states having an energy α bigger than the minimal energy 0 of an amount x. We denote this quantity by:

m (full) t (x) = α θ(x + 0 -α ) = n t ( 0 + x) ( 25
)
This quantity is more complex than n t (x) since it involves the minimal energy 0 . However, following an approach introduced in [23], it is possible to derive an equation for m

(full) t

(x) that we can solve numerically. We sketch here the derivation. The reader uninterested to the technical details can skip directly to Eqs. (32,33).

The first step is to introduce the generating function of m (full) t (x) with |λ| < 1:

χ t (x; λ) = λ m (full) t (x) (26)
This function satisfies the following identity:

χ t (x; λ) = dx λ nt(x +x) δ(x -0 ) = 1 + ∂ x dx λ nt(x +x) θ( 0 -x ) (27)
To prove the validity of (27) we used (∂

x λ n )θ = (∂ x λ n )θ = ∂ x (λ n θ) -λ n (∂ x θ).
The integrand of equation ( 27) satifies a discrete KPP equation ( 16) since it can be recast into a multiplicative form:

λ nt(x +x) θ( 0 -x ) = α λ θ(x +x-α ) θ( α -x ) (28)
In this case the initial condition reads g(x ) = λ θ(x +x) θ(-x ). The average number of states above the minimum can be written as:

m (full) t (x) = ∂ λ χ t (x; λ)| λ=1 ( 29
)
We thus expand for λ = 1 -ε

λ nt(x +x) θ( 0 -x ) = Ω t (x ) -εR t (x ; x) + O( 2 ) ( 30
)
The function Ω t (x) is defined in (17). This expansion leads to:

m (full) t (x) = dx ∂ z R t (x ; x) = dx r t (x ; x) ( 31
)
where we set r t (x ; x) = ∂ x R t (x ; x). Now, by plugging the expansion (30) in (16) and taking the first order in ε, we obtain a linear equation for R t (x ; x)

R t+1 (x ; x) = 2 dτ p(τ )Ω t (x -τ )R t (x -τ ; x) , R t=0 (x ; x) = 1 2 θ(x + x)θ(-x ) . ( 32
)
Moreover because of the linearity of (32), r t (x; z) satisfies

r t+1 (x ; x) = 2 dτ p(τ )Ω t (x -τ )r t (x -τ ; x) , r t=0 (x ; x) = 1 2 δ(x + x) (33)
By numerically solving simultaneously the equation for Ω t (x) and ( 33), one has direct access to the t dependence of m (full) t (x) (31).

Average number of states above the minimum with maximal overlap

In the main text we introduced the quantity m q,t (x) which counts the average number of states above the minimum with a maximal overlap q. Simply by definition lim q→t m q,t (x) = m (full) t (x). On the contrary, here we are interested in the limit m q (x) = lim t→∞ m q,t (x). States of an infinite tree with maximal overlap q can be built by first taking a tree of q levels where the thresholds on each bond are distributed according to p(τ ) and then adding the q + 1 level with thresholds distributed according as χ 0 ∼ -w min (x). Indeed since we will be interested only in differences in energies we can use the distribution of fluctuations around the extensive part of the minimal energy. Then m for this tree corresponds exactly to m q (x) for an infinite tree. The only difference with respect to equations (32,33) is that term Ω t (x) is replaced w min (x + c(β c )q). The equation (33) becomes:

r q+1 (x ; x) = 2 dτ p(τ )w min (x -τ + c(β c )q)r q (x -τ ; x) (34)
Using r q (x ; x) we obtain:

m q (x) = dx r q (x ; x) . ( 35
)
In the main text, in Figure 3 (middle panel), we show that m q (x) for large q converges to e βcx and is different from

m (full) t (x).

Numerical solutions

The directed numerical integration of Ω t (x) (17) is not feasible at large t as the position of the front moves in the backward direction. In order to maitain the integration limits fixed in a window [x min , x max ] one needs to shift the front by fixing its position as t grows. Moreover to ensure a better numerical stability we evolve H t (x) = e -βcx (1 -Ω t (x)) then the equation for (17) becomes:

H t+1 (x) = dτ p(τ )[2H t (x -τ ) -e βc(x-τ ) H t (x -τ )] (36)
with p(τ ) = p(τ )e -βcτ and initial condition H 0 (x) = θ(x)e -βcx . To fix the position of H t (x) we compute numerically (Riemann sum) b t = dxH t (x)e βcx . At each step t if b t -b 0 > 0 we make the substitution H

t (x) → H t (x -(b t -b 0 )).
As we shift we should assign a value of H t (x) to the points in [x min , x min +b t -b 0 ]; to do so interpolate the last 10 points using a spline. We employed p(τ ) = 1 √ 2π e -x 2 /2 and discretized our integration limits [x min , x max ] in N = 150000 points with x min = -10000 and x max = 30. The equation for r t (x; z) (33) is solved in parallel to H t (x) taking into account the aforementioned translations. We use the same procedure to solve for r q (x; z) (34).

Large q limit of maximal overlap states

In the main text (12) we report the convergence of m q (x) to e βcx as q → ∞ . To get further insights about the large q behavior of m q (x), we observe that Eq. ( 34) can be cast into the application of a q-independent linear operator, by going in the co-moving frame. In other words, we introduce the shifted quantity r (s) q (x ; x) = r q (x -c(β c )q; x). In this way, Eq. ( 34) becomes r (s) q+1 (x ; x) = 2 dτ p(τ )w min (x -τ -c(β c ))r

(s) q (x -τ -c(β c ); x) = [L • r (s) q ](x ) = [L q • r (s) 0 ](x ) (37)
where we implicitly defined the linear operator L and in the last equality we iterated the linear equation to obtain a formal solution starting from the initial condition. The advantage of the formulation in terms of the shifted quantity r (s) q (x ; x) is that it involves a linear operator L independent of q. In the limit q → ∞, the operator L acts as a projector on its largest eigenvector (assuming a gap is present). This can be easily identified taking the derivative of Eq. ( 22) with respect to x , which shows

ρ min (x ) = [L • ρ min (x )](x ) . ( 38
)
with ρ min (x) = -w min (x). In other words, the probability distribution function of the fluctuations around the minimum is an eigenvector. However, since the operator L is not self-adjoint, in order to compute the projector on the eigenvector in (38), we need to determine the corresponding left eigenvector min (x). We introduce the standard scalar product between two functions (x) and r(x) as , r = dx (x)r(x) . The explicit form of min (x) cannot be determined in general as it depends on the specific form of the threshold distribution p(τ ). So, it might look surprising that eventually the large q drastically simplify to a universal form, but a subtle mechanism is at play. Indeed, in the limit x → -∞, we assume min (x) ∼ e -βx . Plugging it in Eq. (40) and using that w min (x → -∞) = 1, we can determine the value of β. We obtain 2 dx e -βx p(x -x -c(β c )) = e -βx (44)

and after the change of variables τ = x -x-c(β c ), it coincides with Eq. ( 21) with β = β c . Since ρ min (x )

x→-∞ -→ x e βcx and r (s) 0 (x ; x) = ρ min (x + x), we see that both the numerator and denominator in Eq. ( 12) are formally infinite as the integrand have a finite limit for x → -∞. So the ratio in Eq. (42) needs to be evaluated by a limiting procedure. In order to regularize we introduce a cutoff Λ and set

5 FIG. 3 .

 53 FIG. 3. Left: Exact numerical solution using the algorithm of section B of [26]. The threshold are Gaussian distributed with zero mean and variance σ 2 = 1/12. Hence βc = √ 2 ln 2/σ. Numerics is performed for moderate sizes t = 12, 15, 17, 19, 21, 23. The analytical prediction (solid line), valid when t → ∞, is approached, but finite-size effects are still strong. Middle: m q (x) is obtained by numerical integration (section C of [26]) of Eqs. (9) starting from an initial threshold drawn from the fixed point distribution -w (x). m (full) ∞ (x) (dashed black line) is instead obtained for an initial threshold drawn from ρ(t) and for t = 10000. Right: Average effective permeability κ for the data of Fig. 3 left. Only the first ∼ t channels are sufficient to reach an effective permeability close to the Newtonian one. Inset of right figure: t dependence of n (ch) SAT , defined as the average number of channels needed to saturate at the effective permeability δ • κ Newton .

Explicit derivation flow up to two channels FIG. 4 .

 channels4 FIG. 4. Schematics for the Cayley tree with two open channels.

FIG. 5 .

 5 FIG. 5. Schematics for the Cayley tree with three open channels in the three possible geometrical arrangements.

7 .

 7 Left panel: the dispersion relation of the velocity (black solid line) as a function of β.

FIG. 8 .

 8 FIG. 8. Construction of the states with maximal overlap q.

( 39 )

 39 Then, min (x) must satisfymin (x) = [L † • min ](x) (40)where the adjoint of L satisfies , L • r = L † • , r . In this way, we can formally express lim the x-dependence is hidden in the scalar product in Eq. (42), which depends on min (x ). Explicitly, the dual operator in Eq. (40) takes the form[L † • ](x) = 2w min (x) dx (x )p(x -x -c(β c ))(43)

  (x) = θ(x + Λ) min (x). Then, we havelim q→∞ m q (x) = lim Λ→∞ ∞ -Λ dx min (x )ρ min (x + x ) ∞ -Λ dx min (x )ρ min (x ) = lim Λ→∞ e βcx 0 -Λ dx (x + x) + O(1) 0 -Λ dx x + O(1)= e βcx (45) which gives the expected exponential behavior.

Acknowledgements. -We thank X. Cao for discussions and suggestions. ADL acknowledges support by the ANR JCJC grant ANR-21-CE47-0003 (TamEnt). We acknowledge the support of the Simons foundation (grant No. 454941, S. Franz). This work was partly supported by the Research Council of Norway through its Center of Excellence funding scheme, project number 262644. Further support, also from the Research Council of Norway, was provided through its INTPART program, project number 309139. This work was also supported by "Investissements d'Avenir du LabEx" PALM (ANR-10-LABX0039-PALM).

Case II

The second case corresponds to the opening of the second excited channel from the ground state with an overlap q0α2 > qα1α2 = q0α1 . See figure 5 middle. The pressure P 2 reads:

When t → ∞, the previous argument for the first excited channel sets q0α1 /t ≈ 0. In this limit, the resulting expression for P 2 is:

This expression is identical to equation ( 7) with the substitution α → α 2 , and applying the same arguments of the first channel we arrive at setting q0α2 /t ≈ 0, leading P 2 = α2 .

Case III

The last case is the mirror of the previous one, with the second channels that opens from the first one with overlap qα1α2 > q0α1 = q0α2 . See figure 5 right. The pressure P 2 reads:

When t → ∞, the previous argument for the first excited channel sets q0α1 /t ≈ 0. In this limit, the resulting expression for P 2 is:

This expression is again similar to equation ( 7) and applying the same arguments of the first channel we arrive at setting qα1α2 /t ≈ 0, leading P 2 = α2 .

Exact numerical solution of the flow for the Cayley tree

The Darcy flow on the Cayley tree can be exactly numerically solved for moderate t and the code used in the paper is available here. The solution consists in finding the pressures P 0 , P 1 , P 2 , . . . at which each new channel opens. The


Christian Kankolongo

Didier Lasseux

Marc Prat

Florent Tony Zaouter

F Ledrappier

Imbibition dynamics including corner

Keywords: imbibition, corner flow, sealing

come

Introduction

In the field of sealing, mechanical assemblies commonly consist of two flanges made of hard metallic material and a seal made of soft metallic material as shown in figure 1a. These flanges have surface defects due to the manufacturing process, e.g., face-turning, which produces a characteristic spiralgrooved profile as shown in figure 1b. The residual opening defects thus constitute leak paths for the fluid to be sealed and the modeling of fluid flow in the resulting connected aperture field between the high and low pressure zones is essential in order to predict seal efficiency. To help understand the flow within the connected channels, in particular under two-phase flow conditions, (i.e., when the fluid to be sealed invades the contact initially saturated by air), it is of interest to analyze the case of a rather simple situation made of single spiral-grooved channel, taking into account the presence of sharp corners as shown in figure 1c. In this context, liquid films can develop in the corners of the channel in the case of imbibition, as studied in [1,2]. Indeed, the presence of the corner films may significantly modify the overall imbibition kinetic and can explain, for instance, capillary instabilities leading to gas trapping. For instance, Dong and Chatzis [1] studied the imbibition of a wetting liquid film in the corner without liquid invasion in the whole channel as shown in figure 1d. By making use of similarity solutions, they showed that the extension of the liquid film and the leak volume in the corner of a rectilinear channel of square cross-section are proportional to the square root of time when flow is triggered by the capillary pressure. Amyot [2] made a similar study in a rectilinear channel of triangular cross-section (figure 1d) with a flow rate imposed at the channel inlet. In accordance with experimental observations, it was shown that the dynamics of the central meniscus is disturbed by that of the liquid film in the corner. The present work reports on analytical solutions for the dynamics of the central meniscus, including that of the liquid film in the corners of a spiral-grooved channel, complementing the analysis reported in the latter reference.

Model

The interest is focused on imbibition of a liquid (the β-phase) in a spiral-grooved channel of triangular cross-section as schematized in figure 1c, initially saturated by a gas (the γ-phase). It is assumed that tan ψ ≪ 1 and the gas phase remains at constant pressure. Under these assumptions and on the basis of the lubrication theory, the Reynolds' equations describing mass and momentum conservation for the β-phase occupying the region A β , are given by

∂q βi ∂x i = 0, q βi = -ρ β h(y) 3 12µ β ∂p β ∂x i , i = 1, 2, (1)
In addition, assuming that the capillary number is small enough (Ca = µ β u/σ < tan ψ, u being a characteristic velocity in the β-phase), the boundary condition at the fluid-fluid interface can be shown to reduce to the Young-Laplace relationship p γ -p β = σκ, at x = L(y, t).

(2)

In these equations, q β is the lineic mass flow-rate, p β and p γ are the liquid and gas pressure, respectively, ρ β is the liquid density, µ β is the liquid dynamic viscosity, h(y) is the aperture, κ is the double mean curvature and σ is the surface tension. For completeness, the kinematic condition at the fluid-fluid interface parametrized by y = e(x, t) must be considered, i.e. ∂e ∂t +

q βx (e) ρ β h(e) ∂e ∂x -q βy (e) ρ β h(e) = 0. By integration of the above equations, considering an imposed volume flow rate, Q 0 , at the channel inlet, the expression of the liquid film extension as a function of time, L f (t), can be obtained, as in [2], by making use of similarity solutions. It is given by

L f (t) = 5 24 1/5 tan(ψ) 1/5 Q 1/5 0 σ cos θ e µ β 2/5 t 3/5 . (3)
A more detailed solution can be derived, that reads

L f (t) = 5 12 1/5 cos 4 (θ e + ψ/2) cos 4 (ψ/2) A 1/5 Q 1/5 0 σ µ β 2/5 t 3/5 . (4)
In equations ( 3) and ( 4), θ e is the equilibrium contact angle and A is a constant which depends of ψ and θ e . The evolution of the liquid film extension, L f (t), is reported in Figure 2a according to the two solutions given in equations (3) (ST) and (4) (SA) using values of the different parameters given in [2]. Note that the two solutions are very close to each other. This is confirmed by the ratio, R ST /SA , of the two expressions given in equations ( 3) and (4), represented in Figure 2, which remains close to unity for an angle of the corner, ψ such that tan ψ < 0.5.

Coupling between the dynamics of the liquid film flow in the corner and that of the central meniscus (L(y, t) in figure 1c) will be further discussed in our presentation.

Références

[1] Dong, M. and Chatzis, I., The Imbibition and flow of a wetting liquid along the corners of a square capillary tube, Journal of Colloid and Interface Science, 172, 278-288 (1995).

[2] Amyot, O., Contribution à l'étude des écoulements diphasiques à travers un contact rugueux, Thèse de doctorat de l'Université de Poitiers (2004).

Figure 1 :

 1 Figure 1: Flange -seal interface.

Figure 2 :

 2 Figure 2: (a) Evolution of the extension for an imposed volume flow rate ; (b) Ratio of liquid film extensions as a function of ψ.

FICPM-Rueil-Malmaison, France -October 17-19, 2023.