Benjamin Dollet

Kennedy Nexon

Chagua Encarnación

Romain Gautier

Philippe Marmottant

Drying by pervaporation in elementary channel networks

. In branched networks, experiments reveal velocity discontinuities of the menisci as they split at the nodes. In loops, it is found that the drying rate depends on the number of menisci bounding a given connected water region; when there are two such menisci, a prediction of the dynamics of each of them is proposed, based on the pervaporation-induced hydrodynamics inside the channels. Experiments and model predictions compares favourably for the global drying rate. Some deviations are found for the dynamics of individual menisci, which are ascribed to the sensitivity of the dynamics to small fluctuations in wetting conditions.

Introduction

Pervaporation is the process of drying of a liquid region by mass transport through a permeable medium. It occurs naturally in vascular plant leaves: water (or more precisely sap) is flowing in a networks of veins, diffuses out of the veins through the leaf tissue and then evaporates at the surface of the leaf. This fascinating network of veins, called venation, can be very simple, with only one vein per leaf for Selaginella species, or can be very complex, with multiple connections between veins for flowering plants. A first aspect of this complexity is the hierarchy of vein sizes with main veins, secondary veins, tertiary veins, and so forth. Another important aspect is the connectivity of the network, on which we focus here. When starting from the main vein where water is entering the leaf, veins can either: (i) split into several branches, (ii) reconnect, creating loops and a net-like reticulation of the network. In normal conditions water is flowing through the whole network. However, under severe drought, a cavitation bubble can appear and grow [START_REF] Tyree | Vulnerability of xylem to cavitation and embolism[END_REF][START_REF] Cochard | Cavitation in trees[END_REF], the remaining water in the network then dries out, eventually disabling the leaf function. This mechanism may threaten a growing proportion of plants and trees in the context of global warming [START_REF] Choat | Triggers of tree mortality under drought[END_REF]. In a different context, pervaporation through thin slabs of a water-permeable polymer, polydimethylsiloxane (PDMS) is used in so-called microevaporators in chemical engineering. It is useful because it ensures a homogeneous concentration of solutes or colloids and therefore a proper study of their phase diagrams [START_REF] Walker | An evaporation-based microfluidic sample concentration method[END_REF][START_REF] Eijkel | Osmosis and pervaporation in polyimide submicron microfluidic channel structures[END_REF][START_REF] Leng | Microevaporators for kinetic exploration of phase diagrams[END_REF][START_REF] Merlin | Microfluidic-assisted growth of colloidal crystals[END_REF][START_REF] Ziane | Drying with no concentration gradient in large microfluidic droplets[END_REF][START_REF] Ziemecka | Hydrogen peroxide concentration by pervaporation of a ternary liquid solution in microfluidics[END_REF]. However, despite the importance of these biological and chemical applications, pervaporation-induced flows in channel networks are little known, in contrast with the well-established field of pressure-driven flows in porous media.

To bridge this gap with a physical standpoint, our approach is to understand the drying using artificial channels in a permeable material, following the concept of biomimetic leaves incorporating a man-made microfluidic network in a porous material [START_REF] Noblin | Optimal vein density in artificial and real leaves[END_REF]. Geometry and physical properties of the channels in the network are then well controlled. During drying, air invades the microfluidic network initially filled with water. The interfaces (menisci) separating air and remaining liquid, travel along all the branches of the network, and eventually either reach end points or merge with other ones.

The dynamics of moving menisci and expelled liquid has been observed in 2D network of microchannels in literature, but with different driving mechanisms: first, when a source of air pressure pushes liquid plugs out of a microfluidic network towards openings, for instance to study air-way reopening in lung [START_REF] Song | Local interactions and the global organization of a two-phase flow in a branching tree[END_REF][START_REF] Song | The air-liquid flow in a microfluidic airway tree[END_REF][START_REF] Mamba | Dynamics of a liquid plug in a capillary tube under cyclic forcing: memory effects and airway reopening[END_REF]. Note the behaviour of these liquid plugs is similar to long droplets trapped in a network, causing them to potentially divide [START_REF] Link | Geometrically mediated breakup of drops in microfluidic devices[END_REF][START_REF] Salkin | Passive breakups of isolated drops and one-dimensional assemblies of drops inmicrofluidic geometries: experiments and models[END_REF]. There is also a similarity with the drainage of a large amount of liquid when pushed by another liquid [START_REF] Armstrong | Interfacial velocities and capillary pressure gradients during Haines jumps[END_REF]. Second, the motion of menisci in networks has been studied when wetting drives the liquid, in the case of the spontaneous imbibition of a 2D network by a wetting liquid [START_REF] Sadjadi | Scaling theory for spontaneous imbibition in random networks of elongated pores[END_REF]. Here, the specificity of our study is that water constantly pervaporates all over the network, driving the flow.

The goal of this article is to provide a physical study on the dynamics of drying by pervaporation in channel networks. Such a dynamics has been previously studied and understood in single, straight channels [START_REF] Noblin | Optimal vein density in artificial and real leaves[END_REF][START_REF] Ziane | Drying with no concentration gradient in large microfluidic droplets[END_REF][START_REF] Dollet | Drying of channels by evaporation through a permeable medium[END_REF]. In this paper, we build up on this knowledge to study the most elementary network topologies: simple branches and loops; our aim is to provide a quantitative picture in these simple cases, which could be the " building block" for more complex networks such as those encountered in natural leaves.

Materials and methods

Materials

In all our designs, there is a straight entrance channel into which drying starts. We added to all designs a circle of radius 1 mm at the root of the entrance channel. In plant leaves, embolism starts always in larger conduits (Brodribb et al. 2016a,b). The larger conduits are located on the larger veins, which occur at the base of the leaf in the petiole. Artificial embolism has therefore been proposed by injecting gas in the petiole of plant leaves [START_REF] Hochberg | Visualizing embolism propagation in gas-injected leaves[END_REF], to mimick the development of embolism in natural conditions. In our setup, the entrance channel thus mimics the main vein of leaves.

To study the effect of branching (without reconnections), we start with four designs with a single node where the entrance channel arrives and from where several straight branches, either two, three, four or five, emerge (Fig. 1). We shall henceforth call these designs single-node networks. All branches are of length L b = 4.0 mm. For each singlenode network, all channels (entrance channel and branches) have an equal angular distance from the node.

We then study branched networks with a larger number and a higher hierarchy of nodes, but at which each incoming branch subdivides in only two branches. First, we designed such a "tree" with three levels of nodes (Fig. 2a). At each node, the two branches are separated by an angle of 60 • . Contrary to Fig. 1, not all the branches have the same length, to determine whether length variations have an effect on the drying kinetics. Finally, we study a natural vein network, namely the leaf of the fern Adiantum (Brodribb et al. 2016a). To do this, we adapted a mask drawn from the real leaf [START_REF] Bienaimé | Embolies dans les plantes[END_REF] and modified the channel widths such that they all have the same width w = 100 µm (Fig. 2b), which is the common width of all the branched designs in the current study.

To study simple loops, we created four designs where the entrance channel splits at a first node, henceforth called entrance node, in two channels making a circle of diameter 2R = 4.3 mm, and reconnecting at a second node, henceforth called exit node, from which a single straight terminal channel emerges. One of the loop is symmetric, i.e. the distance between both nodes is equal in the two channels connecting them (Fig. 3a), and the three other ones are asymmetric, with a shortest angular distance β between the nodes, measured from the centre of the circle, equal to 120 • (Fig. 3b andd) or 60 • (Fig. 3c); this angle is obviously 180 • in Fig. 3a. The terminal channel has a length L t = 3.0 mm in Fig. 3a, b and c, and 39 mm in Fig. 3d, to determine whether L t has an effect on the drying dynamics.

We fabricated the channels based on all these designs as described in [START_REF] Dollet | Drying of channels by evaporation through a permeable medium[END_REF]. We created masks from the designs using high-resolution printing. The channels were molded in PDMS using standard soft lithography techniques. From the masks, we created a mold in a photoresistive material (SU8) on a silicon wafer. We mixed liquid PDMS (Sylgard 184, Dow Corning) with a curing agent in volumetric proportions 9:1. This mixture was degassed then spin coated on the mold at 300 revolutions per minute (rpm) for 15 s then 900 rpm for 40 s, to create an imprint of small thickness, and baked at 65 • C for one hour to reticulate the PDMS. This imprint was then bonded to a glass slide as follows. The bare glass slide was exposed for 30 s to a plasma, then the thin, flexible PDMS imprint was immediately deposited on top of the plasma-activated glass surface, carefully avoiding folding the PDMS sheet. This simple procedure was sufficient to firmly bond PDMS and glass. We measured the channel thickness h and the PDMS thickness H using an interferometer. This is crucial because these thicknesses presented significant variations between different designs, despite the use of the same spin-coating parameters. The values of the thicknesses are tabulated in Tab. 1. To speed up the measurements, we gathered the four single-node networks on one slide, the tree and Adiantum networks on a second one, and all loops on a third one, at a distance large enough that diffusive interactions remain negligible. After Fig. 3 from [START_REF] Noblin | Optimal vein density in artificial and real leaves[END_REF], it is the case as soon as channels are separated by a distance of order 5δ, with δ = H -h. In our case, δ ≈ 60 µm (Tab. 1), hence the distance between channels should be larger than 0.3 mm. In practice, the distance between different designs is larger than 1 mm.

Methods

