2.1 Growth of a critical streamer

To arrive at the radial growth rate of the streamer or channel, consider the transport process that move particles into the channel and the reverse process of attachment into the network. The evolution of the number of particles $N$ in the cylindrical channel of height $L$ and radius $R(t)$ is given by the net flux:

(2.1) $$\begin{eqnarray}\displaystyle \frac{\text{d}N}{\text{d}t}=J_{in}-J_{out}=2\unicode[STIX]{x03C0}RL(j_{in}-j_{out}), & & \displaystyle\end{eqnarray}$$

where $J$ and $j$ are the rate and flux per unit area of particles into the growing channel ( $in$ ) and out into the network ( $out$ ), respectively. To the best approximation, for tall enough channels ( $L$ is large compared to the mesh size of the gel network) the density of particles in the interior of the channel will be equal to its bulk value $\unicode[STIX]{x1D719}$ , so that the number of particles is related to the geometry of the channel through $N=3\unicode[STIX]{x03C0}R^{2}L\unicode[STIX]{x1D719}/4\unicode[STIX]{x03C0}a^{3}$ . When the channel spans the height of the network, $L$ is unchanged over time and the one-dimensional growth model of $R(t)$ is:

(2.2) $$\begin{eqnarray}\displaystyle \frac{\text{d}R}{\text{d}t}=\frac{4\unicode[STIX]{x03C0}a^{3}}{3\unicode[STIX]{x1D719}}(j_{in}-j_{out})=\frac{1}{n}(j_{in}-j_{out}), & & \displaystyle\end{eqnarray}$$

where $n$ is the bulk number density of particles.

The flux of particles per unit area into the open channel is proportional the surface number density of particles at the network–channel interface, ${\tilde{n}}$ , and the rate of particle bond breaking, $k_{break}$ :

(2.3) $$\begin{eqnarray}\displaystyle j_{in}=d_{1}k_{break}{\tilde{n}}. & & \displaystyle\end{eqnarray}$$

Note that throughout this paper, $d_{i}$ for $i=\{1,\ldots ,6\}$ are unknown $O(1)$ non-dimensional, scalar constants. The boundary between the interior of the channel at bulk density $\unicode[STIX]{x1D719}$ and the percolated gel network of fractal dimension $d_{f}$ constitutes the ‘wall’ of the streamer. For a fractal structure, the surface density of particles attached to the network at this interface is ${\tilde{n}}=\unicode[STIX]{x1D719}^{d_{f}/3}(d_{2}/a^{2})(R/a)^{(d_{f}/3-1)}$ . As one would expect, for a compact structure with $d_{f}=3$ , ${\tilde{n}}$ is independent of the radius. In the absence of an external force, we hypothesize that $k_{break}$ is set by the Kramers escape time of a particle diffusing out of an attractive well of depth $U$ and width $\unicode[STIX]{x1D6E5}$ (Kramers Reference Kramers1940):

(2.4) $$\begin{eqnarray}\displaystyle k_{break}=\unicode[STIX]{x1D70F}_{K}^{-1}=\frac{D}{a^{2}}\left(\frac{a}{\unicode[STIX]{x1D6E5}}\right)^{2}\frac{U}{kT}\text{e}^{-U/(d_{3}kT)}=\frac{D}{a^{2}}\unicode[STIX]{x1D6FF}^{-2}\unicode[STIX]{x1D716}^{-1}\text{e}^{-1/(d_{3}\unicode[STIX]{x1D716})}, & & \displaystyle\end{eqnarray}$$

where $D$ is the diffusivity of a particle of radius $a$ . The coefficient $d_{3}$ is included to account for the fact that the particles sit in a complex energy landscape in the gel network for which the barrier between bound and unbound states could be smaller than $U$ owing to elastic stresses stored intrinsically in the network during its formation, or could be larger due to a high coordination number. Due to the back flow of fluid through the channel, the particles at the interface will feel a hydrodynamic drag, due to the shear stress, $\unicode[STIX]{x1D70F}$ , at the channel wall. This results in a relative force $d_{4}\unicode[STIX]{x1D70F}a^{2}/kT$ that stretches the inter-particle bonds and accelerates the breaking process. This activated hopping leads to a higher rate of erosion. Assuming simple Poiseuille flow in the cylindrical channel the shear stress at the wall is related to the channel radius through $\unicode[STIX]{x1D70F}=|\unicode[STIX]{x1D735}p|R/2$ and $k_{break}$ becomes:

(2.5) $$\begin{eqnarray}\displaystyle k_{break}=\frac{D}{a^{2}}\unicode[STIX]{x1D6FF}^{-2}\unicode[STIX]{x1D716}^{-1}\text{e}^{-1/d_{3}\unicode[STIX]{x1D716}+d_{4}|\unicode[STIX]{x1D735}p|Ra^{2}\unicode[STIX]{x1D6E5}/(2kT)}=\frac{D}{a^{2}}\unicode[STIX]{x1D6FF}^{-2}\unicode[STIX]{x1D716}^{-1}\text{e}^{-1/(d_{3}\unicode[STIX]{x1D716})+R/R^{\ast }}, & & \displaystyle\end{eqnarray}$$

$R^{\ast }$ is an effective gravitational length, i.e. the characteristic length scale of the erosion process:

(2.6) $$\begin{eqnarray}\displaystyle R^{\ast }=\frac{2kT}{d_{4}a^{2}\unicode[STIX]{x1D6E5}|\unicode[STIX]{x1D735}p|}=\frac{8\unicode[STIX]{x03C0}a}{3d_{4}}\unicode[STIX]{x1D716}\unicode[STIX]{x1D6FF}^{-1}G^{-1}. & & \displaystyle\end{eqnarray}$$

It sets the characteristic pore size, beyond which the activated hopping dominates the bond breaking. The erosive flux of particles into the channel becomes:

(2.7) $$\begin{eqnarray}\displaystyle j_{in}=\frac{4\unicode[STIX]{x03C0}}{3}d_{1}d_{2}\unicode[STIX]{x1D719}^{d_{f}/3-1}R\left(\frac{R}{a}\right)^{(d_{f}/3-2)}\frac{nD}{a^{2}}\unicode[STIX]{x1D6FF}^{-2}\unicode[STIX]{x1D716}^{-1}\text{e}^{-1/(d_{3}\unicode[STIX]{x1D716})+R/R^{\ast }}. & & \displaystyle\end{eqnarray}$$

The flux of individual particles currently in the bulk of the channel back onto the network is driven by diffusion according to:

(2.8) $$\begin{eqnarray}\displaystyle j_{out}=\frac{nD}{x}=\frac{nD}{R}f(Pe), & & \displaystyle\end{eqnarray}$$

where $Pe$ is the advective Péclet number near the wall $Pe=\unicode[STIX]{x1D70F}R^{2}/\unicode[STIX]{x1D702}D=3\unicode[STIX]{x03C0}R^{3}a|\unicode[STIX]{x1D735}p|/kT=(9/4)\unicode[STIX]{x1D716}G(R/a)^{3}$ . Here, $x$ is the thickness of the diffusive boundary layer at the wall of the channel: $f(Pe)=1$ for $Pe\ll 1$ and $f(Pe)\sim Pe^{1/2}$ when $Pe\gg 1$ (Acrivos & Goddard Reference Acrivos and Goddard1965; Goddard & Acrivos Reference Goddard and Acrivos1966). Except for very early times when the channel initially forms, the radius of the channel exceeds the primary particle size, $R\gg a$ , and $Pe\gg 1$ so that to first order:

(2.9) $$\begin{eqnarray}\displaystyle j_{out}=d_{5}R\frac{nD}{a^{2}}\unicode[STIX]{x1D716}^{1/2}G^{1/2}\left(\frac{R}{a}\right)^{-1/2}. & & \displaystyle\end{eqnarray}$$

Note that $j_{in}\gg j_{out}$ for $R\gg a$ so that to first approximation the flux of particles from the channel onto the network can be neglected, i.e. $j_{out}\approx 0$ . Rescaling time $t$ on the pure diffusive Kramers escape time, $\hat{t}=tD/(\unicode[STIX]{x1D6E5}^{2}\unicode[STIX]{x1D716})\text{e}^{-1/(d_{3}\unicode[STIX]{x1D716})}=t/\unicode[STIX]{x1D70F}_{K}$ , the model predicts the following evolution equation for $R(t)$ :

(2.10) $$\begin{eqnarray}\displaystyle \frac{\text{d}R}{\text{d}\hat{t}}=\frac{4\unicode[STIX]{x03C0}}{3}d_{1}d_{2}a\unicode[STIX]{x1D719}^{d_{f}/3-1}\left(\frac{R}{a}\right)^{d_{f}/3-1}\text{e}^{R/R^{\ast }}. & & \displaystyle\end{eqnarray}$$

There are two competing effects that set the erosion of particles and the evolution of the radius of the channel. For fractal structures with $d_{f}<3$ , the quantity $(R/a)^{d_{f}/3-1}$ suggests that the growth rate of $R$ is slowed as $R$ becomes larger as the total number of particles eroded by the back flow is decreasing relative to the number of particles in the channel. In contrast, the term $\text{e}^{R/R^{\ast }}$ , arises because the rate of bond breaking is exponentially dependent on the shear stress exerted by back flow, which itself is linear in $R$ . Hence the activated rate of erosion grows exponentially with the channel radius. The result is that for $R>R^{\ast }$ , the erosion process accelerates exponentially. More generally, such first-order ordinary differential equations, which have a super-linear flux or grow rate, exhibit a finite-time blowup singularity (Ide & Sornette Reference Ide and Sornette2002).

In the case of this erosion model, we predict an exponentially dependent flux, an ultrafast type of super-linear growth, so that the channel radius will grow infinitely large at a critical point in time. This blowup time is defined such that: $R\rightarrow \infty$ as $t\rightarrow t_{blowup}$ . This is the main result of our phenomenological model. At a certain critical point in time in freely settling gels a hydrodynamic instability occurs. The channel radius grows unstably to span the width of the gel. In practice, this corresponds to the streamer being of comparable size to the width of the container, as seen by Starrs et al. (Reference Starrs, Poon, Hibberd and Robins2002), when large scale fluid back flow deconstructs the gel network. The solid is ripped apart and the gel rapidly collapses. Note, that due to the exponential growth profile, the fate of the network is determined long before $R$ reaches the system size. When $R\geqslant R^{\ast }$ , the channel growth rate is exponential and catastrophic collapse is unavoidable. However, macroscopically the radial growth of the channel may not manifest itself in a dramatic increase in the settling velocity of the network’s top interface until $t$ is right near $t_{blowup}$ .

To explain this observation consider that for an incompressible Newtonian suspending fluid the settling velocity of a gel and therefore the velocity of the interface will be proportional to the average fluid back flow velocity, $\langle u_{f}\rangle$ , driven by the constant pressure gradient $|\unicode[STIX]{x1D735}p|$ across the network. Here, $\langle u_{f}\rangle$ is the average volumetric flow rate of fluid back flow through the intact gel and the streamer channel. For simplicity, the intact gel network can be thought of as a porous Darcy medium with permeability $\unicode[STIX]{x1D705}$ , whereas the streamer is a cylindrical channel with Poiseuille flow, as before.

The fluid velocity through the gel network is given by:

(2.11) $$\begin{eqnarray}\displaystyle u_{f}^{Darcy}=-\frac{\unicode[STIX]{x1D705}}{\unicode[STIX]{x1D702}}\unicode[STIX]{x1D735}p. & & \displaystyle\end{eqnarray}$$

The total flow rate through the medium with cross-sectional dimension $L$ and a cylindrical pore of radius $R$ within it is $Q_{f}^{Darcy}=-(L^{2}-\unicode[STIX]{x03C0}R^{2})(\unicode[STIX]{x1D705}/\unicode[STIX]{x1D702})\unicode[STIX]{x1D735}p$ . In the cylindrical channel the velocity profile is: $u_{f}^{channel}=-(1/4\unicode[STIX]{x1D702})\unicode[STIX]{x1D735}p(C-r^{2})$ , where $r$ is the radial position and the constant $C$ is set by the boundary condition at the network–channel interface, $r=R(t)$ . This boundary condition will be highly dependent on the exact fractal structure of the network. At lowest order this can be captured by introducing an inverse Navier’s slip length $\unicode[STIX]{x1D706}$ so that at the boundary (Beavers & Joseph Reference Beavers and Joseph1967): $u_{f}^{channel}(r=R(t))-u_{f}^{Darcy}=-\unicode[STIX]{x1D706}(\text{d}u_{f}^{P}/\text{d}r)$ . Therefore the resulting parabolic flow profile inside the channel is:

(2.12) $$\begin{eqnarray}\displaystyle u_{f}^{channel}=-\frac{1}{4\unicode[STIX]{x1D702}}\unicode[STIX]{x1D735}p(R^{2}-r^{2}+2\unicode[STIX]{x1D706}R+4\unicode[STIX]{x1D705}), & & \displaystyle\end{eqnarray}$$

