Diffusion-driven growth of calcium carbonate polymorphs in microchannels

Abstract

We propose a novel approach to characterize the growth of individual crystals. Calcium chloride and sodium carbonate solutions have been injected into a Y-shaped microfluidic channel at various stoichiometric ratios, and the development of calcium carbonate has been monitored. The formation of calcite and vaterite depends not only on the stoichiometric ratio of the reactants but also on the region of the reactor where they form. From the crystal images, we have mapped the surface growth of the particles and have shown that closer to the confluence of the microchannel the crystal growth is significant. Both morphologies mainly form in the carbonate-rich zone, supported by numerical modeling. Moreover, the side growth of the calcite particles is diffusion-controlled and independent of the crystal orientation and the stoichiometric ratio of the reactants injected.

---

1 Introduction

Crystallization is an abundant phenomenon both in industrial processes and in nature. Precipitation reactions are used industrially to produce proteins,¹ polymers,² and paint pigments³ among other products, however, the formation of solid particles is not always beneficial, for example in the case of limescale⁴ or kidney stone⁵ accumulation. To improve the yield or the selectivity in industrial processes, or to decrease the growth of undesired crystals, a deeper understanding of crystallization is inevitable.⁶

One way to comprehend the formation of the crystalline phase is to study its kinetics, which is typically divided into nucleation and growth processes. The birth of crystals can be quantitatively described by the nucleation time, which can be determined experimentally by monitoring the turbidity with a spectroscope,^7–9 a fast camera,¹⁰ by focused beam reflectance count,¹¹ or by conductivity measurements.¹² There are various techniques, which provide more details on the growth of particles, such as the one developed by Mullin et al. by which the growth of crystal faces can be determined by a simple laboratory crystallization apparatus.¹³ More advanced methods use atomic force microscopy¹⁴ to capture the formation of single layers from solution.

Microfluidic devices can be advantageous when studying precipitation since they require small amounts of chemicals and offer well-defined fluid flow with an efficient heat and mass transfer compared to large-scale conventional techniques.¹⁵ Therefore, microchannels have been widely used in crystallization to synthesize active pharmaceutical ingredients,¹⁶ nanomaterials¹⁷ or to mimic and understand kidney stone formation.¹⁸

We have developed a strategy where it is possible to characterize the growth of individual crystals profoundly. Calcium carbonate has been selected as a model precipitate, since its flow-driven nucleation, growth,¹⁹ and polymorph selectivity²⁰ in microchannels have already been studied. With our method, we aim to map the crystal surface growth profile in the entire channel. Moreover, we determine the growth of crystal faces separately and investigate the effect of flow direction and concentration on the development of particle sides. A numerical modeling study has also been carried out to corroborate the experimental findings.

2 Experimental

The growth of calcium carbonate particles was monitored using the experimental setup shown in Fig. 1. From the reactants, CaCl₂·2H₂O (VWR-Life Science, ACS) and Na₂CO₃ (Honeywell, ACS, ≥99.5%), 7.5 mM and 15 mM solutions were prepared. The pH of sodium carbonate solutions was set to 10 by the dropwise addition of 1 M HCl solution (VWR-Life Science, 37%) using a pH Meter (Thermo Orion). The reactant solutions were filtered with syringe filters (VWR) with a pore size of 0.22 μm. A commercially available microfluidic chip (microfluidic ChipShops, H-shaped Channel Chips) made out of poly(methyl methacrylate) served as the reactor. The cross-section of the reactor where the two reactants were pumped in was 75 μm × 75 μm, which at the confluence widened to 150 μm × 75 μm, and the channel length was 1 cm. The solutions were injected into the microchannel using a syringe pump (kdScientific Legato 180) with 0.5 μL min⁻¹ flow rate per channel, naturally setting a time-independent concentration and pH gradient. To record the experiments, a camera (Nikon LV-TV) mounted on a microscope (Nikon Eclipse Ts2R) was utilized. From the recorded grayscale images various properties of the individual crystals on the bottom of the channel have been determined with an in-house program. The boundaries of particles have been localized from the grayscale minima and the cross section areas have been calculated using the Shoelace formula²¹ on the boundaries. In addition for calcite particles, a quadrilateral was fitted to the boundary points using Hough transformation,²² allowing the temporal monitoring of crystal edge/side evolution. The precise position of the crystals have been determined along the channel from the images using reference points.

