Gianluca Lavalle
email: gianluca.lavalle@limsi.fr

Nicolas Grenier

Sophie Mergui

Georg F Dietze
email: dietze@fast.u-psud.fr

Solitary waves on superconfined falling liquid films

Solitary traveling waves are prominent features covering the surface of a falling liquid film and are known to promote heat/mass transfer. We focus on the little studied case where they are subject to an extremely confined counter-current gas flow, and we identify two novel secondary instabilities. At high gas velocities, a catastrophic instability develops, leading to flooding through wave reversal and liquid arrest. At lower gas velocities, an oscillatory instability occurs, producing a high-frequency periodic modulation of the wave height. Conjunction of this self-sustained oscillatory state and vortices forming in the liquid is shown to enhance mixing. We also show that the gas flow can cause extreme local film thinning, leading to almost dry patches where the liquid thickness is very small.

Flows through miniaturized channels are usually laminar and require active or passive methods to enhance mixing [START_REF] Stone | Engineering flows in small devices: Microfluidics toward a lab-on-a-chip[END_REF]. In the case of stratified two-phase flows, this can be achieved through flow structures resulting from interfacial instabilities [START_REF] Hu | Viscous wave breaking and ligament formation in microfluidic systems[END_REF]. In horizontal pressure-driven configurations, the instability is known to require a strong enough relative motion between the two phases [START_REF] Yih | Instability due to viscosity stratification[END_REF][START_REF] Yiantsios | Linear stability of plane poiseuille flow of two superposed fluids[END_REF]. By contrast, we study the vertical configuration of a gravity-driven falling liquid film in contact with an extremely confined gas, where the inertia-driven Kapitza instability [START_REF] Kapitza | Wave flow of thin layer of viscous fluid (in Russian)[END_REF] generates solitary surface waves that greatly enhance heat/mass transfer [START_REF] Yoshimura | Enhancement of mass transfer into a falling laminar liquid film by twodimensional surface waves-some experimental observations and modeling[END_REF]. In this case, the gas velocity can be used as a free parameter to act on the waves.

Many works have studied the effect of a gas flow on the stability and dynamics of wavy falling liquid films [START_REF] Alekseenko | Primary instabilities of liquid film flow sheared by turbulent gas stream[END_REF][START_REF] Tseluiko | Nonlinear waves in counter-current gas-liquid film flow[END_REF][START_REF] Trifonov | Instabilities of a gas-liquid flow between two inclined plates analyzed using the navier-stokes equations[END_REF]. Such studies have usually concentrated on weakly/moderately confined configurations. In that case, only the tangential interfacial viscous stresses exerted by the gas play a significant role, and turbulent conditions are required to affect meaningfully the liquid film [START_REF] Trifonov | Counter-current gas-liquid wavy film flow between the vertical plates analyzed using the Navier-Stokes equations[END_REF]. By contrast, we focus on extreme confinement levels , where the gas flow and liquid film are of comparable thickness [START_REF] Setu | Superconfinement tailors fluid flow at microscales[END_REF]. Here, gas pressure variations induced by surface waves become relevant [START_REF] Dietze | Wavy liquid films in interaction with a confined laminar gas flow[END_REF], and the film is affected by the gas even under laminar conditions [START_REF] Lavalle | Suppression of the kapitza instability in confined falling liquid films[END_REF]. This configuration is bound to be more prone to flooding, i.e. occlusion of the channel by the liquid [START_REF] Vlachos | Visual observations of flooding in narrow rectangular channels[END_REF], wave reversal [START_REF] Tseluiko | Nonlinear waves in counter-current gas-liquid film flow[END_REF], or liquid arrest [START_REF] Trifonov | Counter-current gas-liquid wavy film flow between the vertical plates analyzed using the Navier-Stokes equations[END_REF]. Yet, it is unclear how this unwanted phenomenon is provoked under extreme confinement. Does it result from a lack of saturated traveling-wave solutions or due to a loss of their stability? Answering such question may produce valuable insights for the design of falling film microreactors [START_REF] Zhang | Hydrodynamics and mass transfer of gas-liquid flow in a falling film microreactor[END_REF][START_REF] Seebauer | Tube bundle falling film microreactor for performing gas liquid reactions[END_REF][START_REF][END_REF].

We have studied numerically solitary traveling waves under the effect of an extremely confined counter-current gas flow. For the conditions investigated, which represent real liquid properties and confinement levels typical of falling film microreactors, we find that flooding can be triggered at high gas velocities by a catastrophic secondary instability leading to wave reversal and liquid arrest. At lower gas velocities, we identify an oscillatory secondary instability, which causes a periodic (in time and space) modulation of the wave height. This self-sustained oscillatory state, in conjunction with vortical structures forming in the liquid, is found to enhance mixing within the solitary waves. It can be controlled through the gas velocity and thus represents an attractive approach to enhance inter-phase heat/mass transfer in small-scale devices. When the liquid volume and not the liquid flow rate is fixed, increasing the gas flow can cause pronounced film thinning, producing almost dry patches where the liquid thickness becomes extremely small.

We consider a two-dimensional liquid film falling down a vertical wall in interaction with a counter-current gas flow confined by a second wall at a distance H ⋆ =1 mm (the star superscript denotes dimensional quantities). The working liquid (index l) corresponds to a 83% by weight aqueous dimethylsulfoxide (DMSO) solution typically used in experiments [START_REF] Dietze | Experimental study of flow separation in laminar falling liquid films[END_REF] and has density ρ l =1098.3 kg/m 3 , dynamic viscosity µ l =3.13 mPas, and surface tension σ=0.0484 N/m. The gas (index g) is air. We focus on solitary waves, i.e. composite structures consisting of a large tear-shaped hump and several capillary precursor ripples [START_REF] Alekseenko | Wave Flow of Liquid Films[END_REF] (see graph 1b). Solitary waves are typically produced by nonlinear coarsening dynamics [START_REF] Chang | Coarsening dynamics of falling-film solitary waves[END_REF] that restructure the initial train of short waves emerging from linear wave selection, as the result of secondary subharmonic and/or sideband instabilities [START_REF] Liu | Onset of spatially chaotic waves on flowing films[END_REF]. They produce maximal heat/mass transfer intensification in unconfined films [START_REF] Rastaturin | Optimal regimes of heat-mass transfer in a falling film[END_REF].

We simulate this flow using the low-dimensional model of Dietze and Ruyer-Quil [START_REF] Dietze | Wavy liquid films in interaction with a confined laminar gas flow[END_REF], which is based on the weighted residual integral boundary layer (WRIBL) approach [START_REF] Ruyer-Quil | Improved modeling of flows down inclined planes[END_REF][START_REF] Kalliadasis | Falling Liquid Films[END_REF], accounts for the full coupling between the two phases, and has been substantially validated [START_REF] Dietze | Wavy liquid films in interaction with a confined laminar gas flow[END_REF]:

