Self-similarity in particle accumulation on the advancing meniscus

Abstract When a mixture of viscous oil and non-colloidal particles displaces air between two parallel plates, the shear-induced migration of particles leads to the gradual accumulation of particles on the advancing oil–air interface. This particle accumulation results in the fingering of an otherwise stable fluid–fluid interface. While previous works have focused on the resultant instability, one unexplored yet striking feature of the experiments is the self-similarity in the concentration profile of the accumulating particles. In this paper, we rationalise this self-similar behaviour by deriving a depth-averaged particle transport equation based on the suspension balance model, following the theoretical framework of Ramachandran (J. Fluid Mech., vol. 734, 2013, pp. 219–252). The solutions to the particle transport equation are shown to be self-similar with slight deviations, and in excellent agreement with experimental observations. Our results demonstrate that the combination of the shear-induced migration, the advancing fluid–fluid interface and Taylor dispersion yield the self-similar and gradual accumulation of particles.


Introduction
Particle-laden flows have been studied extensively both theoretically and experimentally, owing to their relevance in avalanches and mudflows (Savage & Lun 1988;Gray & Ancey 2009;Gray & Kokelaar 2010), three-dimensional printing of complex fluids (Lewis 2006;Roh et al. 2017) and cell migration in biological systems (Vejlens 1938;Zhou & Chang 2005). In particular, there is much interest in understanding hydrodynamic interactions of suspended particles that may lead to non-uniform particle distribution or particle accumulation on a fluid-fluid interface. For instance, Vejlens (1938) experimentally studied blood flowing through capillary tubes and observed a much stronger red † Email address for correspondence: sungyon@umn.edu colour in the region behind an advancing meniscus, which indicates a higher local concentration of red blood cells. More recently, Zhou & Chang (2005) demonstrated that the accumulation of red blood cells at the meniscus causes the penetration failure of blood suspensions into a capillary with a diameter smaller than 100 micrometres, which is a major obstacle in miniaturising blood diagnostic tools. The observation of particle accumulation on the fluid-fluid interface is not limited to blood flows and even extends to geological systems. Bhattacharji & Smith (1964) observed the concentration gradient of minerals in rock-structure fissures, which is due to the accumulation of solid particles.
To rationalise the mineral concentration gradient (Bhattacharji & Smith 1964), Bhattacharji & Savic (1965) hypothesised that the accumulation of solid particles is caused by a recirculation flow near the interface. They developed a mathematical model to predict the velocity field of this recirculation flow in the absence of the particles, later known as a 'fountain flow'. Following their work, Karnis & Mason (1967) experimentally investigated the accumulation of particles behind an advancing meniscus by tracking aluminium trace particles to obtain the fountain flow streamlines. In qualitative agreement with the theoretical prediction (Bhattacharji & Savic 1965), their results indicate that particles entering the fountain flow region at the meniscus move radially away from the axis towards the wall and are directed radially inwards again as they exit the fountain flow region. At the same time, the particles close to the wall are trapped in the slow-moving region and remain at the wall, resulting in the accumulation of particles at the meniscus. In addition, Karnis & Mason (1967) claimed that the particle accumulation must vanish when the wall effects are eliminated.
Despite the apparent success of the fountain flow model, Chapman (1990) proposed that shear-induced migration be an alternative mechanism for particle accumulation, the effect of which does not vanish with wall effects. Shear-induced migration refers to the tendency of particles to move from high-shear to low-shear regions (Leighton & Acrivos 1987;Phillips et al. 1992;Nott & Brady 1994), which has been considered experimentally (Hookham 1986;Abott et al. 1991;Koh, Hookham & Leal 1994) and theoretically (Cook 2008;Ward et al. 2009;Murisic et al. 2011Murisic et al. , 2013Lee, Stokes & Bertozzi 2014;Mavromoustaki & Bertozzi 2014;Wang & Bertozzi 2014;Lee, Wong & Bertozzi 2015;Snook, Butler & Guazzelli 2016). In particular, a number of recent studies have focused on building the theoretical framework for particle transport in the axial direction, while capturing the effects of shear-induced migration normal to the flow direction. Select examples include modelling of particle dispersion in different configurations (Griffiths & Stone 2012;Ramachandran 2013;Christov & Stone 2014), while others considered the effects of particle friction (Lecampion & Garagash 2014) and a yield-stress interstitial fluid (Hormozi & Frigaard 2017).
Then, connecting shear-induced migration to particle accretion, particles migrate towards the channel centreline in pressure-driven flow and assume a higher particle average velocity than that of the fluid upstream of the interface. This velocity differential results in a net flux of particles towards the meniscus and causes particle accumulation (Xu, Kim & Lee 2016;Luo, Chen & Lee 2018). Chapman demonstrated that while the accumulation rate increases with the particle size, there exists a non-zero asymptote when the particle size approaches zero. More recently, Ramachandran & Leighton (2007) supported the mechanism proposed by Chapman and confirmed that particle accumulation is observed even when the wall effects vanish. Therefore, we conjecture that particle accumulation is governed primarily by the shear-induced migration far upstream, while the fountain flow near the meniscus may play a secondary role.
One of the consequences of particle accumulation on the fluid-fluid interface is the emergence of the interfacial instability. For instance, Tang et al. (2000) first observed viscous fingering, when a mixture of particles and viscous fluid displaces air inside a Hele-Shaw cell. This surprising phenomenon, or 'particle-induced viscous fingering', contrasts with the stable case of a pure viscous liquid displacing air. Following their work, Ramachandran & Leighton (2010) observed a similar fingering instability upon squeezing the particle suspension between two parallel plates. More recently, Xu et al. (2016) experimentally characterised fingering patterns for varying particle volume fractions. They also successfully validated the effects of shear-induced migration that leads to the accumulation of particles on the advancing oil-air interface, which was followed by the exploration of the role of channel confinement (Kim, Xu & Lee 2017) and linear stability analysis to predict the critical wavenumber (Hooshanginejad, Druecke & Lee 2019).
Despite the recent development in particle-induced viscous fingering, one of the fundamental phenomena that leads to fingering remains unexplained. While particle accumulation on the interface has been attributed to shear-induced migration, Xu et al. (2016) noted that particles accrete on the advancing interface in such a way that the depth-averaged particle volume fractionφ is independent of time when scaled appropriately. In this paper, we focus on understanding this self-similar behaviour of particle accumulation both experimentally and theoretically. In § 2, we first quantitatively show the self-similar behaviour by defining an inner radius of the suspension, R in , which separates the quasi-steady region far upstream of the interface from the particle accumulation region. In § 3, we mathematically model the suspension flow as a continuum and obtain a depth-averaged particle transport equation, based on the theoretical framework developed by Ramachandran (2013). The numerical solutions to the transport equation are also presented in § 4 and the paper is concluded in § 5.

