# Nonlinear dynamics of air invasion in one-dimensional
compliant fluid networks

Ludovic Jami, François-Xavier Gauci, Céline Cohen, Xavier Noblin, Ludovic Keiser

► To cite this version:

Ludovic Jami, François-Xavier Gauci, Céline Cohen, Xavier Noblin, Ludovic Keiser. Nonlinear dynamics of air invasion in one-dimensional compliant fluid networks. 2025. hal-05291063

HAL Id: hal-05291063

https://hal.science/hal-05291063v1

Preprint submitted on 30 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Nonlinear dynamics of air invasion in one-dimensional compliant fluid networks

Ludovic Jami,¹, * François-Xavier Gauci,¹ Céline Cohen,¹ Xavier Noblin,¹ and Ludovic Keiser¹, †

¹Université Côte d'Azur, CNRS UMR 7010, Institut de Physique de Nice, 17 rue Julien Lauprêtre, 06200 Nice, France

Vascular networks exhibit a remarkable diversity of architectures and transport mechanisms across biological systems. Inspired by embolism propagation in plant xylem — where air invades water-filled conduits under negative pressure — we study air penetration in compliant one-dimensional hydrodynamic networks experiencing mass loss by pervaporation. Using a theoretical framework grounded in biomimetic models, we show that embolism dynamics are shaped by the interplay between network compliance and viscous dissipation. In particular, the competition between two timescales—the pressure diffusion time, $\tau_{diff}$, and the pervaporation time, $\tau_{pv}$—governs the emergence of complex, history-dependent behaviors. When $\tau_{diff} \sim \tau_{pv}$, we uncover previously unreported nonlinear feedback between the internal pressure field and the embolism front, leading to transient depressurization and delayed interface motion. These results offer a minimal framework for understanding embolism dynamics in slow-relaxing vascular systems and provide design principles for soft microfluidic circuits with tunable, nonlinear response.

# I. INTRODUCTION

Vascular networks are essential architectures in living organisms, spanning a wide range of geometries and transport strategies [1–6]. From distributing nutrients and signaling molecules [7] to maintaining homeostasis, their functions rely on designs that often reconcile transport efficiency [8] with robustness [9]. These networks frequently exhibit a coupling between flow and elasticity, termed hydraulic compliance, where the fluid storage capacity depends on local pressure [10].

This rich interplay has inspired the development of artificial fluidic networks, especially in the context of flexible microfluidics [11]. These systems have revealed a diversity of nonlinear dynamics [12, 13], triggered by inertial effects [14], nonreciprocal elements [15, 16], or active boundaries such as cilia [17]. Such mechanisms can lead to complex and sometimes counterintuitive behaviors — e.g., the emergence of the Braess paradox [18] or fluidic memory elements [19] — highlighting the potential of nonlinear fluid-structure interactions. In confined geometries, both elasto-capillary [20–22] and elastohydrodynamic couplings [23–29] have been harnessed to design responsive and/or adaptable transport systems.

Among natural vascular systems, the xylem in plants is particularly intriguing [30–32]. It sustains upward water transport from roots to leaves, driven by negative pressure gradients established via leaf transpiration. Some species tolerate tensions down to –10 MPa, yet drought-induced embolism — the entry and spread of gas into xylem conduits — can severely disrupt water transport, leading to plant decline or death [33–35]. In leaves, embolism propagation has been directly visualized using high-resolution imaging techniques like the “optical vulnerability method” [36, 37], revealing an intermittent pattern, possibly governed by bordered pits acting as capillary barriers. Furthermore, structural failure, such as wall buckling or collapse, under extreme negative pressures has been documented [38–40], pointing to a role of fluid-structure interactions in embolism dynamics. However, despite progress in imaging, the multiscale and rapid nature of embolism propagation makes it difficult to access the full pressure dynamics experimentally.

To bridge this gap, biomimetic approaches have been developed to replicate xylem behavior in controllable artificial settings. Early models used hydrogels or silicone networks to study transpiration [41, 42], cavitation [43–46], and air-seeding processes [47–49]. More recently, researchers have mimicked bordered pits with microfluidic constrictions, imposing Laplace pressure thresholds that temporarily block air invasion until pressure relaxation permits sudden jumps [50]. The elasticity of the channels plays a key role, modulating how pressure evolve dynamically with volume loss by pervaporation (mimicking evapotranspiration).

Recent work has shown that these principles can be extended to channel arrays inspired by *Adiantum* leaves, reproducing intermittent embolism dynamics consistent with real leaves (Figure 1, [51]). Yet, existing studies consider uniform liquid pressure and focus on pervaporation or capillarity, neglecting viscous dissipation. This assumption breaks down in networks combining wide channels (producing large evaporation fluxes) and narrow constrictions (with high hydraulic resistance), where hydrodynamic pressure variations can rival or exceed capillary thresholds.

<footnote>* jami.ludovic@gmail.com</footnote>

<footnote>† ludovic.keiser@univ-cotedazur.fr</footnote>
Embolism spreading in Adiantum leaves

<figure><img src="image_2.png" /><figcaption>FIG. 1. Left: Picture of the embolism propagation in an *Adiantum* leaf (scale bar = 5 mm) and embolism growth as a function of time for different veins. Courtesy of the authors of Ref. [51], where the analysis was carried out based on the experiments from Brodribb *et al.* [36]. For each vein, the dynamics is highly intermittent, with sudden sub-minute propagation events followed by hour-long periods of rest. Right: Embolism dynamics in an one-dimensional pervaporating biomimetic network made with PDMS (seen from the top), where analogous intermittent propagations were evidenced. Courtesy of A. Pellegrin for the picture (whose vertical dimension represents 4 mm), and of the authors of Ref. [51] for the data. The white vessels are embolized (full of air) while the grey ones are still full of water. Note that the curves representing the embolized percentage of the networks exhibit a convex shape. This is a signature of a coherent network - regardless of the number of channels $N_0$ - where pressure variations are rapidly transmitted throughout the structure. The purpose of our current study is to explore configurations where pressure diffusion is slower and where networks are less coherent.</figcaption></figure>

Embolism spreading in 1D 'biomimetic leaves'

<figure><img src="image_3.png" /></figure>

In this work, we develop a model of embolism propagation in compliant, xylem-inspired, one-dimensional fluidic networks with significant hydrodynamic resistance. While the system is deliberately idealized and not intended to quantitatively replicate plant xylem, it contains a minimal set of physical ingredients — compliance, pervaporation, and capillary thresholds — relevant to embolism dynamics in synthetic and biological conduits alike.

By capturing the coupling between pervaporation-induced flows, Laplace thresholds, and viscous pressure drops, we explore how key design parameters govern the dynamics. In particular, we demonstrate that when the pressure diffusion timescale is comparable to or larger than the pervaporation timescale, a nonlinear feedback emerges between the embolism front position and the internal pressure field. This coupling leads to complex, history-dependent dynamics, where embolism progression is no longer solely dictated by local capillary pressure thresholds but also by the evolving pressure distribution across the network. This framework offers interesting perspectives on embolism spreading in plants, where pressure heterogeneity may naturally arise from long transport paths and anatomical complexity. It also informs the design of soft, nonlinear fluidic circuits with emergent behaviors.

The paper is structured as follows. In Sec. II, we describe the physical model and the architecture of the networks. In Sec. III, we present the numerical resolution in both discrete (Subsec. IIIA) and continuous (Subsec. IIIB) settings, and analyze how network topology influences pressure drops (Subsec. IIID). In Sec. IV, we identify asymptotic regimes where pervaporation dynamics are slow (Subsec. IV A) or fast (Subsec. IV B) compared to pressure diffusion. We finally discuss the broader implications for nonlinear transport in engineered fluidic networks (Subsec. VA) and plant physiology (Subsec. VB).

## II. DESCRIPTION OF THE SYSTEM: PERVAPORATION, ELASTO-CAPILLARY AND HYDRODYNAMIC COUPLING OF THE CHANNELS.

We investigate series of $N_0$ microchannels fabricated in PDMS, each of length $l$, width $w$ and height $h$, separated by constrictions of length $l_c$, width $w_c$ and height $h_c$ (Figure 2). In the numerical simulations, we restrict ourselves to constrictions having the same height as the channels ($h = h_c$), but in analytical expressions we retain distinct notations to clearly separate their geometric contributions. The total thickness of the PDMS leaf $H$ is uniform, such that the thickness separating the ceiling of the microchannels and constrictions from the external air is $\delta = H - h$. The configuration is a dead-end series: the last ($N_0^{th}$) channel is closed by a wall, and the first constriction being
<figure><img src="image_4.png" /><figcaption>FIG. 2. Illustration of a series of compliant pervaporating channels connected by constrictions. (a) Illustration of the embolism propagation in the channel series. The width and length of the channels w and l and of the constriction are indicated. In step i, the embolism front is halted at the constriction between channels i − 1 and i . To cross the constriction (step ii), pressure p$_{k}$$_{+1}$ must drop below the Laplace pressure imposed by the constriction (eq. 4). In step iii, the embolism has fully invaded channel i , and the interface reaches the next constriction. Step iv is equivalent to step i but with the embolism front, indexed by k , now located between channels i and i + 1. (b) Scheme of the deformation of the channel when losing water volume by evaporation. In this scheme, the channel dimensions h , w , l and δ are defined, as well as the section of the deformed channel S = S$_{0}$ + Δ S with S$_{0}$ = hw the section at rest. (c) Pressure difference between two successive channels i and i + 1 generates a flow through the connecting constriction.</figcaption></figure>

