Drift-diffusive liquid migration in partly saturated sheared granular media

Abstract We study liquid migration in partly saturated shear bands of granular media where liquid is transported away from the shear-band centre. Earlier studies show that the liquid migration can be modelled as a diffusive process with a shear-rate-dependent diffusion coefficient. Here, we apply this model in a two-dimensional Cartesian split-bottom shear cell with one wide, steady shear band. Initially, a high liquid concentration peak develops at the edges of the shear band, which propagates away from the shear band, splitting the shear cell into a liquid-depleted shear-band region and an outer region not yet affected by the liquid migration. Assuming the liquid transport in the vertical direction is negligible, we simplify the liquid migration model to a one-dimensional form evolving over time. By coordinate transformation, an analytically solvable drift-diffusion model is obtained for liquid migration from the simplified model. From here, we obtain analytical solutions for the liquid concentration as a function of space and time. The significance of the mechanisms is studied in terms of the local Péclet number. While drift enhances drying of the shear band and accumulates the liquid in the peak, diffusion shifts the liquid further away from the shear band. To validate the model, we predict numerically the trajectory of the liquid concentration peak from the continuum model and compare with discrete particle method (DPM) simulations. Our continuum model results give a perfect qualitative and an approximate quantitative agreement with the overall results predicted from the DPM model.


Introduction
Liquid saturation is of tremendous importance to the stability of soil structures. Granular materials generally gain strength with increasing liquid content (Herminghaus 2005;Radjai & Richefeu 2009;Roy et al. 2016;Liefferink et al. 2020) until the liquid saturation reaches a small percentage of the available pore volume. A further increase in liquid saturation in porous soil may cause a dramatic decrease in strength leading, e.g. to landslides or soil collapses (Pailha & Pouliquen 2009;Iverson 2012;Strauch & Herminghaus 2012;George & Iverson 2015;Tomac & Gutierrez 2020). Furthermore, liquid transport or migration induced by shear can lead to a local increase in liquid concentration in the soil pores. Thus, shear-driven liquid migration within soil pores plays an important role for the overall soil properties. Liquid migration is also of great interest in a variety of other applications, such as chemical processing (Rushton 1952), pharmaceutical industries (Cullen et al. 2015), powder technology or in wet granulation processes (Kwant, Prins & Van Swaaij 1995;Jarrett et al. 2016). In wet granulation, grains are mixed with liquid and initial surface wetting is carried out by inducing liquid migration by shearing actions of e.g. the blades inside a rotating device. Thus, understanding the liquid transport phenomena in sheared wet granular media is of great importance for the granular community.
Pore liquids reconfigure in different ways depending on the saturation level of the granular materials. In fully saturated granular materials, liquid is 'sucked' into dilating shear bands (Hicher, Wahyudi & Tessier 1994;Tillemans & Herrmann 1995) with increasing porosity. In contrast, liquid transport at low liquid contents is induced by several different processes. Firstly, it is known that, in shear flows, particles undergo a self-diffusive motion and therefore, liquid which is carried by the menisci will diffuse in space (Campbell 1997;Utter & Behringer 2004). This particle diffusivity has been shown to be proportional to the local shear rate in quasi-static dense flows. Secondly, there is a transport of liquid associated with liquid bridge rupture and formation. Reconfiguration of liquid bridges in the shear band, induced by shear (Long et al. 2019), leads to a local liquid bridge redistribution and liquid transport where liquid is driven out of the shear band (Mani et al. 2012). Both self-diffusion of particles and liquid bridge rupture processes are functions of the shear rate. Thirdly, the equilibrium distribution in the bridge plays a fundamental role in liquid transport (Mani et al. 2015), the liquid transport potential between two capillary bridges on the same grain being proportional to the difference in their capillary pressure. However, we neglect the latter mode of liquid transport in our present study. While the liquid redistribution phenomenon is limited to small shear scale, i.e. happens at the beginning of shearing (Roy, Luding & Weinhart 2018), overall liquid transport is rather a slow process driven by the local shear rate and is the subject matter of this paper.
The focus of our discussion here is to understand the mechanisms of liquid transport in partly saturated granular media at the continuum scale. Liquid migration or transport from the shear band has been understood as a shear-rate-driven diffusion phenomenon with the diffusivity coefficient proportional to the shear rate (Mani et al. 2012;Mani, Kadau & Herrmann 2013). The liquid concentration profile shows remarkable features, particularly in a split-bottom shear cell, where the spatial shear-rate profile is an error function (Ries, Wolf & Unger 2007;Dijksman & van Hecke 2010;Henann & Kamrin 2013) and its width increases with the height in the system (Ries et al. 2007). More precisely, in this set-up, liquid migrates from the zone of high shear rate to the relatively slowly sheared or non-sheared zones. While the shear band gets depleted, a liquid concentration peak is initially observed at the edges of the shear band where the second gradient of the shear rate is largest (Mani et al. 2012). However, whether this liquid concentration peak is stationary over time or behaves differently is still an open question. We address this in the paper by investigating the dynamics of the liquid concentration peak trajectory, using a continuum model for liquid transport. Additionally, we use discrete particle method (DPM) simulations (Athani & Rognon 2019;Vescovi, Berzi & di Prisco 2019;Kuhn & Daouadij 2019;Xiong et al. 2019) in the open source code MercuryDPM (Thornton et al. 2013a,b;Weinhart et al. 2016;Weinhart 2017;Weinhart et al. 2020) to obtain the shear-band width (or velocity profile) imposed in the continuum model, and to compare and validate the liquid transport model. This paper is arranged as follows: the system set-up is explained in § 2 and the continuum model, non-dimensionalisation of the length and time scales and the different features of liquid migration are explained in § 3. The numerical scheme for solving the continuum model is explained in § 3.2. A modification of the liquid migration model to a simplified form is explained in § 5 and the results obtained from this model are compared with that of the full continuum model in § 5.1. We show a suitable transformation of the simplified diffusive liquid migration equation to a drift-diffusion form and the approaches for analytical solutions in § 6.1. In § 6.2, we discuss about the significance of the drift and diffusion processes for liquid migration. We validate the continuum model by comparing the results with the DPM model in § 7. Finally, we draw our conclusions in § 8.