The experiments are performed as described in [START_REF] Dollet | Drying of channels by evaporation through a permeable medium[END_REF]; we only repeat here the main steps. We opened the channels by manually cutting through the PDMS single-node networks (Fig. 1) tree and Adiantum (Fig. 2) loops (Fig. 3 layer across the circle with a scalpel at the junction between the channels and the round part (these cuts are clearly visible in Figs. 1 and2). Once the channels were water filled, we placed them under dry atmosphere by imposing a constant dry air flux above the PDMS layer. Each channel progressively dried out, as menisci separating water from air advanced downstream through the channel from the roots to farther branches, until reaching the ends of the terminal branches. The motion of the menisci was imaged with a CCD camera.

To measure the location of the menisci as a function of time, we used the reslice operation of the freeware ImageJ and a home-made Matlab script. For each design, we draw a series of segmented lines following the middle of each channel, and we generate spatiotemporal diagrams along each line, for which we get the time evolution of the location of the menisci. Obtaining the velocity of the meniscus required to take the discrete derivative of the length of the channel filled with water as a function of time. We chose to approximate the data for each meniscus, between its creation and annihilation (e.g. at a node), by a polynomial function of degree n 9. We chose the degree of the polynomial which gives the best coefficient of determination R 2 between the experimental data and the interpolating polynomial. We performed such an operation on each meniscus as long as it can be unambiguously defined between splitting or annihilation events (see next section).

Qualitative description of the experimental drying dynamics

Branched networks

In branched networks, the single meniscus in the entrance channel first moves until it splits at the entrance node, as shown in Figs. 4. In the single-node networks, each meniscus then moves until it annihilates against the terminal wall of each terminal branch. The dynamics in each branch has the same speed within a few % (Fig. 4a). In the tree and Adiantum (Fig. 4b), there is a "cascade" of consecutive splitting events as deeper nodes are reached by the different menisci, until the terminal branches are reached.

Loops

Drying in loops proceeds as follows. First, a single meniscus moves in the entrance channel until it splits in two menisci at the entrance node. Then we distinguish two cases. For the symmetric loop (Fig. 5a), the two menisci move almost symmetrically, until they merge at the exit node; the resulting meniscus in the terminal channel moves until it annihilates on the terminal wall of this channel, which ends the drying process. For the asymmetric loops (Fig. 5b), where the loop has a short and a long arm, the meniscus in the short arm reaches first the exit node, where it splits in two. We have then two disconnected water regions and three menisci: one water region is in the terminal channel, while the other, bounded by two menisci, is in the long arm, until both regions dry.

Theory

Drying dynamics for a single channel

In this Section, we recall the main ingredients to model the drying dynamics by pervaporation in a single channel; see [START_REF] Dollet | Drying of channels by evaporation through a permeable medium[END_REF] for full details. This case is important because the drying models in networks, derived in Secs. 4.2 and 4.3, elaborate on this simplest case, and also because it describes drying in each terminal branch of the various networks.

According to the experimental configurations, we consider a semi-infinite channel of width w and height h between a glass slide and a PDMS layer of thickness H (Fig. 1e), filled with water over a length L between the channel extremity and an air/water meniscus. The outer air is at controlled relative humidity RH = c out a /c sat a with c out a the water vapour concentration in the outer air and c sat a the water vapour concentration in saturated air; in our experiments, RH = 0. The channels dry because of water diffusion from the water-filled part of the channel to the outer air through water-permeable PDMS. We showed in [START_REF] Dollet | Drying of channels by evaporation through a permeable medium[END_REF] that there are two contributions to the rate of decrease of the liquid volume: one by diffusion from the water-PDMS interface (Q ), and one by evaporation from the meniscus (Q g ). Assuming that L is much greater than H and w, a condition generally met in practice because microfluidic channels are very slender, the water concentration profile in the PDMS between the water-filled part of the channel and the outer air is almost invariant along the channel length; hence, Q = q L with q the diffusive flux per unit length along the channel. We henceforth consider Q and Q g as volumetric fluxes, and not molar fluxes as in [START_REF] Dollet | Drying of channels by evaporation through a permeable medium[END_REF], for simplicity; the main conclusions are unaffected by this choice. Equating the variation of volume of the water-filled part of the channel to the diffusive flux yields hw L = -q L + Q g , hence:

L + L + L g τ = 0, (4.1)
with:

L g = Q g q , τ = hw q . (4.2)
Eq. ( 4.1) has the solution:

L(t) = [L(t = 0) + L g ]e -t/τ -L g , (4.3)
for t τ ln[L(t = 0)/L g ], which is the time where the channel fully dries. In Dollet et al.

(2019), we derived analytical predictions for both q and Q g , hence for L g and τ , by solving first the two-dimensional diffusion problem in the cross-section of the water-filled part of the channel, then the three-dimensional diffusion problem in the PDMS and the air-filled part of the channel away from the meniscus. These calculations yield:

L g = αD a D P hw w δ + ξ , (4.4)
and:

τ = ρ D P M C sat P (1 -RH) hw w δ + ξ , (4.5)
with α = 0.03 the Henry constant quantifying the water affinity in PDMS [START_REF] Harley | Thermodynamic study in dynamic water vapor sorption in sylgard-184[END_REF], D a = 2 × 10 -5 m 2 /s the diffusivity of water vapour in air, D P = 10 -9 m 2 /s the diffusivity of water in PDMS [START_REF] Watson | The behaviour in water in poly(dimethylsiloxane)[END_REF], ρ = 10 3 kg/m 3 the density of water, M = 0.018 kg/mol the molar mass of water, C sat P = 40 mol/m 3 the saturation concentration of water in PDMS [START_REF] Randall | Permeation-driven flow in poly(dimethylsiloxane) microfluidic devices[END_REF], δ = H -h the PDMS thickness between the channel top and the outer air, and ξ a dimensionless factor quantifying diffusion between the channel side walls and the outer air; under reasonable geometrical assumptions quantified in [START_REF] Dollet | Drying of channels by evaporation through a permeable medium[END_REF], it equals:

ξ = 2 π ln (H + δ)h δ 2 + H δ ln H + δ h .
In anticipation to what follows, it is worth reminding that the evaporative flux Q g across the meniscus is related to water vapour in the air inside the channel: by Fick's law, Q g = -hwD a dc a /dx | x =0 , with c a the water vapour concentration field and x the distance from the meniscus. The spatial dependence c a (x ) is itself related to the water transfer from the channel to the outer atmosphere through PDMS, hence it results from the coupled diffusion of water in air and in PDMS. In [START_REF] Dollet | Drying of channels by evaporation through a permeable medium[END_REF] we showed that d 2 c a /dx 2 -c a /L 2 g = 0, with L g given by (4.4). Hence, L g can be interpreted as the typical distance over which air (and also PDMS) dries away from the meniscus (Fig. 1f). This distance is much larger than the typical cross-section dimensions: from (4.4), L g ≈ M √ hδ with an amplification factor:

M = αD a D P , (4.6)
which equals 30 with the aforementioned values for the different parameters. It is much larger than one, mostly because the diffusivity of water through PDMS is much less efficient than in air (D P D a ), which "delays" the drying of air inside the channel. It implies that the evaporation flux per unit area of the meniscus is much stronger (by a factor M ) than the liquid pervaporation flux per unit area of the PDMS/channel interface. Yet, the contribution of liquid pervaporation to the total drying rate dominates

if Q > Q g or if L > L g ≈ M √
hδ, which is the most common case in our experiments. Increasing L g implies weaker concentration gradients, which may at first sight lower Q g and contradict the very definition (4.2) of L g . This apparent contradiction is solved by the fact that increasing L g also implies increasing the prefactor relating Q g to the concentration gradient, relative to the prefactor relating q to the concentration gradient.

The role of the evaporation through the meniscus into the air is somewhat reminiscent of the role of fins in heat transfer, which enhance the heat flux from a heated surface [START_REF] Bejan | Heat Transfer[END_REF]. In some sense, evaporation through the air channel acts as an extra effective "fin" (in the sense of mass transfer, not of heat transfer) which adds up to the pervaporation flux from the PDMS/water interface.

Drying dynamics for branched channels

In branched channels, our experiments always begin with a single meniscus inside the entrance channel; as soon as the meniscus meets the first node of the network, it splits into two or more menisci, depending on the number of branches; and so forth for each of these menisci if they reach further nodes (Fig. 4).

We start our discussion with the single meniscus in the entrance channel at the beginning. Because the water-filled part spans multiple branches, the water concentration profile in the PDMS may differ from the single-channel case. More precisely, in the vicinity of a node, different branches are close enough to "compete" for drying because they are coupled in the diffusion problem. However, as shown in [START_REF] Noblin | Optimal vein density in artificial and real leaves[END_REF], diffusive interactions between two neighbouring channels decay exponentially with H as the characteristic length scale. Hence, except if branches start from a node with a small separation angle, diffusive interactions thus fade away at a distance of order H along each branch from a node. Therefore, provided the length of each branch between two consecutive nodes (or between a node and an extremity) is much larger than H, we can neglect the edge effect due to the nodes and extremities and consider, as in the case of a single channel, that the concentration profile between each water-filled branch and the outer air is invariant along the branch length, and similar to the single-channel case. In this case, the diffusive flux from the water-PDMS interface and the outer air still writes: Q = q L, with L the total length of the water-filled branches ahead from the meniscus (Fig. 6a). Hence, the single meniscus in the entrance channel is predicted to obey (4.1) with solution (4.3).

L 0 L 1 L 2 L l (a) (b) (c) 0 x 0 x 1 x 2 L 0 L 1 L 2
Let us now consider one meniscus inside one of the branches, obtained after a splitting event. From the previous discussion, we still have Q = q L, with L the total length of the water-filled branches ahead from the meniscus. For simplicity, we shall henceforth call this quantity L the "water length". However, if the distance travelled by the meniscus from the node where it has been created (Fig. 6b) is smaller than, or of order of magnitude L g , the water concentration profile behind the meniscus, responsible for the contribution Q g (see Sec. 4.1), is affected by the presence of the node. In other words, menisci issued from a given splitting event interact over a range of order L g , which is much larger than H or w (see Sec. 4.1). Hence in the case L g , we cannot use the single-channel prediction for Q g .

In the Appendix, we derive a prediction for the diffusive flux from the meniscus accounting for this interaction. We show in Sec. 5.1.1 that the discrepancy between this prediction and the single-channel prediction remains small in the parameter range of our experiments. Hence, to simplify the forthcoming discussion of our experiments, we shall simply use (4.1) as a prediction with L the total length of the water-filled branches ahead from the meniscus under consideration, unless explicitly stated.

Drying dynamics for loops

The key difference between branched channels without loops on one hand, and loops on the other hand, is the following. In branched channels, each connected water region is bounded by one meniscus. On the contrary, once the meniscus has split at the entrance node of the designs with loops, there is a connected water region bounded by two menisci (snapshots at t = 15.6 and 39 min in Fig. 5a, and at t = 13.4 min in Fig. 5b). However, the model developed in Sec. 4.2 predicts the rate of drying of a given connected water region, but not how such a rate determines the individual velocity of two (or more generally multiple) menisci bounding that region.

To go one step further, we must delve into the details of the hydrodynamics driven by pervaporation. We start with the simplest case, namely the single channel with a water region of length L bounded by one meniscus and one terminal wall, a case for which pervaporation-induced flow is well documented [START_REF] Verneuil | Permeation-induced flows: Consequences for silicone-based microfluidics[END_REF][START_REF] Randall | Permeation-driven flow in poly(dimethylsiloxane) microfluidic devices[END_REF][START_REF] Leng | Microevaporators for kinetic exploration of phase diagrams[END_REF]. We denote by x the streamwise axis, with x = 0 at the wall and x = L at the meniscus, and u(x, t) the velocity, averaged over the cross-section. Since q is the flux of water leaving the channel per unit length, water conservation inside the channel (assuming incompressible flow) implies that ∂u/∂x = -q /hw. To integrate this equation, we use the classical impermeability condition at the terminal wall:

u(0, t) = 0, (4.7)
to get:

u = - q hw x. (4.8)
It is worth noting that u(L, t) = L(t), because of the evaporation flux Q g leaving the meniscus; this is an illustration of Stefan flows. The correct relationship is:

u(L, t) = L + Q g hw . (4.9)
Moreover, the impermeability assumption (4.7) is, strictly speaking, wrong because of the "leakiness" of the PDMS walls. However, because of the large value of the amplification factor (4.6), the water velocity across the terminal wall remains small in comparison to Q g /hw, which justifies using (4.7) as a first approximation. This also enables to neglect the velocity across the walls everywhere, hence we can use a no-slip boundary condition there, and we can also neglect the transverse pressure gradients compared to the longitudinal pressure gradient. Therefore, we can use the Poiseuille law relating the average velocity with the pressure gradient for a laminar flow in a rectangular channel:

u = - S η ∂p ∂x , (4.10)
with the section S equal to (see e.g. [START_REF] Bruus | Theoretical Microfluidics[END_REF]):

S = 1 12 h 2 - 16 π 5 h 3 w ∞ n=0 1 (2n + 1) 5 tanh (2n + 1)πw 2h . (4.11)
The boundary condition for the pressure field comes from Laplace pressure jump across the meniscus:

p a -p(L, t) = γκ, (4.12)
with p a the atmospheric pressure, γ the surface tension and κ the curvature, counted positive if the centre of curvature is on the air side. The curvature depends on the width and height of the channel, and on the contact angles of water on both glass and PDMS; it is positive with our convention because both contact angles are smaller than 90 • , as can be seen from the meniscus curvature on snapshots. From (4.8) and (4.10), we thus get the pressure field in a single channel:

p(x, t) = p a -γκ - ηq 2hwS (L 2 -x 2 ). (4.13)
We now consider a connected water region with two menisci. If it spans a single channel, in particular the long arm of a loop as in Fig. 5b (see the snapshot at t = 47.5 min), its drying dynamics is very similar to the single-channel case discussed in Sec. 4.1, except that the water flux Q g across a given meniscus must be counted twice. Hence instead of (4.1), we have:

L + L + 2L g τ = 0, (4.14)
provided that the distance travelled by each meniscus from the node where it has been created is smaller than, or of order of magnitude L g (see the discussion of Sec. 4.2). Under this assumption, the dynamics of each of the two menisci must be similar and they must move at the same velocity.

If a connected water region with two menisci spans three "branches" connecting at a node, as in Fig. 5b (see the snapshot at t = 13.4 min), we must consider each branch. For the terminal channel of length L 0 , we simply adapt (4.7) and (4.8): u 0 (L 0 , t) = 0 and ∂u 0 /∂x 0 = -q /hw (see Fig. 6c for the orientation of the axis); in particular, we get U 0 ≡ u 0 (0, t) = q L 0 /hw. For each of the two parts between the node and one meniscus, ∂u i /∂x i = -q /hw (for i = 1, 2), hence:

u i (x i , t) = U i - q hw x i , (4.15)
where U i is the velocity in channel i at the node. One relation between the two unknowns U 1 and U 2 comes from the flux conservation at the node:

U 0 + U 1 + U 2 = 0. (4.16)
To get a second relation, we compute the pressure difference between the node and each of the two menisci:

Li 0 ∂p i ∂x i dx i = p a -γκ -p 0 ,
for i = 1, 2, after (4.12), and denoting p 0 the pressure at the node. We also have:

Li 0 ∂p i ∂x i dx i = - η S Li 0 u i dx i = - η S - q 2hw L 2 i + U i L i ,
from (4.10) and the expression (4.15) of the velocity, whence the second relation between U 1 and U 2 :

- q 2hw L 2 1 + U 1 L 1 = - q 2hw L 2 2 + U 2 L 2 .
(4.17)

Solving the system (4.16) and (4.17), we get:

U 1 = q hw 1 2 (L 1 -L 2 ) - L 2 L 1 + L 2 L 0 , U 2 = - q hw 1 2 (L 1 -L 2 ) + L 1 L 1 + L 2 L 0 .
We can finally predict the velocity of the menisci. Inserting U 1 and U 2 in (4.15), taking into account that Li = u i (L i , t) -Qg hw from (4.9) and using (4.2), we get:

L1 = - 1 τ 1 2 (L 1 + L 2 ) + L 2 L 1 + L 2 L 0 + L g , L2 = - 1 τ 1 2 (L 1 + L 2 ) + L 1 L 1 + L 2 L 0 + L g . (4.18)
These predictions deserve some comments. First, Li < 0, which is of course consistent with the loss of liquid water upon drying. Second, the largest velocity | Li | is predicted to occur in the shortest arm. This is also expected, because the same pressure drop drives the liquid between the node and each meniscus. Third, when L 0 = 0, the equations are compatible with (4.14) when setting L = L 1 + L 2 . Fourth, the dynamics becomes nonlinear, contrary to (4.1) and (4.14).

We now solve the equations (4.18). As expected, the sum of the lengths obeys a simple equation:

L1 + L2 = - L 1 + L 2 + L 0 + 2L g τ , (4.19)
with solution L 1 + L 2 = (L 10 + L 20 + L 0 + 2L g )e -t/τ -(L 0 + 2L g ), denoting the initial conditions: L i (t = 0) = L i0 for i = 1, 2. Inserting this solution in (4.18), we find that the difference of the lengths obeys a linear equation, albeit with a time-dependent coefficient:

L1 -L2 = L 0 τ L 1 -L 2 (L 10 + L 20 + L 0 + 2L g )e -t/τ -(L 0 + 2L g ) .
This is a separable equation which is easy to solve, and which yields:

L 1 -L 2 = (L 10 -L 20 ) (L 10 + L 20 + L 0 + 2L g )e -t/τ -(L 0 + 2L g ) L 10 + L 20 L0/(L0+2Lg) × exp - L 0 L 0 + 2L g t τ .
Of course, this solution holds only as long as both lengths remain positive, which ceases when the length in the shortest arm decreases to zero. This occurs at a critical time which must be solved numerically using the expressions of L 1 and L 2 .

Experiments

Having derived all theoretical tools to predict the dynamics in our elementary networks, we now proceed to a quantitative description of the experiments, and compare them to our predictions.

Branched networks

Single-node networks

We begin by measuring the drying dynamics in single-node networks (Fig. 1). The distance travelled by the menisci is plotted as a function of time in Fig. 7; the time origin here is somewhat arbitrary, corresponding to the beginning of the camera recording. The distance increases with time, as drying proceeds. For a given number of branches, all branches show similar dynamics, with a drying velocity decreasing at increasing distance. The striking feature of the drying dynamics is a velocity jump: the meniscus velocity is larger in the main channel than in the branches, the jump being larger at increasing number of branches (Fig. 7).

The origin of the jump can be qualitatively explained from the discussion in Sec. 4.2, illustrated in Figs. 6a andb. The water length from the meniscus coming from the entrance channel, as it arrives at the node, equals N L b with L b the length of each branch. As it splits into N menisci, each of these has its own water length equal to the branch length L b . More quantitatively, the dynamics of all menisci should obey (4.1). To test this idea, we plot the meniscus velocity as a function of the water length for all menisci in Fig. 8a. It shows indeed that all data follow a common trend, which is well fitted by a straight line. The entrance channels show small deviations from this trend. We have no clear explanation for these deviations, but they may due to the somewhat uncontrolled opening of the channels by scalpel cutting (see Sec. 2.2). However, they remain within 10% from the straight line, suggesting that the simple model (4.1) already captures the dynamics of drying in branched channels with good accuracy. The value of the fitting parameters will be discussed in Sec. 6.1 for all experiments. As discussed at the end of Sec. 4.2, we should in principle not use the single-channel version of the flux Q g from the meniscus for the dynamics inside the branches, but rather the expression (A 7) derived in the Appendix, which accounts for the coupled evaporation from different menisci issued from the same node. Hence, Eq. (4.1) should be replaced by:

L + L + f L g τ = 0, (5.1)
with a factor f given by, after (A 7) and since L b = + L (Fig. 6b):

f = cosh[(L b -L)/L g ] + N sinh[(L b -L)/L g ] sinh[(L b -L)/L g ] + N cosh[(L b -L)/L g ] . (5.2)
To quantify the magnitude of this correction, we focus on the data for the branches, and we fit them by Eqs. (5.1) and (5.2), using τ and L g as free fitting parameters. The results are shown in Fig. 8b. They show indeed a small deviation from a linear relationship between the meniscus velocity and the water length as the latter tends towards the branch length, i.e. as the meniscus is close to the node. This is in qualitative agreement The line is a linear fit of all data, with best fitting parameters in (4.1): τ = 1.8 × 10 3 s and Lg = 1.1 mm. (b) Velocity of the meniscus as a function of the water length, in branches of the network of two ( ), three ( ), four ( ) and five ( ) branches. The curves are fits of the date for two (plain curve), three (dashed curve), four (dotted curve) and five (dash-dotted curve) branches by (5.1) and ( 5.2), with best fitting parameters: τ = 1.7 × 10 3 s and Lg = 1.1 mm. The straight line corresponds to Eq. (4.1) using the latter two parameters.

with our model of evaporation coupling of the menisci, even though the agreement is not fully quantitative, perhaps because the hypotheses underlying the model (notably that w and H are much lower than other lengths under consideration) become questionable in the range of length where evaporation coupling is significant. However, the main result from this analysis is that evaporation coupling has only a small effect on the velocity of the menisci inside the branches, because the difference of velocity is only 10% lower with coupling than without (compare the straight lines and the curves in Fig. 8b). Hence, to simplify the forthcoming discussion of our experiments, we shall simply use (4.1) as a prediction with L the total length of the water-filled branches ahead from the meniscus under consideration.

Multinode networks

Overall, the results on single-node networks suggest a strong correlation between the meniscus velocity and the water length. We further test this idea by studying more complex branched channels. We consider first the "tree" network of Fig. 2a. We first plot the distance travelled by the menisci as a function of time in Fig. 9a. This graph shows that at each node, there is a jump of velocity, consistently with the previous observations on single-node networks. Moreover, there is a signature of the asymmetry of the eight terminal branches: the drying occurs in shorter time, but with a lower velocity, in the two short branches (inset of Fig. 9a). According to the previous discussion, we measure the meniscus velocity and water length in the main channel and in each branch, and we plot the dependency between these two quantities in Fig. 9b. It shows that the data are gathered in four subgroups, depending on the level of the branch considered inside the tree; as in the previous case, the gaps between each subgroup is associated to a velocity jump as a meniscus splits in two at each node. Moreover, the data collapse very well on a straight line, which confirms that the simple model (4.1) describes well the drying dynamics even in cases with multiple nodes.

We finally consider the Adiantum design (Fig. 2b). Since there are many edges (45 in total), we measure only the average velocity in each edge, and we retain the standard deviation of velocity as an error bar. We plot the meniscus velocity as a function of the water length in Fig. 10. All data collapse again very well on a single straight line, with only one significant outlier. This collapse is all the more remarkable that the Adiantum design does not really obey the model approximations; because it comprises many veins, these are not all slender, and they are close enough that diffusive interactions may become significant.

Loops

Meniscus position and water length

In the branched networks without reconnections, we have quantified the position of the menisci by the distance travelled from the origin, which was unambiguous. Conversely, in loops, there are several possible choices to quantify the position x of a meniscus. We proceed as follows, and as illustrated in Fig. 11a; in each case, we also relate x to the water length L or to the quantities L 1 and L 2 previously introduced. For the meniscus in the terminal channel, we use the distance s from the end of the channel (Fig. 11a,i). In this case, we thus have simply: x = s and L = s.

When there is one meniscus per arm, for the meniscus in the short arm, we set x = L t + R θ, where we recall that L t the length of the terminal channel, R the radius of the loop, and where 0 θ β is the angular distance between the meniscus and the exit node (Fig. 11a,ii), and β the small angle between both nodes and s the distance between the meniscus and the entrance node (Fig. 11a,i). For the menisci in the long arm, we set x = L t + R θ with 0 θ 2π -β the angular distance between the meniscus and the exit node (Fig. 11a,ii). We thus have: L = L t + R (θ + θ ). For these menisci, the lengths L 1 and L 2 introduced in Sec. 4.3 are proportional to the angular distances: L 1 = R θ and L 2 = R θ. When there are two menisci in the long arm (Fig. 11a,iii), we define their position x ± = L t = R θ ± (with θ - θ + ), and the water length of the region between them is

L = R (θ + -θ -).
For the meniscus in the entrance channel, we set x = L t + R β + s, and we have:

L = L t + 2πR + s.

Symmetric loop

In the symmetric loop, the position of each meniscus is plotted as a function of time in Fig. 12a. The dynamics of the meniscus in the entrance channel and that in the terminal channel are smooth, and as before, the velocity of each meniscus decreases in time. The two menisci in the arms of the loop do not behave perfectly symmetrically, and their displacement curves show some irregular inflections, but the difference remains small and they reach the exit node simultaneously. Like for the single-node branched networks, there is a velocity jump as the meniscus splits in two at the entrance node, the velocity of each of the two new menisci being smaller. Interestingly, the opposite velocity jump is observed at the exit node: the velocity of the meniscus in the terminal channel is greater than that of the two menisci merging at the exit node.

We now compare this experimental dynamics with our model. The dynamics of the entrance and the terminal meniscus is simple, since each is the only meniscus present in the network during their period of existence; they should obey (4.3), provided the water length is defined as before. When there are two menisci in the symmetric loop, the model of Sec. 4.3 show that they should have the same velocity; it is a special case of the coupled equations (4.18). The fitting parameters τ and L g are taken once and for all the data (and not separately for each curve). Fig. 12a shows a good agreement between the experiments and the predictions, except for the aforementioned small difference between the two menisci in the two arms. Furthermore, the sum of their velocity, which is the drying velocity | L|, is predicted to follow (4.19), which is similar to (4.1) except for the prefactor of L g , which counts the number of menisci bounding the drying water region. Therefore, we plot the drying velocity | L| as a function of the water length in Fig. 12b. Indeed, the data gather in three subgroups, with the same slope as expected from (4.1) and (4.19), but with an offset when there are two menisci. More precisely, the data are very well fitted by (4.1) and (4.19).

Asymmetric loops

For the asymmetric loops, the three cases described for the symmetric loop, namely (i) the meniscus in the entrance channel, (ii) one meniscus in each arm of the loop and (iii) the meniscus in the terminal channel, are supplemented with a fourth case: the existence of a water region bounded by two menisci in the long arm. For all these cases, we plot the experimental dynamics of each single meniscus in Fig. 11b andc for two of the asymmetric loops, and we compare it to the predictions of Sec. 4.3. As discussed in that Section, the dynamics of the menisci in the entrance and terminal channel is predicted to obey (4.1). When there is one meniscus in the short arm and another in the long arm, their dynamics is predicted to obey the coupled equations (4.18), until the meniscus in the short arm reached the exit node. Finally, when there is a water region in the long arm bounded by two menisci, its dynamics is predicted to follow (4.14) where each meniscus contributes to half of the velocity. Fig. 11b andc shows that while the agreement between experiments and theory is good for the meniscus in the entrance and terminal channels, and still rather good for the two menisci bounding a water region in the long arm, there is a significant discrepancy in the phase of drying where there is one meniscus per arm. While our model predicts that the meniscus in the short arm should have the largest velocity, the contrary occurs in Fig. 11c. In Fig. 11b, even if the meniscus in the short arm is faster, the difference of velocities between the two menisci is much larger than the prediction. A likely origin of this discrepancy is discussed in Sec. 6.2.

To go further, we plot the water length of the different water regions as a function of time in Fig. 13a, c and e, and the drying velocity | L| as a function of the water length in Fig. 13b, d andf. It shows that in the spite of the discrepancy for individual menisci, the drying velocity of the water region bounded by a meniscus in the short arm and another in the long arm is well captured by the model. Another striking feature is the significantly larger drying velocity, for a given water length, of the water region in the long arm compared to the water region in the terminal channel. This is perhaps the best illustration of the influence of the contribution of the menisci to the drying of a given water region. For instance, in the asymmetric loop of Fig. 3b, even though it is initially twice shorter, the water region in the terminal channel takes almost the same time to dry as the water region in the long arm. Comparing Fig. 12a and Fig. 13a and c also shows that for the same total channel length, the total drying time is shorter in asymmetric loops compared to the symmetric one. To compare these data with the model, we proceed like for the symmetric loop, with as new feature the water region in the long arm. Since it is bounded by two menisci, its dynamics should obey (4.14). Hence as previously, we plot the velocity of the meniscus, or the sum of the velocity of the menisci if there are two, as a function of the water length, in Fig. 13b, d andf. Here again, the data gather in four subgroups, with the same slope but with an offset correlated to the number of menisci. The data are well fitted by (4.1), (4.14) and (4.19) with τ and L g as fitting parameters taken once and for all the data; the agreement is slightly less good than for the symmetric loop, but the deviations remain relatively small. We notice finally that there are no qualitative differences between the three asymmetric loops. In the most asymmetric loop (Fig. 3c), the duration of the phase of coexistence of two menisci in the two arms decreases, while the water region created in the long arm after meniscus splitting at the exit node is initially longer and takes a longer time to dry, as expected. For a longer terminal channel, the durations of the phase of coexistence and of the water region in the long arm are both reduced.

Discussion

In this section, we focus on two further aspects of our study. First, we must now discuss the values of the fitting parameters τ and L g , especially the first one which gives the typical drying time. Second, we shall discuss the discrepancy observed between the experiments and the model in the specific case of individual menisci in the loops.

Values of the fitting parameters

The fits of all our data are related to the two same parameters: the typical drying time τ , and a length L g related to the contribution of the meniscus to drying. We summarise the values of these parameters in Tab. 2, together with their theoretical predictions using (4.4) and (4.5).

The values of the typical time τ shows some dispersion: it is lowest for the singlenode networks, highest for the tree and Adiantum, and intermediate for the loops; the values for the tree and Adiantum on one hand, and for the various loops on the other hand, are very similar. This apparent dispersion is actually well correlated with the difference in thickness for the different designs (Tab. 1). Indeed, our model predicts the following dependence of τ on the geometrical parameters, see Eq. (4.5): τ ∝ hw/(ξ+w/δ). Therefore, we plot τ as a function of hw/(ξ+w/δ) in Fig. 14. The data then align correctly on an affine trend, although the data for the single-node networks lie slightly below the line. Moreover, the slope of the affine fit, equal to 1.9 s/µm 2 , compares favourably with the prediction, from (4.5): ρ/D P M C sat P = 1.4 s/µm 2 , although it is somewhat higher (likewise, fitted values of τ are larger than theoretical estimates). However, the values of C sat P and D P taken from the literature are known only with poor precision. Hence, to get a better check, we have repeated experiments on single channels, on two different chips, performed with the same protocol as for the various networks. Fits of the resulting curves L(t) (not shown) by (4.3) yield fitting values for τ which are added on Fig. 14 as cross symbols. These single-channel data follow the same trend as the data for the networks within a few percents. This constitutes a further validation of our models, and of their extension from single channels to networks.

Tab. 2 shows that the values of L g are more dispersed than that of τ , and the value for Adiantum is even negative. However, we believe that this dispersion is largely artificial. Indeed, we fit over quite an extended range of water lengths (larger than 50 mm for Adiantum and for the asymmetric loop with a long terminal channel, see Figs. 10 and13f) compared to the expected value of L g (of order 1 mm). Hence, a small uncertainty on the slope of the curves may result in a very large relative deviation on L g . More importantly, this shows that the meniscus contribution is secondary compared to the drying by direct pervaporation between the water-filled part of the channels and the outer atmosphere through the PDMS. 2) as a function of the geometrical parameter hw/(ξ + w/δ). The line is an affine fit of the data for networks (round symbols), with a value of the slope 1.9 s/µm 2 as free fitting parameter. Cross symbols are data for single channels.

