Secondary Nucleation by Interparticle Energies. III. Nucleation Rate Model.

A nucleation rate model for describing the kinetics of secondary nucleation caused by interparticle energies (SNIPEs) is derived theoretically, verified numerically, and validated experimentally. The theoretical derivation reveals that the SNIPE mechanism can be viewed as enhanced primary nucleation, i.e., primary nucleation with a lower thermodynamic energy barrier (for nucleation) and a smaller critical nucleus size, both caused by the interparticle interactions and the associated energy between the surface of a seed crystal and a molecular cluster in solution, as shown in part I of this series. In the case of a sufficiently agitated suspension, the model depends on four parameters: two reflecting primary nucleation kinetics and the other two accounting for the intensity and effective spatial range of the interparticle interactions. As a numerical verification of the model, we show that the nucleation kinetics described by the SNIPE rate model is in quantitative agreement with those given by the kinetic rate equation model developed in part II of this series. A sensitivity analysis of the SNIPE rate model is conducted to present the effect of key model parameters on the nucleation kinetics. Moreover, the SNIPE rate model is validated by fitting the model to the time-resolved data of secondary nucleation experiments as well as to two other, well-known secondary nucleation rate models. Importantly, all of the estimated parameter values for the SNIPE model were consistent with the theoretical estimates, while some of the estimated parameter values for one of the well-known secondary nucleation models deviated from the corresponding theoretical values significantly.

Introduction

Secondary nucleation is the dominant nucleation mechanism in industrial crystallizers,[1] particularly when operated continuously.[2] Despite its importance, the underlying mechanisms of secondary nucleation are still debated.[1−3] Secondary nucleation has been largely attributed to the mechanical attrition of a seed crystal,[1−3] caused by its collision with a stirrer blade, the inner walls of a crystallizer, or another crystal,[4,5] thus generating a secondary nucleus in the form of an attrition fragment whose crystal structure is identical to that of the seed. However, secondary nucleation can occur even without any collision,[6−8] and the resulting secondary nucleus can have a polymorphic/chiral crystal structure different from that of a seed,[9−15] hence the need for other perspectives to complement the attrition mechanism.

A group of other mechanisms explains secondary nucleation without including the attrition of a seed in their description, and some of them even rationalize the formation of a secondary nucleus having a different polymorphic/chiral structure.[16−18] These mechanisms consider secondary nucleation as nucleation induced by the surface of a seed crystal through various phenomena, e.g., a direct contact between the seed surface and a molecular cluster in solution,[19] the detachment of a two-dimensional (2D) nucleus (i.e., a growth unit[20]) from the seed surface,[4] a rapid coagulation of molecular clusters facilitated by the van der Waals forces between the seed surface and the clusters,[16] and the energetic stabilization of a molecular cluster due to its interparticle interactions with the seed surface.[17,18] The last mechanism, called secondary nucleation by interparticle energies (SNIPEs), is the focus of this series.

In part I,[17] the thermodynamics of the SNIPE mechanism has been explained, which reveals that the interparticle interactions between a molecular cluster and a seed surface can lower the thermodynamic energy barrier for nucleation. Moreover, this energetic stabilization effect has been quantified through two parameters, namely, the intensity of the stabilization effect, Est, and the effective range thereof, lst. In part II,[18] the kinetics of the SNIPE mechanism has been demonstrated, which shows that the stabilization effect can enhance nucleation kinetics to the extent that nucleation occurs even when the supersaturation is at such a low level that primary nucleation cannot occur. To address this kinetic aspect, a kinetic rate equation (KRE) model has been developed by including the stabilization effect parameters (i.e., Est and lst) in the conventional KRE model of nucleation (i.e., the Szilard model[21,22]). As part III of the series, this work aims at transforming the developed KRE into a nucleation rate model, which can be used in a population balance equation (PBE) model to describe crystallization processes and to compare its description with experimental data quantitatively. In doing this, we will account also for the effect of mixing in the crystallizer where secondary nucleation may take place. This is something that we did not touch upon in parts I and II, where we assumed that each seed was exposed to the same solution in which they were suspended without any hindrance caused by the presence of the other seeds; this was tantamount to assuming perfect mixing. Therefore, this work allows us both to assess the scientific validity of the SNIPE mechanism and to show its practical applicability.

This work aims at achieving three objectives. The first is to derive the SNIPE-based nucleation rate model (called SNIPE rate model for brevity). The second is to verify the SNIPE rate model by showing that its output matches the nucleation kinetics described by the KRE model developed in part II[18] and by conducting a sensitivity analysis of the SNIPE rate model. The third objective is to validate the SNIPE rate model by using it for fitting experimental data and to assess its goodness-of-fit and theoretical consistency in comparison with those of other secondary nucleation models.

This contribution is structured as follows. Section explains the experimental data to be fitted, the methods for modeling nucleation and growth, the derivation of the SNIPE rate model, and the procedure for fitting. In Section , the verification and sensitivity analysis of the SNIPE rate model are shown, followed by the experimental validation of the model and its comparison with other secondary nucleation models.

Methods

