## Abstract

We examine theoretically the spreading of a viscous liquid drop over a thin film of uniform thickness, assuming the liquid’s viscosity is regulated by the concentration of a solute that is carried passively by the spreading flow. The solute is assumed to be initially heterogeneous, having a spatial distribution with prescribed statistical features. To examine how this variability influences the drop’s motion, we investigate spreading in a planar geometry using lubrication theory, combining numerical simulations with asymptotic analysis. We assume diffusion is sufficient to suppress solute concentration gradients across but not along the film. The solute field beneath the bulk of the drop is stretched by the spreading flow, such that the initial solute concentration immediately behind the drop’s effective contact lines has a long-lived influence on the spreading rate. Over long periods, solute swept up from the precursor film accumulates in a short region behind the contact line, allowing patches of elevated viscosity within the precursor film to hinder spreading. A low-order model provides explicit predictions of the variances in spreading rate and drop location, which are validated against simulations.

## 1. Introduction

The thin liquid film lining lung airways plays an important role in protecting airway tissues from the harmful effects of inhaled particles or aerosol droplets [1–3]. The film is a complex liquid that includes mucins, surfactants and surfactant-associated proteins; its thickness is regulated by osmotic effects driven by ion fluxes across airway epithelial cells and its transport is driven by active motion of cilia on epithelial cells. The film’s rheology is dependent in part on the secretion of mucins from goblet cells distributed across the airway wall; disruption of normal mucin production can lead to harmful effects associated with poor clearance of pathogens. The physical properties of the film in a particular airway of a given individual are therefore subject to considerable uncertainties and intrinsic spatial variability [4,5].

These features motivate this study, in which we seek to relate the spatial heterogeneity of a liquid film to the dynamics of a drop spreading over it. We deliberately focus on a subset of features relevant to airway liquid, neglecting non-Newtonian rheology, osmotic effects, surface-tension gradients, internal stratification and ciliary transport. Instead, we assume that the film’s viscosity is determined by the concentration of a solute (a proxy for mucins, strong determinants of mucus viscosity [6]) that is distributed heterogeneously and diffuses slowly within the film. We wish to establish how spatial variability in solute concentration influences the rate at which an inhaled aerosol droplet might spread over the film. This allows us to address an equivalent, related question: given imperfect knowledge of the film’s properties, what is the likely distribution of spreading rates?

To investigate these questions, we exploit a sequence of approximations. The drop and the film over which it spreads are both assumed to be thin, allowing the flow to be described using lubrication theory, but macroscopic, so that disjoining pressure effects can be neglected. We assume the solute is of an appropriate molecular weight to diffuse across the film during the lifetime of the drop spreading, but not appreciably along it. This allows us to simulate the spreading flow using a pair of coupled transport equations, for the film thickness and the cross-sectionally averaged solute distribution. The initial solute distribution along the film is described as a Gaussian random field with a specified covariance. We can then simulate multiple (Monte Carlo) realizations of spreading dynamics, although this is computationally expensive. Further progress can be made by assuming the drop height significantly exceeds the precursor film thickness. As is well known from numerous studies (reviewed in [7–10]), the drop dynamics is then regulated by the flow in narrow ‘inner’ regions in the neighbourhood of the drop’s effective contact lines. Analysis using characteristics shows that the solute distribution ahead of the drop is swept into each inner region where it concertinas as the drop advances over the film; in contrast, the solute distribution within the remainder of the drop is stretched by the spreading. At any instant, the drop spreading rate is regulated primarily by the viscosity at the rear of each inner region, where it overlaps with the bulk ‘outer’ region. We exploit these observations to derive a set of nonlinear ODEs (a surrogate of the full system) that captures drop spreading rates and which allows the statistical variability in the dynamics to be characterized efficiently. Further simplifications arise when the disorder in the initial viscosity field is weak.

Our study complements numerous previous theoretical studies of drop spreading and contact-line motion. The precursor film regularizes the contact-line singularity [11,12]; we avoid introducing slip or a disjoining pressure, while recognizing that these may be relevant in some applications. While there are numerous potential origins of randomness (thermal fluctuations [13], a rough surface [7,14–18], etc.), we focus here on spatial heterogeneity of the liquid itself, a feature that is particularly relevant to biological applications. The problem is governed by four primary dimensionless quantities (a precursor film thickness; a Péclet number; the variance and correlation length of the initial random field); rather than attempt comprehensive coverage of parameter space, we investigate distinguished limits in which insights are possible through model reduction techniques.

