4.1. Region I

The main part of the bulk drop is in $x_{T}\leqslant x\leqslant x_{M}(t)~(x_{M}\gg 1)$ , where the surfactant concentration is almost uniform, ${\it\Gamma}={\it\Gamma}_{M}(t)$ . The dominant spreading mechanism for the fluid and surfactant is due to the horizontal component of gravity, with Marangoni forces contributing near $x=x_{T}$ and $x=x_{M}$ . We write ${\it\Gamma}={\it\Gamma}_{M}+\hat{{\it\Gamma}}(x,t)$ , taking $|\hat{{\it\Gamma}}(x,t)|\ll {\it\Gamma}_{M}$ . To leading order in $\hat{{\it\Gamma}}(x,t)/{\it\Gamma}_{M}$ , (2.2) becomes

(The dot notation used above represents the time evolution of ${\it\Gamma}_{M}$ and is used, henceforth, to denote the time evolution of any characteristic variable, unless stated otherwise.) Integrating (4.1b ), applying the boundary conditions $h=h_{T}$ , $\hat{{\it\Gamma}}_{x}=0$ at $x=x_{T}$ , gives the surfactant flux,

Substituting (4.2) in (4.1a ) gives the evolution equation of the drop height,

Equation (4.3) has solution

if the following compatibility conditions are satisfied (assuming $(x_{M}h_{M^{-}})_{t}\approx 0$ and ${\dot{x}}_{T}\approx 0$ , i.e. pinning the drop’s trailing edge, which are confirmed by the numerical solution):

where $h_{T}^{(0)}$ is an integration constant. Numerical solutions confirm that the above relationships are approximately satisfied at late times (see figure 6). This solution for $h$ (setting $x_{T}=h_{T}=0$ ) is the same as that derived for gravity-driven spreading of a drop (Huppert Reference Huppert1982; Troian et al. Reference Troian, Herbolzheimer, Safran and Joanny1989; Hocking Reference Hocking1990). On substituting (4.4) in (4.2) and using the relationships in (4.5), we obtain $\hat{{\it\Gamma}}_{x}=0$ , implying ${\it\Gamma}(x,t)={\it\Gamma}_{M}(t)$ in this region. The solutions for $h$ (rescaled by $h_{M^{-}}$ ) and ${\it\Gamma}$ (rescaled by ${\it\Gamma}_{M}$ ) versus $x-x_{T}$ (rescaled by $x_{M}$ ) in this region are shown in figure 7(a,b) (dashed lines). They match the numerical solution closely except that ${\it\Gamma}$ is non-uniform, particularly, near $x=x_{T}$ (where ${\it\Gamma}$ decreases as $x$ increases) and $x=x_{M}$ (where ${\it\Gamma}$ increases as $x$ increases). We first describe the solution behaviour near $x=x_{T}$ . Here, the above linearised analysis is not valid since changes in ${\it\Gamma}$ from ${\it\Gamma}_{M}$ are appreciable (see figure 7 b) and does not allow any meaningful approximation to be made. However, a qualitative understanding can be obtained as follows. Integrating (2.2b ) (neglecting capillary effects and the vertical component of gravity) and applying the boundary conditions $h=h_{T}$ , ${\it\Gamma}={\it\Gamma}_{T}>{\it\Gamma}_{M}$ , ${\it\Gamma}_{x}=0$ at $x=x_{T}$ gives the Marangoni surfactant flux,

The contribution from the unsteady term (first term on the right-hand side of (4.6)) is positive since ${\it\Gamma}_{t}<0$ for all $x$ . For $x$ near $x_{T}$ , this dominates the contribution from horizontal gravity (second term on the right-hand side of (4.6)), which is negative, resulting in a positive (or forward direction) Marangoni surfactant flux. For this scenario to occur, ${\it\Gamma}_{x}<0$ or ${\it\Gamma}$ decreases. As $x$ gets larger, the negative contribution from horizontal gravity increases (due to $h$ increasing), offsetting the positive contribution from the unsteady term, resulting in the Marangoni surfactant flux becoming zero and ${\it\Gamma}={\it\Gamma}_{M}$ , as described by the linearised analysis above. We note that, away from $x=x_{T}$ and $x=x_{M}$ , the cancellation between the unsteady and horizontal gravity terms predicted theoretically via the compatibility conditions in (4.5) for ${\it\Gamma}_{x}=0$ and ${\it\Gamma}={\it\Gamma}_{M}$ is not exactly achieved in the numerical results. This explains the slight deviation from ${\it\Gamma}_{M}$ observed in figure 7(b). We also note from (4.4) that the derivatives of $h$ become large as $x\rightarrow x_{T}$ , implying that the terms neglected here, particularly, the capillary term, need to be included to regularise the solution. This suggests the existence of an inner region near the trailing edge of the drop where gravity, Marangoni and surface tension forces compete. This is not discussed here. Next, we analyse region C to describe the evolution of $h$ and ${\it\Gamma}$ near $x=x_{M}$ .

4.2. Region C