Materials and methods
We hereby describe the experimental method and data that were previously included in (Xu et al. 2016;Luo et al. 2018). The experiments are performed by radially displacing air with a particle suspension inside a Hele-Shaw cell, as depicted in the schematic in figure 1(a). The cell consists of two plates of plexiglass, each with dimensions of 30.5 × 30.5 × 3.8 cm, that are separated by a gap thickness h with standard shims (McMaster) secured at four corners. The suspension is prepared by mixing PDMS silicone oil (United Chemical, viscosity η l = 0.096 Pa · s, density ρ l = 0.96 g cm −3 ) with fluorescent polyethylene particles (Cospheric, diameter d = 125-150 μm, density ρ p = 1.00 g cm −3 ) inside a syringe to an initial volume fraction of φ 0 . The suspension is left to sit for a sufficient time to allow entrapped air bubbles to escape with minimal particle settling.
A complete list of experimental parameters (i.e. h and φ 0 ) is summarised in table 1. A syringe pump (New Era Pump Inc., model NE-1010) is used to inject the suspension into the Hele-Shaw cell through a hole drilled at the centre of the bottom plate via a clear vinyl tube. For all the experiments discussed herein, the volumetric flow rate Q is kept constant at 150 mL min −1 . A light-emitting diode (LED) panel (EnvirOasis, 75 W, 4200 Lumin) is placed under the cell to provide uniform illumination. All the experiments are recorded with a digital camera (Canon 60D, resolution 1920 × 1080, field of view 64 • ) directly from above, at a rate of 30 frames per second.   Table 1. Experimental parameters of the gap thickness, h and particle concentration, φ 0 . The φ 0 range that is denoted as '# -#' increases by an increment of 1 %.
We process the videos collected from experiments in MATLAB with two specific goals: to track the suspension-air interface and to measure the depth-averaged particle concentration,φ, in each experimental image. To track the interface, we use the built-in edge detector in MATLAB, which takes the smoothed images from median filter as an input. The detector utilises the 'Canny' method (Canny 1986) to identify changes in local intensity gradient to find maxima, which are treated as edges. As indicated in figure 2(c), R(t) corresponds to the radius of the circle fitted from all data points on the suspension-air interface at a given time t.
Next, to extractφ, we first systematically subtract the background image from all the experimental images. This allows us to eliminate any inherent non-uniformity in lighting and to obtain the light intensity matrix I from the processed images. Then, the depth-averaged concentrationφ is correlated to I, based on the following relationship (Grasa & Abanades 2001;Xu et al. 2016): where I min and I max are the minimum and maximum light intensity values of the given image, respectively. The empirical parameter, k, is acquired from mass conservation: φ(r) is integrated from the centre to the suspension-air interface such that 2 R δφ r dr = φ 0 (R 2 − δ 2 ). Note that the lower bound of the integral, δ, is chosen to be approximately 6d in order to account for the injection hole.
Based on the extracted values ofφ(r, t), we calculate the radial particle flux, f (r, t), across a given cylindrical surface inside the suspension. As illustrated in the schematic in figure 1(a), we consider a cylindrical control volume whose radius is given by the arbitrary radial position, r, between the inlet and the interface, R(t). Conservation of particles inside the control volume yields f (r, t) = φ 0 Q − dV p /dt, which balances the rate of particle injection at the inlet with the time rate of change of the total volume of particles inside the control volume, or V p = r δ 2πrφ(r)h dr. Hence, to obtain f (r, t), we compute the change in V p for all r (at an increment of δr) between two adjacent times, t 1/2 = t ∓ δt/2. Here, δr is a small radial distance, while δt corresponds to a small temporal deviation from the reference time t. For the given set of experimental data, we set δr ∼ 3d and δt = 1/30 s, which are limited by the spatial resolution of the image and the video frame rate. Figure 1(b) shows the time-elapsed images from two experimental runs at h = 1.4 mm: φ 0 = 0.2 (a) and φ 0 = 0.35 (b). Interfacial deformations associated with viscous fingering are clearly observed at φ 0 = 0.35, while the suspension-air interface remains stable and circular at φ 0 = 0.2. We further quantify this particle accumulation on the interface by experimentally measuringφ, see figure 2(a), which shows an increase with r at φ 0 = 0.25, where r is the radial coordinate defined from the injection centre. Furthermore, a colour map of the experimental image at t = 3 s in figure 2(c) shows a higher colour intensity close to the interface, corresponding to a local increase inφ at φ 0 = 0.25. Both the plot ofφ and the colour map verify the accumulation of particles on the interface even in the stable regime (i.e. φ 0 < 0.3), which is the focus of our current study.

