Drop spreading with random viscosity

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.


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][2][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, surfacetension 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 crosssectionally 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][8][9][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][15][16][17][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.

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 fieldM(X, T) (which for convenience we call the viscosity field), are transported by the cross-sectionally averaged fluid velocityŪ(X, T). As demonstrated in appendix Aa, these equations (in a planar geometry) may be expressed in dimensionless form as 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 HM, which represents the amount of solute per unit length of film.
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 The initial viscosity fieldM(X, 0) is represented as a random field M(X, ω), where ω is an event in an underlying probability space. For fixed X, M is a random variable; for an outcome ω, M is a function of X that we call the sample associated with ω. We assume M = exp(G(X, ω)), where G(X, ω) is a Gaussian random field with zero mean and stationary covariance Here, σ 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,M, etc., with ω although this will be implicit. We wish to establish how the uncertainties in M, represented by σ 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 M (constructed using a Karhunen-Loéve decomposition, see appendix Ab); using this as an initial condition forM(X, 0), we solve (2.1) numerically with the method of lines using fourth-order spatial differences. We collect results from multiple runs to compute statistics (such as mean and variance) of quantities of interest, which we compare with predictions of asymptotic analysis.
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 M, for which the correlation length 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 1a) 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 W(T) and lateral displacement of its mid-point C(T) as the drop spreads, defined by W(T) = a + − a − and C(T) = 1 2 (a + + a − ). (3.1) Because the initial viscosity distribution is heterogeneous (figure 1b), 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     Thus, the linear stretching flow beneath the spreading drop (Ū X > 0, figure 2e) stretches theM field laterally without changing its magnitude. This is illustrated in figure 1b, where the points A ± bound the region in whichŪ X > 0. Ahead of the drop theM field is undisturbed, whereas near the contact line, where the flow is strongly compressive (Ū X < 0, see insets in figure 2e), theM field steepens. This compression is evident from the distributions in figures 1b and 2b. In figure 2b, symbols mark locations at whichM(X, T) =M(±1, 0), demarcating the boundary between stretching and compression of theM field. This is further demonstrated by the pattern of characteristics, which shows uniform stretching of the concentration field beneath the drop (figure 3a), with crowding of characteristics near the contact line (figure 3b), leading to rapid variation ofM in this region. Weak axial diffusion can be expected to suppress such gradients over long times. 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 (σ C , σ W ) at T = 100 and assess the dependence on the variance σ 2 and correlation length l of the initial viscosity field; the means of C and W do not show appreciable dependence on σ or l in this example. For the present, we focus on the square symbols, denoting predictions from simulations of (2.1). Figure 4a,b shows that, as σ increases, σ C and σ W increase. (Simulations for larger σ 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 4c,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, σ C falls; this reflects the effects of axial diffusion in simulations suppressing sharp gradients in the solute field. We seek to quantify the dependence of σ C and σ W on σ and l using an asymptotic model below.

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 1a) 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 η → 0, ensuring that the correlation length of the initial viscosity field exceeds the width of the contact-line regions.

(a) Outer region
For X > a + and X < a − , away from the contact lines, H = η,M =M(X, 0) andŪ = 0. Within the drop, with a − (T) < X < a + (T), we seek a solution of (2.1a) subject to lim where V = 4 3 is the volume of the droplet (and the ∓ symbol here denotes a one-sided limit). Assuming the drop shape to be quasi-static, we write and G 1 is linear inȧ ± , satisfying which we seek to solve subject to Substituting (4.4) into (4.5) and integrating (4.5) once with respect to Y yields The bulk fluid velocity in the outer region is thereforē where H(·) is Heaviside function. The linear stretching flow is illustrated in figure 2c. Thus,Ū out = C + 1 2Ẇ Y ≡ dX/dT, implying thatM = N(Y) exactly satisfies (3.2). The viscosity field beneath the bulk of the drop is stretched linearly, as illustrated in figures 1b and 2b. It is therefore reasonable to identify M ± (T) = lim Y→±1 N =M(a ± (0+), 0), where 0+ denotes the early time at which the asymptotic spreading structure is established (which we take here to be T = 0); the value of M ± may be affected by axial diffusion in practice.
G is singular as Y → ±1 and (4.7) shows its asymptotic behaviour to be Integrating (4.9) twice with respect to Y gives Finding the constants ζ ± 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 Using (4.6) and (4.10) and noting that (provided N is sufficiently smoothly varying), we simplify (4.11) to find, as ± → 0, where the result for ζ − follows analogously from multiplying (4.7) by (1 Here, 14) The asymptotic behaviour of H X (the inner limit of the outer solution) is obtained as (recalling [17,19]), so the cube of the slope may be written in original variables as where α 1 is a constant dependent on the wholeM field in the inner region and φ(X) ≡ ln(−X). Aŝ X → +∞, (4.19) has an asymptotic solution of the form where α 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 towardsX → −∞, we seek the solution satisfyingXH 2XHXX = −r M at the left of the domain, in accordance with the asymptotic solution (4.21). We evaluate α 1 by solvingH 3X + 3r M (ln(−r 1/3 MX ) + 3 2/3 α 1 + 1) = 0 at the domain boundary. WhenM ≡ 1, for example, we find that α 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 theM field. Because the inner-region flow is compressive, a region of high or low viscosity encountered by an advancing contact line will typically manifest as a steep ramp in theM field; the ramp will slowly propagate from the front to the rear of the inner region. To mimic this situation, we chosê for someX r , imposing continuity conditions inĤ,ĤX,ĤXX acrossX r , which we solved on a long domain [−X L ,X R ]. Computed values of α 1 are illustrated in figure 5. For values of r M close to unity, the approximation α 1 ≈ A 1 is robust, particularly whenX r > 0. However, there is a striking difference between cases in which 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 so that in the overlap regions η|ȧ ± | −1/3 |a ± − X| 1, where α ± denotes the values of α 1 at each contact line. Matching (4.16) with (4.24) gives the coupled ODEs   [22], for which I 1 = 1, I ± 2 = 0, M ± = 1; this limit yields a deterministic expression for a ± = ±a 0 , namelyȧ 0 −2 + 3 2/3 A 1 + ln(2a 0 ) − ln η + To solve (4.25) and (4.28), we transform them to a system of differential-algebraic equations by definingȧ ± as two new variables. Figure 2a-c compares the asymptotic predicted drop shape, width W and centre C with simulations of the PDE system (2.1). The ODE model successfully captures the lateral drift of the drop owing to the gradient in the viscosity field. Factors limiting the accuracy of the low-order model are the inclusion of diffusion in the PDE simulations; given the logarithmic (rather than algebraic) dependence on the small parameter η, (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 1b) 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,ȧ , (4.29) indicate how passage of the region of elevated viscosity backwards into the inner region slows the advance of the contact line.