Fig. 1 The experimental setup.

3 Modeling

Toward a deeper understanding of the experimental observations, numerical modeling of the spatiotemporal distribution before precipitation in a microchannel has been carried out mimicking the experiments where two miscible liquids (aqueous solutions of calcium chloride and sodium carbonate) are injected into the microreactor. In water carbonate ions can be protonated via the following reactionswhere K_H1 and K_H2 are the corresponding protonation constants. Calcium ions can form CaOH⁺ complex viawhere β is the stability constant for the complexation reaction. By taking these fast equilibria and the autoprotonation of water

H₂O ⇌ H⁺ + OH⁻, (4)

into account, the mass balances

c_C = [CO₃²⁻] + [HCO₃⁻] + [H₂CO₃], (5)

c_Ca = [Ca²⁺] + [CaOH⁺], (6)

and the charge balanceallow the reduction to three variables: c_C and c_Ca are the total carbonate and calcium concentration, respectively, and [A⁻] is the concentration of the counter ion in the monoprotic acid added to the Na₂CO₃ solution. In eqn (7) the equilibrium concentration of sodium and chloride ions can be substituted with c_C and c_Ca, respectively, as they originate from the Na₂CO₃ and CaCl₂ stock solutions leading to

Since the equilibrium reactions are much faster than the transport processes, the pre-equilibrium approximation can be applied to calculate the concentrations of the 7 species (CO₃²⁻, HCO₃⁻, H₂CO₃, Ca²⁺, CaOH⁺, OH⁻, H⁺) from eqn (1)–(7) as c_C, c_Ca, and [A⁻] are known.

The temporal change in the concentrations in the presence of transport processes is calculated by solving the differential mass balance equation for c_C, c_Ca, and [A⁻]

where the second term on the left side of the equation is the contribution of advection, while on the right side, the diffusion term is presented only as pre-equilibrium approximation holds. Since the flow in the channel is laminar, the Navier–Stokes equation for incompressible fluids (∇·u⃑ = 0) has been usedwhere p is the pressure, u⃑ is the velocity of the fluid flow, ν = 10⁻⁶ m² s⁻¹ is the kinematic viscosity of the liquid, and ρ is the solution density. Dilute solutions are used, therefore the densities are assumed to be equal.

The partial differential equations in eqn (9) and (10) are discretized and solved by OpenFoam software.²³ An operator splitting technique is used where the convective and the diffusive terms are calculated separately, followed by an update of the concentration fields according to eqn (1)–(6) and (8) with the KINSOL module of the SUNDIALS software package.^24,25 Finally, to obtain the spatiotemporal distribution of the supersaturation S at each time step, the concentrations of free calcium and carbonate ions in every cell were used to calculate the supersaturation S as

where K_sp is the solubility product of calcium carbonate.

The reactor was modeled with the same spatial dimensions as the one in the experiments (see Fig. 2). The inlet channels were discretized to (10, 20, 20), while the main one to (200, 40, 20) volume cells in (X, Y, Z) directions. At the walls, Zero Gradient boundary condition was used for the concentration fields and No Slip for the velocity. At the inlets and outlets, fixed values were set for the flow rate and the concentrations in accordance with the experiments. The densities of the species in these dilute solutions were set to ρ = 1000 kg m⁻³. As diffusion coefficients for c_C, c_Ca and A⁻, 0.955 × 10⁻⁹, 0.793 × 10⁻⁹ and 2.03 × 10⁻⁹ m² s⁻¹ were used, respectively, the following equilibrium constants were set: pK_H1 = 10.33, pK_H2 = 6.35, log₁₀(β) = 1.3, pK_w = 14, pK_sp = 8.35.²⁶ For solving the differential equation (eqn (9) and (10)) Euler method was used with a step size of Δt = 1 × 10⁻⁴ s.

Fig. 2 The scheme of the reactor used in the simulations.

4 Results and discussion

4.1 Characterization of crystal formation

