1. Introduction
Colloidal gels, formed through the aggregation of micron-sized particles within a solvent medium (Lekkerkerker et al. Reference Lekkerkerker, Poon, Pusey, Stroobants and Warren1992; Poon Reference Poon2002; Bergenholtz, Poon & Fuchs Reference Bergenholtz, Poon and Fuchs2003; Chen & Schweizer Reference Chen and Schweizer2004), represent a class of soft materials with diverse applications ranging from food and care products (Larson Reference Larson1999; Mezzenga et al. Reference Mezzenga, Schurtenberger, Burbidge and Michel2005) to crop protection formulations (Faers et al. Reference Faers, Choudhury, Lau, McAllister and Luckham2006), and bio-fluids (Darras et al. Reference Darras, Dasanna, John, Gompper, Kaestner, Fedosov and Wagner2022). In any real-world system, a density mismatch between the solvent and colloidal particles is unavoidable (Buscall & White Reference Buscall and White1987). This leads to buoyancy forces and subsequent sedimentation (or creaming) of the colloids, ultimately giving rise to a separation of the suspension into colloid-rich and colloid-poor regions. One of the appealing aspects of using colloidal gels is that they possess the ability to support the particles’ buoyant weight against gravity for a finite, often extended, duration (Allain, Cloitre & Wafra Reference Allain, Cloitre and Wafra1995; Senis, Gorre-Talini & Allain Reference Senis, Gorre-Talini and Allain2001; Starrs et al. Reference Starrs, Poon, Hibberd and Robins2002; Derec et al. Reference Derec, Senis, Talini and Allain2003; Manley et al. Reference Manley, Skotheim, Mahadevan and Weitz2005; Bartlett, Teece & Faers Reference Bartlett, Teece and Faers2012; Harich et al. Reference Harich, Blythe, Hermes, Zaccarelli, Sederman, Gladden and Poon2016; Padmanabhan & Zia Reference Padmanabhan and Zia2018). This time is typically referred to as a delay time and sets the typical shelf life of a gel-based product in many applied contexts (Harich et al. Reference Harich, Blythe, Hermes, Zaccarelli, Sederman, Gladden and Poon2016; Zhou Reference Zhou2018).
The gel's ability to support its own weight can be attributed to its network-like structure (Whitaker et al. Reference Whitaker, Varga, Hsiao, Solomon, Swan and Furst2019) and the internal dynamics (Zaccarelli & Poon Reference Zaccarelli and Poon2009). A colloidal gel forms due to the interplay between thermal/Brownian diffusion and strong, short-ranged attractions between the colloids. The latter cause the colloids to aggregate into a network structure, when the system is quenched into the spinodal region of the phase diagram (Carpineti & Giglio Reference Carpineti and Giglio1992). For attractions that are several times the thermal energy – 
$k_{{B}}T$, with 
$k_{{B}}$ Boltzmann's constant and 
$T$ the temperature – rearrangement of the formed network is slowed down (Foffi et al. Reference Foffi, Dawson, Buldyrev, Sciortino, Zaccarelli and Tartaglia2002). This significantly extends the time it takes the system to relax to thermodynamic equilibrium (Zaccarelli Reference Zaccarelli2007; Royall et al. Reference Royall, Faers, Fussell and Hallett2021). The arrested dynamics also halts (Harich et al. Reference Harich, Blythe, Hermes, Zaccarelli, Sederman, Gladden and Poon2016; Darras et al. Reference Darras, Dasanna, John, Gompper, Kaestner, Fedosov and Wagner2022) or at the very least strongly slows down (Derec et al. Reference Derec, Senis, Talini and Allain2003; Manley et al. Reference Manley, Skotheim, Mahadevan and Weitz2005) buoyancy-mediated separation. The amount by which this separation is suppressed depends sensitively on the way the gel is prepared and the interactions between the colloids. Consequently, the time scales reported for complete collapse can vary significantly, ranging from minutes to months (Harich et al. Reference Harich, Blythe, Hermes, Zaccarelli, Sederman, Gladden and Poon2016).
Accurately characterizing a gel's resistance to gravity in a physical model presents a considerable challenge. This is because experiments reveal that colloidal gels can exhibit several settling behaviours: fast sedimentation, delayed rapid collapse, and slow sedimentation (Harich et al. Reference Harich, Blythe, Hermes, Zaccarelli, Sederman, Gladden and Poon2016). Which of these is present in a given suspension is dependent primarily on the initial volume fraction 
$\phi _0$ of the sample, the strength of attraction between colloids 
$\epsilon _0$, and the size of the particles (Senis et al. Reference Senis, Gorre-Talini and Allain2001; Derec et al. Reference Derec, Senis, Talini and Allain2003; Manley et al. Reference Manley, Skotheim, Mahadevan and Weitz2005; Buzzaccaro et al. Reference Buzzaccaro, Secchi, Brambilla, Piazza and Cipelletti2012; Harich et al. Reference Harich, Blythe, Hermes, Zaccarelli, Sederman, Gladden and Poon2016; Darras et al. Reference Darras, Dasanna, John, Gompper, Kaestner, Fedosov and Wagner2022). The importance of size can be understood by considering the gravitational Péclet number, which gives the dimensionless ratio between sedimentation and diffusion. For a sphere of radius 
$a$, the number is given by
\begin{equation} {{Pe}}_{g} = \frac{4 {\rm \pi}g\,\Delta \rho\,a^{4}}{3 k_{{B}} T} , \end{equation}
where 
$g$ is the acceleration of gravity, and 
$\Delta \rho$ is the density difference between the colloid and the suspending fluid. The strong dependence on size follows from 
${{Pe}}_{g} \propto a^{4}$.
Harich et al. (Reference Harich, Blythe, Hermes, Zaccarelli, Sederman, Gladden and Poon2016) experimentally investigated the settling dynamics of depletion gels comprising 
${\sim }0.65\,\mathrm {\mu }{\rm m}$ colloids. They established a 
$\phi _{0}$–
$\epsilon _{0}$ state diagram for the types of behaviour that were observed. The most striking of these is delayed gravitational collapse (Harich et al. Reference Harich, Blythe, Hermes, Zaccarelli, Sederman, Gladden and Poon2016; Zhou Reference Zhou2018) that occurs for 
$0.15 \lesssim \phi _{0} \lesssim 0.35$. In these situations, the interface between the gel and supernatant phase – a region (nearly) devoid of particles – initially resists sedimentation. However, at a certain point, the system abruptly loses structural integrity and enters a regime of rapid settling, which is characterized by a constant velocity of the interface. This regime concludes when the interface reaches the densified bottom part of the sample, after which the velocity of the interface exponentially decays to zero as the colloid-rich phase compacts. It was suggested in Harich et al. (Reference Harich, Blythe, Hermes, Zaccarelli, Sederman, Gladden and Poon2016) and later shown in Zhou (Reference Zhou2018) that the accumulation of dense debris on top of the gel triggers the catastrophic failure. This debris originates from the curved parts of the meniscus between air and the sample.
Several attempts have been made to theoretically explain the response of colloidal gels to gravitational forces. In 2003, Derec et al. (Reference Derec, Senis, Talini and Allain2003) studied experimentally the rapid collapse of gels formed from strongly aggregating colloidal suspensions, for low values of the initial volume fraction (
$\phi _0<1\,\%$). They showed that in the first resistant regime, the gel–supernatant interface slowly settles at a constant velocity. It was later argued that this follows a power law 
$\phi _0^{(1-D)/(3-D)}$, with 
$D$ the fractal dimension of the gel network (Allain et al. Reference Allain, Cloitre and Wafra1995). In addition, Derec et al. (Reference Derec, Senis, Talini and Allain2003) connect the transition into the second, linear regime to the appearance of fractures within the bulk of the gel. This provides an easy route for the solvent to move up into the supernatant phase, which explains the abrupt increase in the interface's settling velocity. The authors also proposed a simple hydrodynamic model by which the increase in velocity can be estimated, which takes as input the height and radius of a typical crack in the gel.
These insights were built upon in the work of Manley et al. (Reference Manley, Skotheim, Mahadevan and Weitz2005), wherein the collapse dynamics is determined by the balance between gravitational stress and the yield stress of the network (
$\phi _0<1\,\%$). When the former is larger, a gel collapses poroelastically, with a rate of compression that decays exponentially in time. Otherwise, the network eventually yields, leading to rapid linear settling. The authors use a modified Carman–Kozeny relation (Happel & Brenner Reference Happel and Brenner1991) to argue that the characteristic pore size in the gel network is set by the largest length scale in the system, the cluster size. They also show that their experimental data collapse onto a single curve that is well described by their model. However, it is important to note that both studies (Derec et al. Reference Derec, Senis, Talini and Allain2003; Manley et al. Reference Manley, Skotheim, Mahadevan and Weitz2005) focus primarily on colloidal particles with small diameters, typically in the range of tens of nanometres, which can explain some of the differences between the observations in Harich et al. (Reference Harich, Blythe, Hermes, Zaccarelli, Sederman, Gladden and Poon2016) and Zhou (Reference Zhou2018).
More recently, Darras et al. (Reference Darras, Dasanna, John, Gompper, Kaestner, Fedosov and Wagner2022) investigated the sedimentation behaviour of dense soft colloidal gels, formed by the aggregation of red blood cells (
$\phi _{0} \gtrsim 0.25$). Unlike previous models that assumed low initial volume fractions to simplify hydrodynamic descriptions, they developed a theoretical framework capable of accurately reproducing the time evolution of the gel interface height. Their model operates on the assumption that the gel network uniformly compresses as it sediments, resulting in a mostly homogeneous colloid distribution over time. This is a reasonable approximation for the dense and weak gels that they examined. However, while their model successfully captures many aspects of the gel sedimentation process, the initial delay time observed in their gels does not arise spontaneously from the theory and is instead introduced as an adjustable parameter.
The study of collapsing colloidal gels remains an active area of research, particularly in understanding how gravity influences gel structure. Gallegos et al. (Reference Gallegos, Martínez-Rivera, Valadez-Pérez and Castañeda-Priego2023) explored gravitational effects in patchy colloidal systems, identifying a threshold 
${{Pe}}_{g}$ above which gravity enhances particle bonding, promoting clustering and altering the gel network's structure. This leads to the coexistence of dilute and dense phases, or even a crystalline state, depending on 
${{Pe}}_{g}$. Tateno, Wang & Tanaka (Reference Tateno, Wang and Tanaka2024) used confocal microscopy to study gravity-induced collapse at the single-particle level, showing that gel microstructure is determined solely by local volume fraction, while compressive stress forms chain-like structures that support external stress. These insights underscore the role of gravity in gel collapse and help to inform the development of models addressing such phenomena.
Motivated by the absence of theoretical descriptions that capture a delay time, we present here a comprehensive theoretical framework that overcomes this issue. We also account for spatial variation of the gel properties, which should allow us to shed light on phenomena like debris formation. In our model, we describe the stress within the gel as a dilatational viscous response to applied compression (Landau & Lifshitz Reference Landau and Lifshitz1959). Following other analyses, Darcy's law is used to model the dynamic coupling between the solvent and colloidal phases (Carman Reference Carman1939; Gray & Miller Reference Gray and Miller2004; Carrillo & Bourg Reference Carrillo and Bourg2019). The average network pore cross-section is determined by the average surface-to-surface distance between particles in our approach. We show that these elements are sufficient to qualitatively capture features observed in experiments, including the emergence of delay times and the accumulation of debris at the top of the gel.
The remainder of this paper is organized as follows. In § 2, we present the theoretical framework in generic dimensions and boundary conditions. In § 3, we first provide the analytical solution for the instantaneous mean colloidal velocity in samples with homogeneous volume fractions. Subsequently, numerical results depicting the time evolution of the gel–supernatant interface, as well as colloidal density profiles, are presented. We conclude the section by demonstrating how the delay time can be predicted accurately from the characteristics of colloidal flow at the beginning of the collapse. In § 4, we delve into a detailed discussion of the obtained results, addressing the limitations of the models and establishing connections with experimental observations. Finally, in § 5, we draw conclusions and explore potential future directions, including the study of meniscus curvature on the behaviour of a gel.
2. Theoretical model
In this section, we provide a step-by-step derivation of our model. As the primary focus of this paper is the theoretical result, we have included details on the numerical solving of the model in Appendix A.
2.1. Mass conservation equation
Consider a binary mixture of a solid colloidal phase, with volume fraction 
$\phi _{c} \equiv \phi$, and a liquid phase with volume fraction 
$\phi _{l}$. This composition can vary throughout space and time, but we will leave the dependence on position 
$\boldsymbol {r}$ and time 
$t$ of our fields implicit throughout to ease the notation. Volume conservation implies that 
$\phi _{l} = 1 - \phi$. The phases have densities 
$\rho _{c}$ and 
$\rho _{l}$, respectively. Imposing mass conservation for each phase, we obtain the integral expressions
$$\begin{gather} \int_{V} \mathrm{d}V \, \partial_{t} \left[ \rho_{l} (1-\phi) \right] ={-} \int_{\partial V} \mathrm{d}\boldsymbol{A} \boldsymbol{\cdot} \left[ \rho_{l} (1-\phi) \boldsymbol{v}_{l} \right] , \end{gather}$$
$$\begin{gather}\int_{V} \mathrm{d}V \, \partial_{t} \left( \rho_{c} \phi \right) ={-} \int_{\partial V} \mathrm{d}\boldsymbol{A} \boldsymbol{\cdot} \left( \rho_{c} \phi \boldsymbol{v}_{c} \right) . \end{gather}$$
Here, 
$\boldsymbol {v}_{c}$ and 
$\boldsymbol {v}_{l}$ are the colloid and solvent velocities, respectively, ‘
$\boldsymbol {\cdot }$’ indicates the inner product, and 
$\partial _{t}$ is the partial derivative with respect to time. Integrals on the left-hand side are over the volume of a stationary control volume 
$V$, while the right-hand sides integrate the fluxes through that volume's surface 
$\partial V$. Since the control volume is in principle arbitrary, (2.1) and (2.2) give rise to the differential form
$$\begin{gather} \partial_{t} \phi ={-}\boldsymbol{\nabla} \boldsymbol{\cdot} \left( \phi \boldsymbol{v}_{c} \right) , \end{gather}$$
$$\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \left[\phi \boldsymbol{v}_{c} + (1-\phi) \boldsymbol{v}_{l} \right] = 0. \end{gather}$$
Equation (2.3) implies that 
$\phi$ changes only via flow-mediated flux through the boundary, and we will use it to integrate 
$\phi$ in time. Equation (2.4) defines a compressibility criterion on a volume-fraction-weighted flow.
2.2. Momentum balance equation
Next, we derive an expression for the velocity of the colloidal phase, relative to the liquid phase velocity. First, consider the momentum balance equation for the entire binary mixture:
\begin{equation} \frac{{\rm d}}{{\rm d}t} \int_{V} \mathrm{d}V \left[ \rho_{{mix}} \, \boldsymbol{v}_{{mix}} \right] = \int_{V} \mathrm{d}V \left[\rho_{{mix}} \, \boldsymbol{g} + \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\sigma}_{{mix}}\right] . \end{equation}
On the left-hand side, we take the total time derivative of the momentum density, which follows from integrating the total density 
$\rho _{{mix}}(\phi ) = \rho _c \phi + \rho _{l} (1-\phi )$ together with the velocity of the mixture 
$\boldsymbol {v}_{{mix}}$. By adding equations (2.1) and (2.2), we can define a mixture flux 
$\boldsymbol {J}_{mix} = \rho _{mix}(\phi )\,\boldsymbol {v}_{mix}$, and subsequently write the mixture velocity as
\begin{equation} \boldsymbol{v}_{{mix}} = \frac{\rho_{c} \phi \boldsymbol{v}_{c} + \rho_{l} (1-\phi) \boldsymbol{v}_{l} }{ \rho_{{mix}}(\phi) } , \end{equation}
which represents a weighted average by the densities of the two phases. On the right-hand side of (2.5), momentum is introduced into the system by the gravitational acceleration 
$\boldsymbol {g} = (0,0,-g)$, and internally redistributed via (the divergence of) the mixture stress tensor 
$\boldsymbol {\sigma }_{{mix}}$. We define the latter as a sum over a Newtonian term from the liquid solvent, and a non-Newtonian stress 
$\boldsymbol {\sigma }_c$ that is related to the presence of the solid phase in the mixture:
\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\sigma}_{{mix}} ={-}\boldsymbol{\nabla} p + \mu \, \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\mathsf{e}} + \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\sigma}_c. \end{equation}
Here, 
$p$ and 
$\mu$ are the liquid hydrostatic pressure and dynamic viscosity, respectively, and 
$\boldsymbol{\mathsf{e}} = 1/2 [\boldsymbol {\nabla } \boldsymbol {v}_l + (\boldsymbol {\nabla } \boldsymbol {v}_l)^{\rm T}]$ represents the rate-of-strain tensor.
To make progress, we now consider the momentum balance equation for each phase separately. Introducing the momentum exchange terms 
$\boldsymbol {\varSigma }_{ij}$ between phases 
$i,j \in \{l,c\}$, we arrive at
$$\begin{gather} \frac{{\rm d}}{{\rm d}t} \int_{V} \mathrm{d}V \left( \rho_c \phi \boldsymbol{v}_c \right) = \int_{V} \mathrm{d}V \left(\rho_c \phi \boldsymbol{g} + \boldsymbol{\varSigma}_{cc} + \boldsymbol{\varSigma}_{cl}\right), \end{gather}$$
$$\begin{gather}\frac{{\rm d}}{{\rm d}t} \int_{V} \mathrm{d}V \left[ \rho_l (1-\phi) \boldsymbol{v}_l \right] = \int_{V} \mathrm{d}V \left[\rho_l (1-\phi) \boldsymbol{g} + \boldsymbol{\varSigma}_{ll} + \boldsymbol{\varSigma}_{lc}\right]. \end{gather}$$
Clearly, Newton's third law implies 
$\boldsymbol {\varSigma }_{cl} = - \boldsymbol {\varSigma }_{lc}$. Summing the two individual contributions above and subtracting (2.5), we obtain
\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot} \left[ - \mathbb{I} p + \mu \boldsymbol{\mathsf{e}} + \boldsymbol{\sigma}_c + \rho_c \rho_l\,\frac{\phi(1-\phi)}{\rho_m(\phi)}\,\boldsymbol{\delta v}\,\boldsymbol{\delta v} \right] = \boldsymbol{\varSigma}_{ll} + \boldsymbol{\varSigma}_{cc} , \end{equation}
with 
$\mathbb {I}$ the identity operator, and 
$\boldsymbol {\delta v} \equiv \boldsymbol {v}_c - \boldsymbol {v}_l$ the net velocity of the colloids with respect to the background liquid flow. Evaluating the above equation in the absence of a colloidal phase, we can split the momentum exchange terms as
$$\begin{gather} \boldsymbol{\varSigma}_{cc} = \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\tilde{\sigma}}_c, \end{gather}$$
$$\begin{gather}\boldsymbol{\varSigma}_{ll} = \boldsymbol{\nabla} \boldsymbol{\cdot} (-\mathbb{I}p+ \mu \boldsymbol{\mathsf{e}}), \end{gather}$$where we have redefined the non-Newtonian stress
\begin{equation} \boldsymbol{\tilde{\sigma}}_c = \boldsymbol{\sigma}_c + \frac{\rho_c \rho_l \phi(1-\phi)}{\rho_{mix}}\,\boldsymbol{\delta v}\,\boldsymbol{\delta v} , \end{equation}to ease the notation.
Finally, we choose to break up the cross-terms into a hydrostatic contribution and a dynamical (friction) term,
\begin{equation} \boldsymbol{\varSigma}_{cl} ={-}\phi\,\boldsymbol{\nabla} p + \boldsymbol{F}_{{Darcy}} , \end{equation}
the latter of which we choose to model using (a variant of) Darcy's law (Carman Reference Carman1939; Gray & Miller Reference Gray and Miller2004; Carrillo & Bourg Reference Carrillo and Bourg2019). This law states that 
$\boldsymbol {F}_{{Darcy}}$ is proportional to the relative velocities between the two phases, and it represents a coarse-grained contribution of all the microscale dissipative dynamics between the colloids and the solvent. Here, we write
\begin{equation} \boldsymbol{F}_{{Darcy}} = \frac{\mu}{\sigma^2}\,\frac{(1-\phi)(\boldsymbol{v}_{l} - \boldsymbol{v}_{c})}{K(\phi)}. \end{equation}
The dimension-free scalar function 
$K(\phi )$ represents the porosity of the mixture, 
$\sigma = 2a$ the (mean) particle diameter, and 
$\mu$ the viscosity of the suspending fluid, which ensures the correct dimensionality. Note that we must have that 
$K$ diverges in the absence of colloids, such that the Darcy-like term disappears from the equation of motion. Conversely, 
$K$ must go to zero as 
$\phi$ approaches the maximum value allowed, 
$\phi _m$ – depending on the context, this could be random loose packing, random close packing, or another dense arrested state. We will provide an explicit form for 
$K(\phi )$ in § 2.4.
By defining a phase-related convective derivative 
$D^i_t= \partial _t +\boldsymbol {v}_i \boldsymbol {\cdot }\boldsymbol {\nabla }$, and combining (2.8) with (2.10) and (2.12), we can eliminate the pressure 
$p$ from the system. This allows us to arrive at our final relation:
\begin{equation} D^c_t[\rho_c \boldsymbol{v}_c] - D^l_t[\rho_l \boldsymbol{v}_l] + \frac{ \mu\,\boldsymbol{\delta v}} {\sigma^2 \phi K} = \Delta\rho\,\boldsymbol{g} + \frac{\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\tilde{\sigma}}_c}{\phi} - \frac{\mu\,\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\mathsf{e}}}{1-\phi} . \end{equation}2.3. Quasi-hydrostatic limit
We perform a dimensional analysis on (2.14), and conclude that there it is useful to introduce a characteristic velocity 
$v_g$ and length 
$l_R$. The former represents the bulk settling velocity of a single colloid due to buoyant forces, 
$v_g = \Delta \rho \,g \sigma ^2/ \mu$. The latter is the length scale associated with smallest distance that the coarse-grained theory can resolve. That is, 
$l_R^3$ is the representative elementary volume (REV) (Carrillo & Bourg Reference Carrillo and Bourg2019), which is the smallest volume containing sufficient particles to define a local field 
$\phi$. This parameter is system-dependent and can be used to demarcate the scale transition from microscale to macroscale descriptions in porous materials (Gray & Miller Reference Gray and Miller2004).
From our dimensional analysis, we conclude that (2.14) can be rewritten as
\begin{equation} \lim_{{{Re}} \rightarrow 0} \, \frac{ \boldsymbol{\delta v}^*} {K(\phi )} ={-}\phi \hat{\boldsymbol{z}} + \alpha^2 \left[ \boldsymbol{\nabla}^* \boldsymbol{\cdot} \boldsymbol{\sigma}^*_c - \frac{ \phi\,\boldsymbol{\nabla}^* \boldsymbol{\cdot} \boldsymbol{\mathsf{e}}^*}{1-\phi} \right ]. \end{equation}
Here, we have introduced the dimension-free constants 
$\alpha = \sigma / l_R$ – the ratio of the particle diameter to the REV size – and 
${{Re}} = \Delta \rho \,v_g l_R / \mu$ – the Reynolds number. Quantities denoted with an 
$\ast$ superscript are expressed in natural units. The limit 
${{Re}} \rightarrow 0$ emerges naturally, since in real systems this limit holds for the small particles suspended in a viscous medium.
We can further simplify the expression by considering the limit as 
$\alpha \rightarrow 0$. This represents highly correlated systems, where the amount of particles needed to be described at a coarse-grained level is extremely large. However, caution is advised in following this approach. For finite-size systems of typical length 
$L$, 
$\alpha$ is naturally bounded from below by the ratio 
$\sigma /L$. Lower values would make the coarse-grained description meaningless. Moreover, the non-Newtonian contribution to the stress 
$\sigma _c$ may be singular (divergent), particularly at the maximum packing fraction 
$\phi _{m}$. This, for example, prevents the non-physical overlap of hard colloids. This makes taking the limit 
$\alpha \rightarrow 0$ poorly defined in some cases.
The limit wherein 
$\alpha$ and 
${{Re}}$ both tend to zero is commonly referred to as the quasi-hydrostatic regime. Assuming small values of 
$\alpha$ and 
${{Re}}$, we can express
\begin{equation} \boldsymbol{\delta v} \approx K(\phi )\,({-}v_g \phi \hat{\boldsymbol{z}} + \sigma^2 \mu^{{-}1}\, \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{{\sigma}}_c). \end{equation}
At this point, we need to specify the forms of 
$K(\phi )$ and 
$\boldsymbol {\sigma }_c$ in order to solve the problem. This also requires us to combine (2.16) with (2.4) to eliminate the dependence on the liquid velocity, as we will show in § 3.
2.4. Mean porosity
The distribution of pore cross-sections within the colloidal phase is intricately linked to the distribution of particle surface-to-surface distances. This relationship is typically hard to predict (Carman Reference Carman1939; Heijs & Lowe Reference Heijs and Lowe1995; Xu & Yu Reference Xu and Yu2008; Ozgumus, Mobedi & Ozkol Reference Ozgumus, Mobedi and Ozkol2014), as it is influenced by factors such as system temperature, details of interparticle interactions (Ruiz-Franco et al. Reference Ruiz-Franco, Camerin, Gnan and Zaccarelli2020), and potentially, the system's preparation/rheological history (Koumakis et al. Reference Koumakis, Moghimi, Besseling, Poon, Brady and Petekidis2015; Gibaud et al. Reference Gibaud, Dagès, Lidon, Jung, Ahouré, Sztucki, Poulesquen, Hengl, Pignon and Manneville2020; Saint-Michel, Petekidis & Garbin Reference Saint-Michel, Petekidis and Garbin2022; Torre & de Graaf Reference Torre and de Graaf2023). For simplicity, and as a first-order approximation, we opt to characterize the mean pore cross-section using the ansatz
\begin{equation} K(\phi ) = \frac{(1-{\phi} / {\phi_{m}})^3}{k^2_0 \phi} \end{equation}
in this work. This expression adequately reproduces the findings reported in Torquato (Reference Torquato1995), where the author estimates the nearest-neighbour distance in a homogeneous configuration of hard spheres as a function of the volume fraction 
$\phi$. Thus we model the mean pore cross-section simply as the mean surface-to-surface distance between spherical particles. The constant parameter 
$k_0 = 3\sqrt {2}$ is determined by computing the sedimentation velocity of a single sphere in a viscous unbounded fluid. To achieve this, we evaluate (2.16) for 
$\phi =0$, assuming that 
$[K(\phi )\,\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {{\sigma }}_c] |_{\phi =0} = 0$. Note that although 
$\boldsymbol {v_c}$ is non-zero, where the volume fraction vanishes, the colloid flux 
$J_c=\phi \boldsymbol {v_c}$ is zero, as required.
2.5. Dilatational viscosity
We choose to represent the non-Newtonian stress as a simple diagonal operator, which depends on the local 
$\phi$ and the amount of compression in the colloid field:
\begin{equation} \boldsymbol{\sigma}_c = \lambda(\phi)\,\mathbb{I} \left( \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{v}_c \right). \end{equation}
Here, 
$\lambda$ denotes the dilatational viscosity, which physically represents the irreversible resistance to compression or expansion within a fluid. At the microscopic level, it arises from the finite time needed for energy injected into the system to disperse among the various degrees of freedom of colloidal motion (Landau & Lifshitz Reference Landau and Lifshitz1959). Note that (2.18) does not account for shear stress response of the colloidal phase. However, in the present study we limit our focus to one-dimensional (1-D) systems, where the shear component of the stress is not resolved.
We model the dilatational viscosity as the product of a local energy density and relaxation time, 
$\lambda (\phi ) = \epsilon (\phi )\,T_{{r}}(\phi )$. The first can be expressed as the product of particle density 
$n_c = N_c/V$ and the average number of bonds per particle 
$n_b(\phi )$. To estimate the latter, we compute the ratio of incoming 
$j_{{in}}=D_0 / \delta r^2(\phi )$ and outgoing 
$j_{{out}}=D_0 / \varGamma ^2$ fluxes of particles for a single central colloid. Here, 
$D_0 = (k_{{B}} T) / (3 {\rm \pi}\mu \sigma )$ represents the diffusion coefficient of an isolated colloid in an unbounded fluid, and 
$\varGamma \gtrsim \sigma$ denotes the range of the attractive interactions, as measured from the colloid centres. We retain only the fraction of particles that do not escape from the potential well of the central colloid, resulting in
\begin{equation} \epsilon(\phi) = \epsilon_0 n_c\,n_b(\phi) = \epsilon_0\,\frac{3 z_{{max}} }{{\rm \pi} \sigma^3}\,\frac{ \phi \varGamma^2}{\delta r^2(\phi)}\,(1-{\rm e}^{-\epsilon_0 / k_B T }). \end{equation}
Here, 
$z_{{max}}$ represents the maximum number of bonds that a particle can form due to geometric constraints, and 
$\epsilon _0$ is the energy of a single bond. Additionally, and similarly to what we assumed for the mean porosity, 
$\delta r = \sigma [1+ (1-\phi /\phi _m)^{3/2}\phi ^{-1/2}/10 ]$ denotes the mean distance to the closest neighbouring particle in a spatially homogeneous distribution of colloids, approximately reproducing Torquato's original result (Torquato Reference Torquato1995).
Next, we calculate the relaxation time of the colloidal phase. To obtain an estimate, we use an effective diffusion coefficient 
$D(\phi )=D_0\,f(\phi )$ and expand it to first order around 
$\phi =\phi _{m}$ where diffusion is suppressed by crowding, thus yielding 
$D |_{\phi =\phi _m} = 0$:
\begin{equation} T_{{r}}(\phi) =\frac{l^2_R}{D_0 (1-\phi/\phi_m)}. \end{equation}Now we are in a position to combine (2.19) and (2.20) to obtain the dilatational viscosity
\begin{equation} \lambda(\phi) =\mu\,\frac{9 \alpha^{{-}2} z_{{max}} \ U(1-{\rm e}^{{-}U})}{\left[\phi^{1/2}+\tfrac{1}{10} (1-\phi/\phi_m)^{3/2}\right]^2}\,\frac{ \phi^2 }{(1-\phi/\phi_m)}, \end{equation}
where we introduce the dimension-free potential strength 
$U=\epsilon _0 / k_B T$, and assume short-range interactions (
$\varGamma \approx \sigma$). As expected, the stress exhibits a singularity at the maximum packing fraction 
$\phi _m$. This implies that it can serve as a means to dissipate kinetic energy at a rate large enough to suppress applied stresses. In turn, this would enable the (dense) gel to support its own weight.
2.6. Helmholtz decomposition
The above framework is limited to 1-D systems. When extending our analysis to higher dimensions, an additional relationship involving the spatial derivative of the colloid velocity becomes necessary to close the system. One possible approach is to perform a Helmholtz decomposition (Ribeiro, de Campos Velho & Lopes Reference Ribeiro, de Campos Velho and Lopes2016; Klaseboer, Sun & Chan Reference Klaseboer, Sun and Chan2019; Glötzl & Richters Reference Glötzl and Richters2023) on the mixture velocity field. That is, we can express it as
\begin{equation} \boldsymbol{v}_{{mix}} = \boldsymbol{\omega} + \boldsymbol{\zeta}, \end{equation}
where the vector field is divided into a divergence-free component 
$\boldsymbol {\omega }$, and a curl-free component 
$\boldsymbol {\zeta }$. A natural choice for the former is derived from (2.4), namely 
$\boldsymbol {\omega } = \phi \boldsymbol {v}_{c} + (1-\phi ) \boldsymbol {v}_{l}$. Consequently, we can explicitly compute 
$\boldsymbol {\zeta } = \boldsymbol {v}_{{mix}} - \boldsymbol {\omega }$, and obtain the constraint
\begin{equation} \boldsymbol{\nabla} \times \left [ \frac{\Delta \rho}{\rho_{{mix}}}\,\phi (1-\phi)\, \boldsymbol{\delta v} \right ] = 0. \end{equation}
It is worth noting that this decomposition is not always unique and may not be applicable in all scenarios (Glötzl & Richters Reference Glötzl and Richters2023). This is because the smoothness of the mixture velocity field 
$\boldsymbol {v}_{{mix}}$ is a prerequisite for employing this technique.
3. Results
In this section, we explore solutions of our model specifically tailored to 1-D systems. We start our analysis by deriving the exact solution for the initial (
$t = 0$) colloid velocity field 
$\boldsymbol {v}_c$, when there is a homogeneous density profile. This already provides key insights into the origin of the gel's strength, and how it is able to support its own weight (for a finite time). Subsequently, we consider the full time evolution numerically. This clearly demonstrates the gel's resistance to gravitational collapse, and reproduces gel settling regimes qualitatively.
3.1. Homogeneous colloid distribution
Consider a homogeneous isolated system with net volume 
$HL^2$ and height 
$H \ll L$, as well as initial colloid volume fraction 
$\phi (z) |_{t=0} = \phi _0$. Due to the underlying symmetry, the velocity field can be expressed as a function of the 
$z$-coordinate and time only, 
$\boldsymbol {v}_c = v_c(z,t)\,\hat {\boldsymbol {z}}$. Combining (2.4) with (2.16), and defining a characteristic length
\begin{equation} \zeta= \sigma \sqrt{K(\phi_0)\,\lambda(\phi_0)\,(1-\phi_0) / \mu} , \end{equation}we obtain the following equation for the initial colloid velocity field:
\begin{equation} \partial_z {v}_c |_{t=0} = \zeta^2 \, \partial^3_z {v}_c |_{t=0} . \end{equation}
Equation (3.2), when extended to arbitrary dimensions, reveals its nature as a Poisson equation describing the stress 
$\boldsymbol {\sigma }_c$, with a source term linked directly to the gel's compression, expressed as 
$\phi \,\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {v_c}$.
Using mass conservation, we integrate (2.3) in space to obtain the boundary conditions
$$\begin{gather} v_c |^{z=0,H}_{t=0}=0 , \end{gather}$$
$$\begin{gather}\partial^2_z v_c |^{z=0,H}_{t=0} = \frac{v_g \phi_0}{\lambda(\phi_0)}. \end{gather}$$
We can also define a characteristic colloid velocity 
$V_g = v_g \phi _0 (1-\phi _0)\,K(\phi _0)$, representing the colloidal flow in absence of non-Newtonian stress. The presence of a characteristic length scale 
$\zeta$ and velocity 
$V_{g}$ allows us to define a characteristic time 
$\tau _d = \zeta / V_g$. In the next subsection, we will show that the latter can be used to obtain an accurate estimation of the gel collapse delay time.
Finally, combining the above results, the initial colloid velocity can be written as
\begin{equation} {v}_c(z)|_{t=0} = V_g \left [ \frac{\cosh{\dfrac{H/2 - z}{\zeta}}}{\cosh{\dfrac{H}{2\zeta}}} - 1 \right ]. \end{equation}
Here, we recognize 
$\zeta$ as a stress screening length for the gel. In the limit 
$\zeta / H \gg 1$, a strong resistance to compression and expansion is extended to the entire sample, halting the sedimentation process. This can be seen readily by Taylor expanding (3.5) with respect to the small parameter 
$H/(2\zeta )$, and noting that 
${v}_c(z)|_{t=0}$ scales quadratically with it. Note that a similar concept of stress screening length was already introduced in Evans & Starrs (Reference Evans and Starrs2002).
3.2. Time evolution
For our numerical analysis, we set the height of our sample as 
$H=5 \times 10^{-2}\,\mathrm {m}$, and the bulk settling velocity of a single colloid as 
$v_g = 10^{-5}\,\mathrm {m}\,\mathrm {s}^{-1}$, reproducing commonly used experimental set-ups (Senis et al. Reference Senis, Gorre-Talini and Allain2001; Derec et al. Reference Derec, Senis, Talini and Allain2003; Manley et al. Reference Manley, Skotheim, Mahadevan and Weitz2005; Harich et al. Reference Harich, Blythe, Hermes, Zaccarelli, Sederman, Gladden and Poon2016). We also set 
$l_R \approx 1.3 \times 10^{-4}\,{\rm m}$, and incrementally discretize height using a spacing 
$\delta z \in \{ 10^{-3}, 5 \times 10^{-4}, 2.5 \times 10^{-4} \}\,H$, where the smaller values are used for more numerically challenging systems. For the initial volume fraction, we choose 
$\phi _0 \in [0.01, 0.55]$, which spans the range of experimentally considered densities (Harich et al. Reference Harich, Blythe, Hermes, Zaccarelli, Sederman, Gladden and Poon2016; Darras et al. Reference Darras, Dasanna, John, Gompper, Kaestner, Fedosov and Wagner2022), and set 
$\phi _m=0.7$. Finally, we take the (reduced by 
$k_{{B}}T$) potential strength in the range 
$U \in [5, 50]$. Here, 
$U \approx 5$ is a typical lower bound for colloidal gelation (Harich et al. Reference Harich, Blythe, Hermes, Zaccarelli, Sederman, Gladden and Poon2016), while 
$U \approx 50$ represents a very strong gel at intermediate to high 
$\phi _{0}$. With these choices, a single colloid would sediment the entire sample height in a time 
$t_s = 5 \times 10^{3}\,{\rm s}$.
We numerically solve (2.4) and (2.3) for a 1-D system with these choices, to obtain the colloidal density profile as a function of time; see Appendix A for more information. For these solutions, we define the height of the gel interface 
$h(t)$ – between the colloid-rich and colloid-poor parts of the sample – as the tallest point in the system where the local colloid volume fraction 
$\phi (z)$ exceeds 
$\phi _0/2$.
Figure 1 shows 
$h(t)$ for several choices of 
$\phi _0$ and 
$U$. We observe three regimes: an initial delay, a linear ‘collapse’, and an exponential compaction. The former is hard to distinguish for low 
$\phi _0$, while the latter cannot be observed for high 
$\phi _0$, due to the slow settling dynamics for such 
$\phi _0$. Examining the delay more closely, we find that it is significantly influenced by both 
$\phi _0$ and 
$U$, which reflects the monotonic increase of the dilatational viscosity 
$\lambda$ with these parameters. An analytical expression for the delay time will be provided later in this section, in (3.6). In the linear regime, the interface settles at a constant velocity 
$(1-\phi _0/\phi _m)^{-1}V_g$, determined by the porous structure of the gel medium rather than its resistance to deformations. The onset of the third regime occurs when the interface reaches the dense region formed at the bottom of the sample, with the settling of the (largest part of the) gel bulk. In this regime, sedimentation speed is governed predominantly by colloid-solvent hydrodynamic drag, with minimal deviations observed among systems with different potential strengths 
$U$.
Time evolution of the interface height between the settling gel and supernatant region. Lighter colours denote less dense configurations, while darker colours represent denser ones, as labelled, with the symbols indicating various bond strengths; see the legend. The interface height is adjusted by its long-term limit 
$h^{\infty } = H \phi _0 / \phi _m$ (complete separation) and scaled to fit within the range 
$[0,1]$. Dashed black lines are guides to the eye to indicate a linear trend with slope 
$(1-\phi _0/\phi _m)^{-1}\,V_g(\phi _0)$. Time is expressed in units of 
$t_s = H / v_g$, representing the time required for a single colloid to sediment the height of the sample volume 
$H$.
Next, we connect the (early time) behaviour of the interface to the gel's internal structure. Figure 2(a) shows the time evolution of density profiles for a sedimenting gel with 
$\phi _0=0.25$ and 
$U = 30$. In the linear settling regime, the bulk of the gel sediments is almost unperturbed. That is, there is a large section of the gel that remains at 
$\phi _0$. We also see a small (physical) peak at the upper end of this flat range, which represents the accumulation of ‘debris’ on top of the gel. We will see that this originates from the early stages of settling when we turn to figure 2(b).
Density profiles depicting the time evolution of a sedimenting gel with an initial volume fraction 
$\phi _0=0.25$ and dimension-free potential strength 
$U = 30$. Local volume fractions 
$\phi (z)$ are normalized by the maximum volume fraction 
$\phi _m$ (left) or their initial value 
$\phi _0$ (right), with colours indicating temporal evolution from light to dark. The evolution direction is indicated by the black arrow. (a) Colloid density sampled uniformly in time as a function of the reduced height 
$z/H$, with the time interval between each curve given by 
$0.3 t_s$. The initial configuration is represented by a solid line, the linear regime by dashed lines, and the exponential compaction by dotted lines. (b) The flux of colloids 
$J_c = v_c \phi$ (red), normalized by the bare gravitational flux 
$J_g = V_g \phi _0$, is shown alongside normalized local volume fraction (green), plotted against height shifted by the system size 
$H$, and expressed in units of the gel stress screening length 
$\zeta$. Thick lines and symbols denote the system at three representative times: the initial configuration (circles), the creation of debris at the top of the sample (triangles), and during linear collapse (squares). The thin lines give an impression of the behaviour at intermediate times 
$\{0.02,0.04,0.07,0.09,0.13,0.15,0.5, 0.75 \}\,t_s$, and the characteristic time for this initial volume fraction is 
$\tau _d \approx 0.09 t_s$.
Turning our attention to the bottom of the sample, we see that here, the gel begins to compact immediately. That is, for our parameter choices, the gel cannot fully support its own weight; we will return to this in § 4, when we discuss the connection to experiment. The transition from linear collapse to exponential settling occurs at 
$t \approx 2.1 t_{s}$, and can be seen to coincide with the dense region at the bottom meeting with the interface between the gel and supernatant.
In figure 2(b), early-time colloidal flux 
$J_c (z) = \phi \,v_c (z)$ and 
$\phi (z)$ are depicted. As the bulk of the system sediments, it stretches the top part and compresses the bottom part of the gel – referencing (3.5), the top and bottom can be identified to be within 
$\zeta$ of either boundary. The top of the gel tries to resist this stretch until the thinnest part connecting it with the bulk breaks. Unsurprisingly, this is approximately when the bulk has settled at a distance 
$\zeta$. This leaves the top part detached from the rest of the gel. Note that here we did not explicitly put in an attraction between the top of the sample and the gel; the apparent attachment to the top of the gel is caused by the no-flux boundary conditions imposed on the colloid velocity field 
$v_c$. These are necessary to ensure that the system is isolated.
Throughout the evolution of the system, the part that remains at the top decreases in density and creates debris. This will sediment at a faster rate than the denser bulk gel due to the weaker hydrodynamic dissipation that this debris experiences. Ultimately, the debris deposits itself on top of the bulk as it continues to sediment, and causes a local density peak. As this denser material resists settling more, a local density minimum eventually forms in front of it. Revisiting figure 1, we understand that the delayed collapse can be attributed to the top part of the gel resisting collapse for a finite amount of time. Considering 
$J_c(z)$ in figure 2(b), the top layer detaches from the gel bulk, the flux of colloids transforming from a monotonically increasing function to a function with local minima and maxima. The former indicates that the gel is swelling everywhere near the top boundary in a cohesive way. The latter signals that a part of the gel has detached from the bulk and is accumulating on top of the homogeneously collapsing bulk.
Finally, in figure 3 we plot the delay time 
$T_d$ obtained from numerical solutions – defined as the earliest time for which 
$h(t)< H$ – as a function of initial volume fraction 
$\phi _0$. The delay time is scaled by 
$c(U)=\sqrt {U(1-{\rm e}^{-U})}$ to ensure data collapse. This scaling implies that the main contribution of stronger interactions is an increase in the delay time, while the underlying physics is not modified, as confirmed by the data collapse in figure 1. We compare the numerical results with our theoretical prediction
\begin{equation} T_d = a_0 \left ( 1-\frac{\phi_0}{\phi_m} \right ) \tau_d, \end{equation}
where 
$a_0 \approx 5/2$ is a fit parameter likely resulting from the way in which we define the gel interface in our system. We obtained this by computing the time required for an interface moving at a constant velocity 
$V_g(1-\phi _0 / \phi _m)^{-1}$ to cover a distance 
$\zeta$. This result matches our numerical data extremely well across the range of 
$U$ and 
$\phi _0$ considered. Turning to our simple description of the early-stage settling, we conclude that the bulk of the gel must settle at a distance at least 
$\zeta$ to free itself from the non-Newtonian responses of the boundary. Once it does, it is able to settle at a constant velocity 
$V_g$. We will provide further elaboration on this point when we provide additional context to our research in the next section.
Delay time for gravitational gel collapse 
$T_d^{\ast }$ (scaled by 
$c(U)=\sqrt {U(1-{\rm e}^{-U})}$), in units of 
$t_s$, and as a function of initial volume fraction 
$\phi _0$. Numerical results for different values of the potential strength 
$U$ are indicated using symbols; see legend. The dashed line represents the theoretical prediction (3.6), scaled in the same way as the data points.
4. Discussion
We have put forward a theoretical description of a colloidal gel, where we have moved away from previous modelling (Buscall & White Reference Buscall and White1987; Evans & Starrs Reference Evans and Starrs2002; Manley et al. Reference Manley, Skotheim, Mahadevan and Weitz2005; Buzzaccaro et al. Reference Buzzaccaro, Secchi, Brambilla, Piazza and Cipelletti2012; Darras et al. Reference Darras, Dasanna, John, Gompper, Kaestner, Fedosov and Wagner2022) by introducing a dilatational viscosity for the solid phase. This imparts resistance to compression and expansion, which strongly affects the gel's settling behaviour under gravity. In this section, we discuss our results in the context of experimental studies (Derec et al. Reference Derec, Senis, Talini and Allain2003; Manley et al. Reference Manley, Skotheim, Mahadevan and Weitz2005; Harich et al. Reference Harich, Blythe, Hermes, Zaccarelli, Sederman, Gladden and Poon2016; Zhou Reference Zhou2018; Darras et al. Reference Darras, Dasanna, John, Gompper, Kaestner, Fedosov and Wagner2022; John et al. Reference John, Kaestner, Wagner and Darras2024) and other computational (Padmanabhan & Zia Reference Padmanabhan and Zia2018; Varga, Hofmann & Swan Reference Varga, Hofmann and Swan2018; de Graaf et al. Reference de Graaf, Torre, Poon and Hermes2023), and theoretical (Evans & Starrs Reference Evans and Starrs2002; Derec et al. Reference Derec, Senis, Talini and Allain2003; Manley et al. Reference Manley, Skotheim, Mahadevan and Weitz2005; Darras et al. Reference Darras, Dasanna, John, Gompper, Kaestner, Fedosov and Wagner2022) work.
4.1. Layering and length scales
Our 
$t = 0$ analytic expression, (3.5), provides insight into our numerical observations, which we consider valid on the basis of the correspondence in figure 3. A homogeneous initial sample can be approximately divided into three regions: a top and bottom layer, both with a thickness of 
$\zeta$ – the gel stress screening length – and a bulk containing the rest of the sample. As this homogeneous gel begins to sediment, the top and bottom layers resist stretching and compression, respectively, while the bulk settles with a uniform velocity 
$V_g$.
In our calculations, the bottom part typically lacks the strength to withstand compression, leading to the immediate creation of a bottom dense layer that grows without delay (see figure 2a). This indicates that the gel stress screening length 
$\zeta$ is typically too small to propagate the stress response from the bottom of the sample to a height sufficient for forming a resistant layer thick enough to prevent the gel bulk from collapsing. A bottom layer that can be clearly distinguished from the bulk gel for an extended period of time is also observed in certain experimental configurations (Harich et al. Reference Harich, Blythe, Hermes, Zaccarelli, Sederman, Gladden and Poon2016). These are usually associated with gels having strong interparticle interactions, suggesting that their size depends on the strength of the attractions, which aligns with our findings for 
$\zeta$.
Notably, a similar characteristic stress screening length is also reported in the work of Evans & Starrs (Reference Evans and Starrs2002), where a comparable theoretical framework is applied. Their study demonstrated that the characteristic velocity of the collapse could be influenced by the size of this screening length, particularly when the sample size is comparable to it. This suggests that gels with extremely long delay times (Harich et al. Reference Harich, Blythe, Hermes, Zaccarelli, Sederman, Gladden and Poon2016; Zhou Reference Zhou2018; Darras et al. Reference Darras, Dasanna, John, Gompper, Kaestner, Fedosov and Wagner2022) may fall into this category. In such cases, the velocity of the colloids is exponentially suppressed throughout the entire sample, resulting in minimal or no dynamics, and consequently, a significant delay in the onset of collapse.
4.2. Settling and rupture
The initial settling causes an upward solvent flow (backflow), with the velocity dictated primarily by the porous structure of the networked colloids. We believe this to be realistic for unruptured gels in experimental settings (Darras et al. Reference Darras, Dasanna, John, Gompper, Kaestner, Fedosov and Wagner2022), as well as for colloids with small sizes (Derec et al. Reference Derec, Senis, Talini and Allain2003; Manley et al. Reference Manley, Skotheim, Mahadevan and Weitz2005). In these systems, the interface velocity follows a power law that is different from the one in our model, due to the specific choice for porosity 
$K$. The scaling predicted by our theory follows from the assumption that pores are locally homogeneously distributed within an REV, thus recovering the expression for a single isolated sphere sedimenting in an unbounded Newtonian solvent. A more accurate porosity for fractal-like aggregates should also depend on the fractal dimensions of these aggregates (Allain et al. Reference Allain, Cloitre and Wafra1995; Derec et al. Reference Derec, Senis, Talini and Allain2003; Manley et al. Reference Manley, Skotheim, Mahadevan and Weitz2005). However, estimating such a quantity is far from trivial, as it generally depends on the volume fraction, the strength of the interactions, and the shape of the particles (Derec et al. Reference Derec, Senis, Talini and Allain2003), making it space-dependent in our description. Therefore, we choose to use a simpler and well-defined expression as defined in (2.17), which effectively covers the entire range of volume fractions and connects different regimes with a unified qualitative description.
When the centre of mass of the bulk displaces by a height 
$\zeta$, the top part of the gel is maximally extended. This triggers a detachment of the bulk gel, after which the remaining top part effectively becomes debris. This debris subsequently deposits on the gel interface at the top of the bulk area as it continues to settle. This accumulation of debris at the top of the gel interface has been observed in experiments (Harich et al. Reference Harich, Blythe, Hermes, Zaccarelli, Sederman, Gladden and Poon2016; Zhou Reference Zhou2018). In most experiments, there is a curved meniscus to the sample, and the mechanisms by which the debris forms are presumably not identical. However, the inverted-cuvette set-up of Zhou (Reference Zhou2018) comes closest to our 1-D theoretical analysis. Here, a denser layer on top of a uniformly settling gel was observed, though the formation of a thinner region before rupture could not be discerned. It may be that 
$\zeta$ is very different in these experiments, or that the boundary conditions that we considered are not commensurate with those of Zhou (Reference Zhou2018).
Once the interface has detached from the top of the sample, it continues to sediment with constant velocity until it reaches the dense bottom layer. The interface subsequently slows down exponentially as the sample compacts and 
$h$ approaches its asymptotic limit 
$h^{\infty }$. The presence of such a compaction regime is also consistent with experimental observations and previous modelling (Derec et al. Reference Derec, Senis, Talini and Allain2003; Manley et al. Reference Manley, Skotheim, Mahadevan and Weitz2005; Harich et al. Reference Harich, Blythe, Hermes, Zaccarelli, Sederman, Gladden and Poon2016; Razali et al. Reference Razali, Fullerton, Turci, Hallett, Jack and Royall2017; Darras et al. Reference Darras, Dasanna, John, Gompper, Kaestner, Fedosov and Wagner2022; de Graaf et al. Reference de Graaf, Torre, Poon and Hermes2023). However, we do not believe that the constant velocity regime is representative of the ruptured regime that is encountered in experiments (Derec et al. Reference Derec, Senis, Talini and Allain2003; Harich et al. Reference Harich, Blythe, Hermes, Zaccarelli, Sederman, Gladden and Poon2016; Zhou Reference Zhou2018). That is, the transition from delayed to linear collapse is widely accepted to occur when cracks (Derec et al. Reference Derec, Senis, Talini and Allain2003) or streamers (Harich et al. Reference Harich, Blythe, Hermes, Zaccarelli, Sederman, Gladden and Poon2016; Zhou Reference Zhou2018) emerge in a previously homogeneous bulk gel. They could either give preferential pathways for fluid flow or yield the gel network in its entirety (Harich et al. Reference Harich, Blythe, Hermes, Zaccarelli, Sederman, Gladden and Poon2016), due to erosion by fluid backflow (Varga et al. Reference Varga, Hofmann and Swan2018; Zhou Reference Zhou2018; John et al. Reference John, Kaestner, Wagner and Darras2024). Our model gel remains homogeneous, and its bulk stress is relatively unaffected before compaction occurs.
4.3. Extending the delay time
As we do not have cracks or rupture, we might expect our delay times to be longer than in experiment, as also suggested by the inverted-cuvette experiment in Zhou (Reference Zhou2018). However, 
$T_{d}$ is consistently shorter than reported in real systems (Harich et al. Reference Harich, Blythe, Hermes, Zaccarelli, Sederman, Gladden and Poon2016; Zhou Reference Zhou2018; Darras et al. Reference Darras, Dasanna, John, Gompper, Kaestner, Fedosov and Wagner2022). Crucial to understanding the origin of this difference is the form of our dilatational viscosity 
$\lambda (\phi,U)$. Scaling this viscosity by the function 
$c(U)$, as shown in figure 3, preserves its dependency on 
$\phi _0$. This is evidenced by the data collapse in this figure. If we had a different form of 
$c(U)$, then we would expect the collapsed curve to still hold, but the real delay time to change significantly. In proposing a form of 
$\lambda$, we ignored the influence of attractive interactions on the system's relaxation time. These should strongly suppress colloidal diffusion, thereby increasing drastically 
$c(U)$ and 
$T_{d}$. Incorporating a more detailed expression for 
$\lambda$ capable of capturing these effects should yield delay times more closely aligned with experimental values, and could be considered in follow-up studies.
Finally, we note that the computational study by Padmanabhan & Zia (Reference Padmanabhan and Zia2018) predicted both a delay and a collapse without accounting for hydrodynamic interactions. They attribute their observed delay to gravity-enhanced coarsening driven by negative osmotic compressibility. This mechanism is not captured directly in our modelling, which does partially consider colloids diffusion and local rearrangement only in constructing the dilatational viscosity of (2.21). However, their predicted delay times are significantly shorter than those observed in experiments, and the range of gravitational Péclet numbers 
${{Pe}}_{g}$ that they investigated is limited to small values, whereas not all experimental systems fall within this range (Harich et al. Reference Harich, Blythe, Hermes, Zaccarelli, Sederman, Gladden and Poon2016; Zhou Reference Zhou2018). In contrast, our model assumes that in the absence of gravity, there would be no osmotic contribution to the stress affecting local density. We consider the dynamics to be dominated by gravity and the viscous response of the gel network, mediated by hydrodynamics, making the Péclet number effectively infinitely large, as reflected by the missing diffusive term in (2.3).
4.4. Going beyond one dimension
The experiments carried out in the Poon lab (Meeker Reference Meeker1998; Starrs et al. Reference Starrs, Poon, Hibberd and Robins2002; Harich et al. Reference Harich, Blythe, Hermes, Zaccarelli, Sederman, Gladden and Poon2016; Zhou Reference Zhou2018) suggest that a 1-D description may lack necessary features to accurately replicate real gel behaviour. Missing features include the curvature of the suspension–air meniscus, which represents the top boundary of the system, and the interaction between the gel and the meniscus, as well as between the gel and the other solid-confining boundaries (Senis et al. Reference Senis, Gorre-Talini and Allain2001; Evans & Starrs Reference Evans and Starrs2002). That is, tangential slip of colloids along a curved meniscus could lead to local weakening of the structure at the top of the sample, and debris accumulation at specific points, which subsequently yields the gel and triggers the onset of rapid collapse (Zhou Reference Zhou2018).
Additionally, the influence of shear stress response in the gel does not play a role in our 1-D calculations. This response is expected to be relevant for the same reasons that we mentioned, when we introduced the dilatational stress. We expect that shear stress response could potentially produce more resistant gels, as any variation in the flow of colloids in the plane normal to its direction would generate a response. These variations can be caused by microscale factors such as density fluctuations coming from trapped solvent droplets. They can also stem from the applied boundary conditions, e.g. the presence of a liquid–air meniscus (Harich et al. Reference Harich, Blythe, Hermes, Zaccarelli, Sederman, Gladden and Poon2016; Zhou Reference Zhou2018) or an overall incline of the cuvette (Senis et al. Reference Senis, Gorre-Talini and Allain2001).
Finally, our model assumes an instantaneous stress response for the gel network. However, in reality, the gel network undergoes erosion, which changes the morphology of the network and weakens its stress response over time (Varga et al. Reference Varga, Hofmann and Swan2018; de Graaf et al. Reference de Graaf, Torre, Poon and Hermes2023). This erosion also alters the pore distribution, affecting the porosity 
$K$. That is, a gel that is initially much stronger than what we considered in this work will weaken over time, leading to a more pronounced and violent transition from delay to linear collapse. Future work will focus on the implementation of these ‘weakening’ features.
5. Conclusions and outlook
Summarizing, we have developed a theoretical model by which we can study the response of colloidal gels to gravitational stress. Our approach treats the gel as a viscous medium, incorporating a dilatational viscosity for the colloidal phase that varies with local density. Herein, we depart from previous modelling (Buscall & White Reference Buscall and White1987; Evans & Starrs Reference Evans and Starrs2002; Manley et al. Reference Manley, Skotheim, Mahadevan and Weitz2005; Buzzaccaro et al. Reference Buzzaccaro, Secchi, Brambilla, Piazza and Cipelletti2012; Darras et al. Reference Darras, Dasanna, John, Gompper, Kaestner, Fedosov and Wagner2022), which allows us to qualitatively capture several experimentally observed features of colloidal gel sedimentation. The gel viscosity is complemented by Darcy's law, by which we represent flow as a local drag force acting on the solid phase. That is, we work at a level that coarse grains out the complex microscale hydrodynamic flows and interactions.
We computed the full dynamics of the system numerically in an effective 1-D geometry, i.e. only the sample height is relevant. This reveals that our model predicts three distinct regimes: a delay, followed by linear settling, and a final exponential compacting. All of these elements have been reported in experimental systems (Harich et al. Reference Harich, Blythe, Hermes, Zaccarelli, Sederman, Gladden and Poon2016; Zhou Reference Zhou2018; Darras et al. Reference Darras, Dasanna, John, Gompper, Kaestner, Fedosov and Wagner2022). Surprising, solving analytically for the initial (
$t = 0$) behaviour is sufficient to understand the origin of the salient features of the full collapse dynamics in our model. This analysis reveals that there is a natural length scale 
$\zeta$ in the system, representing a gel stress screening length. Thus a sample can be divided into three regions when 
$\zeta < H$: top and bottom layers, and a middle ‘bulk’ part of the gel. The latter remains relatively unaffected during the initial part of the settling, again mirroring experiments (Harich et al. Reference Harich, Blythe, Hermes, Zaccarelli, Sederman, Gladden and Poon2016; Zhou Reference Zhou2018). However, when the former is fully stretched, it triggers the onset of the linear settling regime and the subsequent formation of a dense layer of debris on top of the settling gel. This is reminiscent of some experimental observations that were made on depletion-based colloidal gels (Zhou Reference Zhou2018).
However, the present modelling does not capture all settling behaviours observed in experiment. In part, this can be attributed to the approximations that we made in determining the resistance of the gel to flow. This makes the delay time much shorter than is typically observed in experiment. However, we provide a means to amend this discrepancy through suitable modification of the dissipative term. In part, the choice of modelling the dynamics in 1-D proved useful to gain basic insight, but it fails to capture some of the more essential features of the experiments. For example, we do not account for the the curvature of a liquid–air meniscus, which appears to control much of the onset of rapid collapse (Harich et al. Reference Harich, Blythe, Hermes, Zaccarelli, Sederman, Gladden and Poon2016; Zhou Reference Zhou2018).
Addressing these gaps in our present modelling could bring our results in closer agreement with the experiments. Future work will therefore extend the analysis to higher dimensions and incorporate erosion mechanisms. The present study provides a solid foundation for this and may be readily built upon in other directions.
Acknowledgements
We are grateful to Dr A. Darras, M. van Leeuwen and M. van Schaik for useful discussion.
Funding
The authors acknowledge NWO for funding through OCENW.KLEIN.354.
Declaration of interests
The authors report no conflict of interest.
Data availability statement
The data that support the findings of this study are openly available at https://doi.org/10.24416/UU01-F3SLMJ.
Author contributions
Conceptualization by J.d.G. and K.W.T.; methodology by K.W.T.; numerical calculations by K.W.T.; validation by K.W.T.; investigation by K.W.T.; writing original draft by K.W.T.; writing review and editing by J.d.G.; funding acquisition by J.d.G.; resources by J.d.G.; supervision by J.d.G.
Appendix A. Numerical scheme
In this appendix, we provide the details of the numerical scheme used to solve our model in one dimension. For non-homogeneous density profiles, by combining (2.4) with (2.16), we obtain a non-homogeneous linear differential equation for the gel stress 
$\sigma _c(z)$,
\begin{align} \sigma_c &= \lambda(\phi) \, \partial_z \left [ (1-\phi)\,K(\phi) \left ({-}v_g \phi + \frac{\sigma^2}{\mu}\,\partial_z \sigma_c \right ) \right ] \end{align}
\begin{align} &= A(\phi) + B(\phi)\,\partial_z \sigma_c + C(\phi)\,\partial^2_z \sigma_c , \end{align}where the non-constant coefficients are defined as
$$\begin{gather} A(\phi) ={-} \partial_z \phi \, \lambda(\phi) \, \partial_{\phi} \left [ \phi (1-\phi)\,K(\phi) \right ] v_g , \end{gather}$$
$$\begin{gather}B(\phi) = \partial_z \phi \, \lambda(\phi) \, \partial_{\phi} \left [ (1-\phi)\,K(\phi) \right ] \mu^{{-}1} \sigma^2 , \end{gather}$$
$$\begin{gather}C(\phi) = \lambda(\phi)\,(1-\phi)\,K(\phi)\,\mu^{{-}1} \sigma^2, \end{gather}$$and the boundary conditions are given by
\begin{equation} \left. \partial_z \sigma_c \right |_{z=0,H} = \Delta \rho\,g \, \phi|_{z=0,H}. \end{equation}
The system height 
$H$ is subdivided in 
$N+1$ intervals of size 
$\delta z$, and all fields are represented on this grid using a superscript indicating their evaluation positions. We employ central second-order and first-order finite-difference schemes in the bulk and at the boundaries, respectively, to derive the following finite-difference equations for 
$\sigma _c$:
$$\begin{gather} \left . \sigma^i_c \right |_{0< i< N} = \frac{A^i+ \dfrac{ B^i}{2\,\delta z }\,( \sigma_c^{i+1} - \sigma_c^{i-1} )+ \dfrac{C^i}{\delta z^2}\,(\sigma_c^{i+1} + \sigma_c^{i-1})} {1 + \dfrac{2C^i}{ \delta z^2}} , \end{gather}$$
$$\begin{gather}\left . \sigma^i_c \right |_{i=0} ={-}\Delta \rho\,g\phi^0 + \sigma^1_c , \end{gather}$$
$$\begin{gather}\left . \sigma^{i}_c \right |_{i=N} = \Delta \rho\,g\phi^N - \sigma^{N-1}_c . \end{gather}$$
These equations are solved iteratively to update the stress values until convergence is reached, with the tolerance set to 
$0.001\,\%$.
Once the numerical solution for the stress is computed, it is divided by the dilatational viscosity 
$\lambda$, which is computed using the volume fraction of the current time step, to obtain 
$\partial _z v_c(z,t)$. Next, we integrate this result to find the colloid velocity field,
\begin{equation} \int_0^z {\rm d}z' \, \partial_{z'} v_c(z',t) = v_c(z,t), \end{equation}
where the boundary velocity term 
$v_c(0,t)$ is simplified using mass conservation and the fact that the system is isolated. An approximation of the above integral is computed using the trapezoidal rule, resulting in the discretized colloid velocity:
$$\begin{gather} \left. v_c^i \right |_{0< i< N} = \frac{\delta z}{2} \sum_{j=1}^i \left ( \frac{\sigma_c^{j+1}}{\lambda^{j+1}} + \frac{\sigma_c^{j-1}}{\lambda^{j-1}} \right ) , \end{gather}$$
$$\begin{gather}\left. v_c^i \right |_{i=0} = \left. v_c^i \right |_{i=N} = 0. \end{gather}$$ Now that the colloid velocity field is known, we can advance the colloid volume fraction one time step forwards using (2.3). This is solved numerically using an upwind finite-volume method Versteeg & Malalasekera (Reference Versteeg and Malalasekera1995). Standard finite-difference schemes would fail to converge due to the strong convective nature of (2.3), given the absence of explicit diffusive terms. Thus the change in volume fraction is computed approximately at the centre of each of the 
$N+1$ space intervals using the biased fluxes
\begin{align} J_{{up}}^k &= \begin{cases} \phi^{k} v_c^{i+1}, & v_c^{i+1}>0, \\ \phi^{k+1} v_c^{i+1}, & v_c^{i+1}<0, \end{cases} \end{align}
\begin{align} J_{{down}}^k &= \begin{cases} \phi^{k-1} v_c^{i}, & v_c^i>0, \\ \phi^{k} v_c^{i}, & v_c^i<0, \end{cases} \end{align}
where we have introduced the indices 
$k = i + 1/2$ to relate the velocity field, computed at the 
$N+1$ nodes, with the volume fraction field, computed at the 
$N$ centres of the intervals formed by consecutive nodes. The instantaneous change in volume fraction at a given time step can then be computed using
\begin{equation} \partial_t \phi^{k}(t) ={-} J_{{up}}^k + J_{{down}}^k, \end{equation}
and the updated values of 
$\phi$ are obtained using a first-order forward Euler scheme:
\begin{equation} \phi^{k}(t+\delta t) = \delta t \, \partial_t \phi^{k}(t). \end{equation}
Open access