connected to the ambient air by an opening. The thickness δ is thin enough to enable a relatively fast pervaporation flux. For water-permeable materials like PDMS, the evaporation flux per unit length is given, following Dollet et al. , by [48]:

$$
j = j_0 \left\{ \frac{w}{\delta} + \frac{2}{\pi} \left[ \ln \left( \frac{(H + \delta)h}{\delta^2} \right) + \frac{H}{\delta} \ln \left( \frac{H + \delta}{h} \right) \right] \right\} \quad (1)
$$

with j$_{0}$ = M ρ D$_{p}$c sat p (1 − RH ), where h is the channel height, H = δ + h , M and ρ are, respectively, the molar mass of water and the density of liquid water, D$_{p}$ and c sat p are the diffusivity and saturation concentration of water vapor in the pervaporating material and RH is the relative humidity of the air surrounding the channel network. The second term in the square brackets of equation 1 is a geometrical logarithmic correction term to account for cases when the height of the channel is comparable to the channel width and similar or greater than the thickness of its upper wall (that is when h ∼ w and h ≳ δ ). Pervaporation causes a gradual loss of liquid volume, leading to several possible consequences.

- In the absence of a gas phase, volume loss translates into a pressure drop within the channel. For soft materials like PDMS, this causes an elastic deformation of the upper wall. To describe the relationship between the water volume lost by the channel Δ V$_{i}$ and its pressure decrease p$_{i}$ , we define the compliance C as: Δ V$_{i}$ = Cp$_{i}$ . We here use the notation p$_{i}$ for the difference between the pressure in the channel i with the atmospheric pressure in the air surrounding the PDMS leaf. In the case of rectangular channels with thin upper wall ( δ ≪ w ), thin plate elastic theory enables to express the compliance as a function of the geometric dimensions of the channel and the mechanical parameters of the elastic material [50, 51]:

$$
C = \frac{lw^5 (1 - \nu^2)}{\delta^3 60E} \qquad (2)
$$

with E the Young modulus of the channel material and ν its Poisson’s ratio. In practice, we take ν = 0 . 5, a standard value for nearly incompressible PDMS elastomers.
<table><caption>TABLE I. Physical and geometrical parameters used in the model, grouped by category.</caption><thead><tr><th>Category</th><th>Parameter</th><th>Symbol</th><th>Typical value (unit)</th><th>Fixed / Explored</th></tr></thead><tbody><tr><td rowspan="5">Geometric (channels)</td><td>Number of channels</td><td>N<sub>0</sub></td><td>10 or 64</td><td>Explored</td></tr><tr><td>Channel width</td><td>w</td><td>50–600 µm</td><td>Explored</td></tr><tr><td>Channel length</td><td>l</td><td>0.5 or 5 cm</td><td>Explored</td></tr><tr><td>Channel height</td><td>h</td><td>50 µm</td><td>Fixed</td></tr><tr><td>PDMS ceiling thickness</td><td>δ</td><td>H − h</td><td>Fixed</td></tr><tr><td rowspan="3">Geometric (constrictions)</td><td>Constriction width</td><td>w<sub>c</sub></td><td>1–100 µm</td><td>Explored</td></tr><tr><td>Constriction length</td><td>l<sub>c</sub></td><td>0.3–1000 mm</td><td>Explored</td></tr><tr><td>Constriction height</td><td>h<sub>c</sub></td><td>50 µm</td><td>Fixed</td></tr><tr><td rowspan="2">Material properties</td><td>Young's modulus of PDMS</td><td>E</td><td>1 MPa</td><td>Fixed</td></tr><tr><td>Poisson's ratio of PDMS</td><td>ν</td><td>0.5</td><td>Fixed</td></tr><tr><td rowspan="4">Fluid properties</td><td>Dynamic viscosity of water</td><td>η</td><td>1 × 10<sup>-3</sup> Pa·s</td><td>Fixed</td></tr><tr><td>Surface tension of water</td><td>γ</td><td>72 × 10<sup>-3</sup> N/m</td><td>Fixed</td></tr><tr><td>Density of water</td><td>ρ<sub>w</sub></td><td>1000 kg/m<sup>3</sup></td><td>Fixed</td></tr><tr><td>Molar mass of water</td><td>M</td><td>18 g/mol</td><td>Fixed</td></tr><tr><td rowspan="3">Pervaporation parameters</td><td>Diffusivity of water in PDMS</td><td>D<sub>p</sub></td><td>1 × 10<sup>-9</sup> m<sup>2</sup>/s</td><td>Fixed</td></tr><tr><td>Saturation conc. in PDMS</td><td>c<sub>p</sub><sup>sat</sup></td><td>40 mol/m<sup>3</sup></td><td>Fixed</td></tr><tr><td>Relative humidity</td><td>RH</td><td>0 %</td><td>Fixed</td></tr></tbody></table>

- When a gas bubble is present, pressure differences between adjacent channels drive water flow through the constrictions. The flow rate between channel *i* and *i* + 1 is governed by the hydrodynamic resistance *R* of the constriction $q_i = \frac{p_{i+1} - p_i}{R}$. Assuming $w_c < h_c = h$, the hydrodynamic resistance is expressed as [52]:

$$R = \frac{12\eta l_c}{h_c w_c^3} \frac{1}{1 - 0.63 \frac{w_c}{h_c}}, \quad (3)$$

with $\eta$ the dynamic viscosity of the water.

When embolism propagates across the channel series, the interface can only cross a constriction if the pressure in the downstream channel drops below a Laplace pressure threshold $p_c$. For rectangular constrictions, a completely wetting liquid and assuming $w_c < h_c$, the threshold Laplace pressure is given by [53]:

$$p_c = \gamma \kappa_c \quad (4)$$

with the curvature:

$$\kappa_c = \frac{1}{w_c} + \frac{1}{h_c} + \sqrt{\left(\frac{1}{w_c} - \frac{1}{h_c}\right)^2 + \pi \left(\frac{1}{h_c w_c}\right)}, \quad (5)$$

and $\gamma$ the surface tension of water.

These equations 1, 2, 3 and 4 provide a consistent description of embolism dynamics under a set of reasonable assumptions. Pressure is assumed homogeneous within each channel owing to the high contrast in viscous dissipation between channels and constrictions (Eq. 3). Contact line pinning and additional contact line dissipation are neglected. The flow regime is considered Stokesian, with negligible inertia. Finally, the channel deformation is treated within linear thin plate theory (Eq. 2).

The full range of parameter values used in this study is summarized in Table I. In the following sections, embolism dynamics will be described in general terms, without making further assumptions on the functional forms of *R*, *C*, *j* and $p_c$.
<figure><img src="image_5.png" /><figcaption>FIG. 3. Hydrodynamic resistance initially delays embolism front propagation, followed by a rapid "catch-up" dynamics. Position of the front (channel index) versus time is plotted as hydrodynamic resistance R is increased. For vanishing hydrodynamic resistance (dark magenta), the dynamics is that of a truncated exponential, as demonstrated by Dollet et al. [47]. As R increases (lighter magenta), an inflection point appears in the shape of L vs. t, a signature of the delay induced by the slow diffusion of the pressure in the resistive network. The values of the pressure diffusion timescale $\tau_{\text{diff}} = N_0^2 RC$ are expressed in the legend (see also Eq. 17), and show that the inflection of the shapes appears as $\tau_{\text{diff}}$ becomes comparable to the duration of the network pervaporation. (a) Example with 1D networks of $N_0 = 10$ channels and with $w = 400 \mu\text{m}$, $l = 5$ cm, $w_c = 15 \mu\text{m}$ and $l_c = 0.5$, 15 or 50 mm. The embolism front advances intermittently, following the sequence of steps illustrated in Fig. 2a. (b) Example given with 1D networks of $N_0 = 64$ channels (as in the following figures), with $w = 400 \mu\text{m}$, $l = 0.5$ cm, $w_c = 10$, 15 or 30 $\mu\text{m}$ and $l_c = 10$ mm. See table I for the other parameters considered in this study.</figcaption></figure>

### III. EMBOLISM PROPAGATION AND PRESSURE DYNAMICS IN ONE-DIMENSIONAL NETWORKS

#### A. Embolism dynamics in a discrete network of channels and constrictions in series.