∂ x q l + ∂ t h = 0, ∂ x q g -∂ t h = 0 (1a) {S i ∂ t q i + F ij q i ∂ x q j + G ij q i q j ∂ x h} = (1 -Π ρ ) -Ka ∂ xxx h + (C jl -Π µ C jg ) q j + J j q j (∂ x h) 2 + K j ∂ x q j ∂ x h + L j q j ∂ xx h + M j ∂ xx q j , (1b)
where the indices i and j should be permuted over the phase-indicators l and g using Einstein summation. All quantities have been rendered dimensionless using the scales L=(µ l /ρ l ) 2/3 g -1/3 , U=(µ l g/ρ l ) 1/3 , and T =L/U. Further, q l and q g designate the liquid and gas flow rates per unit width, h is the film height, and the coefficients S i , F ij , G ij , C ij , J j , K j , L j , and M j are known functions of h. The streamwise and crosswise coordinates are denoted x and y, and t denotes time. The dimensionless groups are the Kapitza number Ka=σρ 1/3 g -1/3 µ -4/3 l , the ratios of density and dynamic viscosity Π ρ =ρ g /ρ l and Π µ =µ g /µ l , and the Reynolds numbers Re l =q l0 and Re g =q g0 Π ρ /Π µ , based on the nominal flow rates q l0 and q g0 of the flat-film primary flow. For a counter-current gas flow, q g0 and Re g are negative.

We solve numerically equation ( 1), either on a domain of wavelength Λ subject to periodic boundary conditions, or on an open domain with inlet/outlet conditions, where waves of frequency f are produced through inlet forcing. In both cases, the total flow rate q tot =q l +q g is fixed. We also track traveling-wave solutions of (1) by numerical continuation [START_REF] Allgower | Introduction to Numerical Continuation Methods[END_REF] using the software Auto07P [START_REF] Doedel | AUTO07p: Continuation and bifurcation software for ordinary differential equations[END_REF]. Such solutions are stationary in the reference frame ξ=x-c t moving with the wave speed c. We check their stability by using them as initial conditions for transient periodic computations. Stable solutions will remain unchanged in such computations, whereas unstable ones will evolve toward a different state. Details of the numerical procedures employed for solving the model equations ( 1) are given in [START_REF] Dietze | Wavy liquid films in interaction with a confined laminar gas flow[END_REF].

To study the effect of the oscillatory instability on heat/mass transport (figure 3), we have also performed a direct numerical simulation (DNS) based on the full Navier-Stokes equations, which additionally solves a convectiondiffusion equation for the transport of a passive scalar. The solver employed for the DNS is Gerris [START_REF] Popinet | An accurate adaptive solver for surface-tension-driven interfacial flows[END_REF] and details of the numerical setup are provided in [START_REF] Dietze | Effect of wall corrugations on scalar transfer to a wavy falling liquid film[END_REF].

Figure 1 illustrates how solitary traveling-wave solutions of fixed liquid Reynolds number Re l and wave frequency f =f 0 =0.027 evolve under the effect of an increasingly strong counter-current gas flow. This scenario is most practical from an experimental point of view [START_REF] Kofman | Characteristics of solitary waves on a falling liquid film sheared by a turbulent counter-current gas flow[END_REF]. Graph 1a displays the maximum h max and minimum h min wave height as a function of the gas Reynolds number Re g . The solid (h max ) and dashed (h min ) branches of the displayed curve join at the limit point S, where the Kapitza instability is linearly suppressed and traveling-wave solutions disappear into the flat-film primary flow. The location of S, Re g =-168, can be predicted by linear stability analysis [START_REF] Lavalle | Suppression of the kapitza instability in confined falling liquid films[END_REF].

Graphs 1b, 1c, and 1d represent streamlines in the reference frame of the wave for the solutions marked by arrows 1, 2, and 3 in graph 1a. As the gas velocity is increased, the surface wave is shaped into a flattened elongated liquid hump containing an increasing number of counter-rotating vortices. These result from the subdivision of the single vortex in graph 1b, which is typically observed in roll waves [START_REF] Alekseenko | Application of PIV to velocity measurements in a liquid film flowing down an inclined cylinder[END_REF].

We now check the stability of traveling-wave solutions by advancing them numerically in time through transient periodic computations, as represented in figure 2. In these computations, q tot and the liquid volume (due to periodicity) are fixed, and thus the wave-averaged liquid flow rate ql =Λ -1 Λ 0 q l dx may vary in time w.r.t. the traveling-wave solution used as initial condition. We find that traveling-wave solutions on the black curve segment in graph 1a (Re g >-13.5) are stable. This is shown in graph 2a, representing h max time traces from different transient computations. The black dashed line corresponds to a computation started from the traveling-wave solution at Re g =-1 (arrow 1 in graph 1a), and we see that h max does not change over time. On the blue curve segment in graph 1a, -78<Re g <-13.5, traveling-wave solutions are subject to an oscillatory secondary instability. This produces a high-frequency selfsustained modulation of the wave height, as evidenced by our transient computation in graph 2a (blue solid line), which was started from a traveling-wave solution at Re g =-16 (arrow 2 in graph 1a). We see that h max oscillates in time with a secondary frequency f s ≈ 3f 0 much greater than the nominal frequency f 0 . The amplitude of the oscillation first grows exponentially, proving that oscillations result from linear instability, and then saturates to a constant value.

On the red curve segment in graph 1a, -168<Re g <-78, traveling-wave solutions are subject to a secondary instability that leads to a catastrophic event. This is shown in graph 2b, representing a space-time diagram of the film thickness h(x, t)/H, for a transient computation started from the traveling-wave solution at Re g =-80 (arrow 3 in graph 1a). The profile superimposed on the left represents the initial condition h(x, t=0). In the early stages, tf 0 <1, this profile propagates downstream without changing significantly. Then, for 1≤ tf 0 ≤2.5, a reversal of the wave occurs under the effect of the counter-current gas flow. Finally, for tf 0 >2.5, a new upward traveling wave is formed (profile superimposed on the right). Importantly, the observed wave reversal coincides with liquid arrest, as evidenced by the corresponding time trace of ql , represented by a thick red line in graph 2c, which shows the wave-averaged liquid flow rate to become negative. In a real system, where the liquid inlet flow rate is fixed, wave reversal and liquid arrest are usually associated with flooding (supplementary movie [31] shows how these phenomena unfold in a spatially evolving computation). Thus, based on graph 1a and figure 2, we may conclude that flooding in the studied configuration is triggered by a secondary instability, whereas traveling-wave solutions never cease to exist, up to the linear instability cutoff (point S in graph 1a).

We now concentrate on the oscillatory instability observed for -78<Re g <-13.5. Figure 3 shows how this instability develops spatially, based on transient computations on an open domain with inlet/outlet conditions. Here, Re l is imposed through the time-averaged liquid inlet flow rate ql (x=0)=q l0 , which is subjected to a small harmonic temporal perturbation of frequency f 0 . This allows forcing the desired solitary waves, as is typically done in experiments [START_REF] Kofman | Characteristics of solitary waves on a falling liquid film sheared by a turbulent counter-current gas flow[END_REF].

