Stratification in drying films: a diffusion–diffusiophoresis model

Abstract This research is motivated by the desire to control the solids distribution during the drying of a film containing particles of two different sizes. A variety of particle arrangements in dried films has been seen experimentally, including a thin layer of small particles at the top surface. However, it is not fully understood why this would occur. This work formulates and solves a colloidal hydrodynamics model for (i) diffusion alone and (ii) diffusion plus excluded volume diffusiophoresis, to determine their relative importance in affecting the particle arrangement. The methodology followed is to derive partial differential equations (PDEs) describing the motion of two components in a drying film. The diffusive fluxes are predicted by generalising the Stokes–Einstein diffusion coefficient, with the dispersion compressibility used to produce equations valid up until close packing. A further set of novel equations incorporating diffusiophoresis is derived. The diffusiophoretic mechanism investigated in this work is the small particles being excluded from a volume around the large particles. The resulting PDEs are scaled and solved numerically using a finite volume method. The model includes the chemical potentials of the particles, allowing for incorporation of any interaction term. The relative magnitudes of the fluxes of the differently sized particles are compared using scaling arguments and via numerical results. The diffusion results, without any inter-particle interactions, predict stratification of large particles to the top surface. Addition of excluded volume diffusiophoresis introduces a downwards flux on the large particles, that can result in small-on-top stratification, thus providing a potential explanation of the experimental observations.