Experimental data (see Section ) were used for fitting the output of a population balance equation (see Section ) that describes growth and nucleation in a well-mixed batch reactor, where different nucleation rate models (see Section ) were considered. Three secondary nucleation rate models (see Section ) were employed, including the SNIPE rate model, whose expressions were derived by adopting the approach used in classical nucleation theory (CNT) (see Section ). The parameters of each model were estimated by applying the method of least squares (see Section ).

Data: Benchmark Experimental Study

Fitting a secondary nucleation model to experimental data can be challenging in many ways, including due to the unknown kinetics of concomitant crystallization mechanisms (e.g., primary nucleation and growth), a change in operating temperature and thus temperature-dependent physical parameters (e.g., the specific surface energy[23−25]), and the use of an unconventional crystallizer. To minimize the number of such factors, it is crucial to use a set of experimental data that has been acquired using a widely studied system in the simplest settings, e.g., isothermal seeded batch crystallization of a compound whose primary nucleation and growth kinetics have already been characterized by the literature.

Among experimental investigations on secondary nucleation accompanied by online or ex situ monitoring techniques,[26−31] many were not used in this study due to a different combination of the following reasons: change in operating temperature,[27−29] undetermined seed size distributions,[27,28] use of an unconventional crystallizer,[29] and polymorph formation.[30,31] Finally, ref (26) was chosen as a benchmark study; it concerns isothermal seeded batch crystallization of paracetamol from a 500 mL ethanol solution, with seed size distributions and growth kinetics predetermined by the same researchers in another study[32] as well as with a good understanding of primary nucleation kinetics.[33−35] All reported data measured at 20 °C and 200 rpm were used for estimating the secondary nucleation model parameters (see their main characteristics in Table ). As explained in Section S1 in the Supporting Information, the data of two experiments at 250 and 300 rpm were not used because the effect of stirring intensity on the secondary nucleation rate in those experiments is not consistent with the existing understanding of secondary nucleation.

Table 1

List of Experiments[26] Used for Fitting Secondary Nucleation Modelsa

exp	label in ref (26)	S0	M0seed [g]	sieve size fraction [μm]
E1	S03	1.57	1	120–250
E2	S02	1.42	1	120–250
E3	S05	1.42	3	120–250
E4	S01	1.42	7	120–250
E5	S06	1.42	1	90–125

S0 is the initial bulk supersaturation, and M0seed is the initial mass of seeds.

Population Balance Modeling

The growth and nucleation of a crystal population can be simulated using a population balance equation (PBE) coupled with the solute mass balance. For a well-mixed batch reactor, the PBE can be written aswhere t is the time, L is the characteristic size of the crystal, f(t, L) is the number density function (called PSD), i.e., f(t, L)dL being the number of crystals with length L ∈[L;L + dL] per unit suspension volume,[36]G is the crystal growth rate assumed to be size-independent, and J is the nucleation rate. Neither agglomeration nor breakage is considered in the model. Equation represents the initial condition, with f0(L) being the initial PSD, while eq indicates the boundary condition. The PBE is coupled with the mass balance, which can be written aswhere c(t) is the bulk solute concentration, m = ∫0∞Lf(t, L)dL (I = 0,1,2,···) is the ith moment of the PSD, kv is the volume shape factor, V1 is the molecular volume, κ = 106 m μm–1 is a numerical factor to convert micrometer to meter (this is necessary to be consistent with the units chosen for the different variables, as reported in the Nomenclature at the end of this paper), and c0 is the initial solute concentration.

The PBE was solved numerically using a fully discrete, high-resolution finite volume method[37] with the van Leer Flux limiter[38] and by satisfying the convergence condition by Courant–Friedrichs–Lewy.[39] At the upper bound of the size domain, a numerical outflow boundary condition was applied by employing the zero-order extrapolation method.[37]

Growth Rate Model

The empirical growth rate model used in the benchmark study[26] was employed in this work while using the same parameter valueswhere θG with i = 1, 2, 3 are the growth kinetic parameters, R is the gas constant, T is the temperature, and Δc = c–c is the absolute supersaturation. When Δc is in units of kmol m–3, the values of the growth kinetic parameters are θ1G = 9.979 × 106 μm s–1 (m3 kmol–1)θ, θ2G = 4.056 × 104 J mol-1, and θ3G = 1.602.

Nucleation Rate Models

An overview of a primary nucleation rate model is presented as background (see Section ), followed by a derivation of the SNIPE rate model as well as of the kinetic models of the other secondary nucleation mechanisms (see Section ). In the following, θX denotes the ith parameter of a model identified by the label X.

Background: Primary Nucleation

According to CNT,[22] the primary nucleation rate can be described by the following expression

where θPN (i = 1,2) are the parameters to be estimated and s is the monomer supersaturation defined as[23,40]where Z1 and C1,e are the monomer concentration and the monomer solubility, respectively. In this work, to account for the presence of molecular clusters in solution and its impact on crystallization kinetics, the monomer supersaturation s has been used instead of the bulk supersaturation S = c/ce, where ce is the bulk solubility.[40] For the sake of completeness, a relationship between s and S as well as that between C1,e and ce is reported in Section S2 in the Supporting Information (see also ref (40)).