Graphs 3a and 3b represent computations at Re g =-19 and Re g =-25. The black lines represent snapshots of the film profile (Eulerian view), whereas the red lines trace the path followed by the wave crests as time progresses (Lagrangian view). The oscillatory nature of these paths implies that the wave height of a given solitary wave traveling through the domain undergoes a self-sustained periodic oscillation (in time and space). In graph 3a, the oscillation consists of a single mode with a wavelength about twice smaller than the nominal separation of the solitary waves. In graph 3b, where the gas flow is stronger, the oscillation is enriched by additional modes making it more choppy, but it remains periodic. In both cases, the amplitude of the oscillation is substantial. The green profile in graph 3a corresponds to our DNS, performed on a shorter domain, and supplementary movie [32] shows that this DNS also displays oscillations. The present oscillatory instability differs from the secondary sub-harmonic and side-band instabilities observed for liquid films falling in a quiescent atmosphere [START_REF] Liu | Onset of spatially chaotic waves on flowing films[END_REF]. Those instabilities apply to trains of short periodic waves and result from wave interaction, leading to wave coalescence. The amplitude modulations they induce are of longer and not shorter length than the nominal wavelength. Importantly, these instabilities do not require any interaction between the gas and the wavy liquid film, whereas our instability is precisely caused by such an interaction.

The conjunction of the oscillatory instability (figure 3) and the vortices within the solitary waves (graph 1b) enhances convective heat/mass transport in the liquid film. To show this, we have solved a convection-diffusion equation for the representative scalar Θ in our DNS of the oscillatory regime in graph 3a. We have set Θ=1 at the gas/liquid interface, neglecting the transport resistance in the gas, and Θ=0 at the liquid inlet. The Péclet number, relating convective to diffusional transport, is Pe=UL/α=4590, where α denotes the diffusivity of Θ in the liquid.

The color plots in graphs 3c, 3d, and 3e represent the instantaneous Θ field at successive time points over the course of one nominal wave period 1/f 0 . We focus on the thin scalar plume that is seen to grow from the leading toward the trailing part of the wave in the successive snapshots. It results from the stagnation point flow generated by vortices within the wave (see e.g. graph 1b). This flow destroys the interfacial Θ boundary layer at the front of the wave but compresses the one at the wave back. The latter mechanism is responsible for the significant enhancement of inter-phase heat/mass transfer due to solitary waves [START_REF] Rastaturin | Optimal regimes of heat-mass transfer in a falling film[END_REF], but usually breaks down once the growing plume has reached the rear part of the wave [START_REF] Dietze | Effect of wall corrugations on scalar transfer to a wavy falling liquid film[END_REF]. In our case, the oscillatory instability prevents this from happening, by causing a repeated folding of the plume as it approaches the wave back (graphs 3d and 3e). We had already observed in graphs 1b to 1d the capacity of the gas to rearrange vortical structures within the wave. The supplementary movie [32] shows the plume folding dynamics in action. As a result of this mechanism, the interfacial boundary layer growing into the upstream residual film remains thin, and this increases the overall heat/mass transfer to the liquid.

We shall discuss now how solitary traveling waves behave under an increasing gas flow, when not Re l but the liquid volume is fixed. Such a situation may arise within the pulmonary capillaries [START_REF] Grotberg | Respiratory fluid mechanics[END_REF], where films of secreted mucus are subject to the respiratory air flow. We find that traveling waves can break down due to an extreme local thinning of the liquid film. Figure 4 shows how this occurs for a solitary wave of fixed wavelength Λ=46.8H and liquid volume (per unit width) V l =0.31ΛH. Graph 4a maps film thickness profiles h(x) as a function of the gas Reynolds number Re g (vertical axis), using grey scales. Graph 4b represents the corresponding evolution of the minimal film thickness h min . As the gas velocity is increased, the thin regions of the liquid film are increasingly drained, eventually forming patches that are almost dry. Our numerical continuation breaks down at Re g =-196.3, having reached h ⋆ min =62 nm, which is not far beyond the range of long-range van der Waals forces [START_REF] Israelachvili | Intermolecular and Surface Forces[END_REF]. It is interesting that purely hydrodynamic phenomena lead the film to this critical state, which may be conceived as negative flooding, i.e. the gas "touching" the wet wall instead of the liquid "touching" the dry wall. Of course, to check for actual film rupture, intermolecular interaction properties such as disjoining pressure and wall wettability need to be accounted for [START_REF] Bonn | Wetting and spreading[END_REF]. Integrating such a description [START_REF] Thiele | Sliding drops in the diffuse interface model coupled to hydrodynamics[END_REF] into our model, might be an interesting problem for future work.

In this letter, we have explored a new region of parameter space by considering superconfined wavy falling liquid films in contact with a counter-current gas flow. We have focused on a scenario in which the liquid flow rate and frequency of surface waves are fixed while the gas flow velocity is increased, which is most plausible from an experimental perspective. And, our computations were performed for an existing liquid [START_REF] Dietze | Experimental study of flow separation in laminar falling liquid films[END_REF]. We find that traveling-wave solutions exist no matter how strong the gas flow is, up until the point where the primary instability disappears. However, depending on the gas velocity, they are subject to two different secondary instabilities. At high velocities, the instability is catastrophic and produces flooding through wave reversal and liquid arrest. At lower velocities, the instability is oscillatory, producing a high-frequency periodic (in time and space) modulation of the wave height. This self-sustained oscillation is found to intensify convective transport within the liquid film and thus represents an attractive control mechanism to enhance inter-phase heat/mass transfer. We have also found a scenario in which the effect of a counter-current gas flow generates extremely thin patches within the liquid film.

We conclude by evaluating the feasibility of experiments. Zhang et al. [START_REF] Zhang | Hydrodynamics and mass transfer of gas-liquid flow in a falling film microreactor[END_REF] have produced experimentally stratified water/air flows in falling-film microreactors, where the liquid flowed in 0.3 mm deep parallel grooves of width W ⋆ =1 mm, with its surface pinned on either side to the top of the grooves [START_REF] Al-Rawashdeh | Pseudo 3-d simulation of a falling film microreactor based on realistic channel and film profiles[END_REF], and the gas confined by an upper wall at H ⋆ =0.8 mm. In our case, pinning should be avoided to allow for the motion of the liquid-gas interface on the side walls. Then, the contact angle φ comes into play, and this limits the channel width W in regard to forming a liquid annulus or rupturing the film laterally. Also, W should be sufficiently small to limit formation of 3-dimensional surface waves [START_REF] Kofman | Three-dimensional instabilities of quasi-solitary waves in a falling liquid film[END_REF][START_REF] Dietze | Three-dimensional flow structures in laminar falling liquid films[END_REF]. At the same time, it should be as large as possible, in order to limit the stabilizing effect of the side-wall meniscus on the interfacial instability [START_REF] Pollak | Side wall effects on the instability of thin gravity-driven films -from long-wave to short-wave instability[END_REF]. For water films on PMMA (Polymethylmethacrylate), steel, or PDMS (Polydimethylsiloxane), channel widths of W =5-10 mm should allow producing the desired conditions. For the DMSO-water solution used in our computations, the contact angle on PDMS is φ=93 • [START_REF] Zhang | Binary liquid mixture contact-angle measurements for preciseestimation of surface free energy[END_REF], and this reduces the ability of the liquid to climb the side walls and the stabilizing effect of the meniscus [START_REF] Pollak | Side wall effects on the instability of thin gravity-driven films -from long-wave to short-wave instability[END_REF], relaxing several constraints.