Individual menisci in the loops

Most of our experimental measurements agree very well with our predictions. However significant deviations have been noticed for menisci inside the loops, more precisely when there is one meniscus per arm (Fig. 11b andc).

To discuss a likely origin of this discrepancy, we remind that we have based our model on a pressure argument. In the water regions, the pressure relative to the reference atmospheric pressure is fixed by a Laplace pressure jump across the menisci, ∆p c , and by a viscous pressure drop, ∆p v . We have assumed that the capillary jump is the same across each meniscus. At this stage, it is worth getting an order of magnitude of the two contributions. A reasonable estimate of the meniscus curvature is the channel width: κ ≈ 1/w = 10 4 m -1 , hence ∆p c = γκ ≈ 7 × 10 2 Pa with γ = 70 mN/m the surface tension of water. From (4.13), the order of magnitude of the viscous pressure drop is ∆p v ≈ ηq L 2 /hwS = ηL 2 /τ S, using (4.2). With w = 10 -4 m and h = 4 × 10 -5 m in our experiments, we compute S = 1.1 × 10 -10 m 2 from (4.11). With η = 10 -3 Pa•s the viscosity of water and using τ = 2.5 × 10 3 s (Tab. 2) and L ≈ 10 -2 m as the typical water length in our experiments, we get ∆p v ≈ 0.3 Pa. Even if this estimate is approximative, it shows that the viscous pressure drop is much smaller than the Laplace pressure jump, typically three orders of magnitude smaller. Hence, even a minute change in wetting conditions may easily induce a difference in Laplace pressure between different menisci of the order of the viscous pressure drop, and it can easily slow down or even immobilise one of the menisci while the theory of Sec. 4.3 would predict a smooth dynamics. We believe that this is the main source of uncertainty on the dynamics of individual menisci. It is probably illusory to control the wetting conditions to such an extent that variations of wetting conditions are such that differences of Laplace pressure jump remain below the viscous pressure jump, especially in view of the vegetal applications, and we have checked (data not shown) that redoing the PDMS design of the loops with the same recipe gives different individual dynamics. A crucial point is that, in spite of these uncontrolled variations of the dynamics of individual menisci, the sum of their velocities is perfectly smooth, since it is determined by the global water loss by pervaporation.