Experimental observations
Notably, as shown in figure 2(b), allφ(r, t) profiles approximately collapse on to a single curve, when r is normalised by the instantaneous position of the suspension interface R(t). This plot ofφ with r/R provides direct evidence that particles accumulate on the advancing interface in a self-similar manner, despite some deviations in the normalised curve near the interface. To further analyse this, we presently divideφ(r/R(t)) into three distinct regions (I, II and III), as labelled in figure 2 In Region I near the injection point,φ(r/R(t)) is shown to decrease slightly with increasing r/R. We conjecture that the suspension is initially uniform upon exiting the injection hole (near r = 0), with the volume fraction of φ 0 = 0.25. Then, as the suspension flows radially out inside the channel, suspended particles gradually migrate across streamlines and focus near the centreline (i.e. shear-induced migration). This gradual transition from uniform to non-uniform particle concentration in the z-direction causes the average particle velocity to increase, which leads to a drop inφ from φ 0 . Hence, Region I is characterised by the transient migration of particles across the streamlines owing to the gradient in shear rate.
In Region II, the suspension flow has reached a quasi-fully developed limit, so thatφ is relatively constant, independent of r/R. We presently refer to this value ofφ in Region II asφ up , whereφ up (φ 0 ) < φ 0 . Then, in Region III,φ starts to increase with r/R, as particles start accumulating near the fluid-fluid interface. As evident in figure 2(a), the maximum value reached inφ consistently increases with time. Hence, the close-up plot of Region III below figure 2(b) reveals thatφ does not collapse into a single curve, when r is normalised by R. However, whileφ may not be self-similar over the entire domain of r/R, the relative sizes of Regions II and III and the general trend inφ strongly suggest that particles accumulate in a self-similar way. In addition, following this peak,φ steeply decreases towards the outer edge of the interface. While this could partly be attributed to the curved shape of the meniscus which affects the measurements, the decrease inφ covers a distance that is consistently larger than the size of the meniscus, h, and again follows a self-similar trend. Hence, the gradual increase and a subsequent decrease inφ are important features of our experimental results, which will be explored in our model.
Despite our characterisation of each region, the boundary between Regions II and III is difficult to determine based on the experimental results in figure 2(b), asφ increases gradually when approaching the interface. Hence, instead ofφ, we consider the extracted values of the radial particle flux, f , as a function of r. Figure 3(a) shows the flux profile for φ 0 = 0.25 over r at t = 6 s. Consistent with the plot ofφ in figure 2(b), f also maintains a constant value away from both the inlet and the interface. However, f increases sharply at a certain radial position, which we currently define as R in (t). Then, R in (t) must correspond to the boundary between Regions II and III, as it clearly separates the quasi-steady region and the particle accumulation region downstream of the interface.
Finally, we characterise the temporal evolution of R in (t) for varying φ 0 at h = 1.4 mm. As plotted in figure 3(b), R in (t) increases linearly with t 1/2 for all values of φ 0 considered. From the volume conservation of the mixture (i.e. πR(t) 2 h = Qt), we have R(t) ∝ (Q/h) 1/2 t 1/2 , where Q and h remain constant in the data plotted. Hence, the linear scaling of R in (t) and R(t) with t 1/2 confirm that the size of the particle accumulation zone must also grow proportionately with R(t), which is consistent with self-similarity inφ.