1 FIG. 2 .

 12 FIG. 2. Transient periodic computations started from traveling-wave solutions marked by arrows 1, 2, and 3 in graph 1a. The total flow rate qtot and the liquid volume are fixed. (a) Time evolution of maximal wave height hmax over several nominal wave periods 1/f0. Black dashed: stable regime (arrow 1); blue solid: oscillatory instability regime (arrow 2); (b) space-time plot of film height contours h(x, t)/H for the regime of catastrophic instability (arrow 3), displaying wave reversal. Profiles superimposed on left and right correspond to start and end times; (c) time-traces of the wave-averaged liquid flow rate ql =Λ -1 Λ 0 q l dx for computations from graphs a (black and blue lines) and b (thick red line), showing liquid arrest for the catastrophic instability.

FIG. 3 .

 3 FIG. 3. Regime of oscillatory instability. (a,b) Transient computations on an open domain with inlet/outlet conditions: H ⋆ =1 mm, Re l =15, f0=0.027. Displayed length is 0.8 m. Black lines: snapshots of the film profile y=h(x, t=const) (Eulerian view); red: path followed by the wave maxima, as the oscillating wave humps propagate through the domain (Lagrangian view); green: snapshot from DNS. (a) Reg=-19: single-mode oscillation; (b) Reg=-25: multi-mode oscillation. (c-e) Heat/mass transport data obtained from DNS in graph a: Pe=4590. See also supplementary movie [32]. Instantaneous color plots of the field of the transported scalar Θ over the course of one nominal wave period 1/f0. (c) tf0=0; (d) tf0=0.65; (e) tf0=1.0.

FIG. 4 .

 4 FIG. 4. Extreme film thinning under an increasing counter-current gas flow. Traveling-wave solutions at fixed Λ=46.8H and liquid volume (per unit width) V l =0.31ΛH. (a) Grey scale plot of film surface profiles h(x)/H as a function of Reg (vertical axis); (b) corresponding evolution of the minimal film thickness hmin (vertical axis is logarithmic). The lowest thickness reached is h ⋆ min =62 nm. In both graphs, two representative cross-sectional views are also included.

ACKNOWLEDGMENTS

We appreciate helpful discussions with J. P. Hulin. The work was supported by the ANR wavyFILM project, grant ANR-15-CE06-0016-01 of the French Agence Nationale de la Recherche.


Hayate Nakamura

Tetsuya Ogawa

Motochika Inoue

Takuma Hori

Lizhong Mu

Harunori N Yoshikawa
email: harunori.yoshikawa@univ-cotedazur.fr

Farzam Zoueshtiagh

Georg Dietze

Takahiro Tsukahara

Ichiro Ueno

Pumping effect of heterogeneous meniscus formed around spherical particle

Keywords: dynamic wetting, meniscus, capillary 2010 MSC: 00-01, 99-00

Hypothesis: A disturbance such as a microparticle on the pathway of a spreading droplet has shown the tremendous ability to accelerate locally the motion of the macroscopic contact line (Mu et al., J. Fluid Mech. 830, 2017) [1]. Although this ability has been linked to the particle-liquid interaction, the physical mechanisms behind it are still poorly understood despite its academic interest and the scope of numerous industrial applications in need of fast wetting.

Experiments: In order to better understand the mechanisms behind the particleliquid interaction, we numerically investigate the pressure and velocity fields in the liquid film. The results are compared to experiments assessing the temporal shape variation of the liquid-film meniscus from which pressure difference around

Introduction

Many natural systems exploit the wetting of their surfaces in various aspects. For example, lotus leaves repel water droplets by their microstructures [START_REF] Barthlott | Purity of the sacred lotus, or escape from contamination in biological surfaces[END_REF]; nepenthes rafflesianas profit of quick liquid transport thanks to their unique surface structures [START_REF] Chen | Continuous directional water transport on the peristome surface of nepenthes alata[END_REF]; water-walking insects take postures that deform the free surface and enable them to leap over menisci [START_REF] Hu | Meniscus-climbing insects[END_REF]. Wetting also finds its importance in numerous industrial fields. Efficient heat transfer in heat pipes requires enhanced wetting of their inner surfaces by working fluids. The functioning of a micro-total analysis system (µ-TAS) may fully depend on the wetting phenomena in microfluidic channels and involved elements [START_REF] Jouhara | Heat pipe based systems -advances and applications[END_REF]. For instance, the transport of a small amount of liquid into a small inspection region of µ-TAS is realized by chemical [START_REF] Kataoka | Patterning liquid flow on the microscopic scale[END_REF][START_REF] Kessenich | Connecting wettability, topography, and chem-istry in a simple lipid-montmorillonite system[END_REF] or physical functionalization of surfaces [START_REF] Ang | Graphene-mediated microfluidic transport and nebulization via high frequency rayleigh wave substrate excitation[END_REF] where controlled wetting is of utmost importance.

Owing to its importance, numerous studies on wetting have been reported.

Liquid spreading on a solid surface can be characterized by the behavior of a contact line (CL), which is a junction of solid, liquid, and gas phases. The dynamics of a CL can be characterized by the spreading parameter, S = γ SG -γ SL -γ, where γ SG , γ SL and, γ are the interfacial tensions of solid-gas, solidliquid and liquid-gas interfaces, respectively [START_REF] Ross | The history of the spreading coefficient[END_REF]. A contact line can be pinned onto a solid surface under S < 0 to form a static contact angle predicted by Young's equation [START_REF] Young | An essay on the cohesion of fluids[END_REF]. A CL can advance under S > 0 to extend the wetted area of surface. Focusing on this phenomenon, Hoffman performed experiments of pushing silicone oil in a glass tube with a piston [START_REF] Hoffman | A study of the advancing interface. i. interface shape in liquid-gas systems[END_REF]. The correlation between the contact angle θ and the capillary number Ca = µU/γ was derived; here, µ and U are the dynamics viscosity of the liquid and the velocity of the CL, respectively. The capillary number is an important parameter for estimating the effects of viscous stress on liquid spreading.