The crystal formation and its growth have been investigated in a microfluidic reactor at three different concentration ratios. The results are presented with microscopic images, where at the upper inlet of the Y-shaped channel CaCl₂, while at the lower inlet, Na₂CO₃ solutions are injected. Fig. 3 shows the reactor at the confluence in the case of 2 : 2 (a), 1 : 2 (b), and 2 : 1 (c) stoichiometric ratios. At 2 : 2 composition, small particles form in the middle of the channel and a sharp line appears at the contact line due to the different refractive indexes of the solutions. In general, 10–20 larger and more isolated crystals nucleate in rhombohedral calcite morphology, except in the case of a 2 : 1 ratio, where 2–3 sphere-like vaterite particles also develop as confirmed by Raman microscopy (see Fig. S1 in ESI†). No amorphous calcium carbonate (ACC) or aragonite has been found throughout the experiments at the flow rate and solution composition used.

Fig. 3 Top view images of the microfluidic channel near the confluence at [Ca²⁺] : [CO²⁻₃] = 2 : 2 (a), 1 : 2 (b), and 2 : 1 (c) ratios.

Fig. 4 depicts the top view of the reactor segment 5 mm far from the inlet with various chemical compositions. The contact line is invisible and only the larger crystals remain in the case of the 2 : 2 stoichiometric ratio. Comparing the three images in Fig. 4, we find that, independently of the composition, particles tend to nucleate at the side where the Na₂CO₃ solution is injected into the reactor.

Fig. 4 Top view images of the middle part of the microfluidic channel at [Ca²⁺] : [CO₃²⁻] = 2 : 2 (a), 1 : 2 (b), and 2 : 1 (c) ratios.

Most of the particles close to the outlet, presented in Fig. 5, also nucleate at the carbonate-rich zone but there are some particles in the more acidic, upper regions. The number of crystals is greater than in the other segments of the reactor. Moreover, vaterite crystals appear in the case of 2 : 2 composition, suggesting spatial morphology control besides the stoichiometric composition-driven one.

Fig. 5 Top view images of the microfluidic channel near the outlet at [Ca²⁺] : [CO₃²⁻] = 2 : 2 (a), 1 : 2 (b), and 2 : 1 (c) ratios.

The experimental findings are corroborated with numerical modeling mimicking the experiments. The spatial distribution of the supersaturation S is calculated (see Fig. 6) since nucleation takes place where supersaturation is high. At a 2 : 2 stoichiometric ratio, the highest supersaturation is at X = 0.2 mm from the confluence as shown by a dashed line in Fig. 6a, which explains the formation of small particles at the middle of the reactor close to the confluence in Fig. 3b.

Fig. 6 The spatial distribution of the supersaturation at the initial part of the reactor in the XY-plane (a) in case of 2 : 2 stoichiometric ratio at t = 10 s. The red colored lines indicate the positions where the supersaturation profiles are determined. The supersaturation S transverse to the flow at various stoichiometric ratios at t = 10 s (b). The dashed lines represent the supersaturation profiles at X = 0.2 mm, while the solid ones at X = 1.0 mm.

Farther from the inlet at X = 1.0 mm, supersaturation significantly increases in the region where sodium carbonate is injected, corroborating the experimental findings, where nucleation is found to be dominant at the carbonate-rich region (see Fig. 4). In this region the concentration of deprotonated carbonate ions is higher, hence supersaturation is reached sooner. Furthermore, supersaturation spreads out in the Y direction on moving away from the injection points (see Fig. 6a) due to the mixing effect of diffusion perpendicular to the flow. Therefore, close to the outlet particles can also nucleate at the calcium-rich zones as seen in the experiments in Fig. 5.

4.2 Characterization of crystal growth

Three parallel experiments have been conducted for each stoichiometric ratio, and 4–22 crystals, which are at least 10 μm far from each other, so their growth is independent at the timescale of the experiments, have been selected for characterization. First, the temporal changes in the particle area in the XY-plane have been investigated. Since the cross section areas of calcites and the orthodrome area of vaterites change linearly with time, their growth rates are calculated from the slopes of their temporal evolution.