Theory
3.1. Theoretical formulation: suspension balance model In order to rationalise the self-similar behaviour of theφ profile, we hereby treat the particle-oil mixture as a continuum and develop a thin-film suspension model of our system. We consider the injection of a suspension inside a Hele-Shaw cell from the centre O and define a z-r cylindrical coordinate system as illustrated in the schematics of figure 2(c). The present model is based on the suspension balance model by Nott & Brady (1994), which includes a two-phase formulation expressed for the bulk suspension and particulate phase.
The mass and momentum conservation equations for the mixture are given as where u is the velocity of the suspension and T is the total stress tensor. Similarly, the governing equations for the particle phase correspond to ∂φ ∂t where the superscript 'p' refers to the particle phase; hence, u p and T p represent the particle velocity and stress tensor, respectively. The inter-phase drag force depends on the relative velocity between the particle and the bulk suspension: where the hindrance function corresponds to H(φ) = (1 − φ) 5 (Richardson & Zaki 1954); μ l is the viscosity of the suspending liquid. By combining with (3.1a), we can rewrite ( 3.4) The bulk suspension and particle phase are coupled through the constitutive equations for stress tensors where I is the identity matrix, and E is the bulk suspension rate of strain tensor. We take the following form of the normal stress of the particulate phase in each direction as , ( 3.6) where λ 2 and λ 3 are scalar parameters first defined by Bird, Armstrong & Hassager (1987) independent of φ, which can be determined from rheological measurements of the suspension (Morris & Brady 1998;Zarraga, Hill & Leighton 2000;Boyer, Guazzelli & Pouliquen 2011). Note thatγ = (2E : E) 1/2 is the strain rate; μ s (φ) and μ n (φ) are the effective shear and normal viscosities of the bulk suspension, respectively. Amongst different empirical models for the effective viscosities (Krieger 1972;Morris & Boulay 1999;Zarraga et al. 2000), we presently employ those developed by Morris & Boulay (1999), written as where φ m = 0.6 is the maximum packing fraction of particles. Our experiments are performed by injecting the mixture at a constant volumetric flow rate Q. Hence, the incompressibility of the mixture requires the following constraint for local volume conservation: where the subscript 'r' corresponds to the radial component of the velocity field and h is the channel gap thickness. In addition, the condition of no particle flux at the walls corresponds to where the subscript 'z' is the vertical component.
We hereby introduce dimensionless variables (denoted with an asterisk) to simplify the governing equations: where R 0 is the characteristic radial distance and ≡ h/R 0 , while the characteristic radial velocity is defined as U 0 ≡ Q/(πR 0 h). In particular, since R(t) = (Q/(πh)) 1/2 t 1/2 , the time-dependent location of the interface can be expressed as R * = t * 1/2 based on the dimensionless variables defined.