with a volumetric flow rate: $Q_{f}^{channel}=(-\unicode[STIX]{x03C0}R^{2}/8\unicode[STIX]{x1D702})\unicode[STIX]{x1D735}p(R^{2}+4\unicode[STIX]{x1D706}R+8\unicode[STIX]{x1D705})$ . Thus the average fluid back flow velocity over the pore and porous medium as the gel settles is given by:

(2.13) $$\begin{eqnarray}\displaystyle \langle u_{f}\rangle =\frac{Q_{f}^{channel}+Q_{f}^{Darcy}}{L^{2}}=-\frac{1}{\unicode[STIX]{x1D702}}\unicode[STIX]{x1D735}p\left[\unicode[STIX]{x1D705}+\frac{\unicode[STIX]{x03C0}}{8}\left(\frac{R}{L}\right)^{2}(R^{2}+4\unicode[STIX]{x1D706}R)\right]. & & \displaystyle\end{eqnarray}$$

The channel size will be negligible compared to the overall size of the porous medium for the majority of the growth process, i.e. $R/L\ll 1$ . The rate of collapse of the gel is therefore limited by the back flow of the fluid through the solid network. It is controlled by $\unicode[STIX]{x1D705}$ , yielding a uniform settling profile of the top interface as if there were no growing channel (Manley et al. Reference Manley, Skotheim, Mahadevan and Weitz2005). However, as the settling proceeds and $t\rightarrow t_{blowup}$ , $R$ will grow rapidly and the channel cross-section will become significant, $R/L\sim O(1)$ . The second term in (2.13) will dominate the fluid back flow velocity. Consequently, the interface velocity will appear to increase rapidly as the size of the channel becomes comparable to the size of the gel network. Therefore $t_{blowup}$ sets the time scale for the onset of macroscopic collapse.

2.2 Finite-time singularity in channel growth

A closed form expression for the finite-time singularity, or blowup time predicted by the model can be found by specifying the initial value for the radius of the channel, $R(\hat{t}=0)$ :

(2.14) $$\begin{eqnarray}\displaystyle \hat{t}_{blowup}=d_{6}R^{\ast ^{2-d_{f}/3}}\unicode[STIX]{x1D6E4}(2-d_{f}/3,R(0)/R^{\ast })\unicode[STIX]{x1D719}^{1-d_{f}/3}, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D6E4}(a,x)=\int _{x}^{\infty }s^{a-1}\text{e}^{-s}\,\text{d}s$ is the incomplete gamma function and $d_{6}=(4\unicode[STIX]{x03C0}d_{1}d_{2}/3)^{-1}$ . Our model predicts that the finite-time singularity is intrinsic to all freely settling colloidal gel networks. One interesting feature is that for an initial channel radius of zero, $R(0)=0$ , the blowup time is non-zero. In fact, all initial channel radii $R(0)\ll R^{\ast }$ have nearly the same blowup time: $\hat{t}_{blowup}=d_{6}R^{\ast ^{2-d_{f}/3}}\unicode[STIX]{x1D6E4}(2-d_{f}/3,0)\unicode[STIX]{x1D719}^{1-d_{f}/3}\approx 0.9d_{6}R^{\ast ^{2-d_{f}/3}}\unicode[STIX]{x1D719}^{1-d_{f}/3}$ , when $d_{f}=2$ . In essence, the gel is filled with pores having a heterogeneous size distribution. The fluid will select the largest initial pore and erosion of that pore will be favoured over others.

As discussed in § 4.2 heterogeneities within the initial gel may form from other processes beyond mere thermally driven restructuring (Lu et al. Reference Lu, Zaccarelli, Ciulla, Schofield, Sciortino and Weitz2008). These can have a number of origins, including falling debris that accumulates at an interface (Harich et al. Reference Harich, Blythe, Hermes, Zaccarelli, Sederman, Gladden and Poon2016), external forcing fields (Teece et al. Reference Teece, Hart, Hsu, Gilligan, Faers and Bartlett2014) or included bubbles. Once a streamer is born with initial size $R(0)$ , the growth rate is universally described by (2.10) and the final stages of breakup leading to collapse are identically captured by the model with $R^{\ast }$ being the only controlling parameter.

Because the incomplete gamma function takes on values $[0,1]$ , the behaviour of the blowup time is dictated by the network parameters in the model. When $R(0)=0$ :

(2.15) $$\begin{eqnarray}\displaystyle t_{blowup}\sim \unicode[STIX]{x1D70F}_{D}\unicode[STIX]{x1D719}^{1-d_{f}/3}\unicode[STIX]{x1D6FF}^{d_{f}/3}\unicode[STIX]{x1D716}^{3-d_{f}/3}\text{e}^{1/d_{3}\unicode[STIX]{x1D716}}G^{d_{f}/3-2}, & & \displaystyle\end{eqnarray}$$

where the only unknowns are the scalar coefficient of proportionality and the coefficient associated with Kramers hopping process, $d_{3}$ . This is a clear prediction of how the point in time where the hydrodynamic instability and collapse occur relates to properties of the gel network. These are the adjustable material parameters that are available to engineer stability into the particulate network. Before we proceed to discuss the utility of our model prediction in § 4, we test its validity using extensive computer simulations and comparisons to published experimental data in the next section.

3.1 Simulation methodology

In order to study hydrodynamic instabilities during gel collapse and observe the breakdown of the network, the effects of fluid flows and hydrodynamic forces have to be modelled accurately in large scale simulations. We have recently developed methods for rapid calculation of hydrodynamic interactions in suspensions of mono-disperse spheres (Swan & Wang Reference Swan and Wang2016; Fiore et al. Reference Fiore, Balboa Usabiaga, Donev and Swan2017), where we use the Rotneâ–Prager–â Yamakawa tensor (RPY) to account for long-ranged hydrodynamic interactions with great fidelity (Rotne & Prager Reference Rotne and Prager1969). The positively split Ewald (PSE) algorithm makes the cost of computing Brownian displacements in simulations of colloidal scale particles with hydrodynamic interactions comparable to the cost of computing deterministic displacements in freely draining simulations. The method relies on a new formulation for Ewald summation of the RPY tensor, which guarantees that the real-space and wave-space contributions to the tensor are independently symmetric and positive–definite for all possible particle configurations. Brownian displacements are drawn from a superposition of two independent stochastic samples: a wave-space (far-field) contribution, computed using techniques from fluctuating hydrodynamics and non-uniform fast Fourier transforms; and a real-space contribution, computed using a Krylov subspace method. The combined computational complexity of drawing these two independent samples scales linearly with the number of particles enabling hydrodynamic simulations with system sizes up to $10^{6}$ particles (Fiore et al. Reference Fiore, Balboa Usabiaga, Donev and Swan2017). Higher-order hydrodynamic models are not yet feasible at the scale simulated in this work. The far-field contribution captures much of the effect of bulk flow of freely settling particles within the gel, and we do not expect that the singular lubrication forces between nearly touching particles play a significant role in the dynamics, which is dominated by collective rather than relative particle motion. When superior methods for high-order simulations become available, these assumptions should be checked more thoroughly.

