Thermodynamic Modeling of Multicomponent Rare Earth Nitrates Aqueous Systems

Solubility constants of rare earth (RE) nitrates crystalline hydrates are determined in a wide temperature range (−30 to 120°C), salts solubilities and phase diagrams of water–RE nitrate systems are calculated. For multicomponent (n > 5) solutions of RE nitrates the assessment of solution properties as well as phase diagrams are shown to be feasible within experimental uncertainty. In case of mixtures of RE nitrates with similar hydrodynamic radii of ions, the parameters of RE1–RE2 interparticle interaction can be ignored without losing accuracy of thermodynamic modeling.


INTRODUCTION
The production and consumption of rare earth (RE) elements in different areas of science and technology are key economic indicators of industrial countries. High-purity RE elements are of particular interest due to their active use as catalysts in oil refinement, luminescence activators, components of hightemperature superconductors, permanent magnets, laser crystals, etc. The available production capacities are nevertheless gradually becoming insufficient to satisfy the growing demand for high-purity RE elements. The technology of liquid-phase extraction from acidic water-organic solutions with subsequent deposition of the product from the raffinate or aqueous re-extract is currently mostly used [1]. Yet, the conditions for these multistep processes are generally adjusted empirically by obtaining huge experimental datasets for certain stocks, that is obviously quite challenging. Thermodynamic modeling can be a reasonable alternative which allows to reduce considerably the time and labor for optimization of the separation process. The aim of this work was therefore to construct a thermodynamic model for describing the properties of phases and their equilibria in multicomponent systems of water and RE nitrates in a wide range of temperatures.
Binary systems (aqueous solutions of nitrates of particular RE nitrate) have been extensively investigated in a number of studies [2][3][4][5][6][7][8][9][10] where the osmotic coefficient and solvent activity were mostly determined isopiestically (relative to potassium and cal-cium chlorides) at 25°C. For yttrium, cerium, samarium, europium, gadolinium, terbium, dysprosium, holmium, thulium, and lutecium nitrates, the operating range of salt molality reached saturated solutions (~4 to 5 mol/kg). For lanthanum, praseodymium, neodymium, erbium, and ytterbium, even supersaturated mixtures we studied. The calorimetric properties of aqueous lanthanide nitrate solutions were considered in [11,12]. Spedding et al. [11] measured the heats of dilution and then proposed the dependences of components partial enthalpies vs. solution composition and ionic radii of the RE element. The most comprehensive review of the experimental data on the solubility of RE nitrates in a wide range of temperatures (−25 to 120°C) and the formation of their crystalline hydrates was presented in [13] (Tables 1 and 2). The authors of [14] performed differential scanning calorimetry (DSC) measurements of the liquidus temperature for aqueous solutions of yttrium nitrate, complementing existing data on the solid-liquid equilibria in this system.
The ice liquidus was additionally studied in [15][16][17], as well as the boiling temperatures were determined for the systems water-nitrates of lanthanum, cerium, praseodymium, neodymium, samarium, and europium. The existence range of crystalline hydrates were also examined and compared to early results. In [18], it was shown that experimental investigations of the solubility of RE nitrates penta-and tetrahydrates are not always possible, since they are metastable towards the hexahydrate. Using these data, phase dia-