Fig. 7 illustrates the location of the selected particles in the XY-plane with symbol sizes proportional to their growth rate. In all scenarios, the growth rate is the highest closer to the confluence and it is decreasing farther from it, dropping to only 1–2% at the outlet. This can be explained by the changes in the supersaturation, as S is the highest at the inlet of the channel and it gradually decreases towards the outlet. At high supersaturation areas not only nucleation but also growth are promoted. The surface growth rate of calcites varies on a large scale: the highest is in the range of 0.239 ± 0.002 μm² s⁻¹, while the smallest ones are in the range of 0.0071 ± 0.0001 μm² s⁻¹. In the case of vaterite, the growth rate of the orthodrome area is between 0.454 ± 0.002 μm² s⁻¹ and 0.023 ± 0.001 μm² s⁻¹. All the experimentally determined surface growth rates for both calcite and vaterite crystals are summarized in Tables S1–S4 of the ESI.†

Fig. 7 The positions of the calcite crystals in the channel in case of 2 : 2 (a), 1 : 2 (b), 2 : 1 (c) stoichiometric ratios and the vaterite polymorph at 2 : 1 ratio (d). The symbol size is proportional to the growth rate of the crystalline particles. The three different colors represent three parallel experiments. The X = 0 cm corresponds to the confluence of the Y-shaped channel, while X = 1 cm to the outlet. Sodium carbonate solution is injected at 0 > Y ≥ −75 μm, whereas at 0 < Y ≤ 75 μm the calcium chloride solutions are injected into the reactor.

In the case of the calcite crystals, not only the surface growth rate but also the temporal change in each face size can be determined. The growth of the particle sides in the XZ- and YZ-planes can be substituted with the edge growth in the XY-plane. Two calcite crystals are presented in Fig. 8 as typical examples: in the left one (Fig. 8a) the edges grow significantly, while in the right one (Fig. 8d) barely. The former grows approximately 1.5 μm under 200 s, while the latter only 0.2 μm during the same time interval. Moreover, the growth profiles of their sides are distinctly different by comparing the two crystals. The sides grow together as shown in Fig. 8b which results in a regular rhomboidal shape, and the temporal change in the length of the edges of the crystal is proportional to the square root of time (Fig. 8c), which indicates diffusion control.²⁷ In the case of the other particle (Fig. 8d) the faces grow linearly and in pairs as illustrated in Fig. 8f, therefore, the crystal shape is more elongated. The linear growth is only observed when particle growth is monitored from a larger initial size (l > 10 μm). However, the linearity may be apparent, since the final part of the square root growth seems linear on a short timescale.

Fig. 8 Microscopic images of two calcite crystals at two-time instances (a) and (d), temporal changes in the contour lines of the crystals (b) and (e), and the growth of their sides as a function of time (c) and (f) with t₀ denoting the starting time of monitoring.

As seen from Fig. 8, the growth of the particle edges is not uniform (see the growth rates of all four different sides of the investigated particles summarized in Tables S5–S7 in the ESI†) as the faces can grow either together or in pairs resulting in rhomboid or less regular crystalline shapes, respectively. Therefore, we have also investigated what differentiates the growth of the edges and governs the final crystal shape. Since from previous studies,²⁸ it had been shown for lithium phosphate precipitates grown at 2–4 μL min⁻¹ flow rates, that in microfluidic channels the direction of the crystal growth is influenced by the flow of the solutions, we have created Fig. 9 to study this effect. The relative standard deviations of the intercept are plotted as a function of the relative standard deviations of the slope in Fig. 9a, which have been obtained from the fitting of the growth of the four sides of the particle. The slopes of the fittings have been used as ordinates in Fig. 9b, while on the abscissa the angle between the fluid flow direction and the crystal face is shown. The symbols corresponding to the three experimental compositions are scattered with no formation of any groups in either presentation suggesting that the growth of the crystal faces in the XY-plane is independent of their orientation compared to the flow direction and of the stoichiometric ratios at the flow condition applied. Moreover, most of the crystal edge growth is below 0.3 μm s⁻¹.

Fig. 9 The relative standard deviations of the intercept as a function of the relative standard deviation of the slope obtained from the fitting of the growth of the four crystal edges (a). The slope of the side growth as the function of the angle between the fluid flow direction and the crystal side (b).

5 Conclusion