The correlation between the spreading radius R and the spreading time t of liquid film on a smooth horizontal substrate can be expressed as R ∼ t a [START_REF] Lopez | Spreading kinetics of liquid drops on solids[END_REF][START_REF] Tanner | The spreading of silicone oil drops on horizontal surfaces[END_REF][START_REF] Cazabat | Dynamics of wetting: effects of surface roughness[END_REF]. Here, a is an exponent that depends on the dominant force acting on the liquid film; a = 1/10 and 1/8 in the regimes of the interfacial-force dominant [START_REF] Tanner | The spreading of silicone oil drops on horizontal surfaces[END_REF] and the gravity-dominant [START_REF] Lopez | Spreading kinetics of liquid drops on solids[END_REF], respectively. In particular, the parameter becomes a = 1/7 in the one-dimensional spreading when surface tension is dominant. In the case of rough surfaces, it was indicated that the roughness enhances the wettability under S > 0 [START_REF] Wenzel | Resistance of solid surfaces to wetting by water[END_REF]. Microstructures of solid surfaces change the wetting behavior significantly [START_REF] Courbin | Imbibition by polygonal spreading on microdecorated surfaces[END_REF][START_REF] Xiao | Microscale liquid dynamics and the effect on macroscale propagation in pillar arrays[END_REF][START_REF] Obara | Imbibition of a textured surface decorated by short pillars with rounded edges[END_REF][START_REF] Wang | Surface structure determines dynamic wetting[END_REF][START_REF] Kim | Dynamics of hemiwicking[END_REF]. The coupled effects of surface chemical properties and surface topology on the exponent were investigated [START_REF] Kessenich | Connecting wettability, topography, and chem-istry in a simple lipid-montmorillonite system[END_REF].

Recently, the dynamics of an advancing CL interacting with spherical particles placed on a horizontal substrate was investigated [START_REF] Mu | Sharp acceleration of a macroscopic contact line induced by a particle[END_REF][START_REF] Kondo | Acceleration of the macroscopic contact line of a droplet spreading on a substrate after interaction with a particle[END_REF]. It was shown that the CL can undergo a large acceleration up to 10 times faster than in the absence of particles. Mu et al. [START_REF] Mu | Sharp acceleration of a macroscopic contact line induced by a particle[END_REF] attributed this behavior to the curvature of the meniscus formed around spherical particles and showed that the observed large velocities can be explaneed by their theoretical model based on the Cox's law [START_REF] Cox | The dynamics of the spreading of liquids on a solid surface. part 1. viscous flow[END_REF]. Mu et al. [START_REF] Mu | Control of local wetting by microscopic particles[END_REF] extended the experimental studies of Mu et al. [START_REF] Mu | Sharp acceleration of a macroscopic contact line induced by a particle[END_REF] and Kondo et al. [START_REF] Kondo | Acceleration of the macroscopic contact line of a droplet spreading on a substrate after interaction with a particle[END_REF] but with different particles in cylindrical and triangular prism shapes. They found that spherical particles caused the largest acceleration.

They also showed that the sequential arrangement of microstructures enables a control over the direction of liquid transport, which would constitute a novel efficient alternative to existing liquid transport techniques [START_REF] Mu | Quick liquid propagation on a linear array of micropillars[END_REF].

In this paper, we perform numerical simulations of the interaction between a CL and a single particle by a numerical code based on the InterFoam solver of the open-source software OpenFOAM. The goal is to elucidate the details of the local flow characteristics in the wetting process and unravel the physical mechanisms behind the acceleration of the CL. The behavior of CL predicted by the simulations are first compared with experiments to validate our numerical code. The velocity field and pressure distribution obtained by the simulations are then examined. Based on the computational results, the physical mechanism responsible for the CL acceleration is unraveled.

Methods

Numerical simulations

The open-source solver for computational fluid dynamics, OpenFoam (ver. 5.0), is used to carry out the present numerical simulations. The free-surface modeling technique of the volume of fluid [START_REF] Hirt | Volume of fluid (vof) method for the dynamics of free boundaries[END_REF] is used to track the gas-liquid two-phase flow. The governing equations for the simulations are:

∇ • U = 0 (1) ∂ρU ∂t + ∇ • (ρU U ) = -∇p + ∇ • τ + ρg + F SV (2) ∂α ∂t + ∇ • (αU ) + ∇ • (U r α (1 -α)) = 0. (3)
Equations ( 1 values of 1 and 0 when an computational cell is fully occupied by liquid and gas, respectively, and takes in-between values when an interface is present in the cell. Solving the advection equation enables to track the interface. Note that in order to counteract the numerical diffusion, a compression velocity is used.

We define α = 0.5 isosurface as the liquid-gas interface. The stress due to the interfacial tension is described by F SV = γκn = γκ∇α/∥∇α∥ by continuum surface force method [START_REF] Brackbill | A continuum method for modeling surface tension[END_REF], where κ and n are the curvature of the interface and the unit normal vector of the interface pointing to the gas phase, respectively.

The computational domain is shown in The advancing/receding contact angles for the substrate and particles are set to be 0 • and 10 • , and 10 • and 30 • , respectively. We set the physical properties of liquid and gas as those of 2-cSt silicone oil and air used in the experiments; the surface tension, density and kinematic viscosity of the liquid are γ = 1.83×10 -2 N/m, ρ l = 873 kg/m 3 , ν l = 2.0 × 10 -6 m 2 /s, and the density and kinematic viscosity of the gas are ρ g = 1 kg/m 3 and ν g = 1.48 × 10 -5 m 2 /s, respectively.

The gravitational acceleration is set to be constant at g = 9.81 m/s 2 . At time t = 0 s, the liquid film starts to flow into the computational domain under the condition of P inlet = 0 from the inlet of height h inlet = 20 µm at the x = 0 boundary. The capillary number of the liquid flow is small, e.g., Ca = 3.01 × 10 -4 just before the particle start to be wetted, but larger than in the experiments (Sec. 2.2 and Ref. [START_REF] Mu | Sharp acceleration of a macroscopic contact line induced by a particle[END_REF]) in order to reduce the time for flow computation. Kondo et al. [START_REF] Kondo | Acceleration of the macroscopic contact line of a droplet spreading on a substrate after interaction with a particle[END_REF] indicated that the ratio of the maximum velocity of the contact line to the velocity just prior to the interaction with the particle on the substrate become drastically smaller as increasing Ca. It was indicated that the maximum velocity decreases accompanying with the particle sucked to the bulk of the droplet after the contact with the contact line of the droplet. It must be noted that such trend was obtained in the system with the mobile particle just settled on the substrate. In the present simulation, we employ the particle fixed to the substrate, thus the trend observed by Kondo et al. [START_REF] Kondo | Acceleration of the macroscopic contact line of a droplet spreading on a substrate after interaction with a particle[END_REF] never takes place. Corresponding Reynolds number Re = U 0 D p /ν and Ohnesorge number Oh = √ ρν 2 /(γD p ) = √ Ca/Re for the present simulation are Re = 7.89×10 -2 and Oh = 6.18×10 -2 , respectively. Here U 0 is the velocity of the macroscopic contact line just prior to the contact with the particle. We limit ourselves to deal with the phenomenon under which the effect of inertia is negligible.

Experiment

Experiments are conducted to assay the similarity and the validity of the computational results. The setup, identical to that of Mu et al. [START_REF] Mu | Sharp acceleration of a macroscopic contact line induced by a particle[END_REF], is illustrated in Fig. 2. For detailed description of the setup refer to Mu et al. [START_REF] Mu | Sharp acceleration of a macroscopic contact line induced by a particle[END_REF] The test fluid Co., Ltd.) of diameter D p = 200 µm are used. The particle size here is larger than those used in the previous experiment [START_REF] Mu | Sharp acceleration of a macroscopic contact line induced by a particle[END_REF] and in our numerical simulations in order to capture the gas-liquid interface with higher accuracy. The horizontal distance between the centers of the particle and the needle tip is 4.0 mm. The interaction between the liquid and the particle is monitored simultaneously from top and side by high-speed cameras (Fastcam Mini WX100, Photron Ltd.). The resolution in top and side views are 2.17 µm/pix and 0.77 µm/pix, respectively.