Model assumptions and multiple-time-scale expansion
Our model rests on a few assumptions based on the experimental set-up and our observations. First, we reasonably assume = h/R 0 1, as R 0 ∼ O(10 −1 ) m and h ∼ O(10 −3 ) m based on our experiments. Next, we consider two important time scales that characterise the current system: τ z , the time scale of particle migration across the thin gap and τ r , the time scale of particle advection in the r-direction. While by examining (3.4) and (3.3). Alternatively, τ z = h 2 /D, where D is the shear-induced diffusivity that scales as d 2γ (Leshansky, Morris & Brady 2008;Griffiths & Stone 2012). Then, the ratio of the time scales corresponds to 1. Based on the characteristic experimental parameters, the typical value of χ ranges from 0.6 to 1. Nonetheless, we consider the limit where χ is a small parameter, such that χ 1, for the following reasons. First, assuming χ 1 allows us to simplify our present analysis and to systematically add the effects of Taylor dispersion into our transport equation, in the manner of Ramachandran (2013). Second, the disparity of the time scales is evident in the experimental plot ofφ (see figure 2b), in which the size of Region I -the transient region near the injection centre -is consistently small compared to the overall radial distance. This implies that the time it takes for particles to migrate across the thin gap must be small compared to that of particle transport in the radial direction. Hence, it is reasonable to assume τ z τ r , or χ 1, as a starting point in the present analysis. Finally, we neglect the transient region near the injection centre (Region I) and only focus on Regions II and III in our model. Given the assumptions above, we first employ the lubrication approximations and reduce the governing equations based on 1. Second, we follow the multiple-time-scale expansion by Ramachandran (2013) and break up the time variable into the short (i.e. t * 1 = t * ) and long (i.e. t * 2 = χt * ) time scales, such that Then, the dependent variables can be expanded in the order of χ. We denote the leading order term with the superscript '(0)' and the term linear in χ with '(1)' (e.g. φ = φ (0) + ( χ )φ (1) + O(( χ ) 2 )). Note that all dependent and independent variables henceforth are presented in their dimensionless form unless otherwise noted; the asterisk ( * ) is subsequently dropped for brevity. Specifically, we model our current system in two steps. First, we compute the particle concentration profiles (e.g. φ (0) , φ (1) ) and radial velocity profiles (e.g. u (0) r , u (1) r ) in the order by order fashion. This yields the expression for the particle flux with respect to the local depth-averaged particle concentration, up to O( χ ). Second, we utilise the expression for the particle flux and derive an evolution equation forφ based on the particle volume conservation, which is solved numerically with appropriate initial and boundary conditions. Overall, our model procedure closely follows the derivation of Ramachandran (2013), which we frequently refer to for more details.