System set-up
We consider a common experimental device to study shear bands, the split-bottom shear cell, which consists of two straight 'L'-shapes sliding past each other, as shown in figure 1. We use Cartesian coordinates where the x-direction is perpendicular and the y-direction parallel to the slit, and the z-direction is perpendicular to the bottom plates (Depken, van Saarloos & van Hecke 2006;Depken et al. 2007;Ries et al. 2007). The left and right 'L'-shapes move along the y direction in opposite directions with speeds −V/2 and V/2, respectively. They are separated by a slit that passes through the origin O. The gravitational acceleration g acts in the negative z-direction. The particle bed consists of particles of uniform diameter d p . The width of the shear cell and the height of the particle bed are denoted L and H, respectively. The interstitial space between particles is filled with liquid with an initial homogeneous liquid concentration Q | t=0 = Q 0 . In steady state, the flow is uniform in the y-direction and a shear band propagates from the split position O upwards. We follow the observations of Singh et al. (2014) and assume that the shear band has a Gaussian velocity profile of width with W top the shear-band width at the surface of the flow, and the exponent 0 < α < 1. We further assume that the shear in the z-direction can be neglected and thus the magnitude of the local shear rate,γ , is approximated by the gradient of the y-velocity of the particle phase,γ (2.2) Since the system is symmetric around the y-axis, we study only the right half of the shear cell.