Alternatively, eq can be written in a form that reveals the physical mechanism of CNT[22,41]where B(n) is the frequency of monomer attachment to an n-sized cluster (i.e., a molecular cluster consisting of n monomers), z is the Zeldovich factor, C is the equilibrium concentration of the n-sized clusters, and n* is the size of the critical clusters, i.e., the critical (nucleus) size. Considering that the rate of monomer attachment is often limited by the surface-integration step,[23,42−44] the monomer attachment frequency B(n) is given by[22]where k0 is a lumped coefficient to be estimated from experimental data. This coefficient reflects the mass transport of solute molecules to the surface of an n-sized cluster, and its value is independent of mixing conditions. The Zeldovich factor z is defined as[22]Here, F = ΔG/kBT is the dimensionless Gibbs free energy for forming an n-sized cluster, ΔG is the dimensional Gibbs free energy, and kB is the Boltzmann constant. The equilibrium cluster concentration C is defined as[22]where C0 = ce/∑∞m exp(−Ωm2/3) is the concentration of nucleation sites.[40] The critical size n* can be calculated fromwhere, under the capillary approximation,[22] the Gibbs free energy F is given bywhere Ω = bγ/kBT is the dimensionless surface energy. Here, b = ks(V1/kv)2/3 is the surface area of a molecule, γ is the specific surface energy of a cluster, and ks is the surface area shape factor. Note that substituting eq into eq yields an explicit expression for the critical size, namely, n* = (2Ω/3 ln s)3. Finally, a comparison between eq and the set consisting of eqs –14 results in the theoretical values of the two parameters in eq (θ1PN and θ2PN):Note that, typically, the lumped coefficient k0 and the dimensionless surface energy Ω are the unknown parameters to be estimated by fitting the nucleation kinetic model (eq ) to relevant experimental data (i.e., a set of primary nucleation rates characterized at different values of supersaturation).

Secondary Nucleation

The kinetics of secondary nucleation can be described using either empirical rate expressions or first-principle rate models, as illustrated in Figure . The empirical rate expressions can be fitted to experimental data adequately[26,28,30,31,45,46] but they rarely advance the fundamental understanding of secondary nucleation mechanisms. Alternatively, the first-principle models can be used to fit experimental data as well as to gain a deeper insight into the physical mechanisms of secondary nucleation. The first-principle models can be grouped into two categories: the models based on the attrition mechanisms[4,5,47,48] and those based on the surface-induced nucleation mechanisms, namely, the SNIPE mechanism,[17,18] the surface nucleation (SN) mechanism,[4] and the embryo coagulation secondary nucleation (ECSN) mechanism.[16]

Figure 1

Overview of the secondary nucleation rate models: empirical models[26,28,30,31,45,46] and first-principle models reflecting the attrition mechanism[4,5,47,48] or the surface-induced mechanisms.[4,16−18] The latter includes the SNIPE mechanism,[17,18] the surface nucleation (SN) mechanism,[4] and the embryo coagulation secondary nucleation (ECSN) mechanism.[16]

In this study, the attrition models were not used for fitting due to two main reasons. First, the researchers of the benchmark study indicate that the experimental data were fitted better to a model whereby a secondary nucleation rate is proportional to the total surface area of the crystal population rather than to its total volume,[26] thus suggesting that the surface-induced secondary nucleation mechanisms are likely more suitable for describing these data. Second, unlike other compounds (e.g., l-ascorbic acid,[29] potassium alum, and sodium chlorate[49]) and crystals with high aspect ratios,[50] the effect of attrition on paracetamol crystals (i.e., the model compound in the benchmark study) in a stirred tank appears to be insignificant,[51] likely due to the paracetamol crystals’ mechanical properties. Accordingly, the experimental data of the benchmark study[26] were fitted with all of the three surface-induced secondary nucleation models, whose expressions are presented in the following, including the relevant derivation for the SNIPE rate model.

SNIPE Model

The SNIPE model describes that secondary nucleation occurs because molecular clusters near the surface of seeds can become secondary nuclei thanks to the stabilization effect originating from the interparticle interactions and associated energy between the clusters and the seed surface.[17,18] As explained in part I and illustrated in Figure a, due to the short-range nature of the interparticle energies, the stabilization effect is localized in the stabilization volume (i.e., part of solution surrounding seed crystals) while being negligible in the bulk solution volume (i.e., part of the solution being far away from the seeds). To account for the localized stabilization effect and its influence on the kinetic of the SNIPE mechanism, three compartments and the corresponding volume fractions were introduced in part II,[18] as illustrated in Figure b: the solid volume fraction, ϕs, the stabilization volume fraction, ϕst, and the bulk solution volume fraction, ϕb, with ϕs + ϕst + ϕb = 1. The solid volume fraction ϕs and the stabilization volume fraction ϕst are given by[18]respectively, where lst (≥0) is the model parameter determining the thickness of the stabilization volume and m2seed is the second moment of the PSD of seeds, which measures their total surface area. Note that, in eqs and 18, the unit conversion factor κ is used for the same reason as in eq and that ϕst given by eq is an approximation of its exact geometric value, which is justified as discussed in part II.[18] In eq , it is assumed that very small crystals (e.g., nuclei) do not cause the stabilization effect and thus secondary nucleation;[18] a similar concept has also been applied in other secondary nucleation models.[5,52] Note that the bulk solution volume fraction ϕb can be calculated from the other two volume fractions as ϕb = 1 – ϕs – ϕst.

Figure 2