Near $x=x_{M}$ , between the main bulk drop (region I) and the fluid hump region (region II), there is an abrupt positive jump in the drop height ( $h_{M^{-}}$ to $h_{M^{+}}$ ) and the surfactant concentration gradient ( $k_{M^{-}}$ to $k_{M^{+}}$ ) (see figure 2 at late times). This is due to a reverse (or backward direction) Marangoni flux slowing down the forward direction horizontal gravity flux without reversing the flow. This results in fluid accumulating towards this end and the formation of a shock-like structure. This structure develops over a length scale long enough for surface tension and diffusive effects to be negligible at leading order. Continuity of fluid flux across the region is ensured by the Rankine–Hugoniot condition

Continuity of surfactant flux across the region is ensured by imposing

In deriving these conditions, we assume that ${\it\Gamma}$ is continuous across the region and $k_{M^{-}}\approx 0$ as $x\rightarrow x_{M^{-}}$ . Integrating (2.2b ) (neglecting capillary effects and the vertical component of gravity) and applying the boundary conditions $h=h_{M^{-}}$ , ${\it\Gamma}={\it\Gamma}_{M}$ , ${\it\Gamma}_{x}=k_{M^{-}}=0$ at $x=x_{M^{-}}$ , gives the surfactant flux,

Substitution in (2.2a ) gives

We are unable to solve (4.9) and (4.10) analytically. However, (4.9) can be considered as a quadratic equation for $h$ , solving which we obtain

To make some analytical progress, we assume that, near $x=x_{M}$ , ${\it\Gamma}_{t}\approx \dot{{\it\Gamma}}_{M}$ . The surfactant flux (4.9) can be approximated by

Figure 8(a) compares the surfactant flux obtained from the numerical solution (solid lines) with the approximation in (4.12) (dashed lines) versus $x-x_{M}$ at $t=10^{5},5\times 10^{5},10^{6}$ . We observe that this linear approximation is accurate around $x=x_{M}$ . Further away from $x=x_{M}$ , this approximation breaks down, suggesting that ${\it\Gamma}_{t}\not\approx \dot{{\it\Gamma}}_{M}$ . This is considered in region II. Substituting ${\it\Gamma}_{t}\approx \dot{{\it\Gamma}}_{M}$ in (4.11) and using the compatibility condition (4.5b ), we can rewrite (4.11) as

Figure 8(b) compares the numerical solution for $h$ (solid lines) with its approximation in (4.13) (dashed lines) versus $x-x_{M}$ at $t=10^{5},5\times 10^{5},10^{6}$ . We use the numerical solution for ${\it\Gamma}$ and ${\it\Gamma}_{x}$ while computing the approximation for $h$ . We observe that (4.13) is a very good approximation for $h$ around $x=x_{M}$ . Further away from $x=x_{M}$ , the approximation, although overestimating $h$ , captures the gross behaviour, particularly the jump in $h$ , even though (4.12) breaks down here. This suggests that, although the flow and surfactant transport in this region are unsteady, it is the positive jump in ${\it\Gamma}_{x}$ that causes the positive jump in $h$ as shown in (4.13). It is worth noting that this region (and region II described below) has close analogies with the travelling wave structure described by Levy and co-workers (Levy & Shearer Reference Levy and Shearer2006; Levy et al. Reference Levy, Shearer and Witelski2007) in which the fluid layer height and surfactant concentration gradient are constant on either side of $x=x_{M}$ when the surfactant concentration has compact support. Assuming that (4.11) admits a travelling wave solution and $({\it\Gamma}_{M},k_{M^{-}})\approx 0$ , we have $h(x,t)=h_{M^{-}}$ for $x<x_{M}$ and $h(x,t)=h_{M^{+}}=k_{M^{+}}+\sqrt{k_{M^{+}}^{2}+2s}$ for $x>x_{M}$ , where $s$ denotes the travelling wave speed. Solving the latter equation for $k_{M^{+}}$ , we obtain $k_{M^{+}}=(h_{M^{+}}^{2}-2s)/(2h_{M^{+}})$ . Writing (4.10) in travelling wave coordinates and integrating, we obtain $-6sh+h^{3}=12C$ , where $C$ is an arbitrary constant. Assuming $C\approx 0$ and $s\approx h_{M^{+}}^{2}/3$ (which is the case when the precursor thickness ${\it\delta}\ll 1$ ), we obtain the non-zero solution $h_{M^{+}}=\sqrt{2}h_{M^{-}}$ . These relationships have been identified by Levy et al. (Reference Levy, Shearer and Witelski2007). The influence of the unsteady fluid flow and surfactant transport in (4.9) and (4.10) precludes any similar analytical relationships being derived for our problem.

Figure 8. Evolution of (a) surfactant flux $(h^{2}{\it\Gamma})/2-Mh{\it\Gamma}{\it\Gamma}_{x}$ and (b)  $h$ in region C using data shown in figure 2(a,b) at the times shown. Dashed lines show the approximate solution (4.12) and (4.13) (we use the numerical solution for ${\it\Gamma}$ and ${\it\Gamma}_{x}$ at the corresponding times while computing the approximation for $h$ in (4.13)).

4.3. Region II