## 2. The model problem

We consider the evolution of a thin liquid film having spatially heterogeneous viscosity. The film lies on a flat plane and spreads under the action of surface tension alone. The liquid wets the plane (with zero equilibrium contact angle) and satisfies the no-slip condition at its lower surface; its upper surface is free of external stress. The liquid is assumed to have Newtonian rheology but contains a chemical species, which is transported passively, such that the liquid’s viscosity is linearly proportional to the chemical concentration (the surface tension being unaffected). Provided the film is sufficiently thin, lubrication theory can be used to derive a nonlinear evolution equation for the film thickness *H*(*X*,*T*), as a function of distance *X* along the plane and time *T*. Molecular diffusion is assumed sufficiently strong to suppress transverse but not axial concentration gradients of the chemical species, so that its cross-sectionally averaged concentration, and thus the cross-sectionally averaged solute field *a*and
*b*The evolution equation for *H* describes how fluid is transported by surface-tension-induced pressure gradients associated with gradients of interfacial curvature, at a rate modulated by the local viscosity field; this field is transported by bulk advection and can spread along the film via molecular diffusion. The Péclet number, *Pe*, measuring the strength of advection to diffusion, is chosen to be sufficiently large for axial diffusion to appear only as a weak singular effect. In the absence of diffusion, (2.1b) can be expressed in conservative form for the transported variable