Leading order solutions
Combining (3.1b) and (3.5a), the momentum conservation equations of the bulk suspension can be obtained to the leading order in and in χ as Note that we have assumed the flow to be axisymmetric; hence, the equation in the θ-direction is neglected. Also, at the leading order, the particle transport equation (3.4), combined with (3.2b), (3.3), (3.5b) and (3.6), yields while the no-flux boundary condition (3.9) reduces to ∂ ∂z λ 2 μ n (φ (0) ) ∂u (0) r ∂z z=±h/2 = 0. (3.14) By integrating (3.13) subject to (3.14), we obtain Similarly, based on (3.12b) and the condition ∂ z u r (z = 0) = 0, we can integrate (3.12a) with respect to z to obtain Next, by combining (3.16) and (3.15), we can obtain an implicit expression for φ (0) (z) as where g(φ (0) ) ≡ μ n (φ (0) )/μ s (φ (0) ) and C 0 is a constant independent of z but may have a dependence on r. Then, we integrate (3.16) again and apply no-slip boundary conditions (i.e. u (0) r (z = ±1/2) = 0) to obtain an expression for the mixture velocity, can be determined from the volume conservation of the suspension as 1 = 2r (3.20) Finally, we compute the resultant local particle concentration profile φ (0) (z) from (3.17), by treating the depth-averaged particle concentration,φ = 1/2 −1/2 φ (0) dz, as a known parameter. By systematically varying the value ofφ, we obtain a complete set of solutions for u (0) r (z) and φ (0) (z), examples of which are plotted in figure 4. The plot of φ (0) (z) in figure 4(a) shows that there is a higher concentration at the centreline, z = 0, owing to the shear-induced migration of particles. Then, the corresponding velocity profiles, u (0) r , are normalised by the local mean velocity,ū = 1/2 −1/2 u (0) r dz, and plotted for varyingφ. As shown in figure 4(b), the velocity profile becomes more blunt at z = 0 for increasinḡ φ. This is consistent with the physical picture of particle aggregation near the centreline. With u (0) r (z) and φ (0) (z) thus known, we subsequently compute the dimensionless particle flux, f (0) , as a function ofφ: (3.21)

O( χ ) solutions Within the lubrication framework, we consider the governing equations at O( χ ) to derive u
(1) r and φ (1) . At O( χ ), the momentum conservation of the mixture leads to while (3.4) reduces to . (3.23) We then integrate (3.23) from z = −1/2 to z = 1/2. With no particle flux condition at the walls (i.e. (3.9) at O( χ )), the right-hand side of (3.23) vanishes upon integration, such that On the left-hand side, we apply the chain rule and integration by parts to the second and third terms, respectively, so that Based on the impermeability boundary condition and continuity, we derive the following particle transport equation at O( χ ): By combining (3.23) with (3.26) and employing chain rules (e.g. ∂ r = ∂ rφ ∂φ), we re-write the left-hand side of (3.23) as − ∂φ (0) ∂φ r dz based on continuity. After some algebraic manipulation, (3.23) becomes Finally, integrating (3.28) and applying appropriate boundary conditions leads to where S 7 is a function of z for givenφ. In addition, combining (3.30) with the momentum equations (3.22a) and (3.22b) yields the expression for the perturbation velocity, u (1) r : (3.31) The detailed derivations of φ (1) and u (1) r , as well as the expressions of S 7 and S 11 , are included in the appendix.
Notably, φ (1) and u (1) r depend on the gradient ofφ in the radial direction and onφ, as distinct from the leading order terms. To probe the physical meaning of φ (1) and u (1) r more closely, we plot the functions, S 7 and S 11 , in figure 5. The plot of S 7 in figure 5(a) shows that φ (1) is positive near the walls (z = 1/2) and negative near the centreline before sharply vanishing at z = 0; these effects are increasingly less pronounced at higherφ. This suggests that whenφ increases with r (i.e. ∂ rφ > 0), the effects of shear-induced migration may become mitigated, as the particle concentration is reduced near the centreline and increases near the walls. Similarly, as shown in the plot of S 11 in figure 5 r is positive near z = 0 and becomes negative for increasing z before reaching zero at the walls. While this trend is consistent for allφ, it is reduced at largerφ. This implies that for ∂ rφ > 0, the velocity profile becomes less blunt near the centreline, where the corresponding particle concentration decreases.
This dependence of φ (1) and u (1) r on ∂ rφ will directly impact the dynamics of particle accumulation through the corresponding perturbation in the particle flux, or f (1) . With φ (0) , u (0) r , φ (1) and u (1) r known, we compute f (1) as As shown in figure 6(b), S 12 (φ) is negative and its magnitude generally decreases with φ, while f (0) monotonically increases withφ. Hence, the particle flux must increase more slowly owing to this perturbation in the flux, as particles accumulate on the interface, or as ∂ rφ > 0. We conjecture that this mitigation in the particle flux may lead to the gradual accumulation of particles, which will be fully examined in the next section.