This research is motivated by the desire to control the solids distribution during the drying of a film containing particles of two different sizes. A variety of particle arrangements in dried films has been seen experimentally, including a thin layer of small particles at the top surface. However, it is not fully understood why this would occur. This work formulates and solves a colloidal hydrodynamics model for (i) diffusion alone and (ii) diffusion plus excluded volume diffusiophoresis, to determine their relative importance in affecting the particle arrangement. The methodology followed is to derive partial differential equations (PDEs) describing the motion of two components in a drying film. The diffusive fluxes are predicted by generalising the Stokes-Einstein diffusion coefficient, with the dispersion compressibility used to produce equations valid up until close packing. A further set of novel equations incorporating diffusiophoresis is derived. The diffusiophoretic mechanism investigated in this work is the small particles being excluded from a volume around the large particles. The resulting PDEs are scaled and solved numerically using a finite volume method. The model includes the chemical potentials of the particles, allowing for incorporation of any interaction term. The relative magnitudes of the fluxes of the differently sized particles are compared using scaling arguments and via numerical results. The diffusion results, without any inter-particle interactions, predict stratification of large particles to the top surface. Addition of excluded volume diffusiophoresis introduces a downwards flux on the large particles, that can result in small-on-top stratification, thus providing a potential explanation of the experimental observations. typically preferentially accumulate at the top surface, but it is not fully understood why (Routh 2013). Understanding this could have applications across a wide range of industries, from the surface of catalyst pellets, to coating surfaces with biocides (Mardones et al. 2019). It may have significance in the development of cheaper and more sustainable materials, as the drying process can be engineered such that expensive components are only located where they are required.
As solvent evaporates from a film, the top surface, initially at height H, descends at velocityĖ. The situation and coordinate system used are depicted in figure 1, where z is the spatial coordinate, ξ is the transformed spatial coordinate (see § 3.2) and t is time. Particles which are unable to diffuse away from the top are collected by the top surface. Whether or not the particles can diffuse sufficiently quickly away from the surface is characterised by the Péclet number, Pe = 6πηRĖH/kT, where R is the particle radius, η is the solvent viscosity, k is the Boltzmann constant and T is the temperature (Routh & Zimmerman 2004). Hence it is expected that drying dispersions with different size particles will lead to different particle arrangements in the dried film. On the basis of diffusional arguments alone, it would be expected that larger particles stratify to the top surface.
In light of experimental observations of small particles at the top surface, it is thought that a simple diffusional model may not suffice. The terminology small-on-top stratification is used in this work to describe a layer, including just a thin layer, of small particles at the top surface, and likewise for large-on-top stratification. This work derives a fluid mechanics model for both the simple diffusion case, and the case that adds a type of diffusiophoresis: an excluded volume effect. The aim is to use the numerical results to determine whether diffusiophoresis could explain these experimental observations. It will be shown that diffusiophoresis may result in small-on-top stratification, including a thin layer of small particles at the top surface.
A literature review is carried out in § 1.1 to explain the context of the work and the existing hypotheses. In § 2, the relevant theory is derived, then in § 3, the solution methodology is explained. The results of the simulations are presented and discussed in § 4. Asymptotic solutions are outlined in § 5, before drawing conclusions in § 6.
1.1. Literature review Different approaches have been taken in the literature to model the drying of particulate films and attempt to explain the particle arrangements seen experimentally. Key works which have used thermodynamic models include (i) Routh & Zimmerman (2004), who derived equations for the drying of a film containing just one species of particles and solved them numerically; and (ii) Trueman et al. (2012b), who extended Routh & Zimmerman's work to a film containing particles of two different sizes.
In the expression for the solute chemical potential, μ i , Trueman et al. (2012b) included only an entropic contribution, μ i = μ 0 i + kT ln φ i , where μ 0 i is the reference chemical potential and φ i is the solute volume fraction. Zhou, Jiang & Doi (2017) took a similar approach to Trueman et al., but also included a binary interaction term in the chemical potential. This cross-interaction term in the chemical potential accounts for enthalpic effects, i.e. inter-particle interactions. Zhou et al. (2017) include these up to the order This can predict that the small particles end up on top of the larger particles. However, it does not explain why some images suggest that only a thin layer at the top is made up of small particles (Trueman et al. 2012a). A further shortcoming is that the form of the interaction terms used by Zhou et al. does not account for the change in interactions as the film approaches close packing, i.e. as the solution becomes more concentrated. Zhou et al. consider this to be acceptable, arguing that most of the stratification is formed before the film approaches close packing. Following the approach of Russel, Saville & Schowalter (1989) and Routh & Zimmerman (2004), Trueman et al. (2012b) form particle chemical potential expressions that are valid even in non-dilute systems by using an expression for the solvent chemical potential that diverges at close packing. Having an accurate expression for the solvent chemical potential allows the chemical potentials of the two particles, μ 1 and μ 2 , to be related through the Gibbs-Duhem equation. Therefore, only μ 1 or μ 2 needs to be specified to fully describe the system. In this work we choose to specify ∇μ 1 /∇μ 2 . The approach of Russel et al. will be adopted in this present work, so that the model can be run to close to the end of drying. An alternative to the thermodynamic approach described above is using molecular dynamics. Simulations of this type include diffusiophoresis because they simulate every particle, and do not allow the particles to overlap. In this way, excluded volume, including the small particles being excluded from around the large ones, is part of these models. These simulations can be divided into two types, implicit solvent models, where only the non-solvent particles are directly simulated, and explicit solvent models, where the solvent particles are included. This allows the hydrodynamic interactions to be simulated, resulting in greater accuracy, but at significant computational expense. This present work models volume fraction evolution, rather than every colloidal particle. In doing so, diffusiophoresis is modelled at a coarser level than, for example, molecular dynamics, at which level implicit/explicit solvent methods become important. This is intended to allow different insights to be obtained with the benefit of analytic expressions and without the computational expense.
Considering implicit solvent models, Fortini et al. (2016) used Langevin dynamics with a short-range repulsive interaction. This method predicts layers of small particles on top, but not a single layer of small particles. Fortini et al. argue that this stratification is due to the gradient in osmotic pressure. Similar results were obtained by the Langevin dynamics simulations of Howard, Nikoubashman & Panagiotopoulos (2017a). Note that they obtain film profiles with a thin layer of large particles at the very top surface, which is present in their equilibrium film profiles before commencing drying. They argue that this is due to volume exclusion of the large particles at the interface. This is a well-known entropic effect near a boundary, as, for example, shown by Roth & Dietrich (2000) in comparing Rosenfeld density functional calculations with simulations. Directly beneath is a stratified layer of small particles, which grows over time. In addition to the Langevin dynamics simulations, Howard et al. outline, and subsequently simulate (Howard, Nikoubashman & Panagiotopoulos 2017b), a model based on density functional theory. This can be thought of as an extension to the approach of Zhou et al. (2017), using a high-accuracy equation of state to calculate the chemical potentials. This would be valid for concentrated solutions. Tang, Grest & Cheng (2018) carried out molecular dynamics simulations that model the solvent explicitly, obtaining a modified threshold for the occurrence of small-on-top stratification. Statt, Howard & Panagiotopoulos (2018) argue that hydrodynamic interactions which are modelled in an explicit solvent model are important for reliably predicting stratification: small-on-top stratification is predicted using the implicit solvent method, but not with the explicit solvent one. In contrast, Tang, Grest & Cheng (2019) also compare implicit and explicit solvent models, although with limitations on the film thicknesses that can be assessed with the computationally expensive explicit solvent models. For the systems they studied, the implicit solvent model was found to be sufficient, as both the implicit and explicit solvent models predicted small-on-top stratification.
Adequate modelling of the solvent is important for assessing the effect of diffusiophoresis. A recent hypothesis for explaining the accumulation of small particles at the top surface concerns a type of diffusiophoresis. Diffusiophoresis is the migration of particles along a concentration gradient of a different solute species, which could be either ionic or non-ionic (Anderson & Prieve 1984). Whilst diffusiophoresis can be used as a general term, covering different mechanisms which result in such particle migration, including electrophoresis and chemiphoresis, the hypothesis in this work concerns a specific non-ionic type of diffusiophoresis. Other phoretic terms could be easily added to the model in this work if desired: this model requires an expression for the total flux, so other flux contributions can be added.
The non-ionic type of diffusiophoresis considered is an excluded volume effect, relevant in a mixture of small (component one) and large particles (component two). The small particles are excluded from a layer of solvent of thickness R DP around each of the large particles. For the case of hard spheres, R DP is expected to equal the radius of a small particle, R 1 . However, this work will explore what happens as R DP varies, which may correspond to, for example, non-spherical particles. The situation is depicted in figure 2. The radius of each large particle is R 2 . Note that there will also be excluded volume effects between particles of the same type but, by definition, these are not diffusiophoresis, and are excluded in the present study.
Sear & Warren (2017) outline a simple model for this type of diffusiophoresis in a drying film. It is an idealised model, assuming that particles do not interact and that the particle size ratio is infinite. Diffusion of the larger particles is ignored, on the assumption that it is small compared with the diffusiophoretic flux. Diffusion of the smaller particles is included, by finding the concentration gradient of an ideal gas diffusing between two impenetrable walls, and then calculating the diffusiophoretic flux due to this concentration gradient.
Whilst Sear & Warren (2017) obtain elegant approximate analytical solutions, which are useful for gaining understanding, their model's applicability is limited to dilute solutions and very large size ratios. To aid insight into less idealised situations, particularly where the particle size ratio is not infinite, diffusion of both particles should be included. This allows identification of criteria for when diffusiophoresis is important. To allow the entire drying process to be modelled, not just the early stages, the model could be adapted to be valid up to close packing, again using the approach of Russel et al. (1989).
Diffusiophoresis with different relative strengths of diffusion of the large and small particles was modelled by Ault et al. (2017). The geometry that they chose, particle injection into or withdrawal from a semi-infinite or finite domain, is simpler than the drying geometry and has static boundaries. They formed advection-diffusion equations which are valid for dilute solutions. For the case where the diffusivity of the larger particles is neglected, Ault et al. obtain an analytical solution using the method of characteristics. This is compared with the numerical results which include diffusion of the larger particles. They found that ignoring diffusion of the larger particles accurately predicted the trajectories of the particle fronts, but not the particle concentrations around the fronts. Diffusiophoresis has also been modelled and imaged in a dead-end channel geometry (Shin et al. 2016). This present work will seek to establish the relative importance of diffusion and diffusiophoresis in the geometry of a drying film.
Experimentally, it is difficult to determine whether the top layer is a monolayer of small particles, or thicker. It requires the ability to take measurements throughout the film depth, not just at the top surface. Small-angle X-ray scattering (SAXS) measurements have shown enrichment of small particles in the top several hundred layers of particles (Liu et al. 2019). Existing works which include diffusiophoresis (Fortini et al. 2016;Howard et al. 2017a,b) predict a growing layer of small particles at the top surface, in agreement with this. Experimental observation of a monolayer of small particles would support a different mechanism, such as interfacial effects. Atmuri, Bhatia & Routh (2012) added a surface interaction term to the model of Trueman et al. (2012b) to demonstrate this theoretically.
The review article by Schulz & Keddie (2018) compared predicted stratification regimes based on the interaction potential model of Zhou et al. (2017), and the diffusiophoresis models of Sear & Warren (2017) and Sear (2018), with experiments from the literature. Sear's model considers a jammed layer of small particles at the top surface, with diffusiophoretic drift of the large particles occurring just beneath it. This produces a simple criterion for small-on-top stratification: the speed of the jammed layer must be smaller than that of diffusiophoresis. Good agreement was found between experimental data for drying dilute dispersions and Zhou et al.'s predictions, but less so for concentrated ones. This suggests that a model which is valid for more concentrated solutions but allows input of Zhou et al.'s chemical potential expressions, could improve the agreement. Such a model is developed in this present work. Noting that there was a lack of experimental data in the range required to test the predicted stratification regimes of Sear & Warren and Sear, Schulz et al. (2020) subsequently carried out experiments with films containing colloid and polymer. The data obtained agree qualitatively with the predictions of Sear & Warren and Sear, but not with the exact position of the transition to small-on-top stratification.
It is clear from this survey that several hypotheses could explain the small-on-top observations, although existing simulations seem to predict a thicker final layer than is sometimes seen experimentally (Atmuri et al. 2012). Whilst diffusiophoresis has been suggested as an explanation, it has not been directly modelled in a drying film. This work seeks to address this, adopting a fluid mechanics partial differential equation (PDE) model, since its simple formulation allows insight to be gained.

Diffusion
To make explicit how the theory outlined here compares with previous work, whilst (2.7) in § 2.2 was also used by Trueman et al. (2012b), (2.12)-(2.14) correct the hydrodynamics implementation. These equations are given here in their most general forms. Also, (2.8)-(2.11) add clarification regarding diffusion reference frames. Equations (2.18)-(2.23) are reproduced from the work of Trueman et al., but they are subsequently substituted into the general conservation equation derived in § 2.2. § § 2.3, 3.1 and 3.2 follow the same approach as Trueman et al., although again following through with the general equations.