This region is part of the bulk drop in $x_{M^{+}}(t)\leqslant x\leqslant x_{L}(t)$ where the fluid upwells due to competing reverse Marangoni and horizontal gravity fluxes resulting in fluid accumulation and a hump-like region being formed. There is a large variation in ${\it\Gamma}$ from ${\it\Gamma}_{M}$ to ${\it\Gamma}_{L}$ (see figure 2 b,d), $h$ varies between $h_{M^{+}}$ and $h_{L^{-}}$ , and ${\it\Gamma}_{x}$ varies between $k_{M^{+}}$ and $k_{L^{-}}$ . Using the continuity of surfactant flux at $x=x_{M}$ given by (4.8), the evolution equations in this region are given by (4.9)–(4.11). We use a linear approximation for

Substituting (4.14) in (4.9) and (4.11) gives

Figure 9(a) compares the surfactant flux obtained from the numerical solution (solid lines) with the approximation in (4.15a ) (dashed lines) versus $x-x_{M}$ at $t=10^{5}$ , $5\times 10^{5},10^{6}$ . We observe that (4.15a ) is a much better approximation of the surfactant flux compared to the linear approximation in (4.12). Figure 9(b) compares the numerical solution for $h$ (solid lines) with its approximation in (4.15b ) (dashed lines) versus $x-x_{M}$ at $t=10^{5},5\times 10^{5},10^{6}$ . We use the numerical solution for ${\it\Gamma}$ and ${\it\Gamma}_{x}$ while computing the approximation for $h$ . Although this approximation underestimates the numerical solution, it captures the gross behaviour of the evolution of $h$ in this region. We again note the analogy with the work by Levy and co-workers (Levy & Shearer Reference Levy and Shearer2006; Levy et al. Reference Levy, Shearer and Witelski2007) in which the fluid layer height and surfactant concentration gradient are constant in this region. We observe from figure 9(b) that the fluid height is not constant even at very late time, suggesting the long-lived unsteady influence of the fluid flow and surfactant transport in (4.14) and (4.15). Evaluating (4.15b ) at $x=x_{L}$ gives

Figure 9. Evolution of (a) surfactant flux $(h^{2}{\it\Gamma})/2-Mh{\it\Gamma}{\it\Gamma}_{x}$ and (b)  $h$ in region II using data shown in figure 2(a,b) at the times shown. Dashed lines show the approximate solution (4.15a,b ) (we use the numerical solution for ${\it\Gamma}$ and ${\it\Gamma}_{x}$ at the corresponding times while computing the approximation for $h$ in (4.15b )).

Figure 10. Evolution of (a)  $h$ and (b)  ${\it\Gamma}$ in region A using data shown in figure 2(a,b). Dashed lines show the solution of (4.20) for $\hat{V}=0.18485$ , ${\hat{h}}_{L^{+}}=0.1$ , $\hat{k}_{L^{+}}=-1.797$ and $\hat{k}_{L^{-}}=0.3153$ (this corresponds to a particular solution curve in the family of solutions).

4.4. Region A

A thin transition region at $x=x_{L}(t)$ of height $O(h_{L^{-}})$ connects the fluid hump (region II) to the fluid front (region III). This effective contact-line region is controlled by competing surface tension, horizontal gravity and Marangoni forces. Balancing surface tension and horizontal gravity fluxes gives the width of the region to be $O((\mathit{Ca}\,h_{L^{-}})^{1/3})$ . Balancing convective and horizontal gravity fluxes, ${\dot{x}}_{L}h_{x}\sim (h^{3})_{x}$ , implies ${\dot{x}}_{L}\sim h_{L^{-}}^{2}$ . Balancing horizontal gravity and Marangoni fluxes implies that changes in ${\it\Gamma}$ across this region are of size $(\mathit{Ca}\,h_{L^{-}}^{4})^{1/3}/M$ and are assumed small compared to ${\it\Gamma}_{L}$ . We then set