Extensive simulations of freely settling, attractive, hydrodynamically interacting, colloidal particles of size $a$ in a solvent of viscosity $\unicode[STIX]{x1D702}$ are performed. The simulations contain $N_{sim}=216\,000$ particles having a volume fraction $\unicode[STIX]{x1D719}$ in a cubic simulation box of length $L_{sim}$ with periodic boundary conditions. $N_{sim}$ has been selected to avoid any system size effects (Varga et al. Reference Varga, Wang and Swan2015) and to be able to resolve large scale structural changes (Varga & Swan Reference Varga and Swan2018a ). Other choices of aspect ratio, including stretched and flattened gel columns, have been investigated and do not affect the simulation results that follow. Any interactions with the container, the meniscus or other potential wall effects are neglected in the simulations. We seek to use the simulations to probe only the freely settling region of the network. The interactions with walls can influence the mode of collapse (Poon Reference Poon2002) and this effect is discussed later. However, the hydrodynamic instability we seek to model occurs far from the boundaries in the feely settling region of the gel where the leading-order effects of the sample geometry are thought to be insignificant (Secchi et al. Reference Secchi, Buzzaccaro and Piazza2014). Zero volume flux boundary conditions on the simulation box ensure that sedimentation models the free falling zone (Buscall & White Reference Buscall and White1987), where the gel network is freely settling. Note that as a consequence, both the bottom of the container and the compacting cake region, shown previously to play no role in collapse, are ignored. Furthermore the processes at the top surface of the gel and the role of the meniscus in inducing collapse are not under study. The simulation can be viewed as modelling micron sized colloids and a millimetre sized gel cross-section deep inside a sedimenting network that is in a state of free fall.

We introduce a short-ranged attraction, mimicking the polymer induced depletion attraction in experimental systems (Russel et al. Reference Russel, Saville and Schowalter1989; Poon et al. Reference Poon, Pirie, Haw and Pusey1997; Lu et al. Reference Lu, Conrad, Wyss, Schofield and Weitz2006), and model it through an Asakura–Oosawa form (Asakura & Oosawa Reference Asakura and Oosawa1958):

(3.1) $$\begin{eqnarray}\displaystyle U_{A}(r)=-U\frac{2(2a(1+\unicode[STIX]{x1D6FF}))^{3}-3r(2a(1+\unicode[STIX]{x1D6FF}))^{2}+r^{3}}{2(2a(1+\unicode[STIX]{x1D6FF}))^{3}-6a(2a(1+\unicode[STIX]{x1D6FF}))^{2}+(2a)^{3}}, & & \displaystyle\end{eqnarray}$$

for particle separations $r$ in the range of $2a<r<2(a+R_{g})$ . The polymer radius of gyration, $R_{g}$ , relative to the colloid particle size is varied, $\unicode[STIX]{x1D6FF}=R_{g}/a$ , and the pairwise depletion strength at contact sets the network strength, $\unicode[STIX]{x1D716}=kT/U$ , from the athermal limit, $\unicode[STIX]{x1D716}=0$ , to a hard-sphere dispersion, $\unicode[STIX]{x1D716}\rightarrow \infty$ . The Heyes–Melrose potential-free algorithm ensures hard-sphere repulsion at contact (Heyes & Melrose Reference Heyes and Melrose1993). The uniform gravitational load on all particles is tuned through the gravitational Mason number, $G$ . We study the settling systems over a range of attraction ranges: $\unicode[STIX]{x1D6FF}=0.075{-}0.15$ , strengths: $\unicode[STIX]{x1D716}=0.01{-}0.2$ , Mason numbers: $G=0.1{-}1.0$ and volume fractions: $\unicode[STIX]{x1D719}=5{-}50\,\%$ (see table 1 for definitions of all dimensionless quantities). Initially, the dispersion is allowed to gel for 500 bare diffusion steps, $\unicode[STIX]{x1D70F}_{D}=6\unicode[STIX]{x03C0}\unicode[STIX]{x1D702}a^{3}/kT$ in the absence of gravity. Use of the box-counting method determines the Minkowski–Bouligand dimension (Falconer Reference Falconer2004), $d_{f}$ , for each gel, which characterizes the meso-scale structure and tortuous nature of the random network. Note that due to the finite system size and finite particle size the gels are not true fractal objects over all length scales. At time $t=0$ we introduce a finite gravitational Mason number $G$ and observe the process of free settling over a simulation time of $2500G/\unicode[STIX]{x1D716}\unicode[STIX]{x1D70F}_{D}=2500\unicode[STIX]{x1D70F}_{G}$ , where $\unicode[STIX]{x1D70F}_{G}$ is the characteristic settling time, the time it takes a single particle to settle its own radius in bulk fluid. As the network moves downward, fluid back flow, particle erosion and the eventual failure of the network are observed. All reported simulation results are averaged over three independently generated samples for each combination of $\unicode[STIX]{x1D6FF}$ , $\unicode[STIX]{x1D716}$ , $G$ and $\unicode[STIX]{x1D719}$ .

One assumption made in the simulations, which is not relevant for the theory of streamer growth, is that the initial gelled state emerges from a diffusion limited aggregation process in the absence of any gravitational load. In reality, experimental gels form while settling and it is well known that with large gravitational loading, $G\gg 1$ , the structure of the gel or the gel point can be altered. However, the initial gel structure only has an impact on the streamer erosion theory through the dimensionless coefficients, $d_{i}$ , and the experiments for which data are available reside in the $G\ll 1$ regime modelled well by the initial states chosen for these simulations.

3.2 Comparison with simulations of freely settling colloidal gels