Modeling the stabilization effect of the SNIPE mechanism based on three compartments (a and b) and two compartments (c). The three compartments consist of three volume elements, namely, the solid volume, the stabilization volume, and the bulk solution volume, with their volume fractions denoted by ϕs, ϕst, and ϕb, respectively. Here, C is the equilibrium concentration of n-sized clusters and F is the Gibbs energy for forming an n-sized cluster, with their superscripts (e.g., “st” and “b”) indicating the corresponding volume. The variable with the superscript “eff” represents an effective quantity, whose value is an average of the same quantities in the stabilization and in the bulk volume.

Due to the localized stabilization effect, nucleation in the stabilization volume is thermodynamically more favorable than in the bulk solution. This can be described using a different Gibbs free energy for the formation of an n-sized cluster in the different compartments,[17,18] denoted by Fst and Fb for the stabilization volume and the bulk solution volume, respectively;[17,18] these are defined aswhere Est (≥1) is the model parameter determining the intensity of the stabilization effect. Replacing F in eq with Gibbs free energies in each volume (eq ) results in different equilibrium concentrations of the n-sized clusters in the stabilization volume and the bulk solution volume (denoted by Cst and Cb, respectively), thus yielding the effective equilibrium concentration of the n-sized clusters in the system, Ceff, as a volume-weighted average of Cst and Cbor with the help of eqs and 19where the effective Gibbs energy Feff, which is not a thermodynamic quantity, is defined as

Equations and 22 suggest that the nucleation process with the localized stabilization effect (illustrated in Figure b) can be described as a primary nucleation process, enhanced by the fact there is a lower energy barrier (i.e., Feff), as illustrated in Figure c. The corresponding nucleation rate, J∞, can be calculated using the framework of CNT (eqs –13), with F replaced by FeffHere, zeff is the effective Zeldovich factor calculated from eq and neff* is the effective critical size calculated numerically from eq , while F and n* are replaced by Feff (eq ) and neff*, respectively. Equation contains four parameters to be estimated, namely, two primary nucleation kinetic parameters (k0 and Ω) and two SNIPE parameters (Est and lst) that account for the stabilization effect. It is worth noting that eq constitutes a novel result, which allows using the SNIPE model within a population balance equation (eqs –3).

Effect of Agitation

It is worth noting that, in surface-induced secondary nucleation mechanisms, it is implicitly assumed that the entire surface area of all seed crystals is available for secondary nucleation and that secondary nuclei are swiftly removed by shear from the surroundings of seed crystals, as demonstrated in several experimental studies.[6−8,53−55]

In other words, the following two hydrodynamic conditions must be fulfilled for surface-induced secondary nucleation to be effective. The first condition is a good solid/liquid mixing, through which all seed crystals are well suspended and all of the surface areas of the seed crystals are available for secondary nucleation via interaction with newly formed clusters. The second condition is the presence of sufficient fluid shear around the seed crystals, as a means to shear off secondary nuclei from the stabilization volume around them, as demonstrated in various experimental studies.[6−8,53−55] In practice, as agitation weakens, less surface area of the seeds would be available for secondary nucleation and a fewer number of secondary nuclei would be released from the seeds, thus yielding a secondary nucleation rate smaller than J∞. The secondary nucleation rate becomes eventually negligible in the absence of agitation.

To account for the effect of agitation intensity on the secondary nucleation rate in the SNIPE model, the effective nucleation rate can be described as the product of the secondary nucleation rate in a well-mixed suspension, J∞, and an efficiency factor, ϵ(r), whose value depends on the stirring rate r

We introduce an empirical model for the efficiency factor ϵ(r) since developing a mechanistic model is beyond the scope of this work. The model is based on the physical intuition mentioned above; i.e., the effectiveness of the SNIPE mechanism is monotonically related to the degree of suspension of the seed crystals, while the latter depends on the stirring rate.

As far as r is concerned, there are two limit behaviors, whereby ϵ(0) ≃ 0 and ϵ(r → ∞) ≃ 1, and two threshold values, namely, the just suspended and the uniformly suspended cases. It is worth noting in passing that for very intense stirring, i.e., for r → ∞, secondary nucleation will be controlled by the attrition mechanism,[48] whereby the role of SNIPE will be negligible.

A suspension is considered to be just suspended when no particle remains on the bottom of the agitated vessel for more than 1–2 s; hence, most of the surface of the seeds is exposed to the solution and hence accessible.[56] This occurs for values of the stirring rate around the threshold, rjs, which is defined through the properties of the solution and of the particles, as well as through the agitation intensity, by the Zwietering correlation[56−58]with the reference stirring rate, r0, defined aswhere v is the kinematic viscosity of the liquid, g is the gravitational acceleration constant, ρs and ρ1 are the density of the particle and that of the liquid, respectively, X is the ratio of the mass of suspended solid to the mass of the liquid, L̅v = m4/m3 is the volume-weighted mean size of particles, and dimp is the diameter of the impeller. The parameter pjs is an adjustable parameter whose value depends on the type, geometry, and location of the impeller and typically ranges between 3.4 and 7.1; in fact, it is often estimated from experiments.[58]