In fluidic networks as represented in Figure 2, the pressure in a channel is coupled to the one in the other channels by the constrictions that connect them and enable the fluid to flow between channels. The dynamics of the pressure ($p_i$) in the channel i away from the embolism front is governed by the pressure in the neighboring channels and its pervaporation flux. One can write the following local water volume conservation equation for the channel i:

$$ \frac{d}{dt}p_i = \frac{(p_{i+1} - p_i) - (p_i - p_{i-1})}{RC} - \frac{jl}{C}, \qquad (6) $$

that has structure of a diffusion equation with a discrete diffusion coefficient $1/RC$ associated to a diffusion time RC. The first term on the right-hand side accounts for pressure diffusion between neighboring channels, while the second represents the pressure drop induced by pervaporation. When the embolism front is stopped at a constriction, the pressure in the channel (k) next to the constriction at which the interface is stopped evolves as:

$$ \frac{d}{dt}p_k = \frac{p_{k+1} - p_k}{RC} - \frac{jl}{C}. \qquad (7) $$

When the pressure $p_k$ in the channel k reaches the Laplace pressure of the constriction, $p_c$, the interface front moves forward in the channel. We consider here that the mechanical relaxation of the channel section is instantaneous after the constriction jump leading to a reduction of the filled length of the channel to $\frac{V_k}{wh} = l + \frac{Cp_c}{wh}$. This assumption is consistent with the assumption made that the flow in constrictions is the dominant dissipative process. Then, the water volume $V_k$ in the draining channel decreases as:

$$ \frac{d}{dt}V_k = \frac{p_{k+1} - p_k}{R} - j\frac{V_k}{wh}. \qquad (8) $$

We solved numerically equations 6, 7 and 8 to capture the intermittent dynamics of the embolism propagating through the series of channels and constrictions. We first illustrate the effect of an increase of the hydrodynamic
<figure><img src="image_6.png" /><figcaption>FIG. 4. Pressure dynamics in the channel series. (a) Pressure evolution in successive channels (color-coded) as a function of time, for $N_0 = 10$. (b) Depression in the successive channels, coded by color, as a function of time, for $N_0 = 64$. Only a subset of 12 channels are shown here. (c) Zoom on the pressure minimum ($p_{min}$ at time $t^*$) observed in (b), where the series has significant hydrodynamic resistance. (d) Pressure dynamics near a single embolism event (steps i–iv in panel a of Fig. 2).</figcaption></figure>

resistance on advancing front dynamics in Figure 3. We can see the initial propagation of the embolism is greatly delayed when hydrodynamical resistance increases. This delay results from increased viscous dissipation in the network, which weakens the pressure gradient pulling the interface forward.

Figure 4, shows the typical evolution of the pressure with time in the different channels of a relatively short ($N_0 = 10$, panel a)) and relatively long ($N_0 = 64$, panel b)) series. In this figure, we illustrate the successive steps of the pressure dynamics in the networks as the embolism front stops at a constriction and then progresses in the successive channel. Pressure dynamics in the channels exhibit pressure dips associated with the passage of the embolism front through the constrictions. On top of this impulse train dynamics, a transient pressure minimum is observed associated with significant viscous dissipation occurring in the system.

## B. A continuum model for the embolism spreading

We now derive a continuous description of embolism propagation from the discrete model, describing the propagation of the liquid length $L(t)$ in the $x$ direction, as well as the pressure field $p(x, t)$. The continuum description is valid when considering long series and looking at their large size and long time dynamics, that is when :

$$L_0 = N_0 l \gg l \quad \& \quad \left( \begin{array}{l} dL, dx \gg l \\ dt \gg \Delta t_{tot} \end{array} \right) \qquad (9)$$
where $L_0$, the cumulated length of all the channels, is a finite quantity and $\Delta t_{tot}$ corresponds to the time that the embolism front takes to pass through a constriction and its following channel. Following this procedure, Eq. 6 becomes:

$$ \partial_t p = D \partial_x^2 p - \frac{j}{c} \qquad (10) $$

where $D = \frac{1}{\tau_c}$ with $\tau = R/l$ and $c = C/l$. This equation is a diffusion equation with a sink term ($-j/c$). Concerning the boundary conditions, one have: in $x = L_0$, at the end of the channel, $\partial_x p = 0$, since no flow occurs at the closed end of the network; and, in $x = L_0 - L$, at the interface, $p = p_c = \gamma\kappa_c$ [54]. Eq. 10 is valid in the liquid-filled domain, whose size decreases as much as the liquid pervaporates and the embolism propagates. Thereby, it is more convenient to rewrite the equation in a fixed domain by doing the following appropriate change of spatial variable $x' = 1 - \frac{L_0 - x}{L}$. We then find:

$$ \partial_t p + \frac{1-x'}{L} \dot{L} \partial_{x'} p = \frac{D}{L^2} \partial_{x'}^2 p - \frac{j}{c} \qquad (11) $$

which holds between $x' = 0$ (corresponding to $x = L_0 - L$, and where the boundary condition $p = p_c$ is imposed) and $x' = 1$ (corresponding to $x = L_0$ and where the no-flux boundary condition $\partial_{x'} p$ is imposed). This variable change introduces an advection-like term: $(1 - x')/L \cdot \dot{L} \partial_{x'} p$.

The system of equations can be closed by expressing volume conservation at the interface:

$$ r \dot{L} S_c = \partial_x p (L_0 - L, t) \qquad (12) $$

with $S_c = S_0 + C p_c / l$ the section of the channel when the pressure equals the Laplace pressure of the constriction (the condition to jump the constriction). From integration of 10, one obtains the following equation for $L$:

$$ S_c \dot{L} + (j + \overline{\partial_t S}) L = 0 \qquad (13) $$

with $\overline{\partial_t S} = c \int_0^1 \partial_t p(x', t) dx'$. Note that Eq. 13 is equivalent to the global volume conservation condition $\dot{V} = -jL$ (see Appendix A).

### C. Comparison with non-deforming systems

We now want to explicit the link between our continuous model given by equations 11 & 13 with the simpler model of drying channels having negligible compliance hydrodynamic resistance. Let first note that, even if $\partial_t S$ does depend on $L$ (via $p$), one can integrate this equation for $L$ as if $\overline{\partial_t S}$ was a prescribed function of $t$. Thereby, one gets:

$$ \begin{aligned} L &= L_0 \exp \left( - \int_0^t \frac{j + \overline{\partial_t S}}{S_c} dt' \right) \\ &= L_0 \exp \left( - \frac{t}{\tau_{pv}} \right) \exp \left( \frac{1}{S_c} \int_0^1 c[p_c - p(x', t)] dx' \right), \end{aligned} \qquad (14) $$

where we introduced $\tau_{pv}$, the characteristic pervaporation time of the channel series:

$$ \tau_{pv} = \frac{S_c}{j} \qquad (15) $$

This equation compares directly with the case of infinitely rigid channels with negligible viscous dissipation reported by Dollet et al. [47]:

$$ L(t) \approx L_0 e^{-t/\tau_0}, \text{ with } \tau_0 = S_0/j, \qquad (16) $$

in the limit of negligible meniscus evaporation. In the latter case, the interface progression is straightforwardly forced by pervaporation, in the absence of compliance of the network. For compliant systems, looking at the right hand side of equation 14, one can see that, in the case where all the volume lost by pervaporation is translated into channel depressurization, that is when $c(p_c - p(x', t)) = jt$, no movement of the front is observed and one has $L = L_0$. Thus, in our configuration, the channel section modification delays the effect of pervaporation on the retraction of the interface, as the water loss is compensated by depressurization and section reduction. Conversely, when the section relaxes, it pulls on the interface. This effect renders interface dynamics more intricate than in the rigid case [47, 48, 55].
<figure><img src="image_7.png" /><figcaption>FIG. 5. **Parametric analysis of the discrete model.** Resolution of the equations 6, 7 and 8 varying the design parameters of the series. First panel line represent the temporal evolution of the embolism front. Second panel line plot the dynamics of the pressure in the channels. Last panel line represent the pressure profile $P_n$ in the series where $P_n$ is the pressure at the channel $n$ when the embolism front just finished the draining of a channel. The series considered here have $N_0 = 64$, $l_c = 1$ cm, $l = 5$ mm and $h = 50$ µm.</figcaption></figure>

D. Influence of the network design on the embolism dynamics

We conducted a parametric analysis of the discrete model described by Eqs. 6, 7 and 8, varying the width of the channels $w$ and the width of the constrictions $w_c$. The main outcomes of this parametric analysis is summarized in the figure 5, where three archetypal cases are presented. For each column, the successive rows present respectively, the embolism front trajectory (first row), the channel depressurization dynamics (second row) and the evolution of the pressure $P_n$ (third row). The first column of Fig.5.A) shows the results for the case of a large constriction width ($w_c = 30$ µm) and small channel width ($w = 100$ µm). This design implies a low hydrodynamic resistance $R$ coupled with a low compliance $C$ of the channel. This case is analogous to the one described in a recent study [51], where pressure tends to homogenize rapidly after each jump of the interface that passes a constriction. The second column of Fig.5.B) combines small compliance ($w = 100$: µm) with large hydrodynamic resistance ($w_c = 10$ µm). A large pressure gradient appears quasi-instantaneously, with a pressure minimum in the terminal channel that largely exceeds the capillary pressure of the constrictions. However, the trajectory of the embolism front in the series is not affected, such that cases A) and B) exhibit the same exponential dynamics. Moreover, cases A) and B) appear to exhibit self-similar pressure profiles (third row on Fig. 5). Finally, the third column of Fig.5.C) corresponds to a network with a large compliance ($w = 400$ µm) coupled to a large hydrodynamic resistance ($w_c = 10$ µm). This design exhibits a large diffusion timescale $RC$ for the pressure. In this case, the trajectory of the embolism front is significantly modified as compared to earlier cases (first row). A large depressurization peak appears (second row) in the channel but only after a significant delay corresponding to about one third of the total pervaporation time of the network. The pressure profile appears no more self-similar (third row), a consequence of the strong coupling between the pressure field and the embolism front dynamics. Initially, the pressure profile is completely flat over a significant length away from the embolism front. As the embolism front advances the pressure gradient expands further over the network. It is only when the pressure gradient expanded towards the terminal channel that the latter starts to re-pressurize. The
pressure typically diffuses on the length $N_0$ of the series with a characteristic time $\tau_{\text{diff}}$ which expresses as:

$$ \tau_{\text{diff}} = N_0^2 RC = \frac{L_0^2}{D}. \qquad (17) $$

Overall, the parametric analysis showed the strong influence of the hydrodynamic resistance $R$ and the compliance $C$ on the onset of nonlinear dynamics and the appearance of a strong pressure minimum.

As we considered here long channel series ($N_0 = 64$), the main features observed from the discrete model (intricate interface dynamics, pressure minimum and pressure gradient profiles) can be reproduced from the continuous model of equations 11 and 13. We conducted a parametric sweep of the continuous model varying several geometric design parameters, whose output is represented in Fig. 6 as a phase diagram. The nondimensionalized phase parameters used in the diagram are $\tau_{\text{diff}}/\tau_{\text{pv}}$, which increases when pressure diffusion delay increases, and $\Delta S_c/S_0$, which characterizes the deformation of the channels owing to the capillary pressure of the constrictions [56]. Note that the continuous model was solved using a simplified version of the pervaporation rate, only keeping the linear dependence of $j$ with the channel width and neglecting the corrective terms of Eq. 1. The maximal deformation of the terminal channel, $\Delta S_{\text{min}}/S_0$, is observed to increase with $\tau_{\text{diff}}/\tau_{\text{pv}}$ (color shifting from blue to red). The ratio $p_{\text{min}}/p_c$ of the pressure minimum and the capillary pressure of the constriction is observed to decrease (marker size reducing) as $\Delta S_c/S_0$ increases. It is also noteworthy that for a ratio $\tau_{\text{diff}}/\tau_{\text{pv}}$ larger than one, the depressurization of the terminal channel is so intense that a total collapse of the terminal channel is observed, with outcomes falling outside the validity range of the current model. Channel collapses are marked with on the diagram.

In the following section, we explore how to rationalize analytically the results obtained from the parametric analysis. In particular, we look for prediction of the amplitude and time of appearance of the pressure minimum (respectively, $p_{\text{min}}$ and $t^*$), as a function of the design parameters of the fluid network.

## IV. ANALYTICAL DEVELOPMENTS

In the previous part, we demonstrated that the collective behavior of pervaporating channels can give rise to complex, nonlinear dynamics. These nonlinear patterns are not only observed at the interface, through constriction-induced jumps [51], but also in the global pressure distribution, which becomes inhomogeneous when the hydrodynamic resistance $R$ and compliance $C$ are large. Two main regimes can be distinguished depending on the ratio of the pressure diffusion timescale $\tau_{\text{diff}}$ to the pervaporation timescale $\tau_{\text{pv}}$. In the following, we provide analytical solutions for the pressure dynamics for the first case ($\tau_{\text{diff}} \ll \tau_{\text{pv}}$). For the second case ($\tau_{\text{diff}} \sim \tau_{\text{pv}}$), we also provide analytical estimates. We do not consider here the regime $\tau_{\text{diff}} \gg \tau_{\text{pv}}$, as it leads to uncontrolled channel collapse, where pressure relaxation is too slow to compensate for pervaporation-driven volume loss. Overall, the aim of this section is to rationalize the dynamics to expect for any chosen design to provide predictive guidelines for the design and analysis of both artificial and biological vascular networks.

### A. Networks with rapid pressure diffusion: $\tau_{\text{diff}} \ll \tau_{\text{pv}}$

#### 1. Discrete Model

For fully coherent networks where pressure relaxation diffuses instantaneously on the entire network, the discrete model of equations 6, 7 and 8 can be solved analytically. Indeed, one can consider that, after a constriction jump and a quick transient time, the water flux and the pressure field reach a (temporary) stationary state during the drainage of the channel. This stationary state imposes that pressure suction at each constriction equals the pervaporative loss of the upward channels. Thus, during the stationary state of the drainage, one has:

$$ \frac{p_i - p_{i+1}}{R} = (N_0 - i)jl \qquad (18) $$

In particular, this equation shows that pressure difference between successive channels should decrease with $i$ as the cumulated pervaporation flux of the upward channels decreases. Notably, equation 18 gives that the pressure gradient $\frac{p_i - p_{i+1}}{l} = (N_0 - i)jR$, at a given channel position $i$ is independent of time for a given network. This explains the observed self-similar pressure profiles in Figure 5. The speed of the interface, $v = \frac{l}{\Delta t_{\text{drain}}}$, $= \frac{(N_0 - k)jl}{hw} = \frac{p_k - p_{k+1}}{Rhw}$, noting $k$ the index of the draining channel, is itself not stationary; it decreases, as well as the pressure gradient near the interface $\frac{p_k - p_{k+1}}{l}$, when the interface moves from a channel to the next (as $k$ increases).
<figure><img src="image_8.png" /><figcaption>FIG. 6. Parametric analysis of the continuum model of equations 11 and 13. The different design parameters of the series considered are plotted in the non-dimensional phase space ($\tau_{diff}/\tau_{pv}$, $\Delta S_c/S_0$). Color and size of the scatter dots represent, respectively, the maximal deformation of the terminal channel ($\Delta S_{min}/S_0$) and their minimal pressure ($p_{min}/p_c$, log scale) over the course of the embolism propagation. Crosses indicate cases where the terminal channel undergoes full collapse owing to insufficient pressure compensation. The data correspond to a wide variety of design parameters with $w \in [50, 600]$ µm, $l_c \in [0.3; 10^3]$ mm and $w_c \in [1; 100]$ µm.</figcaption></figure>

Equation 18 can be integrated to yield $p_i$ explicitly:

$$p_i = -\gamma\kappa - jlR \left( N_0 - \frac{i+k-1}{2} \right) (i-k) \quad (19)$$

Here $-\gamma\kappa$ is the “remanent” Laplace pressure of the channel $k$ that is draining (see eq. 4 considering the dimensions $h$ and $w$ of the channels). Equation 19 predicts that the self-similar pressure gradient profile of figure 5 is quadratic in $i$.

Finally, one can notice that considering series without pressure relaxation delay ($\tau_{diff} \ll \tau_{pv}$) is equivalent to considering small number of channels $N_0$, so that during draining all channels are influenced by the pressure relaxation provided by the advancement of the interface. $N_0$ has thus to be small compared to the number of channels over which the pressure diffuses during drainage $\sqrt{\frac{1}{RC}\Delta t_{drain}} = \sqrt{\frac{1}{N_0}\frac{hw/j}{RC}}$ where $hw/j$ is the characteristic time of pervaporation of an individual channel and RC is the characteristic time of diffusion of the pressure between two channels. This is typically fulfilled when $N_0^2 \ll N^{*2} = \frac{hw/j}{RC}$, that is, when $\tau_{diff} \ll \tau_{pv}$. We can thus express, as a function of the geometric design parameters of the series, the characteristic number of channels $N^*$ for which the pressure relaxation delay begins to occur. Assuming $w_c < h_c = h$ (one has to inverse the order of $w_c$ & $h_c$ in the following expression otherwise) and considering only the first order dependence over the geometric parameters, one gets:

$$N^* \propto \sqrt{\frac{h\delta^4 w_c^3 h_c}{l_c l w^5}} \quad (20)$$

$N^*$ can be considered as the critical network size, above which pressure relaxation cannot keep pace with pervaporation, and nonlinear effects emerge. This equation indicates how the design parameters of the network ($N_0$, $l$, $w$, $h$, $l_c$, $w_c$, $h_c$ and $\delta$) should be selected in order to observe similar embolism dynamics.