Continuum model
The nature of the transport equation governing liquid migration is given as (Mani et al. 2012 where Q is the liquid concentration, or volume fraction of liquid, expressed in dimensionless form as the volume of liquid per unit volume,Q is the rate of change of liquid concentration Q andγ is the shear rate given by (2.2). The description of the liquid migration model originally comes from Mani et al. (2012), where the diffusion mechanism is explained in terms of a theoretical model. The main mode of liquid transport in this model happens via rupture of individual capillary bridges. The bridge rupture rate is proportional to the shear rateγ and the number of contacts N. According to Mani et al. (2012), C liq is proportional to the number of contacts, and thus weakly depends on the pressure, or depth z and is also proportional to a geometrical proportionality factor which measures the average volume of liquid leaving the control volume after each rupture event. However, for simplicity, we assume here that the prefactor C liq (which is not the diffusion coefficient) is constant. Thus, the physical significance of C liq is based on the geometric configuration as well as on the packing fraction of the granular materials.
3.1. Non-dimensionalisation It is convenient to redefine the governing equation (3.1) in terms of dimensionless length and time scales. Therefore, we scale the spatial x-and z-coordinates by the particle diameter d p , We further scale the time t by the shear rate at the initial liquid concentration peak location near the free surface,γ c s (evaluation shown in appendix A), All the variables henceforth are non-dimensionalised and we omit the superscript * subsequently.
3.2. Numerical scheme Numerical methods suitable for the solution of the fluid transport equations are a matter of extensive research in computational fluid dynamics. Their application is predominant in various research fields such as geophysical fluid dynamics (Durran & Mobbs 2001; Baumgarten & Kamrin 2019), hydrological processes (Igboekwe & Achi 2011), reactor flow (Hastaoglu & Abba 1996) etc. In the Eulerian solution of equations, difficulties arise because of the dual advective-diffusive nature of the transport equation. When the transport is advection or drift dominated, the equation behaves as a first-order hyperbolic equation, but when the transport is diffusion dominated, the equation behaves as a second-order parabolic equation. To accurately model the drift-diffusion transport, the numerical scheme must be able to handle the mixed parabolic-hyperbolic character of the systems. Eulerian models that have grids fixed in space have a number of difficulties when transport is drift dominated. These include numerical diffusion, oscillations, instabilities and peak clipping because of the numerical representation of advection terms in the transport equation. To avoid instabilities, we choose the finite volume method (FVM) with semi-implicit time stepping (Patankar 1980) to solve the liquid transport equations in this paper. The details of the numerical methods and discretisation are described in appendix B.
The grid sizes are chosen as 400 in the x-direction and 100 in the z-direction, which is equivalent to a grid spacing of dx = dz ≈ 0.08. The solutions are checked for different grid sizes (the finest resolution tested is dx = dz ≈ 0.008) and the trend is maintained both qualitatively and quantitatively. The time step is chosen as dt = 10 −4 . These values meet the necessary Courant-Friedrichs-Lewy (CFL) condition for the stability of the solutions with CFL = 0.59 < CFL max . The details of the calculation of the CFL number is elaborated in appendix C. It is important to emphasise that this is not a sufficient condition for stability and other stability conditions are generally more restrictive than the CFL condition.

Characteristics of liquid migration
Shearing an unsaturated granular system causes a re-distribution of the interstitial liquid.
While an initial transient behaviour shows random local re-distribution of the liquid, a larger shear leads to transport of liquid from the shear zone (Roy et al. 2018). The resultant liquid concentration profile in the shear cell geometry is worth describing in detail.

Liquid concentration profile
Initially, the liquid concentration Q | t=0 = Q 0 = 6.9 × 10 −3 is homogeneous; thus, the time derivative of the liquid concentration Q is proportional to the second spatial derivative ofγ i.e. ∇ 2γ . Thus, the liquid concentration initially decreases in the shear band, where the second derivative ofγ is negative, and a liquid concentration peak is formed at the edge of the shear band, where the second derivative ofγ is largest.

Liquid concentration peak
The location of the liquid concentration peak could be defined as the location of the maximum of the liquid concentration profile. However, this has the following disadvantages: as the grid resolution is finite, the liquid concentration peak can only be located with a limited accuracy, resulting in an undesirable stepwise definition. To avoid these effects, the location of the liquid concentration peak in the right half of the shear cell is defined as the centroid of the liquid concentration profile above the initial value Q 0 . Thus the definition of x c is given as where L = 60 is the width of the shear cell in non-dimensional form. The location of the liquid concentration peak is well approximated by this definition for the continuum model. However, if Q 0 is more scattered, as in the case of data from DPM simulation, we definẽ Q = max(Q − (1 + )Q 0 , 0), where 1 is the standard deviation of Q 0 in the data.

Accumulated liquid concentration
The concentration of liquid accumulated in the edge of the shear band is given as the integral of the liquid concentration profile lying above the value Q 0 as follows: Ideally, there is no loss of liquid from the system (e.g. due to vaporisation). The conservation of liquid volume requires that the volume of liquid accumulated in the edge of the shear band equals the volume of liquid drained from the centre of the shear band. Figure 2 shows a typical liquid concentration profile Q as a function of x at a fixed height, z = 3.6 (W = 3) for the full liquid migration model given by solving (3.4) after t = 3.6. We distinguish the two main features of the liquid concentration profile, namely the peak concentration location x c and the accumulated liquid concentration φ. Further, the dynamic characteristics of these features are the subject of our discussion as and when we simplify the governing equation for liquid migration in the following sections § § 5, 6 and 7.

Simplified model neglecting vertical diffusion
In order to do a detailed theoretical analysis of the mechanisms of the liquid migration process and for simplicity, we reduce (3.4) to a simplified form, neglecting the diffusion in the z-direction. Thus, (3.4) is simplified as The original equation of liquid migration is represented by (3.4) and its simplified form is given by (5.1) . We denote these equations as the full model and the simplified model, respectively. Both equations are solved using the finite volume scheme described in § 3.2. In the following subsection, we make a comparison of the results obtained from the full model with those of the simplified model. Figure 3(a,b) shows contour plots of liquid concentration as a function of space x and z after t = 14.4 solved for the full model and the simplified model, respectively. As observed from the figures, both the models have a minimum liquid concentration at the shear-band centre, i.e. at x = 0, corresponding to the region with dark blue colour. A high liquid concentration is developed at the edges of the shear band for both the models, corresponding to the narrow region with dark red colour. The region close to the boundary, represented by the cyan colour, is not yet affected by the liquid migration and the liquid concentration is unchanged. A height-wise gradient of the liquid concentration is observed inside the peak location for the full model, represented in figure 3(a). This is due to the liquid diffusion in the z-direction. Unlike the full model, the simplified model shows a rather uniform liquid concentration inside the peak location, represented in figure 3(b). Starting to shear from a uniform concentration of liquid Q 0 , the initial location of the liquid concentration peak after a single time step t = dt, is given by x 0 c . This location is obtained analytically for the simplified model from (5.1) as x 0 c = √ 1.5W where ∂ 2γ /∂x 2 is maximum. Note that, for the full model, x 0 c is at the location where ∇ 2γ is maximum. The red solid line in figure 3(b) shows the locus of x 0 c at different heights for the simplified model. The liquid concentration propagates away from the shear band with time and the red dashed lines in figure 3(a,b) represent the location of the peak after t = 14.4 for the two models. Note that this is an intermediate time chosen to show the liquid concentration profile when the initial liquid redistribution phase has ended and the liquid migration phenomenon has started.

Comparison of the full and simplified models
Next, we do a quantitative analysis of the liquid concentration peak x c and the accumulated liquid concentration φ as a function of time at different heights picked from figure 3(a,b). The results of x c and φ at z = 3.6 (W = 3) (blue lines) and z = 5.4 (W = 3.5) (red lines) are shown in figure 4(a,b), respectively. The key parameters x c and φ both increase with time, indicating that the process has not reached a steady state. The solutions of the full and the simplified models are represented by the solid and the dash-dotted lines, respectively. While there is only a difference of less than 5 % in x c between the full and the simplified models, φ is significantly affected by the diffusion in the z-direction. The accumulated liquid concentration increases by 33 % more for the simplified model as compared to the full model after t = 72, closer to the base of the shear cell (blue lines), where the vertical shear gradient is stronger. The location of the liquid peak position x c is insignificantly affected by diffusion in the z-direction as only the horizontal diffusion shifts the liquid peak away from the shear band. Thus, x c is captured with up to more than 95 % accuracy through the simplified model and is the key parameter that we want to focus on by further analysis in this paper. Thus, an approach towards a simplified model solution is likely to deviate quantitatively from the full model in the first place, especially close to the bottom of the shear cell. However, the location of liquid concentration peak x c can be readily captured, which in itself is an important feature of the liquid concentration profile. The location of the liquid concentration peak deviates for the simplified model as compared with the full model by less than 5 %. Although the horizontal x-diffusion is primarily shifting the liquid concentration peak, the component of the z-diffusion that is perpendicular to the liquid concentration profile is also contributing to the process. Hence, a difference is observed between the location of the liquid concentration peak between the full and the simplified models.
Next, we investigate the velocity of propagation of the liquid concentration peak v c as a function of the peak location x c . Figure 5 shows the dependence of v c on x c at three different heights for the full model (solid lines) and simplified model (dashed lines). The propagation velocity of the peak position decreases with increasing x c as the peak moves away from the shear band. Note that the initial peak locations x 0 c are different for the full and the simplified models and hence the x c ranges are also different for the two models. Initially, the liquid peak is not well developed (for small x c ), resulting in a subtle peak whose location is difficult to analyse. These data for the initial time steps are eliminated from our analysis. The velocity of the simplified model (dashed lines) at a lower height, close to the split position, deviates from the behaviour of the full model (solid lines) by up to 20 % at z = 2.7. This is due to the absence of diffusion in the z-direction, which is more prominent at a lower height, close to the split position. The effect of diffusion in the z-direction is weak near the surface and thus the velocity profiles for the full and the simplified models almost collapse close to the free surface at z = 5.4.
It is observed that the trajectory of the location of the liquid concentration peak x c can be predicted from the simplified model with an accuracy of 95 % as compared with that of the full model. The change of velocity of propagation of the liquid concentration peak location v c as a function of the local shear rate is closely predicted by both the full and the simplified models near the free surface, but deviates significantly near the bottom of the shear cell. The variation of the accumulated liquid concentration for the full and the simplified models as a function of time is closer near the free surface, but deviates by approximately 20 % near the bottom where the shear rate is higher. To summarise, the deviation of predictions of accumulated liquid concentration given by the simplified model from the full model is higher where the local shear rate is higher. Also, the liquid peak propagation velocity is proportional to the local shear rate, with a zero velocity corresponding to zero shear rate, indicating that the liquid migration is a dynamic process which is solely shear driven.
The simplification of the model allows further analysis and the development of analytical solutions. We propose to show analytical solutions for the simplified model with suitable transformations of the equation in § 6.1. We choose an intermediate height z = 3.6, where the effect of the local shear rate is moderate and thus the results of analytical predictions are closer to the results of the full continuum model.

Transformation of equation
The fundamental challenge is to understand and predict the transport of interstitial liquid in sheared, partly saturated granular materials. While the simplest picture of diffusive transport with a constant diffusivity cannot explain the dynamics of liquid transport, a model with a variable, shear-rate-dependent coefficient of diffusion can. However, multiple effects happening at the same time, some of which lead to drift-like rather than diffusive transport features (rapid build up and narrowing of the liquid front), make the basic understanding difficult. Therefore, by transforming the variables one can enforce a diffusion term with constant diffusivity D c , which yields a drift term with a variable drift coefficient. This split allows us to study the two terms separately and (with some further simplifications) solve them analytically. Furthermore, the mechanisms of liquid transport can then be separated and understood one by one, where diffusion is a randomising driving term, but the drift, i.e. drift-like transport, is generated by the shear rate due to the shear banding and flow profile.
The details of the transformation of one-dimensional (5.1) and the resulting analytical solutions are explained in this section. We choose an intermediate height z = 3.6 for transforming the simplified model (5.1) and at this height the width of the shear band is W ≈ 3.

Transformation of the equation: drift and diffusion
Next, we aim to transform (5.1) into a form that is analytically more treatable. We follow the approach of Risken (1989) and apply a coordinate transformation, The liquid distribution in the transformed coordinate, Q (t, ξ), is thus given by Applying this change of variables to (5.1) yields a diffusion and a drift term, This equation has a constant diffusion coefficient (equal to 1) and a variable drift coefficient, Applying the coordinate transformation used in the paper has the advantage of a constant diffusion coefficient, which in our opinion results in a split that is better to analyse and understand. Moreover, this form of decomposition adopted by us to analyse the liquid migration phenomenon is rather a novel approach compared with the conventional chain rule decomposition. In order to see the individual contributions of the drift and diffusion terms on the overall liquid transfer, we separate the drift and diffusion processes and show analytical solutions for each in § § 6.1.1 and 6.1.2, respectively.

Drift
In order to measure the contribution of the drift term to the liquid transport, we neglect the diffusion term in (6.3), obtaining the simpler equation The general solution of (6.7) is given by where A is an anti-derivative, and R 0 is defined such the initial condition, Transforming back to the original variable Q drift , we obtain the analytic solution Note, this analytic solution is valid for any shear-rate profileγ (x). Substitutingγ from (2.2) into the definitions of A and E, (6.9) and (6.12a,b), we obtain The plots of Q drift (t, ξ) in the transformed coordinate and Q drift (t, x) in the original coordinate are shown in figures 6(a) and 6(b), respectively. The initial condition, Q drift = Q 0 C liqγ (x) is simply a Gaussian function of x. The initial peak location of Q drift at ξ = 0 is as shown in figure 6(a). One can see from the inset of figure 6(b) that, initially, For large values of t, x 0 ≈ 0 for small x (hence Q ≈ 0) and x 0 ≈ x for large x (hence Q ≈ Q 0 ); in between, the facts that x 0 < x and E(x) has a maximum ensures there is a peak value. As observed from figure 6(b), drift induces the liquid concentration peak x c to move further away from the shear-band centre, leading to complete rapid drying of given by (6.13).
the shear band and surroundings by pushing the liquid to the peak region. As a result, we observe the liquid concentration profile forming a dry shear band that is surrounded by a wet region, with a sharp liquid concentration peak between the dry and wet regions, as shown in figure 6(b). We have verified that the analytical solution Q drift (t, x) given by (6.13) agrees with the numerical solutions of (6.7) (results not shown here). Note that (6.1)-(6.5) are valid for any shear-rate profileγ (x). We only apply the specific value oḟ γ (x) in (6.6)-(6.12a,b) to obtain an analytical solutions. Based on this analysis, we can derive an estimate, x e (t), of the peak position. First, we define we define x 0 e as the peak location of E(x), thus E (x 0 e ) = 0, with denoting the derivative with respect to the variable x. Next, we define x e (t) for t > 0 such that x 0 (t, x e ) = x 0 e . Note that x e > x 0 e , since x 0 is monotonically increasing in x, see the inset of figure 6(b). Further, note that E (x 0 e ) < 0, since E is monotonically decreasing for x > x 0 e , see the inset of figure 6(a). We will now show that x 0 c < x e < x c : substituting x = x e into (6.13) we get (6.14) Since E(x 0 e ) is the maximum value of E, we get Q drift (t, x e ) > Q 0 , and thus x e > x c 0 . Next, we show x e < x c : Differentiating (6.13) and substituting x = x e , we get The first term in the numerator is zero, Thus, x e yields a (lower) estimate of the peak location, in particular for large times, as we observe that the peak narrows over time. Note Thus, x e (t) has the shape of the inverted exponential integral, shifted to the right by a constant. A plot of x e is given in figure 7. As observed in the figure, the peak x e moves away from the shear band with increasing time, leading to drying of the shear band. Thus, the final behaviour of the solution for the drift operator is to push the dry front to infinity, resulting in a completely dry domain everywhere, however, this is a very slow process, as apparent from the figure.

Diffusion
In order to measure the contribution of the diffusion term to the liquid transport, we neglect the drift term in (6.3) An analytical solution of (6.16) is given by the convolution of the initial condition, Q 0 = Q 0 C liqγ , with a kernel function l, where l is defined as To prove that (6.17) satisfies (6.16), it is sufficient to show that (6.20) We solve (6.17) numerically in Matlab to get the final solution. The plots of Q diffusion (t, ξ) in the transformed coordinate and Q diffusion (t, x) in the original coordinate are shown in figures 8(a) and 8(b), respectively. Figure 8(a) shows a typical solution of the diffusion equation in the transformed coordinate, where the peak of the Gaussian solution decreases and broadens with increasing time. The corresponding solution in the original coordinate is shown in figure 8(b). Unlike drift, figure 8(b) shows that diffusion induces a smooth liquid concentration peak x c instead of a sharp interface and slowly leads to the drying of the shear band and its surroundings. As a result, we observe a smooth liquid concentration profile with a nearly dry shear band and its periphery and a subtle liquid concentration peak that moves away from the shear band, shown in figure 8(b). We have verified that  (t, x) in the transformed coordinate system. the analytical solution Q diff (t, x) given by (6.17) agrees with the numerical solutions of (6.16) (results not shown here). So far, we have simplified the full model for the liquid migration process and obtained analytical solutions for the simplified forms. In § 6.2, we use the analytical solutions for the drift and diffusion terms and find their significance, individually, in the overall liquid transfer process.

Significance of drift and diffusion
In this section, we explore the significance of the transformed drift and diffusion processes, respectively, as a part of the overall liquid transport process as given by the one-dimensional simplified model equation. The results of the drift and diffusion contributions are obtained from (6.13) and (6.17), respectively, by analytical solutions which are considered as new mathematical tools of great significance. While both drift and diffusion contributions change over time, they behave qualitatively the same, i.e. having a minimum at the centre x = 0 in the original coordinate system, a propagating liquid concentration peak at the edges of the shear band x c and a constant liquid concentration near the boundary region. Furthermore, the liquid concentration peak location as well as the accumulated liquid concentration changes over time for the solutions of all the aforementioned models. The new mathematical tools are used to investigate further the relative significance of drift and diffusion in the liquid migration process. In the following subsections, we explore the contribution of drift and diffusion processes, individually, to the velocity of propagation of the liquid concentration peak v c and to the accumulated liquid concentration φ at a given height of the shear cell z = 3.6 (W = 3).

Local Péclet number and velocity of propagation of liquid concentration peak
Among the associated dimensionless numbers, we identify the Péclet number for the study of liquid transport in granular media. The Péclet number is a class of dimensionless numbers which measures the rate of advection of a physical quantity by the flow to the rate of diffusion of the same. This dimensionless number is relevant for the study of transport phenomena in fluid flows, signifying the transport by drift and diffusion processes. We define the local Péclet number Pe for liquid transport, obtained from (6.3) in the transformed coordinate, as Pe = D (ξ )ξ , where D (ξ ) is the drift coefficient, the diffusion coefficient is constant (equal to 1) and ξ is the characteristic length in the  .17), respectively, at z ≈ 3.6 (W ≈ 3). The cyan line represents the sum of the drift and the diffusion velocities. The relative significance of the drift and diffusion velocities in (b) is comparable with the corresponding local Péclet number in (a) at a given height z = 3.6.
transformed coordinate system. Physically, this dimensionless number signifies the ratio of mass transfer by drift to that by diffusion processes in the ξ -direction. Note that, although the Péclet number is defined here in the transformed coordinate, we analyse the corresponding results in the original coordinate system. Figure 9(a) shows the variation of the local Péclet number in the domain. The local Pe is independent of time since it is only varying with the shear rateγ and space x, which are independent of time. There is an inner region in the vicinity of the shear band centre where Pe 1, thus diffusion dominates liquid transport and drift is insignificant. Simultaneously, there is an adjoining outer region where the effects of drift and diffusion become comparable (Pe ≈ 1).
In figure 9(b), we compare the liquid concentration peak velocity for the simplified model, drift and diffusion, respectively, at a given height z = 3.6 (W = 3). The red line corresponds to the velocity of the simplified model, v csimplified model . The blue and the green lines correspond to the propagation velocity of liquid concentration peak for the drift and diffusion models, v cdrift and v cdiffusion , respectively. The cyan line represents the sum of the velocities of propagation for the drift and diffusion models, individually, i.e. v cdrift + v cdiffusion . It is observed from figure 9(b) that the red line collapses with the cyan line. This indicates that net contribution to the velocity of propagation of the liquid concentration peak from drift and diffusion is approximately equal to the velocity given by the simplified model throughout the range of x c and therefore v csimplified model = v cdrift + v cdiffusion . (6.21) Thus, we observe that the combination of drift and diffusion processes has an additive effect on the velocity of propagation of the liquid concentration peak for the simplified model. This also shows the validity of the simplified liquid migration model. We also investigate the relative contribution of drift and diffusion to the overall liquid transfer and make a comparison with the local Péclet number. Focusing on the local Péclet number at a given height z = 3.6 in figure 9(a), one observes an inner dark blue region where Pe < 1. Simultaneously, there exists an outer dark red region in the vicinity where Pe > 1. This profile of the Péclet number can be related to the velocity profile described in figure 9(b). The contribution by diffusion is approximately 80 % of the total contribution by drift and diffusion at x c = 4.5. This corresponds to the dark blue region in figure 9(a) where Pe ≈ 0.5. The contribution by diffusion decreases with increasing x c and thus the Péclet number increases. The contributions by drift and diffusion are approximately equal at x c = 6, corresponding to Pe ≈ 1. On further increasing x c , diffusion becomes less dominant and its contribution is only 33 % of the total contribution at x c = 8. This corresponds to the dark red region in figure 9(a) where Pe ≈ 1.3. The mass transfer is dominated by diffusion and is restricted to the inner region where the velocity of liquid concentration peak propagation v c is two orders of magnitude higher than that in the outer region. This is in similar philosophy with the steady-state mass transfer at low Pe (Pe 1) (Jones 1973;Neale & Nader 1974;Srivastava & Srivastava 2006;Bell et al. 2014), a typical example of Brinkman and Darcy flow. Under conditions of low Reynolds number, for transport of liquid past granular materials, in the quasistatic state, the transverse diffusion flow is more important than the drift flow. Thus, we understand the significance of the Péclet number with respect to the propagation velocity of the liquid concentration peak. This shows how the simplified model is used to understand the essence of the liquid migration process.
In this subsection, we investigated the relative significance of drift and diffusion in the velocity of propagation of the liquid concentration peak. Thereby, we validated the simplified model, which is an important mathematical tool for studying the physics of the liquid migration process. The accumulated liquid concentration associated with the increment of the liquid concentration peak is also an important feature of the liquid concentration profile. Thus, in the following § 6.2.2, we investigate the relative contributions of drift and diffusion to the accumulated liquid concentration required for increase of the peak of the liquid concentration.