Explanation of hydrodynamics correction
Using the thermodynamic approach in this present work, it is found that the extension from the one-component case of Routh & Zimmerman (2004) to the two-component case is non-trivial. The Stokes-Einstein diffusion coefficient, D SE , was originally derived for a single particle in infinite solvent, so it needs to be determined how to use D SE in a system that is no longer infinitely dilute. According to Batchelor (1976), D SE applies to a force-free solvent. Hence instead of deriving that the diffusive flux of particles in a one-component film, j 1 , is where −∇μ 1 is the force acting on component one due to its chemical potential, μ 1 , Batchelor derives where φ 1 is the volume fraction of component one. The factor of 1/(1 − φ 1 ) results from a force correction. This force correction could be expressed as where (4/3)πR 3 1 is the volume of a particle of component one, ν s is the volume of a particle of the solvent and μ s is the chemical potential of the solvent. The subscript eff denotes the effective driving force for diffusion. Since the ratio of particle volumes is a constant, (2.3) could be equivalently written as This correction appears in the equations of Routh & Zimmerman (2004). Trueman et al. (2012b) incorrectly extended this by writing where φ 2 is the volume fraction of component two, and μ 2 is its chemical potential. Whilst this appears to be a logical extension of (2.4), (2.3) is still correct for a film with any number of components, following the reasoning of Batchelor (1976 where B ij is the bulk mobility coefficient, which is generally a function of all φ j and ratios of R j . Substituting in expressions for B ij allows an equation of the form of (2.12) to be derived (Batchelor 1983). This present work will use the methodology of Trueman et al. (2012b), but with the corrected (2.12). The numerical results of Trueman et al. solved diffusion fluxes of the same leading order in Pe 1 and Pe 2 as the equations in this present work, and so the results in § 4.1 are expected to show the same qualitative behaviour.

Governing equations
The conservation equation for one type of particles can be written as where N 1 is the total flux of particles of component one. Following the approach of Cussler (2009), the total flux can be arbitrarily split up into a convective and diffusive flux where the volume average velocity of type one particles is U 1 , and their diffusive flux j a 1 relative to the reference velocity U a is The chemical potential of component i in the two-component film is denoted by μ i , and the diffusion coefficient of particles of type one in particles of type i, using a reference velocity U a , is denoted by D a 1i . It is convenient to choose the volume average velocity U v as U a , since U v = 0 in a static film. Hence (2.7) becomes and diffusion coefficients using a volume-fixed reference frame are required. For a volume-fixed reference frame, so expressions for D v 11 and D v 12 are required. It is desired to relate these to the Stokes-Einstein diffusion coefficient of component one, D SE,1 = kT/6πηR 1 .
Hydrodynamic effects are taken into account via K 11 (φ 1 , φ 2 ) and K 12 (φ 1 , φ 2 ), which are referred to as sedimentation coefficients by analogy with sedimentation (Batchelor 1976). It is expected that these will be functions of R 1 and R 2 as well as φ 1 and φ 2 . This gives (2.12) Defining the sedimentation coefficients such that U 1 is a velocity relative to the volume average velocity, and so can be used in (2.10) and (2.11), the conservation equation for component one becomes Hence, D v 11 and D v 12 can be found in relation to D SE,1 by comparison with the sedimentation coefficients. Note that, in this work, the diffusive fluxes are written in terms of chemical potential gradients, as opposed to concentration gradients. This is why we divide D SE,1 by kT in the flux expressions, (2.12) and (2.13), and why a factor of φ 1 appears in front of the ∇μ 1 term in (2.13), and likewise for the φ 2 term in front of the ∇μ 2 .
Expressions for K 11 (φ 1 , φ 2 ) and K 12 (φ 1 , φ 2 ) for rigid spheres with zero interaction potential in dilute solution can be found in Batchelor's work (1983). However, it is desired to use generalised forms of these expressions to run the model up to close packing. A factor of φ 1 is included in front of the ∇μ 2 term in (2.13) in order to match the leading-order coefficients in (φ 1 , φ 2 ) for dilute solution. Batchelor's (1983) dilute expressions will need to be adapted such that the sedimentation coefficients fall to zero as the mixture approaches close packing, due to hydrodynamic hindrance (Russel et al. 1989). Steric effects could be incorporated into the chemical potential terms if desired.
We consider the case of one-dimensional drying with the top surface decreasing normally to the substrate. By scaling withẑ = z/H andt = tĖ/H, (2.13), in one dimension, becomes (2.14)

Evaporation rate
For the purpose of (2.14),Ė needs only to be a constant, characteristic evaporation velocity used for non-dimensionalising. When evaporation at the top surface is included in the boundary condition ( § 3.2), it is assumed that the evaporation rate is approximately constant. Evaporation is driven by the thermodynamic driving force between the vapour pressure of the solvent at the top of the film and the vapour pressure of the gas above the film. The vapour pressure of the gas above the film is affected by the time scale for diffusion from the layer of saturated vapour directly above the film to the bulk gas. This time scale is orders of magnitude smaller than the evaporation time (Popòv 2005), leading to a quasi-static problem (Routh 2013). Routh & Russel (1998) derived an expression forĖ, where k m is the gas-side mass transfer coefficient, p vap is the vapour pressure of the solvent at the top of the film and p ∞ vap is the vapour pressure of the bulk gas. 928 A15-8 The chemical potential of the solvent is related to the osmotic pressure, Π , as where μ 0 s is the solvent reference potential and Z(φ 1 , φ 2 ) is the compressibility, which accounts for the non-ideality of the osmotic pressure at high volume fractions. Relating the chemical potential to the solvent pressure through μ s − μ 0 vap is the solvent reference vapour pressure, and extending the approach of Routh & Russel (1998) to a two-component mixture, (2.16) is substituted into (2.15), givinġ Since ν s (4/3)πR 3 1 , and ν s (4/3)πR 3 2 , as long the film has not yet reached close packing, we obtainĖ ∼ k m ν s ( p 0 vap − p ∞ vap )/kT, which is independent of film composition. In other words, the driving force is dominated by the humidity of the gas. In this one-dimensional drying model, a large surface area of film is assumed, so geometric edge effects are not applicable (Routh 2013). When the film is close packed, Z(φ 1 , φ 2 ) diverges.
This work models drying up until the film is close packed throughout. There are stages of drying beyond this (deformation and aging) that lead to film formation (Sonzogni et al. 2018), but it is the first part of drying, where the particles can move throughout the film, that affects their arrangement in the dried film. If the particles are colloids, as they would be in many common examples of films (Routh 2013), then they are assumed to be colloidally stable. However, the model would also be valid for particles smaller than the colloidal range until the composition at which they precipitate.

Derivation of the spatial derivative of the chemical potential
The chemical potential of component one can be written as a function of φ 1 and φ 2 , such that ∂μ 1 /∂ẑ = f (φ 1 (ẑ), φ 2 (ẑ)), and likewise for the chemical potential of component two, μ 2 . These expressions could be directly inputted into (2.14). However, it is difficult to find expressions for the chemical potential of the solutes that are valid as the solution approaches close packing. One means of resolving this is to relate μ 1 and μ 2 to the chemical potential of the solvent, μ s , for which an expression that diverges at close packing can be found.
For a system which contains n 1 = φ 1 /(4/3)πR 3 1 particles of component one, n 2 = φ 2 /(4/3)πR 3 2 particles of component two and n s = (1 − φ 1 − φ 2 )/ν s particles of solvent per unit volume, conservation of volume can be expressed as 4 3 πR 3 1 n 1 + 4 3 πR 3 2 n 2 + ν s n s = 1. (2.18) The Gibbs-Duhem equation relates the chemical potentials in the system at constant temperature and pressure as n 1 ∇μ 1 + n 2 ∇μ 2 + n s ∇μ s = 0. (2.19) The Gibbs-Duhem equation (2.19), in one dimension, can be rearranged, giving and (2.21) Using the spatial derivative of (2.16), plus the relationship R 1 /R 2 = Pe 1 /Pe 2 , (2.20) and (2.21) become and The purpose of rewriting the chemical potential gradients as in (2.22) and (2.23) is to obtain the advantages, discussed below, of expressing them as functions of ∂μ s /∂ẑ and (∂μ 1 /∂ẑ)/(∂μ 2 /∂ẑ). The chosen model formulation intends for the effects of concentrated solution to be addressed via the divergence in solvent chemical potential, μ s , and the particle interactions to be input via (∂μ 1 /∂ẑ)/(∂μ 2 /∂ẑ). Hence, this extends approaches that are valid only for more dilute solution, for example, Zhou et al. (2017). Since the ratio of the chemical potential gradients is determined physically by the inter-particle interactions, the chemistry of the types of particles that are used to form the film, such as their surface charge, will be relevant (Atmuri et al. 2012).
Note that the Gibbs-Duhem equation (2.19) relates μ 1 , μ 2 and μ s . Equation (2.16) gives an expression for μ s which diverges at close packing. Hence, only one of μ 1 and μ 2 , or a relationship between them, can also be specified. Since (2.20) and (2.21) give ∂μ 1 /∂ẑ and ∂μ 2 /∂ẑ only as implicit functions to be inserted into the conservation equation, different rearrangements of the Gibbs-Duhem equation could have been chosen instead. The particular rearrangements in (2.20) and (2.21) were chosen since their right-hand sides depend on ∂μ s /∂ẑ and the ratio (∂μ 1 /∂ẑ)/(∂μ 2 /∂ẑ). The solvent chemical potential is measurable up to close packing, where it diverges. This means that ∂μ 1 /∂ẑ and ∂μ 2 /∂ẑ also diverge at close packing, in an unknown fashion. An expression chosen for the ratio (∂μ 1 /∂ẑ)/(∂μ 2 /∂ẑ) is more likely to be valid approaching close packing than expressions that could be put forward for each particle separately. Equations (2.20) and (2.21) are also of the same form as each other, as desired.
Equation (2.16) for the osmotic pressure is used to generate an expression for μ s only, rather than being used as an equation of state to also generate the other chemical potentials. The thermodynamic consistency of this approach with respect to the Maxwell relations is further explained in § S6 of the supplementary material (SI) available at https://doi.org/10. 1017/jfm.2021.800.

Diffusiophoresis
The governing conservation equation (2.7) requires an expression for the total flux. In this section, the aim is to investigate the effect of diffusiophoresis. The total flux is written as the sum of two terms, one described as the diffusion term, and the other described as the diffusiophoresis term.
The deviation of osmotic pressure from ideality, and hence Z(φ 1 , φ 2 ) from unity, accounts for the total excluded volume (Russel et al. 1989). However, this model uses the same form of Z(φ 1 , φ 2 ) for the diffusion term, for both the flux of component one and two. In not differentiating between the expressions used for the two fluxes, this form of equation will not give rise to diffusiophoresis. It will give a diffusional model in the sense of the fluxes of both particle types remaining inversely proportional to their particle radius. The flux contribution that is due to diffusiophoresis can be found from first principles for a dilute solution, and written as a separate term. The language of this paper refers to this as the diffusiophoresis term.
Excluded volume effects are more general than diffusiophoresis resulting from excluded volume, since the latter specifically refers to one component moving in response to the concentration gradient of another, whereas in general a component can have volume excluded by both itself and other components. An alternative approach could write the total flux as one term, with a form of Z(φ 1 , φ 2 ) from an equation of state that would include the diffusiophoretic flux. However, the approach taken here allows the effect of diffusiophoresis to be identified separately.
With this approach, a model including diffusiophoresis is constructed, assuming R 2 R 1 . The particles of type one are modelled using an extension of the Asakura-Oosawa model to concentrated solution. Under the Asakura-Oosawa model, the type one particles are assumed to be excluded from a layer of solvent of thickness R DP around each particle of type two. Using this model, in a dilute solution, the drift velocity of the large particles due to diffusiophoresis, U P , in a static film is given by where Γ 1 is the diffusiophoretic drift coefficient (Anderson & Prieve 1984;Sear & Warren 2017). Following the approach of Marbach, Yoshida & Bocquet (2017), this is extended to concentrated solution via invoking the osmotic pressure (SI, § S1.2), to obtain This is taken to be the extra slip velocity, in addition to that due to diffusion, between the particles of component two and the solvent. For the case of hard spheres, which is considered here, R DP = R 1 . Since U P gives the relative diffusiophoretic velocity between type two particles and the solvent, this needs to be converted to a velocity relative to the volume average velocity, U v . Note that the result in (2.24) is the correct diffusiophoretic drift velocity for dilute solution, expected to agree with explicit solvent models. This contrasts with the result of implicit solvent methods which neglect the solvent dynamics, and consequently may overpredict the diffusiophoretic velocity (Sear & Warren 2017).
The velocity of each type of particle i is split up into a component due to diffusion (denoted by the superscript D), which has been found in § 2.2, and a component due to diffusiophoresis (denoted by the superscript P) The local average velocities of the type one, type two and solvent particles are denoted by U 1 , U 2 and U s , respectively. It is necessary to obtain relationships between U P 1 , U P 2 and U P s . First, note that for a static film, U v = 0, so (2.27) The components of the velocities due to diffusion already obey (2.28) Hence, subtracting equation (2.28) from (2.27), provides The velocity of component two and the velocity of the solvent are related by Under the assumptions of this model, diffusiophoresis does not affect the velocity of component one relative to the solvent, This is the case in an idealised model and is therefore used as a starting point for a non-ideal model, with finite size ratios. Although deviation from this will be expected when the particle size ratio is finite (Howard & Nikoubashman 2020), it would not be expected to be large. Solving (2.29)-(2.31) simultaneously, and including sedimentation coefficients, K P1 (φ 1 , φ 2 ) and K P2 (φ 1 , φ 2 ), obtains and To satisfy continuity, it is required that Note that this is not a general result; it is just a result of how this model chooses to write the total flux as two separate terms. Addressing the more general question of whether K 12 (φ 1 , φ 2 ) = K 21 (φ 1 , φ 2 ), these are not equal, as can be seen from analytical expressions for dilute solution (Batchelor 1983). However, equating these expressions would be a reasonable approximation for generating example solutions. K ij (φ 1 , φ 2 ) and K P (φ 1 , φ 2 ) are both denoted with 'K' since they are both sedimentation coefficients, i.e. they represent the same physical phenomenon of hydrodynamic hindrance. The different subscripts identify the terms for which these are the coefficients. Equation (2.32) can be interpreted as showing that component one and the solvent move backwards compared with the diffusiophoretic motion of component two, in order to conserve volume. Using (2.12), which gives U D 1 , and the analogous expression for 928 A15-12 U D 2 , the complete expressions for the component velocities including both diffusion and diffusiophoresis are and (2.35) Including diffusiophoresis, the conservation equation for particles of type one is now and that for particles of type two becomes (2.37) Note that the divergence of the total flux should reflect the total excluded volume, which is incorporated into Z(φ 1 , φ 2 ), hence all the flux terms, including the diffusiophoresis terms, should diverge in the same form as Z(φ 1 , φ 2 ). Since all of the flux terms contain a factor of ∇μ 1 or ∇μ 2 , which diverge via Z(φ 1 , φ 2 ) through the formulations for (2.22) and (2.23), (2.36) and (2.37) obtain suitable divergence at close packing.

Methodology
3.1. Scaling Equations (2.22) and (2.23), the spatial derivatives of the chemical potentials, are substituted into (2.36), and scaled in one dimension, to form the conservation equation for component one, A key observation is that the only dimensionless groups appearing in (3.1) and (3.2) are Pe 1 and Pe 2 i.e. no further dimensionless groups appear due to diffusiophoresis, as would be expected from dimensional analysis.

Coordinate transform
In order to create a static top boundary, the conservation equations are transformed using ξ =ẑ/(1 −t) and τ =t. This results in and This is the coordinate transform that is classically used in solving problems with this geometry, as in Routh & Zimmerman (2004), Trueman et al. (2012b) and Howard et al. (2017b). The boundary conditions are no flux of particles at both the top and bottom boundaries (see Section S1.1 of the SI). For the top boundary condition, the evaporation rate is assumed to be constant, as explained in § 2.2. Setting K P (φ 1 , φ 2 ) = 0 in (3.3)-(3.4) gives a system of PDEs to describe a diffusion-only model that neglects diffusiophoresis. This can be compared with the diffusion-diffusiophoresis model, with non-zero K P (φ 1 , φ 2 ), to observe the effect of diffusiophoresis. Discussion regarding satisfying the Onsager reciprocal relations with each model is provided in the SI, § S3.3.

Numerical method
The system of (3.3)-(3.4) can be solved if φ 1,t=0 , φ 2,t=0 , Pe 1 , Pe 2 and the maximum volume fraction, φ m , which could be a function of φ 1 and φ 2 , are specified. Appropriate forms of K(φ 1 , φ 2 ) and Z(φ 1 , φ 2 ) can be taken from the literature. Also needing to be specified is the ratio of the chemical potential gradients, as explained in § 2.3.
The compressibility is taken to be Z(φ 1 , φ 2 ) = φ m (φ m − φ 1 − φ 2 ) −1 (Trueman et al. 2012b). The canonical concentrated extension of the dilute form of K(φ) for the one-component case, 1 − 6.55φ, is (1 − φ) 6.55 (Russel et al. 1989). For the two-component case, this is extended to K ij|i,j=1 or 2 (φ 1 , φ 2 ) = (1 − φ 1 − φ 2 ) 6.55 (Trueman et al. 2012b). This is adopted for the example results in this work. The particle concentration up to which this is applicable is discussed in § S3.1 of the SI. It is shown in § S3.2 of the SI that the stratification results are minimally affected by how the form of K is extended.
A code was written in MATLAB to solve the PDEs numerically. The governing PDEs are of a similar form to diffusion-advection equations, where ξ/(1 − τ ) is the effective advection velocity, but the effective diffusion coefficient is a complicated function of φ 1 and φ 2 . The PDEs are solved numerically using a finite volume method for the transformed diffusion term. This conserves mass by construction and is an appropriate choice for results containing sharp fronts. Due to the ξ -dependence of its coefficient, the quasi-convective term cannot be formulated using the finite volume method, so finite differences are employed instead. Since the ∂φ/∂ξ term represents advection due to the stretching of the coordinate frame, this term is coded using backward finite differences. This ensures that information travels from upstream to downstream, avoiding potential stability issues that might result from using forward or central finite differences instead. The Euler method is used for time-stepping. The numerical results are verified by checking the final mass balance for each component, and by comparison with asymptotic solutions for high Pe ( § 5). Figure 3 shows the diffusion-only model results for Pe 2 /Pe 1 = 2, as the geometric mean of the Péclet numbers is increased from less than one, to one, to greater than one. The initial volume fractions are φ 1,t=0 = φ 2,t=0 = 0.10, and the maximum packing fraction is taken as φ m = 0.64, the value for randomly packed monodisperse spheres. Any drying film, regardless of its initial volume fraction, will go through higher concentrations in the course of drying. It is most interesting to show examples where the film is initially not too concentrated, since then there is time for the particles to rearrange themselves before movement becomes too hydrodynamically hindered. Hence why φ 1,t=0 = φ 2,t=0 = 0.10 is chosen for these examples. Drying an initially more concentrated film involves the same phenomena as starting from a more dilute film: what changes (when diffusiophoresis is included) is the proportions of the drying time when different phenomena are dominant (see § S4 in the SI). § S5 in the SI demonstrates that there is minimal effect on the stratification results as φ m is varied. From a mass balance on the particles, the end of drying would be at τ f =t f = 0.6875.

Diffusion
The model allows input of any form of K ij . To run an example with numerical ease, and due to not being able to analytically derive the cross-terms for concentrated solution, K ij|i / = j (φ 1 , φ 2 ) = 0 is used. The same model could be rerun with K 12 (φ 1 , φ 2 ), .55 would be expected to give similar results. This is because the fluxes for the cross-terms are small.
In addition to modifying for hydrodynamic hindrance via K(φ 1 , φ 2 ), the diffusion-only model, as run in these examples, also differs from that of Zhou et al. (2017) in the expressions used for the chemical potential: (i) not including the cross-interaction terms in the chemical potential expressions; and (ii) writing the chemical potentials via the osmotic pressure so that they are valid for concentrated solution. For the ratio of the chemical potential gradients, (∂μ 1 /∂ẑ)/(∂μ 2 /∂ẑ), a non-interacting expression is chosen: (∂(ln φ 1 )/∂ẑ)/(∂(ln φ 2 )/∂ẑ). Whilst this expression is chosen for these example numerical results, the model is general, and can be run with input of any form of (∂μ 1 /∂ẑ)/(∂μ 2 /∂ẑ), for example, from different equations of state (EoS), including the Boublík-Mansoori-Carnahan-Starling-Leland EoS for hard sphere mixtures (Heyes & Santos 2016).
The parameter values used are collected in table 1. To give an idea of a physical system that could be represented with these Péclet numbers, Liu et al. (2018) estimated evaporation rates ofĖ = 2.0 × 10 −8 m s −1 at T = 296 K and 60 % humidity. For a film of initial height H = 1.0 mm, and a solvent of viscosity η = 0.001 Pa s, particles of size {7.5, 15} nm would correspond to Péclet numbers of {0.7, 1.4}. Silica particles, for example, of this size can be readily obtained.
The temporal and spatial resolutions, in (ξ =ẑ/(1 −t), τ =t) coordinates, are selected in order to give mass balance errors of < 1 % in each component by near the end of drying. For example, figure 3(a) used ξ = 0.01 and τ = 10 −8 , giving a mass balance error of 0.17% in φ 1 and 0.42% in φ 2 at τ = 0.683. Increasing Pe requires smaller ξ in order to resolve the increasing sharpness of the curves, and correspondingly smaller τ for stability. It is checked that the change on further improving the resolution is acceptably small. For example, halving the ξ used for figure 3(a) roughly halves the mass balance errors to 0.08% in φ 1 and 0.21% in φ 2 . Halving the τ used in figure 3(a) produces a negligible change in the mass balance errors. The resolutions used to generate each figure are tabulated in § S1.3 of the SI.  Large-on-top stratification is predicted, shown by the blue lines being above the red lines in the upper part of the film, as expected for a diffusion-only model; (b) Pe 1 = 0.70 and Pe 2 = 1.40. As the Péclet numbers now straddle one, there is a greater difference between the volume fractions of the two components at the top surface, compared with (a); (c) Pe 1 = 2.8 and Pe 2 = 5.6. Péclet numbers greater than one cause there to be a sharp transition in volume fraction between the lower and upper parts of the film. (a) Pe 1 , Pe 2 < 1(Pe = 6πηRĖH/kT), (b) Pe 1 < 1, Pe 2 > 1 and (c) Pe 1 , Pe 2 > 1.

Diffusion only
In figure 3, the films become more concentrated over time. The volume fraction profiles show the large particles stratifying to the top surface, as expected for a purely diffusional  In figure 3(a), since Pe 1 , Pe 2 < 1, the film remains fairly uniform over the course of drying, but with the film concentration increasing towards the top surface. As √ Pe 1 Pe 2 is increased from 0.25 (figure 3a) to 1.0 (figure 3b), the degree of stratification increases, and the concentration profiles sharpen. Increasing √ Pe 1 Pe 2 again from 1.0 (figure 3b) to 4.0 (figure 3c), further sharpens the profiles, starting to form a sharp transition between the volume fractions at the top and bottom of the film. This sharpening is due to slower diffusion relative to the evaporation rate, such that the lower parts of the film do not 'see' the effect of the evaporating top surface.

Diffusiophoresis
The same sets of parameters are run as in figure 3, but now with the addition of the diffusiophoresis term. The form used for K P (φ 1 , φ 2 ) is taken as (1 − φ 1 − φ 2 ) 6.55 , which is the same as for K ii (φ 1 , φ 2 ). The exact functional forms used for K ii (φ 1 , φ 2 ) and K P (φ 1 , φ 2 ) do not significantly impact the numerical results. Interestingly, if K ii (φ 1 , φ 2 ) and K P (φ 1 , φ 2 ) are equal, then the functional form does not affect the relative fluxes, which is the important factor when determining the importance of diffusiophoresis, relative to diffusion.
The resulting volume fraction profiles are shown in figure 4.

Non-enhanced diffusiophoresis
In contrast to the diffusion-only results, at all Pe modelled, including diffusiophoresis results in the film being almost uniform in composition. Byt = 0.60, the films are now only slightly large-on-top stratified, with the exception of the film with low Pe (figure 4a), which has just transitioned from slightly large-on-top stratified to slightly small-on-top stratified. As is discussed in § S4 of the SI, the films in figures 4(b) and 4(c) will also transition to slightly small-on-top stratified even later than the latest time shown on these graphs. These results therefore suggest that diffusiophoresis acts against large-on-top stratification, such that the effects of the diffusion and diffusiophoresis terms largely counteract each other. As in the diffusion-only case, increasing √ Pe 1 Pe 2 results in sharper profiles, as shown in figures 4(b) and 4(c). There is little difference between the red and the blue curves throughout drying, indicating that the film is nearly uniform in composition; (b) Pe 1 = 0.7 and Pe 2 = 1.4. The film profiles are sharper than in (a), due to the higher Péclet numbers, but there is still little stratification between the two components; (c) Pe 1 = 2.8 and Pe 2 = 5.6. As in (b), little stratification between the two components develops, although the film profiles are sharper than (b) due to the Péclet numbers being increased again. (a) Pe 1 , Pe 2 < 1(Pe = 6πηRĖH/kT), (b) Pe 1 < 1, Pe 2 > 1 and (c) Pe 1 , Pe 2 > 1.

Increased diffusiophoresis
To verify that it is diffusiophoresis in the model which is promoting small-on-top stratification, the model is run with an increased diffusiophoretic drift coefficient. From (2.25), this corresponds to a larger value of R DP , i.e. a larger excluded volume.
For example, non-spherical particles, such as polymers, may be excluded from a distance equal to their radius of their gyration. The radius R 1 in this case would correspond to the radius of a sphere with the equivalent volume. The ratio R DP /R 1 would therefore be expected to be greater than unity for a non-spherical particle. However, in order to accurately model non-spherical particles, the other parameters, K ij (φ 1 , φ 2 ) and Z(φ 1 , φ 2 ), would also need adapting. As merely an example of increasing diffusiophoresis alone, results for the case R DP /R 1 = (4π/3) 1/2 are shown in figure 5.
At small Pe (figure 5a), the stratification increases over time, as small particles accumulate at the top surface. As the Pe numbers are increased (figure 5b), the stratification develops at a greater rate, until the top surface is almost entirely small particles. At large Pe numbers (Pe 1 , Pe 2 > 1, figure 5c), a thin layer of small particles at the top surface rapidly develops at early drying times. The thickness of this layer grows over time. Directly beneath this layer, φ 1 falls sharply to a much smaller value (∼0.15), as φ 2 , at ∼0.49, dominates this region of this film. Over the lower region of the film, where the effect of evaporation from the top surface is yet to significantly manifest, the volume fractions of both components fall to their initial values.
Continuing this present model to the end of drying predicts a film of two regions: an upper layer of almost entirely small particles, at an approximately constant concentration, and a lower layer with a majority of large particles, at an approximately constant concentration. A drawback of these results is that the form of K ii used allows particle rearrangement within the close-packed region (SI, § S3.1). This, numerically, allows the layer of small particles to grow into the close-packed region beneath it. However, physically, such a degree of rearrangement is unlikely.

Key mechanism
The results can be understood by considering the scaling of (3.3) and (3.4). The diffusiophoretic flux of component two has an extra order of (φ 1 , φ 2 ) in its coefficient, compared with the diffusive flux of component two, which is only just compensated for by the prefactors (Pe 2 /Pe 1 ) and 9/4. Hence the results for hard spheres with diffusiophoresis in figure 4 show minimal stratification. Noting this concentration dependence of the relative importance of diffusiophoresis, in the SI, § S4, the effect of varying the initial concentration on the film structure is presented.
For figure 5, the ratio R DP /R 1 can be considered as a geometric factor, relating the excluded distance of the small particles around the large particles, to the radius of the small particles. With the value of R DP /R 1 adopted for this figure, the extra order of (φ 1 , φ 2 ) in the coefficients of the diffusiophoretic flux of component two is more than compensated for by its other coefficients. This is sufficient to make the magnitude of the total diffusive and diffusiophoretic flux of component two greater than that of component one, resulting in small-on-top stratification. Note that the scaling of the ratio of the diffusiophoretic flux to the diffusive flux of component two increases with φ 2 . Hence just beneath the front of small particles, where φ 2 is at its largest value in the film, there is a strong driving force for more small particles to travel to the top surface. This is seen most clearly in figure 5(c).
To illustrate the relative magnitudes of the different fluxes across the film, figure 6 is plotted. The fluxes are plotted againstẑ for Pe 1 = 0.7 and Pe 2 = 1.4, so the flux profiles can be compared against the concentration profiles given in figures 4(b) and 5(b). An early time,t = 0.17, when the profiles are forming, is chosen for this illustration. The cases with non-enhanced and enhanced diffusiophoresis are shown for comparison. In figure 6(a), the non-enhanced case, the total fluxes for components one and two are similar throughout the  as seen by the increasing difference between the volume fractions of the two components at the top surface ast increases; (b) Pe 1 = 0.7 and Pe 2 = 1.4. Small-on-top stratification develops more rapidly than in (a), as can be seen by the greater distance between the red and blue lines at the top surface. This is due to the increase in the Péclet numbers; (c) Pe 1 = 2.8 and Pe 2 = 5.6. Significant stratification between the two components develops rapidly, as can be seen from the large volume fraction difference at the top surface att = 0.17. The depth of the layer enriched in small particles at the top surface grows over time. (a) Pe 1 , Pe 2 < 1(Pe = 6πηRĖH/kT), (b) Pe 1 < 1, Pe 2 > 1 and (c) Pe 1 , Pe 2 > 1. film, hence a nearly uniform film results in figure 4(b). The total flux for each component is in the negativeẑ-direction, consistent with the solvent evaporating from top surface. The magnitude of the diffusive flux of component one (the smaller particles) is larger than that of component two, despite the gradients being similar, since this is proportional to 1/Pe. In figure 6(b), at the uppermost part of the film, which corresponds to the small particles-enriched section in figure 5(b), the magnitude of the total flux of component two is less than that of component one. These total fluxes are both in the negativeẑ-direction. However, the mechanism that causes the layer of small particles to continue to build, occurs below this layer, where the magnitude of the total flux of component two is greater

Scenario comparison
For ease of comparison, the results from figures 3-5 are summarised in figure 7. The degree of stratification, defined by (4.1), is plotted against the geometric mean of the component Péclet numbers at selected times throughout drying, for the diffusion-only, diffusiophoresis and enhanced diffusiophoresis cases.
There are different possible stratification measures that can be used, but the one below, similar to the average separation measure used by Tang et al. (2018), is deemed to be appropriate here From the summary figure, it can be clearly seen that the diffusion-only model leads to large-on-top stratification, whilst adding diffusiophoresis promotes β in the direction of small-on-top. In the diffusion-only and enhanced diffusiophoresis cases, the greatest magnitudes of β are reached as the Péclet numbers straddle one, and β increases with time. The non-enhanced diffusiophoresis results remain negligibly stratified throughout.

Effect of cross-interactions
Similarly to figure 5, Howard et al. (2017a) obtain a growing layer of small particles, although they did not run both Pe 1 , Pe 2 < 1. The particle size ratios explored by Howard et al., Pe 2 /Pe 1 = {4, 6, 8}, were larger than the Pe 2 /Pe 1 = 2 studied here, but based on scaling arguments, similar trends would be expected. Howard et al. found faster rate of growth of the layer of small particles at higher Pe 2 /Pe 1 , and that would also be found from this model. Fortini et al. (2016) obtain a growing layer of small particles, although they only present simulations for Pe 2 /Pe 1 = 7. The large Péclet numbers tested may explain why they predicted such a sharp stratification.
Whilst the work of Howard et al. and Fortini et al. give similar small-on-top stratification to figure 5 with enhanced diffusiophoresis, they differ from the approximately uniform films predicted in figure 4, for hard spheres. The work of Zhou et al. (2017) demonstrates that cross-interaction terms alone, in the absence of diffusiophoresis, can also give rise to small-on-top stratification. Note also that Howard et al. (2017a,b) use models accurate to higher order, and in doing so account for cross-interaction terms. Since their models combine these EoS with diffusiophoresis, it suggests that both diffusiophoresis and cross-interaction terms contribute to small-on-top stratification. From figure 4, which did not include cross-interaction terms, in order to see the effect of diffusiophoresis separately, diffusiophoresis alone seems not to be sufficient to give rise to significant small-on-top stratification. Figure 5 suggests that in cases with enhanced volume exclusion, diffusiophoresis can cause small-on-top stratification in the absence of cross-interactions. Cross-interaction terms could be incorporated into this model by using (1.1) to form an expression for (∂μ 1 /∂ẑ)/(∂μ 2 /∂ẑ), to replace the example in table 1 which only includes entropy.
Further comparing these results with the work of Zhou et al. (2017), the term in the Zhou et al. model which is responsible for the small-on-top stratification is the cross-interaction term; diffusiophoresis could also be thought of as a cross-interaction term. Zhou et al. argue that the flux of large particles due to the cross-interaction term must be greater than the flux due to diffusion of the large particles down their own concentration gradient, in order to observe small-on-top stratification. A more general explanation of the condition, applicable to this diffusiophoresis model also, would be that the cross-interaction term must be sufficiently large to make the total flux of the large particles exceed that of the small particles. Zhou et al.'s cross-interaction term disproportionately affected the flux of the larger particles, and that is also the case with the diffusiophoresis term in this model.
One similarity between this model and the Langevin dynamics simulations of Fortini et al. (2016) and Howard et al. (2017a) is that diffusiophoresis results from a short-range repulsive interaction between the two particle types, whilst Fortini et al. use the Yukawa interaction, and Howard et al. adopt the Weeks-Chandler-Andersen potential. As explained in § 2.4, an alternative approach to adding diffusiophoresis to the diffusion model in (2.13) would be to include the exclusion potential in the chemical potential term. 4.4.2. Finite large particle R 2 size effects Note also that the diffusiophoresis term is idealised, assuming infinite size ratio between the two particle sizes, yet is it being used to model a finite size ratio (Sear & Warren 2017). The results in figure 4 therefore give an upper bound for the predicted effect of diffusiophoresis. This is because the surface of a component two particle of finite radius will not all be able to get as close to the surface of a component one particle, compared with the case of infinite-sized component two particles. Whilst this upper bound will provide insight into the importance of diffusiophoresis, and it is likely to be most important when R 2 R 1 , an improved model would include the R 2 -dependence for the case where R 2 is finite. In this work, Pe 2 /Pe 1 = 2 was run as an example. The model can be run with a larger size ratio, such as Pe 2 /Pe 1 = 10, in which case the infinite size ratio assumption becomes more accurate, as long as the dispersion remains colloidally stable.
There are various approaches to correcting the model for when R 2 is finite. The simplest correction is that by Anderson, Lowell & Prieve (1982), which still models the small particles as a solute in a fluid but corrects the geometry to include the finite size of R 2 . Including the correction term to first order in Pe 1 /Pe 2 results in the diffusiophoretic velocity being 2/3 of that given in (2.24). This suggests the order of magnitude of the error that is introduced from using the Asakura-Oosawa model. Brady (2011) derives equations for a colloidal system that no longer requires the small particles to be much smaller than the large ones, showing that the continuum and colloidal approaches agree in the limit where this requirement is satisfied. Whilst being valid for any particle size ratio, Brady notes that the approach is analytically limited to dilute solution. Marbach, Yoshida & Bocquet (2020) derive a similar expression to Brady, by considering Stokes flow around a large particle and the force-free particle condition for diffusiophoretic motion. In addition, for the case of finite particle size ratio, the present model could be improved to accurately model the remaining excluded volume/diffusiophoresis effects. Besides small particles being excluded from large particles (small-large), there are also large-large, small-small and large-small exclusions to consider.

Other model comments
Other possible improvements to the model include making φ m a function of Pe 2 /Pe 1 , φ 1 and φ 2 , to model how particles of two different sizes pack. The deviation from φ m = 0.64 will be more significant as Pe 2 /Pe 1 is increased, but assuming a constant value should be a reasonable approximation for Pe 2 /Pe 1 = 2, as in figures 3-5.
The dilute limit of the equations can be inferred by construction of the model: the leading order of the diffusion velocities is just −(kT/6πηR i )K ii ∇ ln φ i , and the leading order of the diffusiophoretic drift velocity for the large particles is −(3φ 1 kT/8πηR 1 )K P ∇ ln φ 1 .

Asymptotic solution for large Pe
As Pe 1 and Pe 2 are increased, the shapes of the curves appear to attain a sharp transition front with constant shape. This suggests that an asymptotic solution could be sought in the regime of high Pe 1 and Pe 2 , with Pe 2 > Pe 1 . Physically, this corresponds to a film dried with a fast evaporation rate, such that the diffusive processes do not have time to flatten the shape of the transition front. The method of Russel et al. (1989), for one particle type, also adopted by Trueman et al. (2012b) for two particle types, is followed.
It is expected that the total particle volume fraction at the top of the film will be φ m , as fast evaporation will lead to this being quickly obtained, whilst the volume fractions at the bottom of the film will be φ 1 (τ =t = 0) and φ 2 (τ = 0). As the Péclet numbers are increased, there will be a sharper transition between these values. The position of this transition is denoted as P(τ ).
The method is outlined here, with further details given in § S2 of the SI. The general strategy is to use the change of variable X = [ξ − P(τ )]Pe 2 to expand the region around the transition. In the case of the diffusion-only model, there is only one transition front, so the solution method is that found in the supplementary material of Trueman et al. (2012b), but with the governing hydrodynamics equations corrected. In the case of the diffusion-diffusiophoresis model with strong diffusiophoresis, a second transition appears, which complicates the solution, as will be outlined in § 5.2.
To leading order in Pe 2 , the PDEs in (3.3)-(3.4) reduce to ordinary differential equations (ODEs) for dφ 1 /dX and dφ 2 /dX, describing the shape of the transition front. The position of the front P(τ ) is found from integrating the ODEs across the discontinuity from X = −∞ to X = ∞ and summing This is the same result as for the one component case, and as would be obtained from a mass balance across the drying film, assuming a sharp discontinuity at P(τ ). Interestingly, it was necessary to assume large Pe 1 as well as large Pe 2 , for the assumption of a sharp discontinuity in φ 1 at P(τ ) to be valid.
The resulting asymptotic solution for an example set of Péclet numbers is shown in figure 8. In the transition region, there is excellent agreement between the asymptotic and numerical solutions. This validates the numerical code used to generate figures 3-5. When φ 1,τ =0 = φ 2,τ =0 , the asymptotic solution predicts that φ 1,X→∞ = φ 2,X→∞ (SI). Whilst there is therefore disagreement in the values at the top surface, these do tend to the values given by the asymptotic solution over time. The small excess of large particles at the top surface, a diffusive effect, counterbalances the excess of small particles in the rest of the film (see (5.2)), to satisfy the mass balance.

Diffusiophoresis
As for the diffusion-only results, the diffusiophoresis results at high Pe 1 and Pe 2 appear to be curves of constant shape moving through the film. The numerical code is verified using the equations with enhanced diffusiophoresis. This is because the main purpose of this asymptotic solution is to specifically verify the coding of the diffusiophoresis term. The approach to deriving the asymptotic solution for the non-enhanced diffusiophoresis case is very similar to that described in § 5.1. However, there are two key transitions in the curves for the enhanced diffusiophoresis case: The first is the position of the particle front i.e. the front of total volume fraction. This is denoted by P(τ ), as in the diffusion-only case. The second is between the region of constant high φ 1 and low φ 2 values and the region of constant high φ 2 and low φ 1 . The position of this transition is denoted by Q(τ ). The transitions are annotated on figure 9. This section seeks to construct an asymptotic solution for all parts of the film: is used to expand around this transition. It is assumed that Q(τ ) takes the form which is a front travelling at constant speed v in (ẑ,t) coordinates. The speed v at which Q(τ ) travels can be found by integrating the equation for dφ 1 /dY or dφ 2 /dY across the transition (SI, S2.2.6). This is equivalent to a mass balance on φ 1 or φ 2 across Q(τ ). Physically, both the diffusion and diffusiophoretic fluxes work to propel this front, changing (φ 1 , φ 2 ) from the values on one side of the front to those on the other. Note that it is expected that v < φ m /[φ m − (φ 1,τ =0 + φ 2,τ =0 )], since Q(τ ) is located higher in the film than P(τ ). The ODEs for dφ 1 /dY and dφ 2 /dY can be integrated, and are combined for the case where K 12 (φ 1 , φ 2 ) = K 21 (φ 1 , φ 2 ) = 0, K 11 (φ 1 , φ 2 ) = K 22 (φ 1 , φ 2 ) = K P (φ 1 , φ 2 ) and the particle chemical potentials are purely entropic, dμ 1 /dμ 2 = d ln φ 1 /d ln φ 2 . In the upper region of the film, φ 1 + φ 2 ≈ φ m . Substituting in φ 2 = φ m − φ 1 gives a quadratic equation for the values of φ 1 either side of Q(τ ).
There is reasonably good agreement between the numerical and asymptotic solutions, although the asymptotic solution slightly overpredicts the sharpness of the transition around P(τ ). Using the numerical results for φ 1 to reconstruct the φ 2 curves with the asymptotic solution, and vice versa, gives excellent agreement. Therefore, the small discrepancy that does arise is mostly from the step of finding X(φ 1 ).
At fixed values of φ 1,τ =0 , φ 2,τ =0 and φ m , the asymptotic solution for v shows minimal variation with Pe 1 and Pe 2 . This allows a schematic diagram to be drawn for the regions  Figure 10. Schematic of the regions of the asymptotic solution for large Pe, when φ 1,t=0 = φ 2,t=0 = 0.10 and φ m = 0.64. The unreachable region is shaded in grey.
in the film, as exemplified in figure 10 for the set of φ 1,t=0 , φ 2,t=0 and φ m values that were used in figure 9. Aside from verifying the numerical code, the asymptotic solutions shed insight into the large Pe regime. Remembering the definition of the Péclet number, Pe = 6πηRĖH/kT, this corresponds to evaporation being fast compared with diffusion. Highly stratified films are predicted in this regime, with the asymptotic solutions allowing prediction of the composition of each region of the film. Figure 10 shows the depths of both (i) the region between the top surface and Q(τ ), and (ii) the region between Q(τ ) and P(τ ), increasing linearly witht. The depth of the region between the top surface and Q(τ ) is the predicted thickness of the layer of small particles at the top surface. It can be seen in figure 10 that this small particle-rich layer is thinner than the large particle-rich layer beneath it, as a result of the relative speeds at which the top surface, Q(τ ) and P(τ ) recede. These fronts recede at speeds of 1, v and φ m /[φ m − (φ 1,τ =0 + φ 2,τ =0 )], respectively, in (ẑ,t) coordinates. The speeds are presented as the magnitudes of the gradients of the diagonal lines separating the regions in figure 10. Hence the schematic diagram may provide an explanation for why a thin layer of small particles can form at the top surface of a drying film.

Conclusions
PDE models were successfully derived to describe the volume fraction evolution of two components in a drying film, for both diffusion-only and diffusion-diffusiophoresis cases. The diffusion term is written in a general form, allowing input of sedimentation coefficients and the compressibility, and can be used with different chemical potential expressions.
Both the scaling and the numerical results of the diffusion-only model predict large-on-top stratification. This model was run without interactions, but these could be included in future work. In contrast, the numerical results from the diffusion-diffusiophoresis model demonstrate that diffusiophoresis is a feasible explanation of the experimental observations of a thin layer of small particles at the top surface of a dried film, but it is dependent on the size of the excluded volume. The results suggest that diffusiophoresis combined with another effect, such as cross-interactions, may be necessary. Whether it is diffusiophoresis causing this will therefore need to be experimentally verified. Diffusiophoresis is typically investigated by setting up a gradient of salt and observing the motion of colloidal particles along this gradient, for example in dead-end microfluidic channels (Shin et al. 2016;Singh et al. 2020). This could be adapted with two colloidal particle sizes. In addition, the idealised diffusiophoresis term in this model will need to be adapted to predict whether diffusiophoresis is still expected to be significant with finite R 2 . Nevertheless, this simple PDE model, that allows the diffusiophoresis term to be switched on and off, gives clear insight into what aspect of more complex models -namely, a repulsive cross-interaction, with short range being sufficient -leads to their predictions of small-on-top stratification.
In the enhanced diffusiophoresis model, the layer of small particles at the top surface is predicted to grow over time, similarly to the results of Langevin dynamics simulations in the literature (Fortini et al. 2016;Howard et al. 2017a), and at faster rates at higher Pe numbers. This suggests that in order to obtain a top surface film of a desired component, such as a biocide, fast drying rates should be used. Asymptotic solutions were found for this regime. Declaration of interests. The authors report no conflict of interest.