To illustrate the impact of heterogeneous viscosity, we consider spreading of a droplet sitting on a precursor film. The drop has an initial parabolic profile with (dimensional) height *h*_{0} and half-width *l*_{0}, from which we define an aspect ratio *ϵ*=*h*_{0}/*l*_{0}≪1; the precursor film surrounding the drop has thickness *ηh*_{0}, where *η*≪1. We do not attempt to model the impact of the drop with the film or subsequent mixing of material. The initial condition on *H* is simply
*ω* is an event in an underlying probability space. For fixed *X*, *ω*, *X* that we call the sample associated with *ω*. We assume *σ*^{2} is the variance of the Gaussian random field and *l* is the correlation length of the initial viscosity distribution. The squared exponential covariance function (2.3) yields smooth samples of the Gaussian random field, facilitating numerical simulations. For convenience, we do not label variables *H*, *ω* although this will be implicit. We wish to establish how the uncertainties in *σ*^{2} and *l*, propagate through (2.1), when the precursor film is vanishingly thin (*η*≪1) and diffusion is weak (1≪*Pe*≪*ϵ*^{−1}).

To close the problem we impose no-flux conditions at |*X*|=*L* for some *L*≫1, ensuring that the film sufficiently far from the drop remains undisturbed as the drop spreads. To perform numerical simulation, we draw a sample of

Results of simulations are presented in §3 and in figures 1–4. Figures 1–4 also include approximations from a low-order model, derived in §4. We restrict our attention to times over which the drop remains significantly thicker than the precursor film.

## 3. Simulations

Figure 1 shows an example of drop spreading given a sample of *l* is shorter than the drop width, and the variance *σ* is sufficiently large to ensure that variations in the film’s viscosity span an order of magnitude. The drop retains a parabolic profile as it spreads. Insets near each contact line (figure 1*a*) show a characteristic dimple in the film thickness where the drop connects to the precursor film. We use the local minimum to identify the contact-line locations, defining *X*=*a*_{±}(*T*) to be the locations at which *H* reaches its first minimum as *X* increases (decreases) from the drop centre, where *a*_{−}(*T*)<*a*_{+}(*T*). We use these variables to characterize the drop width *b*), the two contact lines travel at slightly different speeds: in this example, the left-hand contact line has travelled a little further than the right-hand contact line (|*a*_{−}(10)|>*a*_{+}(10)); the region of high viscosity near *X*=1.6 appears to restrain the motion of the right-hand contact line.

The simulation in figure 2 demonstrates how the drop behaves when the correlation length *l* is large compared with the drop width. The viscosity is larger on the right-hand side of the drop, leading to slight leftward displacement of the drop centre as it spreads (figure 2*d*). As in the majority of cases investigated, the bulk velocity of the spreading flow *e*) is approximately linear beneath the drop, falling abruptly to zero (with small flow reversal) near each contact line.

Transport of the *e*) stretches the *b*, where the points *A*_{±}′ bound the region in which *e*), the *b* and 2*b*. In figure 2*b*, symbols mark locations at which *a*), with crowding of characteristics near the contact line (figure 3*b*), leading to rapid variation of

Figure 4 presents statistics describing drop spreading over multiple realizations of the initial viscosity field. We use 1000 samples to estimate the standard deviation of the drop centre and width (*T*=100 and assess the dependence on the variance *σ*^{2} and correlation length *l* of the initial viscosity field; the means of *σ* or *l* in this example. For the present, we focus on the square symbols, denoting predictions from simulations of (2.1). Figure 4*a*,*b* shows that, as *σ* increases, *σ* were limited by the difficulty of resolving very large viscosity gradients that accumulated in the contact-line region, for the chosen value of *Pe*.) The standard deviations show notable dependence on the correlation length (figure 4*c*,*d*): for small *l*, the viscosities at the left and right contact lines are uncorrelated, whereas they become increasingly similar as *l* increases. Consequently, fluctuations in drop width increase in magnitude as *l* increases (if one contact line is, say, hindered, then the other is also likely to be), whereas there is less tendency for the drop to drift sideways (the mean drift remains very close to zero). This behaviour is consistent with studies of drops spreading on random surfaces [18]. As *l* becomes very small, *σ* and *l* using an asymptotic model below.

## 4. Derivation of a low-order model

With *η*≪1, the drop motion is slow and is dominated by the flow in the neighbourhood of the contact lines. We now investigate the impact of readjustment of the viscosity field on this motion, initially neglecting the influence of axial diffusion in (2.1). We divide the flow into an outer region in which the drop adopts an equilibrium shape to leading order, with narrow regions at each contact line (illustrated by insets to figure 1*a*) governed by a modified form of the Landau–Levich equation. While the overall structure of the flow follows the uniform-viscosity case [10,17,19], we seek to identify how variations of film properties modify the drop spreading rate. We follow previous authors [17,20] in matching the cube of the interface slope between inner and outer regions, rather than invoking an intermediate region. Formally, we assume that *σ* and *l* are *O*(1) as

### (a) Outer region

For *X*>*a*_{+} and *X*<*a*_{−}, away from the contact lines, *H*=*η*, *a*_{−}(*T*)<*X*<*a*_{+}(*T*), we seek a solution of (2.1a) subject to
*T*=*O*(1)), we expand *G* as *G*=*G*_{0}+*G*_{1}+⋯ , where *G*_{0} is the quasi-static solution
*G*_{1} is linear in *Y* yields
*c*. Thus, *b* and 2*b*. It is therefore reasonable to identify *T*=0); the value of *M*_{±} may be affected by axial diffusion in practice.

*G* is singular as *Y* gives
*ζ*_{±} requires use of conditions (4.6), which cannot be carried out analytically for arbitrary *N*(*Y*). However, following [17] and multiplying (4.7) by (1−*Y*)(1+*Y*)^{2} and integrating by parts with respect to *Y* from −1+*ϵ*_{−} to 1−*ϵ*_{+} (*ϵ*_{−} and *ϵ*_{+} are small and positive), we have
*a*and
*b*(provided *N* is sufficiently smoothly varying), we simplify (4.11) to find, as *ϵ*_{±}→0,
*ζ*_{−} follows analogously from multiplying (4.7) by (1+*Y*)(1−*Y*)^{2}. Here,

The asymptotic behaviour of *H*_{X} (the inner limit of the outer solution) is obtained as

### (b) Inner region

For convenience, we restrict attention initially to the inner region at the right-hand contact line. Here, we stretch *X* and enlarge *H* using
*η*≪1,

Matching to the outer region requires *M* help us to construct far-field asymptotic solutions of (4.18). One boundary condition can be simplified by writing

Denoting *α*_{1} is a constant dependent on the whole *α*_{2} is also a constant. The solution (4.22) represents a one-parameter family of solutions of (4.19); only one member of that family satisfies (4.21). We solve (4.19) numerically to determine *α*_{2}, and hence *α*_{1}. Shooting towards *α*_{1} by solving

When *α*_{1}≈−0.63≡*A*_{1}, say, in accordance with prior studies (Peng *et al.* [21], for example, report a value equivalent to −0.61 using a three-term expansion over an unspecified domain). We now consider how *α*_{1} is influenced by spatial variations of the *α*_{1} are illustrated in figure 5. For values of *r*_{M} close to unity, the approximation *α*_{1}≈*A*_{1} is robust, particularly when *r*_{M}≪1 (representing a drop spreading into a high-viscosity region) and *r*_{M}≫1 (a drop encountering a low-viscosity region): *α*_{1} becomes large and positive in the former case (with *α*_{1}∝1/*r*_{M} for fixed *X*_{r}≤0) but remains relatively small and negative in the latter. The jump in viscosity therefore has greatest influence on the magnitude of *α*_{1} when the jump extends into regions where the film is thicker and when the contact line is encountering a region of elevated viscosity. We explore the consequences of these variations below.

### (c) Matching

Written in outer variables, the outer limit of the inner problem (4.21) is
*α*_{±} denotes the values of *α*_{1} at each contact line. Matching (4.16) with (4.24) gives the coupled ODEs
*a*where
*b*The system (4.25) constitutes a simplified model for the slow spreading of a drop over a heterogeneous film, and recalls similar descriptions of drop-spreading on homogeneous films [22], for which *M*_{±}=1; this limit yields a deterministic expression for *a*_{±}=±*a*_{0}, namely
*t*^{1/7} in a planar geometry).

Equation (4.25) indicates that the contact-line speeds are coupled via hydrodynamic effects within the bulk drop, and are dependent primarily on the viscosities *M*_{±} upstream of each contact line. This becomes evident on taking leading-order terms as *α*_{±} in (4.25b): if variations are modest, we adopt the approximation
*α*_{±} and, potentially, *M*_{±}.

To solve (4.25) and (4.28), we transform them to a system of differential-algebraic equations by defining *a*–*c* compares the asymptotic predicted drop shape, width *η*, (4.25) provides notably greater quantitative accuracy than (4.27).

The limitations of the assumption *α*_{±}=*A*_{1} are illustrated in figure 1. In this example, *M*_{+}<*M*_{−} (compare the viscosities at *A*_{±}), leading to initial rightward drift of the drop. However, the peak in viscosity near *X*=1.6 causes the right-hand contact line to slow and the drop to then drift left. In this example, the viscosity field within the inner region (between *A*′_{+} and the open circle in figure 1*b*) has a ramp with magnitude *r*_{M}≈0.2. As it propagates into the contact-line region, *α*_{+} can be expected to increase as indicated in figure 5. The dominant terms in (4.25) when *η*≪1 and *α*_{+}≫1,

### (d) The bulk velocity field

To understand solute transport in more detail, we use the asymptotic approximation to describe the bulk velocity field. In the inner region, the velocity is
*e* and is used to construct the characteristics along which *a*_{±}(*T*) from (4.25) and (4.28), and the numerical solution of *M*′=1. It is clear from figure 3*b* that *X* derivative of *a*_{+}, placing it formally within the overlap region between the inner and outer solutions. This defines the boundary between expansive and compressive regions (figure 3). Thus, in the absence of diffusion, all the solute swept up by the contact line remains confined to a narrow zone immediately behind the contact line, within which the asymptotic inner region is confined.

### (e) Weak disorder

When *σ*≪1, we can use a perturbation method to quantify the variability in solutions of (4.25) and (4.28) explicitly. We write the random variables *M*_{±} and *a*_{±} as sums of their mean and a small random perturbation
*a*_{0} satisfies (4.26). We anticipate that perturbations are *O*(*σ*) smaller than leading-order terms and *M*_{±1}〉=0. For simplicity, we restrict attention to leading-order terms as *a*_{±1} give, in terms of the drop displacement *a*and
*b*using (2.3) directly to evaluate covariances.

Predictions of the PDE system (2.1), the ODE system (4.25) and (4.28) and the explicit expressions (4.35) are compared in figure 4. The linear dependence of *σ* is reflected by simulations, except for larger variance where the assumption *l* where again the effects of axial diffusion are likely to become significant. The present predictions could be refined by incorporating *a*_{0} in (4.26), which would also include the influence of the initial viscosity distribution across the bulk of the drop through the integrals *M*_{±}. However, we focus instead on incorporting the effects of axial diffusion.

## 5. An approximate model for axial diffusion

Compression of the viscosity field at the contact line limits the time over which computations can be pursued (for fixed *Pe*≫1), particularly when initial solute gradients are large. Naturally, axial diffusion can be expected to have a growing influence in each compressive region as time increases. We now develop a simple model to describe drop motion when diffusion is sufficient to homogenize the solute in the short compressive regions behind each contact line.

We model the wedge behind the contact line, in which the flow is compressive, by taking *H*≈*θ*_{+}(*a*_{+}−*X*)+*η* for *b*_{+}<*X*<*a*_{+} and *H*=*η* for *X*>*a*_{+}, where *b*_{+} represents the rear boundary of the wedge, defined below. (For clarity, we initially consider only the right-hand contact line.) In this simple compartmental approach, *θ*_{+} represents the approximate contact angle at the edge of the outer region. The inner region near the front of the wedge at *X*=*a*_{+} has length *η*/*θ*_{+}, which is short compared with *a*_{+}; we introduce *ε*_{+}=*η*/(*a*_{+}*θ*_{+})≪1. Within the inner region *H*=*η* and *ε*_{+}→0, of the evolving solute concentration in the wedge as
*α*_{1}=*A*_{1} and *M*_{±} replaced by