The frame rates of both cameras are 50 frames per second (fps). A 100 W cold light (MegaLight 100, Moritex Corp.) is used for light source.

Results and discussion

Experiment

The spatio-temporal profiles of the liquid film interacting with a particle are derived from the videos captured by the high-speed cameras. The snapshot views in Fig. 3 (a,b) show the spreading behavior of the droplet near the particle in both views. We here define as t 0 [s] the moment at which the meniscus starts to envelop the particle observed from top. The velocity of the CL, U CL , is obtained from video analysis. In order to evaluate U CL , we apply the process as indicated in Fig. 4 in [START_REF] Mu | Sharp acceleration of a macroscopic contact line induced by a particle[END_REF] to the successive images observed from the side.

That is, we construct a space-time diagram along the substrate surface for tracking the position of the macroscopic contact line, X CL (t), from the successive side-view images. The velocity U CL is then the computed from the temporal evolution of X CL (t). The capillary number just before interacting with particle corresponds to Ca ∼ 3.7×10 -6 . Under the present condition for the experiment, the Reynolds and Ohnesorge numbers are Re ∼ 3.9 ×10 -3 and Oh ∼ 3.1 ×10 -2 , respectively. We thus neglect the inertia effect in the present system as [START_REF] Mu | Sharp acceleration of a macroscopic contact line induced by a particle[END_REF]. The profiles of liquid film in the symmetry plane (i.e., the plane z = z p ) extracted from side-view images are shown in Fig. 3 (c). By monitoring the temporal evolution of the spatial profiles of the meniscus formed around the particle, we obtain the same tendency as that in the previous study using smaller spherical particles of diameter D p = 50 µm [START_REF] Mu | Sharp acceleration of a macroscopic contact line induced by a particle[END_REF]. That is, when the spreading liquid film makes a contact with the particle, a concave meniscus, i.e., a meniscus with negative curvature in the x-y plane, forms around the particle foot, as seen in Fig. 3 (b,c).

From meniscus profiles we evaluate the pressure difference ∆p gl between liquid and gas phases based on the stress balance at interface,

∆p gl := p g -p l = γ ( ξ yy ( 1 + ξ 2 y ) 3/2 - 1 ξ ( 1 + ξ 2 y ) 1/2 ) , ( 4
)
where p g , p l represent the pressures on the gas and liquid sides of the interface, respectively, and x = ξ(t, y) is the meniscus profile in the symmetry plane. ξ y we have also assumed that the azimuthal variation of meniscus shape around the particle is locally negligible at the upstream and downstream sides. This assumption is supported by the fact that the CL downstream of the meniscus remains circular throughout meniscus development as shown in Fig. 3 (a). We estimate ∆p gl for an observed meniscus by seeking a value of ∆p gl that produces the solution of Eq. ( 4) fitting the best the experimental profile. This evaluation of Laplace pressure is validated from numerical results (Appendix).

The temporal variation of velocity of the CL, U CL are shown in Figure 4 (b).

The net pressure difference (Fig. 4 (c)) between upstream and downstream of the particle is assessed from ∆p ′(ud) := p where superscripts (up) and (down) indicate 'upstream' and 'downstream' with respect to the particle, respectively. These figures show a rapid increase of U CL at the beginning of meniscus formation (at 0 ≤ t -t 0 [s] ≤ 2.6), before a decrease (at t -t 0 [s] > 2.6) as the pressure difference drops. Here, U CL hits maximum value slightly after the highest pressure difference ∆p ′(ud) is reached.

This observation suggests that the pressure inside the meniscus is an important factor for the acceleration of the .

Simulation

Figure 5 shows numerical results of the spreading distance X CL and the U CL on the substrate with or without the particle. The spreading distance X CL is defined as the CL position along the line of (x, y = 0, z = z p ), and U CL is approximated by central differences. As shown in Fig. 5, on the flat substrate without the particle, the simulations correctly reproduce the liquid film spreads by following Tanner's law [START_REF] Tanner | The spreading of silicone oil drops on horizontal surfaces[END_REF], X CL ∝ t 1/7 and U CL ∝ t -6/7 . Similar to the experiments, the presence of the particle significantly alters the spreading of the liquid film: the velocity increases to reach a maximum value and then gradually decreases.

The theoretical analysis in the previous study [START_REF] Mu | Sharp acceleration of a macroscopic contact line induced by a particle[END_REF] found that the meniscus formed around the particle plays an important role in the acceleration phenomenon. We now investigate in detail what happens within the meniscus based on our numerical simulations, as shown in Fig. 6. In our simulation, the liquid film rises up around the particle right after the interaction, and the meniscus with negative curvature of the interface is spontaneously formed, which is similar to our experimental results. At the initial stage of meniscus formation (i.e., just after the liquid film contacts with the particle), the velocity on the meniscus surface exhibits an increase locally at the contact part (t = 15.0 ms).

As the liquid film rises around the particle, the CL on the substrate spreads forward (t = 17.0 ms), and then decreases (t = 19.0 ms). Here, the interface shape of the liquid film in the x-y cross section is shown in Fig. 7 to discuss the detail of its time evolution. The inset shows the comparison of the CL angular position φ CL around the particle and X CL . Because the CL of the spreading liquid film keeps advancing downstream of the particle, the distance over which the meniscus connects the substrate to the particle becomes longer and longer, and this causes the meniscus curvature (in the x-y plane) to decrease in magnitude. Such trend qualitatively agrees with the experimental results [START_REF] Mu | Sharp acceleration of a macroscopic contact line induced by a particle[END_REF]. We note the time scale of meniscus formation O(10 -2 s) is different from that in the experiments O(1 s), since the capillary number just in prior to the liquid contact with the particle in our simulation is larger than that in the experiment.

We monitor the temporal variation of the pressure distribution in the meniscus during its formation process, which is difficult to be observed in the experiment. Figure 8 illustrates the spatio-temporal variation of the pressure around the particle at different y. At t = 15.0 ms, the contact line makes a contact with the particle to form the meniscus. The lines that are almost parallel to the z axis at this instant are the contours of α = 0.5 at each y position. Only Vectors indicate the velocity on the meniscus surface. Note that the particle periphery at the same level in y as the front of the meniscus is illustrated in (a) for the sake of visibility.

the liquid near the substrate touches the particle, the contour at y = 0 is deformed, and those at y = 5 and 10 exhibit almost straight shape. Due to the negative curvature of the liquid film to form the meniscus as aforementioned, the pressure in the meniscus around the partciel foot becomes negative (y = 0).

As time elapses, the contact line along the particle surface rises to form larger meniscus, the contour becomes deforms at higher positions (see y = 5 at t = 18, and y = 5 and 10 at t = 21). As the meniscus develops around the particle, the absolute value of negative pressure becomes smaller in the meniscus because the curvature radius becomes larger as illustrated in Fig. 7. One finds that there exists negative gradient of the pressure around the particle in the x direction; the pressure in the upstream region against the particle becomes higher than the downstream region. This pressure gradient drives the liquid inside the meniscus toward the downstream side against the particle like a pump. It is obvious that the distribution of the pressure is symmetrical about the plane of z = z p . In the following we will discuss the pressure distribution in that plane.

As shown in Fig. 9 (a), the pressure inside the meniscus is lower than the atmospheric pressure at an initial stage (t = 15.0 ms). As the liquid film rises around the particle, the pressure inside the meniscus increases towards the atmospheric pressure because of the diminishing curvature caused by the ongoing advancement of the CL (t = 15.0 to 21.0 ms) as shown in Fig. 7 as well.

At the same time, the pressure difference between upstream and downstream sides is generated inside the meniscus at the particle foot (t = 18.0 ms), and it eventually diminishes (t = 21.0 ms). These results can be also explaneed by assuming that the Young-Laplace equation ( 4) is applicable to dynamical wetting. As aforementioned, the pressure difference between the menisci on the upstream and downstream sides drives the liquid forward along the particle foot.

The time evolution of the pressure difference ∆p = p (up) -p (down) inside the menisci on the upstream and the downstream sides and the CL velocity U CL are shown in Fig. 9 (b). Here, the pressures p (up) and p (down) are extracted at x up and x down , corresponding to (x, y, z) [µm] = (x p ± 15, 2.5, z p ). The CL velocity after the encounter with the particle (marked by the gray region in panel 9 (b) can be roughly classified into three regions of (I) acceleration, (II) deceleration, and (III) convergence. The volume function α up and α down at the two points x up and x down are also shown in the upper graph of Fig. 9 (b), in order to show when the advancing contact line has reached the downstream evaluation point x down . This happens at t = 14.5 ms. Before that instant, the pressure p (down) at x = x down corresponds to the ambient gas pressure and thus the time trace of ∆p carries the signature of p (up) . This pressure decreases when the meniscus starts to form on the upstream side of the particle, producing the minimum in the ∆p trace. Then, as the meniscus forms around the particle foot and reaches x down , ∆p becomes positive and reaches a maximum. This is due to the asymmetric nature of the meniscus, which displays a greater curvature on the downstream side. The CL velocity on the substrate reaches the maximum value slightly after it. From these results, we conclude that the meniscus plays a role as a 'pump' to suck a liquid from the bulk of the film and supply it behind the particle due to the pressure difference inside the meniscus between the front (upstream) and rear (downstream) of the particle.

In order to investigate the detailed behavior of the liquid inside the meniscus, the pressure, the absolute velocity U , and the velocities in x-, y-and z-directions in the liquid film around the particle on (x, y, z) = (x p , y, z) cross section are shown in Fig. 11. Initially the velocity in the x direction of the liquid film, U x , develops in the vicinity of the gas-liquid interface near the contact line on the particle, and then the large U x region transit to the z direction as the liquid film rises around the particle (region (i) → region (ii)). The small region near the particle foot (region (iii)) prevents the liquid to supply forward because of the viscous effect. The variation of the position where the velocity becomes maximum corresponds to the transition of the dominant driving force to 'pump' the liquid toward downstream region from the pressure difference within the meniscus in the main flow direction to the interfacial tension acting 275 on the interface with the curvature (red square in Fig. 6 (a)). One finds that the distribution of U x is almost the same as that of U , thus the pumping effect contributes the liquid transport almost along the x direction.

viscous dissipation??

Concluding remarks

The physical mechanism behind the local acceleration of the contact line of a spreading droplet due to its interaction with a microparticle as reported by Mu et al. [START_REF] Mu | Sharp acceleration of a macroscopic contact line induced by a particle[END_REF][START_REF] Mu | Control of local wetting by microscopic particles[END_REF][START_REF] Mu | Quick liquid propagation on a linear array of micropillars[END_REF] was investigated numerically and experimentally. The aim of the numerical simulations was to access the physical data such as the pressure and velocity fields in the liquid film forming around the particle that were not accessible experimentally. The experiments were used to validate the numerical simulations through the time-dependent shape taken by the meniscus around the particle's foot. Numerical simulations reproduce the experimental observations of an accelerating microscopic contact line around the particle. In particular, the simulations similarly show that the particle-liquid film interaction enables the microscopic contact line of a spreading of droplet to accelerate and overcome locally the commonly known power laws of Tanner or Lopez [START_REF] Lopez | Spreading kinetics of liquid drops on solids[END_REF][START_REF] Tanner | The spreading of silicone oil drops on horizontal surfaces[END_REF][START_REF] Cazabat | Dynamics of wetting: effects of surface roughness[END_REF]. This result is of primary importance as it provides a numerical tool to predict the wetting dynamic of a substrate by a liquid with respect to topological conditions.

The simulations were limited to small capillary numbers Ca = µU CL /γ ∼ O (10 -4 ) to O (10 -6 ). In order to unravel the physics behind the contact line acceleration, a particular focus was made in studying the process of meniscus formation right after the contact of the contact line with the particle. Three regimes of behaviour for the contact line were identified: an accelerating, a decelerating, and a converging regime [START_REF] Kondo | Acceleration of the macroscopic contact line of a droplet spreading on a substrate after interaction with a particle[END_REF]. It is found that the characteristic time for the meniscus formation around the particle becomes of the order of O(10 -2 s) which is much longer than that reported in the literature for the fully wetting of a smooth substrate by a liquid [START_REF] Biance | First steps in the spreading of a liquid droplet[END_REF]. Furthermore, the relaxation of the Laplace pressure occurs within a long characteristic time of the order of µD p /2γ ∼ O(10 -6 s). This long characteristic time is linked to the significant limitation of the liquid within a liquid film. Indeed, the pressure difference between the upstream and downstream regions of the meniscus around the particle is behind the pumping mechanism that delivers the liquid from the bulk to the contact line behind the particle. Once the pressure difference vanishes the local contact line around the particle decelerates till it is matched by the main contact line of the spreading droplet that follows the Tanner's law [START_REF] Tanner | The spreading of silicone oil drops on horizontal surfaces[END_REF]. Those findings leads a comprehensive description of the enhancement of the droplet spreading on the microstructured solid surfaces [START_REF] Courbin | Imbibition by polygonal spreading on microdecorated surfaces[END_REF][START_REF] Xiao | Microscale liquid dynamics and the effect on macroscale propagation in pillar arrays[END_REF][START_REF] Obara | Imbibition of a textured surface decorated by short pillars with rounded edges[END_REF][START_REF] Wang | Surface structure determines dynamic wetting[END_REF][START_REF] Kim | Dynamics of hemiwicking[END_REF]. We will apply these findings and the numerical simulations to optimize the wetting of a microstructured surface. It will, in particular, be used to unravel the physical mechanism behind the dynamics of successive accelerations experienced by a contact line through an array of pillars reported by [START_REF] Mu | Quick liquid propagation on a linear array of micropillars[END_REF][START_REF] Spruijt | Liquid filmification from menisci[END_REF] where the local acceleration of the contact line can reach up to 100 times that of a spreading droplet on a smooth substrate. Such understanding and data are of primary importance for designing economical and low-energy consuming liquid transport devices, for example, in optimal designing of hemiwicking structures for heat-pipe wicks for improving the performance of lithium-ion batteries by structuring electrode surfaces [START_REF] Pfleging | A review of laser electrode processing for development and manufacturing of lithium-ion batteries[END_REF].

with a constant Laplace pressure is a solution of the Eqs. 4. Note that ∆p gl in Eqs. 4 is the average of the pressures at the meniscus surface ξ (y) given by the following:

  )-(3) are, respectively, the continuity equation of incompressible flows, the Navier-Stokes equation, and an advection equation of volume function α. The symbols U , p, τ , F SV , U r represent velocity, pressure, viscous stress, force induced by the interfacial tension, and relative velocity between liquid and gas phases, respectively. The volume fraction α of the liquid phase takes

Fig. 1 .

 1 The length and grid number of the domain are (L x , L y , L z ) [µm] = (300, 80, 300) and (N x , N y , N z ) = (150, 40, 150), respectively. The origin of the coordinate system is set on the substrate surface at the upstream edge of the computational domain. No-slip and no-penetration conditions are imposed at the substrate surface (x, y = 0, z). The top plane (x, y = L y , z) and the side planes ((x, y, z = 0); (x, y, z = L z )) are set to mirror boundaries. A spherical solid boundary representing a particle of diameter D p = 50 µm is positioned at the center of the substrate (x p = 150, y p = 25, z p = 150) µm. The wettability of the substrate and that of the particle are prescribed by contact angles θ s = 5 • and θ p = 20 • , respectively.

Figure 1 :

 1 Figure 1: Computational domain of spreading liquid film on horizontal substrate with a single particle.

Figure 2 :

 2 Figure 2: Experimental apparatus for observation of droplet spreading on substrate. Simultaneous observation from above and side is realized with synchronized two high-speed cameras.

Figure 3 :

 3 Figure 3: Meniscus formation around a particle of a diameter Dp = 200 µm observed in an experiment. Images in top view (a) and in side view (b) and meniscus profiles (c) extracted from the side view images.

Figure 4 :

 4 Figure 4: memo by IU: I plan changing the figure for (a); temporal variations of the meniscus are also illustrated as Mu et al.(JFM 2017) in order to reflect the comment 4 by the reviewer1. What do you think? (a) Example of evaluating profiles of meniscus (red) formed around particle from experimental result. Inset: temporal variations of the contact-line positions along the particle, ϕ, and that along the substrate, X CL . Temporal variations of (b) CL velocity U CL and (c) evaluated pressure difference between menisci on upstream and downstream sides ∆p ′(ud) l :=

Figure 5 :

 5 Figure 5: Temporal variations of the CL position X CL and velocity U CL on substrates. The CL position is derived on the line parallel to the x direction passing under the particle foot on the substrate as illustrated in inset. Note that the CL position is not well defined in the gray area, in which the CL is traveling around the particle foot.

Figure 6 :

 6 Figure 6: Meniscus formation around a particle of a diameter 50 µm observed in a simulation.

Figure 7 :

 7 Figure 7: Temporal variation of meniscus profile around a particle on plane of z = zp by numerical simulation. First profile corresponds to the meniscus surface at t = 0.015 s and the time interval between profiles is of 0.0002 s. Inset illustrates temporal variations of angular position φ CL of the meniscus front on the particle and the position of CL on the substrate X CL .

Figure 8 :

 8 Figure 8: Spatio-temporal distributions of pressure at different height y.

Figure 9 :

 9 Figure 9: (a) Temporal variation of pressure around a particle on plane of z = zp. (b) Temporal variations of pressure difference between upstream (at xup as shown in (a)) and downstream (x down ) sides in meniscus, and CL velocity U CL . The monitoring points inside the meniscus on the upstream and downstream sides are defined as xup = (135, 2.5, 150) [µm], x down = (165, 2.5, 150) [µm], respectively. Top of (b) indicates the temporal variation of the volume fraction α measured at xup (αup) and at x down (α down ). The position of CL on the substrate cannot be well defined in the gray area.

Figure 10 :

 10 Figure 10: Spatial-temporal variations of pressure, absolute velocity, and velocity component in each direction inside the meniscus formed around the particle on plane of (x, y, z) = (x, y, zp).

Figure 11 :

 11 Figure 11: viscous dissipation??

Figure A. 1 : 1 2

 11 Figure A.1: (a) Profiles of a downstream side of meniscus developing around a spherical particle at t = 15, 15.2, 15.4, 15.6, 15.8, 16, 16.2, and 16.4 ms. Profiles determined in a numerical simulation and theoretical ones are compared. The latter profiles are used to estimate Laplace pressure ∆p gl . (b) Temporal variation of the pressure difference ∆p gl at gas-liquid interface of downstream side of meniscus by estimation based on Laplace pressure and simulation.

  Figure A.1 (b) shows the comparison between the pressure at the gas-liquid interface by estimated from the profiles of Fig. A.1 (a) and computed by simulation. A slight decrease of ∆p gl with y is observed and trends are consistent. From the above, this Laplace pressure measurement is valid enough to indicate the relationship between the CL velocity and the pressure difference inside the meniscus, which are estimated by the meniscus profiles, in the experiment.

Acknowledgment

This work was partially supported from Japan Society for the Promotion of Science (JSPS) by Grant-in-Aid for Challenging Exploratory Research (grant number: 16K14176) and by Grant-in-Aid for Scientific Research (grant number: 19H02083). Authors acknowledge a financial support from MAE (French Ministry of Foreign Affairs and International Development) as well as from JSPS through the PHC Sakura program 'Wetting dynamics in the presence of particles.' This work benefited from the support of the project FEFS ANR-CE21-2018 of the French National Research Agency (ANR).

Appendix

Validation of Laplace pressure measurement

Laplace pressure ∆p gl of a meniscus developing around a spherical particle is estimated from its profile in the x -y plane, forms around the particle foot, as seen in Fig. 3 (b,c). The profile x = ξ (y) of an axisymmetric meniscus