Figure 10(a,b) shows a family of quasi-steady solutions when the PDE data are scaled using these variables. These solutions are parametrised by $\hat{V}$ (the effective contact-line speed, representing the ratio of the advective to horizontal gravity fluxes), ${\hat{h}}_{L^{+}}$ (the ratio of the downstream film thickness to that upstream), and $\hat{k}_{L^{+}}$ and $\hat{k}_{L^{-}}$ (representing the ratio of the Marangoni to horizontal gravity fluxes downstream and upstream of the contact line, respectively). Assuming ${\it\xi}$ , $H$ , $G$ , $\hat{V}$ , ${\hat{h}}_{L^{+}}$ , $\hat{k}_{L^{+}}$ and $\hat{k}_{L^{-}}$ are all $O(1)$ , (2.2) reduces at leading order to subject to the boundary (or matching) conditions where $\hat{D}=(Dh_{L^{-}}^{2/3})/\mathit{Ca}^{1/3}$ . Equation (4.19) assumes that $h$ is constant far upstream and downstream, where it matches onto regions II and III, respectively. This is reasonable considering that $h$ in both these regions is only very weakly linear (see figure 10 a). The assumption of constant $h$ matching onto region III is valid only for sufficiently small values of precursor film thickness ${\it\delta}$ (see description of region III below). Also ${\it\Gamma}$ is assumed linear, matching onto regions II and III with slope $k_{L^{-}}~({>}0)$ and $k_{L^{+}}~({<}0)$ , respectively. While this matching condition is approximate for the fluid hump region, it is exact when matching onto the fluid front ahead, since ${\it\Gamma}$ is linear there (see region III description below). Integrating (4.18) and applying the boundary conditions (4.19), we obtain In addition, continuity of fluid and surfactant flux are ensured across this region by the following conditions (in original variables): Equation (4.20a ) is the same as the evolution of $h$ in the case of only gravity-driven spreading (Troian et al. Reference Troian, Herbolzheimer, Safran and Joanny1989; Hocking Reference Hocking1990; Bertozzi & Brenner Reference Bertozzi and Brenner1997). The effective contact-line speed ${\dot{x}}_{L}$ in (4.21a ) shows the competition between gravity and the Marangoni effect, with Marangoni slowing down the contact-line speed. The numerical solution of (4.20) subject to the boundary conditions in (4.19) is shown in figure 10(a,b) (dashed lines) for $\hat{V}=0.18485$ , ${\hat{h}}_{L^{+}}=0.1$ , $\hat{k}_{L^{+}}=-1.797$ and $\hat{k}_{L^{-}}=0.3153$ (the constant arising from integrating (4.20b ) is chosen such that $G(0)=0$ , since by definition ${\it\xi}=0$ is where $G$ has a maximum). It matches the numerical solution closely and corresponds to a particular solution curve in the family of solutions. Equation (4.20) is similar to the differential equations derived by Levy et al. (Reference Levy, Shearer and Witelski2007) for the fluid layer height and surfactant concentration assuming a travelling wave solution structure in this region. Moreover, their equations are also coupled, unlike (4.20). We assume ${\it\Gamma}_{t}$ to be small in this region compared to ${\it\Gamma}_{L}$ , which allows the equations to be decoupled. Retaining this unnecessarily complicates the equations without significantly influencing the solution structure.

Figure 11. Evolution of (a)  $h$ and (b)  ${\it\Gamma}$ in region III using data shown in figure 2(a,b). Dashed lines show the solution in (4.24a,b ).

4.5. Region III

This region contains the spreading fluid front in $x_{L}<x<x_{SL}$ of height $O({\it\delta})$ . The flow in this region is predominantly driven by Marangoni forces and exhibits the self-similar structure similar to that identified in Jensen & Grotberg (Reference Jensen and Grotberg1992). Balancing convective and Marangoni fluxes, ${\dot{x}}_{SL}h_{x}\sim (Mh^{2}{\it\Gamma}_{x})_{x}$ , implies ${\it\Gamma}_{x}\sim {\dot{x}}_{SL}/(M{\it\delta})$ . Hence, ${\it\Gamma}\sim {\dot{x}}_{SL}(x_{SL}-x_{L})/(M{\it\delta})$ ( $x_{SL}>x_{L}$ for all $t$ ). We set

Figure 11(a,b) shows the numerical solution in this region to collapse when the PDE data are scaled using these variables. With error $O(1/(x_{SL}-x_{L}))$ , (2.2) becomes subject to $H=2$ and $G=0$ at ${\it\xi}=1$ (see Jensen & Grotberg Reference Jensen and Grotberg1992). Assuming ${\it\gamma}=1+({\it\alpha}{\it\alpha}_{1})/2$ (figure 12 a confirms this; noise in the data is due to numerical approximation in computing the derivatives involved in the relevant variables) and integrating (4.23), we obtain

where $A=1+2/({\it\alpha}{\it\alpha}_{1})$ . The time evolution of $A$ seen in figure 12(b) shows it to be large and negative, with $A\approx -50$ (the noise in the data is attributed to the numerical approximation in computing the derivatives involved). The dashed lines in figure 11(a,b) show $H$ and $G$ given by (4.24), which match the rescaled numerical solution very well. We note that, since $A$ is large and negative, so the slope of $h$ is small and positive, hence, $h\approx 2{\it\delta}$ in this region, which makes it appear to look like a ‘step’. This is true for sufficiently small values of the precursor film thickness ${\it\delta}$ . However, for larger values of ${\it\delta}$ , $A$ decreases in magnitude and hence the slope of this region increases and it no longer looks like a step. Hence, for sufficiently small ${\it\delta}$ and $x\rightarrow x_{L}^{+}$ , we obtain

It is worthwhile comparing (4.24) with the travelling wave solution derived by Levy & Shearer (Reference Levy and Shearer2006), where the fluid layer height and surfactant concentration gradient in this region are constant. In the travelling wave frame of reference, a balance between the convective, Marangoni and horizontal gravity fluxes admits constant solutions in fluid layer height and surfactant concentration gradient if and only if the travelling wave speed $s=-h_{L^{+}}k_{L^{+}}+h_{L^{+}}^{2}/2$ (the surfactant transport speed). However, in our problem the dynamics of this region is controlled by convective and Marangoni fluxes, which dominate that due to horizontal gravity. This balance admits a linear solution in the fluid height (and surfactant concentration), which distinguishes it from that derived by Levy & Shearer (Reference Levy and Shearer2006).

Figure 12. (a) Numerical validation of ${\it\gamma}=1+({\it\alpha}{\it\alpha}_{1})/2$ and (b) the time evolution of $A=1+2/({\it\alpha}{\it\alpha}_{1})$ , showing it to be large and negative and $A\approx -50$ . The noisy data are due to numerical approximation in computing the derivatives involved in the relevant variables.