The weak-disorder limit of this model is particularly revealing. Writing *a*_{±}=±*a*_{0}+*a*_{±1}+⋯ and *a*_{0}(0)=1 and *a*_{0} increases the latter component dominates the former, showing a fading memory of the initial condition. Because *a*
*b*
*c*
*d*

Figure 4 shows how the modified system captures the reduction in *l* falls to zero. For *l*=*O*(1), (5.7) recovers (4.35) for *a*_{0}≫1 at leading order. Likewise, both variances are proportional to *l*≪1 with *a*_{0}=*O*(1). The increase in variance with *l* when the correlation length is short compared with the drop radius can be explained as follows: diffusion acting within each wedge will tend to suppress the effect of fluctuations in the accumulated viscosity field; however, increasing the correlation length suppresses this effect, in each contact line independently, promoting variation in drop location and width. The approximations (4.35) and (5.7) indicate how variances change with time as the drop spreads (through their dependence on *a*_{0}), as long as the drop stays thick compared with the precursor film.

## 6. Discussion

Complex liquids in natural environments can have spatially heterogeneous properties that influence, and are transported by, a flow. Although diffusion can be expected to suppress spatial gradients over long timescales, heterogeneity will persist in liquids containing large molecular-weight structures with low mobility. In practical applications, the heterogeneity can often be quantified at best at a statistical level, requiring flow outcomes to be described in terms of distributions. The example we present here illustrates some of the challenges of this task. Regions of strong compression quickly generate large spatial gradients in the transported material, far narrower than the physical boundary layers within the flow, which rapidly inflate computational cost; this cost is magnified by the requirement to simulate multiple realizations of the problem.