One can similarly express the maximal depression in the terminal channel in leading order of the design parameters with: $jlR^{\frac{N_0(N_0+1)}{2}} \propto \frac{wll_c}{\delta w_c^2 h_c} \frac{N_0(N_0+1)}{2}$. Therefore, assuming $w \gg w_c$, the ratio between the minimal pressure in the terminal channel $p_{min}$ and the Laplace pressure imposed by the constriction goes as:

$$\frac{p_{\min}}{p_c} \propto \frac{wll_c}{\delta w_c^2 h_c} \frac{N_0(N_0+1)}{2} \quad (21)$$
Assuming a stationary pressure profile during drainage of the channels leads to an analytical solution for the pressure profile in the series. This enables to rationalize the depressurization peak observed at the terminal channel of series without pressure diffusion delay in the middle column of figure 5 and to relate it to the series design parameters. In this simplified case, the pressure gradient observed is merely the suction profile needed to sustain the pervaporation flux, given the viscous dissipation in the series.

## 2. Continuum limit

In the continuous description, the absence of delay in the diffusion of the pressure relaxation implies that pressure gradient always drains as much water volume as it is lost by pervaporation, i.e. that:

$$ -\frac{1}{\tau} \partial_x p(x, t) = j(L_0 - x), \qquad (22) $$

which is the continuous equivalent of equation 18. Combined with the volume conservation condition 12, it integrates into, $L = L_0 e^{-t/\tau_{pv}}$ for the liquid length (i.e. finding back the exponential decay from Dollet et al. [47]), and for the pressure field:

$$ \begin{aligned} p(x,t) &= p_c - \tau j(x - (L_0 - L)) \frac{L_0 + L - x}{2} \\ &= p_c - \tau j(x - L_0(1 - e^{-t/\tau_{pv}})) \frac{L_0(1 + e^{-t/\tau_{pv}}) - x}{2} \end{aligned} \qquad (23) $$

The absence of pressure diffusion delay implies a clear timescale separation between the evolution of $L$ (having the characteristic timescale $\tau_{pv}$) and $p$ (transiently adapting to any change in $L$ with characteristic timescale $\tau_{diff}$). Here, $p$ is a fast variable ("slaved variable"), instantaneously adapting to the slow evolution of $L$. In fact, equation 22 can be obtained from an integration of the governing equation 11, supposing $p$ stationary and $L$ constant so that the evolution of the pressure field is quasi-static (quasi-stationary) and overall evolves with the characteristic time $\tau_{pv}$.

## B. Networks with large pressure diffusion delay $\tau_{diff} \sim \tau_{pv}$

We now turn to the opposite regime, where pressure diffusion becomes comparable to pervaporation, leading to a delayed depressurization and a stronger coupling between pressure and embolism front dynamics.

### 1. pressure minimum

For networks presenting large pressure relaxation delays, the progression of the embolism front in the channel series owing to pervaporation proceeds over similar timescales as pressure diffusion in the network. Therefore, the embolism dynamics proceeds from a tight coupling between these two effects; such that both diffusion and interface advection terms play a role in the dynamics of the pressure field. These processes have respectively a typical timescale $\tau_{diff}$ and $\tau_{pv}$, and pressure in the network typically evolves with a timescale given by the composition of both processes:

$$ \tau^* = \frac{\tau_{pv} \tau_{diff}}{\tau_{pv} + \tau_{diff}}. \qquad (24) $$

In the early times of the embolism spreading, the gradient profile spreads slowly over the network. The pressure and the section of terminal channels only decrease owing to the volume they lose by pervaporation, regardless of the movement of the interface that does not influence them. A minimum for the pressure $p_{min}$ is then reached at $t^* \sim \tau^*$, when a significant number of channels has been embolized and the relaxation of the pressure has started to diffuse from the front to end of the channel. As shown in figure 7.a, $t^*$ is indeed found to be of the order of $\tau^*$, with $t^*$ varying between $\tau^*$ and $2\tau^*$ for $\tau_{diff}/\tau_{pv}$ varying from $10^{-2}$ to 8. Note that, close to $\tau_{diff}/\tau_{pv} = 1$, we have $t^* = \tau^*$. Moreover, data for $t^*/\tau^*$ all collapse on a unique curve which varies symmetrically around $\tau_{diff}/\tau_{pv} = 1$.

All the values of the minimal pressure $p_{min}$ for the different design parameters are also found to collapse onto a single curve when rescaling the depression $p_{min} - p_c$ by its value in the quasi-stationary regime $j\tau L_0^2$ and looking at the dependence of this depression with the ratio $\tau_{diff}/\tau_{pv}$, as shown in figure 7.b. This collapse suggests that the nonlinear interplay between embolism front dynamics and pressure field evolution is fully captured by the dimensionless ratio
<figure><img src="image_9.png" /><figcaption>FIG. 7. Collapse of the pressure minimum values is obtained by quasi-static and coupled rate process normalisation. a) Time at which the pressure minimum is reached in the channel series $t^*$ normalized by $\tau^* = \frac{\tau_{pv}\tau_{diff}}{\tau_{pv}+\tau_{diff}}$ as a function of the ratio of the diffusion and the pervaporation time. The parametric exploration of the continuous model confirms that $t^*$ is close to $\tau^*$ over a large range of $\frac{\tau_{diff}}{\tau_{pv}}$ values. The red dotted-line is provided as a guide for the eyes, its equation is $\frac{t^*}{\tau^*} = 0.87 \left( \frac{0.49+(\tau_{diff}/\tau_{pv})^2}{\tau_{diff}/\tau_{pv}} \right)^{0.22}$. b) Nondimensionalized pressure minimum of the channel series versus $\tau_{pv}/\tau_{diff}$. The dashed red fitting curve is given by the equation $f(\tau_{diff}, \tau_{pv}) = 1/(1+0.7(\tau_{diff}/\tau_{pv})^{0.63})$.</figcaption></figure>

$\tau_{diff}/\tau_{pv}$. Following this rescaling, noting $\tilde{p} = \frac{p-p_c}{j\tau L_0^2}$ and $\tilde{t} = t/\tau_{pv}$, and further rescaling the liquid length as $\tilde{L} = L/L_0$, we can rewrite the governing equation 11 and 12 as:

$$ \partial_t \tilde{p} + \frac{1-\tilde{x}}{\tilde{L}} d_t \tilde{L} \partial_{\tilde{x}} \tilde{p} = \frac{\tau_{pv}}{\tau_{diff}} \frac{\partial^2 \tilde{p}}{\tilde{L}^2} - \frac{\tau_{pv}}{\tau_{diff}} \quad (25) $$

and:

$$ d_t \tilde{L} = \frac{\partial_{\tilde{x}} \tilde{p}(0, \tilde{t})}{\tilde{L}} \quad (26) $$

with initial conditions $\tilde{L} = 1$ and $\tilde{p} = 0$ and boundary conditions $\tilde{p}(0, \tilde{t}) = 0$ and $\partial_{\tilde{x}} p(1, \tilde{t}) = 0$. Equations 25 and 26 show that any parametric exploration of the embolism dynamics can be reduced to variation of the unique parameter $\frac{\tau_{diff}}{\tau_{pv}}$ for the dimensionless variables $\tilde{p}$, $\tilde{L}$ and $\tilde{t}$. Figure 7.b shows the collapse of all the minimal pressure data of the parametric exploration previously present in Figure 6). Thus, $\tilde{p} = \frac{p-p_c}{j\tau L_0^2}$ is found to be a decreasing function of $\frac{\tau_{diff}}{\tau_{pv}}$. Interestingly, this non-dimensionalization may mask non-monotonic dependencies of $p_{min}$ on geometric design parameters. For instance, $p_{min}$ is non-monotonic with the width of the channels $w$ owing to the competitive effect of $w$ on $\frac{\tau_{diff}}{\tau_{pv}} \propto w^5$ and on $j\tau L_0^2 \propto w$.

## 2. Embolism length dynamics

We now try to predict the dynamics of the liquid length $L$. The speed of retraction of the interface is governed by the pressure gradient at the interface as outlined by equation 12. Slow diffusion confines the water injection, from the advancing of the interface, and pressure relaxation is confined near the interface, preventing the establishment of a stationary gradient over the full system. Absence of pressure relaxation on the terminal channels implies that all the water injected by the moving interface is removed from the system by the pervaporation over the length $\sim \sqrt{Dt}$ where the gradient establishes within time $t$. This hypothesis of limitation by the diffusion and the compensation of the water fluxes, only valid for $t \ll \tau_{diff} \sim \tau_{pv}$, can thus be translated into:

$$ S_c \dot{L} \approx -j \sqrt{Dt}, \quad (27) $$

which can be further solved into:

$$ L = L_0 - \frac{2}{3} \sqrt{D} \frac{j}{S_c} t^{3/2}, \quad (28) $$
<figure><img src="image_10.png" /><figcaption>FIG. 8. Embolism dynamics of resistive and deforming channel series collapse on a piecewise master curve. Non-dimensionalized embolism length is plotted as a function of time. In (a) all the data of figure 6 are represented and length is plotted at $t^*$, the time at which the pressure minimum is reached. Inset is the log-log scale view. Red continuous curve has equation $f(\tau_{pv}/\tau_{diff}) = 2/3 (\tau_{pv}/\tau_{diff})^{3/2}$. In (b) a representative selection of the length data is plotted. Red dashed lines plots the piecewise function of equation 30. This semi-analytical model appears to rather correctly fit the data.</figcaption></figure>

or:

$$ \left(1 - \frac{L}{L_0}\right) \sqrt{\frac{\tau_{diff}}{\tau_{pv}}} = \frac{2}{3} \left(\frac{t}{\tau_{pv}}\right)^{3/2}. \qquad (29) $$

Figures 8.a and 8.b show that this behavior for $L$ is indeed observed for $t < \tau_{pv}/4$ universally in all the trajectory obtained in the parametric exploration of the continuous model. For $t > \tau_{pv}/4$, $L(t)$ starts to deviate from this behavior. Indeed, the depression accumulated in the terminal channels starts to diffuse and pull on the interface so that the catching-up of pressure relaxation accelerates the movement of the interface. Surprisingly, during this phase, the embolism trajectories of the different series are still found to collapse in figure 8.a and b. During this transition regime, the interplay between the accumulated depressurization in distal channels and the advancing interface results in an accelerating front, even though pervaporation remains constant. This regime is thus driven by delayed elastic recoil and restored pressure communication. Eventually, when the liquid length $L$ has been significantly reduced, the pressure diffusion over the remaining length become very fast and no more pressure relaxation delays affect the embolism dynamics, such that quasistatic regime for $p$ and exponential decreasing of $L$ are found back. Based on the previous observations, we build a semi-analytical solution for $L(t)$, smoothly transitioning between the diffusion-limited regime (eq. 28) and the quasistatic regime (eq. 23), between $t = \tau^*$ and $t = 1.5\tau^*$:

$$ L(t) = \begin{cases} L_{dl} = L_0 \left(1 - \frac{2}{3} \frac{t^{3/2}}{\tau_{pv} \sqrt{\tau_{diff}}}\right) & \text{for } t < \tau^* \\ \sqrt{L_{dl}^2 \left(1 - \frac{t-\tau^*}{0.5\tau^*}\right) + L_{qs}^2 \frac{t-\tau^*}{0.5\tau^*}} & \text{for } \tau^* < t < 1.5\tau^* \\ L_{qs} = L^+ e^{-\frac{t-\tau^*}{\tau_{pv}}} & \text{for } t > 1.5\tau^*. \end{cases} \qquad (30) $$

Here, $L^+ \approx L_0 \left(1 - \frac{\tau^*}{\tau_{pv}} + \frac{4}{15} \frac{\tau^*^{5/2}}{\tau_{pv}^2 \tau_{diff}}\right)$ is the liquid length ensuring volume conservation if the embolism dynamics were to switch instantaneously from a diffusion-limited regime ($L_{dl}$) to a quasi-static regime ($L_{qs}$) at $\tau^*$ while conserving the water volume. As shown in Figure 8.b, this expression of $L(t)$ gives a rather correct prediction of the complex trajectory of the embolism depending on the design parameters of network. This analytical form provides a predictive tool for designing microfluidic systems with targeted embolism dynamics by tuning resistance and compliance.

# V. PERSPECTIVES

## A. Fluidic networks

Our results reveal how the interplay between compliance, viscous dissipation, and pervaporation can drive rich and nonlinear embolism dynamics in compliant microfluidic networks — dynamics which were previously unexplored.

Beyond embolism, nonlinear behaviors in related systems, such as memristor-like hydrodynamic resistance [19] or fluidic Braess paradoxes [18], suggest that pervaporating networks may also exhibit unexpected emergent phenomena.
In our system, nonlinearity arises not from flow-rate-dependent resistance but from the evolving functional topology of the network — as the advancing embolism front dynamically alters hydraulic connectivity. This mechanism bears similarities to frangible porous media, where structural changes dynamically influence transport [57, 58].

Extending this framework from one-dimensional to two-dimensional networks represents a promising next step, as branched and looped architectures are expected to introduce new collective behaviors and complex transport regimes. Previous investigations in networks with fast pressure diffusion evidenced embolism trajectories solely controlled by capillarity [48, 55]. In networks with slow pressure diffusion, our work suggests more subtle and intricate bubble propagation, and we can expect the induction of preferential embolism propagation paths, spatially localized collapse zones, or even front bifurcations, effects absent in one-dimensional networks. These investigations could be relevant both for fundamental studies of nonlinear dynamics [59] and for applications in adaptive fluidic networks [60].

On the experimental side, future studies will be guided by the design principles established here, leveraging recent pervaporation-driven methods [51, 61, 62]. Experimental setups using PDMS microchannels under negative pressure [63] offer promising avenues for validation. Furthermore, the possible triggering of cavitation events in regions far from the embolism front [64] could significantly modify embolism progression patterns, opening new questions in the physics of drying fluidic systems, and in particular for plant-inspired networks.

## B. Embolism propagation in plants

While our model is designed as a minimal physical description of embolism dynamics in linear fluidic networks, it may offer useful analogies for interpreting embolism propagation in plant vascular systems — particularly in species with simple or weakly branched architectures, such as Adiantum leaves or non-dividing stems [36, 37, 51].

In natural xylem networks, properties such as hydraulic resistance (R) and compliance (C) are difficult to access experimentally [65], and additional processes — including lateral water storage and interactions with living tissues — further complicate the picture [66]. Nonetheless, our model isolates a few key mechanisms — such as pressure diffusion delays and compliance-driven buffering — that could be relevant under certain conditions.

Embolism propagation in plants has been reported to occur intermittently and over a range of time scales [36, 67], the origin of which remains unclear. In this context, pressure diffusion delays could contribute an additional dynamical ingredient. The emergence of large transient pressure drops, as predicted in our study, might also help interpret sudden cavitation events observed in high-resolution optical or acoustic recordings.

While we model pit membranes as idealized constrictions with a capillary threshold, their actual structure is more complex. Intervessel pits in angiosperms contain porous and often deformable membranes [68], and gymnosperm pits include specialized torus-margo valves [16, 69]. These elements likely influence embolism propagation in ways not captured by our model. Yet, the idea that the pit geometry sets both a resistance and a capillary threshold — and that these two properties modulate embolism dynamics — remains qualitatively relevant.

Our findings suggest that large pressure drops may exceed capillary thresholds. Extending this approach to reticulate networks could help explore preferential embolism routes [37] or structural adaptations that maintain connectivity in vascular systems. Finally, the model opens the possibility that the parameter $\tau L_0^2$, which governs pressure field evolution in the continuous limit, could be of functional significance. Whether such a quantity is constrained or tuned across species remains an open question, but it may offer a useful perspective for interpreting hydraulic behavior in terms of global network properties.

## VI. CONCLUSION

In this article, we theoretically investigated air invasion in one-dimensional networks of compliant microchannels undergoing pervaporation. While such systems were previously explored in our earlier work [51], we extended the analysis by introducing a sufficiently large hydrodynamic resistance in the constrictions, which induces significant pressure and deformation inhomogeneities along the biomimetic veins.

We identified two key timescales that govern the dynamics of this system: the intrinsic pervaporation timescale $\tau_{pv}$ and the pressure diffusion timescale $\tau_{diff}$. When these timescales become comparable — or when $\tau_{diff} > \tau_{pv}$ — pressure within certain parts of the network, particularly near its distal end, can drop well below the capillary threshold $p_c$ required for embolism propagation.

By combining discrete network simulations and a continuum theoretical approach, we developed a predictive framework for the embolism dynamics. Within the continuous description, we showed that the system behavior could be fully captured by a single nondimensional parameter, $\tau_{diff}/\tau_{pv}$ (see Eqs. 25 and 26). Furthermore, we demonstrated that a diffusion-limited approximation captures essential features of the nonlinear dynamics and provides a simple, semi-analytical understanding of interface propagation.
Overall, the system analyzed here offers a fundamental illustration of how nonlinear behavior can emerge in fluidic networks owing to the coupling between pressure field and interfacial motion. Our model provides a basis for further studies of nonlinear elasto-fluidic dynamics, both in the context of sap flow in plants and in engineered microfluidic systems. Extending this approach to branched or looped architectures, or to networks exhibiting instabilities and multistable flow paths, could reveal richer emergent phenomena.

ACKNOWLEDGEMENT