(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 Noting that lim X→a ± ∓ U out =ȧ ± = limX →−∞ U in± , we can construct a composite approximation ofŪ using (4.8) asŪ This approximation is illustrated in figure 2e and is used to construct the characteristics along whichM is transported in figure 3, using solutions of a ± (T) from (4.25) and (4.28), and the numerical solution ofH from (4.19) with M = 1. It is clear from figure 3b thatŪ com <ȧ + in the neighbourhood of the contact line, implying that characteristics move smoothly through the inner region. Evaluating the X derivative ofŪ com reveals that the velocity maximum lies a distance of order [a + η(ȧ + M + ) −1/3 ] 1/2 (the geometric mean of the inner and outer lengthscales) behind 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 M ± = 1 + M ±1 + · · · , a ± = ±a 0 + a ±1 + · · · , (4.32) where a 0 satisfies (4.26). We anticipate that perturbations are O(σ ) smaller than leading-order terms and M ±1 = G(±1, ω) with M ±1 = 0. For simplicity, we restrict attention to leading-order terms as η → 0, expanding (4.27) using (4.32). Expressions for a ±1 give, in terms of the drop 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 σ C and σ W on σ is reflected by simulations, except for larger variance where the assumption M ± =M(±1, 0) breaks down, because the effects of axial diffusion may also be significant. The dependence on the correlation length is also captured well for l 1, but not at smaller l where again the effects of axial diffusion are likely to become significant. The present predictions could be refined by incorporating O(1/ ln η) corrections to 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 I 1 , I ± 2 in (4.25) and their correlations with M ± . However, we focus instead on incorporting the effects of axial diffusion.

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Ū =ȧ + (1 − (η/H)) (from (4.30)); combining this with the stretching flow in the outer region gives a composite expression for the velocity across the wedge as Thus,Ū X = 0 at b + = a + (ε + + 1) − √ a + ε + W, defining the rear of the wedge. The fluid volume in the wedge, V + = a + b + H dX, satisfies, from (2.1a), The mass of solute in the wedge, At the front of the wedge, H = η andŪ = 0; at the rear, H = ηW/ √ a + ε + W andŪ =ȧ + (W + a + ε + − 2 √ a + ε + W)/W. We then assume that the compressed viscosity field is mixed by diffusion within the wedge, so that the integral in (5.3) may be approximated as N + ≈M + V + andM| b+ =M + in that case (5.2), (5.4) give the leading-order approximation, as ε + → 0, of the evolving solute concentration in the wedge as˙M (treating the left-hand contact line analogously). The reservoir of solute in the wedge is fed by a source in the film ahead, and diluted by expansion of the wedge. Our candidate model for spreading, accounting for the effects of diffusion where the flow is strongly compressive, therefore uses (4.25) with α 1 = A 1 and M ± replaced byM ± (t), supplemented with (5.5).

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 molecularweight 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 4e 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. 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 G, which then gives one sample of M via exponentiation. Given the spatial grid {X i }, i = 0, 1, . . . , n, on the computational domain [−L, L], the covariance function k G produces a covariance matrix K = {k G (X i , X j )}, 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 G := [G(X 0 , ω), . . . , G(X n , ω)] T can be generated as where ξ 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 The quality of the approximation of G ≈ G is determined by the sizes of the neglected eigenvalues λ n +1 , . . . , λ n . Here, we choose smallest n such that λ n /λ 0 < 10 −3 .
(c) Compressive wedge approximation in the weak disorder limit Cov(Y(X 1 ), Y(X j ))