4.6. Region B

The structure of the kinematic shock near $x=x_{SL}$ has been described in detail previously (Jensen & Grotberg Reference Jensen and Grotberg1992; Jensen Reference Jensen1994; Jensen & Halpern Reference Jensen and Halpern1998) and is not discussed here. In the parameter regime relevant here, the discontinuity in the film thickness (equal to $2{\it\delta}$ ) is smoothed by surface tension effects and the jump in the surfactant gradient (equal to $k_{L^{+}}$ ) is smoothed by surface diffusion (see figure 12 a,b near ${\it\xi}=1$ ), but the region remains dynamically passive. Ahead of region C, the film in region IV is undisturbed and equals the precursor thickness.

4.7. DAE model

Equations (4.5b,c ), (4.7), (4.8), (4.16), (4.21a,b ) and (4.25a−c ) provide a system of differential–algebraic equations for the characteristic variables ${\it\Gamma}_{M}$ , $h_{T}$ , $x_{M}$ , $k_{M^{+}}$ , $h_{L^{-}}$ , $x_{L}$ , $k_{L^{-}}$ , $h_{L^{+}}$ , ${\it\Gamma}_{L}$ and $k_{L^{+}}$ , respectively. The total fluid volume and surfactant mass in the region $x_{T}\leqslant x<\infty$ are conserved and can be approximated as

where $\mathscr{V}$ and $\mathscr{M}$ are the total fluid volume and surfactant mass, respectively. Since we do not have analytical expressions for $h$ and ${\it\Gamma}$ in region II, we approximate this region as a trapezium when evaluating the fluid volume and surfactant mass. Equation (4.26) provide equations for the characteristic variables $h_{M^{+}}$ and $x_{SL}$ , respectively. The equation for the remaining characteristic variable $h_{M_{-}}$ is derived as follows. Since $h$ in region I is similar to that due to spreading under gravity alone with no surfactant, we can use the solution for $h$ in (4.4) (with $h_{T}$ , $x_{T}\approx 0$ ) and similarity to write $h_{M_{-}}/{\hat{h}}_{L}=\sqrt{x_{M}/\hat{x}_{L}}$ , where $\hat{x}_{L}$ and ${\hat{h}}_{L}$ are the location and height of the leading edge of the drop spreading under gravity alone. Equation (4.26a ) for gravity-driven spreading with no surfactant reduces to $\mathscr{V}=2{\hat{h}}_{L}\hat{x}_{L}/3$ . The effective contact-line speed (4.21a ) for this case reduces to $\dot{\hat{x}}_{L}={\hat{h}}_{L}^{2}/3$ (setting $M=0$ and $h_{L^{+}}={\it\delta}\ll {\hat{h}}_{L}$ ). Solving these, we obtain $\hat{x}_{L}=[(9/4)\mathscr{V}^{2}(t-t_{0})+\hat{x}_{L}^{0^{3}}]^{1/3}$ and ${\hat{h}}_{L}=3\mathscr{V}/\{2[(9/4)\mathscr{V}^{2}(t-t_{0})+\hat{x}_{L}^{0^{3}}]^{1/3}\}$ , where $\hat{x}_{L}^{0}$ is the value of $\hat{x}_{L}$ at $t=t_{0}$ . Hence,

The system of DAEs are parametrised by $\mathscr{V}$ , $\mathscr{M}$ , $M$ and ${\it\delta}$ . The PDE simulations used $\mathscr{M}=4/3$ and $\mathscr{V}=2$ . Initial values for the differential equations, ${\it\Gamma}_{M}^{0}$ , $h_{T}^{0}$ , $x_{M}^{0}$ , $x_{L}^{0}$ , $x_{SL}^{0}$ and $\hat{x}_{L}^{0}$ , are fitted from the numerical solution at $t_{0}=10^{5}$ . These are used to determine consistent initial conditions for the algebraic equations. The DAEs are solved numerically, and their solution provides a good approximation for the behaviour in $t>10^{5}$ (dashed lines in figure 5 a–d). The unsteady nature of the flow and surfactant transport in regions C and II prohibits a uniformly asymptotic approximation and limits the accuracy of the predictions of $h_{M^{+}}$ , $h_{L^{-}}$ , $k_{M^{+}}$ and $k_{L^{-}}$ , but nevertheless illustrates an important feature of the spreading dynamics. Once $x_{M}\gg x_{T}$ , further simplification is possible: $h_{M_{-}}$ and $x_{M}$ are similar to gravity-driven spreading, so $h_{M_{-}}\propto {\hat{h}}_{L}\propto t^{-1/3}$ and $x_{M}\propto \hat{x}_{L}\propto t^{1/3}$ (Hocking Reference Hocking1990). The majority of the surfactant is contained in regions II and III and ${\it\Gamma}_{M}\ll {\it\Gamma}_{L}$ , so (4.26b ) can be approximated as $\mathscr{M}\approx {\it\Gamma}_{L}(x_{L}-x_{M})/2+{\dot{x}}_{SL}(x_{SL}-x_{L})^{2}/(4{\it\delta}M)$ . This implies ${\dot{x}}_{SL}(x_{SL}-x_{L})^{2}\sim 1$ and ${\it\Gamma}_{L}(x_{L}-x_{M})\sim 1$ , which suggests $x_{SL}\propto x_{L}\propto t^{1/3}$ and ${\it\Gamma}_{L}\propto t^{-1/3}$ (consistent with (4.25b )). Equation (4.25c ) gives $k_{L^{+}}\propto {\dot{x}}_{SL}\propto t^{-2/3}$ . A balance of terms in (4.21b ) then implies that $h_{L_{-}}\propto k_{L^{+}}^{1/2}\propto t^{-1/3}$ and $k_{L^{-}}h_{L_{-}}\propto k_{L^{+}}$ , hence $k_{L^{-}}\propto t^{-1/3}$ . A balance of terms in (4.8) gives $h_{M_{+}}\propto h_{M_{-}}\propto t^{-1/3}$ and $k_{M^{+}}h_{M_{+}}\propto h_{M_{-}}^{2}$ , hence $k_{M^{+}}\propto t^{-1/3}$ . Finally, (4.5c ) gives $h_{T}\propto t^{-1/2}$ . These scalings are all borne out by the numerical solutions in figure 5(a–d).