A suspension is considered to be uniformly suspended. Hence, particles are not only suspended but also homogeneously distributed in the suspension; thereby, the suspension can be indeed considered well mixed, for stirring rates beyond a threshold, which can be expressed in terms of the just suspended threshold, i.e.[56]where the parameter pus is an adjustable parameter, whose value can be assumed to range between 1.3 and 1.7, with a rather large level of uncertainty.[59]

Based on the physical evidence above, as quantified by eqs and 27, and in consideration of the plausible ranges of values of the empirical parameters pjs and pus given above, we propose to describe the dependence of the efficiency factor ϵ on the stirring rate r using a logistic function, i.e., a sigmoid function. The transition from poor to good suspension (inflection point of the curve) occurs at a value r/r0 = (r/r0)infl, i.e., an adjustable parameter that might have a value around 6 based on the feasible values of the parameters pjs and pus mentioned above. The value of ϵ is essentially 0 (very poor suspension and SNIPE negligible) and 1 (uniform suspension, hence ϵ ≈ 1) for r/r0 < ((r/r0)infl – 3) and r/r0 > ((r/r0)infl + 3), respectively. The functional form of ϵ iswhere pϵ is an empirical parameter that defines the slope of the logistic function at the inflection point and should be larger than or equal to 1.

Surface Nucleation Model

The surface nucleation (SN) rate model[4] suggests that secondary nucleation can occur because part of the 2D nuclei (i.e., growth units formed on the surface of seeds) can be detached and dispersed into the bulk solution, thus becoming secondary nuclei. The corresponding nucleation rate is given by the following expression[3,27,60]where θSN (i = 1, 2) are the parameters to be estimated. The theoretical values for these two parameters are given by[3,27,60]where ψ (0 < ψ < 1) represents the fraction of 2D nuclei being detached from the surface of the seeds, d1 is the diameter of a solute molecule, and D is the diffusion coefficient of the solute molecule in the solution. Note that the value of ψ may depend on the stirring rate, thus allowing the SN rate model to predict a higher secondary nucleation rate at a higher stirring rate in a manner similar to that proposed for the SNIPE model through eqs and 28.

Embryo Coagulation Secondary Nucleation Model

The embryo coagulation secondary nucleation (ECSN) model[16] suggests that secondary nucleation can occur through the rapid coagulation of molecular clusters near the surface of seed crystals caused by the van der Waals forces between the clusters and the seed surface. The corresponding nucleation rate is given by the following expression[16]where θEC (i= 1, 2, 3) are the parameter to be estimated. The parameter θ1EC represents the initial size of the subcritical clusters that coagulate; hence, its value is constrained between the smallest possible size (called the cutoff size in the original article[16]) and the largest possible size (i.e., the critical nucleus size). The parameter θ2EC is a lumped parameter defined aswhere A131 is the Hamaker constant determining the intensity of van der Waals energy between the clusters and the seed, and x is the surface-to-surface distance between a seed crystal and a group of clusters that coagulate. Finally, the theoretical value of the parameter θ3EC is given bywhere η is the dynamic viscosity of the solvent. It is worth noting that, although the original ECSN model cannot account for the effect of the stirring rate on the secondary nucleation rate, this effect can be considered in the ECSN model using the efficiency factor ϵ(r), as done in the case of the SNIPE model (see eqs and 28).

Parameter Estimation Procedure: Least Squares

As explained in Section S3 in the Supporting Information, the parameters of the secondary nucleation rate models were estimated from experimental data by applying the method of least squares, while the corresponding confidence intervals were determined using a linear approximation of the models. The time evolution of bulk supersaturation S(t) was used as a response variable in parameter estimation, as commonly done in the literature.[32,40,61,62] To assign equal weight to each experiment independent of its time duration, experimental data with a longer duration were downsampled, thus fixing the number of the data points of each experiment used for fitting a model; this is a data preprocessing method applied also elsewhere.[40,63]

Results and Discussion

Unless otherwise specified, the results reported in this work are based on the physicochemical properties of the system used in the benchmark study[26] (paracetamol in ethanol) at 20 °C (for a summary of these properties, see Table 1 in part II of ref (18)).

Model Verification and Sensitivity Analysis

To verify the SNIPE rate model (eq ) derived in Section , the SNIPE nucleation rate, JSNIPE, and the corresponding KRE nucleation rate, JKRE, calculated from the simulations of the KRE model presented in part II[18] were compared at different conditions. For this comparison, we assumed perfect mixing, thus setting the efficiency factor ϵ equal to unity. To this aim, the nucleation rates JKRE were calculated from KRE simulations performed under various conditions, as explained in Section S4 in the Supporting Information, and their dimensionless version, ζKRE, was compared with the dimensionless version of JSNIPE, denoted by ζSNIPE, with the dimensionless nucleation rate defined as ζ = J/k0C1,e2. With this definition, the parameter k0 drops from the SNIPE model, thus leading to a model with only three parameters: Ω, Est, and lst.

In the parity plot of Figure , the nucleation rates ζKRE and ζSNIPE are plotted against each other. All of the points in the (ζKRE, ζSNIPE) plane fall on the diagonal line ζSNIPE = ζKRE, thus verifying that the description of the SNIPE rate model and that of the KRE model developed in part II[18] are consistent. Note that, as indicated by an arrow in Figure , the minimum nucleation rate was obtained when the stabilization effect was not considered in the models by setting either Est = 1 (zero intensity) or lst = 0 (zero thickness) or both; such a minimum value corresponds to the kinetics of primary nucleation.