This work has been supported by the French government, through the grant ANR-19-CE30-0010 (PHYSAP) and through the UniCA Investments in the Future project managed by the National Research Agency (ANR) with the reference number ANR-15-IDEX-01. The authors want to thank Benjamin Dollet and Philippe Marmottant for fruitful discussions, as well as Virgile Thiévenaz, José Torres-Ruiz, Eric Badel and Hervé Cochard.
# APPENDIX A: GLOBAL VOLUME CONSERVATION CONDITION IN THE CONTINUUM MODEL

We here demonstrate the equivalence of equation 13 with the global conservation of the water volume that must be fulfilled in the system. Noting $V$ the total volume of water in the system, global volume conservation can be expressed as follows:

$$ \dot{V} = -jL \qquad (31) $$

Let first note that $V(t) = \int_{L_0-L}^{L_0} S(x,t)dx$ and $S(x,t) = S_0 + cp(x,t)$, we have

$$ V(t) = S_0 L(t) + c \int_{L_0-L}^{L_0} p(x,t)dx \qquad (32) $$

To get the derivative of the volume, one has to use the Leibniz formula to account for the time-dependency of the integral boundaries. This gives:

$$ \begin{aligned} \dot{V} = S_0 \dot{L} + c \dot{L}_0 p(x = L_0, t) & \\ - c \frac{d}{dt}(L_0 - L)p(x = L_0 - L, t) & \\ + c \int_{L_0-L}^{L_0} \partial_t p(x, t) dx & \end{aligned} \qquad (33) $$

As $L_0$ is a constant, the first term vanishes and noting that $\frac{d}{dt}(L_0 - L) = -\dot{L}$, $p(x = L_0 - L) = p_c$ and $S_c = (S_0 + \frac{c p_c}{l})$, one can reexpress the global volume conservation as:

$$ \dot{V} = \dot{L} S_c + c \int_{L_0-L}^{L_0} \partial_t p(x) dx = -jL \qquad (34) $$

A change of variables $x' = \frac{x-L_0-L}{L}$ in the integral enables to simplify the expression and to find back equation 13.
[1] D. M. Wootton and D. N. Ku, Fluid mechanics of vascular systems, diseases, and thrombosis, Annual review of biomedical engineering 1, 299 (1999).

[2] S. Marbach, N. Ziethen, L. Bastin, F. K. Bäuerle, and K. Alim, Vein fate determined by flow-based but time-delayed integration of network architecture, eLife 12, e78100 (2023).

[3] M. Le Verge-Serandour and K. Alim, *Physarum polycephalum*: Smart Network Adaptation, Annual Review of Condensed Matter Physics 15, 263 (2024).

[4] L. Oyarte Galvez, C. Bisot, P. Bourrianne, R. Cargill, M. Klein, M. Van Son, J. Van Krugten, V. Caldas, T. Clerc, K.-K. Lin, F. Kahane, S. Van Staalduine, J. D. Stewart, V. Terry, B. Turcu, S. Van Otterdijk, A. Babu, M. Kamp, M. Seynen, B. Steenbeek, J. Zomerdijk, E. Tutucci, M. Sheldrake, C. Godin, V. Kokkoris, H. A. Stone, E. T. Kiers, and T. S. Shimizu, A travelling-wave strategy for plant-fungal trade, Nature 639, 172 (2025).

[5] T. Ogawa, S. Koyama, T. Omori, K. Kikuchi, H. de Maleprade, R. E. Goldstein, and T. Ishikawa, The architecture of sponge choanocyte chambers maximizes mechanical pumping efficiency, bioRxiv 10.1101/2024.02.22.581376 (2024), https://www.biorxiv.org/content/early/2024/02/22/2024.02.22.581376.full.pdf.

[6] K. Jensen, K. Berg-Sørensen, H. Bruus, N. Holbrook, J. Liesche, A. Schulz, M. Zwieniecki, and T. Bohr, Sap flow and sugar transport in plants, Rev. Mod. Phys. 88, 035007 (2016).

[7] K. Alim, N. Andrew, A. Pringle, and M. P. Brenner, Mechanism of signal propagation in *Physarum polycephalum*, Proceedings of the National Academy of Sciences 114, 5136 (2017).

[8] S. Liese, L. Mahadevan, and A. Carlson, Balancing efficiency and homogeneity of biomaterial transport in networks, EPL 135, 58001 (2021).

[9] L. Pereira, L. Kaack, X. Guan, L. D. M. Silva, M. T. Miranda, G. S. Pires, R. V. Ribeiro, H. J. Schenk, and S. Jansen, Angiosperms follow a convex trade-off to optimize hydraulic safety and efficiency, New Phytologist 240, 1788 (2023).

[10] M. Heil and A. L. Hazel, Fluid-structure interaction in internal physiological flows, Annu. Rev. Fluid Mech. 43, 141 (2011), publisher: Annual Reviews.

[11] H. Fallahi, J. Zhang, H.-P. Phan, and N.-T. Nguyen, Flexible microfluidics: Fundamentals, recent developments, and applications, Micromachines 10, 830 (2019).

[12] S. Battat, D. A. Weitz, and G. M. Whitesides, Nonlinear Phenomena in Microfluidics, Chemical Reviews 122, 6921 (2022).

[13] H. M. Xia, J. W. Wu, J. J. Zheng, J. Zhang, and Z. P. Wang, Nonlinear microfluidics: device physics, functions, and applications, Lab on a Chip 21, 1241 (2021).

[14] D. J. Case, J.-R. Angilella, and A. E. Motter, Spontaneous oscillations and negative-conductance transitions in microfluidic networks, Science Advances 6, eaay6761 (2020).

[15] M. Brandenbourger, A. Dangremont, R. Sprik, and C. Coulais, Tunable flow asymmetry and flow rectification with bio-inspired soft leaflets, Phys. Rev. Fluids 5, 084102 (2020).

[16] K. Park, A. Tixier, M. Paludan, E. Østergaard, M. Zwieniecki, and K. H. Jensen, Fluid-structure interactions enable passive flow control in real and biomimetic plants, Phys. Rev. Fluids 6, 123102 (2021).

[17] E. Jambon-Puillet, Dense array of elastic hairs obstructing a fluidic channel, arXiv (2025), arXiv:2501.01875 [cond-mat.soft].

[18] D. J. Case, Y. Liu, I. Z. Kiss, J.-R. Angilella, and A. E. Motter, Braess's paradox and programmable behaviour in microfluidic networks, Nature 574, 647 (2019).

[19] A. Martínez-Calvo, M. D. Biviano, A. H. Christensen, E. Katifori, K. H. Jensen, and M. Ruiz-García, The fluidic memristor as a collective phenomenon in elastohydrodynamic networks, Nat Commun 15, 3121 (2024).

[20] A. T. Bradley, F. Box, I. J. Hewitt, and D. Vella, Wettability-Independent Droplet Transport by *Bendotaxis*, Physical Review Letters 122, 074503 (2019).

[21] A. T. Bradley, I. J. Hewitt, and D. Vella, Bendocapillary instability of liquid in a flexible-walled channel, Journal of Fluid Mechanics 955, A26 (2023).

[22] C. Ushay, E. Jambon-Puillet, and P.-T. Brun, Interfacial flows past arrays of elastic fibers, Physical Review Fluids 8, 044001 (2023).

[23] D. Pihler-Puzović, P. Illien, M. Heil, and A. Juel, Suppression of complex fingerlike patterns at the interface between air and a viscous fluid by elastic membranes, Phys. Rev. Lett. 108, 074502 (2012), publisher: APS.

[24] A. Juel, D. Pihler-Puzović, and M. Heil, Instabilities in blistering, Annu. Rev. Fluid Mech. 50, 691 (2018).

[25] J.-F. Louf, J. Knoblauch, and K. H. Jensen, Bending and Stretching of Soft Pores Enable Passive Control of Fluid Flows, Physical Review Letters 125, 098101 (2020).

[26] J. Cappello, O. Du Roure, F. Gallaire, C. Duprat, and A. Lindner, Fiber Buckling in Confined Viscous Flows: An Absolute Instability Described by the Linear Ginzburg-Landau Equation, Physical Review Letters 129, 074504 (2022).

[27] G. Guyard, F. Restagno, and J. D. McGraw, Elastohydrodynamic relaxation of soft and deformable microchannels, Phys. Rev. Lett. 129, 204501 (2022).

[28] M. V. Paludan, B. Dollet, P. Marmottant, and K. H. Jensen, Elastohydrodynamic interactions in soft hydraulic knots, J. Fluid Mech. 984, A55 (2024).

[29] H. Garg, P. G. Ledda, J. S. Pedersen, and M. Pezzulla, Passive Viscous Flow Selection via Fluid-Induced Buckling, Physical Review Letters 133, 084001 (2024).

[30] M. T. Tyree and M. H. Zimmermann, *Xylem Structure and the Ascent of Sap* (Springer Science & Business Media, 2013).

[31] A. D. Stroock, V. V. Pagay, M. A. Zwieniecki, and N. Michele Holbrook, The Physicochemical Hydrodynamics of Vascular
Plants, Annu. Rev. Fluid Mech. **46**, 615 (2014).