5.1. Small-wavenumber analysis

In this section we examine instabilities with wavelength intermediate between the width of region A and regions II and III, and assume exponential time dependence so that ${\hat{H}}_{{\it\tau}}$ in (5.2a ) becomes ${\it\beta}{\hat{H}}$ . Numerical dispersion relations (see figure 14) suggest that the growth rate is linear at small wavenumbers, so we write $[{\hat{H}},{\hat{G}},{\it\beta}]=[H,G,{\it\beta}]_{0}+q[H,G,{\it\beta}]_{1}+q^{2}[H,G,{\it\beta}]_{2}+\cdots \,$ . Substituting in (5.2), at leading order in $q$ , we obtain $(H_{0},G_{0})=(-H_{s{\it\xi}},-G_{s{\it\xi}})$ and ${\it\beta}_{0}=0$ , which represents a translation of the base state $(H_{s},G_{s})$ . At $O(q)$ , we obtain

We assume that $H_{1}$ and all its derivatives decay to zero as ${\it\xi}\rightarrow \pm \infty$ and $G_{1{\it\xi}}$ tends to a constant $A$ as ${\it\xi}\rightarrow \infty$ and $G_{1{\it\xi}}$ tends to a constant $B$ as ${\it\xi}\rightarrow -\infty$ . The last two conditions are required to balance the unsteady and the Marangoni terms (first and third terms) in (5.3a ) and are also motivated in part by the fact that the leading-order solution $G_{0}$ tends to a constant value $-\hat{k}_{L^{-}}$ and $-\hat{k}_{L^{+}}$ as ${\it\xi}\rightarrow -\infty$ , $\infty$ , respectively. Using $H_{s}\rightarrow {\hat{h}}_{L^{+}}$ , $G_{s{\it\xi}}\rightarrow \hat{k}_{L^{+}}$ , $H_{s{\it\xi}},H_{s{\it\xi}{\it\xi}{\it\xi}}\rightarrow 0$ as ${\it\xi}\rightarrow \infty$ and $H_{s}\rightarrow 1$ , $G_{s{\it\xi}}\rightarrow \hat{k}_{L^{-}}$ , $H_{s{\it\xi}},H_{s{\it\xi}{\it\xi}{\it\xi}}\rightarrow 0$ as ${\it\xi}\rightarrow -\infty$ , we obtain the following conditions after integrating (5.3) between $-\infty$ and $\infty$ :

To determine the growth rate ${\it\beta}_{1}$ , we need to determine either $A$ or $B$ , which is done as follows. When $q\ll 1$ , we assume that ${\hat{H}}$ and ${\hat{G}}$ have a three-region structure: an inner region near the drop’s contact line in which ${\it\xi}=O(1)$ and ${\hat{H}}=-H_{s{\it\xi}}$ , ${\hat{G}}=-G_{s{\it\xi}}$ to leading order, and two outer regions of length scales $O(1/q)$ . We suppose that the drop’s contact line (when viewed from the outer regions) lies along ${\it\xi}={\it\epsilon}\exp (\text{i}q{\it\eta}+{\it\beta}{\it\tau})$ , so that to leading order, at the outer limits of the inner region, we have

We now show how perturbations to $G$ affect the surfactant gradient ${\hat{G}}_{{\it\xi}}$ ahead of and behind the contact line.

To describe the outer regions, we rescale with ${\it\xi}=z/q$ , ${\hat{G}}=\tilde{G}/q$ , ${\it\beta}=q\tilde{{\it\beta}}$ (anticipating $\tilde{{\it\beta}}=O(1)$ ), so that (5.2) becomes (after dividing by $q$ )

At leading order in $q$ as $z\rightarrow \infty$ , we can write (5.6) as which reduces to