Figure 3

Dimensionless nucleation rates obtained from the KRE model simulations ζKRE (abscissa) and those from the SNIPE rate model ζSNIPE (ordinate) under different conditions, summarized in Section S4 in the Supporting Information. The diagonal line (ζSNIPE = ζKRE) indicates the agreement in the nucleation kinetics described by the two models.

Furthermore, for a sensitivity analysis of the SNIPE rate model (eq ), the supersaturation and three model parameters were varied: the supersaturation s in the range between 1.3 and 2.4, the dimensionless surface energy Ω among the values {1.5, 3.0, 4.5, 6.0}, the intensity of the stabilization effect Est among the values {1.00, 1.10, 1.17, 1.2, 1.25}, and the range of the stabilization effect lst among the values {0, 1, 10, 100}. The results of the sensitivity analysis are presented in terms of the dimensionless nucleation rate ζSNIPE in Figure , where the solid lines without markers indicate a base case with a set of parameters representing primary nucleation (i.e., no stabilization effect) with a specific value of the dimensionless surface energy: Est = 1, lst = 0, and Ω = 4.5.

Figure 4

Dimensionless nucleation rate ζSNIPE given by the SNIPE rate model (eq ) with different values of the supersaturation s, the intensity of the stabilization effect Est, the range of the stabilization effect lst, and the dimensionless surface energy Ω: (a) Est = 1 and lst = 0, (b) Ω = 4.5 and lst = 10, and (c) Ω = 4.5 and Est = 1.2.

The behavior of the model without the stabilization effect (Est = 1 and lst = 0) is shown in Figure a, where the nucleation rate ζSNIPE increases with increasing supersaturation, s, and decreasing surface energy, Ω. Note that a decrease in the surface energy Ω increases the nucleation rate ζSNIPE by orders of magnitude while making the dependence of the nucleation rate ζSNIPE on the supersaturation s flatter.

The main characteristics of the model that includes the stabilization effect (Est > 1 and lst > 0) are illustrated in Figure b and 4c, where increasing the intensity, Est, and range, lst, of the stabilization effect enhances the nucleation rate ζSNIPE, mainly at a low level of supersaturation s. Note that the intensity of stabilization (Est) has a major effect on the SNIPE nucleation rate (see Figure b), while the range of stabilization (lst) has a minor effect (see Figure c), as also shown in Section 4.2 of part II.[18]

Regardless of the values of Est and lst, the nucleation rate given by the SNIPE rate model converges to that of the base case (i.e., primary nucleation) as the supersaturation s increases. This unique and useful feature of the SNIPE rate model allows us to describe both primary and secondary nucleation using the same nucleation rate model, which transitions continuously from the secondary nucleation rate to the primary nucleation rate with increasing supersaturation s. Note that this feature is absent in the other secondary nucleation models, which require the use of both primary and secondary nucleation models to describe crystallization processes.[26,28,64]

Model Validation and Comparison

The objective of this section is to validate the SNIPE rate model (eq ) by proving its capability to describe experimental data through fitting and to assess its goodness-of-fit and theoretical consistency in comparison with the other surface-induced secondary nucleation models. To this aim, a comparison of the secondary nucleation models is presented (see Section ), followed by a discussion on the fidelity of the fitted SNIPE model to the experimental data of the benchmark study (see Section ).[26]

It is worth noting that, under most experimental conditions in the benchmark study, the reference stirring rate r0 is 31 rpm (calculated from eq ); hence, according to eq , the threshold stirring rate for establishing the just suspended condition, rjs, ranges between 107 and 223 rpm. Since this upper bound of rjs is very close to the stirring rate employed in the reference experiments (i.e., 200 rpm), it is very likely that the suspensions in the reference experiments were under a good mixing condition. Accordingly, in the following analysis, the values of the efficiency factor ϵ are varied between 0.5 and 1. For the sake of brevity, results based on ϵ = 1 are discussed unless otherwise mentioned.

Comparison of Secondary Nucleation Models

As explained in Section , the SNIPE mechanism can be viewed as primary nucleation enhanced by the stabilization effect; so, in the case of a well-mixed suspension, its rate expression (eq ) depends on four parameters: two parameters (Est and lst) to quantify the stabilization effect and the other two (k0 and Ω) to reflect primary nucleation kinetics. In this study, the latter (k0 and Ω) were determined independently from the former (Est and lst) because the parameters k0 and Ω can be estimated from a set of primary nucleation rates, thus allowing us to reflect the primary nucleation kinetics of the studied system correctly. To consider the uncertainty in the estimated values of k0 and Ω originating from the inherent uncertainty in the measurement of nucleation rates,[65] two different sets of nucleation rate data were employed, thus yielding two sets of primary nucleation parameter values (see Section S5 in the Supporting Information). Since the parameter values of the two sets are very similar, only the first set was used in the following unless otherwise stated. Besides, the parameter lst was set to value 4 because the value of lst calculated in part I[17] ranges between 2 and 7 and because, in this interval of values of lst, the nucleation rate barely changes (recall from Figure c that the parameter lst has a minor effect on the nucleation rate). Consequently, the parameter Est was the only parameter of the SNIPE rate model estimated by fitting the experimental data of the benchmark study.[26]