Conclusions

We have unravelled the dynamics of two elementary events during the drying of a hydraulic networks by pervaporation: (i) the splitting of a meniscus in several branches and (ii) the annihilation of two menisci.

Concerning the splitting, our main result is that further away from a node the meniscus velocity is always proportional to the remaining liquid length plus a constant axial diffusion flux. Just after splitting events, these fluxes interact before tending to a steady value. Because of the creation of new interfaces and thus a multiplication of axial diffusion fluxes, there is overall a slightly enhanced total drying rate compared to a drying process with only one meniscus.

Concerning the annihilation of menisci in loops around a liquid plug, we found that menisci velocities may present an irregular behaviour, owing to the extreme sensitivity of the system to slight inhomogeneities in wetting properties. However the sum of the velocities, thence the total drying rate, is remarkably regular.

This work opens a lot of perspectives. First, the present findings show that drying is enhanced by the multiplication of menisci and patterns with lot of branches should therefore accelerate drying compared to networks with the same water length but less branches. This principle should be helpful to design networks.

Second, an indetermination of velocities is often observed when a liquid region is bounded by several menisci before annihilation, even with homogeneously designed channels. We thus expect very original behaviour of menisci in networks with controlled geometrical inhomogeneities, such as varying channel width. This proves important since in real networks, such as the ones found in plant tissues, where the channel size is actually highly heterogeneous. with L g given by (4.4). In contact with the meniscus, air is saturated with water vapour, which yields the boundary condition:

c a (x = 0) = c sat a . (A 3)
The other boundary condition is:

c a (x = ) = c , ( A 4)
with c the concentration at the node, yet to be determined. Solving (A 2) with boundary conditions (A 3) and (A 4) yields: Let us now consider the entrance channel on the other side of the node, considered infinitely long. Let us denote x the distance from the node. The concentration field inside the entrance channel also obeys (A 2) with x instead of x , with boundary conditions: c a (x = 0) = c by continuity of the concentration field, and lim x →∞ c a (x ) = 0, hence: c a (x ) = c e -x /Lg .

(A 6)

To close the problem, in the limit where the volume of the node is negligible, the total diffusive flux from all channels issued from it must vanish, hence:

dc a dx x =0 = N dc a dx x = .
Computing these fluxes from (A 5) and (A 6), we obtain that c = N c sat a /[sinh( /L g ) + N cosh( /L g )]. From this value and (A 5), we can compute dc a /dx | x =0 and, from (A 1), the flux at the meniscus:

Q g = Q g0
cosh( /L g ) + N sinh( /L g ) sinh( /L g ) + N cosh( /L g ) , (A 7)

with Q g0 = hwD a V m c sat a /L g the flux at the meniscus in a single channel. The quantity Q g /Q g0 is plotted as a function of /L g , for various numbers of branches, in Fig. 15. At given N , it increases from 1/N for = 0 to 1 for /L g → ∞. At fixed /L g , it is a decreasing function of N . It is significantly smaller than one only if /L g remains close to one or lower.

Figure 1 .

 1 Figure 1. Snapshots of single-node networks, with (a) two, (b) three, (c) five and (d) four branches. (e) Transverse side view of the channels of width w and height h between a glass slide and a PDMS sheet of thickness H. (f) Longitudinal side view of a branch. The meniscus, sketched as a vertical line inside the channel, separates a water-filled region on its right side and an air-filled region on its left side. The varying blue shade depicts the decreasing water vapour concentration ca in the air-filled region, also sketched by the red dashed curve which shows the decrease of ca over the characteristic length Lg. The green dashed arrows show the water flux leaving the channel by pervaporation through PDMS. The red dashed arrow shows the evaporation across the meniscus. The blue arrows inside the water-filled region represent the liquid flow induced by pervaporation.

Figure 2 .

 2 Figure 2. Snapshots of branched networks with more than one hierarchical level: (a) network with three hierarchical levels, and (b) network from an Adiantum leaf (Brodribb et al. 2016a).

Figure 3 .

 3 Figure 3. Snapshots of the four designs with loops. In (a), (b) and (c), the terminal channel is inside the circle merely to spare some space. (e) Sketch of the different elements of the loop.

Figure 4 .

 4 Figure 4. Snapshots of the progressive drying in (a) a single-node network with four branches, and (b) the Adiantum design. Darker areas in the networks represent dried parts of the channels, while lighter areas are still filled with water. Horizontal bars represent a length of 1 mm.

Figure 5 .

 5 Figure 5. Snapshots of the progressive drying in (a) the symmetric loop, and (b) an asymmetric loop (the dynamics is qualitatively similar for the other asymmetric loops). Horizontal bars represent a length of 1 mm.

Figure 6 .

 6 Figure 6. In branched networks, (a) for a meniscus before a node, L is the total length of the water-filled branches ahead from the meniscus, here L = L0 + L1 + L2. (b) Meniscus created after splitting in a node, having travelled a distance from the node where it has been created. (c) Sketch of the topology corresponding to the phase of the drying dynamics in loops shown in Fig. 5a (snapshots at t = 15.6 and 39 min) and b (snapshot at t = 13.4 min).

Figure 7 .

 7 Figure 7. Drying dynamics in branched networks. Displacement of the meniscus from the origin as a function of time for networks with (a) two, (b) three, (c) four and (d) five branches. Symbols correspond to different branches as sketched on each subfigure: round symbols correspond to the main channel, other symbols correspond to branches. The origin is marked with a filled black circle.

Figure 8 .

 8 Figure 8. (a) Velocity of the meniscus as a function of the water length, in the main channel of the network of two (•), three ( * ), four (+) and five (×) branches, and in the branches (triangles).The line is a linear fit of all data, with best fitting parameters in (4.1): τ = 1.8 × 10 3 s and Lg = 1.1 mm. (b) Velocity of the meniscus as a function of the water length, in branches of the network of two ( ), three ( ), four ( ) and five ( ) branches. The curves are fits of the date for two (plain curve), three (dashed curve), four (dotted curve) and five (dash-dotted curve) branches by (5.1) and (5.2), with best fitting parameters: τ = 1.7 × 10 3 s and Lg = 1.1 mm. The straight line corresponds to Eq. (4.1) using the latter two parameters.

Figure 9 .Figure 10 .

 910 Figure 9. Drying dynamics in the network of Fig. 2a. (a) Displacement of the meniscus from the origin as a function of time (inset: zoom on the terminal branches). (b) Velocity of the meniscus as a function of the water length, and sketch of all branches with their symbols used in this figure. The line is a linear fit of all data, with best fitting parameters giving in (4.1): τ = 2.9 × 10 3 s and Lg = 0.8 mm.

Figure 11 .

 11 Figure 11. (a) Snapshots of the different menisci in a loop, superimposed with the symbols used to plot their positions in panels (b) and (c). (i) Meniscus coming from the entrance channel, and splitting at the entrance node into (ii) two menisci in each branch of the loop. Then the meniscus in the short arm splits at the exit node and (iii) the remaining water splits in two disconnected regions, one in the terminal channel and the other, bounded by two menisci, in the long arm. (b,c) Dynamics of individual menisci in two asymmetric loops: time evolution of the different menisci (b) in the loop of Fig. 3b and (c) in the loop of Fig. 3c; the menisci are distingued by different symbols and colours as shown by in panel (a). The lines are the predictions of the dynamics of each meniscus by the model of Sec. 4.3, with best fitting parameters: (b) τ = 2.6 × 10 3 s and Lg = 1.7 mm, and (c) τ = 2.5 × 10 3 s and Lg = 1.7 mm.

Figure 12 .

 12 Figure 12. Drying dynamics in the symmetric loop: (a) time evolution of the position of the menisci, in the entrance channel (•), in both arms (♦ and ) and in the terminal channel ( ). The curves are fits by (4.3) and (4.18), with best fitting parameters: τ = 2.6×10 3 s and Lg = 2.1 mm. (b) Drying velocity | L| as a function of the water length in the different phases of drying in the symmetric loop. The new symbol ( * ) as compared to panel (a) corresponds to the water length of the region bounded by two menisci, one in each arm. The lines are fits by (4.1) and (4.19).