CHEMICAL THERMODYNAMICS AND THERMOCHEMISTRY
grams of water-RE nitrate systems (RE = La, Ce, Pr, Nd, Sm, and Eu) were plotted from the eutectic temperature to the boiling point. A summary of the literature experimental data for the H 2 O-RE(NO 3 ) 3 binary systems is presented in Table 2.
Published efforts of thermodynamical modeling of H 2 O-RE(NO 3 ) 3 systems where RE = Y, La, Ce, Pr, Nd, Sm, Eu, Gd, Tb, Dy, Ho, Er, Tm, Yb, and Lu, have been based mainly on vapor-liquid equilibria and employed the Pitzer model [20]. The authors of [2][3][4] managed to describe their results with acceptable accuracy only for a narrow range of concentrations with a maximum at 1.9 mol/kg. Employing a modified Pitzer model allowed to expand it. Pérez-Villaseñor et al. [19] succeeded in describing the properties of H 2 O-RE(NO 3 ) 3 systems up to a salt content of 6 mol/kg with an accuracy better than 1%. In [20], a five-parameter version of the Pitzer model was used to calculate not only the osmotic coefficients above the saturation point but the solubility of salt at 25°C as well. The molar Gibbs energies of RE(NO 3 ) 3 ⋅nH 2 O crystalline hydrate formation, where n = 5 and 6, were estimated simultaneously. The complete dissociation of RE nitrates in water without the formation of intermediate species or associates was assumed in [2,19,20].
Dealing with the same, Chatterjee et al. [21] wrongly assumed that the Pitzer model describes osmotic coefficients well but not the mean ionic activity coefficients at ionic strengths greater than 4 mol/kg. It was nevertheless shown in [22] that an error was made in the calculations of [21] (the dependence of parameter C cx of the Pitzer model on concentration was not considered). Introducing it allows to describe the activity coefficients up to an ionic strength of 20 mol/kg. In addition, a correlation between the Pitzer model parameters and the ionic radius of the RE element was revealed in [21] that can be used to estimate unknown parameters of the H 2 O-Pm(NO 3 ) 3 system.
Modeling RE solutions is not limited to the Pitzer model. For the comparison, Chatterjee et al. [21] employed an expanded version of the semiempirical Bromley model [23]. The electrolyte generalized local composition model (eGLCM) developed at the Laboratory of Chemical Thermodynamics of the Chemistry Department of Lomonosov Moscow State University, which accurately describes binary and multicomponent systems in a wide range of temperatures and concentrations, was used in [24].
Experimental data on ternary systems water-RE nitrate 1-RE nitrate 2 are nonsystematic (Table 3). Mutual solubilities of several RE nitrates were given in [13]  Thermodynamic properties of ternary systems containing a nitrate anion and two RE cations were studied in [25][26][27]: 3 . The isopiestic osmotic coefficients were determined experimentally at 25°C in a wide range of concentrations, from infinitely dilute solutions to saturated. A similar situation was also observed for quaternary and quinary systems of water-RE nitrates with Y 3+ , La 3+ , Pr 3+ , Nd 3+ , and Er 3+ cations [27][28][29] ( Table 4). It was noted in [25] that the absence of interaction between salts makes the empirical Zdanovskii rule [30][31][32][33] applicable with quite high accuracy. The authors of [25][26][27][28][29] used a seven- parameter Pitzer model with a modified parameter C for more rigorous modeling. Investigating the properties of multicomponent aqueous solutions of RE nitrates is therefore generally limited to a temperature of 25°С (with the exception of binary systems). The main type of data are osmotic coefficients, salts solubility, and (rarely) thermochemical properties. Thermodynamic modeling is generally performed at one temperature only (25°С) using different versions of the Pitzer model. CALCULATION PROCEDURE eGLCM model [24] was used for thermodynamic modeling of the properties of aqueous solutions of RE nitrates and phase equilibria in the considered systems. It has proven its reliability with respect to extraction systems, as it allowed to calculate the binary diagram of the water-tributyl phosphate system (the extracting agent), which was unsuccessful with any other models [34].

THERMODYNAMIC MODEL OF THE LIQUID PHASE
The expression for the molar excess Gibbs energy in the eGLCM model contains three terms: where is the contribution of long-range interactions, is the contribution of middle-range interactions, and is the contribution of short-range interactions. A symmetric reference state was used to normalize the properties, and the concentration was in mole fractions. Activity coefficient when for any species in a solution (such state is hypothetical for ions). The dissociation of electrolytes in a solution was assumed to be complete.
The mole fraction of the kth species can be calculated as where summation is made over all species in the solution.
In the eGLCM model, the Pitzer-Debye-Hückel equation is used to consider long-range interactions in a symmetric reference state for each component: where is the ionic strength expressed in the scale of mole fractions; is the ionic strength in a hypothetical single-component system consisting of an i-th component; and is a parameter of the closest approach of ions with diameter = 5.4671 × 10 -10 m:   The density and dielectric permittivity of a dilute solution are normally assumed to be equal to the corresponding values of the solvent; in this study, however, the indicated values were determined using empirical equations (mixing rules) to achieve a better description: where , is the molar volume of The molar weight of a solution is calculated from molar weight of its components: where and are the molar fraction and molar weight of the ith component, respectively.
The contribution from long-range interactions to the activity coefficient is expressed as Middle-range interactions in the eGLCM model are those involving charged species and which are not considered in the Pitzer-Debye-Hückel theory. Corresponding excess Gibbs energy is expressed as where is the interaction parameter between the -th and -th species. It depends on the ionic strength and has the form of a symmetric matrix: and . , since the th and th species are uncharged. Parameter is fitted with the empirical dependence where and are the matrices of interaction between the components, for interaction between uncharged species and −1 for interaction between ions, and . The expression for the corresponding contribution to the activity coefficient has the form where The excess Gibbs energy for short-range interactions in the eGLCM model are written as where where and are structural parameters (the relative surface area and van der Waals volume) of the i-th species; and are temperature-dependent parameters of binary interactions between the i-th and j-th species ( , , The contribution from short-range interactions to the activity coefficient takes the form The properties of pure water were taken from [35,36]; structural parameters and , from [24] (Table 5).

OPTIMIZING MODEL PARAMETERS AND CALCULATING THE PHASE DIAGRAMS OF BINARY AND MULTICOMPONENT SYSTEMS
Parameters of the eGLCM model were optimized using the Levenberg-Marquardt algorithm [37] in the MATLAB ® R2017a programming environment. The sum of squared deviations between experimental and analytical data was minimized. The general form of the target function was where x i exp is the experimental value, x i calc is the value calculated using the model, is the vector of model parameters, and is the number of experimental points. The standard error of regression was estimated as where is the number of experimental points, and is the number of the model parameters.
x a x f x a n = σ = − 

Binary Water-RE Nitrate Systems
The parameters of middle-range interactions b ij and c ij were optimized in [24]   NO are hypothetical activities of RE and nitrate ions in a non-symmetric reference state, and is the activity of water in a symmetric reference state. Since the eGLCM model deals with a symmetric one by default, the following expression was used to determine the activities at infinite dilution ( at ): where is the activity coefficient of the th component in the non-symmetric reference state; is the activity coefficient of the th component in the symmetric reference state; and is the activity coefficient in the symmetric reference state at infinite dilution of component in water (i.e., when ).
The dissociation constants of the salt crystalline hydrates were described by three-parameter empirical dependences in the form where A, B, and C are coefficients determined from least squares, based on deviations of the liquidus points measured experimentally and calculated by the eGLCM model.
The corresponding equilibrium constant for ice was determined using the known stability parameters of pure water: where J/mol and are the enthalpy and temperature of melting of water, and J/(mol K) is the change in the heat capacity of water at the temperature of melting [38].

Multicomponent Systems
According to [13], solid solutions of two or more RE nitrates are precipitated from aqueous solutions, instead of pure crystalline hydrates. Depending on the combination of the nitrates, such solid solutions can be either continuous or limited. The former were considered to be ideal, and the latter were described using the Margules formalism.    The calculations and construction of the phase diagrams sections were done using the TernAPI program developed at Lomonosov Moscow State University's Laboratory of Chemical Thermodynamics [39].

RESULTS AND DISCUSSION
Two-Component Systems The dissociation constants of RE crystalline hydrates were assessed using literature data on their solubilities. Obtained parameters of the temperature dependences are given in Table 7. The values of the constants at 25°С coincide with those listed in [24]. Binary systems containing Ho, Tm, and Lu were not considered, since there were no measurements for them at temperatures other than 25°С. Further in the text, we consider only those investigated multicomponent systems that were studied most comprehensively (i.e., those containing Y, La, Pr, Nd, and Sm [13][14][15][16]) (Figs. 1a-1e).
Note that we were able to describe the available experimental data by introducing a temperature dependence only for dissociation constants of RE crystalline hydrates while keeping parameters b ij and c ij of the eGLCM model constant. Since experimental studies were performed only for the osmotic coefficients, heats of dilution, and heat capacities of the solutions at 25°С, we considered it unnesessary to make b ij and c ij temperature dependent, which also increased the number of optimized parameters considerably.
Looking on the ice liquidus, one can notice that the calculated curve is a bit steeper than the experimental one in almost all cases. Any efforts to introduce the temperature dependence for b ij and c ij made no noticeable improvements here either. We associate this deviation with the features of differential thermal analysis (DTA) used by the authors of [15][16][17]. The incomplete deconvolution of the DTA signal can lead in particular to the increased results for the melting point, mostly due to the delay in the recorded system response. Additional experiments are nevertheless required to clarify the reasons.
Note also that there are almost no solubility data for lower RE crystalline hydrates, probably because of the high equilibrium temperatures. In cases where the metastable solubility is known, the eGLCM model with the optimized parameters describes it very well (dashed lines in Figs. 1a-1e).

Three-Component Systems
A thermodynamic model of ternary systems containing water and two RE nitrates (Y, La, Pr, Nd, Er, and Sm) was constructed on the basis of binary subsystems with earlier determined parameters b ij and c ij .

Parameters of -interparticle interaction
were optimized using experimental data on the water activity. They proved to be statistically significant in 3 , and allowed to describe better the properties of solutions (Figs. 2a and 2b). Their necessity is likely due to the strong difference between the hydrodynamic radii of RE 3+ ions (the average distance between an ion and the center of the nearest water molecule) [40]. Table 8 lists the optimized values of the parameters for all of the studied systems. In all cases, c 23 = 0, which can be considered as the independence of interparticle interactions on a solution's ionic strength. Differences between the calculated and the experimentally measured water activities are listed ln K ln C T in Table 9. The root-mean-square error does not exceed 0.25%, and in most cases it does not go beyond the limit of 0.09%.
The mutual solubilities of RE nitrates were predicted using optimized values of binary and ternary parameters of interparticle interactions. Since there are no data for phase equilibria between the lower RE crystalline hydrates, in this work we confined to constructing the phase diagrams sections only up tõ 70 wt % salts. Results are presented in Fig. 2 (Fig. 2a) in combination with the eGLCM model shows its high predictive ability.

Multicomponent Systems
For calculations in multicomponent systems we used the model parameters determined at the previous step, assuming there were no higher order interactions in the solutions of several salts or they could be ignored relative to others. Unfortunately, there is no informa-tion on liquid-solid equilibria in multicomponent aqueous mixtures of RE nitrates. Still the water activity in quaternary and quinary systems is quite predictable,   Table 9. Resulting accuracy of describing the activity of water in RE nitrate solutions using the eGLCM model (δ is the root-mean-square error) confirming the quality of constructed eGLCM model. As was expected, the accuracy of describing the properties of a system rises along with its dimension. The comparison of calculation results with the experimental data for the studied systems are given in Table 9. The rootmean-square error does not exceed 0.09%. CONCLUSIONS Constructed eGLCM model for multicomponent RE nitrates solutions describes highly accurate the different types of phase equilibria as well as the properties of aqueous solutions of RE nitrates from −30 to 120°С within the experimental error using a fairly small set of parameters. Accounting the binary interactions only allows to extend the eGLCM model's range of applicability to systems containing three or more RE elements. In case of the mixtures of RE nitrates with similar sizes of ions in the solution, parameters of -interaction can be ignored without losing in accuracy of modeling. FUNDING This work was performed as a part of State Assignment no. 121031300039-1, "Chemical Thermodynamics and Theoretical Materials Science" with partial financial support of RFBR (grant no. 18-29-24167).

OPEN ACCESS
This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.