We have studied the formation and have characterized the growth of individual calcium carbonate particles in a microfluidic reactor into which calcium chloride and sodium carbonate have been injected at [Ca²⁺] : [CO²⁻₃] = 2 : 2, 1 : 2, and 2 : 1 stoichiometric ratios. Mostly calcite crystals appear in the microchannel along the flow direction, except at a 2 : 1 stoichiometric ratio, where vaterite particles are also present as confirmed by Raman microscopy. The majority of the crystals form on the side of the carbonate-rich zone where supersaturation is at maximum and corroborated by numerical calculations.

With the use of our experimental setup and in-house evaluation software, we have quantified the surface growth of the particles. From the positions of the crystals in the reactor we have found that the growth is linear and it is dominant close to the confluence. The growth rates of each side of calcite particles scale with the square root of time, i.e., they are diffusion-driven. Furthermore, at the applied injection rates the growth is independent of both the crystal orientation and the solution composition. Further studies would be necessary to investigate the possible effect of shear on the crystal growth for individual particles.

We believe that the developed method is useful in other experimental setups to determine the surface growth of crystals or to obtain the growth rates of individual faces in the case of even irregularly shaped crystals. These data can be further implemented in theoretical dynamic studies as valuable constant parameters.

Data availability

The data supporting this article have been included in the main text and in the ESI.†

Author contributions

R. M. Á.: investigation, formal analysis, visualization P. P.: investigation, software, formal analysis, writing – original draft D. H.: conceptualization, software, methodology, supervision Á. T.: conceptualization, methodology, supervision, funding acquisition. All authors edited and reviewed the text.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

The authors thank the financial support from the National Research, Development, and Innovation Office (K138844 and TKP2021-NVA-19). P. P. acknowledges the support from the New National Excellence Program of the Ministry for Culture and Innovation from the source of the National Research, Development, and Innovation Fund (ÚNKP-23-3-SZTE-492). The authors are also grateful to Ádám Balog for carrying out the Raman microscopy measurements. The authors thank the University of Szeged Open Access Fund (7391) for support. We acknowledge the KIFÜ's HPC supercomputing infrastructure.

References