The parameter estimation results for all three secondary nucleation models are summarized in Table , including the estimated values of all model parameters (type “est.”), their theoretical values (type “calc.”), and the corresponding values of the objective function H(q) (eq S4 in Section S3 in the Supporting Information) that have been minimized by the method of least squares. For all three models, the minimized values of the objective function are similar, thus suggesting that the statistical goodness-of-fit of all models is comparable. The theoretical consistency of each model is assessed by comparing the estimated parameter values with the corresponding theoretical estimates.

Table 2

Parameter Estimation Results for All Surface-Induced Secondary Nucleation Modelsa

model	quantity	type	unit	value	95% confidence interval
SNIPE	H(θ̑)			0.115
	Est	est.		1.23	±0.01
		calc.		1.17
SN	H(θ̑)			0.117
	θ1SN × 10–5	est.	m–2 s–1	1.2	0.2–8.6
		calc.	m–2 s–1
	θ2SN × 102	est.		1.3	0–7.9
		calc.		6.7
ECSN	H(θ̑)			0.133
	θ1EC	est.		1	1–49
		calc.		≥500
	θ2EC	est.		9.2	0–255.8
		calc.		1.6
	θ3EC	est.	m–2 s–1	1.2 × 1010	9.7 × 10–12 to 1.5 × 1011
		calc.	m–2 s–1	5.3 × 1016

H(θ̑) is the optimized value of the objective function (eq S4 in Section S3 in the Supporting Information).

In the case of the SNIPE rate model, the parameter estimation result was Ȇst = 1.23 ± 0.01, with the narrow confidence interval (±0.01) indicating a high level of confidence in the estimated value. The value of Êst changed only about ±5% even when the value of lst was increased by a factor of 10 (from the value of 4) and when the second set of the estimated primary nucleation parameter values was used. Analogously, halving the value of the efficiency factor ϵ to 0.5 increased the value of Ȇst by around 3%. We underline that all of the values of Ȇst, spanning from 1.1 to 1.3 at different conditions, are physically reasonable, as elaborated in part I of this series.[17]

Likewise, the estimated values of the two parameters for the surface nucleation (SN) model appear to be consistent with the theory: the efficiency factor ψ (=5 × 10–24) calculated from the estimated θ̑1SN using eq is of a comparable order of magnitude as the corresponding estimates reported in the literature,[27,60] let alone within the physically acceptable range (i.e., 0 < ψ < 1).[4] Furthermore, the confidence interval of the estimated θ̂2SN ([0, 0.079]) includes its theoretical value (0.067) given by eq .

On the contrary, in the case of the ECSN model, two of the three estimated parameters significantly deviate from their theoretical values: the estimated θ̑1EC (=1) is at least 2 orders of magnitude smaller than the theoretical minimum value of about 500 (called the cutoff size in ref (16), where the ECSN model was originally developed). Furthermore, the estimated θ̑3EC (=1.2 × 1010) is 16 orders of magnitude smaller than the theoretical value θ3EC (=5.3 × 1026) calculated using eq . Only θ̑2EC (=9.2) is of the same order of magnitude as the theoretical estimate calculated from eq , where the physically reasonable values of the Hamaker constant A131 = 2 × 10−20 J (see ref (66)) and of the coagulation distance x = d1/2 (see ref (16)) are employed. This clearly suggests that the physical mechanism described by the ECSN model is not consistent with the experimental data of the benchmark study.[26]

It is remarkable that the SNIPE rate model with only one free parameter (Est) fitted the benchmark experimental data[26] and that the surface nucleation model with two free parameters (θ1SN and θ2SN) did the same. Moreover, it is notable that the resulting estimate (Ȇst) has a physically meaningful value with a narrow confidence interval. This was possible mainly because the primary nucleation kinetic parameters (k0 and Ω) of the SNIPE model were estimated independently from Est using a set of primary nucleation rate data from the literature.[35] In cases where relevant primary nucleation rate data are absent, a set of three parameters {k0, Ω, Est} needs to be estimated simultaneously from secondary nucleation experiments. However, increasing the number of parameters to be estimated makes the estimates more uncertain and the parameter estimation problem more challenging (eq S4 in Section S3 in the Supporting Information). Accordingly, to estimate physically meaningful parameter values with high certainty, it is highly recommended to estimate the primary nucleation parameters (k0 and Ω) from primary nucleation rate measurement and then to estimate the parameter Est from secondary nucleation experiments.

Validation of the SNIPE Rate Model

In Figure , the concatenated time series of the experimental data (circle markers) and the corresponding predictions of a PBE using the estimated SNIPE rate model parameters (lines) are shown, with wnuclei = m3nuclei/(m3nuclei + m3seed) representing the mass fraction of nuclei. Since the fits of all three secondary nucleation rate models to the experimental data are comparable with each other (see Section ), only the fitting results of the SNIPE are illustrated for brevity. It is worth understanding from Table that the effect of the initial supersaturation on the crystallization process can be seen from a comparison between experiments E1 and E2, the effect of the initial mass of seeds from experiments E2–E4, and the effect of the seed size from experiments E2 and E5.

Figure 5