[32] E. Katifori, The transport network of a leaf, *Comptes Rendus Physique **19***, 244 (2018).

[33] B. Choat, S. Jansen, T. J. Brodribb, H. Cochard, S. Delzon, R. Bhaskar, S. J. Bucci, T. S. Feild, S. M. Gleason, U. G. Hacke, *et al.*, Global convergence in the vulnerability of forests to drought, *Nature **491***, 752 (2012).

[34] B. Choat, T. J. Brodribb, C. R. Brodersen, R. A. Duursma, R. López, and B. E. Medlyn, Triggers of tree mortality under drought, *Nature **558***, 531 (2018).

[35] J. M. Torres-Ruiz, H. Cochard, S. Delzon, T. Boivin, R. Burlett, M. Cailleret, D. Corso, C. E. Delmas, M. De Caceres, A. Diaz-Espejo, *et al.*, Plant hydraulics at the heart of plant, crops and ecosystem functions in the face of climate change, *New Phytologist **241***, 984 (2024).

[36] T. J. Brodribb, D. Bienaimé, and P. Marmottant, Revealing catastrophic failure of leaf networks under stress, *Proc. Natl. Acad. Sci. **113***, 4865 (2016).

[37] T. J. Brodribb, M. Carriqui, S. Delzon, and C. Lucani, Optical measurement of stem xylem vulnerability, *Plant Physiology **174***, 2054 (2017), publisher: American Society of Plant Biologists.

[38] H. Cochard, F. Froux, S. Mayr, and C. Coutand, Xylem wall collapse in water-stressed pine needles, *Plant Physiol. **134***, 401 (2004).

[39] Y. J. Zhang, F. E. Rockwell, A. C. Graham, T. Alexander, and N. M. Holbrook, Reversible leaf xylem collapse: a potential 'circuit breaker' against cavitation, *Plant Physiol. **172***, 2261 (2016).

[40] A. R. Chin, P. Guzmán-Delgado, S. C. Sillett, L. P. Kerhoulas, A. R. A. A. R. McElrone, and M. A. Zwieniecki, Tracheid buckling buys time, foliar water uptake pays it back: Coordination of leaf structure and function in tall redwood trees, *Plant Cell Environ. **45***, 2607 (2022).

[41] T. D. Wheeler and A. D. Stroock, The transpiration of water at negative pressures in a synthetic tree, *Nature **455***, 208 (2008).

[42] X. Noblin, L. Mahadevan, I. A. Coomaraswamy, D. A. Weitz, N. M. Holbrook, and M. A. Zwieniecki, Optimal vein density in artificial and real leaves, *Proc. Natl. Acad. Sci. **105***, 9140 (2008), publisher: National Acad Sciences.

[43] C. Duan, R. Karnik, M.-C. Lu, and A. Majumdar, Evaporation-induced cavitation in nanofluidic channels, *Proc. Natl. Acad. Sci. **109***, 3688 (2012), publisher: National Acad Sciences.

[44] O. Vincent, P. Marmottant, P. A. Quinto-Su, and C.-D. Ohl, Birth and growth of cavitation bubbles within water under tension confined in a simple synthetic tree, *Phys. Rev. Lett. **108***, 184502 (2012), publisher: APS.

[45] O. Vincent, D. A. Sessoms, E. J. Huber, J. Guioth, and A. D. Stroock, Drying by Cavitation and Poroelastic Relaxations in Porous Media with Macroscopic Pores Connected by Nanoscale Throats, *Phys. Rev. Lett. **113***, 134501 (2014).

[46] M. A. Bruning, M. Costalonga, J. H. Snoeijer, and Marín, Turning drops into bubbles: cavitation by vapor diffusion through elastic networks, *Phys. Rev. Lett. **123***, 214501 (2019), publisher: APS.

[47] B. Dollet, J. F. Louf, M. Alonzo, K. H. Jensen, and P. Marmottant, Drying of channels by evaporation through a permeable medium, *J. R. Soc. Interface **16***, 20180690 (2019).

[48] B. Dollet, K. N. C. Encarnación, R. Gautier, and P. Marmottant, Drying by pervaporation in elementary channel networks, *J. Fluid Mech. **906***, A6 (2021).

[49] O. Vincent, T. Tassin, E. J. Huber, and A. D. Stroock, Tunable transport in bidisperse porous materials with vascular structure, *Phys. Rev. Fluids **9***, 064202 (2024), publisher: American Physical Society.

[50] L. Keiser, P. Marmottant, and B. Dollet, Intermittent air invasion in pervaporating compliant microchannels, *J. Fluid Mech. **948***, A52 (2022), publisher: Cambridge University Press.

[51] L. Keiser, B. Dollet, and P. Marmottant, Embolism propagation in *Adiantum* leaves and in a biomimetic system with constrictions, *Journal of The Royal Society Interface **21***, 20240103 (2024).

[52] H. Bruus, *Theoretical Microfluidics* (Oxford University Press, 2008).

[53] H. Wong, S. Morris, and C. Radke, Three-dimensional menisci in polygonal capillaries, *Journal of Colloid and Interface Science **148***, 317 (1992).

[54] The latter condition ensures the continuity, between the discrete and the continuum model, of the volume injection in the remaining wet channel series due to the movement of the interface. By definition of the continuum approximation, noting $\langle \cdot \rangle$ the average for many channel draining observations, the water volume flux coming from the embolization of the channels is: $\langle S_c l / \Delta t_{tot} \rangle$ with $S_c = S_0 + C p_c / l$. Given that $\langle l / \Delta t_{tot} \rangle = \dot{L}$, we indeed find that we have to impose the boundary condition $p = p_c$ at the interface to conserve the volume flux between the discrete and continuous model.

[55] K. N. C. Encarnación, P. Marmottant, and B. Dollet, Pervaporation-induced drying in networks of channels of variable width, *Microfluid. Nanofluid. **25***, 71 (2021).

[56] Note that the continuous model was solved using a simplified version of the pervaporation rate, only keeping the linear dependence of $j$ with the $w$ and neglecting the corrective terms of Eq. 1.

[57] N. J. Derr, D. C. Fronk, C. A. Weber, A. Mahadevan, C. H. Rycroft, and L. Mahadevan, Flow-Driven Branching in a Frangible Porous Medium, *Phys. Rev. Lett. **125***, 158002 (2020).

[58] A. Zareei, D. Pan, and A. Amir, Temporal evolution of erosion in pore networks: From homogenization to instability, *Phys. Rev. Lett. **128***, 234501 (2022).

[59] S. Fancher and E. Katifori, Mechanical response in elastic fluid flow networks, *Phys. Rev. Fluids **7***, 013101 (2022).

[60] A. Winn and E. Katifori, Operating principles of peristaltic pumping through a dense array of valves, *Journal of Fluid Mechanics **989***, A18 (2024).

[61] P. Bacchin, J. Leng, and J. B. Salmon, Microfluidic evaporation, pervaporation, and osmosis: From passive pumping to solute concentration, *Chem. Rev. **122***, 6938 (2022), publisher: ACS Publications.
[62] F.-X. Gauci, L. Jami, L. Keiser, C. Cohen, and X. Noblin, Channel deformations during elastocapillary spreading of gaseous embolisms in biomimetic leaves, *accepted in Interface Focus*, arXiv preprint arXiv:2501.05832 (2025).

[63] M. Bruning, M. Costalonga, J. Snoeijer, and A. Marin, Turning Drops into Bubbles: Cavitation by Vapor Diffusion through Elastic Networks, Phys. Rev. Lett. **123**, 214501 (2019).

[64] H. Pingulkar, C. Ayela, and J.-B. Salmon, Pervaporation-driven electrokinetic energy harvesting using poly(dimethylsiloxane) microfluidic chips, Lab on a Chip **24**, 5328 (2024), arXiv:2411.15226 [cond-mat].

[65] C. J. Blackman and T. J. Brodribb, Two measures of leaf capacitance: insights into the water transport pathway and hydraulic conductance in leaves, Functional Plant Biology **38**, 118 (2011), publisher: CSIRO Publishing.

[66] T. Hölttä, H. Cochard, E. Nikinmaa, and M. Mencuccini, Capacitive effect of cavitation in xylem conduits: results from a dynamic model, Plant, Cell & Environment **32**, 10 (2009).

[67] V. Thievenaz, J.-M. Torres-Ruiz, H. Cochard, E. Badel, and X. Noblin, Unpublished data, .

[68] Y. Zhang, L. Pereira, L. Kaack, J. Liu, and S. Jansen, Gold perfusion experiments support the multi-layered, mesoporous nature of intervessel pit membranes in angiosperm xylem, New Phytologist (2024), publisher: Wiley Online Library.

[69] J. Pittermann, J. S. Sperry, U. G. Hacke, J. K. Wheeler, and E. H. Sikkema, Torus-Margo Pits Help Conifers Compete with Angiosperms, Science **310**, 1924 (2005).