Accumulated liquid concentrations
The differential equation for the transport of liquid from the shear band is composed of two terms in the transformed coordinate, namely, diffusion and drift. While the drift is the directed flow of liquid by the bulk motion, diffusion leads to random spreading of liquid under the influence of shear from higher to lower shear rate. However, in spite of different significance of drift and diffusion on the process, the accumulated liquid concentrations contributed by the two are additive. It is expected that, if the two accumulated liquid concentrations are separately integrated over time, the resulting accumulated liquid concentrations should be comparable with those of the simplified model, at least for a small incremental time scale, when other effects, such as coupling, are negligible. We verify this from our results.
In this section, we analyse the incremental accumulated liquid concentrations Δφ tot Δt from the two contributions by drift and diffusion separately. The objective is to see if the resultant accumulated liquid concentrations from the two components of the drift and diffusive terms, together, contribute to the same incremental accumulated liquid concentrations as obtained by integrating (5.1). Moreover, we check the validity of this assumption for different initial conditions. We start running the drift and diffusion models from different initial conditions, such that they have an initial accumulated liquid concentration equal to that of the simplified model φ simplified model at time t. The net accumulated liquid concentration φ tot Δt in time Δt contributed by the drift and diffusion processes is given by  where φ drift Δt and φ diffusion Δt are the contributions to the incremental accumulated liquid concentration by drift and diffusion, respectively. All the aforementioned accumulated liquid concentrations are obtained as described in (4.2). Figure 10(a) shows the comparison of the sum of the accumulated liquid concentrations contributed by drift and diffusion, indicated as φ tot Δt , with the accumulated liquid concentration for the simplified model, φ simplified model . We observe that, within a very small time interval Δt = 0.72, φ tot Δt ≈ φ simplified model Δt , signifying that the effects of drift and diffusion are additive within a very short duration of time. However, with further progress in time, the net contribution from the two effects φ tot Δt over-predicts the accumulated liquid concentration φ simplified model Δt . The deviation of the accumulated liquid concentrations decreases at longer time (or shear) when the incremental drift and diffusion fluxes become weaker. This is shown in figure 10(b) where the simulations are run, starting at t = 288 and for Δt = 72. The trends show that the two accumulated liquid concentrations φ tot Δt and φ simplified model Δt coincide, indicating that the incremental accumulated liquid concentrations are equal even after a long time. The contributions by drift and diffusion to the accumulated liquid concentrations are shown by the blue and the green lines respectively. While drift leads to a positive (increasing) contribution to accumulated liquid concentration, diffusion leads to a negative (decreasing) contribution and the net accumulated liquid concentration is constant relative to the simplified model. This shows another method of validating the simplified model and the mathematical tools used for predicting the drift and diffusion behaviour of liquid migration.
In this subsection, we have seen that the liquid volume required for the increment of the liquid concentration peak from any initial condition is contributed from the sum of the liquid volumes from the drift and diffusion processes. This observation is valid for a small time interval and the sum deviates from the liquid volume of the simplified model for a longer time interval. Further, the deviation slows down as the local shear rate is weaker when the peak moves away from the shear-band centre. Depending on the initial condition, the contribution by diffusion to the accumulated liquid concentration might be incremental or decremental. In § 6.2.3, we analyse in detail the sign of the incremental accumulated liquid concentrations for drift and diffusion processes, which determines whether the shear band wets or dries.