Concatenated time series of the experimental data used for fitting secondary nucleation models (circle marker) and the corresponding predictions of a population balance equation with the fitted SNIPE rate model (lines): (a) the bulk supersaturation S, (b) the second moment of seeds m2seed, (c) the nucleation rate J, and (d) the mass fraction of nuclei wnuclei. The conditions of each experiment, indicated by labels (E1–E5), are summarized in Table . The experimental data are from ref (26).

In Figure a, there is a reasonable agreement between the experimental data and the model predictions in all five experiments, thus proving the capability of the SNIPE model to describe crystallization processes quantitatively. In Figure b, the second moment of the seeds, m2seed, increases monotonically over time in all experiments simply due to the growth of the seeds.

Figure c shows the time evolution of the nucleation rate J, with some cases having a maximum in the middle of the experiment, as also observed in the benchmark study (see Figure 14 in ref (26)). The maximum in the nucleation rate J can appear because the nucleation rate J is positively correlated to the second moment of seeds m2seed and the supersaturation S and because m2seed increases and S decreases over time. Besides, (i) decreasing the initial supersaturation (S0) reduces the initial nucleation rate (compare E1 and E2), (ii) adding a larger amount of seeds (i.e., a higher value of M0seed) increases the initial nucleation rate while shortening the time to reach a maximum in the nucleation rate J (compare E2–E4), and (iii) reducing the size of seed crystals increases the initial second moment of seeds and thus enhances the initial nucleation rate (compare E2 and E5).

In Figure d, the mass fraction of nuclei wnuclei increases monotonically over time in all cases due to the growth of nuclei. Note that the final value of wnuclei in each experiment highlights the importance of secondary nucleation in the corresponding experiment: a small final value of wnuclei indicates that secondary nucleation has been suppressed in the corresponding experimental conditions. For instance, lowering the initial supersaturation S0 suppressed secondary nucleation (compare E1 and E2) by slowing down nucleation kinetics significantly, which is consistent with the findings in the benchmark study (see Section 5.1 in ref (26)). Analogously, increasing the initial mass of the seeds (M0seed) suppressed secondary nucleation (compare E2–E4), which is also reported in the reference study (see Figure 10 in ref (26)). Lastly, decreasing the size of seeds (i.e., increasing the total surface area of the seeds) suppressed secondary nucleation (compare E2 and E5) because the increased surface area accelerates the consumption of supersaturation by growth so that a smaller amount of nuclei can form. This result is also in agreement with the benchmark study[26] and with results reported by other researchers.[67]

In summary, the predictions of a PBE with the fitted SNIPE model (illustrated in Figure ) are physically sound and they are consistent with the findings in the benchmark study,[26] thus validating the SNIPE rate model. It is worth noting that the validity of the fitted model is restricted to the range of experimental operating conditions used for fitting, which is admittedly narrow, as only two levels of initial supersaturation and two different seed populations are considered. Moreover, a cross validation remains to be performed by utilizing the fitted model to describe a new set of experimental data that have not been used for fitting, as done elsewhere;[63,68] this is, however, beyond the scope of this work.

Conclusions

A SNIPE rate model (secondary nucleation rate model based on the SNIPE mechanism) has been derived theoretically, verified numerically, and validated experimentally. The theoretical derivation showed that the SNIPE mechanism can be described as enhanced primary nucleation, with the magnitude of the enhancement depending on the intensity and range of the stabilization effect originating from the interparticle energies between the surface of seeds and nearby molecular clusters. Consequently, the derived model has four parameters in the case of a sufficiently agitated suspension: two parameters (k0 and Ω) reflecting primary nucleation kinetics and the other two (Est and lst) accounting for the stabilization effect.

The SNIPE rate model was verified by showing that the nucleation kinetics predicted by the model quantitatively agrees with those described by the kinetic rate equation model developed in part II.[18] Furthermore, a sensitivity analysis of the rate model demonstrated that the model predictions for the nucleation rate depend strongly on the dimensionless surface energy (Ω) and the intensity of the stabilization effect (Est) but rather weakly on the effective range of the stabilization (lst).

The experimental validation of the SNIPE rate model was performed by fitting the model to time-resolved experimental data of secondary nucleation. When fitting the SNIPE model to these data, only one model parameter (Est) was varied and its optimal value was estimated by the method of least squares, while the values of the other three parameters were determined independently. Remarkably, the SNIPE model with only one free parameter (Est) fitted the experimental data as well as the other two surface-induced secondary nucleation models (i.e., surface nucleation model[4] and the embryo coagulation secondary nucleation model[16]) that require two or three parameters, and the estimated Est was of a comparable magnitude to the theoretical values presented in part I.[17]

This article highlights the practical applicability and theoretical consistency of the SNIPE theory, whose thermodynamic and kinetic aspects have been addressed in parts I[17] and II[18] of this series, respectively. We underline that the SNIPE rate model can be useful in describing the kinetics of nucleation in various crystallization processes, for instance, by explaining secondary nucleation in a reactor where the attrition mechanism cannot prevail (e.g., in a plug-flow reactor or in an air-lift crystallizer[29,60]), by describing both primary and secondary nucleation with a single nucleation model, and by rationalizing the formation of a secondary nucleus whose polymorphic/chiral crystal structure is different from that of a seed.[17,18]