The example we consider here, motivated by an application in respiratory physiology, illustrates the benefits (and limitations) of low-order approximations of the flow, which can be used to predict outcomes and their variability. When heterogeneity is weak, drop spreading rates are determined primarily by conditions near each contact line at the start of the spreading process; the drop ‘samples’ restricted features of the initial viscosity distribution and these have long-lived influence. In this case, we were able to derive explicit expressions for the mean and variance of variables describing the drop’s motion, in terms of parameters describing the structure of the initially heterogeneous liquid. A more complex picture emerges for a strongly heterogeneous liquid, for which spreading is inhibited by patches of elevated viscosity encountered by the advancing contact lines. In this case, we derived an ad hoc model that shows how a reservoir of solute immediately behind the contact line regulates spreading rates over long timescales.

We have focused attention on a parameter range that is accessible to analysis. A slender geometry allowed the use of lubrication theory and the assumption of a fully wetting fluid interacting with a pre-wetted surface obviated the need to include disjoining pressure. We assumed a linear relationship between viscosity and the distribution of a passively transported solvent, and assumed that the solvent had a sufficiently large molecular weight for diffusion to suppress gradients across but not along the liquid layer. We also assumed that the correlation length of the initial solute distribution exceeded the film thickness. The resulting system of coupled nonlinear hyperbolic/parabolic PDEs generates solutions with large localized gradients requiring careful numerical treatment. To gain physical insights, we derived asymptotic approximations exploiting the difference between the height of the drop and the depth of the precursor film over which it spreads. This yielded a hierarchy of algebraic/ODE systems for the location of the two contact lines, which can predict the mean and variance of drop width and lateral displacement. Naturally, many features arising in applications should be addressed in future studies, not least spreading in two spatial dimensions.