A30-17
One of the major results found in this subsection is that after reducing to a quasi-one-dimensional system and transforming variables, one can separately solve the equation without a diffusion term or without a drift term. The two solutions, when added together, agree with solutions of the unseparated equation when proper initial conditions are used. The fundamental approximation of this additive splitting can be shown mathematically for an incremental time step and this is described in appendix D.

Drift vs diffusion: drying and wetting
In this section we analyse the increase in accumulated liquid concentration by drift and diffusion individually, relative to the accumulated liquid concentration at any time t. We start simulations corresponding to drift and diffusion with an initial condition of an accumulated liquid concentration φ simplified model corresponding to time t. The relative increase in accumulated liquid concentration by drift in time Δt with respect to the initial condition is given by and we express the relative increase in accumulated liquid concentration by diffusion in a similar way. Since the total liquid is conserved in the system, an increment in accumulated liquid concentration at the peak location indicates a decrement of equal amount of accumulated liquid concentration in the shear band. Thus, a positive value of η (Δφ > 0) signifies drying of the shear band with respect to the initial condition and a negative value of η (Δφ < 0) signifies wetting of the shear band. Figure 11 shows the relative increment in the accumulated liquid concentration for the drift and diffusion processes as a function of the initial accumulated liquid concentration φ simplified model . Here, η decreases for both drift and diffusion processes with increasing φ simplified model until it reaches a condition such that the change in accumulated liquid concentration ratio η becomes very small. As the liquid concentration peak propagates away from the shear band with increasing time, the local shear rate becomes smaller. With decreasing shear rate, the driving forces for liquid transport by both drift and diffusion mechanisms become slow enough, such that it requires very long time to transport liquid. Diffusion leads to wetting of the shear band under certain initial conditions when liquid is transported back to the shear band such that η becomes negative. However, drift always leads to drying of the shear band only. This is yet another aspect of the physics of liquid migration which we are able to understand by disintegrating the simplified model into drift and diffusion components.
Thus, the mechanism of liquid transport can now be separated and understood one by one, under varied initial conditions. While drift always leads to drying of the shear band (with rapid build up and narrowing of the liquid front), diffusion can sometimes lead to wetting of the shear band, depending on the initial condition. We studied in details the role of drift and diffusion in the liquid migration process in § 6.2. The origin of this theoretical studies begins with the model given by (3.4) which was originally proposed and validated by Mani et al. (2012). We compare this model once again with the DPM simulations as described in § 7.   (Mani et al. 2012). To validate the proposed model given by (3.4), we simulate a simple linear split-bottom shear cell as described in § 2 using DPM. The so-called DPM is used for modelling of particulate materials as an approach towards the macroscopic understanding of microscopic behaviour. Contact models are at the physical basis of DPM simulations. We perform DPM simulations using the open source code MercuryDPM (Thornton et al. 2013a,b;Weinhart et al. 2016;Weinhart 2017;Weinhart et al. 2020). To model the Cartesian split-bottom shear cell (see figure 12), small particles are glued to the sidewalls and bottom to make the surface rough and the shear cell is filled in with particles. A periodic boundary condition is applied in the y-direction. All the parameters for particles and the contact model for the DPM simulations, which are used in their dimensional form, are given in table 1. All the other dimensions of the shear cell geometry and parameters of the particle and contact model are scaled by the particle diameter d p , gravity g and particle mass m p . Note that, although we use d p andγ c s as original scaling parameters in this paper, we scale the input parameters for our DPM simulations by d p , g and m p . Therefore, the mentioned scaling of all the parameters for DPM simulations are shown in table 2. The non-dimensional value of gravity g in terms of the original scaling parameters d p andγ c s is equal to 3.5 × 10 4 (equivalent to g = 10 ms −2 in dimensional form). All the non-dimensional values of other parameters in table 2 in terms of the original scaling parameters d p andγ c s can be obtained by substituting g = 3.5 × 10 4 , d p = 1 and m p = 1. The details of the contact model are given in Roy et al. (2016) and we describe the mechanism of liquid bridge formation and rupture in appendices 1 and 2, respectively.  The left colour bar represents the magnitude of the liquid bridge radius and the right colour bar represents particle velocity v y in the y-direction. Note that the colours and scales of the liquid bridge radius and the particle velocity v y are adjusted to a suitable range for better visualisation. Figure 12 shows a snapshot from the DPM simulation showing the liquid migration from the shear band. The particles are coloured according to decreasing particle velocity v y from red to blue. The liquid bridges are shown in the form of cylinders with their length signifying the radius equivalent to liquid bridge volume at the contact. The liquid bridges are coloured according to decreasing liquid bridge radius from red to blue. Note that the colours and scales are adjusted to a suitable range for better visualisation. It is evident from the figure that the liquid bridge concentration is lowest inside the shear band, highest near the edges of the shear band and has an intermediate concentration near the walls.