Figure 13 .

 13 Figure 13. Drying dynamics in the asymmetric loops, of Fig. 3b (top), 3c (middle) and 3d (bottom). Left panels: time evolution of the water length, for the meniscus in the entrance channel (•), for two menisci with one in each arm of the loop ( * ), for the meniscus in the terminal channel ( ) and for the water region in the long arm ( ). Right panels: drying velocity | L| as a function of the water length in the different phases of drying. The lines are fits by (4.1), (4.14) and (4.19), with best fitting parameters: (b) τ = 2.6×10 3 s and Lg = 1.7 mm, (d) τ = 2.5×10 3 s and Lg = 1.7 mm and (f) τ = 2.6 × 10 3 s and Lg = 2.7 mm.

Figure 14 .

 14 Figure14. Value of the fitting parameter τ (see Tab. 2) as a function of the geometrical parameter hw/(ξ + w/δ). The line is an affine fit of the data for networks (round symbols), with a value of the slope 1.9 s/µm 2 as free fitting parameter. Cross symbols are data for single channels.

Figure 15 .

 15 Figure 15. Plot of Qg/Qg0 as a function of /Lg, for N = 2 (plain curve), 3 (long dashed curve), 4 (short dashed curve) and 5 (dash-dotted curve).

  c a (x ) = c sinh x /L g + c sat a sinh( -x )/L g sinh /L g . (A5)

Table 1 .

 1 Values of the channel thickness h and PDMS thickness H measured for all the channels.

	)

Table 2 .

 2 Values of the fitting parameters τ and Lg for all the channels, and their theoretical predictions.

		single-node networks tree Adiantum	loops
	Fig.	1	2a	2b	3a 3b 3c 3d
	τ (×10 3 s)	1.8	2.9	3.0	2.5 2.6 2.5 2.6
	theoretical τ (×10 3 s) from (4.5)	1.20		1.57	1.37
	Lg (mm)	1.1	0.8	-0.4	1.9 1.7 1.6 2.7
	theoretical Lg (mm) from (4.4)	0.72		0.82	0.77

† Email address for correspondence: benjamin.dollet@univ-grenoble-alpes.fr

We thank Mathieu Alonzo for preliminary experiments, and Danièle Centanni for experimental support. We acknowledge ANR for financial support, under the grant PHYSAP ANR-19-CE30-0010-02. The authors report no conflict of interest.

Appendix A. Computation of Q g in the case of multiple menisci

In this Appendix, we will quantify how the diffusive flux Q g from a meniscus is modified in the presence of multiple menisci. We consider the situation where N menisci are issued from the same node and are a distance from it; to simplify, we assume that this distance is the same for all menisci, a condition which is reasonably well met in experiments (Fig. 7).

Let us consider one of the N menisci. The diffusive flux Q g issued from it is given by:

where V m is the gas molar volume and c a (x ) is the water vapour concentration field in the air-filled part 0 < x < of the branch where the meniscus belongs, where x designs the distance from the meniscus. Eq. (A 1) relies on the assumption that the concentration field does not vary significantly in the cross-section, which as shown in [START_REF] Dollet | Drying of channels by evaporation through a permeable medium[END_REF] holds provided the amplification factor (4.6) is much larger than one, a condition met in practice.

Predicting Q g requires to predict the concentration field. In [START_REF] Dollet | Drying of channels by evaporation through a permeable medium[END_REF], we showed that it obeys the simple differential equation:


Sumanta Pattanaik

David Vanderhaeghe

Pascal Barla

Joëlle Thollot

François X Sillion

Dynamic point distribution for stroke-based rendering

Keywords: Preserved distrib, New distrib, Coherent distrib, Coherence params Importance map & distribution params

We present a new point distribution algorithm that is well adapted to stroke-based rendering systems. Its main characteristic is to deal efficiently with three conflicting constraints: the distribution of points should retain a good repartition in 2D; their motion should tightly follow the target motion in the underlying scene; and as few points as possible should be added or deleted from frame to frame. We show that previous methods fail to meet at least one of these constraints in the general case, as opposed to our approach that is independent of scene complexity and motion. As a result, our algorithm is able to take 3D scenes as well as videos as input and create non-uniform distributions with good temporal coherence and density properties. To illustrate it, we show applications in four different styles: stippling, pointillism, hatching and painterly.

Introduction

Creating animations by hand requires the artist to paint or draw each frame independently. In complex styles such as painterly, watercolor, airbrush, stippling or hatching, this process does not take into account frame-to-frame coherence of "brush strokes". This can lead to vibrating, evocative, and rich visual experiences when the technique is mastered (see the work of George Schwizgebel, Alexander Petrov, Frederic Back or Bill Plimton for instance). But this noise may also become annoying for the viewer, and artists may want to avoid it while keeping the richness of a complex style. However, ensuring frame-to-frame coherence is a very tedious and difficult task to achieve by hand.

One of the goals of non-photorealistic rendering is precisely to assist the user in complicated and tedious tasks, and automate most of the non-creative process. It then found a natural application in the creation of animations in various styles. Among these, stroke-based rendering [START_REF] Hertzmann | A survey of stroke-based rendering[END_REF] has received much attention, since it consists in assisting the artist in creating animations where images are composed of many small 2D drawing primitives; a task that is typically time-consuming when done by hand. Previous techniques rendered animations in specific styles, either from an input animated 3D scene or from a video. In this paper, as we consider both types of input, we will refer indiscriminately to a 3D scene or a video as the underlying scene.

By considering only small drawing primitives, a natural approach, taken by most of previous work, is to attach them to anchor points in the underlying scene. Hence, when put in motion, these anchor points will guide the primitives in the picture and give a sense of movement. While the definition of a primitive might vary depending on the chosen style, defining the anchor point distribution from frame-toframe is common and essential to any stroke-based system. It also raises multiple, conflicting constraints: the distribution of anchor points should retain a good repartition in 2D; their motion should tightly follow the target motion in the underlying scene; and as few points as possible should be added or deleted from frame to frame, a property often referred to as temporal coherence.

Naive distributions cannot satisfy all of the above requirements. For instance, if one distributes new drawing primitives at each frame, the first constraint might be easily satisfied, but since no motion is matched, and no coherence is ensured, the remaining two constraints are violated. Another approach would be to use a single fixed distribution for the whole animation (good repartition and of course perfect coherence), but it will obviously violate the motion constraint. More sophisticated approaches have been proposed in the literature, but as explained in Section 2, they fail to satisfy all the constraints in general cases. This motivates the development of a more flexible distribution method for the creation of stroke-based animations. In this paper, we present a new distribution technique that works for any kind of input, and achieves interactive frame-rates independently of the scene complexity. Regarding the afore-mentioned constraints, our algorithm, presented in Section 3, provides the following contributions:

• Generation of non-uniform Poisson-disk distributions with blue-noise properties, necessary for applications such as stippling or pointillism; • Matching of any motion (rigid-or soft-body deformations, vector flows, occlusions), making the approach tractable for arbitrary 3D scenes and videos; • Control of temporal coherence via the management of coarse-to-fine transitions when anchor points need to be added or deleted.

To achieve these goals, we sample stroke positions in the image plane and project them to the underlying scene to apply the scene motion. At each frame, previous samples points are reprojected to the image plane and updated with respect to an importance function. New points are inserted if needed. Our main results, presented in Section 4, show the advantages of our distribution over previous approaches to primitive distribution. However, for the purpose of illustration, we also show in Section 5 how our distribution algorithm can be used to render animations in various, albeit simple, styles.

Previous work

As our goal here is to compare with other distribution techniques, we first present classic methods employed in the sampling literature, before reviewing previous approaches encountered in stroke-based rendering.

Sampling techniques usually rely on an importance map to create non-uniform point distributions (with more points distributed in areas of higher importance). When applied to time-varying importance maps, either for video environment map sampling [HSK * 05] or for primitive distribution [START_REF] Secord | Fast primitive distribution for illustration[END_REF], the samples follow the motion of the peaks of the time-varying importance. This approach is well-suited for environment map sampling for instance, but in the case of primitive distribution, it results in points floating over the surface as the viewpoint or importance change.

In contrast, we are interested in having points follow a given motion, while their distribution matches a distinct importance map: for instance in stippling, points should follow object motion, with more points in dark regions of the image. The recent tile-based distribution method of Kopf et al. [START_REF] Kopf | Recursive wang tiles for real-time blue noise[END_REF] shows such a behavior: they propose a realtime blue-noise point distribution technique that matches very accurately any importance map. However, because their method relies on pre-computed tileable point distributions, it can only directly deal with 2D rigid motions (panning and zooming inside an image). Adapting the method to more general motions, for instance via non-uniform warping, is not trivial. In this paper, we thus take an alternative approach that does not rely on any precomputation.

Unlike sampling techniques, stroke-based rendering systems precisely focused on ways to convey motion. Different techniques have been proposed, depending on whether a 3D scene or a video is given as input. Because of the amount of strokes we seek to distribute, we only focus here on automatic distributions. In comparison, in the work of Kalnins et al. [KMM * 02] or Daniels [START_REF] Daniels | Deep Canvas in Disney's Tarzan[END_REF], the user has to manually place each stroke onto 3D surfaces and provide levels-ofdetail for each of them.

For 3D scenes, anchor points are usually directly distributed onto 3D object surfaces and projected in the picture plane to guide primitives motion. Introduced by Meier [START_REF] Meier | Painterly rendering for animation[END_REF] for painterly rendering using a static 3D distribution (hence limited in the range of possible viewpoints), the method has been extended by various authors to adapt to viewpoint changes [START_REF] Cornish | Viewdependent particles for interactive non-photorealistic rendering[END_REF][START_REF] Pastor | Real-time animated stippling[END_REF][START_REF] Nehab | Stratified point sampling of 3d models[END_REF][START_REF] Haller | Real-time painterly rendering for MR applications[END_REF]. They all first construct a hierarchy of 3D anchor points over an object's surface in pre-process. Then they choose a "cut" in the hierarchy at runtime so that the resulting projected points best match the required distribution (e.g. specified via an importance map). While this approach has the advantage of accurately following objects motion with a very good temporal coherence, it suffers from limitations that make it impractical for general scenes. First, the hierarchy construction, built in pre-process, might be time-consuming for large models and some knowledge of "how close" any object can come to the point of view is needed to control its precision. Second, as each object is considered individually, distribuc The Eurographics Association 2007.

Figure 2:

The different steps of our dynamic distribution algorithm: pruning, distribution, blending&sliding. A coherent nonuniform distribution is dynamically updated from frame to frame regardless of how the scene moved: 3D objects can be deformed, viewpoint motion is not restricted, and motion in video can be provided via its optical flow.

tions cannot overlap multiple objects (although it might be necessary, for instance with distant objects); they can only follow individual object motion with restricted deformations (no topology change, nor vector flows over surfaces for instance); and they become increasingly time-and memoryconsuming with scene complexity. Finally, the repartition constraint cannot be ensured since the points are only guaranteed to appear roughly in the middle of neighbor points when projected to the picture plane [START_REF] Pastor | Frame-Coherent 3D Stippling for Non-photorealistic Computer Graphics[END_REF]. Moreover, distributions obtained this way sometimes present artifacts such as patterns, accumulation or grouping motifs (especially at grazing angles). As a result, 3D-based techniques, while very efficient at ensuring temporal coherence, fail to fully satisfy motion and distribution constraints.

For videos, on the other hand, since there are no objects to rely on, no pre-process is needed. Rather, as in the compelling painterly renderer of Hays and Essa [START_REF] Hays | ESSA Image and video based painterly animation[END_REF], anchor points follow video optical flow. Provided the optical flow is a good approximation of motion, the resulting anchor point distributions satisfy the motion and temporal coherence constraints. However, this approach needs to be adapted to work with 3D scene input, as it does not take into account occlusions: for instance, if we use Hays and Essa's system on a motion flow extracted from a 3D scene, anchor points may move from one object to the other. Moreover, their method for dealing with distributions assumes a thick, opaque paint medium: multiple uniform distributions are painted in sequence, from coarse to fine. Not only this approach makes the system very slow (around 80 seconds per frame), but it will fail to adapt to other styles such as stippling or pointillism, since there are no continuous variations of anchor points density, and distributions layers are built independently so that two anchor points from different layers may be arbitrary close to each other.

The "animosaics" method of Smith et al. [START_REF] Smith | Animosaics[END_REF] provide a distribution approach similar to [START_REF] Hays | ESSA Image and video based painterly animation[END_REF], even if used another context: they take a vector based animation as input and mainly focus on packing mosaics with a uniform distribution. Therefore, regarding our goals, the approach shares the same drawbacks as Hays and Essa's work. Videobased methods hence have important limitations concerning motion and distribution constraints.

To overcome the limitations of previous work we make use of an hybrid approach that combines advantages of 3D-and video-based techniques, while having distribution properties sufficiently close to the ones obtained with sampling techniques, at least in the case of stroke-based rendering.

Distribution algorithm

Our goal is to distribute drawing primitives in screen-space in such a way that they follow a given motion in the depicted scene without exhibiting temporal artifacts. To permit non-uniform densities, we also allow the user to specify an importance map I, a function that assigns to any point p in screen-space an importance value I(p) ∈ [0, 1]. This function, which typically evolves during the animation, can either be defined by 3D scene properties (object IDs, depth, painted textures, etc) or screen-space properties (luminance, optical flow magnitude, etc).

We choose to ground our dynamic distribution on the static Poisson-disk distribution approach of McCool and Fiume [START_REF] Mccool | Hierarchical poisson disk sampling distributions[END_REF]. As argued in Section 3.1, it is best adapted to the particular constraints of stroke-based rendering. We then explain in Section 3.2 how we extend this algorithm to deal with dynamic settings. The underlying idea, summarized in Figure 2, is that at each step of our algorithm, points are moved to the next frame according to available scene motion; they are then pruned using a Poisson-disk criterion, and holes in the distribution are filled. The main contribution of this method is that through this process, the updated distribution retains the properties of the static algorithm.

Static distribution

We first present the guidelines of the Poisson-disk distribution technique of McCool and Fiume [START_REF] Mccool | Hierarchical poisson disk sampling distributions[END_REF], also called relaxation dart throwing. Standard dart throwing [START_REF] Cook R | Stochastic sampling in computer graphics[END_REF] randomly distributes points in the image plane and keeps a point only if no other point is found within a disk of radius r around it. McCool and Fiume extended this algorithm with a hierarchical formulation: initially, points are randomly distributed with a large radius r = rmax; once no more space has been found after n trials the radius is multiplied by α ∈ (0, 1). The algorithm stops as soon as r ≤ r min . The quality of the distribution can be controlled by the parameters n and α, allow us to run our algorithm at interactive frame rates for previewing (n = 100, α = 0.9), or create high quality animations for off-line rendering (n = 1000, α = 0.99). Apart from having good statistical properties (see the survey on Poisson-disk distributions by Lagae and Dutre [START_REF] Lagae | A Comparison of Methods for Generating Poisson Disk Distributions[END_REF]) and being simple to implement, the algorithm is well suited to our problem: it enables non-uniform densities, provides a hierarchical distribution and can be built incrementally. In particular, the point density can be made to follow an arbitrary importance map I defined in screenspace, via the use of an importance criterion: p is kept only if 1 -r(p)-rmin rmax-rmin > I(p), with r(p) ∈ [r min , rmax] the radius used to distribute it. This way, the user is able to finely control the distribution of points by manipulating r min and rmax.

Dynamic distribution

Up until now, we explained how points are distributed for a single frame. We start any animation with such a point distribution (Figure 3(a)). Then, we move points to the next frame. This produces clusters of points that we clean out by pruning unwanted points; but it also creates holes that we fill by distributing new points.

Depending on the input given to our system, motion is performed in different ways. To deal with complex 3D scene motions, we project distributed points from screen-space onto object surfaces, move objects (with any kind of motion, including skinning and deformations), and project points back to screen-space. In practice, we take as input triangle meshes, and the projection onto a surface is done by rendering at each frame buffers that hold mesh IDs, face IDs, and triangle barycentric coordinates; then for each point we obtain its corresponding IDs and coordinates simply by reading the buffers. In the next frame, the point is located in 3D and perspectively projected. This approach allows points to stick to their triangle even if they are severely distorted. For videos we simply translate points along an optical flow provided in input (we rely on a classical gradient based method available in Adobe R After Effects R ).

Another important aspect of motion is occlusion: while for videos, points are considered hidden only when they fall out of the image dimensions, they can disappear behind other surfaces when dealing with 3D scenes. We perform a simple depth test to detect occlusions, and tag hidden points that will be treated specifically in the next step. Moreover, points falling onto the background in 3D scenes are considered hidden. Figure 3(b) points after they have been moved to a subsequent frame, with hidden points in yellow.

In both videos and 3D scenes, points can be either visible or hidden after they have been moved to the next frame. In the following, we only consider points that are still visible. In order to clean out clutters, we re-insert them in the order they have been previously distributed using a Poissondisk criterion to reject unwanted samples. For each radius r L corresponding to a level L of the hierarchical distribution, we check all the previously inserted points against the importance criterion of Section 3.1. Accepted points are reinserted with radius r L and we try to re-insert the remaining points at level L -1. The process stops as soon as all the points have been re-inserted or the finest level has been checked, usually rejecting a subset of the previous points. These rejected points are tagged as deleted (red points in Figure 3), while accepted points (in blue) are preserved so that the distribution is coherent from frame to frame.

Finally, holes still remain in the distribution due to appearing surfaces or optical flow dilation. We thus need to distribute new points to fill these holes with the same properties as in the static case of Section 3.1. Considering all the re-inserted points as part of the current distribution, we reuse McCool and Fiume's algorithm: each level L of the hierarchical distribution is processed in a top-down approach; points are randomly distributed with a radius r L , until no more space has been found for n trials, discarding points that not meet the importance criterion. The algorithm stops after the finest level, corresponding to r ≤ r min , has been processed. This approach allows inserted points to appear in the highest possible level of the hierarchy, hence resulting in a better repartition.

Temporal coherence

Using the algorithm described in Section 3.2 guarantees that most of the previous points are kept from frame-to-frame. However, we cannot avoid inserting or deleting some points, as motion and density are usually conflicting constraints. Hence, in order to further improve the temporal coherence of our approach, we propose two mechanisms: blending and sliding.

Blending: Inserting or deleting points as soon as this is requested by our distribution algorithm often results in disturbing popping artifacts. As in previous work [PFS03, HE04, KCODL06], we rather blend points in or out using smooth transitions. A blending duration is assigned to each inserted and deleted point, that is incremented or decremented at subsequent frames. Therefore, the distribution algorithm only tags points for being inserted or deleted, and leaves the task of performing the operation to the blending mechanism.

The main difference in our approach compared to previous work is that we only focus on blending duration, so that the way blending is actually performed can use of this value, but vary depending on the application (see Section 5). Moreover, we propose to distinguish between two different cases where blending occurs: density and occlusion events. For occlusion and dis-occlusion events, points are deleted (resp. inserted) directly. In practice, detecting which points become occluded in the current frame is easy (using the zbuffer), but finding out what points are dis-occluded requires to store some more information (z-buffer and transformation matrices from the previous frame). For points that appear or disappear by density, we developed another mechanism that assigns blending durations proportional to the radius of the inserted or deleted point. In practice, it produces smooth transitions where "details" are added progressively, from coarse to fine.

Sliding: Another problem may appear when we remove a point due to the Poisson-disk constraint, and it leaves enough place to add another point in its neighborhood: even with blending enabled, it can result in a flickering artifact during animation. We avoid this problem by matching a pair {p d , p i } of deleted and inserted points, where p d is the closest deleted points to p i in screen-space, using the radius of p i as a threshold. Then instead of removing p d , we allow it to slide linearly towards the position of the newly inserted one. We only perform sliding for density events though, otherwise hidden points would tend to slide across occlusion boundaries.

In practice, we track the original points p d and p i and linearly interpolate between their projected positions. This results in a projected distance dp(t) between the sliding point and its target position, with t ∈ [0,ts] where ts is a sliding duration set by the user. While this would work for simple motions, more complex movements such as rotations can stretch Relative radius statistics α/n RDT ours 0.9/100 0.69 0.67 0.99/1000 0.80 0.74 out the projected pair, hence increasing dp(t); which is opposite to the goal of the sliding mechanism. Therefore, we also clamp dp(t) to the distance that would be obtained with a straightforward 2D linear interpolation, dp(0)(tts)/ts. This way, the interpolated position gets progressively closer to the inserted point, while still conveying object motion.

Results

Figure 1 shows stills from dynamic distributions created with our algorithm, either from 3D scenes or video input. The full animations can be viewed in the accompanying video (available at http://artis.imag.fr/ Publications/2007/VBTS07a/). The point distribution tightly matches required motion and density in all cases, while exhibiting good temporal coherence.

The sampling times are independent of scene complexity, and our method requires no pre-process. The performance of our distribution algorithm thus only depends on image dimensions and the choice of the n and α parameters that drive points density. Note that as in [START_REF] Mccool | Hierarchical poisson disk sampling distributions[END_REF], we use a uniform grid to speed up distance tests during distribution. For interactive manipulation (n = 100/α = 0.9) we get frame rates between 5 and 10 fps, while for off-line rendering (n = 1000/α = 0.99), we are able to render point distributions at approximately 1 fps. These parameter settings are the ones advocated by McCool and Fiume [START_REF] Mccool | Hierarchical poisson disk sampling distributions[END_REF] and in both cases, we render images at 720 × 540 resolution, and get performances similar to them. now present statistical and visual evaluations of our distribution algorithm in comparison with previous work.

Statistical evaluation

We first want to check that our dynamic distribution has qualities similar to that of McCool and Fiume's original algorithm [START_REF] Mccool | Hierarchical poisson disk sampling distributions[END_REF]. A study by Lagae and Dutré [START_REF] Lagae | A Comparison of Methods for Generating Poisson Disk Distributions[END_REF] proposed to use two measurements to evaluate Poisson-disk distributions: relative radius that measures the amount of "packing" in the distribution; and power spectrum that gives infor- mation about its frequency content. Regarding these measurements, that are only valid for uniform distributions, Mc-Cool and Fiume's static distribution method is one of the best available techniques. We found no similar evaluation tool for non-uniform distributions.

We thus measured the power spectrum and relative radii for the frames of a test animation using the torus model of Figure 3, and a dynamic uniform distribution that fills the image (using r min = 2, rmax = 10, α = 0.99 and n = 1000). We expected that our approach, by introducing coherence of points distribution, would result in worse relative radii, or reveal regularities in the power spectrum. However, it did not show significant influence (see Figure 4): the spectrum clearly exhibits blue-noise properties (we removed the central peak for clarity as in [START_REF] Lagae | A Comparison of Methods for Generating Poisson Disk Distributions[END_REF]) and relative results are only slightly better in the static case. Note that as our distribution is not periodic, horizontal and vertical lines appear in the spectrum. While these measures do not outperform static Poissondisk distribution methods, they show that our technique has comparable statistical qualities while it satisfies additional constraints, i.e. scene motion and temporal coherence.

Visual evaluation

previous section compared our approach to previous techniques using uniform distribution statistics. However, such comparison is sometimes infeasible, either because no statistics are available, or because we deal with non-uniform distributions. In this section, we thus resort to visual evaluation, to compare our approach with object-based methods and illustrate its behavior with non-uniform densities.

The pioneering work of Meier [START_REF] Meier | Painterly rendering for animation[END_REF] produced remarkably stable object-based distributions: this is because her distribution is static, and done once in pre-process. A similar image-based approach, obtained using a static 2D distribution, would introduce the so-called "shower-door effect". Thus the object-based approach of Meier gives much better results, albeit for the limited range of viewpoints where the distribution is dense enough, or not too packed. But as soon as we need to insert or remove points to allow camera and object motions to be unconstrained, we lose stability and need blending. This is a problem common to object-and image-based approaches: in the first case, points are distributed on a surface and selected/deselected depending on a refinement function; while in our case, points are added/deleted using a Poisson-disk criterion. While both approaches have similar problems, object-based techniques have the additional burden of finding a good refinement function. And as shown in Figure 5(a), there is no simple way to get a good 2D distribution starting from a 3D distribution, even after introducing randomization as in Pastor et al. [START_REF] Pastor | Real-time animated stippling[END_REF]; while our method, shown in Figure 5(b), simply avoids this problem by directly distributing points in the picture plane.

Non-uniform distributions created with our approach depend on the values of r min , rmax and α. While r min and rmax control the density in regions where the importance map is equal to 1 and 0 respectively, the value of α the continuity of the non-uniform distribution. As shown in Figure 6, a low value for α has an effect similar to quantization, while a reasonable value results in smooth gradations. In practice, we found that using α > 0.9 gives very good results in most cases.

Applications

We now present some applications of our distribution algorithm that work with 3D scenes or video input (with optical flow). Each of them target a different visual style: stippling, pointillism, hatching or painterly. Note that we are concerned here with a very low-level notion of style: stipples are dots of varying size, pointillism introduces color, hatching needs to deal with orientation (e.g. for cross-hatching), and painterly refers to the use of small paint strokes of varying color and orientation. There is much more research to do to characterize how a style can be specified, but this falls out of the scope of this paper.

Our stroke-based rendering prototype system is similar in spirit to previous painting systems [Mei96, DOM * 01, HE04, PFS03] in that it separates low-level strokes drawing from high-level style control. At each frame, stroke parameters are Stroke density is a special parameter that deserves more attention since it is common to all styles. In our system, it is controlled by a time-varying importance map that is input to the distribution algorithm at each frame. As it plays an important role in anchor points distribution, we implemented a variety of importance mappings in the form of GPU shaders, so that new mappings can be added easily. Each style uses a different importance mapping, as shown in Figure 7 and in the accompanying video.

Stippling: For stippling, we use an approach similar to Pastor et al. [START_REF] Pastor | Real-time animated stippling[END_REF]: stipples are simply black dots of varying size; their density follows the luminance information coming from the scene (our importance map); and stipples are blended in or out by making them grow or shrink.

Pointillism: Pointillism is a variant of stippling where we make use of color information to color the stipples, and the point distribution is close to uniform. The points are added or removed by alpha blending. Moreover, as in [START_REF] Meier | Painterly rendering for animation[END_REF], we jitter each point's color in the Lab color space.

Hatching: Hatching introduces two new parameters: strokes length and orientation. In order to make the distribution visible in our results, we intentionally chose a simple style: stroke length and orientation are uniform across the image. To create cross-hatching, we simply create two independent distributions of strokes with different orientations. Moreover, as in previous work [HP00, HE04], we cut strokes at object boundaries, that we measure as depth discontinuities.

Painterly: Finally, painterly rendering combines all the characteristics of previous styles at once (color, length and orientation) and adds another parameter: stroke thickness. We again set the parameter values to uniform to better show the underlying distribution. Paint strokes are blended using the half-toning technique of Durand et al. [DOM * 01].

Discussion and Future Work

The work presented in this paper proposes an new approach to the dynamic distribution of anchor points for stroke-based rendering: it consists in a hybrid method that combines the advantages of distributing points in image space and the richness of object-space motions.

In the near future we plan to adapt our algorithm to work with more complex motions, such as specular reflexions motion, or natural phenomena (e.g. fluid, smoke, etc). Indeed, the only requirement in our case is to be able to project a point to the scene, and after it has been moved, to project it back to the picture plane.

c The Eurographics Association 2007.

Another area of research resides in using more complex strokes, such as long paint strokes, or mosaics. Already, when dealing with stippling, we cannot ensure an accurate tone reproduction, even if we get satisfying results in practice. While with hatching or painterly styles, one interesting question would be to adapt the distribution to primitives direction (e.g. using "Poisson-ellipses"). But this will not be enough as soon as long strokes will be involved. In this case, one possibility would be to use our approach as an underlying structure, and attach long strokes to multiple anchor points. And in the case of mosaics, another issue is raised: we need to meet boundary constraints between primitives. Our approach already offers good packing properties, but additional mechanisms are needed to make mosaic primitives "touch" each other without inter-penetrating.

Finally, while our distribution, by construction, makes use of the coherence of points from frame to frame, we additionally proposed two temporal coherence mechanisms: blending and sliding. Our simple assumption being that by reducing the number of insertions and deletions, and making the remaining ones less noticeable, we effectively enhance coherence from frame to frame. Interestingly, while temporal coherence is a recurrent problem in non-photorealistic rendering, there is no standard way to measure it to our knowledge, making a formal study an interesting avenue of future work.

Figure 1 :

 1 Figure 1: Our algorithm allows to distribute points on arbitrary scenes with non-uniform density and temporal coherence. (a) A rigid 3D model. (b) a 3D scene composed of different objects. (c) A deformable model.

c

  The Eurographics Association 2007.

Figure 3 :

 3 Figure 3: (a) Points are initially distributed in screen-space, here using luminance to drive density. They are then projected onto the 3D object. (b) After they have been moved to the next frame, points are reprojected to screen-space. Red points are rejected based on a Poisson-disk criterion; blue points from the previous frame are still valid; green points are added to fill in the blanks; and yellow points are hidden. (c) The distribution after it has been updated. The red, green and blue frames locate corresponding points that have been respectively removed, added or preserved.

Figure 4 :

 4 Figure 4: Left: Power spectrum of our point distribution method obtained by averaging periodograms for each frame of a test animation (central peak removed). Right: Relative statistics for static relaxation dart throwing (RDT) and our dynamic point distribution method.

Figure 5 :

 5 Figure 5: Comparison with 3D-based distributions. (a) The technique of Pastor et al. [PFS03] (figure taken from their paper) cannot ensure a good repartition (left half), even after randomization (right half). (b) Our approach, in comparison, produces distributions with blue-noise properties.

Figure 6 :

 6 Figure 6: Non-uniform density. (a) Using α = 0.5 results in quantification effect while (b) with α = 0.9 the importance gradation is faithfully reproduced.

Figure 7 :

 7 Figure 7: Our distribution can be used as a basis for various styles: (a) Stippling rendition of a deformable hand model, (b) a hatching style applied to a rigid shoe model, (c) pointillist rendering of a scene consisting in multiple objects and (d) a painterly rendering of a video that follows optical flow.computed for the whole image and stored in a style buffer. Stroke anchor points are then distributed using our dynamic algorithm. And finally, primitives are drawn by taking their parameters at their anchor point location in the style buffer. Note however that nothing is done to make stroke parameters evolve coherently, hence stroke orientation, length or thickness may change suddenly in the animations. Such temporal artifacts are specific to every new style and we plan to study temporal coherence for individual styles in future work.

c The Eurographics Association 2007.

Acknowledgments

We would like to thank Florent Moulin and Laurence Boissieux for the 3D models; and Lionel Baboud, Cyril Soler, Kartic Subr and Kaleigh Smith for their suggestions.