1. H. Huettmann, S. Zich, M. Berkemeyer, W. Buchinger and A. Jungbauer, Design of industrial crystallization of interferon gamma: Phase diagrams and solubility curves, Chem. Eng. Sci., 2015, 126, 341–348.
2. W. Zhang, H. Zhang, L. Yang, Y. Tang and P. Tang, Enhanced crystallization and properties of poly(ethylene terephthalate) nanocomposites with zeolites from 3D to 2D topologies, J. Ind. Eng. Chem., 2022, 109, 510–520.
3. S. Švarcová, E. Kočí, P. Bezdička, S. Garrappa, L. Kobera, J. Plocek, J. Brus, M. Šťastný and D. Hradil, Uncovering lead formate crystallization in oil-based paintings, Dalton Trans., 2020, 49, 5044–5054.
4. D. Georgiou, D. Bendos, M. Kalis and C. Koutis, Removal and/or prevention of limescale in plumbing tubes by a radio-frequency alternating electric field inductance device, J. Water Process Eng., 2018, 22, 34–40.
5. F. L. Coe, A. Evan and E. Worcester, Kidney stone disease, J. Clin. Invest., 2005, 115, 2598–2608.
6. J. Ulrich and M. Jones, Industrial crystallization: Developments in research and technology, Chem. Eng. Res. Des., 2004, 82, 1567–1570.
7. M. Emmanuel, P. Papp, G. Schuszter, A. Deák, L. Janovák, A. Tóth and D. Horváth, Nucleation kinetics of lithium phosphate precipitation, CrystEngComm, 2022, 24, 4447–4453.
8. N. P. Das, R. Zahorán, L. Janovák, A. Deák, A. Tóth, D. Horváth and G. Schuszter, Kinetic characterization of precipitation reactions: Possible link between a phenomenological equation and reaction pathway, Cryst. Growth Des., 2020, 20, 7392–7398.
9. W. Guo, W. Xia, K. Cai, Y. Wu, B. Qiu, Z. Liang, C. Qu and R. Zou, Kinetic-controlled formation of bimetallic metal–organic framework hybrid structures, Small, 2017, 13, 1702049.
10. R. Zahorán, A. Kukovecz, A. Tóth, D. Horváth and G. Schuszter, High-speed tracking of fast chemical precipitations, Phys. Chem. Chem. Phys., 2019, 21, 11345–11350.
11. J. Schöll, L. Vicum, M. Müller and M. Mazzotti, Precipitation of L-glutamic acid: Determination of nucleation kinetics, Chem. Eng. Technol., 2006, 29, 257–264.
12. T. Maqbool, P. Srikiratiwong and H. S. Fogler, Effect of temperature on the precipitation kinetics of asphaltenes, Energy Fuels, 2011, 25, 694–700.
13. J. W. Mullin and A. Amatavivadhana, Growth kinetics of ammonium- and potassium-dihydrogen phosphate crystals, J. Appl. Chem., 1967, 17, 151–156.
14. S. Verma and P. J. Shlichta, Imaging techniques for mapping solution parameters, growth rate, and surface features during the growth of crystals from solution, Prog. Cryst. Growth Charact. Mater., 2008, 54, 1–120.
15. A. Saikia, R. Newar, S. Das, A. Singh, D. J. Deuri and A. Baruah, Scopes and challenges of microfluidic technology for nanoparticle synthesis, photocatalysis and sensor applications: A comprehensive review, Chem. Eng. Res. Des., 2023, 193, 516–539.
16. R. Fortt, R. Tona, P. M. Martin-Soladana, G. Ward, D. Lai, J. Durrant and N. Douillet, Extractive crystallization of cabotegravir in droplet-based microfluidic devices, J. Cryst. Growth., 2020, 552, 125908.
17. P. R. Makgwane and S. S. Ray, Synthesis of nanomaterials by continuous-flow microfluidics: a review, J. Nanosci. Nanotechnol., 2014, 14, 1338–1363.
18. S. Bourg, K. Rakotozandriny, I. T. Lucas, E. Letavernier, C. Bonhomme, F. Babonneau and A. Abou-Hassan, Confining calcium oxalate crystal growth in a carbonated apatite-coated microfluidic channel to better understand the role of Randall's plaque in kidney stone formation, Lab Chip, 2024, 24, 2017–2024.
19. L. Li, J. R. Sanchez, F. Kohler, A. Røyne and D. K. Dysthe, Microfluidic control of nucleation and growth of CaCO₃, Cryst. Growth Des., 2018, 18, 4528–4535.
20. Y. Zeng, J. Cao, Z. Wang, J. Guo and J. Lu, Formation of amorphous calcium carbonate and its transformation mechanism to crystalline CaCO₃ in laminar microfluidics, Cryst. Growth Des., 2018, 18, 1710–1721.
21. B. Braden, The surveyor's area formula, Coll. Math. J., 1986, 17, 326–337.
22. R. O. Duda and P. E. Hart, Use of the Hough transformation to detect lines and curves in pictures, Commun. ACM, 1972, 15, 11–15.
23. H. G. Weller, G. Tabor, H. Jasak and C. Fureby, A tensorial approach to computational continuum mechanics using object-oriented techniques, Comput. Phys., 1998, 12, 620–631.
24. D. J. Gardner, D. R. Reynolds, C. S. Woodward and C. J. Balos, Enabling new flexibility in the SUNDIALS suite of nonlinear and differential/algebraic equation solvers, ACM Trans. Math. Softw., 2022, 48, 1–24.
25. A. C. Hindmarsh, P. N. Brown, K. E. Grant, S. L. Lee, R. Serban, D. E. Shumaker and C. S. Woodward, SUNDIALS: Suite of nonlinear and differential/algebraic equation solvers, ACM Trans. Math. Softw., 2005, 31, 363–396.
26. S. Kotrly and L. Šucha, Handbook of Chemical Equilibria in Analytical Chemistry, Ellis Horwood Ltd, 1985.
27. K. Ananth, J. P. Kavanagh, R. C. Walton and P. N. Rao, Enlargement of calcium oxalate stones to clinically significant size in an in-vitro stone generator, BJU Int., 2002, 90, 939–944.
28. M. Emmanuel, D. Horváth and A. Tóth, Flow-driven crystal growth of lithium phosphate in microchannels, CrystEngComm, 2020, 22, 4887–4893.

---

Footnotes

† Electronic supplementary information (ESI) available. Raman spectra of the polymorphs, the surface and the side growth rates of the particles. See DOI: https://doi.org/10.1039/d4ra07137a

‡ These authors contributed equally to this work.

---