Discrete to continuum
We use the coarse-graining post-processing tool MercuryCG (Thornton et al. 2012;Weinhart et al. 2012aWeinhart et al. ,b, 2013Tunuguntla, Thornton & Weinhart 2016;Tunuguntla, Weinhart & Thornton 2017a,b) to translate DPM data from discrete to continuum scale, averaged over time and over the y-direction, at x and z positions on a 400 × 100 grid over the whole shear cell. We use a Gaussian spatial coarse-graining function with a coarse-graining width (standard deviation) of one particle diameter. Since the liquid concentration profile is dynamic in nature, we temporally averaged over a short period of time of Δt = 0.25 (5 snapshots) in order to get the dynamic behaviour of liquid transport. Thus, we obtained average values for the liquid concentration Q present in the liquid bridges and liquid films of the DPM simulation. The details of the coarse-graining formulation are given in appendix E.3. Likewise, we obtained other continuum fields, such as the local pressure, concentration, velocity and the shear rate, in order to draw the shear-band profile. The velocity of particles v y in this geometry is typically fitted by an error function of the spatial length x at a given height (Fenistein, van de Meent & van Hecke 2004;Ries et al. 2007). Figure 13(a) shows the particle velocity v y as a function of x at z = 3.6 and z = 5.4. The red solid and the dashed lines represent the fitting of the velocity profile at the two heights, respectively. We obtain the width of the shear band W by fitting the particle velocity v y as a function of x at different heights z. Thereby, we obtain the shear-band width W as a function of z, as shown in figure 13(b). The red solid line in this figure corresponds to the fitting of the data by (2.1) with the fitting parameters W top ≈ 4 as the width of the shear band near the free surface, H = 8 as the filling height and α = 0.88 as the power obtained by fitting DPM data as detailed by Singh et al. (2014). Substituting (2.1) in (2.2), we obtain the shear-rate profile in the continuum models. Figure 14(a,b) shows a comparison of the results from the continuum full model with those from the coarse-grained DPM simulations. The location of the peak liquid concentration of the continuum model is well aligned with the peak liquid concentration of the DPM model, as observed in figure 14(b). Liquid gets depleted from the region near the shear band and a liquid concentration peak appears in the neighbouring region. The liquid concentrations close to the right boundary remain unchanged for both the DPM and continuum models. The location of the peak liquid density for the continuum model, denoted by the magenta-coloured markers in figure 14(b) is well aligned with the liquid density peak for the DPM model. The peak liquid concentration is decreasing with height for both the continuum and DPM models. However, certain differences still exist in the liquid concentration profiles of the two models: the curvature of the liquid concentration peak locus is nearly linear in the DPM simulations, but convex in the continuum model results, as observed in figure 14(a,b). Also, the width of the liquid concentration peak for the DPM model is slightly larger than that of the continuum model, as observed in figure 14(a,b). Figure 15(a) shows a comparison of the liquid concentration profile as a function of space x at two different heights z = 3.6 (W = 3) and z = 5.4 (W = 3.5) for the DPM (scattered points) and continuum models (solid lines) after time t = 72. The comparison shows that the liquid concentration profiles are in good qualitative agreement for the two models, although the liquid concentration peak is slightly lower for the continuum model at z = 5.4. This is also evident from the contour plots of the two models in figure 14(a,b). The peak location x c is well aligned for the two models at both heights. As observed in figure 15(b), the trajectories of the peak location x c for the two models are in good agreement. There are various possible explanations for some quantitative differences between the two models: firstly, we fitted the continuum model with the DPM data by using a constant value of C liq . However, our soft particles in the DPM simulation are subjected to a confining pressure under gravity. The particles near the base are therefore more compressed than the particles at the free surface. Thus, the number of contacts per particle near the bottom is slightly higher than at the free surface (not shown). Secondly, in the continuum theory, we assume a steady-state shear-rate profile at the beginning of the simulation. However, the DPM simulations show an evolving shear rate inside the shear band, until it reaches a steady state. Thirdly, the liquid concentration profile in the DPM is not developed during the transient and is more reliable after a longer simulation time. So there is a subtle difference in the liquid concentration profile between the two during the initial transient phases.