The ODE model revealed the mechanism whereby drop motion is arrested when encountering a patch of elevated viscosity. The viscosity field is initially steepened within an inner layer near the contact line (of approximate width *ηt*^{2/7} at large times—this is short compared with the drop width of approximate order *t*^{1/7}, assuming *t*≪*η*^{−7}), leading to a ramp-like distribution passing slowly towards the rear of the inner region. (The shear rate in the inner region is order 1/(*t*^{8/7}*η*), sufficient to lead to exponentially rapid compression of the viscosity field.) As the more viscous liquid invades the inner region, the matching parameter *α*_{±} increases in magnitude (figure 5), causing slowing of the contact line’s motion. Compression of the viscosity field extends into the overlap region between the inner and outer problems (a distance of order (*ηt*^{3/7})^{1/2} behind the contact line). Thus, in the absence of longitudinal diffusion, characteristics are confined within this ‘compressive wedge’. In practice, diffusion can be expected to suppresses gradients in this wedge over large times. In this case, spreading rates are increasingly regulated by features of the viscosity field encountered during the spreading process. The increasing influence of material encountered during spreading is neatly illustrated by (5.6).

In terms of the application motivating this study, the primary insight concerns the role of the precursor (mucus) film in regulating the spreading of an inhaled drop of a different material. Assuming the liquids are fully miscible, the drop will slowly accumulate endogenous material at its leading edge and fluctuations in its motion can therefore be associated directly with initial heterogeneities in the precursor film. As figure 4*e* shows, fluctuations in drop location (relative to drop radius) are most pronounced when the drop radius is comparable to the correlation length of the viscosity distribution. From a methodological perspective, our study shows how low-order physical models, combined with weak disorder expansions, can be effective in quantifying the statistical variability in flow outcomes.

## Data accessibility

Simulation codes of (2.1) and (4.19) with (4.23) are deposited in Dryad: http://dx.doi.org/10.5061/dryad.t2v1b.

## Authors' contributions

O.E.J. and F.X. jointly conceived the study; F.X. performed numerical simulations; O.E.J. and F.X. jointly undertook the model development and analysis; O.E.J. and F.X. jointly drafted the manuscript.

## Competing interests

We declare we have no competing interests.

## Funding

This work was supported by EPSRC grant no. EP/K037145/1.

## Acknowledgements

We are grateful to a referee for spotting an error in an earlier version of this work.

## Appendix A

**(a) Model derivation**

We consider a liquid layer bounded below by a horizontal solid substrate and above by a free surface. We introduce Cartesian coordinates (*x*,*z*) such that the solid substrate lies on *z*=0 and the free surface occupies *z*=*h*(*x*,*t*), where *t* is time. The liquid is incompressible, Newtonian and has dynamic viscosity *μ*(*x*,*z*,*t*), constant density *ρ* and uniform surface tension *σ*. The Reynolds number is assumed sufficiently small for inertia to be neglected, so that the flow field (*u*,*w*) and pressure *p* satisfy the Stokes equations, which when accounting for spatially varying viscosity may be expressed as
*a*
*b*
*c*We assume the viscosity *μ*=*μ*(*c*(*x*,*z*,*t*)) is determined by the distribution of a chemical species with concentration *c*(*x*,*z*,*t*), which satisfies the transport equation
*γ* is a constant diffusivity. No-flux conditions are imposed on *c* at *z*=0 and *z*=*h*.