where ${\it\alpha}=[({\hat{h}}_{L^{+}}^{2}/2)+({\hat{h}}_{L^{+}}\hat{k}_{L^{+}}/2)-{\hat{h}}_{L^{+}}^{2}+\hat{V}]$ . Thus ${\hat{H}}={\hat{H}}_{0}\exp (\tilde{{\it\beta}}z/{\it\alpha})$ for some ${\hat{H}}_{0}$ . For $\tilde{{\it\beta}}>0$ and ${\it\alpha}>0$ (observed numerically), we require $\tilde{H}_{0}=0$ for perturbations to remain bounded as $z\rightarrow \infty$ (the contact line moves faster than any growing disturbances to the fluid front ahead). Thus $\tilde{G}=\tilde{G}_{0}\exp (-z)$ as $z\rightarrow \infty$ for some $\tilde{G}_{0}$ . To match with (5.5b ) we take $\tilde{G}_{0}=-q\hat{k}_{L^{+}}$ . Thus towards the contact line as $z\rightarrow 0^{+}$ , ${\hat{H}}=0$ and

Hence, as ${\it\xi}\rightarrow \infty$ , $G_{1{\it\xi}}\rightarrow A=\hat{k}_{L^{+}}$ . Similarly, (5.6) at leading order in $q$ as $z\rightarrow -\infty$ can be written as

which reduces to

where ${\it\alpha}_{1}=[(1/2)+(\hat{k}_{L^{-}}/2)-1+\hat{V}]$ . Thus ${\hat{H}}={\hat{H}}_{0}\exp (\tilde{{\it\beta}}z/{\it\alpha}_{1})$ for some ${\hat{H}}_{0}$ . For $\tilde{{\it\beta}}>0$ and ${\it\alpha}_{1}<0$ (observed numerically), we require $\tilde{H}_{0}=0$ for perturbations to remain bounded as $z\rightarrow -\infty$ . Thus $\tilde{G}=\tilde{G}_{0}\exp (z)$ as $z\rightarrow -\infty$ for some $\tilde{G}_{0}$ . To match with (5.5a ) we take $\tilde{G}_{0}=-q\hat{k}_{L^{-}}$ . Thus towards the contact line as $z\rightarrow 0^{-}$ , ${\hat{H}}=0$ and

Hence, as ${\it\xi}\rightarrow -\infty$ , $G_{1{\it\xi}}\rightarrow B=-\hat{k}_{L^{-}}$ . Using (5.4a ) we obtain the $O(q)$ growth rate,

Figure 15. Numerical validation of the relationships (a)  $\hat{k}_{L^{-}}=-\hat{k}_{L^{+}}{\hat{h}}_{L^{+}}$ , (b)  $\hat{k}_{L^{-}}=(1-{\hat{h}}_{L^{+}}^{2})/4$ and (c)  $\hat{k}_{L^{+}}=-(1-{\hat{h}}_{L^{+}}^{2})/(4{\hat{h}}_{L^{+}})$ .

Noting that $\hat{k}_{L^{-}}>0$ , $\hat{k}_{L^{+}}<0$ and ${\hat{h}}_{L^{+}}<1$ , the growth rate ${\it\beta}_{1}>0$ if $\hat{k}_{L^{-}}/2>-\hat{k}_{L^{+}}{\hat{h}}_{L^{+}}^{2}/2$ , i.e. the $O(q)$ Marangoni flux behind the contact line dominates that ahead of the contact line. If this condition is satisfied, a mechanism for the flow to become linearly unstable at $O(q)$ is due to the forward Marangoni flux behind the contact line pulling more fluid into the contact-line region than the forward Marangoni flux ahead of the contact line dragging fluid out. This results in growth of perturbations in $h$ , consequently destabilising the flow and the contact line (consistent with figure 13, where ${\hat{H}}$ is positive and growing immediately behind the contact line, coinciding with ${\hat{G}}_{{\it\xi}}$ , which is negative there). We now verify whether the above condition is satisfied. Using (5.4b ) implies $\hat{k}_{L^{-}}=-\hat{k}_{L^{+}}{\hat{h}}_{L^{+}}$ (figure 15 a shows that this relationship is approximately satisfied; the slight difference is due to the ambiguity in determining $k_{L^{-}}$ from the numerical solution), using which the above condition can be rewritten as $\hat{k}_{L^{-}}(1-{\hat{h}}_{L^{+}}^{2})>0$ , which is always satisfied. We also know using (4.21b ) that $\hat{k}_{L^{-}}=(1-{\hat{h}}_{L^{+}}^{2})/2+{\hat{h}}_{L^{+}}\hat{k}_{L^{+}}$ , which, when combined with the above relationship, gives $\hat{k}_{L^{-}}=(1-{\hat{h}}_{L^{+}}^{2})/4$ and $\hat{k}_{L^{+}}=-(1-{\hat{h}}_{L^{+}}^{2})/(4{\hat{h}}_{L^{+}})$ (figure 15 b,c shows that these relationships are approximately satisfied; the slight difference in figure 15 b is due to the ambiguity in determining $k_{L^{-}}$ from the numerical solution). Equation (5.13) can now be rewritten as