Discussion and conclusion
We have studied wide, partly saturated shear bands in a split-bottom shear cell geometry using continuum theory and discrete particle simulations. Just as in experiments, the liquid content decreases in the shear band, and a peak of liquid concentration, located initially on the inflection point of the shear-rate profile, propagates away from the shear band, making the fluid-depleted region wider and drier with time. A simple diffusion-driven model for liquid transport explains this phenomenon with a variable, shear-rate-dependent diffusivity. Being diffusion driven, the peak velocity decreases exponentially with distance from the shear band, and thus no stationary state is reached in the tails of the shear band (figure 9b). We tracked the trajectory of the liquid concentration peak location from DPM and continuum theory and showed that the two numerical solutions are in good agreement.
By transforming the spatial coordinate, the continuum model is transformed from a shear-rate-dependent diffusion model to a model with constant diffusion and shear-rate-dependent drift. This approach is better to analyse since one-dimensional analytical solutions can be obtained for the liquid concentration profile, individually, for both drift and diffusion driven liquid migration. For short time intervals, the accumulated liquid concentration can be obtained from the superposition of the accumulated liquid concentrations of the drift and diffusion processes. Further, the velocities of propagation of the liquid concentration peak for the simplified model can be predicted as the sum of the velocities of propagation of liquid concentration peak of the drift and diffusion processes.
The relative dominance of the drift and diffusion related liquid transport can be understood via the local Péclet number. It is observed that diffusion dominates in the centre of the shear band, but both drift and diffusion are significant for the depletion of the shear band and the overall liquid transport in the tails. There is an inner region in the vicinity of the shear-band region where Pe 1 and diffusion dominates drift. Simultaneously, there is an adjoining outer region of Pe ≈ 1 where the effects of drift and diffusion become comparable (Pe ≈ 1). While drift enhances the drying of the shear band and accumulates the liquid in the peak, diffusion enhances the shift or transport of the liquid away from the shear band. In this way, both processes are important in determining the overall liquid migration features.
We calibrated the diffusivity and velocity profiles in the continuum model using DPM simulations. The results of the continuum model are comparable with the DPM simulations. The simulation time is as less as 30 minutes for the continuum model compared to 6 to 7 days for the corresponding DPM simulations of same system size, showing the effectiveness of a continuum approach. The comparison of the continuum model with the DPM shows that gravity plays a significant role in the liquid transport process too. However, this is not taken into account in the continuum model. Liquid migration in unsaturated granular media is a known phenomenon which is validated experimentally and numerically by other authors. Our contribution has been to simplify the continuum model and obtain analytical solutions for the simplified liquid migration model. The new tool is able to capture certain features such as the liquid concentration peak and the velocity of propagation of the peak within accuracies of 95 % and 80 %, respectively. Moreover, one can analyse the effects of drift and diffusion separately using the new tools for liquid migration. This gives us a clear understanding of the physics of the liquid migration process. Certain other features of the liquid concentration profile, e.g. the accumulated liquid concentration in the peak region, are also obtained by the new tool to within 60 % accuracy, but this prediction works for short time intervals. In the long term, the new tool over-predicts accumulated concentration due to the cumulative error in the prediction. The prediction time for the new tool is less than 4 minutes as compared to 30 minutes for the numerical simulations for the continuum models and 5 to 6 days for the expensive DPM simulations. Thus, the new mathematical tool is advantageous and useful for a quick prediction and understanding of the liquid migration process. Certain features of the liquid migration process are predicted well by this mathematical tool, e.g. the peak liquid concentration position is predicted to within 95 % accuracy and the velocity of liquid migration is predicted with up to 80 % accuracy. Other features, e.g. the accumulated liquid concentration, are well predicted for short time intervals at up to 60 % accuracy. However, one needs to use the full continuum model to predict the long-term behaviour of this feature.
Nevertheless, there are certain scopes for improvement of our approach in simulating the original continuum model. Firstly, our results on tracking the liquid concentration peak location showed no effect of grid size resolution. However, a high minimum resolution is necessary to have a smooth velocity profile. Thus, a numerical algorithm with adaptive spatial resolution would be able to solve the continuum model more accurately. Secondly, an effective contribution term due to gravity in the overall liquid transport flux is worth adding to the continuum model for better agreement with the DPM results. These will be taken into consideration in our future work. the peak locations is given aṡ The above equation is evaluated to findγ c s ≈ 0.36 s −1 .