We introduce a length scale *l*_{0}, height scale *h*_{0} (defining *ϵ*≡*h*_{0}/*l*_{0}≪1), viscosity scale *μ*_{0}, velocity scale *u*_{0}=*ϵ*^{3}*σ*/*μ*_{0} and concentration scale *c*_{0}, and non-dimensionalize variables using (*x*,*z*,*h*)=*h*_{0}(*X*/*ϵ*,*Z*,*H*), (*u*,*w*)=*u*_{0}(*U*,*ϵW*), *t*=(*h*_{0}/*ϵu*_{0})*T*, *p*=(*μ*_{0}*u*_{0}/*ϵh*_{0})*P*, with *μ*=*μ*_{0}*M*(*C*) and *c*=*c*_{0}*C*, where the function *M*(*C*) is to be specified. We define the Péclet number *Pe*=*l*_{0}*u*_{0}/*γ*.

The governing equations and boundary conditions become
*a*
*b*
*c*
*d*
with *M*=*M*(*C*). Here we have eliminated terms of *O*(*ϵ*^{2}) from the flow equations (as is standard in lubrication theory) but have retained all terms in the solute transport equation. It follows that *P*(*X*,*T*)=−*H*_{XX}, and integration of the horizontal momentum equation yields

We define an averaging operator on the function *Φ*(*X*,*Z*,*T*) as

We introduce the decomposition

We seek the limit in which *Pe*≫1 as *ϵ*→0. Anticipating the dominant balance in (A 12) to be *C*′=*O*(*ϵ*^{2}*Pe*), the terms in (A 11) fall into four categories with magnitude *O*(1) (advection), *O*(1/*Pe*) (diffusion), *O*(*ϵ*^{2}*Pe*) (Taylor dispersion) and *O*(*ϵ*^{2}). Thus for 1/*ϵ*≫*Pe*≫1, the approximation of (A 11) up to *O*(1/*Pe*) is
*O*(1/*Pe*) contribution of diffusion in (A 13) to facilitate numerical simulations. Assuming that the viscosity *M* linearly depends on *C* yields (2.1b) and *O*(1/*Pe*,*ϵ*^{2}*Pe*,*ϵ*^{2}), as in (2.1a).

A similar formulation has been adopted by Karapetsas *et al.* [23] in a study of thin-film suspension flow, for which a nonlinear relation between viscosity and particle concentration was retained. As in that study, we assume here that the solute field does not influence the interfacial tension; for suspensions, linearization of *M*(*C*) is appropriate at low volume fractions [24].

**(b) Karhunen–Loéve decomposition**

We use the Karhunen–Loéve decomposition to sample the Gaussian random field *X*_{i}}, *i*=0,1,…,*n*, on the computational domain [−*L*,*L*], the covariance function *i*,*j*=0,1,…,*n*, which can be factorized as *K*=*V* *ΛV* ^{T}, where *Λ* is the (*n*+1)×(*n*+1) diagonal matrix of eigenvalues, λ_{0}≥λ_{1}≥⋯≥λ_{n}≥0, of *K*, and *V* =[*v*_{0},…,*v*_{n}] is the (*n*+1)×(*n*+1) matrix whose columns *v*_{i} are the eigenvectors of *K*. To resolve the *M* field properly, we choose the grid width such that it is smaller than one fifth of the correlation length *l*. As in [25], the discrete random field *ξ*_{i} are independent and identically distributed Gaussian random variables with zero mean and unit variance. For large *n*, we further truncate the sum in (A 15) after *n*′ (≪*n*) terms and use the approximate discrete random field
_{n′+1},…,λ_{n}. Here, we choose smallest *n*′ such that λ_{n′}/λ_{0}<10^{−3}.

**(c) Compressive wedge approximation in the weak disorder limit**

Writing *X*=(*a*_{0}−1)/*N*, *X*_{j}=1+(*j*−1)Δ*X* and *N* a large positive integer. Clearly,

- Received April 16, 2016.
- Accepted September 15, 2016.

- © 2016 The Authors.

Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited.