1. Introduction
Coating flows are ubiquitous in nature and industrial applications. Nature provides astonishing examples of the capability of coating flows to modify the topography of the substrate via chemical and thermodynamic reactions. The structures that can be observed in limestone caves, known as speleothems, are characterized by a morphogenesis that is related to the hydrodynamic instability of a coating flow (Short et al. Reference Short, Baygents, Beck, Stone, Toomey III and Goldstein2005; Meakin & Jamtveit Reference Meakin and Jamtveit2010; Camporeale Reference Camporeale2015; Bertagni & Camporeale Reference Bertagni and Camporeale2017). These fascinating structures originate from the interaction between hydrodynamics and chemistry. The control of the instability related to coating processes is an important task in industrial applications as many fabrication processes involve the presence of a thin film flowing on a substrate (Weinstein & Ruschak Reference Weinstein and Ruschak2004). Thin elastic shells of constant thickness can be fabricated by polymerization of the film, as performed in Lee et al. (Reference Lee, Brun, Marthelot, Balestra, Gallaire and Reis2016). Marthelot et al. (Reference Marthelot, Strong, Reis and Brun2018b) showed a remarkable example of control of the flow instability to produce textured surfaces, by rotation of a cylindrical substrate.
The Rayleigh–Taylor instability is a phenomenon that occurs when a heavier fluid is placed above a lighter one. When a horizontal flat interface is considered, under the only effect of gravity, all wavelengths are unstable (Rayleigh Reference Rayleigh1882; Taylor Reference Taylor1950). The introduction of capillary effects bounds the range of unstable wavelengths (Chandrasekhar Reference Chandrasekhar2013). When an upper wall confines the overhanging fluid, the resulting pattern is characterized by lenses arranged in hexagonal or square arrays (Fermigier et al. Reference Fermigier, Limat, Wesfreid, Boudinet and Quilliet1992). The lenses may saturate for small enough initial thickness (Marthelot et al. Reference Marthelot, Strong, Reis and Brun2018a), or algebraically grow (Yiantsios & Higgins Reference Yiantsios and Higgins1989; Lister, Rallison & Rees Reference Lister, Rallison and Rees2010), eventually resulting in dripping droplets.
The problem of the dynamics of a thin film is usually studied within the context of the lubrication approximation. The model assumes much larger characteristic lengths in the directions which lay along the substrate than in the normal-to-the-substrate direction (Ruschak Reference Ruschak1978; Wilson Reference Wilson1982; Babchin et al. Reference Babchin, Frenkel, Levich and Sivashinsky1983; Weinstein & Ruschak Reference Weinstein and Ruschak2004).
In the case of an inclined substrate, the route from a flat film towards dripping drops still needs to be analysed. When the substrate is tilted with respect to the horizontal direction, the gravity component parallel to the substrate creates a flow. In this work, we consider a configuration with a permanent influx, in opposition to the case of cylindrical and spherical substrates in which a transient release of fluid is studied (Balestra et al. Reference Balestra, Kofman, Brun, Scheid and Gallaire2018a; Balestra, Nguyen & Gallaire Reference Balestra, Nguyen and Gallaire2018b). A strong modulation of the thickness along the direction perpendicular to the flow (spanwise direction) is identified as rivulet formation (Charogiannis et al. Reference Charogiannis, Denner, van Wachem, Kalliadasis, Scheid and Markides2018). The presence of a predominant rivulet pattern when the inertia of the fluid is negligible was experimentally observed by Charogiannis et al. (Reference Charogiannis, Denner, van Wachem, Kalliadasis, Scheid and Markides2018). Similar rivulet patterns were observed by Rietz et al. (Reference Rietz, Scheid, Gallaire, Kofman, Kneer and Rohlfs2017), in an experimental set-up where gravity was replaced by centrifugal acceleration. Lerisson et al. (Reference Lerisson, Ledda, Balestra and Gallaire2019) showed that a state characterized by lenses travelling on the rivulets may emerge, depending on the inclination angle and flow rate.
The stability analysis was performed by linearizing the flow equations around a constant thickness, revealing that the flat film solution is always unstable to perturbations (Brun et al. Reference Brun, Damiano, Rieu, Balestra and Gallaire2015). These authors found experimentally a link between dripping and the absolute instability of the flow, modelled with the one-dimensional lubrication equation. The model was refined introducing inertial and viscous extensional stresses (Scheid, Kofman & Rohlfs Reference Scheid, Kofman and Rohlfs2016; Kofman et al. Reference Kofman, Rohlfs, Gallaire, Scheid and Ruyer-Quil2018). These authors showed that the occurrence of the absolute instability does not predict the dripping satisfactorily.
In Lerisson et al. (Reference Lerisson, Ledda, Balestra and Gallaire2020) an experimental set-up able to continuously feed an inclined planar substrate with fluid was presented. Using a very viscous fluid such that inertial effects are negligible, the natural emergence of elongated, streamwise-oriented, steady patterns was observed. A detailed analysis of the appearance of these so-called rivulets was then performed, both when a spanwise-periodic forcing is imposed at the inlet and when the rivulets emerge naturally from the lateral boundaries of the experiment. The forced dynamics revealed that there is a narrow range of attainable spacings of rivulets. The nonlinear simulations agreed with the thickness measured in experiments, observing steady and streamwise-invariant rivulet states, periodic along the spanwise direction. The one-dimensional and saturated rivulet profile was recovered by simple static arguments, i.e. the equilibrium between capillary effects and hydrostatic pressure gradient (Roman, Gay & Clanet Reference Roman, Gay and Clanet2020; Zaccaria et al. Reference Zaccaria, Bigoni, Noselli and Misseroni2011; Duprat & Stone Reference Duprat and Stone2015). The correct shape was obtained imposing the local flow rate along the direction transverse to the rivulet profile.
In this work, we aim at rationalizing the observations of steady rivulet patterns by investigating the intrinsic rivulets selection and their stability.
The paper is organized as follows. We first introduce an experimental visualization for the evolution of the film when the inlet is steadily forced along the spanwise direction. A numerical study for an initial condition that mimics these experimental conditions, namely a regular pattern of sinusoidal perturbations in the spanwise direction, is performed. Periodic boundary conditions in all in-plane directions are imposed. The experimental and numerical results are then rationalized by a secondary stability analysis. We perturb the one-dimensional rivulet profile along the streamwise direction with a normal mode expansion and obtain a dispersion relation characterizing the secondary growth of lenses. We thus present a comparison of the secondary stability study with experimental measurements of the spatial amplification of disturbances over steady rivulets. The last section is devoted to the study of the emergence of rivulets from a flat film when the film is impulsively perturbed. We introduce a qualitative experimental visualization when the film is excited by a localized perturbation in the thickness, the results of which are numerically reproduced. A weakly nonlinear model is eventually proposed to rationalize these observations.
2. Experimental apparatus
The experimental apparatus is the same as that described in Lerisson et al. (Reference Lerisson, Ledda, Balestra and Gallaire2020) (see figure 1). The substrate is an orientable glass plate of length 
$\hat {L}_x=600$ mm and width 
$\hat {W}_i=300$ mm, whose angle with respect to the vertical is varied from 
$\theta =20^{\circ }$ to 
$\theta =80^{\circ }$. The fluid is silicon oil (Bluestar Silicons 47V1000) of density 
$\rho =974$ kg m
$^{-3}$, viscosity 
$\mu =1089$ mPa s and surface tension coefficient 
$\gamma =21$ mN m
$^{-1}$. The oil is injected through a horizontal rectangular opening of a reservoir and flows beneath the substrate. The flow rate is driven by the height difference with another reservoir that creates a hydrostatic pressure gradient. The flow rate can be varied by changing the height difference of the two reservoirs. The system is designed in such a way that it is possible to steadily modify the inlet condition in the spanwise direction by adding a sinusoidal or a comb-like blades (see sketches in figure 1a). The sinusoidal blade is placed below the inlet with an angle of 
$30^\circ$ with respect to the substrate, and the fluid fills the gap between the glass and the blade. Systematic measurements of the thickness give a thickness perturbation amplitude of 
$\simeq 250\ \mathrm {\mu }$m. The comb-like blade presents teeth of thickness 
$\hat {t}_t=1$ mm, streamwise size of 
$\hat {l}_{\textrm {d} t}=5$ mm and spanwise size of 
$2$ mm. The teeth occlude the inlet and the fluid covers them by capillarity.
Figure 1. (a) Sketch of the experimental apparatus with the detail of the sinusoidal and comb-like blades for the steady forcing at the inlet along the spanwise direction. (b) Photo of the experimental apparatus.
The volumic flow rate 
$q$ is measured by weighing the oil leaving the substrate for 180 s. We define the equivalent Nusselt thickness 
$h_N$ as well as the reduced capillary length 
$l_c^*$:
\begin{equation} h_N=\left(\frac{3\nu q}{\hat{W}_i g \cos(\theta)}\right)^{1/3},\quad l_c^*=\frac{l_c}{\sqrt{\sin(\theta)}},\end{equation}
where 
$l_c=\sqrt {\gamma /(\rho g)}$ is the capillary length, g the gravity acceleration and $\nu=\mu/\rho$ the kinematic viscosity. We define a coordinate system 
$(\hat {x},\hat {y}, \hat {z})$, where 
$\hat {x}$, 
$\hat {y}$ and 
$\hat {z}$ are respectively the streamwise, spanwise and normal-to-the-substrate directions.
We employ a qualitative visualization technique based on shadowgraphs that are constructed looking at the distortion of the rays coming from a point light source through the liquid film. The surface deformation will focus or defocus the initially homogenous light and forms patterns that are highly sensitive to slight deformations. The combination of small and large deformations (Settles Reference Settles2001; Moisy, Rabaud & Salsac Reference Moisy, Rabaud and Salsac2009) within the same experiment makes the visualization impossible to relate to quantitative measurements of the thickness amplitude. However, the experiment gives access to the phases of perturbations, and thus to qualitative observations of the emerging pattern.
We measure the film thickness using the confocal chromatic sensor STIL-CCS located on the upper part of the glass plate. We choose an acquisition rate of 100 Hz. The position of the sensor can be adjusted in the normal-to-the-substrate and spanwise directions once the streamwise location is selected.
3. Observations of the secondary stability and instability of rivulets
3.1. Experimental observations
In this section, we briefly present selected results from the study of Lerisson et al. (Reference Lerisson, Ledda, Balestra and Gallaire2020) in the presence and absence of the spanwise inlet forcing devices shown in figure 1(a). Figure 2(a) shows a film thickness distribution obtained using an absorption technique (reproduced from Lerisson et al. Reference Lerisson, Ledda, Balestra and Gallaire2020). The inlet spanwise thickness profile is amplified, and streamwise-saturated and steady rivulets are observed downstream. The saturated rivulets are periodic along the spanwise direction. There is a narrow range of attainable spacings, when the inlet is forced, around the value 
$\hat {L}_r=2\pi \sqrt {2}l_c^*$ (value shown in figure 2a), i.e. the most amplified wavelength in the dispersion relation of the flat film. Interestingly, even in the absence of the spanwise inlet forcing devices, the predominant spacing of the emerging rivulet structures is 
$\hat {L}_r$ (see figure 2b).
Figure 2. (a) Film thickness for 
$\theta = 39^{\circ }$ and 
$h_N = 1515\ \mathrm {\mu }\text {m}$ (
$u=1.5$), steady inlet forcing with the sinusoidal blade at the wavelength 
$\hat {L}_f = 2\pi \sqrt {2}l_c/\sqrt {\sin \theta }$, extracted from Lerisson et al. (Reference Lerisson, Ledda, Balestra and Gallaire2020). The thickness is measured with the absorption method and normalized by the flat film thickness 
$h_N$. The in-plane distances are normalized by the reduced capillary length 
$l_c/\sqrt {\sin \theta }$. (b) Typical rivulet pattern in the absence of the inlet forcing devices (figure 1a), 
$\theta =20^\circ$.
However, far downstream in figure 2(a), oscillations appear on the rivulet profiles. These oscillations are amplified and rivulets carrying travelling lenses are observed, for these values of angle and flow rate.
3.2. Numerical observations
The aim of this section is to numerically study the emerging patterns for an initial condition that mimics the experimental conditions described in the previous section.
We consider a gravity-driven thin film of viscous Newtonian fluid flowing under a planar substrate inclined with respect to the vertical with an angle 
$\theta$. We introduce the following non-dimensionalization:
\begin{equation} x=\hat{x}/l_c^*; \quad y=\hat{y}/l_c^*; \quad z=\hat{z}/h_N; \quad t=\hat{t}/\tau^*, \end{equation}
where 
$\tau ^*={\nu l_c^2}/{h_N^3 \sin ^2(\theta ) g}$ is the characteristic time scale of the Rayleigh–Taylor instability. The numerical model for the evolution of the film thickness 
$h$ is the lubrication equation in which the complete expression of the curvature is retained (Ruschak Reference Ruschak1978; Wilson Reference Wilson1982; Weinstein & Ruschak Reference Weinstein and Ruschak2004):
\begin{equation} \partial_t h +u h^2 \partial_x h +\tfrac{1}{3} \boldsymbol{\nabla} \boldsymbol{\cdot} [ h^3 \boldsymbol{\nabla} h + h^3 \boldsymbol{\nabla} \kappa]=0,\end{equation}
where 
$\boldsymbol {\nabla }$ operates in the 
$(x,y)$ directions, 
$u= \cot {(\theta )} \tilde {l}_c^*$ and 
$\tilde {l}_c^*= {l}_c^*/h_N$. The linear advection velocity 
$u$ corresponds to the surface film velocity at which linear interface thickness perturbations with respect to a flat condition are advected downstream (Brun et al. Reference Brun, Damiano, Rieu, Balestra and Gallaire2015). In physical quantities, an increase of the parameter 
$u$ implies a decrease of the flow rate (since 
$u$ is inversely proportional to 
$h_N$) or 
$\theta$. The curvature 
$\kappa$ reads
\begin{equation} { \kappa}=\frac{\partial_{{x}{x}} {h}(1+(\partial_{{y}}{h})^2) +\partial_{{y}{y}} {h}(1+(\partial_{{x}}{h})^2) - 2\partial_{{x}{y}}{h}\partial_{{x}}{h}\partial_{{y}}{h}}{(1+(\partial_{{x}}{h})^2 +(\partial_{{y}}{h})^2)^{3/2}}. \end{equation} The two-dimensional equation is implemented in COMSOL Multiphysics. We use the built-in finite elements method solver, exploiting cubic elements with Lagrangian shape functions and a fully implicit time solver. The largest mesh element size is half of the reduced capillary length 
$\tilde {l}_c^*$. The domain size is 
$L_x\times L_y$, where 
$L_x=231$ and 
$L_y=106$, leading to approximately 
$50\,000$ elements. A convergence analysis was performed, showing that convergence is achieved for this characteristic size of the elements. This characteristic element size was also validated by the experimental and numerical comparisons in Lerisson et al. (Reference Lerisson, Ledda, Balestra and Gallaire2020). The equations are solved for the variables (
$h,\kappa$). For all the considered cases, periodic boundary conditions are used.
Experimentally, in the absence of the spanwise inlet perturbation device described in figure 1(a), the rivulet spacing is the one dictated by the most amplified mode in the flat film dispersion relation, i.e. 
${L}_r=2\pi \sqrt {2}$ (Lerisson et al. Reference Lerisson, Ledda, Balestra and Gallaire2020). We numerically study the nonlinear time evolution when a streamwise-invariant sinusoidal initial condition is considered. We choose as initial condition a sinus of wavelength 
${L}_r=2\pi \sqrt {2}$:
\begin{equation} h(x,y,t=0)=\bar{h}_N\left(1+{A}\cos\left(\frac{2\pi y}{L_r}\right)\right), \end{equation}
where 
${A}=10^{-2}$, and 
$\bar {h}_N=0.54$ is the initial value of the thickness that gives, for a pure streamwise saturated structure, the same local flow rate in the streamwise direction as a flat film of thickness 
$h=1$ (§ 5.3 in Lerisson et al. Reference Lerisson, Ledda, Balestra and Gallaire2020).
We introduce the moving reference frame at the linear advection velocity 
$u$
$(\xi =x-u t,y)$. Figure 3 shows the evolution of the thickness with time for (a) 
$u=5.45$ and (b) 
$u=1.5$. For visualization purposes, we focus in the region 
$\xi \in [-8\pi \sqrt {2},8\pi \sqrt {2}]$ and 
$y \in [-6\pi \sqrt {2},6\pi \sqrt {2}]$. In both cases, the streamwise invariant initial condition is amplified and reaches, at 
$t=800$, a saturated state in the streamwise direction. For (a) 
$u=5.45$, we do not observe any further evolution of the pattern for 
$t>800$. For (b) 
$u=1.5$, at 
$t=800$ the rivulet profiles saturate. For 
$t>800$, however, streamwise thickness perturbations grow, and at 
$t=1200$ the flow is characterized by lenses travelling on the rivulets.
Figure 3. Nonlinear response in the case of a streamwise-invariant sinusoidal initial condition, for (a) 
$u=5.45$ and (b) 
$u=1.5$. From left to right: 
$t=1000$, 
$t=1200$. Results are reported in the moving reference frame at the linear advection velocity 
$(\xi =x-u t,y)$.
The streamwise-invariant sinusoidal initial condition is amplified leading to a rivulet pattern saturated in space and time, periodic along the spanwise direction. The absence (respectively presence) of observable streamwise perturbations on the rivulet profiles at high (respectively low) values of 
$u$ suggests that the stability of the rivulet profile to streamwise perturbations may be directly related to the linear advection velocity.
The experimental observations of predominant spanwise-periodic rivulet patterns and the occurrence of lenses on the rivulets are confirmed by the nonlinear simulations with periodic boundary conditions. In the following, we aim at rationalizing the emergence of predominant rivulet structures and their destabilization.
4. Secondary stability analysis of rivulets
In § 3.1 it was experimentally shown that rivulet structures grow in the domain and saturate to a steady and spanwise-periodic state, invariant along the streamwise direction. However, for low values of 
$u$ and at large distances from the inlet, the rivulet profile undergoes an instability and travelling lenses emerge on the rivulet structures, as shown in figure 2(a). The saturation of the rivulet structures and the occurrence of lenses was also observed in the nonlinear numerical simulation of figure 3. No lateral interactions between rivulets are observed as the lenses grow. Here, we study the robustness of the saturated rivulet pattern via a secondary stability analysis. We first introduce the steady, streamwise-invariant and spanwise-periodic rivulet profile 
$H_r(y)$, and then we focus on its local stability properties when perturbed along the streamwise direction, 
$x$. The validity of the local stability analysis is limited to the regions where steady and one-dimensional rivulets are observed.
4.1. Base flow
In this section, we define the saturated rivulet profile 
$H_r(y)$, serving as base flow for the local stability analysis. The numerical base flow is the large-time solution (
$t=10\ 000$) of the one-dimensional model presented in § 5 of Lerisson et al. (Reference Lerisson, Ledda, Balestra and Gallaire2020). The profile, of periodic wavelength 
${L}_r$, is given by a one-dimensional model in which the flow rate in the streamwise direction coincides with the one of a flat film of thickness 
$h= 1$, leading to a mean value 
$\bar {h}_N=0.54$ of the thickness of the rivulet. The numerical procedure revealed that the rivulet profile is slowly saturating to a steady state 
$H_r(y)$. In figure 4, we report the numerical periodic profile at 
$t=10\ 000$ (solid line) used for the stability analysis. The rivulet is characterized by a central lobe of large thickness that saturates to a steady profile described by the pendulum equation (red circles in figure 4), while the side lobes (of very low thickness) are slowly draining with a power law 
$t^{-1/2}$ (Lister et al. Reference Lister, Rallison and Rees2010). It is remarkable that, with the considered non-dimensionalization, the profiles are independent of 
$u$, i.e. there is a unique rivulet shape (Lerisson et al. Reference Lerisson, Ledda, Balestra and Gallaire2020). The numerical profile agrees well with the experimental results (dots in figure 4) and can therefore be safely used as base flow 
$H_r(y)$ for the stability analysis.
Figure 4. Periodic rivulet profile (black line) used for the stability analysis, compared with the results of the pendulum equation of § 5.4 of Lerisson et al. (Reference Lerisson, Ledda, Balestra and Gallaire2020) (red circles), and with the experimental results for three central rivulets (grey dots), from Lerisson et al. (Reference Lerisson, Ledda, Balestra and Gallaire2020), for 10 transverse measurements at two different streamwise locations, at 
$\theta =39^\circ$ and different 
$h_N$. The red dashed line denotes the mean thickness 
$\bar {h}_N$ of the rivulet.
4.2. Dispersion relation
Following the classical approach of the local stability analysis, we consider as a base state the single, spanwise-periodic and steady rivulet 
$H_r(y)$ described in § 4.1. The quasi-steadiness of the rivulet profile allows us to neglect the slow evolution of the side lobes at long times and thus to consider a normal mode expansion in time and along the direction in which the base state is invariant, i.e. the streamwise direction 
$x$ (Schmid, Henningson & Jankowski Reference Schmid, Henningson and Jankowski2002). The spanwise periodicity governing the base state 
$H_r(y)$ is also enforced on the perturbation. The following normal mode decomposition is therefore used:
\begin{equation} h(x,y,t)=H_r(y)+\varepsilon \tilde{\eta}(x,y,t)=H_r(y)+\varepsilon \eta(y)\, \textrm{e}^{\mathrm{i}(k_x x -\omega t)}, \quad {\varepsilon \ll 1,} \end{equation}
where 
$\tilde {\eta }$ is the thickness perturbation with respect to the base flow profile 
$H_r(y)$. By considering the two-dimensional nonlinear equation (3.2) and introducing the normal mode decomposition (4.1), one obtains, up to 
$ {O}(\varepsilon )$
\begin{align} &\varepsilon \partial_t \tilde{\eta} + \varepsilon u H_r^2\partial_x \tilde{\eta} + \frac{1}{3} \partial_y \left[ H_r^3\left(\frac{\mathrm{d}H_r}{\mathrm{d} y} +\frac{\mathrm{d}\kappa_{(0)}}{\mathrm{d} y}\right)+ \varepsilon H_r^3 \partial_y\tilde{\eta} \right.\nonumber\\ &\quad \left. +\,\varepsilon H_r^3 \partial_y \tilde{\kappa}_{(1)} + 3\varepsilon H_r^2\left(\frac{\mathrm{d}\kappa_{(0)}}{\mathrm{d} y} + \frac{\mathrm{d}H_r}{\mathrm{d} y}\right) \tilde{\eta} \right] + \frac{\varepsilon}{3} \partial_x \left[ H_r^3 (\partial_x \tilde{\kappa}_{(1)} +\partial_x \tilde{\eta}) \right] =0, \end{align}
where 
$\kappa _{(0)}$ is the base flow curvature, i.e. (3.3) evaluated for the base flow 
$H_r(y)$, 
$\kappa _{(0)}= ({\mathrm {d}^2H_r}/{\mathrm {d}y^2})/{(1+({\mathrm {d}H_r}/{\mathrm {d}y})^2)^{3/2}}$. Furthermore, 
$\tilde {\kappa }_{(1)}$ is the first-order term of the curvature, i.e. the Jacobian of the curvature evaluated in the base flow and applied to 
$\tilde {\eta }$ (
$\tilde {\kappa }_{(1)}= [\partial _{\tilde {\eta }}\kappa (H_r)]{\tilde {\eta }}$). The full expression of the operator 
$\partial _{\tilde {\eta }}\kappa (H_r)$ is reported in appendix A. Deriving this expression with respect to 
$x$ and 
$y$, we obtain 
$\partial _x \tilde {\kappa }_{(1)}=\mathrm {i}k_x {\kappa }_{(1)}(y)\exp ({\mathrm {i}(k_x x -\omega t)})$ and 
$\partial _y \tilde {\kappa }_{(1)}=({\mathrm {d} {\kappa }_{(1)}}/{\mathrm {d}y})(y)\exp ({\mathrm {i}(k_x x -\omega t)})$.
At 
$ {O}(1)$ the base flow equation is recovered, while at 
$ {O}(\varepsilon )$ one obtains the following evolution equation for the perturbation:
\begin{align} &-\mathrm{i} \omega \eta+\mathrm{i}k_xu H_r^2\eta + \frac{1}{3} \frac{\mathrm{d}}{\mathrm{d} y} \left[ 3H_r^2 \left(\frac{\mathrm{d}H_r}{\mathrm{d} y}+ \frac{\mathrm{d}\kappa_{(0)}}{\mathrm{d} y}\right)\eta \right.\nonumber\\ &\quad +\left. H_r^3 \left(\frac{\mathrm{d}{\kappa}_{(1)}}{\mathrm{d} y}+ \frac{\mathrm{d}\eta}{\mathrm{d} y}\right) \right] - \frac{1}{3} k_x^2\left[ H_r^3 \left({\kappa}_{(1)} +\eta\right) \right] =0, \end{align}
which is the dispersion relation 
$D_r(\omega , k_x)=0$. The base flow 
$H_r(y)$ can be perturbed by (i) imposing the streamwise wavenumber 
$k_x\in \mathbb {R}$ and looking at the temporal evolution through the complex frequency 
$\omega \in \mathbb {C}$ (temporal stability analysis) or (ii) imposing a temporal forcing of real frequency 
$\omega$ and looking at the spatial amplification of the perturbation, embodied by the complex spatial wavenumber 
$k_x\in \mathbb {C}$ (spatial stability analysis).
The numerical implementation of (4.3) is performed in MATLAB by a spectral collocation Fourier method. Once discretized, the eigenfunction problem (4.3) becomes an eigenvalue problem. The temporal and spatial stability analyses are respectively solved using the built-in MATLAB functions eig and polyeig. Numerical convergence is achieved for 100 collocations points. A preparatory analysis on the numerical rivulet profile 
$H_r(y)$ used as base flow for the stability analysis revealed a variation of the eigenvalues of the order of the numerical discretization, as long as 
$t>5000$.
4.3. Temporal stability analysis
In this section, we report the results for the temporal stability analysis. Positive (respectively negatives) values of the temporal growth rate 
$\mathrm {Im}(\omega )=\omega _i$ denote unstable (respectively stable) wavenumbers. A preliminary analysis on the spectrum revealed that all the eigenvalues have negative 
$\omega _i$ for all 
$k_x$, except one that is analysed in the following.
In figure 5(a) we report the variation of 
$\omega _i$ with 
$k_x$, for different values of 
$u$. The dispersion relations are characterized by a local maximum associated with the dominant wavenumber, and by a value of the wavenumber beyond which the temporal growth rate is negative (the cutoff wavenumber), i.e. perturbations with wavenumber larger than the cutoff are damped. Rivulets are strongly stabilized as the value of 
$u$ increases. For 
$u=1$ the growth rate 
$\omega _i$ presents its maximum value at a dominant wavenumber close to 
$k_x=0.56$, while the cutoff wavenumber 
$k_x^{cut}=0.8$. An increase of 
$u$ quickly quenches large wavenumbers. Both the dominant growth rate and the cutoff wavenumber decrease. At 
$u=5$, 
$k_x^{cut} \sim 10^{-2}$, with 
${\rm max}(\omega _i) \sim 10^{-3}$. For these values of 
$u$, the unstable wavelengths are of the order of one hundred reduced capillary lengths. The real frequency 
$\mathrm {Re}(\omega )=\omega _r$ increases slightly less than linearly with 
$k_x$ (figure 5b). The resulting phase velocities 
$\omega _r/k_x$ increase as 
$u$ increases.
Figure 5. (a) Temporal growth rate 
$\omega _i$ and (b) real frequency 
$\omega _r$ as functions of the streamwise wavenumber 
$k_x$, from the temporal stability analysis, for 
$u=1$ (blue line), 
$u=1.5$ (red line), 
$u=2$ (yellow line), 
$u=2.5$ (purple line), 
$u=3$ (green line), 
$u=5$ (light blue line).
In figure 6(a) we show the real (dashed-dotted line) and imaginary (dashed line) parts of the mode 
$\eta (y)$ for the dominant wavenumber 
$k_x=0.5$, normalized by the maximum modulus 
${\rm max}(|\eta |)$, for 
$u=1.5$. The mode is non-zero only in the steady central lobe region. For the same value of 
$u$, in figure 6(b) we report a three-dimensional plot of the linear combination of the base flow 
$H_r(y)$ (extended in the 
$x$ direction along which it is invariant) with the mode at the dominant wavenumber (normalized by the maximum modulus), i.e. 
$h(x,y)=H_r(y)+{A}\,\mathrm {Re}(\eta (y) \exp (\mathrm {i} k_x x))$, with 
${A}=0.25$ an arbitrary amplitude for visualization purposes. The resulting pattern is characterized by rivulets that carry lenses. The temporal dependence of the mode, which is not represented in figure 6(b), is characterized by a growing amplitude 
$\exp (\omega _i t)$ and by an oscillating behaviour 
$\exp (\mathrm {i} \omega _r t)$. The presence of a non-zero real part of 
$\omega$ (figure 5b) implies that the perturbations are oscillating in time at fixed locations. This effect is related to the advection as lenses are travelling along the streamwise direction.
Figure 6. Temporal stability analysis, 
$u=1.5$. (a) Real (solid line) and imaginary (dashed line) parts of the eigenvector 
$\eta (y)$, for the dominant wavenumber 
$k_x=0.5$, normalized by the maximum modulus. (b) Linear combination of the base flow 
$H_r(y)$ (extended in the 
$x$ direction along which it is invariant) with the mode at the dominant wavenumber (normalized with the maximum modulus), i.e. 
$h(x,y)=H_r(y)+{A}\,\mathrm {Re}(\eta (y) \exp (\mathrm {i} k_x x))$. 
$A=0.25$ is an arbitrary amplitude for visualization purposes.
The stability analysis reveals the occurrence of a secondary instability of the saturated and one-dimensional rivulets, which is located in the steady central lobe and leads to a pattern characterized by lenses that travel on the rivulets. Nevertheless, an increase in the advection 
$u$ induces a very strong stabilization and only very large wavelengths remain slightly unstable. The stabilization is related to the advection term. In particular, perturbations in regions of different thickness experience different advection velocities, proportional to 
$u H_r^2$ (Kalliadasis et al. Reference Kalliadasis, Ruyer-Quil, Scheid and Velarde2012). Regions of higher thickness travel faster than regions of lower thickness, leading to a steepening of the interface profile and eventually to a capillary levelling of perturbations. This steepening–levelling mechanism is all the more pronounced as 
$u$ is large. Small wavelengths, which present high interface gradients, are progressively stabilized with 
$u$, leading to a cutoff wavelength of the order of 
$10^2 {l}_c^*$ at 
$u=5$. In the numerical simulation of figure 3(a) the resulting pattern does not show any appreciable streamwise perturbations since the cutoff wavelength (
$L_c =2\pi /k_x^{cut} \approx 2 \times 10^2$) is of the order of the maximum acceptable wavelength fitting in the domain. These results are consistent with the experimental observations of Lerisson et al. (Reference Lerisson, Ledda, Balestra and Gallaire2020) when large values of 
$u$ are considered. For 
$u>3$, only very large wavelengths are unstable and they are eventually suppressed because of the size of the experiment (
$2 \times 10^2 l_c^* < L_x < 3 \times 10^2 l_c^*$).
4.4. Spatial stability analysis
In this section, we study the spatial stability properties of the rivulet base flow 
$H_r(y)$ introduced in § 4.1. The saturated rivulet profile is perturbed with a temporal harmonic perturbation of real frequency 
$\omega =\omega _r$ and we look for the spatial evolution of the perturbation, in terms of spatial growth rate 
$-\mathrm {Im}(k_x)$ and streamwise wavenumber 
$\mathrm {Re}(k_x)$ through the dispersion relation 
$D_r(k_x,\omega )$ ((4.3)). Positive values of the spatial growth rate denote unstable configurations associated with downstream propagating waves (Huerre & Rossi Reference Huerre and Rossi1998; Schmid et al. Reference Schmid, Henningson and Jankowski2002; Gallaire & Brun Reference Gallaire and Brun2017). The spectrum is characterized by only one unstable mode associated with downstream propagating waves, which is described in the following.
In figure 7(a) we report the spatial growth rate 
$-\mathrm {Im}(k_x)$ as a function of 
$\omega$. The spatial growth rate presents a behaviour similar to the temporal growth rate of § 4.3, i.e. characterized by a maximum (dominant) value and a cutoff frequency beyond which perturbations are damped. The dominant value of 
$-\mathrm {Im}(k_x)$ strongly decreases with 
$u$, while its associated dominant frequency presents a non-monotonic behaviour. The same non-monotonic behaviour is observed in the cutoff frequency. The streamwise wavenumber 
$\mathrm {Re}(k_x)$ (figure 7b) shows, to a good approximation, a linear dependence with 
$\omega$. For fixed 
$\omega$, the value of 
$\mathrm {Re}(k_x)$ decreases with 
$u$.
Figure 7. (a) Spatial growth rate and (b) streamwise wavenumber as functions of 
$\omega$, from the spatial stability analysis, for 
$u=1$ (blue line), 
$u=1.5$ (orange line), 
$u=2$ (yellow line), 
$u=2.5$ (purple line), 
$u=3$ (green line). The circles identify the values of the spatial growth rate obtained by the Gaster transformation.
The results of the spatial stability analysis are compared with those of the temporal stability analysis, suitably rescaled by the Gaster transformation (Gaster Reference Gaster1962), valid for strongly convectively unstable systems (see appendix B for details). Within this approximation, from the temporal stability analysis of § 4.3 (labelled with 
$(T)$) we retrieve the spatial stability analysis properties (labelled with 
$(S)$) through the relations
\begin{equation} \omega_r(S)=\omega_r(T), \quad \mathrm{Re}(k_x(S))=\mathrm{Re}(k_x(T)), \quad \mathrm{Im}(k_x(S))=-\frac{\omega_i(T)}{\dfrac{\partial \omega_r}{\partial k_x}(T)}. \end{equation}
The results of the Gaster transformation (4.4a–c) (circles) are in good agreement with the spatial stability analysis results in figure 7(a), for 
$u>1$. In appendix B we report the results for 
$u<1$, where the Gaster transformation prediction deviates from the spatial stability analysis results.
In the following, we experimentally investigate the link between the spatial stability analysis and the observable dynamics.
5. Experimental measurements of the rivulet secondary instability
5.1. Methods
As described in § 3.1, steady rivulets invade the experiment and saturate along the streamwise direction (figure 2). At a certain distance from the inlet, streamwise oscillations on the rivulet profiles grow and evolve into travelling lenses. We investigate the dependence of the overall dynamics and the amplitude of lenses with the parameters, by exploring different angles 
$40^{\circ }<\theta <80^{\circ }$ and thicknesses in the range 
$0.12<h_N/l_c< 1$ (related to the flow rate by (2.1a,b)). Note that 
$u= {\cot {(\theta )} {l}_c}/{\sqrt {\sin {\theta }}h_N}$, i.e. high values of the linear advection velocity correspond to low values of the flow rate or 
$\theta$. We modify the inlet condition using the spanwise comb-like blade (figure 1a) with the optimal spacing predicted by the flat film linear dispersion relation, i.e. 
$\hat {L}_r=2\pi \sqrt {2}l_c^*$. The requirement of a reasonably small and constant error in a large range of the flow parameters, exempted from a case-dependent calibration procedure, makes the confocal chromatic sensor STIL-CCS a suitable candidate. The latter is placed at the end of the plate to measure the variation of the amplitude of lenses oscillations as a function of 
$\theta$ and 
$h_N/l_c$.
The procedure is the following. We place the comb in position, and we wait the time necessary for rivulets to invade the whole domain. We then measure the central rivulet maximum thickness 
$\hat {h}_{max}(t)$ for 20 periods. This leads to a registration time 
$T$ that goes from 20 to 2000 s, depending on the angle. Once the data are registered, the flow rate is increased. We wait the time necessary to advect all the transient effects away from the glass plate; this varies from one minute, for 
$\theta =40^{\circ }$, to one hour, for 
$\theta =80^{\circ }$. Assuming the saturated rivulet profile (figure 4), we transform the point measurement of the maximum thickness in an estimate of the integral flux (i.e. 
$h_N$ in (2.1a,b)) by introducing the average thickness 
$\hat {\bar {h}}$ as follows:
\begin{equation} \hat{\bar{h}}=\left(\frac{3}{T} \int_0^T \frac{\hat{h}_{max}^3(t)}{3} \, \textrm{d} t\right)^{1/3} , \quad h_N=\hat{\bar{h}}/1.71, \end{equation}
being 
$\hat {h}_{max}=1.71h_N$ for a steady and saturated one-dimensional rivulet (Lerisson et al. Reference Lerisson, Ledda, Balestra and Gallaire2020). The deviation 
$\hat {\varDelta }$ from the average thickness value is computed as
\begin{equation} \hat{\varDelta}=\sqrt{\frac{1}{T}\int_0^T (\hat{h}_{max}(t)-\hat{\bar{h}} )^2 \,\mathrm{d}t}, \end{equation}
which is non-dimensionalized using the capillary length, i.e. 
$\varDelta =\hat {\varDelta }/l_c$. Two typical measurements are reported in figure 8.
Figure 8. Registered maximum height 
$\hat {h}_{max}(t)$ of the rivulet, for (a) 
$\theta =40^{\circ }$ and 
$h_N=1190\,\ \mathrm {\mu }$m , (b) 
$\theta =40^{\circ }$ and 
$h_N=1418\,\ \mathrm {\mu }$m. The black dashed line denotes 
$\hat {\bar {h}}$, and the red dashed lines 
$\hat {\bar {h}}\pm \hat {\varDelta }$. (a) 
$u=1.85$, 
$\hat {\varDelta }=0.004$ and (b) 
$u=1.57$, 
$\hat {\varDelta }=0.132$.
5.2. Results
In figure 9 we report the deviation 
$\varDelta$ as a function of 
$h_N/l_c$, for different angles 
$\theta$. At low values of 
$h_N/l_c$, 
$\varDelta$ is constant at a plateau value around 
${\varDelta } \sim 10^{-3}$. The plateau corresponds to the resolution of the optical sensor and is of order 
$1\,\ \mathrm {\mu }$m. At higher values of 
$h_N/l_c$, 
$\varDelta$ increases with 
$h_N/l_c$. We measure an increase of 
$\varDelta$ of two decades.
Figure 9. Values of 
$\varDelta$ (blue dots) as a function of 
$h_N/l_c$, for different values of 
$\theta$. The black horizontal line denotes the plateau value due to the resolution of the optical sensor. The red lines denote the amplification estimated using the spatial stability analysis of § 4.4 and the size of the plate, i.e. 
$\varDelta ={\varDelta }_0\exp (-\mathrm {Im}(k_x)L)$, with an initial amplitude chosen to obtain a good fit of the experimental data, (a) 
$\varDelta _0=2 \times 10^{-4}$, (b) 
$\varDelta _0=3 \times 10^{-5}$, (c) 
$\varDelta _0=7.5 \times 10^{-6}$, (d) 
$\varDelta _0=1 \times 10^{-6}$, (e) 
$\varDelta _0=1.5 \times 10^{-5}$, (f) 
$\varDelta _0=3 \times 10^{-6}$.
The amplitude of the oscillations at the end of the plate is compared with the theoretical findings of the spatial stability analysis. The spatial amplification at a distance 
$x$ of a temporal perturbation 
$\varDelta _0$ on a fully developed rivulet profile reads
\begin{equation} {\varDelta}/{\varDelta}_0=\exp(-\mathrm{Im}(k_x)x). \end{equation} We assume that the observable disturbances are the inlet ones as they are amplified on the largest distance, i.e. 
$x=L$. The perturbation amplitude 
$\varDelta _0$ originates from background noise that is below the sensitivity of our measurement sensor (
${\sim }1\ \mathrm {\mu }$m). We assume that the noise triggers the dominant mode described in §§ 4.3 and 4.4, and that 
$\varDelta _0$ is constant for a fixed angle. Note that the dominant spatial growth rate changes with 
$h_N/l_c$ since the value of 
$u$ is varied.
In figure 9 the red lines denote the theoretical values of 
${\varDelta }$ for an inlet perturbation amplitude 
${\varDelta }_0$ chosen to obtain a good fit of the experimental data. The measurement is then not a direct measure of the spatial growth rate, but of the variation of the spatial growth rate with the parameters. The variation of the deviation with the parameters well agrees with the linear prediction.
In figure 10 the experimental measurements of 
$\varDelta$ (coloured dots) are summarized and reported together with the spatial amplification 
${\varDelta }/{\varDelta }_0$ obtained by the spatial stability analysis (red dashed lines). At low values of 
$h_N/l_c$ the experimental values of 
$\varDelta$ are below the resolution of the optical sensor. As 
$h_N/l_c$ increases, 
$\varDelta$ emerges from the measurement resolution and we observe an increase of two orders of magnitudes in the considered range of parameters. This strong increase can be correlated to the theoretical amplification curves. At very low values of 
$h_N/l_c$ and inclination angles the theoretical amplification is of order 
${\varDelta }/{\varDelta }_0 \sim 10^{0}$. Low values of the flow rate (
$h_N/l_c$) or 
$\theta$ imply high values of 
$u$. In particular, the iso-level with value 
$\exp (-\mathrm {Im}(k_x)L)=1.3$ roughly corresponds to the case 
$u=3.5$. As 
$h_N/l_c$ and 
$\theta$ are increased the theoretical amplification rapidly grows.
Figure 10. Results of the analysis in the 
$(\theta ,h_N/l_c)$ plane: experimental measurements of 
$\varDelta$ (coloured dots) and inlet disturbance amplification 
$\varDelta /\varDelta _0=\exp (-\mathrm {Im}(k_x)L)$ evaluated by the spatial stability analysis of § 4.4 (red iso-contours).
Our analysis suggests that the occurrence of streamwise oscillations on the rivulet profile is strongly related to the advection. The measured deviations strongly vary with 
$u$. When high values of 
$u$ are considered, the occurrence of a steady and saturated rivulet pattern is observed (figure 11a). For low enough values of 
$u$, a state characterized by lenses which travel on rivulets is observed (figure 11b), as shown in (Lerisson et al. Reference Lerisson, Ledda, Balestra and Gallaire2019). Small variations in the advection lead to dramatic effects on the overall pattern dynamics. A change in the inclination of the plate of 10
$^\circ$, e.g. from 
$\theta =60^\circ$ to 
$\theta =50^\circ$ at 
$h_N/l_c=0.55$, is enough to pass from a state characterized by large amplitude lenses to a rivulet pattern.
Figure 11. Representative patterns at 
$\theta =45^{\circ }$ for (a) 
$h_N=623\,\ \mathrm {\mu }$m, i.e. 
$h_N/l_c=0.42$ and 
$u=2.83$, characterized by rivulets, and for (b) 
$h_N=1352\ \mathrm {\mu }$m, i.e. 
$h_N/l_c=0.92$ and 
$u=1.29$, characterized by rivulets which carry lenses.
In the route to dripping, the formation of lenses can be interpreted as a secondary instability of steady and streamwise-saturated rivulets, in which the role of the advection is essential.
6. Linear and nonlinear impulse response: breaking of isotropy and emergence of rivulets
In the previous sections, we numerically and experimentally studied the stability of steady and streamwise-saturated rivulet structures with respect to streamwise perturbations, and the link with the growth of travelling lenses. As observed in figure 2, the instability of rivulets and the consequent emergence of lenses is preceded by the formation of rivulet structures that invade the whole domain. Hereafter, we aim at giving a physical insight into the predominance of rivulet structures by studying the response of the flat film to an impulsive perturbation localized in space and time, i.e. the impulse response.
6.1. Experimental observation
In this section, we introduce a qualitative visualization of the evolution of a localized perturbation in the film thickness. The experimental apparatus is set without any inlet perturbation devices shown in figure 1(a). When high inclination angles and low flow rates are considered (i.e. high values of 
$u$), we experimentally observe a large region characterized by a uniform flat film where thickness perturbations from the lateral boundaries of the experiment do not penetrate (Lerisson et al. Reference Lerisson, Ledda, Balestra and Gallaire2020). In this region, we trigger the destabilization with a thickness perturbation by blowing a puff of air with a syringe. The whole field is then projected on a screen via the shadowgraph technique and captured with a camera.
In figure 12 we show the evolution of the perturbation with time. The initially localized perturbation is advected away in the streamwise direction with a constant velocity and spreads in the domain. The perturbation phase lines are concentric circles in the upstream part of the response. Nevertheless, the isotropy disappears in the downstream part. The shadowgraph reveals that the phase lines tend to be parallel to the streamwise direction, the effect becoming more and more evident as the time increases.
Figure 12. Shadowgraph visualization of an experimental impulse response, for 
$\theta =20^{\circ }$ and 
$h_N=1292\ \mathrm {\mu }$m, i.e. 
$u=5.45$. Time increases going to the right and each snapshot is separated by 
$15$ s. (a) 
$t=0$ s, (b) 
$t=15$ s, (c) 
$t=30$ s and (d) 
$t=45$ s.
The presence of phase lines aligned with the streamwise directions suggests the existence of a wavefront characterized by streamwise structures, i.e. rivulets, when the flat film is perturbed using an impulse thickness perturbation. The selection of a streamwise wavefront is not related to the boundaries of the thin film in the experiment, i.e. the rivulets selection is intrinsic.
6.2. Numerical observation
Inspired by this experimental observation, in this section we numerically simulate the impulse response, via (3.2), for the same values of angle and flow rate used in the shadowgraph of figure 12, i.e. 
$u=5.45$, in a double-periodic domain. The initial condition is taken in the form
\begin{equation} h(x,y,0)=1+{A}\exp\left(-\frac{x^2+y^2}{2}\right), \end{equation}
where 
$A=10^{-2}$. In figure 13 we plot the time evolution of the response in the moving reference frame 
$(\xi =x-u t,y)$, from 
$t=0$ to 
$t=90$. In the moving reference frame, the response progressively invades the domain from the initial impulse location. At 
$t=30$, we observe circular phase lines. At 
$t=60$ the response loses its isotropy in the downstream part. At 
$t=90$ streamwise structures are dominant in the downstream front of the response and they are also observable upstream.
Figure 13. Impulse response, 
$\theta =20^{\circ }$ and 
$h_N=1292\ \mathrm {\mu }$m (
$u=5.45$). The time increases from left to right and the time step is 
$30$. Results are reported in the moving reference frame at the linear advection velocity 
$(\xi =x-u t,y)$. (a) 
$t=0$, (b) 
$t=30$, (c) 
$t=60$, (d) 
$t=90$.
In the moving reference frame, the response spreads from the initial impulse location, meaning that in the fixed reference frame, the response is advected downstream at the linear advection velocity 
$u$. The numerical evolution qualitatively agrees with the experimental observation of § 6.1. We first observe the evolution of the impulse response into an isotropic pattern. At large times, the response mostly evolve towards streamwise structures. However, the complicated form of (3.2), including nonlinear advection, hydrostatic pressure distribution and capillary effects, does not allow one to identify the physical mechanisms that lead to the emergence of streamwise structures observed in figures 12 and 13. Lerisson et al. (Reference Lerisson, Ledda, Balestra and Gallaire2020) furthermore observed that the rivulet propagation and growth are well described by the linear stability analysis of the flat film even at large amplitudes of the thickness perturbation, beyond the expected validity of the linear theory. Hereafter, we study the origin of the selection of rivulet structures by the linear and weakly nonlinear dynamics.
6.3. Linear response
Upon introduction of the decomposition 
$h=1+\varepsilon \eta$ (
$\varepsilon \ll 1$) in (3.2), the linearized equation at 
$ {O}(\varepsilon )$ for the evolution of the thickness perturbation 
$\eta$ with respect to the flat film reads
\begin{equation} \partial_{t} \eta+u \partial_{x} \eta+\tfrac{1}{3}[\boldsymbol{\nabla}^{2} \eta+\boldsymbol{\nabla}^{4} \eta]=0.\end{equation}
The dispersion relation is recovered by introducing the normal mode decomposition 
$\eta \propto \exp [\mathrm {i}(\boldsymbol {k} \boldsymbol {\cdot } \boldsymbol {x}-\omega t)]$ , with 
$\boldsymbol {k}=(k_{x}, k_{y})$, where 
$k_x$ and 
$k_y$ denote respectively the streamwise and spanwise wavenumbers
\begin{equation} \omega=u k_{x}+\frac{\mathrm{i}}{3}({k}^{2}-{k}^{4}),\end{equation}
where 
$k=\sqrt {k_x^2+k_y^2}$. The dispersion relation 
$D(\omega ,k_x,k_y)=0$ is characterized by an isotropic temporal growth rate 
$\omega _i$, as shown in figure 14(a). The temporal frequency 
$\omega _r$ is linear in 
$k_x$ and does not depend on 
$k_y$.
Figure 14. (a) Temporal growth rate 
$\omega _i$ as a function of 
$k_x$ and 
$k_y$. (b) Linear impulse response, 
$u=5.45$. The vertical and horizontal axes are respectively the streamwise and spanwise directions. From left to right: 
$t=30$, 
$t=60$, 
$t=90$. Results are reported in the moving reference frame at the linear advection velocity 
$(\xi =x-u t,y)$.
The initial condition for the numerical simulation is the thickness perturbation 
$\eta (x,y,0)={{A}} \exp (-x^2/2 -y^2/2)$, where 
${A}=10^{-2}$. The linear numerical simulation results for 
$u=5.45$, in the moving reference frame 
$(\xi =x-u t,y)$, are presented in figure 14(b). As time increases, the perturbation spreads in concentric circles from the initial impulse location. Similarly to the nonlinear simulation of figure 13, the response is advected away at the linear advection velocity 
$u$, in the fixed reference frame. The results can be rationalized considering the dispersion relation of (6.2). The wavepacket is non-dispersive since 
$\omega _r$ is linear in 
$k_x$. This means that there is no distortion of the wavepacket. Since the growth is isotropic, concentric circles invade the domain and at the same time are advected downstream with constant velocity 
$\omega _r/k_x=u$. Higher values of 
$u$ imply faster advection velocities. In the moving reference frame 
$(\xi ,y)$, the (6.2) reads
\begin{equation} \partial_t \eta +\tfrac{1}{3}[ \boldsymbol{\nabla}_{\xi y}^2 \eta + \boldsymbol{\nabla}_{\xi y}^4 \eta ]=0,\end{equation}
where 
$\boldsymbol {\nabla }_{\xi y}$ operates in the reference frame 
$(\xi ,y)$. In this reference frame, the response spreads in perfectly isotropic concentric circles without being advected away. The linear dynamics agrees well with the early-time evolution of the nonlinear simulation shown in figure 13, when the amplitude of the perturbations is still very small. However, since the linearized dynamics is not able to capture the anisotropy of the pattern observed in the nonlinear simulation, we propose next a weakly nonlinear study.
6.4. Weakly nonlinear response: the Nepomnyashchy equation
We consider a weakly nonlinear model for the flow of a thin film on the underside of an inclined planar substrate. Following Kalliadasis et al. (Reference Kalliadasis, Ruyer-Quil, Scheid and Velarde2012), the derivation is based on a multiple scale approach combined with an asymptotic expansion. Under the assumption of small interfacial disturbances and 
$u= {O}(1)$, the weakly nonlinear dynamics for a thickness perturbation 
$\eta$ with respect to the flat film reads
\begin{equation} \partial_t \eta +2u \eta \partial_\xi \eta+\tfrac{1}{3} [ \boldsymbol{\nabla}_{\xi y}^2 \eta + \boldsymbol{\nabla}_{\xi y}^4 \eta ]=0, \end{equation}
where 
$\boldsymbol {\nabla }_{\xi y}$ operates in moving the reference frame 
$(\xi ,y)$. The equation is formally analogous to the Nepomnyashchy equation (Kalliadasis et al. Reference Kalliadasis, Ruyer-Quil, Scheid and Velarde2012). We consider the evolution of the thickness perturbation 
$\eta$ starting from a Gaussian impulse 
$\eta (\xi ,y,0)={A} \exp (-\xi ^2/2 -y^2/2)$ (
${A}=10^{-2}$), in analogy with the linear simulation.
In figure 15(a) we report the thickness perturbation evolution. The initial localized perturbation spreads in the domain and is always centred in the vicinity of the initial impulse position, because of the moving reference frame. At 
$t=30$ the perturbation has spread isotropically in the domain. Nevertheless, at 
$t=60$, streamwise structures arise. At 
$t=90$, the streamwise structures have invaded most of the perturbation region.
Figure 15. Case 
$u=5.45$. (a) Impulse response in the moving reference frame 
$(\xi =x-u t,y)$ from the weakly nonlinear model and (b) its two-dimensional Fourier energy spectrum. From left to right: 
$t=0$, 
$t=30$, 
$t=60$, 
$t=90$.
In figure 15(b) we show the two-dimensional Fourier energy spectrum of 
$\eta$, normalized by its maximum value. Since we are considering a real signal, the Fourier spectrum is symmetric with respect to the 
$k_x$ and 
$k_y$ axes. We thus report only the values in the first quadrant (
$k_x>0,k_y>0$). At 
$t=0$, we observe the Fourier spectrum of a Gaussian impulse, which is a Gaussian centred around (
$k_x=0,k_y=0$), i.e. the initial spectrum is isotropic. At 
$t=30$ the energy is located in a region around 
$\sqrt {k_x^2+k_y^2}=1/\sqrt {2}$. As time progressively increases, the energy concentrates towards 
$(k_x=0,k_y=1/\sqrt {2})$.
Initially, the response is characterized by an isotropic pattern, reminiscent of the linear growth that is experienced in the first stages of the perturbation growth. As the amplitude becomes sufficiently large, the spectrum shows that the energy is focusing on the axis 
$k_x=0$, i.e. streamwise structures are selected. The emergence of streamwise structures agrees well with the results of the fully nonlinear simulation and with the experimental observation. Moreover, the spectrum is localized around 
$k_y=1/\sqrt {2}$, the most amplified wavelength predicted by the flat film dispersion relation ((6.3)), and rivulet structures are growing exponentially. Thus, the dynamics of pure streamwise structures stays linear, even in the weakly nonlinear regime.
The origin of the selection of rivulet structures is identified in the weakly nonlinear advection term 
$2u \eta \partial _\xi \eta$, which acts in indirect manner to favour rivulet structures while damping all other orientations. The weakly nonlinear model of (6.5) is formally analogous to the linear model of (6.4), except for the weakly nonlinear advection term. It should be noticed that this term influences the dynamics of streamwise-inhomogeneous structures only, on which it is seen to have a damping effect. The nonlinear advection term embodies the difference in the perturbation advection velocity in regions of different thickness and is known to create wave steepening (Babchin et al. Reference Babchin, Frenkel, Levich and Sivashinsky1983). The emerging steep gradients are damped by surface tension effects, levelling therefore the non-streamwise structures. In conclusion, the most unstable solution in the weakly nonlinear regime is the one in which the capillary damping is reduced the most, as the term 
$2u \eta \partial _\xi \eta$, responsible of wave steepening, vanishes.
When only streamwise structures are present, the advection term disappears and the weakly nonlinear model is formally analogous to the linear equation in the moving reference frame (6.4). Consequently, the response of streamwise structures is linear up to second order in the perturbation.
In conclusion, the weakly nonlinear dynamics gives an insight into the origin of the emergence of rivulet structures: the latter are the only ones screened from the action of the difference in the advection. The dynamics of pure streamwise structures remains linear even in the weakly nonlinear regime, thus explaining the agreement between the linear prediction and the experimental measurements at large amplitudes observed in Lerisson et al. (Reference Lerisson, Ledda, Balestra and Gallaire2020). At late times, rivulets eventually invade the perturbation region. In the case of steady inlet forcing (figure 2) rivulets invade the whole domain and steady and streamwise saturated rivulet structures emerge downstream, as a result of the weakly nonlinear dynamics. As seen in the previous sections, rivulets may eventually destabilize through a secondary instability, resulting in travelling lenses. In both the emergence and the stability of rivulets, the differences in advection in regions of different thickness is crucial.
7. Conclusions
In this paper, we studied the selection and stability of rivulet structures in a thin film flowing under an inclined planar substrate. When the inlet is steadily forced along the spanwise direction, predominant rivulet structures were experimentally observed, which may destabilize at some distance from the inlet through the development of travelling lenses. Inspired by this experimental observation, we performed a nonlinear simulation with periodic boundary conditions, starting from an initial condition that mimicked the experimental forcing. The response to a streamwise-invariant sinusoidal initial condition confirmed the emergence of a persistent pattern of saturated rivulets, which may destabilize.
We then focused on the study of the mechanisms that may explain the behaviours observed in the experiment and numerical simulations, by studying the secondary stability of one-dimensional and saturated rivulets when perturbed in the streamwise direction. As the relative importance of advection increases, short wavelengths are progressively stabilized and only very large wavelengths remain slightly unstable. We relate their stabilization to the different advection of thickness perturbations on the rivulet profile. An increase in the advection results in steeper gradients for the same perturbation wavelength. Capillary forces counteract the wave steepening and eventually damp the perturbation, for high enough values of the advection. We compared the theoretical results for the spatial amplification of disturbances of the inlet flow rate with extensive experimental measurements of oscillations on rivulets, and confirmed the observation of a steady and saturated rivulet state when high values of 
$u$ are considered.
Finally, we gave an insight into the early-stage selection of streamwise-aligned structures, as observed in Lerisson et al. (Reference Lerisson, Ledda, Balestra and Gallaire2020), by studying the evolution of a localized impulse in the flat film. The experimental response showed that the wavefront selects mostly streamwise structures. The numerical impulse response also showed an initial isotropic growth followed by the selection of predominant rivulet structures. The numerical results were rationalized using a weakly nonlinear model, which showed the same selection of rivulets. The strength of the weakly nonlinear model was to identify one source of nonlinearity as the selection mechanism of streamwise structures, i.e. the weakly nonlinear advection. The latter is known to create wave steepening, counteracted by capillary terms. The evolution leads to levelling of all but streamwise structures. We concluded that the departure from a flat film towards streamwise structures is the solution in which the wave steepening and capillary damping effects are reduced the most. As a consequence, the selection of streamwise structures is due to the difference in the advection of perturbations in regions of different thickness, which acts to level all but pure streamwise perturbations (rivulets), while the dynamics of the latter remains linear even in the weakly nonlinear regime, thus rationalizing the results of Lerisson et al. (Reference Lerisson, Ledda, Balestra and Gallaire2020).
Our work aimed at laying rigorous foundations in the study of coating flows on the underside of planar substrates, interpreting the route to dripping as a destabilization of the flat film towards rivulets followed by a secondary instability. Nevertheless, several open questions are left. In complement to the spatio-temporal impulse response studied in this work, the response to a permanent in time but localized in space defect was considered briefly in Lerisson et al. (Reference Lerisson, Ledda, Balestra and Gallaire2020). However, a more detailed study to properly quantify the evolution of the response, e.g. in terms of asymptotic properties of the linear response, still needs to be performed.
Despite the predominance of streamwise-oriented structures, for some conditions, lenses appear on rivulets. While in this work a first analysis was performed in terms of spatial growth, a complete analysis of the precise evolution of perturbations along the streamwise direction remains to be pursued. In particular, a weakly non-parallel approach combined with a global resolvent technique could be suitable in this case. Furthermore, although the rivulet configuration shown in figure 11(b) may seem regular, we sometimes observe catastrophic events: lenses can merge in the streamwise direction, and eventually drip. While this work and the one of Lerisson et al. (Reference Lerisson, Ledda, Balestra and Gallaire2020) were focused on the emergence and stability of steady structures, further investigations focused on the dynamics of the travelling lenses are crucial to understand the route to dripping.
Acknowledgements
We acknowledge the anonymous referees for the valuable comments that helped to improve the manuscript. We acknowledge the Swiss National Science Foundation under grant 
$200021 \_ 178971$.
Declaration of interests
The authors declare no conflict of interest.
Appendix A. Expression of the Jacobian of the curvature
The operator 
$\partial _{\tilde {\eta }}\kappa (H)$, in general form, reads
\begin{align} \partial_{\tilde{\eta}}\kappa(H)&=\frac{(1+(\partial_{y}H)^2){\partial_{xx}} +(1+(\partial_{x}H)^2){\partial_{yy}}-2\partial_{x}H\partial_{y}H\partial_{xy}} {(1+(\partial_{x}H)^2+(\partial_{y}H)^2)^{3/2}} \nonumber\\ &\quad+\frac{2((\partial_{x}H)^2\partial_{y}H-\partial_{x}H(\partial_{y}H)^2)}{(1 +(\partial_{x}H)^2+(\partial_{y}H)^2)^{3/2}}{\partial_x} +\frac{2((\partial_{y}H)^2\partial_{x}H-\partial_{y}H(\partial_{x}H)^2)}{(1 +(\partial_{x}H)^2+(\partial_{y}H)^2)^{3/2}}{\partial_y}\nonumber\\ &\quad -3\frac{\partial_{xx}H(1+(\partial_{y}H)^2)+\partial_{yy}H(1+(\partial_{x}H)^2) -2\partial_{x}H\partial_{y}H\partial_{xy}H}{(1+(\partial_{x}H)^2 +(\partial_{y}H)^2)^{5/2}}(\partial_{x}H{\partial_x}+\partial_{y}H{\partial_y}). \end{align}The operator is evaluated for a base flow 
$H(x,y)=H_r(y)$ (i.e. 
$\partial _xH=0$) and we impose 
$\partial _x \tilde {\eta }=\mathrm {i}k_x \tilde {\eta }$.
Appendix B. Absolute-convective transition of the saturated rivulet profile
The purpose of this appendix is to verify the application of the Gaster transformation used in § 4.4. The Gaster transformation is applied in the context of strongly convectively unstable systems.
To verify the convective nature of the instability of the one-dimensional and steady rivulet profile, we evaluate the value of 
$u$ at which the absolute-convective transition occurs. We thus apply the Briggs–Bers criterion (Briggs Reference Briggs and Bers1964; Bers Reference Bers, DeWitt and Peyraud1975; Huerre & Monkewitz Reference Huerre and Monkewitz1990; Schmid et al. Reference Schmid, Henningson and Jankowski2002) to the dispersion relation 
$D_r(\omega ,k_x)=0$, (4.3). We look for the saddle points in the complex 
$k_x$ plane 
${\partial \omega }/{\partial k_x}=0$ and evaluate the imaginary part of 
$\omega$ at the saddle point 
$\mathrm {Im}(\omega _{0})$. The absolute-convective transition occurs when 
$\mathrm {Im}(\omega _{0})=0$. A spectral code is implemented in MATLAB, and saddle points are searched for with the built-in function fsolve. We identified a single saddle point in the complex-
$k_x$ plane. The absolute-convective transition occurs at 
$u_{0}=0.56$ (figure 16), which is much lower than the values of 
$u$ used throughout this work. Interestingly, the convective-absolute transition for the flat film takes place at 
$u_0=0.54$ (Brun et al. Reference Brun, Damiano, Rieu, Balestra and Gallaire2015), very close to the saturated rivulet value.
Figure 16. (a) Imaginary and (b) real parts of the complex frequency 
$\omega _0$ for the absolute-convective stability analysis. The absolute-convective transition occurs at 
$u_0=0.56$.
In figure 17 we report the comparison between the spatial stability analysis and the Gaster transformation, for 
$u<1$. As 
$u$ approaches the value for the absolute-convective transition, the prediction of the Gaster transformation deviates from the spatial stability analysis results.
Figure 17. Spatial growth rate given by the spatial stability analysis (solid lines) and by the Gaster transformation (circles), for 
$u=0.6$ (blue), 
$u=0.7$ (orange), 
$u=0.8$ (yellow), 
$u=0.9$ (purple).
