C O 2 , g ↔ K 1 C O 2 , a q (2-1) H 2 O g ↔ K 2 H 2 O a q (2-2)

C O 2 , a q + 2 N H 3 , a q ↔ K 3   N H 2 C O O − + N H 4 + (2-3) C O 2 , a q + O H − ↔ K 4 H C O 3 − (2-4) O H − + H + ↔ K 5 H 2 O (2-5) N H 4 + ↔ K 6 N H 3 + H + (2-6) C O 3 2 − + H + ↔ K 7 H C O 3 − (2-7) C a N O 3 + ↔ K 8 C a 2 + + N O 3 − (2-8) C a C O 3 , a q + H + ↔ K 9 C a 2 + + H C O 3 − (2-9) C a H C O 3 + ↔ K 10 C a 2 + + H C O 3 − (2-10) C a O H + + H + ↔ K 11 C a 2 + + H 2 O (2-11)

C a 2 + + C O 3 2 − ↔ K 12 C a C O 3 , s (2-12)

The geochemical equilibrium package EQ3/6 was used to identify the relevant speciation reactions; reactions which result in ionic species at very low concentrations have been neglected.

The following assumptions have been made:

• N₂ is assumed to be insoluble, and thus it does not absorb in the aqueous solution.

The reactions in the bulk liquid (Eqs 2-3–2-11) are considered to be instantaneous hence they can be assumed to be at equilibrium. The corresponding chemical equilibrium conditions must be fulfilled, which can be written as: K i = ∏ j a j , p r o d u c t ∏ k a k , e d u c t (2-13) a j =   γ j c j (2-14) a j and c j are the activity and the concentration of the ion or molecule i ; K i is the equilibrium constant associated to the corresponding equilibrium (the equilibrium constants are reported in the Supplementary Appendix SA); γ i is the activity coefficient of the ion or molecule i . Literature correlations were used to calculate the activity coefficient for NH₃ (Maeda and Kato, 1995). The activity of water and the activity coefficient of CO_2,aq were calculated according to Wolery (2002). The activity coefficient of CaCO_3,aq was assumed to be equal to one. The B-dot equations, which consider the dominat long-range interactions between the ions and molecules, were used to calculate the activity coefficient of charged species. A complete description and discussion of these equations can be found elsewhere (Wolery, 2002).

In addition, the electro neutrality condition has to be met: ∑ i c i z i = 0 (2-15) where the sum extends on all the ionic (charged) species and z i corresponds to the charge of ion i . The concentration of the ionic species and molecules in solution can be used to define the following four overall concentrations in solution, which in turn will have to satisfy the overall material balances (see below): c C = c C O 2 , a q + c H C O 3 − + c C O 3 2 − + c C a H C O 3 + + c C a C O 3 , a q + c N H 2 C O O − (2-16) c C a = c C a 2 + + c C a O H + + c C a N O 3 + + c C a H C O 3 + + c C a C O 3 , a q (2-17) c N O 3 = c N O 3 − + c C a N O 3 + (2-18) c N H 3 = c N H 4 + + c N H 3 + c N H 2 C O O − (2-19)

The system of equilibrium reactions is solved by making use of the geochemical equilibrium software package EQ3/6, v8.0a (Wolery, 2002). EQ3/6 is a software package for modeling geochemical interactions between aqueous solutions, solids, and gases, following principles of chemical thermodynamics and chemical kinetics, which is maintained and kept updated at the Lawrence Livermore National Laboratory ¹ , and has been used in earlier mineral carbonation studies (Hariharan and Mazzotti, 2017a; Hariharan and Mazzotti, 2017b; Hariharan et al., 2016; Hariharan, 2017; Werner et al., 2014a; Werner et al., 2014b). The equilibrium constant for Eq. 2-3 is provided by Thomsen and Rasmussen (1999).

The equations above involve the concentrations of the 14 molecules and ionic species in solution: c C O 2 , a q ,   c H C O 3 − ,   c C O 3 2 − ,   c C a H C O 3 + ,   c C a C O 3 , a q ,   c N H 2 C O O − ,   c C a 2 + ,   c C a O H + ,   c C a N O 3 + ,   c N O 3 − ,   c N H 4 + ,   c N H 3 ,   c H +   and   c O H − . If the total concentrations of carbon c C , of calcium   c C a , of ammonium c N H 3 , and of nitrate c N O 3 are known, for instance because they are determined through the material balances for the whole system (see below) the liquid speciation equations (the nine Eqs 2-3–2-11) together with the definitions of the overall concentrations (Eqs 2-16–2-19) and the electro neutrality condition (Eq. 2-15) can be solved numerically to obtain the 14 unknown concentrations in solution.

The material balance for water in the reactor can be written as: d V s ρ d t = Q i n − Q o u t (3-1)

The product of the liquid phase volume V s and the liquid phase density ρ corresponds to the amount of liquid present in the reactor. The system is operated at steady state, thus Q i n is equal to Q o u t , which from now are called simply Q. The mass balances for Ca, C, NH₃ and N O 3 − in the aqueous solution are: V s ρ d ( c C a ) d t = Q ( c C a , i n − c C a ) − r c V s (3-2) V s ρ d ( c C ) d t = Q ( c C , i n − c C ) + V s ( f M T − r c ) (3-3) V s ρ d ( c N H 3 ) d t = Q ( c N H 3 , i n − c N H 3 ) (3-4) V s ρ d ( c N O 3 ) d t = Q ( c N O 3 , i n − c N O 3 ) (3-5)

The initial conditions for the system have to be assigned accordingly in case the calculation of the transient between start up and steady state is needed; they are not necessary for steady state simulations. The quantity c i represents the overall concentration of component i in the aqueous solution, i.e., corresponding to the left hand side of Eqs 2.15–2.18.

The gas phase mass balance is d n t o t d t = F i n − F o u t − f M T V s (3-6)

The total number of moles of gas in the reactor’s overhead space, n t o t , is calculated by the ideal gas law n t o t = P R T ( V R − V s ) (3-7) P corresponds to the system pressure, R to the ideal gas constant and T to the reactor temperature. The CO₂ mass balance in the gas phase is d ( n t o t y C O 2 ) d t = F i n y C O 2 , i n − F o u t y C O 2 , o u t − f M T V s (3-8) y C O 2 corresponds to the CO₂ mole fraction in the overhead space. The outlet gas stream is assumed to be dry, thus the CO₂ outlet concentration is calculated as follows: y C O 2 , o u t = y C O 2 y C O 2 + y N 2 (3-9) y N 2 corresponds to the N₂ mole fraction in the overhead space of the reactor, whereby the sum of the molar gas fractions for CO₂, N₂ and H₂O in the overhead space must be equal to one. ∑ i y i = 1 (3-10)

The quantity f M T that appears in Eqs 3-3, 3-7, 3-9 couples the equations for the gas phase and for the solution and needs to be defined based on an underlying mass transfer model; this is discussed in detail in Section 3.4.

The solid phase mass balance is d n C a C O 3 d t = r c V s − k V Q ρ c ρ M C a C O 3 μ 3 (3-11) where n C a C O 3 is the number of moles of CaCO₃ in the reactor, k V the volume shape factor of the crystals, equal to π/6 for spheres, ρ c the mass density and M C a C O 3 the molar mass of CaCO₃, and μ 3 is the 3^rd moment of the crystal size distribution, which must be calculated from the solution of the population balance equation, as discussed in Section 3.5 below. The quantity r c represents the number of moles of CaCO₃ that are transferred to the solid phase from solution per unit time and unit volume of the solution.

A model for the gas-liquid mass transfer rate of CO₂ is adopted, which is based on the two-film theory. The rate of gaseous CO₂ absorbed by the liquid phase can be written as Levenspiel (1999a): f M T =   k L   a e f f   E   ρ ( c C O 2 , a q ∗   − c C O 2 , a q ) (3-12) c C O 2 , a q ∗ = f C O 2 H C O 2 (3-13) k L    is the liquid phase mass transfer coefficient, a eff corresponds to the effective mass transfer area per unit volume of the suspension, E the enhancement factor, f CO 2   is the arithmetic mean bulk gas phase fugacity of CO₂ of the incoming and exiting gas stream, H C O 2 is the Henry constant of CO_2, c CO 2 , a q is the molecular CO₂ concentration in the bulk liquid and c C O 2 , a q ∗ is the concentration of molecular CO₂ at the gas—liquid interface, which is in equilibrium with the gas phase. Note that Eq. 3.14 is the isofugacity for the CO₂ gas-liquid equilibrium referred to by Eq. 2-1. As a consequence of its definition, f MT is given in moles of CO₂ per unit time per unit volume of the solution.

A kinetic expression for the precipitation of CaCO₃ is developed. Calcium carbonate can precipitate in three main anhydrous forms: calcite, aragonite and vaterite, listed here by decreasing stability (Flaten et al., 2009). The two dominant mechanisms during precipitation are nucleation and crystal growth. In the case of the continuous process modeled here, nucleation is described according to the classic nucleation theory. To model the precipitation process of CaCO₃, we assume that the suspension is well mixed and that nucleation and growth are the dominant precipitation mechanisms, i.e., we have decided for the sake of simplicity to neglect the effect of agglomeration and breakage. The corresponding population balance equation can be written as (Davey, 2000; Ramkrishna, 2000; Alvarez et al., 2011): ∂ f ∂ t + G ∂ f ∂ L + 1 τ ( f − f i n ) = 0 (3-21) where f corresponds to the number based particle size distribution in the reactor, G to the size independent growth rate, L to the characteristic length of a particle,   τ to the residence time of the suspension in the reactor, i.e.,   τ = V s ρ s Q , and f i n ( L ) is the number based particle size distribution of the population of particles entering the reactor with the feed stream.

Since the system operates at steady state, and the feed is free of particles, Eq. 3-21 can be simplified to: G d f d L + 1 τ f = 0 (3-22)

With the boundary condition: f ( L = 0 ) = J G (3-23) where J is the rate of particle nucleation and G is the rate of crystal growth. It is worth noting that for the sake of simplicity but without loss of generality we have decided not to distinguish between different polymorphs precipitating; as a consequence the driving force for nucleation and growth is expressed with respect to the crystallization of the stable polymorph, i.e., calcite.

The nucleation rate is given by the classical nucleation theory relationship (Jun et al., 2016; Reis et al., 2018): J = A ⁡ exp ( − 16 π σ 3 v m 2 3 k B 3 T 3 ( log ( S ) ) 2 ) (3-24) where S corresponds to the supersaturation with respect to calcite, σ is the surface energy, v m represents the molecular volume of C a C O 3 , both listed in the Supplementary Appendix SB. k B is the Boltzmann. Considering the boundary conditions, the population balance, Eq. 3-22, can be integrated to obtain the particle size distribution: f ( L ) = J G exp ( − L G τ ) (3-25)

One can see, that as the characteristic length of the particle approaches the size of the smallest particle (L = 0), the corresponding boundary condition is fulfilled. Moreover, the regularity condition is fulfilled, since the value of f(L) becomes zero as the characteristic length approaches infinity (Ramkrishna, 2000).

The growth rate was assumed to be continuous in the entire range of calcite saturation values, meaning that the energy required to integrate a molecule to the solid lattice is low and that every growth unit arriving finds a site in the lattice (Davey, 2000). Thus, the growth rate can be expressed as: G = k c ( S − 1 ) α (3-26) where k c is the rate constant for growth, α corresponds to an empirical exponent, and S is the supersaturation with respect to calcite, which is defined as: S =   a C O 3 2 − a C a 2 + K S P (3-27) where K S P is the solubility product of calcite. Different growth rate models can be found in the literature, some of them are more of a theoretical nature while others more of an empirical one. Brečević (2007) developed a correlation for k c at 25°C and for S ranging from 1 to 3 as a function of the ionic strength assuming α equal to 2. However, the discussed system operates at much higher values of supersaturation. Chen et al. (1997) claims that for gas—liquid systems, α is equal to 2.9. In addition, the surface integration is supposed to be rate limiting (also supported by others, see Nancollas and Reddy, 1971; Kazmierczak et al., 1982). Wolthers et al. (2012) develop an empirical model for k c valid for ionic strengths ranging from 0.001–0.7 M: k c = I − β p H − γ ( a C a 2 + a C O 3 2 − ) − Δ (3-28)

They estimated β to be equal to 0.004, γ to be equal to 10.71 and Δ to be equal to 0.35. Furthermore, they assume α equal to 2. Nehrke et al. (2007) identified experimentally that the growth rate of calcite has a maximum at a ratio of the a C a 2 + to a C O 3 2 − equal to one.

Finally, the jth moment of the crystal size distribution μ j is defined as: μ j = ∫ 0 ∞ L j f ( L ) d L (3-29)

Using Eq. 3-26 the integral can explicitly be calculated in the case of interest as: μ 2 = 2 J τ 3 G 2 (3-30) μ 3 = 6 J τ 4 G 3 (3-31)

Thus, the moments of the crystal size distribution are function of the residence time and of the composition of the solution, through the values of J and G that depend in turn on the supersaturation S . Moreover, Eq. 3-11 can be solved at steady state, thus yielding the following expression for r c : r c = 6 k V ρ c M C a C O 3 J τ 3 G 3 (3-32)

It is worth noting that the quantity r c can be calculated independently of Eq. 3-12 by considering that CaCO₃ uptake by the solid phase occurs through crystal growth only, since new nuclei are assumed to have L = 0 hence zero mass. To this aim, one must calculate the product of the growth rate and the total crystal surface, which is proportional to the second moment of the particle size distribution, μ 2 : the proper expression is r c = 0.5 k A ρ c M C a C O 3 G μ 2   , which yields Eq. 3-32 when substituting Eq. 3-31 and considering that k A k V = 6 for particles that do not change shape as they grow.

In the following, the KOPs and the KPIs, which have been developed in the scope of a previous experimental study (Tiefenthaler and Mazzotti, 2021) will be introduced as they are used in the Section 4 of the manuscript, particularly for the comparison between simulations and experiments. The stoichiometric ratio in the feed to the reactor,   ψ F , is defined as: ψ F = F i n y C O 2 , i n Q i n c C a , i n (3-33)

The CO₂ absorption efficiency η C O 2 , a b s . is defined and calculated as: η C O 2 , a b s . =   f M T V s F i n y C O 2 , i n (3-34)

The product of the feed stoichiometry ψ F and the CO₂ absorption efficiency η C O 2 , a b s . determines the actual ratio of the two reactants that enters the aqueous solution, ψ , which is thus: ψ = f M T V s   Q i n   c C a , i n = η C O 2 , a b s .     ψ F (3-35)

The efficiency of the precipitation of CO₂ from the aqueous solution, η C O 2 , p r e c . , is defined as: η C O 2 , p r e c . = r c   f M T (3-36)

Finally, the CO₂ mineralization efficiency of the process with respect to CO₂, η C O 2 , is defined as: η C O 2 = η C O 2 , a b s .   η C O 2 , p r e c .   = r c   V s F i n   y C O 2 , i n (3-37) whilst the overall efficiency of the process with respect to calcium, η C a , is calculated as: η C a = r c V s Q i n   c C a , i n (3-38)

Note that η C a =   η C O 2 ψ F (3-39)

The process model consists of a thermodynamic model, for the speciation, of material balances and of a population balance equation, as well as of kinetics constitutive equations for mass transfer and for nucleation and growth in the crystallization process. By solving the model equations, the key quantities f M T , r c and L 50 are calculated as function of the operating conditions. It is worth noticing that all these quantities are coupled, i.e., they are interdependent. As these quantities have been measured in our recent experimental study, the empirical model parameters k L  and a eff as well as A and α have been estimated by minimizing the sum of the squares of the difference between the calculated and the measured values of these quantities, or in other terms the sum of the mean relative deviations (MRD), which are calculated for each measured quantity as follows: M R D = 1 N exp   ∑ i = 1 N exp ( x i mod e l − x i exp ⁡ e r i m e n t ) 2 x i exp ⁡ e r i m e n t ; x i = f M T   ; r C ; L 50 ; (4-1) where N e x p represents the number of experiments which are considered for estimating the parameters (exp = 16 in this work). The set of parameters reported in Table 2 has been obtained.

The two mass transfer related parameters, k L  and a eff , fall in the typical ranges observed in the literature. The estimated k L  is in the range for absorption processes (Last and Stichlmair, 2002). Furthermore, the interfacial surface area a eff for stirred bubble tanks varies usually between 20 and 200 m²/m³ (Levenspiel, 1999b). The values of the two parameters in the nucleation and growth rate, A and α , must be interpreted as empirical, while the corresponding values reported in the literature have a degree of uncertainty. Therefore, we consider the values reported in the table accurate enough and acceptable (Nehrke et al., 2007).

Parity plots, shown in Figure 2, are used to illustrate the quantitative agreement between modelling and experimental results. They illustrate the experimental values (x-axis) of f M T   (Figure 2A), r C (Figure 2B) and the volume based average particle size L 50 (Figure 2C) against the output of the model (y-axis). In case the experimental results and the model were in perfect agreement, the points would lie on the solid line with an angle of 45°. The dashed and the dotted line enclose the area with a 10% and a 25% relative deviation. Note, that f M T , r c and L 50 are obtained through three independent measurements. In case a point falls above the solid line, the model overestimates the corresponding quantity, while on the other hand, points falling below the solid line are underestimated by the model.

The experimental and modelling results for f M T show a very high quantitative agreement—as visualized by Figure 2A. Many points fall on the solid line, none of them exhibits a MRD exceeding 10%. This holds true for the whole range, from low (5.4 mmol min⁻¹ L⁻¹) to high (11.5 mmol min⁻¹ L⁻¹) absorption rates.

From an experimental viewpoint, the rate of CO₂ absorption is the quantity exhibiting the highest accuracy, since it is measured directly online with a previously calibrated mass spectrometer. The empirical parameters k L  and a eff are usually a function of the properties of the aqueous solution, the reactor geometry and the mixing of the aqueous solution (Laakkonen, 2005). The effect of mixing is not considered by the modelling equations. The stirring speed was kept constant in the experimental campaign. Despite the variation of the inlet gas flowrate between 12.8 and 18.5 mmol min⁻¹, the effects on the mass transfer parameters seem to be small.

Moreover, Figure 2B shows that the model tends to underestimate the precipitation rate in many cases by 25% or more, possibly due to the fact that the model describes growth and nucleation in a rather empirical way, for instance without accounting for the effect of the solution composition on the growth rate; these effects might be important, as described in literature (Jung et al., 2000). Despite its relative simplicity the model manages to give a quantitative description of the precipitation rate, at least of the right order of magnitude.

The third quantity, which is measured independently, i.e.,   L 50 , exhibits for experiments f, n, p, t, k, o, and d a very good quantitative agreement, while the other experimental points exhibit a deviation of mostly below 25% (Figure 2C). This good quantitative agreement is remarkable, since certain mechanisms, such as agglomeration, which affect the particle size in the experimental campaign (Tiefenthaler and Mazzotti, 2021), are not considered by the model.

The values of MRD obtained for each measured variable support the observations made based on Figure 2. The quantitative agreement between experimental and modelling results is the highest for f M T , which exhibits a MRD of only 2.3%, followed by L 50 with a MRD of 16%, and finally by r c with a MRD of 27%.

It is worth noting that we have achieved a good accuracy between simulations and experiments using a primary nucleation rate expression. This might appear to be in contradiction with the fact that the model describes a continuous crystallization process, where new CaCO₃ crystals are formed in the presence of already existing CaCO₃ crystals, i.e., under conditions where secondary nucleation should be favored. We have tried to describe the experimental results using a secondary nucleation rate equation, with the specific feature that the nucleation rate is proportional to a moment of the crystal size distribution [the third moment, as in a recent paper (Chen et al., 1997)], but we have failed, i.e., always obtaining a worse fit than what shown in Figure 2 (results not presented here for brevity). We have rationalized this outcome by considering that primary particles are very small, thus making secondary nucleation by attrition, i.e., the dominant mechanism, unlikely, and that supersaturation is very high, thus favoring primary nucleation.

The model allows to compute process performance for any set of operating conditions, and thus helps to gain additional insights, far beyond the experimental evidence. For this reason the model was evaluated in hundreds of operating points (about 400), covering calcium concentrations from 0.2 to 0.35 molal, liquid feed rates from 9 to 36 g per minute and gas feed rates from 12.9 to 24 mmol per minute. The temperature, the CO₂ inlet concentration, and the NH₄NO₃ concentration in the feed solution were kept constant at 25 °C, 25% CO₂ and 0.7 molal, respectively.