The colloidal gel networks in the simulations all exhibit gravitational collapse during free settling. While the exact time point at which this occurs strongly depends on the network parameters $\unicode[STIX]{x1D6FF}$ , $\unicode[STIX]{x1D716}$ , $G$ and $\unicode[STIX]{x1D719}$ , all gels eventually experience a hydrodynamic instability and fail catastrophically as shown in figure 2 (see the supplementary movies available at https://doi.org/10.1017/jfm.2018.725). Initially, the network settles uniformly under its own weight and fluid flows upwards through the gel pores. After an initiation period, a single streamer nucleates in the gel. This streamer grows radially as individual particles and small clusters are eroded and swept upwards with the back flow. The streamer then grows and spans the height of the settling gel. There now exists a continuous channel for fluid back flow and the entire network is destabilized as the streamer rapidly expands in the radial direction. Eventually, portions of the gel compact, network integrity loss occurs and large domains of the gel move independently both up and downwards. The gel is no more and the network has undergone catastrophic collapse.

Figure 2. In simulations we observe freely sedimenting gels in a gravitational field $\boldsymbol{g}$ and study the effects of fluid back flow over time (see § 3.1 for simulation details). These are snapshots of a simulation with $\unicode[STIX]{x1D6FF}=0.1$ , $\unicode[STIX]{x1D716}=0.05$ , $G=0.5$ and $\unicode[STIX]{x1D719}=20\,\%$ (see also supplementary movies of the collapse). The differently coloured layers are purely for illustrative purposes, indicate initial particle positions in the gel, and are meant to guide the eye through the breakdown of the network during free settling shown for a particle depth of $30a$ . The dispersion gelled over $500\unicode[STIX]{x1D70F}_{D}$ in the absence of gravity and has a structure characterized by $d_{f}=2.05$ . (a) At $t=0\unicode[STIX]{x1D70F}_{D}$ gravity is turned on in the simulation and the network begins to settle. (b) After initial uniform settling, a single cylindrical streamer nucleates (bottom centre). (c) The streamer is both growing radially and its height spans the bottom half of the network. (d) At the onset of settling rate increase ( $t=t_{crossover}$ ) the streamer spans the height of the gel. (e) The streamer continues to grow radially, filling the entire sample, destabilizing and changing the uniform settling of the network. (f) There is complete loss of network integrity as entire aggregates are ripped off the gel.

We quantify this process using two independent metrics, the evolution of the network settling velocity and the growth of the streamer volume over time. The average network settling velocity, $U(t)$ , is measured by computing the velocity of the centre of mass of the gel in the frame of an external observer, the laboratory frame as opposed to the frame of zero volume flux, the Lagrangian perspective, for all gels under study (note that we exclude the velocity of particles that are not attached to the percolated structure). Figure 3 shows this network settling velocity normalized on its initial value, $U(0)$ , as a function of dimensionless simulation time for increasing values of $G$ and constant $\unicode[STIX]{x1D6FF}$ , $\unicode[STIX]{x1D716}$ and $\unicode[STIX]{x1D719}$ . The data exhibit three distinct regimes of settling. As seen, the network begins to settle with a constant initial uniform velocity. Then at a certain point, the velocity grows as a power law in time until it enters the third regime, where it reaches a new plateau of the settling velocity with $U_{final}/U(0)\sim O(10)$ , consistent across all gels studied. Interestingly, Starrs et al. (Reference Starrs, Poon, Hibberd and Robins2002) also observed a tenfold increase in their velocities from initial settling to final collapse.

Figure 3. The average network settling rate $U(t)$ normalized by the initial rate $U(0)$ as a function of simulation time in units of the characteristic single-particle settling time $\unicode[STIX]{x1D70F}_{G}$ for increasing $G$ with $\unicode[STIX]{x1D6FF}=0.1$ , $\unicode[STIX]{x1D716}=0.05$ and $\unicode[STIX]{x1D719}=20\,\%$ . The onset of rapid growth in settling rate, $t_{crossover}$ , is the point in time where the power-law growth at a given $G$ intersects the plateau of uniform settling $U(t)=U(0)$ , marked by the two dashed lines.

In agreement with visual observations, a stronger gravitational force results in an earlier onset of the power-law growth in the settling velocity. We term the transition to power-law growth the crossover time, $t_{crossover}$ , and identify it in each simulation as the point where the best-fit power-law line intersects the initial settling velocity, as illustrated in figure 3 by the dashed lines. $t_{crossover}$ changes by orders of magnitude depending on the values of the network parameters, but eventually a transition to an increasing settling velocity and final plateau is observed for even the strongest gels studied here.

In parallel with the measurements of the network settling velocity, we track the growth of the overall streamer volume for each gel, employing an approach similar to the box-counting method. At each point in time in the simulation the simulation box is divided into a three-dimensional grid of cubic boxes of size $L_{box}$ . The number of particles in each box, $N_{box}$ is counted and the probability distribution, $P(N_{box})$ , is computed. For a nascent randomly percolated gel structure of fractal dimension $d_{f}$ , this distribution should resemble a Gaussian with a mean related to bulk volume fraction $\unicode[STIX]{x1D719}$ . The width of the distribution will be a function of $d_{f}$ . In contrast, when a streamer forms this distribution is altered significantly. The streamer is largely void of particles so that when a streamer forms, $P(N_{box})$ will show an increase as $N_{box}\rightarrow 0$ . We expect that the observed distribution of $N_{box}$ will be a superposition of a contribution from a streamer and due to the bulk of the gel. For the reported results the box size was chosen as $L_{box}=10a$ , however we found that the measured distribution was not sensitive to the exact value of box size when $L_{box}$ was bigger than the number density correlation length within the gel.

The computed distribution of box occupancy is fit to a Gaussian in the neighbourhood of the bulk gel peak, then the occupancy distribution of the streamer is inferred from the difference $P(N_{box})-P_{fit}$ . The streamer volume is computed as the expected value of the free volume (box volume minus particle volume) under the streamer occupancy distribution. That is, the streamer volume is equated to the volume of the fraction of boxes that do not belong to the intact network. As settling proceeds and the network evolves, the occupancy distribution shifts as the contribution due to the streamer overwhelms that due to the gel. Near $t_{crossover}^{volume}$ , the fraction of boxes belonging to the streamer region rapidly increases and the free volume of the streamer exhibits power-law growth in time. While the peak value of the Gaussian fit to the bulk gel decreases with time, its mean value and variance remain essentially the same over the course of the simulation. The structure of the bulk of the gel is unaltered by erosion of the streamer.

Figure 4(a) plots the evolution of the streamer volume $V(t)$ normalized on the simulation box volume for the same set of gels as in figure 3. Initially there is no noticeable streamer present. At a certain point in time, which decreases with increasing $G$ , $V(t)$ exhibits a power-law growth, very similar to the behaviour exhibited by $U(t)$ . Again, it is possible to extract a crossover time, $t_{crossover}^{volume}$ , using the same method as with the settling rate. In figure 4(b) we plot $t_{crossover}^{volume}$ versus $t_{crossover}$ and find that the two time scales are identical to within the measurement errors. The acceleration of the settling network and its ultimate collapse is directly correlated with the nucleation and subsequent growth of the streamer in the gel in accordance with the mechanism described by the model.

Figure 4. There is a direct correlation between the increase in streamer volume and the accelerating settling network, which supports the model’s premise that streamer nucleation and growth are causes for loss of network integrity in settling gels. (a) The growth of streamer volume normalized on the simulation box volume plotted as a function of simulation time for increasing $G$ with $\unicode[STIX]{x1D6FF}=0.1$ , $\unicode[STIX]{x1D716}=0.05$ and $\unicode[STIX]{x1D719}=20\,\%$ . (b) The onset of rapid growth in streamer volume $t_{crossover}^{volume}$ plotted versus the onset of power-law growth in network settling velocity, $t_{crossover}$ .

The crossover time measured in simulations marks the beginning of the collapse of the gel, after which the network rapidly loses its integrity. While the dynamic simulations exhibit the same qualitative behaviour as described by the model, it remains a question whether (2.15) can predict how $t_{crossover}$ depends on the network parameters. Figure 5 examines the dependence of $t_{crossover}$ on $G$ , $\unicode[STIX]{x1D719}$ , $\unicode[STIX]{x1D6FF}$ and $\unicode[STIX]{x1D716}$ separately.

Figure 5. The onset of rapid network collapse, $t_{crossover}$ as measured in dynamic simulations is plotted individually as a function of the network parameters $G$ , $\unicode[STIX]{x1D719}$ , $\unicode[STIX]{x1D6FF}$ and $\unicode[STIX]{x1D716}$ respectively, while all others are held constant. The collapse dynamics exhibit the power-law behaviour predicted by the micromechanical model for networks with $d_{f}\approx 2$ . Each data point is the average of three independent simulations and error bars represent 95 % confidence intervals. The slope of the best-fit line in (d) yields the scaling coefficient $d_{3}$ that sets the relevant strength of attraction, $d_{3}\approx 5.6$ .

Figure 5(a) shows the effect of increasing gravitational Mason number for three different network strengths and volume fractions. Also shown is the expected scaling behaviour as given by (2.15) for a network with fractal dimension $d_{f}=2$ : $t_{blowup}\sim G^{d_{f}/3-2}=G^{-4/3}$ . Indeed for a range of gravitational Mason numbers for all three conditions the crossover time exhibits the scaling that the model would predict. Small deviations from the $-4/3$ scaling are expected as the measured $d_{f}$ for the gels in these simulations range between 1.7–2.3. For large gravitational forces ( $G>1$ ) a different mechanism controls the network collapse. The value of $t_{crossover}$ appears independent of the gravitational Mason number and the situation is one of weakly aggregated clusters settling freely (Huh et al. Reference Huh, Lynch and Furst2007).

Next the dependence of $t_{crossover}$ on the volume fraction is analysed for two different combinations of $\unicode[STIX]{x1D6FF}$ , $\unicode[STIX]{x1D716}$ and $G$ along with the expected scaling behaviour from the model, assuming $d_{f}=2$ : $t_{blowup}\sim \unicode[STIX]{x1D719}^{1/3}$ . The crossover time exhibits roughly the model behaviour for a large range of volume fractions, as shown in figure 5(b). As the fractal dimension is especially sensitive to $\unicode[STIX]{x1D719}$ , it is no surprise that $t_{crossover}$ deviates slightly from the $\unicode[STIX]{x1D719}^{1/3}$ scaling. The model prediction breaks down for $\unicode[STIX]{x1D719}>30\,\%$ , at which point the starting material cannot be considered a gel, but instead approaches the properties of a colloidal glass. Clearly, other restructuring processes dominate settling in these dense structures (Zaccarelli Reference Zaccarelli2007).

Similarly, the critical time exhibits the model predicted dependence on the attraction range, $t_{crossover}\sim \unicode[STIX]{x1D6FF}^{d_{f}/3}$ as shown in figure 5(c). For the values of $\unicode[STIX]{x1D6FF}$ investigated here, the law of corresponding states suggests that the thermodynamic restoring forces for these dispersions will all be very similar (Noro & Frenkel Reference Noro and Frenkel2000). To maintain the short-range nature of the colloidal bonds however, $\unicode[STIX]{x1D6FF}$ could not be increased above this narrow range since it known that the evolution of the collapse dynamics is markedly different for gels with long-range attractions (Teece et al. Reference Teece, Faers and Bartlett2011). So, agreement with the model should be considered tentative.

Finally, we consider the effect of the network strength $\unicode[STIX]{x1D716}$ on $t_{crossover}$ . Since (2.15) suggests both a power-law and exponential dependence on $\unicode[STIX]{x1D716}$ , figure 5(d) plots $\log (\unicode[STIX]{x1D716}^{d_{f}/3-3}t_{crossover}/\unicode[STIX]{x1D70F}_{D})$ versus $\unicode[STIX]{x1D716}^{-1}$ . Note, in the model, the proportionality constant modulating the strength of attraction, $d_{3}$ is left undetermined. Indeed there appears to be a linear relation between the time scale of collapse of the network and the Kramers escape rate, as previously observed experimentally by Teece et al. (Reference Teece, Hart, Hsu, Gilligan, Faers and Bartlett2014). The slope of the best-fit line for all three combinations of $\unicode[STIX]{x1D6FF}$ , $G$ and $\unicode[STIX]{x1D719}$ indicates that the best choice for the constant is $d_{3}=5.6$ and the crossover time exhibits the relationship with $\unicode[STIX]{x1D716}$ that the model predicts. While we cannot provide a physical explanation for this particular value, it is a result of the simplified approach to approximate the network erosion in terms of the bond breaking rate or escape probability of individual particles from the attractive well. For very large attractions, $\unicode[STIX]{x1D716}\leqslant 10^{-2}$ a small deviation from the predicted scaling is observed. Section 4.2 discusses possible reasons for this and ways to improve the model.

Figure 6. The onset of rapid network collapse, $t_{crossover}$ as measured in dynamic simulations is plotted as a function of the model prediction $t_{blowup}$ , computed for the combination of $\unicode[STIX]{x1D6FF}$ , $\unicode[STIX]{x1D716}$ , $G$ , $\unicode[STIX]{x1D719}$ and $d_{f}$ used in the respective simulation run. Each data point is the average of three independent simulations with the same unique combination of parameters in the parameter space and error bars represent 95 % confidence intervals. The model predicts the collapse dynamics to within a scalar coefficient. A linear best fit through the data (dashed line) is employed to identify the coefficient.

It appears that the phenomenological model adequately predicts the dependence of the critical time scale in dynamic simulations on the dimensionless parameters and captures the essential features of collapse. Equation (2.14) gives a quantitative prediction for the critical blowup time to within a scalar constant, provided only the values of the five dimensionless network parameters. In figure 6 we plot the measured $t_{crossover}$ versus $t_{blowup}$ and find a direct parity between the two quantities. Here, every data point is an individual combination of $\unicode[STIX]{x1D6FF}$ , $\unicode[STIX]{x1D716}$ , $G$ and $\unicode[STIX]{x1D719}$ and the calculated $d_{f}$ , averaged over three independent simulation runs. Given this linear relationship, the slope of the best-fit straight line through the data can be used to obtain the missing scalar constant for the model, arriving at the final result:

(3.2) $$\begin{eqnarray}\displaystyle t_{blowup}=540\unicode[STIX]{x1D70F}_{D}\unicode[STIX]{x1D719}^{1-d_{f}/3}\unicode[STIX]{x1D6FF}^{d_{f}/3}\unicode[STIX]{x1D716}^{3-d_{f}/3}\text{e}^{1/(5.6\unicode[STIX]{x1D716})}G^{d_{f}/3-2}. & & \displaystyle\end{eqnarray}$$

For a colloidal gel (3.2) can be applied to calculate the point in time where the hydrodynamic instability occurs, leading to rapid collapse.

3.3 The effect of finite initial channel radius

During the analysis of the computational model results presented so far, when relating $t_{crossover}$ to $t_{blowup}$ it was implicitly assumed that $R(0)=0$ and that a streamer is nucleated immediately after the start of the simulation. However, as will also be discussed in § 4.2 an initial vertical channel might already be present in the gel network at $t=0$ . Or, on observing a gel undergoing collapse, the initial starting point itself may not be well defined in an actual system, especially since in experiments colloidal dispersions do not gel independently from the influence of the gravitational field they are in. Regardless, the model permits a solution for $t_{blowup}$ when assuming an initial condition $R_{0}\neq 0$ for (2.10). Intuitively, it is expected that catastrophic collapse in a gel with a channel present will occur sooner than in an unperturbed gel. Indeed, using the result in (2.14) the model prediction for the shortened finite-time singularity relative to the unperturbed case is found to be:

(3.3) $$\begin{eqnarray}\displaystyle \frac{t_{blowup}|_{R_{0}}}{t_{blowup}|_{0}}=\frac{\unicode[STIX]{x1D6E4}(2-d_{f}/3,R_{0}/R^{\ast })}{\unicode[STIX]{x1D6E4}(2-d_{f}/3,0)}. & & \displaystyle\end{eqnarray}$$

The ratio of gamma functions is guaranteed to be less than unity, regardless of the value of $d_{f}$ , in agreement with the intuitive expectation. The only controlling parameter is the value of the initial channel radius $R_{0}$ relative to $R^{\ast }$ and this ratio is predicted to be independent of $\unicode[STIX]{x1D719}$ .

To test this prediction we conduct an additional set of simulations of freely settling gels where a vertical cylindrical channel of height $L_{sim}$ and radius $R_{0}$ devoid of particles is seeded at the centre of the gel at $t=0$ . The parameter $R_{0}/R^{\ast }$ is varied systematically. This is achieved both, by increasing $R_{0}$ relative to the overall size of the gel, and separately, by decreasing $R^{\ast }$ . Remember, the gravitational length depends on $\unicode[STIX]{x1D6FF}$ , $\unicode[STIX]{x1D716}$ and $G$ so that the effect of varying network parameters is directly incorporated. Because of the exponential growth, in order to resolve the decrease in blowup time, we explore a parameter range of two orders of magnitude in $R_{0}/R^{\ast }$ . As long as the initial pore is small relative to the system size, the assumptions of our model should hold and the onset and dynamics of collapse should only be determined by the radius of the streamer. In each simulation the crossover time with a seeded pore is measured and compared to the corresponding crossover time in an unseeded gel at the conditions. We plot their ratio in figure 7 along with the prediction of (3.3).

Figure 7. The measured values of $t_{crossover}$ in the presence of a seeded hole at $t=0$ relative to the time scale of an unperturbed sample are plotted as a function of the relative seeded channel radius $R_{0}/R^{\ast }$ for different $\unicode[STIX]{x1D719}$ . $R^{\ast }$ is set by a combination of $\unicode[STIX]{x1D6FF}$ , $\unicode[STIX]{x1D716}$ and $G$ . Each data point is the average of three independent simulations with the same unique combination of parameters in the parameter space and error bars represent 95 % confidence intervals. The measured values are in good agreement with the model prediction of (3.3) using a value of $d_{f}=2$ (dashed line).

We find an excellent agreement between (3.3) and the simulation results. For small seeded channel sizes, $R(0)\ll R^{\ast }$ , the collapse times are nearly the same with only a negligible decrease, as stated in § 2.2. In contrast, for $R_{0}/R^{\ast }~O(1)$ , $t_{crossover}|_{R_{0}}$ falls off due to the activated hopping of particles off the network, as expected from the theory. Simulations show that the seeded channel provides a free path for fluid back flow when the gel is settling. Particles on the network–streamer interface are immediately ripped off the network and swept upwards. This provides further evidence that the growth of a streamer in a settling gel is indeed the cause for the rapid instability leading to collapse. As we have shown, the predictions of the simple model of streamer growth agree well with the collapse dynamics observed in simulations for a large range of parameter values and initial gel states.

3.4 Comparison with experimental observations

So far it was shown that the model correctly captures the mechanism of collapse seen in dynamic simulations. Here we show that our theory also describes the gel collapse behaviour in two vastly different experimental systems by fitting available data of the evolution of the observed streamer radius as a function of time to our model.

Starrs et al. (Reference Starrs, Poon, Hibberd and Robins2002) studied a system of polydispserse poly(methyl methacrylate) (PMMA) particles that were induced to gel in the presence of polystyrene depletants in a tetralin and cis-decalin solvent blend. In so-called ‘weak’ gels streamers formed, and the gel exhibited catastrophic collapse following a hydrodynamic instability. The weak gels were differentiated from ‘strong’ gels considered in the study by the lack of hydrodynamic instability and the observation of steady, poro-elastic compression rather than a dramatic collapse. Secchi et al. (Reference Secchi, Buzzaccaro and Piazza2014) looked at an aqueous dispersions of spherical particles of MFA, a copolymer of tetrafluoroethylene and perfluoromethylvinylether. Depletion interactions were induced by the addition of a surfactant that forms globular micelles leading to a short-ranged attraction. Here, the onset and subsequent radial growth of streamers destabilizing the network were also observed. In both cases, snapshots depicting the increasing size of the streamer with time were included by the authors, providing two experimental data sets for the evolution of the streamer radius, $R(t)$ . Table 2 provides the parameters for the two experimental systems.

Table 2. Experimental parameters of collapsing gels as found in Starrs et al. (Reference Starrs, Poon, Hibberd and Robins2002) and Secchi et al. (Reference Secchi, Buzzaccaro and Piazza2014).

Assuming an initial pore size, $R(0)=0$ , equation (2.10) has the solution:

(3.4) $$\begin{eqnarray}\displaystyle \frac{t_{blowup}-t}{t_{blowup}}=\frac{\unicode[STIX]{x1D6E4}(2-d_{f}/3,R(t)/R^{\ast })}{\unicode[STIX]{x1D6E4}(2-d_{f}/3,0)}, & & \displaystyle\end{eqnarray}$$

which has three parameters: $t_{blowup}$ , $R^{\ast }$ and $d_{f}$ . As seen earlier, the model predictions are not sensitive to the value of $d_{f}$ in the range of 1.7–2.3 and so in the absence of any further information from the experiments, we assume that $d_{f}=2$ for both networks, a typical value for dilute colloidal gels (Zaccarelli Reference Zaccarelli2007). Therefore, it remains to obtain a best fit of (3.4) to the two experimental data sets to extract the parameters $t_{blowup}$ and $R^{\ast }$ . Figure 8 compares the experimentally observed streamer radius to the best-fit model predictions. The measured dynamics for both experimental systems follows the model quite well, and table 3 provides the best-fit values for $t_{blowup}$ and $R^{\ast }$ . While $t_{blowup}$ and $R^{\ast }$ are not reported directly in either paper, these best-fit values match well the final collapse time for the gel in either experiment, as shown in table 2. Similarly, the characteristic radius in either case is an order of magnitude smaller than the width of the experimental gel columns, which supports the claim that rapid collapse begins suddenly as the streamer rapidly expands beyond $R^{\ast }$ . For reference, the supplementary section contains images of the experimental gels with outlines demarcating the pore at different points in time.

Figure 8. Time evolution of streamer radius $R$ as extracted from experimental data published in the literature and the best-fit predictions of our model. (a) Comparison of data published by Starrs et al. (Reference Starrs, Poon, Hibberd and Robins2002) (open squares) and the best-fit prediction of the model (dashed line). (b) Comparison of data published by Secchi et al. (Reference Secchi, Buzzaccaro and Piazza2014) (open circles) and the best-fit prediction of the model (dashed line).

The quantities $t_{blowup}$ and $R^{\ast }$ can be estimated from experimental parameters independent of the model fit. Using the values of the experimental quantities in table 2, $t_{blowup}$ and $R^{\ast }$ are computed directly using (2.6) and (3.2) assuming again $d_{f}=2$ and taking the coefficient $d_{4}=1$ . These estimates are reported in table 3. In the case of the work by Starrs et al. (Reference Starrs, Poon, Hibberd and Robins2002) the best-fit values for both parameters are in good agreement with what can be computed a priori. This would suggest that the proposed model subsumes all essential factors contributing to streamer growth and can accurately describe the dynamics of network erosion and collapse. In contrast, the estimates of $R^{\ast }$ and $t_{blowup}$ for the other example are an order of magnitude larger than the corresponding best-fit values. The value of $R^{\ast }$ exceeds the dimensions of the gel and $t_{blowup}$ exceeds the observation window reported (Secchi et al. Reference Secchi, Buzzaccaro and Piazza2014). Since the dynamics is well captured by the best fit, this would indicate that there is no seeded pore in the gel. Instead, we conclude that there is some uncertainty in the experimental parameters. While we have assumed that the characteristic building blocks of the gel in the experiments are individual MFA particles, Secchi et al. (Reference Secchi, Buzzaccaro and Piazza2014) note that there is strong evidence to suggest that the particles move in aggregated clusters. As the model displays a high sensitivity to the size of particles within the gel (see § 4.1) the value of the hydrodynamic radius $a$ can significantly impact the quality of the predicted blowup time. Indeed, using $a=3\times 90$  nm as the characteristics size, we recover estimates of $R^{\ast }$ and $t_{blowup}$ , shown in the last row of table 3, that are both physically reasonable and agree with the best-fit values.

Table 3. Values of $t_{blowup}$ and $R^{\ast }$ as found from the model best fit to the experimental data and their estimates based on details of the experimental parameters in table 2. In the case of Secchi et al. (Reference Secchi, Buzzaccaro and Piazza2014) estimates are based on the primary particle radius and three times the size for $a$ .

As a result of the exponential growth rate of the streamer radius, the value of $R(t)$ and the dynamics in the vicinity of blowup are very sensitive to $t$ and change rapidly, as seen in figure 8. While $t_{blowup}$ and $R^{\ast }$ depend on the experimental parameters and are unique to each study, (3.4) suggests that $\unicode[STIX]{x1D6E4}(2-d_{f}/3,R(t)/R^{\ast })$ will exhibit a universal linear dependence on the distance to blowup, $(t_{blowup}-t)$ . Figure 9 plots the experimentally measured streamer radii as a function of distance to blowup, and we observe indeed a one-to-one parity as expected from theory. This indicates that the model for streamer growth contains the necessary dynamics to understand and track the hydrodynamic instability leading to the collapse of the network in both colloidal gels. As a final note, the experimental values of $G$ that we estimate from the experimental parameters are two orders of magnitude smaller than what we have studied via simulations $O(10^{-3})$ . Recalling that the blowup time scales approximately as $G^{-4/3}$ , simulations in which pores emerge and grow at such small $G$ would be too time consuming to be performed at present. However, we believe that strong connection between the simulations and the model support application of the model to these experimental results.

Figure 9. Parity plot of the predicted and the measured distance to blowup (see text for details) for the two very different experimental systems. The theory is able to capture the underlying collapse dynamics well, independent of the experimental details.