Particle transport equation and self-similarity
In order to computeφ, we must consider the depth-averaged particle transport equation that combines terms up to O(( χ ) 2 ). Hence, we start by integrating (3.4) over the thin gap at O(( χ ) 2 ). By utilising 1/2 −1/2 φ (1) dz = 0 and taking the same steps as O( χ ) (e.g. (3.25a) and (3.25b)), we obtain Finally, we combine (3.26) and (4.1) to derive the following depth-averaged particle transport equation: To find the similarity solution to (4.2), we seek a solution with the following structure: where α and ψ can be found. Substituting it into (4.2), we have (4.4) Note from the plot of f (0) in figure 6(a), we may reasonably approximate f (0) as a linear function inφ, so that f (0) (φ) = C 1 for a constant C 1 . If we assume S 12 (φ) = C 2φ γ for some constant C 2 and γ , and choose α such that αγ + 1 2 = 0, then (4.4) reduces to which provides a rule to determine ψ. As a result,φ = t 1/(2γ ) ψ(r/ √ t) is a self-similar solution to (4.2). In figure 6(b), we fit S 12 (φ) with −0.0022 ×φ −1.399 and show that the latter approximates the former with a slight deviation. This deviation explains why our solution is not exactly self-similar. However, as shown in figure 6(b), the deviation is so small that particles are expected to accumulate in an approximately self-similar manner. Hence, this approximate self-similar nature of (4.1) is consistent with the experimental measurements ofφ in figure 2(b).

Simulation results: comparison with experiments
We solve (4.2) numerically using the upwind finite difference scheme, by applying the following initial and boundary conditions: The initial condition (4.6a) indicates that there are no particles inside the domain at t = 0, while (4.6c) imposes zero particle concentration downstream of the immiscible interface, R(t). The latter simultaneously satisfies the condition of zero particle flux for r > R(t). In addition, we imposeφ at the inlet to beφ up , the constant depth-averaged particle concentration in Region II. In the experiments, the depth-averaged concentration is φ 0 at the inlet and gradually transitions toφ up downstream (Region I). As the current model considers only Regions II and III, it is reasonable to useφ(r = 0) =φ up as our boundary condition, such that the flux at the inlet is equal to the upstream particle flux Hence, we first calculateφ up by revisiting the solution of φ (0) (z) in § 3.3. Specifically, we compute φ (0) (z) from (3.17), subject to the particle volume conservation in the fully developed limit (Region II): φ 0 = 2r Note that the solutions in this regime reduce to the steady-state results of a one-dimensional channel flows, as ∂ rφ = 0 ensures that φ (1) = 0. Finally, we obtainφ up by integrating the resultant φ (0) (z) across the thin gap. Figure 7 shows the plot of φ 0 −φ up versus φ 0 from theory (solid line), together with the experimental measurements (symbols). The results show a qualitative match between theory and experiments, although the theoretical prediction of φ 0 −φ up is persistently larger than the experimental measurements. We subsequently utiliseφ up (φ 0 ) from our model as the boundary condition in (4.6b).
We examine some characteristic solutions to (4.2) for varying values of χ. Note that the second term in (4.2) characterises the radial advection of particles, while the last term corresponds to the radial diffusion, or Taylor dispersion, of particles that acts to smooth out the changes inφ in the r-direction. Hence, the value of χ determines the relative importance of Taylor dispersion versus particle advection in the particle transport equation. Figure 8 shows the plot ofφ(r) at t = 14 s for φ 0 = 0.2 with values of χ ranging from 0.4 to 1. As expected, the increase inφ becomes less steep, as the value of χ is systematically increased. Following this increase,φ is shown to rapidly decrease to zero at the interface, largely independent of the value of χ. Next, we compare the numerical solutions to (4.2) with our experiments for φ 0 = 0.2 and 0.25, respectively, which fall in the regime of clear particle accumulation and no viscous fingering. Here we empirically set the values of χ to best fit the experimental results ofφ, as shown in figure 9(a) and (b). The fitted value corresponds to χ = 0.9 at both concentrations. While the actual values of χ used do not strictly meet the assumption of χ 1, we note that the characteristic size of Taylor diffusion in (4.2) is still an order of magnitude smaller than that of the advection term, as |S 12 | ∼ O(0.01) and f (0) ∼ O(0.1) in figure 6. The numerical solutions ofφ are plotted as a function of dimensional r (in solid lines) at six different times, t. The markers indicate the experimental measurements of the depth-averaged particle concentrations.
The numerical solutions agree well with the experimental observation for both concentrations except for the over-prediction ofφ near the interface particularly for φ 0 = 0.25. This over-prediction is expected given that our current model setsφ(r = 0) =φ up and neglects Region I in whichφ gradually decreases from φ 0 at r = 0 toφ up downstream. We have also plotted the numerical solutions ofφ as a function of r/R in Region III in the inset of figure 9(a). It clearly shows that our theoretical results are also self-similar, despite small deviations near the interface, which is consistent with the experimental observations in figure 2(b). In addition, the numerical solutions of the dimensional flux, f , are plotted and compared with the experimental results in figure 9(c,d) for both concentrations at varying t. The theoretical results are in qualitative agreement with the experiments at all times.
Overall, our reduced model has demonstrated the shear-induced migration of the particles far upstream and the presence of the fluid-fluid interface are barebones physical ingredients that can lead to particle accretion on the interface in the current geometry. Then, the gradual accumulation of particles is achieved through the inclusion of the perturbation in the particle flux that depends on both localφ and the gradient ofφ in the radial direction. This perturbation appears in the particle transport equation as the axial diffusion term, which mitigates how steeplyφ increases towards the interface. The resultant particle transport equation yields solutions that are approximately self-similar and in remarkable agreement with the experimental data.