Appendix B. Discretisation by FVM
The governing equation of diffusion can easily be derived from the general transport equation for property Q, the liquid concentration. Considering a transient state diffusion process in the x and z-directions in the full model, the equation is given by The key step of the FVM is the integration of the governing equation over the control volume. Thus, the above equation is integrated as and thus, where C liq A x (d(γ Q)/dx) and C liq A z (d(γ Q)/dz) are the diffusion fluxes and the subscripts E, W, N and S correspond to fluxes from the east, west, north and south directions, respectively, as illustrated in figure 16. Here, V denotes the volume of the domain,n denotes the normal to the surface, dV and dS denote the volume and surface area of the control volume, respectively; A x and A z represent the area across the east/west and north/south faces of the control volume, respectively. In case of diffusion in the simplified model, we have fluxes in the x-direction only and (B3) reduces to We approximate the diffusion fluxes with a simple first-order central difference scheme as where the subscript P corresponds to the amount of the quantity in each control volume. The discretised equation for each control volume is solved by semi-implicit method and is given by where, and Equation (B9) is solved semi-implicitly by using tridiagonal solution method in Matlab to obtain Q. We use a no-flux boundary condition in both the x and z-directions.
The velocity v 0 = dx 0 /dt thus satisfies Similarly, we can show Substituting (D5) and (D10a,b) into (D11), we obtain v tot This will now help us solving for the time derivative of φ(t), which satisfies Differentiating with respect to time, using Leibniz's rule, yields Substituting (D5) into ∂Q/∂t and (D11) into v 0 , we obtain Thus, we have shown that, for an incremental time step, v c andφ are additive. The fact that it works for longer time intervals is not guaranteed, but is observed, and this is a result of the paper.

Appendix E. Liquid migration model for DPM
In our present study, we use a liquid migration model as proposed by Mani et al. (2012Mani et al. ( , 2013. The methodology is quite straightforward according to the given reference: liquid is transferred locally whenever contacts are formed or broken. The particles and the liquid are considered as two different entities in the system. Liquid is either associated with a particle as a thin liquid film of volume V f i , or with a contact as a liquid bridge of volume V b ij . We describe the liquid migration model in the following subsections. E.1. Liquid bridge formation When two particles come into contact (i.e. overlap), a new liquid bridge is formed from the liquid contained in the particles' film. Since there can be some liquid of volume V min trapped in the roughness of the grains (Herminghaus 2005;Scheel et al. 2008), V f i must be larger or equal to V min . Therefore, the available liquid for bridge formation is V f i − V min . Usually, the length scale of roughness is small compared to the particle size and film volume V f i , so V min is often very small. Moreover, V min is fixed and trapped in the particles; thus, without loss of generality, we assume V min = 0 for our simulations. Thus, the bridge volume V b ij,new is given as where, V max is the maximum liquid bridge volume which is imposed in our simulations. We denote here the film volumes available on the interacting particles i and j as V f i and V j f , respectively, and liquid bridge volume as V b ij . The available liquid volume for bridge formation is the sum of the film volumes available on the individual particles V f i + V f j . This model is designed for small liquid contents and large contact angles with fast and easy transport of fluid on the surface. Figure 17 shows a schematic figure of liquid bridge formation.

A30-29
To avoid clustering of liquid by coalescence, the maximum bridge volume is restricted to V b ij < V max with the maximal V max = βr p 3 . The excess volume V max − V b ij remains as the film volume in the interacting particles. Our liquid bridge model is valid for small enough liquid contents only. Therefore, as a model simplification, we do not allow the formation of liquid clusters via coalescence by using the maximal V max of the bridge volume, which must not be exceeded. We use the capillary force model from Willett et al. (2000), which is limited to the maximal liquid bridge volume of V max = β. The appropriate value for the maximal liquid bridge volume β can be estimated by considering that for monodisperse spheres from Scheel et al. (2008) as β = 0.007.

E.2. Liquid bridge rupture
Bridge volumes are bound to contacts until they rupture. When the distance between two particles with a liquid bridge in between exceeds the rupture distance of the liquid bridge, the liquid bridge ruptures and the bridge volume is distributed to the neighbouring contacts where, p ∈ i, j and n ∈ neighbouring particles in contact and N c is the number of neighbouring contacts associated with the particles i and j. Figure 18 shows a schematic representation of liquid bridge rupture. The total liquid volume conservation is ensured in the system.

E.3. Coarse graining of the liquid distribution
We denote by V b ij the volume of the liquid bridge between particles i and j. The volume of the liquid film on particle i is denoted by V f i . We denote the Gaussian coarse-graining function with a width of one particle diameter by Φ(x, z) i , centred around the position of particle i. We denote the normalised integral of the Gaussian coarse-graining function along the branch vector between particles i and j by ψ(x, z) ij , see Weinhart et al. (2012a). Then the liquid concentration at position (x, z) is given by Weinhart et al. (2012b) for more details and the notation.