Equation (5.14) shows that the $O(q)$ growth rate is always positive and independent of the angle of inclination. This suggests that the flow is linearly unstable for all angles of inclination, at least those within the validity of our model. This supports the speculation by Edmonstone et al. (Reference Edmonstone, Matar and Craster2005b ). We also note from numerics that ${\hat{h}}_{L^{+}}$ increases gradually with time, so ${\it\beta}_{1}$ decreases, confirming the observation in figure 14(a) that, as time increases, the base state becomes less linearly unstable. The analysis also shows the Marangoni effect to be dominant at this order, with gravity having no influence. Figure 16(a,b) shows the $O(q)$ approximation ${\it\beta}\approx q/8$ (dashed lines) along with the numerical dispersion relation (solid lines) for inclination angles ${\it\theta}=90^{\circ }$ and 3°, respectively. We observe that this approximation slightly underestimates the numerical dispersion relation for ${\it\theta}=90^{\circ }$ and overestimates it for ${\it\theta}=3^{\circ }$ .

We now determine the $O(q^{2})$ growth rate ${\it\beta}_{2}$ . At $O(q^{2})$ we obtain

We assume that $H_{1}$ , $H_{2}$ and their derivatives decay to zero as ${\it\xi}\rightarrow \pm \infty$ . Equations (5.9) and (5.12) suggest that $G_{2{\it\xi}}\rightarrow -\hat{k}_{L^{+}}{\it\xi}$ as ${\it\xi}\rightarrow \infty$ and $G_{2{\it\xi}}\rightarrow -\hat{k}_{L^{-}}{\it\xi}$ as ${\it\xi}\rightarrow -\infty$ , respectively. Using $H_{s}\rightarrow {\hat{h}}_{L^{+}}$ , $G_{s{\it\xi}}\rightarrow \hat{k}_{L^{+}}$ , $H_{s{\it\xi}},H_{s{\it\xi}{\it\xi}{\it\xi}}\rightarrow 0$ as ${\it\xi}\rightarrow \infty$ , $H_{s}\rightarrow 1$ , $G_{s{\it\xi}}\rightarrow \hat{k}_{L^{-}}$ , $H_{s{\it\xi}},H_{s{\it\xi}{\it\xi}{\it\xi}}\rightarrow 0$ as ${\it\xi}\rightarrow -\infty$ , $H_{0}=-H_{s,{\it\xi}}$ , $G_{0}=-G_{s{\it\xi}}$ and ${\it\beta}_{0}=0$ , we obtain the following condition after integrating (5.15a ) between $-\infty$ and $\infty$ :

Figure 16. Comparison of late-time numerical dispersion relation (solid lines) with the small-wavenumber approximation, ${\it\beta}={\it\beta}_{1}q$ (dashed lines) and ${\it\beta}={\it\beta}_{1}q+{\it\beta}_{2}q^{2}$ (dot-dashed lines) for inclination angles (a)  ${\it\theta}=90^{\circ }$ and (b)  ${\it\theta}=3^{\circ }$ . Values ${\it\beta}_{1}=1/8$ and ${\it\beta}_{2}=0.0207$ are used for ${\it\theta}=90^{\circ }$ , and ${\it\beta}_{2}=-0.0847$ for ${\it\theta}=3^{\circ }$ .

Equation (5.16) shows that the $O(q^{2})$ growth rate ${\it\beta}_{2}$ is controlled by competing Marangoni, capillary and vertical gravity fluxes. The first two terms on the right-hand side of (5.16) cancel out the divergent contribution to the Marangoni flux integral (third term on the right) far ahead of and behind the contact line, respectively. The capillary and vertical gravity flux integrals (fourth and fifth terms on the right, respectively) are convergent. For spreading without surfactant, (5.16) reduces to that determined by Bertozzi & Brenner (Reference Bertozzi and Brenner1997). We can rewrite (5.16) using (4.20) as

The integrand in the second term in brackets in (5.17) is positive if $H_{s}>1$ for a considerable part of the domain. Hence, a large capillary ridge is necessary for this to happen, as identified by Bertozzi & Brenner (Reference Bertozzi and Brenner1997) for gravity-driven spreading. We numerically approximate the integrals involved to compute the growth rate ${\it\beta}_{2}$ . For large inclination angles, the second integral in (5.17) is positive and dominates the first, which is negative, resulting in a positive growth rate. The negative contribution from the first integral is due to the steep drop in surfactant concentration observed immediately ahead of where the base state surfactant concentration has a maximum (see inset in figures 2 d and 10 b). This results in a net downward flow out of the contact-line region, which has a stabilising influence on the growth of perturbations at this order. As the inclination angle decreases, the positive contribution from the second integral gradually decreases, becoming negative below some threshold value of ${\it\theta}$ , and the contribution from the second integral decreases in magnitude while still being negative. Hence, below this threshold value, ${\it\beta}_{2}$ is always negative. Figure 16(a,b) shows the $O(q^{2})$ approximation ${\it\beta}\approx q/8+{\it\beta}_{2}q^{2}$ (dot-dashed lines) along with the numerical dispersion relation (solid lines) for inclination angles ${\it\theta}=90^{\circ }$ and 3°, respectively, with ${\it\beta}_{2}=0.0207$ for ${\it\theta}=90^{\circ }$ and ${\it\beta}_{2}=-0.0847$ for ${\it\theta}=3^{\circ }$ . We observe that the $O(q^{2})$ correction provides a much better approximation of the dispersion relation at small wavenumbers.