Summary and conclusions
In summary, our work has focused on understanding the self-similar behaviour of particle accumulation observed in suspension flow, through experiments and theoretical modelling.
This particle accumulation at the advancing interface is attributed to the shear-induced migration of particles far up stream, which causes the particles to assume a faster average velocity than the mixture. We first carry out the experiments by injecting the suspension into the Hele-Shaw cell with varying initial particle concentrations, φ 0 , at a constant volume flow rate Q. We extract the depth-averaged particle volume fractionφ and the particle flux f (r) via image processing. The evolution ofφ as a function of r collapses on to a single curve when r is normalised by the instantaneous position of the interface, R(t). We divide the collapsed curves ofφ(r/R(t)) into three distinct regions: Regions I, II and III. In particular, Region II refers to the region where the suspension reaches a quasi-fully developed flow with constantφ =φ up , whileφ increases towards the interface in Region III. We experimentally quantify the boundary between Region II and Region III as R in (t). The linear scaling of both R in (t) and R(t) with respect to t 1/2 also confirms the self-similar growth of particle accumulation.
We develop a thin-film model to rationalise the self-similarity ofφ based on the suspension balance model. By combining the lubrication approximations and multiple-time-scale expansion, we obtain the depth-averaged upstream concentration,φ up , as well as the particle flux (i.e. f (0) and f (1) ) for given localφ(r) and the gradient of φ(r). Then, we derive a depth-averaged particle transport equation by considering the conservation of the particle volume. In particular, the perturbation in the flux appears in the particle transport equation as the axial diffusion term that mitigates how steeplyφ increases towards the interface. The resultant particle transport equation yields solutions that are approximately self-similar, and in remarkable agreement with the experimental data.
The good match between theory and experiments demonstrates that our reduced model successfully captures the self-similar behaviour of particle accumulation on the advancing fluid-fluid interface. However, the model has a number of limitations that need to be addressed. For instance, the validity of the model for χ ∼ O(1) is questionable, despite the good agreement with the experiments. In addition, the present model results suggest that the presence of the fountain flow and the shape of the meniscus play a minimal role in the dynamics of particle accretion, which is not well understood. Finally, resolving the particle dynamics near the fluid-fluid interface is only the first step in predicting the threshold particle concentration necessary to trigger particle-induced viscous fingering. In particular, previous experimental measurements have shown thatφ needs to attain a critical slope with respect to r before miscible fingering is observed (Luo et al. 2018). However, the physical and mathematical basis for this observation is currently not well understood. Hence, predicting the onset of fingering requires building on the current model forφ(r, t) to derive a fully nonlinear two-dimensional model of the suspension flow inside a thin gap.
Since g(φ (0) )z is also independent of z according to (3.17), the second term on the right-hand side of (A6), [−g(φ (0) )(μ s (φ)∂ z u r ) (1) ], must be a function of r and t only, which will be combined with B 0 (r, t) as a single constant B(r, t). Then, (A6) can be finally simplified as The constant B(r, t) can be determined using the constraint Note that S 6 (φ) = 1/2 −1/2 1 S 4 (φ, z) dz.