An account is made of the experimental, theoretical, and computational developments that led to our current understanding of the colloidal aggregation problem when a gravitational field is present. Starting with unaggregated colloids, a review is made of the advances that led to the founding of the barometric equation for the distribution of colloidal particles in a suspension, noticing that for large bodies, like large colloidal aggregates, their final fate in equilibrium is to be at the bottom of the container. Then, we briefly review the aggregation of colloids in the absence of gravity that has been amply studied by both experiments and simulations. For this purpose, the paradigmatic case of the DLVO interaction is taken as an example. Next, a brief revision is made of the seminal experimental work of C. Allain and collaborators on the colloidal aggregation problem when an external gravitational field is present, centering our study in the nongelling situations, that is, for dilute colloidal suspensions, when only sedimentation and deposition of single clusters occur. Afterward, the development of different computer simulations that treat this case of single cluster sedimentation and deposition is reviewed, and note how the different improvements of the algorithms lead to better correspondences with the experimental systems. We finally discuss further possible improvements of the algorithms and end with proposals for future work.
- colloidal particles
- Brownian motion
- colloidal sedimentation
- colloidal aggregation
- cluster sedimentation
- cluster deposition
Small colloidal particles—of the order from 10 nm up to 10 μm in size—immersed in a liquid experience a motion, which is the combination of a steady gravitational drift downward (assuming that their density is higher than that of the surrounding liquid, otherwise the drift would go upward) and a Brownian random movement. This is due to its translational kinetic energy coming from the principle of equipartition of energy, which, as realized by Einstein, does not distinguish the thermal motion of a solvent molecule from that of a suspended colloid or any other types of particle. The kinetic energy of a particle with mass m, translating with a speed v is
The equipartition principle guarantees that in thermal equilibrium, all components of a solution (solvent molecules as well as colloids or any other particles) have the same average translational kinetic energy, which is fixed by the absolute temperature T:
Thus, the root-mean-square speed of a particle is
showing that at a given temperature, a colloid with a large mass moves much slower than a molecule. This movement of the colloidal particles occurs in random (erratic) steps, whose finite size comes from the damping of their movement by the surrounding fluid. In the free diffusion of a particle in a liquid, far away from other particles or a wall, the diffusion coefficient D is given by Einstein’s equation :
where f is the friction coefficient of the particle. The simplest case is a spherical particle with radius a in a Newtonian liquid with viscosity η. In this case, f is the so-called Stokes friction factor: f = 6πηa. The combined result
is called the Stokes-Einstein diffusion coefficient for translational sphere diffusion. With this definition for the diffusion coefficient, Einstein reported the equation for the average quadratic displacement of a particle as follows:
If in addition to this Brownian movement a gravitational field along the vertical z direction is present in the system, the particles experience a gravitational force leading to a final equilibrium distribution of the particles which is dependent on the z coordinate. With the use of the methods of statistical mechanics, this equilibrium average density profile is well known to be given by
where z is the vertical position measured from the free surface of the suspension fluid, considered as the origin of the upward z coordinate. Here c is the colloid number density. This barometric distribution, which only holds for noninteracting particles, has a thickness characterized by the “gravitational length”: lg = kBT/mg. It is easily verified that for common gas molecules as nitrogen, lg is of the order of 9 km, while for colloidal spheres of 0.1 μm in diameter with a density of, say, ρ = 2 g/cm3, it is only around 1 mm, reducing its mass by the corresponding mass of the displaced fluid. For a given system of initial, homogeneous number density co = Ntot/Ah inside a prismatic cell of area A and height h, c(z = 0) is given by
Perrin started the use of “well-defined” colloids to study the density profile on a spatial scale accessible to an optical microscope [2–6]. It should be noted that Eq. (7) assumes that all colloidal particles have the same mass m. He used fairly monodisperse spheres, obtained from laborious fractionation procedures by centrifugation of solutions of the natural gum “gamboge.” In Figure 1, it is shown a Perrin’s microscopic image of the sedimentation-diffusion equilibrium of these resin spheres of about half a micron in diameter immersed in water. His method to obtain the exponential density profile was to compare the average number of particles in a given layer of the sample to the same average number in another layer some tens of microns below. For example, from Eq. (7), we can write
In one series of measurements, he used gamboge particles of around 0.24 μm in diameter with a density of ρ = 1.35 g/cm3 (which would imply a reduced density of ρ = 0.35 g/cm3). This would proportion the quantity inside the exponential with the value of −0.15, which was relatively not very difficult to measure. He used cells of a very short height (100 μm). In one of Perrin’s works , he warned very clearly: “… One must use capillary tubing to avoid convection currents which would muddle the surface of the cloud, and place them vertically in a thermostat. These currents are produced with extreme ease in large tubes, coming from the temperature gradients.”
As a final remark on the experimental work, if we combine Eqs. (7) and (8) and use the values of the density and radius of Perrin’s particles given above, the quantity α = mgh/kBT becomes as large as 600 for a usual laboratory vial of around 10 cm in height. Hence, in equilibrium the vial should be mostly empty of particles, except close to the bottom where the number density should follow the law c = co α exp(−α Δ/h) ; here Δ is the vertical distance from the bottom to the measuring point.
Subsequently, Mason and Weaver  investigated theoretically the manner in which the system approaches the equilibrium exponential distribution, not only for a constant uniform initial distribution co but for an arbitrary initial distribution co = f(z). A consideration of the combined effect of downward drift and diffusion led them to obtain a partial differential equation for the number density of particles as a function of depth and time. Their solution showed that no matter what the initial distribution was, the final, long-time distribution was the one shown in Eq. (7). Perrin’s exponential distribution was questioned by several researchers, among them Burton and Bishop , since Perrin performed his measurements very close to the free surface of the suspension (of the order of a hundred microns), arguing that the distribution was merely a surface phenomena that hold for only very small distances below the free surface. They used instead a tube of a height of around 1 m  containing copper colloids, waiting for about 50 days for the equilibrium to be established and concluded that “… these suspensions of fine particles always reach a uniform distribution throughout the body of the liquid, and the variation expressed in Perrin’s distribution law must be confined to a very small distance at the surface.” The theoretical expression derived by Mason and Weaver  was used by Weaver  to try to settle the issue. He gave an estimate of the duration of the transient state for Burton’s particles and tubing, before the final equilibrium distribution was reached. He found that one must wait for the time that a drifting particle sedimenting downward would travel of the order of two tube lengths, in order for the transient term in their expression to become negligible. However, he concluded that the time of 50 days spent by Burton was long enough to attain the equilibrium distribution, according to his calculations. To make a long story short, while the kinetic-molecular theory became widely accepted, the pesky colloids refused to settle in the laboratory bottles and many wondered why. Svedberg  pointed out that “… the discrepancy was probably due to the mixing effect of convection currents. For very small particles, the velocity of fall is extremely small and even slight temperature differences within the sol would probably disturb the equilibrium considerably.” The solution to the puzzle came years later and was described in a paper by McDowell and Usher , using a very fine tubing and temperature-controlled apparatus, in which the temperature variations did not exceed 0.001°C. They found that the exponential law held for not only a tenth of a mm, as found by Perrin, but up to almost 1 cm, the length of their tubing apparatus, and no sensible departure from this ideal behavior was apparent. Finally, it should be mentioned that for very large particles, like the colloidal aggregates, we are going to study (sometimes consisting of up to 10,000 colloids) the quantity α becomes very large, around 6 × 106. Needless to say that in equilibrium those aggregates would not be found in the bulk of the suspension but at the bottom.
If, furthermore, the colloids experience not only the Brownian motion and the gravitational settling but also they are allowed to aggregate or coagulate, as many colloids do under the appropriate conditions, the situation grows in complexity. Although the aggregation of colloids has been studied for over a century [12, 13], it has been in the last three decades—after Witten and Sander  advanced a single aggregate diffusive model now called DLA, but more appropriate for our study of colloidal aggregation, after the proposal of a diffusive model for cluster-cluster aggregation, made concurrently by Meakin  and Kolb et al. —that the number of investigations dealing with the clustering of particles has increased considerably. This growing was enhanced by the discovery that the colloidal aggregates present a spatial (scale invariance) fractal behavior . Although there are different ways in which the colloidal particles can aggregate, the paradigmatic and more common way to start the subject is with the so-called DLVO theory (after Derjaguin, Landau, Verwey, and Overbeek) of colloidal stability [18, 19]. In this theory, the colloids aggregate due to the strong attractive interactions when they are close to contact each other. This attraction comes from the London-van der Waals dipole-dipole forces between the molecules composing the two colloidal particles. Stability occurs when the repulsion of the two colloids, due to charges of the same sign on their surfaces, prevents the particles to get close enough to allow the dipole-dipole forces to take over. However, this stability can be disrupted by the addition to the colloidal suspension of a good amount of ionizable salt, which imposes the formation of a so-called double layer of ions around the macroion (colloidal particle). This double layer screens very effectively the electric repulsion between the colloids and permits their close approach, in such a way that they can get hooked by the London-van der Waals forces. In Figure 2, the form of the interaction potential between two of these macroions is presented, for different amounts of the added salt. Note how in some of the curves there are two minima: The first and deeper minimum (where coagulation occurs) and a secondary minimum, separated by a potential barrier. In the curve marked (D), the barrier has almost disappeared and the two particles can coagulate as soon as they collide; it is said by the computer simulators that the sticking probability after contact is one. This is the rapid coagulation regime. On the contrary in the cases (A) and (B), particularly in case (A), many collisions between the particles are necessary in order that a thermal fluctuation may help the two colliding particles to overcome this barrier. In this case, we say that the sticking probability after encounters is very low (very close to 0). We are now in the slow coagulation regime. In the olden literature on this subject (see, e.g., Ref. ), this is considered by the Fuchs’ stability factor, which is very roughly inversely proportional to the sticking probability used by computer simulators. The first case is called diffusion-limited cluster aggregation (DLCA), while the second is termed reaction-limited cluster aggregation (RLCA). One last thing about the DLVO interaction potential: since the primary minimum where coagulation occurs is very deep and of very short range, once two particles coagulate, the formed bond stays fixed with only a slight possibility and range for the particles to roll over on top of each other due to the roughness of their surfaces, that is, the real interaction potential should be actually anisotropic very close to contact. This forces the formed bond to be rigid, up to a certain degree, and this degree allows for the elasticity and deformability of a large colloidal cluster subjected to some stresses, which depends on the nature of each particular colloid we are dealing with. This deformability (including cluster fracture or breakup) has been studied in a number of works by several researchers, who considered the aggregation of clusters subjected to shear flow [21–25].
Let us call the number of aggregates in the sample containing s colloidal particles at time t as Ns(t); this would proportion the number density of such clusters as cs(t) = Ns(t)/V, where V is the volume of the sample. Therefore, the total number density of the monomeric units will be given by . Smoluchowski [12, 13] presented the first theoretical description for the rate of change of the quantities cs(t) (s = 1, 2, 3 …) in an equation that resembled the law of mass action in chemical reactions, named after him:
The quantities ki,j are called the reaction rate constants or simply kernels. Furthermore, he was able to derive an expression for the rate constant k1,1 that proportions the conversion of monomeric colloidal particles to dimers, on the assumption that this conversion strictly came from the collision rate of monomers diffusing (performing random walks) before the collision. His expression is k1,1 = 8kBT/3η. It is customary to define the Brownian collision time as tB = 2/cok1,1 = 3 η / 4 kBT co , which is the average time elapsed between two consecutive collisions of monomeric colloids in the sample, at the beginning of the aggregation. As it stands, the Smoluchowski equation with the k1,1 given above would be valid for particles sticking at first contact after diffusing, that is, for the DLCA mechanism. However, many researchers have proposed many other forms for the rate constants that may be valid for other aggregation mechanisms, including RLCA. In the present paper, we will mainly concentrate and discuss the DLCA case coupled with sedimentation.
The colloidal aggregation problem when no gravitational field is present is now relatively well understood, profiting from many results coming from the computer simulations, for both DLCA and RLCA [26–30]. In this case, the experimentalists need to resort to a number of tricks, like flipping up and down the coagulation cell quite often during the aggregation process or considering a suspension fluid whose density closely matches that of the colloidal particles (e.g., with a mixture of water and heavy water at the correct proportions), in order to eliminate the effect of gravity. Let us now briefly review our current knowledge of the DLCA process in three dimensions [31, 32]. (a) In the dilute limit, the fractal dimension df of the clusters is about 1.80. For more concentrated suspensions, it increases from this value with an added term depending on the square root of the particle concentration. (b) The number-average and weight-average cluster sizes grow linearly with time in the dilute limit Sn(t) ~ tz´ and Sw(t) ~ tz, with z = z´ = 1. Once again, the exponents z and z´ increase from one as a square root of the particle concentration. (c) The cluster size distribution Ns(t) is asymmetrically bell-shaped in the late stages of the aggregation process.
2. Colloidal aggregation coupled with sedimentation
When cluster diffusion and aggregation are coupled to the sedimentation experienced mainly by the large clusters, there appear much richer, nontrivial phenomena that proliferate in many natural and industrial processes, as opposed to the vastly investigated problem of aggregation of colloids induced merely by diffusion. Out of the numerous instances, we can mention the clarification of liquids, the aggregation and deposition of asphaltenes in crude oil, the settling of bacteria clusters in quiet water, and a number of precipitation techniques employed by the chemical industry. As quoted, studies of situations where aggregation is due to the combined action of Brownian motion and an imposed external field like gravity were scarce more than two decades ago [33, 34]. It has been only in the last 20 years that researchers have paid an increasing attention to such phenomenon. As noted before, with the addition of salts or flocculants, which screen the electrostatic repulsion between the particles or which establish bridges between them, the aggregation of particles into clusters is induced. These clusters in turn collide with other clusters, stick together, and become larger. At the beginning, the aggregates are small and essentially move by diffusion, with only a very slight sedimentation velocity. As the aggregation proceeds, the settling velocity of the large aggregates becomes significant, and the mechanism for the movement of these large aggregates by sedimentation intensifies.
2.1. The first meaningful experimental studies
A significant advance on the understanding of the process of aggregation coupled with sedimentation came after the experimental work by Allain and collaborators [35–41]. One interesting result appeared in Ref. , in which it was found that the clusters settling in a high cylindrical cell of length 0.8 m and 4.4 cm in diameter can grow up to a certain maximum cluster size smax that cannot be exceeded, with a linear size of around lmax = 0.5 mm. This was explained in the following way: at the very beginning of the growth process, the aggregates are too small to be sensitive to hydrodynamic stresses; however, as soon as they grow beyond some limiting size and their settling velocity becomes important, the hydrodynamic force exerted by the solvent exceeds the colloidal binding force. In this way, when two large clusters come into contact, the aggregates are pulled apart and go on moving separately. In all their works, they used the same aqueous suspensions of calcium carbonate particles with a density of ρ = 2.7 g/cm3 (which implies a reduced density of 1.7 g/cm3) and radius of 0.035 μm. Enough salt was added to the suspension in order for them to be in the rapid coagulation regime (DLCA). A theoretical estimation of the lmax for these aggregating particles was made and found to be around lmax = 0.1 mm. Furthermore, by measuring the settling velocity of the large clusters, they found that the fractal dimension of the aggregates grows from its DLCA value of 1.8 to a value of df ≅ 2.2–2.3. Another very interesting result was published in Ref. , where they reported a graph separating the different regimes of growth. This plot is shown in Figure 3, where on the vertical axis, it showed the logarithm of the colloidal volume fraction (φ), while on the horizontal axis, it displayed the logarithm of the time.
Allain et al. [37, 38] explained the large fractal dimension of the big, sedimenting clusters as coming from a restructuring mechanism due to the hydrodynamic stresses felt by those large clusters, forcing their long branches to elastically bend and form loops, which would make them more compact. Furthermore, they found experimentally for their system that φ* as of the order of φ* ≅ 3 × 10− 3. Let us try to obtain a rough, simpler estimate than the one obtained by them for φ*, from their data. As the cluster size distribution is bell-shaped for DLCA (sticking probability equal to one), with a typical cluster size s corresponding to the position of the maximum, we will assume momentarily that all the clusters have this size. In the following equations, we are going to make a scaling analysis, using always the sign ~ to indicate dependence and order of magnitude only, forgetting about constant factors like (4/3)π. Therefore, the results are expected to give the order of magnitude only. The effective volume fraction of all the clusters filling up the aggregation cell is , where Rs is the typical radius. By the definition of the fractal dimension df, we have . Therefore, . To obtain the volume fraction at which the system starts to gel, we need to make φe ~ 1, which implies . Now, (A) if there was no gravity involved in the problem, gelation would start when there is a single aggregate filling up the whole cell of linear size L, that is, Rs ~ L. Hence, . Assuming the size of their particles (a = 0.035 × 10−6 m) and a typical linear size of the coagulation cell of L = 0.1 m and taking the df = 1.8 for the DLCA case, we get φ*~ 1.8 × 10− 8. This is of the order of the very small φg concentration found in Ref. . (B) If the system is under the action of the gravitational field, there is not going to be a single aggregate filling up the whole system, because the clusters cannot grow beyond an Rmax, which would be taken as half of the linear size of the maximum size cluster they found. That is, Rmax ~ lmax/2 ~ 0.25 mm, which would imply . Assuming in this case that the fractal dimension is df = 2.2, as found in Ref. , we get φ* ~ 0.83 × 10− 3, a value which is close to the φ* ≈ 3 × 10− 3 found by the authors. In the rest of the present article, we will focus mainly on concentrations below the φ*; that is, the colloidal aggregation coupled with sedimentation will be studied for single clusters diffusing and sedimenting downward and depositing at the bottom. Nevertheless, if in some cases a concentration was picked to reach gelation, the simulations were stopped at the gel point.
2.2. The first computer simulation studies
One of the first attempts to study the colloidal aggregation coupled with sedimentation was presented in Ref. , although with some nonessential errors but also with one essential misassumption. It started by considering the sedimentation velocity vs experienced by a cluster of N spherical particles of radius a and mass mo
where ρo is the particle density, ρ is that of the suspension fluid, f = 6πηRg is the cluster´s friction coefficient, D(=kBT/f) is its diffusion coefficient, Rg is its radius of gyration, η is the solvent viscosity, and T is the temperature. The spirit of the simulation rests on finding out, for each cluster of a certain size, by how much the cluster settles downward as compared to how much it diffuses around. Let us define to as the time for a cluster to diffuse a particle diameter (d = 2a), that is, to = 2a2/D. This was a nonessential error which will not change the nature of the results. The correct expression in the three-dimensional case is to = 2a2/3D that has been used in the future simulations (see below). During the same time, the cluster drifts a distance , where is the Peclet number of the individual colloidal particles in the fluid. It is found that Pe is of order unity if the particles are 1 μm in diameter, ρo − ρ is of order 1 g/cm3, and T is room temperature. On the other hand, if the diameter is 0.1 μm, such quantity is of the order of 10− 4, while if the diameter is 10 μm, Pe goes up to 104. As can be seen, the transition size between diffusive and drifting behavior for individual particles with density different from that of the medium occurs around 1 μm.
The physical system was simulated on a simple cubic lattice with periodic boundary conditions on the three spatial directions and with lattice spacing equal to d (=2a, the diameter of the particles). Each of the cells of the lattice can be occupied by a colloidal particle or can be empty, in which case is considered to be occupied by the suspension fluid. As the aggregation proceeds, we deal with a collection of clusters made of nearest neighbor lattice cells that are diffusing randomly and sedimenting downward. The Monte Carlo (MC) step time was picked in order to avoid a cluster to move farther than one lattice spacing. Define Δtdif≡d2/6Dmax as the time taken by the most mobile cluster to diffuse one lattice spacing. Also, Let be the time taken by the largest cluster to sediment downward one lattice spacing. The algorithm used is as follows.
If Δtdrif < Δtdif, then
pick a cluster in a cyclic way. That is, the clusters are numbered, and after having chosen a cluster, the next cluster to choose will be the following in the list, given that sedimentation is a deterministic process.
The time is increased by Δtdrif/Nc(t) where Nc(t) is the number of clusters in the sample at time t. This comes from the fact that after a cycle, when all the clusters are picked once (assuming that no merging took place), the time increase is Δtdrif.
The quantity Δtdrifvs/d is calculated, and the result is added to a real variable associated to that cluster, that is, called Zsed. If the sum is greater than one, the cluster is moved downward one lattice spacing, and the new value of the variable becomes the remainder of the sum modulo one.
If the cluster is moved, we check for overlapping with other clusters, in which case the moved cluster is taken back to the original position and the overlapping clusters are merged. Afterward, we go back to the starting situation to calculate Δtdif and Δtdrif.
The cluster now moves one lattice spacing on a random direction, with probability .
We then check for overlapping with other clusters, as in point (4).
If, on the other hand, Δtdif < Δtdrif we follow exactly the same procedure (a), with the exception that everywhere we find Δtdrif in a formula, it is replaced by Δtdif.
The simulations were stopped just before gelation. Only one volume fraction was considered, at the value φ = 0.01, for which 50 simulations of 156,250 particles were done for each of the following Peclet numbers: Pe = 10− 1 (the drift is felt predominantly by clusters of size ~10 or larger), 10− 3 (the drift is felt by clusters of ~1000 or larger), and 10− 8 (DLCA simulation, no drift). The results from those simulations were quite encouraging. For the structure of the clusters, the radius of gyration of each of the formed clusters was plotted vs. the cluster size as a point in a log-log graph. The inverse of the slope would proportion the fractal dimension of the formed clusters. In Figure 4, we are showing this plot for the case Pe = 10− 3. It is seen how the plotted points lie on not one but two straight lines, with a breaking point around the cluster size s = 1000. Two straight lines were fitted, one for the sizes below 1000 and another for the sizes above this value. For the small clusters, a fractal dimension df1 ≅ 1.87 was obtained, while for the large ones, it was found df2 ≅ 2.27. The df1 is in accord with the results for the df in the DLCA case for this concentration [31, 32], while the df2 coincides with the value obtained by Allain et al. [37, 38]. For the other two cases, only one straight line was seen. The first, for Pe = 10− 8, proportioned the value df ≅ 1.89, which again is in good accord with the known results for pure DLCA [31, 32], while the second is in very good agreement with the values obtained by Allain et al. The kinetic behavior was also studied, and it was seen that the weight-average cluster size presented the well-known behavior Sw(t) ~ tz, for the DLCA case (Pe = 10−8) only. For the other two Peclet numbers studied, the curve for Sw(t) departed from the straight line in the log-log plots, this departure being more pronounced in the Pe = 0.1 case, looking more like an exponential growth rather than a power law of the time. In other words, there was a sharp speeding up of the aggregation rate. Finally, when plotting the cluster size distribution Ns(t) as a function of s for different times during the aggregation process, it was seen a widening of this curve particularly at the later times becoming algebraically decaying, indicating the coexistence of very large clusters with smaller ones. Those three effects—the increase in the fractal dimension of the large clusters, the sharp increase in the aggregation kinetics, and the widening of the cluster size distribution—were explained by the same mechanism: there is an apparent instability in which the big clusters start to grow bigger by sweeping the smaller ones on their fast settling rate downward, increasing in this way their growth rate, capturing the smaller ones that fit in their holes and cavities (making the bigger ones more compact, with a higher fractal dimension), and forcing the very large clusters to coexist with smaller ones, not yet being trapped by the larger aggregates. It should be pointed out that it is not necessary to invoke the restructuring of the clusters—given that this is a rigid aggregation model—in order to achieve the higher values of the fractal dimension but only the sweeping of the small clusters and particles by the large drifting clusters, getting into the holes of the large ones which forces them to compactify.
One of the shortcomings of the model in Ref.  was that it was not possible to obtain the much higher initial volume fraction required to have gelation. This was due to the fact that there was no cluster deposition built-in in the model. We should point out that in the model, if a large cluster is rapidly settling downward, when it reaches the bottom of the simulation cell, it enters to the cell on the top of it, due to the periodic boundary conditions, but never deposits at the bottom of a real cell. To remedy this situation, an improvement of the model was made in Ref. . To include the height H of any prismatic cell, a cubic lattice was considered with periodic boundary conditions on the three spatial directions, which represented in an average way any cubic slice of the prism of height H (≅20L in Ref. ). In this sense, the model is a sort of mean field simulation. Whenever a cluster moves on the z direction a lattice spacing d, a quantity Zs belonging to each cluster is updated that proportions the total distance in units of d which the cluster has moved downward (taking into account the movement of all its ancestors, i.e., the largest of the two quantities Zs of the “parent“ clusters is inherited to the “child” cluster after a merging). If the movement is downward, the cluster is taken out of the cubic box (which means that it has been deposited on the bottom) with probability 1/(H − Zs), because of all the possible H − Zs rows where the center of the cluster can be, only at lowest row the cluster can leave the sample by a sedimentation step. Thus, the bigger clusters with the larger Zs are taken out of the sample with higher probability. As in the previous case, the simulations were stopped just before gelation for the gelling systems; for the nongelling simulations, they were stopped after all the clusters had been deposited. One volume fraction φ = 0.01 was considered, for which five values of Pe were chosen; Pe = 10− 8 (DLCA simulation), 10− 3, 10− 2, 10− 1, and 100. For the first three Peclet numbers, the system was gelling, while for the last Peclet number, it was nongelling; for Pe = 0.1, some of the simulations percolated the box, while the others did not. The size of the cubic box was taken as L = 270—which implies 196,830 particles—and H = 10 L. Another improvement made in Ref.  was in calculating the df. In the previous work , every cluster formed during the aggregation process was plotted as a point in the log-log plots of Rg vs. size. It was difficult to see where the breaking points of the straight lines (if any) occurred in the graph. Now it was plotted the average radius of gyration vs. the number of particles in the cluster, where the average was performed over all clusters formed lying on segments of constant size on the logarithmic size scale. In Figure 5, on the left panel, it is shown this plot for the case Pe = 0.01. We now clearly see a breaking point at around the size of about 200. On the left and right of the breaking point, there are well-defined straight lines between the arrows labeled “a” and “b” that give the two fractal dimensions: dfa = 1.868 for the small clusters and dfb = 2.270 for the large ones. Again, the dfa is in accord with the value for DLCA at that concentration, while the dfb is close to the Allain et al. value (2.2–2.3). In the right panel of Figure 5, we now see the average of Sw over all the simulations made for each Peclet number, except for Pe = 0.1 for which the crosses correspond to the average over the gelling simulations, while the triangles are for the nongelling simulations. As we can see, for the gelling simulations, we again find that the Sw separates from the DLCA simulation, the highest the Peclet number, the sooner the curve starts to take off. For the nongelling simulations, we now see a very interesting phenomenon: After the speeding up of the Sw, it reaches a maximum, and the average cluster size diminishes all the way to zero. This indicates that the dynamics of the aggregation accelerates and then, afterward, it slows down. This is due to the depletion of mass in the bulk as the clusters deposit on the bottom. In the same Ref.  and with more details in Ref. , it was estimated the critical concentration to have gelation for the different Peclet numbers studied, which was around the value found by Allain et al. Finally, in Ref. , the scaling of the cluster size distribution was studied, and it was shown that, for very low Peclet numbers, the cluster size distribution was asymmetrically bell-shaped as in DLCA. As the Peclet number increased, it was seen how that distribution was transforming into an algebraically decaying distribution (with an exponent τ ≅ 1.73), implying the coexistence of small and large clusters.
Subsequently Leone et al.  performed an off-lattice simulation with 10,000 particles inside a square prism with periodic boundary conditions on the x and y directions only, but with no periodicity on the vertical z direction, in an effort to consider the dependence of the aggregating quantities on the vertical z direction. The mechanisms for moving the particles were similar as those in Ref. . However, this z dependence was obtained only for three heights: the authors divided the prismatic cell in three different equally sized regions—top, middle, and bottom. The aggregating quantities were obtained by calculating their averages separately in each of those regions. Unfortunately, due to the small number of particles and the very low concentration considered, they were not able to see very clearly the breaking point on the Rg vs. size plot that differentiates the two cluster fractal dimensions. For example, for the top region they found for DLCA while df ≅ 1.88 for Pe = 0.1, while for the middle region, they found df ≅ 1.74 in the DLCA case, while for Pe = 0.1, it was obtained df ≅ 1.89. As in Ref. , they found the slowing down of the aggregation kinetics due to the depletion of mass in the bulk, but with no speeding up before, which is probably due to the very low concentration studied which prevented the sweeping mechanism to act thoroughly. Subsequently, with this off-lattice model, the case of a sticking probability lower than one was investigated . It was found as expected that the fractal dimension increased to values close to the Allain et al. values of 2.2, but as a result of lowering the sticking probability, not actually due to the sedimentation effect. Finally, they also studied the nature of the deposited clusters and their fractal dimension, distinguishing two cases: for low Peclet numbers, the cluster reactions were occurring mainly in the bulk before deposition; hence, their fractal dimension was higher than that for higher Peclet numbers because, in this last case, the deposited clusters were small and their reactions were mainly occurring on the base of the prism (assuming no sticking of the clusters and particles with the base of the prism), leading to a two-dimensional aggregation that is known to produce fractal clusters with df ≅ 1.45.
The low values of the df for the sedimenting clusters for this very low value of the volume fraction (φ ≅ 6.7 × 10−5), as compared to the high values obtained in Refs. [43, 44], prompted us to consider the concentration effects for the problem . The two computer models, Refs. [43–45], were used to include the effects of concentration, spanning three orders of magnitude in volume fraction. Recall that for the high concentrated case [43, 44], there was a crossover from a DLCA type of regime to another regime, for larger clusters, in which the fractal dimension took a value close to the Allain et al. value of 2.2–2.3 and a speeding up followed by a slowing down of the aggregation rate that we called the sweeping scaling regime. As the concentration was lowered from this high value of φ = 0.01 all the way up to φ = 0.0001, it was found that the crossover to the sweeping regime remained as long as the height of the cell was drastically increased. For example, for φ = 0.01 and H = 2700, there was a crossover to the sweeping scaling regime for almost all the Peclet numbers studied (except for the DLCA case, of course), with a high value of the of the fractal dimension—of the order of the Allain et al. value of 2.2–2.3—as in Refs. [43, 44]. Similarly, there was a speeding up of the aggregation rate followed by the slowing down for the nongelling simulations. For the φ = 0.001 case, there was also a crossover to a sweeping scaling regime for the two heights studied, which were H = 4500 and 45,000. Similarly, the speeding up and, afterward, the slowing down of the aggregation rate were found. However, in this case it was found for both heights that the fractal dimension for the sweeping scaling regime did not increase as high as the Allain et al. value, but it stayed around df ≈ 2.05. For the φ = 0.0001 and H =7000 case, it was not possible to find a sweeping scaling regime for the large clusters. In other words, in the log-log plots of Rg vs. size, only one straight line was found that proportioned the value df ≈ 1.83, that is, close to the usual DLCA value. However, by increasing drastically the height (and the computer time) up to H = 70,000, the sweeping scaling regime was recovered again, with a df ≈ 2.05 as in the previous φ = 0.001 case. This corroborates the assertion that it is possible to find a sweeping scaling regime for low concentrations if the height of the cell is highly increased. Nevertheless, it was intriguing to find that the fractal dimension for the sweeping scaling regime decreased to the value of df ≈ 2.05 for the two low concentrations studied. This could indicate that there may be a universal value for this fractal dimension, as long as the concentration is not too high. Concerning about the clusters deposited at the bottom of the prism, the new regime besides those found from the bulk was confirmed. This is the restricted quasi-two-dimensional aggregation of those clusters that have been deposited on the bottom of the container (assuming again no sticking of the clusters to the base of the prism). If the sediment clusters have been grown essentially in the bulk, which occurs basically for the more concentrated systems, one or the two fractal dimensions can be found in the sediments. However, if the system is dilute enough and/or for a high Peclet number, we may favor the arrival of not very large clusters that can aggregate in a quasi-two-dimensional way.
2.3. Cluster restructuring vs. cavity occlusion to compact a cluster
So far, so good. Maybe not as good, since now it is necessary to explain the lower fractal dimension obtained in the simulation (df ≈ 2.05) for the sweeping scaling regime at the fairly low concentrations used in the experiments, as compared to the Allain et al. value (2.2–2.3). Besides that, in Ref. , the authors presented a two-dimensional simulation on a rectangular lattice with periodical boundary conditions, that is, they considered the z coordinate on the direction of the gravitational field and only one perpendicular direction. They show that for aggregation combined with sedimentation, their clusters did not present spatial scale invariance, that is, they are not self-similar fractal objects. This came from the fact that the mean width perpendicular to the field direction scaled as a power law with the cluster size, while the mean height did not. The method they used to treat the mechanisms for diffusion and sedimentation of clusters differed from the previous simulations [42–47], in that they were using three probabilities for moving a cluster, p↓ ↓ for a cluster moving downward, p↑ ↑ for that cluster moving upward, and p⊥ for the cluster to move on the perpendicular direction. However, these probabilities were calculated from the Stokes-Einstein formalism, and this should not present an essential problem. The only concern is about the deterministic nature of sedimentation, but it may not be very relevant. Another interesting result found in that paper was that their clusters were initially elongated along the field direction, while as time proceeds, the situation reversed and the clusters became oblate. Although this lack of spatial scaling was somewhat annoying, the decisive event that led us to consider a more realistic simulation, not the mean field one that we have shown, was the observation of a figure in Ref. , which is presented below in Figure 6.
In Figure 6, where the horizontal lines were drawn by the present author, we can see that the behavior of the growing and settling clusters is different in each of the layers, labeled with the numbers 1 through 7. For the top layers, we can see that the clusters grow up to a certain maximum size, and then the layer starts deplenishing of colloidal matter, which forces the mean cluster size to become smaller and smaller. For the bottom layers, a similar situation happens with the difference that this maximum size occurs at later times and, additionally, is much larger. In other words, the aggregation of colloids coupled with sedimentation is stratified, in the sense that the aggregating quantities behave differently from layer to layer. Our mean field simulation, although it gives results similar to the experimental ones, does not consider this stratification, since the simulation box represented in an average way any cubic slice of the aggregating prism of height H. Maybe by considering the precise, proper algorithm, it would be possible to give an answer to this lack of spatial scaling of the clusters. Furthermore, we also may be able to find a fractal dimension in the sweeping regime similar to the experimentalist value. A quick estimate to count the number of lattice cells that would be required to consider a square prism similar in volume (in units of particle volumes) to the one used by the experimentalists proportioned the number of lattice cells around 1018 (we recall that the long tube used in Ref.  was 4.4 cm in diameter, while the tube length was 80 cm and the particle diameter was 0.07 μm). To perform a simulation of that size is an impossible task, even with present-day computers. To circumvent the problem, the search for a simulation that could proportion the behavior inside a slice of arbitrary vertical thickness Δ and centered at any given depth Z (in lattice spacings) of the aggregation tube was started. In Refs. [49, 50], such simulation was presented, quite similar to the one described in Refs. [43, 44]. We considered a cubic lattice of lattice spacing d (the diameter of the particles) and lateral size L = 650 d, with periodic boundary conditions on the x, y, and z directions that represent a layer centered at the integer depth Z (in lattice spacings), measured from the top of a lattice prism, and of thickness ∆ (which is taken multiple of two). This differs from the previous algorithm, in that the cubic lattice there was representing the whole prism of height H. As before, whenever a cluster moves in the z direction a lattice spacing, an integer quantity Zs is updated that proportions the total distance that the cluster has moved downward. Again, every time a cluster is picked at each Monte Carlo step of time ∆t, we check whether it diffuses and/or sediments downward a lattice spacing. As we want to simulate the clusters inside the abovementioned layer of the prism, of all the Z + Δ/2 − Zs ≥ 1 possible rows of cells where the center of this cluster can be, only at the lowest row—with the position Z + Δ/2—will be taking out of the layer by a sedimentation step, while if it is at the row at Z + Δ/2, the sedimentation step will take it inside the layer. Therefore, if Zs has not grown enough for such cluster such that Z − Δ/2 − Zs ≥ 1, the cluster is taken out of the box with probability zero (not taken out), because the cluster has the same probability to be on either of those two rows. On the other hand, if Z − Δ/2 − Zs ≤ 0, the cluster cannot be at the row Z − Δ/2 and is therefore taken out of the box with probability 1/(Z + Δ/2 − Zs). This is the manner used to simulate the depletion of colloidal matter inside the different layers in the aggregating prism, due to sedimentation. Finally, as before, at each diffusion or sedimentation movement that does not take the cluster out of the box, we check for overlapping with other clusters. If that happens, the cluster is taken back to its original position and the overlapping clusters are merged. The volume fraction was fixed at the value φ = 0.001. Four depths were considered at Z = 500, 5000, 50,000, and 500,000, and one single ∆ = 100. For each of those depths, six Peclet numbers were considered: Pe = 0.0001, 0.001, 0.01, 0.1, 1.0, and 0.0 (DLCA simulation). In turn, for each Peclet number, a series of 10 simulations of 274,625 particles were made, in order to have enough statistics to evaluate the structural and dynamical quantities. To investigate the lack of spatial self-similarity of the formed clusters, not only the log-log plots of the radius of gyration vs. size were made but also of the x, y, and z components of the radius of gyration tensor, where , etc. In addition, the following anisotropy measures were defined: Axz ≡ 〈Rgx/Rg〉/〈Rgz/Rg〉 and Ayz ≡ 〈Rgy/Rg〉/〈Rgz/Rg〉, which were also calculated as a function of size. By analyzing these quantities , it was found a whole variety of behaviors for the structure of the clusters in this problem, depending on the sedimentation strength (Pe), the layer depth (Z), and the region of sizes considered. Generally speaking the four radii (Rgx, Rgy, Rgz and Rg) scale as a power law with N, with the same scaling power, for the small, non-settling clusters, except for the cases with a high sedimentation strengths: Pe = 1.0 and 0.1 for all depths. This was called the quasi-DLCA regime, where the “quasi” means that, as there is still some sweeping of even smaller clusters, the fractal dimension is a little bit higher than the usual DLCA df. For the large, settling clusters, it was found cases for which the four radii scale again as a power law with N, with the same scaling power, making it possible to define a settling-cluster fractal dimension, a behavior that was called the sweeping scaling regime. There are, however, cases for which the scaling powers for the horizontal and vertical directions differ, obtaining therefore self-affine settling clusters. There are still some cases for which only the mean width or the mean height scale as a power law with N and even further cases for which no scaling as a power law of the four radii is obtained. With respect to the anisotropy measures as a function of the size, it was generally found that for the small clusters, the measure remained either close to one or diminishes somewhat from this value, making the clusters prolate along the field direction. However, as the clusters grow in size and start sweeping smaller clusters, they become oblate and in some cases quite oblate, reaching values for the anisotropy measures of up to 3. For more details about this, the reader is referred to Ref. . With respect to the fractal dimension for the sweeping scaling regime for large clusters (after the quasi-DLCA regime for the smaller ones) which generally occurs for low Peclet numbers and big depths, it was found as in Ref.  that the df was generally around 2.05. This reinforces the idea that this may be a universal value for perfectly rigid clusters. However, this needs to be tested with more simulations for many different cases. It should be mentioned that we did not allow the clusters to rotate as in the experimental systems which, in principle, may be able to do so, canceling in this way the anisotropies found (prolate or oblate clusters). Additional anisotropies besides the prolate or oblate shape may be obtained in this problem. A moment of reflection shows that, even for intermediate drift strengths, the clusters are more compact on the bottom than at the top. This is due to the fact that the bottom is the part of the cluster that is catching up the small clusters and single particles, becoming in this way more compact. It may be conceivable that the new anisotropy may hamper the rotation of the clusters, because a higher gravitational force per unit volume would be experienced at the bottom than at the top, while the medium frictional resistance would be the same everywhere around the cluster, under a non-draining assumption. Hence, the possibility exists to find the prolate or oblate shape of the experimental clusters after all, in the case that they may not be allowed to rotate; this, however, needs to be tested experimentally.
Concerning the kinetics of the aggregation, the speeding up with respect to the DLCA case, followed by the slowing down due to the depletion of mass in the layer, was confirmed. Perhaps more interesting, the experimental behavior found in Wafra´s thesis was obtained  that the aggregates in a layer reach a maximum size at a certain time, the deeper the layer is, the longer it takes to reach this maximum which, however, is of a much larger size. In Figure 7, all these kinetic behaviors are shown. Note that in those cases for which no decrease in the Sw is shown, particularly for low Peclet numbers and big depths, indicates that the system gelled at the point where the curve was cut.
As mentioned, the low value of the fractal dimension (~2.05) for the sweeping scaling regime, as compared to the Allain et al. value (2.2–2.3), still needs to be explained. There is much evidence, both theoretical [22, 25, 51, 52] and experimental [21, 23, 24, 53], that a tenuous fractal object in a shear flow would experience hydrodynamic stresses that can change its structure. At the beginning of the growth process, the aggregates are too small to be sensitive to those stresses. Nonetheless, as soon as they grow beyond some limiting size, they may undergo structural rearrangements due to the stresses felt when they sediment faster. It has been shown that the total force exerted by the fluid on the cluster is concentrated into the most peripheral particles . The cluster-cluster aggregates are loopless and have an elastic behavior similar to that of a contorted chain. The elasticity is dominated by the bending elasticity which, when a pair of forces applied at the tips of the cluster bends the corresponding branches, makes the branches touch and stick. Because of this, the cluster could become more compact, with a higher fractal dimension. This is the explanation advanced by Allain et al.  for obtaining a high fractal dimension for the settling clusters. However, there is also much simulational evidence for obtaining a high fractal dimension of large rigid clusters in the sweeping scaling regime, coming from the occlusion by the small clusters and particles that enter into the holes and cavities of those large clusters. The correct explanation is probably a combination of both effects, which would make the large settling clusters to push the value of the apparent fractal dimension (“apparent” for the cases in which the scaling behavior is not the same along the horizontal and vertical directions) to the order of 2.2–2.3. Another factor that may increase the fractal dimension of settling clusters is the inclusion of a sticking probability less than one. As it is well known, without a gravitational field acting on the system, the df increases from its DLCA value of ~1.8 to an RLCA value of ~2.1. With the inclusion of sedimentation, this could increase the value to around 2.2–2.3. However, as stated by Allain et al., they believe their system belongs to a regime with a sticking probability equal to one. There is also the influence of the hydrodynamic interactions between clusters, whose effect on the df is unclear. The general belief, nonetheless, is that by working with sufficiently dilute systems (φ ≤ 0.001), their influence would be negligible. Needless to say that more research on this interesting subject would be worthwhile to explore, due to its practical application in several aspects of science.
3. Some other relevant simulational results
Hawick [54, 55] performed some simulations in 2D and 3D with the major emphasis in the visualization of the clusters formed during their aggregation and sedimentation. As in  he used rectangular in 2D or prismatic boxes in 3D, with their longest length along the field, vertical direction. Similar to , periodic boundary conditions were considered on the direction(s) perpendicular to the field, but not for the vertical direction with fixed vertical boundaries. For the movement of the cluster, the author used the probabilistic mechanisms proposed in Ref. . In Ref. , for the 2D clusters, he found that the large clusters have an inverted V or W shape or more generally an inverted tree-like structure. However, for the 3D clusters, this was not the case, and they were not so anisotropic. One of the shortcomings of this work is the smallness of the systems used. For example, in Ref.  for the 3D case, a box of size 128 × 642 was used. One of the interesting results shown and visualized was the elongated shape along the field direction of the cluster for a Peclet number of 0.01, while by changing this number to 0.1, wider clusters were found, which is in line with the results presented here. In another work, Whitmer and Luijten  consider a molecular dynamics simulation for the sedimentation of diffusing and aggregating colloidal clusters. The potential of interaction between the colloids was taken as the Asakura-Oosawa-Vrij model [57, 58] of depletion-induced attraction, while the Brownian motion was controlled via a Langevin thermostat. They were mainly interested in the structure of the sediment at the bottom of the container and noticed that, for low attractions incapable of inducing cluster formation during sedimentation, this structure is similar to that expected for equilibrium colloidal suspensions. Modest attraction strengths, however, are found to enhance the growth of crystalline phases in the sediment. Furthermore, high attraction strengths induce the formation of colloidal clusters during sedimentation, which settled into loosely packed disordered structures.
4. Discussion and future work
One of the criticisms that can be made to the abovementioned simulational models is that do not proportion quantities that are usually measured by experimentalists. For example, the absence of a cluster fractal dimension due to the anisotropies found is not currently available to the experimentalists. In several experimental works, researchers are mainly focused on sedimentation rates and sedimentation halftimes, something that they can measure in their aggregation vials . The structure of the sediment, as studied in Ref. , is another possible topic that merits a further consideration both by computer simulators and experimentalists. The sedimentation halftimes, for example, can be studied computationally by measuring the movement of the center of mass of all the colloidal particles in the clusters, as they settle downward. As it was stated above, in this review, we were mainly focused on the sedimentation and deposition of single clusters, that is, for concentrations lower than the volume fraction φ* define by Allain et al. . We should not be afraid to go beyond the φ* and try to simulate the formation and the collective settling of a gelled suspension [39–41, 60]. For this purpose, perhaps we would need to do an off-lattice simulation with non-perfectly rigid bonds between the colloidal particles, but that can have some mild tolerance for bending and also for breaking under a high stress. As we can see, the field is wide open to exciting, future research.
The author is grateful to the Academic Supercomputing Committee at the National Autonomous University of Mexico for the computational resources assigned through DGTIC-UNAM (project SC16-1-IR-52).