Studying the speciation and mineral-fluid partitioning of the rare earth elements (REE) allows us to delineate the key processes responsible for the formation of economic REE mineral deposits in natural systems. Hydrothermal REE-bearing calcite is typically hosted in carbonatites and alkaline rocks, such as the giant Bayan Obo REE deposit in China and potential REE deposits such as Bear Lodge, WY. The compositions of these hydrothermal veins yield valuable information regarding pressure , temperature , salinity, and other physicochemical conditions under which the REE can be fractionated and concentrated in crustal fluids. This study presents numerical simulation results of the speciation of REE in aqueous NaCl-H2O-CO2-bearing hydrothermal fluids and a new partitioning model between calcite and fluids at different -- conditions. Results show that, in a high CO2 and low salinity system, bicarbonate/carbonate are the main transporting ligands for the REE, but predominance shifts to chloride complexes in systems with high CO2 and high salinity. Hydroxyl REE complexes may be important for the solubility and transport of the REE in alkaline fluids. These numerical predictions allow us to make quantitative interpretations of hydrothermal processes in REE mineral deposits, particularly in carbonatites, and show where future experimental work will be essential in improving our modeling capabilities for these ore-forming processes.

1. Introduction

The economic significance of the REE is expected to increase in the coming years due to their varied uses in high technology and green industries [1, 2]. The REE commonly occur in trace concentrations in the earth’s crust [3], while economic REE deposits typically are associated with carbonatite and alkaline/peralkaline igneous systems at intraplate tectonic settings. Prominent examples include Bayan Obo in China [2], Strange Lake and Nechalacho in Canada [47], Mountain Pass, CA [8], and Bear Lodge, WY, in the USA [9].

Igneous processes such as partial melting, crustal assimilation, melt immiscibility, and fractional crystallization, along with crustal metasomatism, can lead to the formation of REE-enriched carbonatite, alkaline, and peralkaline melts [2, 1015]. Numerous studies have also indicated the importance of hydrothermal processes for the mobilization and concentration of the REE in the late magmatic evolution stages of carbonatite and alkaline/peralkaline systems [4, 5, 1619]. Secondary REE enrichment can be significant in the late stages of pluton emplacement, as evidenced by the Strange Lake REE-Zr-Nb deposit in Canada, where metasomatic processes led to hydrothermal mobilization and mineralization of the REE at the periphery of the pluton [4, 5, 16, 17]. Similar metasomatic processes were observed in the Tamazeght alkaline HFSE/REE-enriched pluton in Morocco [20] and other deposits [2, 4, 79].

Quartz-hosted fluid inclusions are generally used within mineralized carbonatites and peralkaline igneous complexes to determine the salinity, metal concentrations, and temperature of such ore-forming fluids [4, 16, 2126]. Another approach that can be used to understand the evolution of fluids in REE mineral deposits is the study of fluid-mineral trace element partitioning. Calcite and fluorite are common gangue minerals, and trace elements are more compatible in these minerals than in quartz. Calcite has proven to be a useful geochemical tracer for the behavior of the REE in different types of ore-forming fluids and also for determining fluid reservoir properties [27]. The partitioning of the REE results in characteristic patterns which can be used in evaluating ore-forming processes in hydrothermal mineral deposits [27]. The REE are compatible within the calcite structure, where substitution of REE3+ for Ca2+ can be achieved due to their similar ionic radii; this has also been observed in fluorite [2730].

Many studies have analyzed REE and other trace elements in natural calcite (Table 1), but interpreting these signatures is difficult because it requires an understanding of the interplay of REE speciation in hydrothermal fluids, the solubility of the minerals, and the incorporation mechanisms of REE in these minerals as a function of pressure , temperature , and fluid compositions . Recent experiments enable us to more accurately predict the formation of REE complexes with the major ligands Cl, F, and at hydrothermal conditions [3335]. However, there is a lack of thermodynamic data at elevated temperatures for most REE hydroxyl and carbonate/bicarbonate complexes, which are mostly based on theoretical extrapolations from low temperature [3639]. Because of the complexity of natural systems, numerical modeling and experimental studies provide the means to explore the link between the compositional changes of these fluids and the trace element signatures of associated gangue minerals such as calcite.

Experimental studies of REE partitioning between fluid-calcite [28, 30, 40] and fluid-fluorite [29] have been carried out at temperatures <100°C. In the study of van Hinsberg et al. [29], the partitioning behavior of the REE between fluid-fluorite yielded predictable results according to the REE ionic radii, and their experimental data were fit to the lattice strain model [41, 42]. This model was later applied to the partitioning of REE between fluid-calcite for a set of selected REE by Voigt et al. [30]. Rimstidt et al. [28] used available distribution coefficient data from the literature to evaluate trace element partitioning between fluid-carbonates at temperatures up to 100°C. These data were fit to a semiempirical equilibrium partitioning model [43]. In their study, they observed a similar behavior of the REE with ionic radii. While these experimental data will be useful to develop future thermodynamic models of trace element partitioning, they are not currently applicable to model natural system using available numerical modeling codes. In contrast, Curti et al. [40] studied the thermodynamics of Eu uptake into calcite under varying pH and pCO2 conditions at room temperature and found that the partitioning of Eu could be predicted using equilibrium thermodynamics and different types of solid solution models. In their study, they used the dual thermodynamic approach, which is based on Gibbs energy minimization and implemented in the GEM-Selektor code package [44, 45].

In this study, we present numerical simulations for the speciation and partitioning of REE between calcite and crustal fluids at hydrothermal temperatures (~100–400°C) and pressures up to 1 kbar. The simulated conditions aim to approximate the chemistry of simple aqueous-carbonic fluids, as well as their driving mechanisms for REE mobilization associated with REE deposit formation in carbonatites and alkaline systems. We simulate the REE speciation in CO2-NaCl-poor near-surface crustal fluids and in deeper CO2-NaCl-rich fluids to approximate conditions relevant to important mineralization stages in carbonatites, such as the Bayan Obo REE deposit. The simulated conditions aim to delineate the effects of changing temperature, pH, and salinity on the REE speciation and partitioning between calcite-fluid. Additionally, we briefly explore current methods for implementing fluid-mineral partitioning data into the GEM-Selektor code package [45] and define where the model may need more experimental data for accurately predicting the partitioning of REE in natural systems. To our knowledge, the simulation of the REE partitioning into calcite at hydrothermal conditions has not been undertaken, and here we apply it to ore-forming processes for the first time.

2. Theoretical Background

2.1. The Lattice Strain Model

Semiempirical partitioning models, such as the lattice strain model [41, 42], aim to determine a set of partition coefficients for ion where, for example, the partitioning of a trace ion between fluid and mineral is described byand the partitioning of a major cation bywhere [REE3+] and [Ca2+] indicate the concentrations of the respective ions in aqueous solution and in the calcite crystal. More details on how to relate partition coefficients and equilibrium constants can be found in Wood and Blundy [42]. In the rhombohedral crystal structure of calcite, Ca2+ has a 6-fold coordination, which corresponds to a Shannon radius of 1.00 Å, whereas the REE3+ have ionic radii varying between 1.032 Å (La) and 0.861 Å (Lu). The lattice strain model is based on the Brice equation [41], and it considers the partitioning of an element as a function of its ionic radius, charge, and the elastic properties of the mineral lattice site (i.e., Young’s modulus). When the REE partition from an aqueous solution into a mineral like calcite, strain occurs on the crystal lattice because of this slight size and charge mismatch between REE3+ and Ca2+.

Initial low temperature partitioning data by van Hinsberg et al. [29], Voigt et al. [30], and Rimstidt et al. [28] indicate that the lattice strain model can be successfully applied to the partitioning of REE between aqueous fluids and hydrothermal minerals, including calcite and fluorite. These studies have shown that partition coefficients of REE can be related to their ionic radii and fit to the lattice strain model according to

where is the partition coefficient for a lattice site with a cation of ideal ionic radius and charge, is the partition coefficient of a trace ion , is the Avogadro constant, is Young’s modulus, and and are the ionic radii of the crystallographic site (i.e., Ca2+ for calcite) and the substituting ion (i.e., REE3+), respectively. Young’s elastic modulus can be verified by experimental methods.

An experimentally based lattice strain fit can then be used to determine the REE3+ concentrations of an aqueous fluid at equilibrium with a mineral and even predict partition coefficients for elements with the same charge, which were not determined experimentally. A drawback of this model is that currently it cannot be implemented in numerical modeling codes for modeling equilibrium thermodynamics. The partition coefficients will have to be determined experimentally as a function of temperature and fluid chemistry to obtain the necessary parameters for fitting REE data to the lattice strain model. The Brice equation [41] will either need to be modified or coded for implementing the various parameters of (3).

2.2. Multicomponent Solid Solution Models

To simulate the partitioning of REE3+ between fluid and mineral, the substitution mechanisms of the ion in the crystal lattice and the energy involved need to be determined. There are two types of possible models that can be used, including the “forward modeling” and the “inverse modeling” approaches, which define what type of thermodynamic data have to be implemented in a thermodynamic database.

The “forward modeling” approach uses knowledge of either the solubility products or Gibbs energy of formation of known REE mineral endmembers (i.e., REE2(CO3)3, NaREE(CO3)2, REEOHCO3, and REE(OH)3), and it can be readily implemented in a numerical modeling code if thermodynamic data are available for them. This approach assumes an ideal solid solution model for calcite and the corresponding endmembers, since trace elements can be assumed to occur in concentrations small enough where activities of cations in the solid solutions can be assumed to be equal to mole fractions. There are several possible substitution mechanisms in the structure of calcite responsible for the partitioning of the REE3+:where, for example, in Reaction (4), 2REE3+ ions replace 3Ca2+ and create a site vacancy, while Reaction (5) corresponds to a coupled substitution. Currently, only high temperature thermodynamic data are available for the REE(OH)3(s) endmembers, and therefore only substitution mechanism (7) was tested in the present numerical modeling study.

“Inverse modeling” is an alternative approach based on experimental partitioning data and the dual thermodynamic approach [44]. Using this method, it is possible to determine the Gibbs free energy of “fictive” endmembers and build a multicomponent solid solution model that may involve more complex stoichiometries. The advantage of this method is the possibility of analyzing several experimental datasets at different pH and salinities. This method can also be used to determine possible substitution mechanisms to be implemented in the model. This approach was used to interpret room temperature experimental data for the uptake of Eu into calcite [40]. However, implementation of this model and determination of the substitution mechanisms of REE in calcite will require additional experimental data to be applicable at hydrothermal conditions in ore-forming systems. These calculations also need to be complemented by spectroscopic measurements to verify that the calculated thermodynamic model can be tied to observed substitution mechanisms in the crystal structure [46].

2.3. Numerical Modeling Codes: LMA and GEM Methods

Current numerical modeling codes can implement various solid solution models. Law of mass action (LMA) codes such as PHREEQC [47], TOUGHREACT [48], and Geochemist’s Workbench® [49] use databases involving equilibrium constants for a set of mineral dissolution reactions and master species, and they utilize the Newton-Raphson method to solve for chemical equilibria. These codes are robust and fast but not adequate to simulate the partitioning of REE between fluids and minerals since they are limited to simple binary and/or tertiary solid solution models. Gibbs energy minimization (GEM) codes such as GEM-Selektor [45] are more ideal to solve complex multicomponent and multiphase systems, which can involve any number of solid solution endmembers and mixing on different crystallographic sites. The GEM-Selektor code package includes several types of nonideal mixing models [44, 5052]. This can be partly related to the larger thermodynamic input data (i.e., Gibbs energy , enthalpy , entropy , heat capacity functions , and volume ) of their thermodynamic databases and to their numerical solver [44, 45] in comparison to LMA codes. Both code types have their advantages and disadvantages, and the recent development of the Reaktoro code package [53, 54], which combines the strengths of both LMA and GEM methods, has seen many promising improvements for simulating complex chemical systems and reactive transport problems. The present study focuses on the GEM method.

3. Methods

3.1. Thermodynamic Dataset and Activity Models

The thermodynamic dataset used in this study is the MINES database (http://tdb.mines.edu), which can be downloaded and managed with the GEM-Selektor code package [45]. A detailed list of minerals, aqueous species, and gases used from this database for the Ca-C-Na-Cl-O-H-REE components system can be found in Table 2. This dataset comprises mineral data from Holland and Powell [55] and aqueous species from SUPCRT92 [36, 57, 63, 64], with the recently updated internally consistent dataset of Miron et al. [58] using the GEMSFITS code package [65]. Experimental data for REE chloride aqueous species were implemented using data from the study of Migdisov et al. [35], the HCl dissociation constants of Tagirov et al. [59], and the properties of REE(OH)3(s) compiled in the study of Navrotsky et al. [56]. Data for aqueous REE carbonate and bicarbonate species were included using the theoretical predictions of Haas et al. [36], which are currently the only available high temperature thermodynamic data for these species. We acknowledge the high level of uncertainty associated with predicting the speciation of REE hydroxyl complexes using the data of Haas et al. [36], but the lack of other internally consistent datasets for these species necessitates its use for the current modeling example until more experimental data become available. For an in-depth review on available experimental data for REE aqueous species and limitations, the reader is referred to the recent study by Migdisov et al. [39].

The standard state adopted in this study was unit activity for pure H2O and for the aqueous species in a hypothetical one-molal standard solution referenced at infinite dilution at any pressure and temperature. The standard state for minerals was that of a pure phase, assuming Raoultian behavior where activity is equal to mole fraction for an ideal solid solution. All activity models and equations of state were calculated using the TSolMod library class implemented in GEM-Selektor [50]. The activity model used for aqueous species was the extended Debye-Hückel according to Robinson and Stokes [66], using the extended Debye-Hückel parameters of Helgeson et al. [67], assuming NaCl as the background electrolyte. The equation of state used for real gases (CO2-CH4-CO-H2O) was the Peng-Robinson-Stryjek-Vera (PRSV) model [62]. The solid solution model adopted for simulating the partitioning of REE into calcite was an ideal solid solution between CaCO3 and the REE hydroxide endmembers. This model includes correction for the configurational entropy resulting from the assumption of random mixing on the crystallographic Ca2+ site of calcite.

3.2. Simulation Setup

Numerical simulations of calcite-fluid interaction were carried out using the GEM-Selektor code package (http://gems.web.psi.ch). In these thermodynamic calculations, we aim to determine the speciation of REE in aqueous and carbonic fluid mixtures and their partitioning behavior between fluid-calcite as a function of temperature, salinity, and pH. The simulated models assume various scenarios of hydrothermal calcite vein formation to study the chemistry of the fluids responsible for REE transport in ore-forming systems, with a particular focus on carbonatites. The model uses the components Ca-C-Na-Cl-O-H-REE, that is, the chloride, hydroxyl, and carbonate/bicarbonate REE complexes. Other fluid components including F, S, P and their associated REE-bearing minerals (e.g., monazite-(Ce) and bastnäsite-(Ce)), as observed in Bayan Obo [68] were not included in the present simulations. The main reason these were omitted was to allow us to better distinguish the processes responsible for REE fractionation solely associated with the solubility of calcite in aqueous-carbonic fluids. Omission of these components was also necessary in order to reduce the overall number of variables in the model. The starting fluid compositions are reported in Table 3, and these were modeled from generalized fluid inclusion data taken from various natural REE-enriched carbonatite and alkaline systems (Table 4).

Three speciation models with different initial fluid compositions were assumed (Models Ia–Ic), along with a fluid-calcite partitioning model (Model II). The three main subsystems were modeled with variable temperature (150 and 350°C), initial XCO2 (0.01–0.1), and salinity (0.5 and 20 wt.% NaCl). The initial REE concentration was set to 100 ppm, within the upper range of concentrations retrieved from natural fluid inclusion data [10, 23, 25, 39, 69, 70], and used in previous numerical modeling studies [71]. Variable pH values were fixed in the speciation model by titrating different amounts of NaOH into the starting fluids to reach saturation with calcite. In these models, the mobility of the REE was fixed by the stability of solid REE(OH)3.

Model Ia (CO2-NaCl-poor system) was conducted at temperatures of 150 and 350°C and at saturated water vapor saturation pressure . Speciation Model Ia represents a low pressure system, which is a simplified proxy for near-surface evolved crustal fluids with low XCO2 and salinity. Such a fluid can be considered as an analogue for a magmatic fluid, which later evolved due to CO2 loss from immiscibility, neutralization/fluid-rock interaction, and/or dilution from meteoric water input. These different processes can account for the low salinity and XCO2 of the fluid in shallow ore-forming environment.

Models Ib and Ic were simulated using a simplified fluid inclusion dataset (i.e., omitting F, S, and P) from the Bayan Obo REE deposit. Model Ib represents a NaCl-poor, CO2-rich system, relating to the early (monazite) mineralization stage for this deposit, while Model Ic represents a CO2-NaCl-rich composition related to the main stage mineralization event at Bayan Obo [70] (Tables 3 and 4). Model Ib (CO2-rich and NaCl-poor system) was conducted at temperatures of 150 and 350°C, a pressure 1 kbar, and a salinity of 0.5 wt.% NaCl. Model Ic (CO2-rich and NaCl-rich system) was simulated at the same conditions but with a salinity of 20 wt% NaCl.

Model II is a REE partitioning model between fluid-calcite as a function of temperature from 25 to 400°C, with conditions similar to Models Ib and Ic. This model aims to approximate the REE trace element concentrations in calcite expected in a hydrothermal calcite vein from proximal to distal from a pluton heat source. The initial REE concentration was set to 100 ppm, with salinities of 0.5, 5, and 20 wt.% NaCl and XCO2 of 0.1. The partitioning of REE was modeled by an ideal solid solution model between calcite and REE(OH)3(s), assuming the substitution mechanism of Reaction (7).

4. REE Speciation Model

4.1. Model Ia (CO2-NaCl-Poor): Shallow System

Speciation Model Ia simulates the formation of a calcite vein in a near-surface (low pressure), evolved crustal fluid; results are displayed in Figure 1. These simulations are modeled with low initial CO2 and NaCl concentrations (Table 3). The REE chloride species are generally dominant at low pH up to ~4.25 to 5.25 and hydroxyl REE species at higher pH values. This transition from REE chloride to REE hydroxyl dominant species is shifted to lower pH values with increasing atomic number (i.e., Ce and Nd, followed by Er). The REE carbonate and bicarbonate species are minor at all simulated conditions of this model. With increasing temperature, the activity of REE3+ becomes insignificant over the activity of the REE chloride species.

In the simulations at 150°C, the dominant Ce species is to pH values of ~5.25 (Figure 1(a)). Above this pH, Ce-hydroxyl species become dominant in order of and then with increasing pH. Simulations of Nd-bearing species at 150°C show similar trends to those of the Ce-bearing species. From low pH values to ~4.5, is dominant (Figure 1(c)). Above these pH values, the hydroxyl species become dominant in order of Nd(OH)2+, Nd, , and with increasing pH. The dominant Er-bearing species at 150°C (Figure 1(e)) is ErCl2+ to pH values of ~4.25 followed by the hydroxyl species in order of Er(OH)2+, , and finally , displaying a large predominance field at pH values of >6.5.

In contrast to the simulations at 150°C, the simulated Ce-hydroxyl species at 350°C display an increase in the stability of over (Figure 1(b)). Additionally, with increasing pH, hydroxyl species evolve in predominance from to to . The dominant Nd-bearing species at 350°C are NdCl+2 and to pH values of ~4.75 (Figure 1(d)) followed by at pH of ~4.75 to ~6.75, and to higher pH, the latter showing a large field of predominance. The dominant Er-bearing species at low pH is , rather than ErCl2+, to a pH of 4.5 (Figure 1(f)), followed by which is dominant to pH values of ~5.75, and at higher pH. This is in contrast to the 150°C simulations, where Er(OH)2+ is the first Er-hydroxyl species to be dominant. The dominance of hydroxyl-bearing species at pH > 6 is due to the increased availability of free OH ions at elevated pH and overall lack of enough available carbonic species due to the low XCO2 in the simulated fluid and the formation of calcite.

Calcite is stable in the 150°C models from a pH of 6.02 to >10, whereas, in the 350°C models, it is stable from a pH of 5.92 to >10 (Figure 1). Solid Ce(OH)3 is stable from a pH of 5.34 to 10 and a pH of 4.81 to 8.83 at 150 and 350°C, respectively (Figures 1(a) and 1(b)). Solid Nd- and Er(OH)3 become stable at lower pH with increasing atomic number (Figures 1(c), 1(d), 1(e), and 1(f)). Solid Nd(OH)3 is stable from a pH of 5.28 to 9.98 and 4.71 to 8.23 at 150 and 350°C, respectively (Figures 1(c) and 1(d)). Solid Er(OH)3 is stable at 150°C from pH 4.28 to 9.70 and from pH 3.60 to 8.30 at 350°C (Figures 1(e) and 1(f)). The smaller stability field of the solid REE hydroxides at 350°C can be related to the stabilization of the species with increasing pH.

4.2. Model Ib (CO2-Rich, NaCl-Poor): Deeper Carbonatite System

Speciation Model Ib simulates conditions of calcite vein formation with an initial fluid composition corresponding to the early mineralization stage of Bayan Obo [70] (Tables 3 and 4). In this model, XCO2 was set to 0.1 and the salinity is low (0.5 wt.% NaCl). As shown in Figure 2, the simulated REE speciation is considerably modified in comparison to Model Ia due to the increased amount of CO2 within the system. The stability of the REE chloride species and REE3+ is generally restricted to pH values < 4 at 150°C, due to the predominance of REE carbonate and bicarbonate species, whereas REE hydroxyl species are dominant at pH values > 6–6.5. This transition from REE3+/REE chloride to REE carbonate/bicarbonate is shifted to lower pH values with increasing atomic number (i.e., Ce and Nd, followed by Er). Simulation results from the 350°C models show that, at this elevated temperature, REE bicarbonate and –carbonate species display very low activities, and only REE chloride and hydroxyl species are dominant.

In the simulations at 150°C, is the dominant REE chloride, and the dominant REE bicarbonate species to a pH of ~6.5 (Figure 2(a)). At higher pH, is dominant to pH values of 9.25, followed by at higher pH. Cerium and Nd formed stronger bicarbonate complexes than , as evidenced by their larger pH range of predominance (Figures 2(a), 2(c), and 2(e)). Simulations of Nd-bearing species at 150°C show similar trends to those of the Ce-bearing species but the major chloride species are to a pH of ~3.75. This is followed by and a small pH range of ~5.75–6 where is stable (Figure 2(c)). Similar to Model Ia, and are stable to high pH. Simulations of Er-bearing species at 150°C showed some differences to Ce in that ErCl2+, rather than , is the dominant Er chloride species to a pH value of 3.75 (Figure 2(e)). The stability of and is only over a small pH range in comparison to Nd and Ce species, and is dominant at pH > 6.25.

The simulations carried out at 350°C display a much larger stability of the REE chloride and hydroxyl species in comparison to the 150°C simulations (Figures 2(b), 2(d), and 2(f)). Bicarbonate and carbonate REE species are unimportant at these conditions. Instead of as the dominant Ce-chloride species at 150°C, CeCl2+ is dominant at 350°C up to pH values of ~4.5, followed by a small pH interval where is stable. Dominance then shifts to at pH of ~5.25 to 6.75. This is followed by the dominance of to higher pH (Figure 2(b)). The overall simulated Nd- and Er-bearing species at 350°C are similar to Ce except that the pH transitions between REE chloride and hydroxyl species are shifted to lower pH values with increasing atomic number of the REE (Figures 2(d) and 2(f)).

Calcite is stable in the 150°C models above a pH value of 5.07, whereas, in the 350°C models, it is stable above a pH value of 5.37 (Figure 2). This shift in the calcite solubility with temperature can be explained by the temperature dependence of the ionization of to form . Solid Ce(OH)3 is stable above pH 6.11 and from a pH of 4.96 to 8.09 at 150 and 350°C, respectively (Figures 2(a) and 2(b)). Solid Nd- and Er(OH)3 become stable at lower pH with increasing atomic number. Solid Nd(OH)3 is stable above a pH of 6.03 and from pH 4.96 to 7.45 at 150 and 350°C, respectively (Figures 2(c) and 2(d)). Solid Er(OH)3 is stable at 150°C at pH values from 4.46 to 9.85 and from pH values of 4.01–7.41 at 350°C (Figures 2(e) and 2(f)). The smaller stability field of the solid REE hydroxides at 350°C can be related to the stabilization of the species with increasing pH, similar to Model Ia.

4.3. Model Ic (CO2-NaCl-Rich): Deeper Carbonatite System

Speciation Model Ic simulates conditions of calcite vein formation with an initial fluid composition corresponding to the main mineralization stage at Bayan Obo [70] (Tables 3 and 4). The results of the simulations of these saline (20 wt.% NaCl) and CO2-rich (XCO2 of 0.1) fluids are shown in Figure 3. In the simulations at 150°C, the REE chloride species are dominant to pH values of ~5.75–6.5, with a considerable increase in dissolved REE activities in comparison to Models Ia and Ib. The onset of solid REE hydroxide stability is linked to the change from the dominance of REE chloride species at low pH to that of REE carbonate or bicarbonate species at 150°C or to REE hydroxyl species at 350°C at higher pH. In the 350°C simulations, this transition occurs at pH of 4.75–5.25, though only hydroxyl-bearing species were important in these models. As in the 150°C simulation, this transition pH at 350°C appears to shift to lower values with increasing atomic number of the REE. The transition from chloride to hydroxyl species predominance is due to the increase in available OH ions with increasing pH, while the overall lowered activities are due to the precipitation of solid REE hydroxides making the REE unavailable to form aqueous complexes. The REE bicarbonate and carbonate species are less important under the simulated conditions at elevated temperature due to the increase in salinity in this model which favors higher activities of the chloride-bearing complexes.

The simulations conducted at 150°C show that is the principal species to pH values of 6.5 (Figure 3(a)). From pH 6.5 to 9.25, is dominant, followed by to higher pH values (Figure 3(a)). Neodymium- and Er-bearing species formed at the simulated conditions are similar to Ce, albeit occurring at slightly different pH; is dominant at pH values up to 6 (Nd) and 5.5 (Er), followed by and to pH values ~6.75 and 6, respectively. Overall the predominant REE species then shifts to the hydroxyl species. From pH of ~6.75 to 8 and 6 to 6.25, and are dominant, respectively, followed by and at higher pH values (Figures 3(c) and 3(e)). The effect of decreasing transition pH with increasing atomic number is evident in these results (Figures 3(a), 3(c), and 3(e)).

At 350°C, Ce-bearing species behaves the same as in the model at 150°C, though with pH transitions shifted to slightly lower values (Figure 3(b)). The Nd- and Er-bearing simulations show a definite change in the speciation results at this higher temperature; REE carbonate species are, as seen in Model Ib, only important in high CO2 and lower temperature (150°C) systems, regardless of salinity. The carbonate species are not dominant for these REE at elevated temperature. Rare earth element hydroxyl species become dominant at this temperature at pH values as low as 4.6 for (Figures 3(d) and 3(f)).

Calcite is stable in the 150°C models at pH above 5.26, whereas, in the 350°C models, it is stable above a pH of 5.65 (Figure 3). Solid Ce(OH)3 is stable at pH above 6.19 and pH values of 5.29 to 7.96 at 150 and 350°C, respectively (Figures 3(a) and 3(b)). Solid Nd(OH)3 is stable at pH values above 6.07 at 150°C and at a pH between 5.19 to 7.29 at 350°C (Figures 3(c) and 3(d)). Solid Er(OH)3 is stable at pH values between 4.89 to 9.84 at 150°C and from 4.19 to 7.26 at 350°C (Figures 3(e) and 3(f)).

5. Calcite-Fluid REE Partitioning (Model II)

5.1. Model IIa: REE Fractionation in Calcite

In this model, we simulate the solubility of calcite and the partitioning of all the REE between fluid-calcite using a solid solution model according to Reaction (7). As shown in Figure 4, we investigate the role of temperature and variable salinities of 0.5, 5, and 20 wt.% NaCl on the mole fraction of precipitated REE in calcite ( calcite).

The overall partitioning behavior of the REE is similar at all simulated conditions, where the heavy (H)REE, including Er, Dy, Yb, and Ho, are preferentially incorporated into calcite at low temperature and the light (L)REE, including La, Ce, and Pr, at higher temperature. An overall decrease with temperature of simulated in calcite is also observed for all the REE, which correlates with an increase of calcite precipitation at higher temperature. The simulated display more significant fractionation at low temperature (<250°C); this becomes less pronounced at elevated temperature (>350°C). This behavior is exemplified by the simulations with 0.5 wt.% NaCl and 5 g of added calcite (Figure 4(a)). In these simulations, Dy and Er are preferentially incorporated into calcite at ~100°C, and the LREE at temperatures >150°C until in calcite does not change significantly for either the LREE or HREE at >250°C.

Comparison of the simulations with 5 g of added calcite (Figures 4(a), 4(c), and 4(e)) versus the simulations with 10 g of added calcite (Figures 4(b), 4(d), and 4(f)) displays an overall shift of in calcite with temperature, which correlates with the onset of calcite saturation. Simulations with 20 wt.% NaCl (Figures 4(e) and 4(f)) illustrate this partitioning behavior, where peaks of and occur at temperatures of ~200°C in the model with 5 g of added calcite, whereas these peaks are shifted to ~150°C in the model with 10 g of added calcite, respectively.

Simulations carried out at 0.5, 5, and 20 wt.% NaCl display an overall decrease in calcite saturation with increasing salinity. This decrease in calcite solubility results in an important shift in the fractionation behavior of HREE versus LREE with increasing temperature. More importantly, the fractionation behavior of the REE becomes more pronounced with increasing salinity. For example, in the simulations with 5 g of added calcite (Figures 4(a), 4(c), and 4(e)), the fractionation between LREE and HREE in calcite is less significant at temperatures of >400°C for the simulations with an initial salinity of 20 wt.% NaCl, whereas this temperature threshold drops to ~250°C for the simulation with 0.5 wt.% NaCl. These observations suggest that salinity and calcite saturation both play a major role in the partitioning behavior of the REE in natural systems.

5.2. Model IIb: REE Speciation as a Function of Temperature

The stability of aqueous REE complexes is simulated in Model IIb for Ce, Nd, and Er as a function of temperature and salinity, where the REE were allowed to fractionate into calcite (Figure 5). Carbonate and bicarbonate are the major ligands for REE speciation at temperatures ranging between <200–250°C in the low salinity system (0.5 wt.% NaCl; Figures 5(a), 5(c), and 5(e)). The stable LREE (i.e., Ce and Nd) species are dominantly bicarbonates (i.e., and ) at the simulated conditions, whereas the HREE display a larger stability for carbonate species (i.e., ). An overall decrease in simulated REE ligand activities is observed with increasing temperature due to an increase in REE-bearing calcite precipitation related to its retrograde solubility. The simulated activities of REE chloride species increase with temperature, level off, and decrease again. These species are relatively unimportant, however, due to the elevated pH conditions at calcite saturation. The simulated REE species are dominantly REE hydroxyl complexes with very low overall REE activities in the aqueous fluid above a temperature of 200°C.

The dominant REE species are different in the simulations with a salinity of 20 wt.% NaCl. The carbonate and bicarbonate species are replaced by REE chlorides (i.e., REECl2+ and ) in the saline system at temperatures up to 300–350°C for Ce and Nd and up to 250°C for Er (Figures 5(b), 5(d), and 5(f)). At elevated temperatures, a higher stability of LREE over the HREE chloride species is observed, where the overall activity of REE decreases sharply between 100 and 400°C. It follows that the REE hydroxyl species play a minor role for the mobility of the REE in the saline system.

6. Implications for REE Partitioning in Natural Systems

6.1. Uncertainties and Assumptions of the Model

The numerical simulations presented in this study are meant to represent several scenarios of calcite vein formation. It must be emphasized that these models assume a simplified closed-system under varying hydrothermal conditions to assess in more detail the chemistry of REE-bearing fluids saturated with calcite. In an effort to better understand the speciation and partitioning of the REE, simplifications of the simulated fluids required the exclusion of some components sometimes found in natural systems, such as F, P, and S. While these exclusions represent a limitation in the way of directly applying these models to certain natural systems, doing so allowed us to selectively assess the effects of REE hydroxides and calcite solubility on the REE partitioning behavior in these fluids, decoupled from other possible controls. For example, F-bearing minerals such as bastnäsite-(Ce) and P-bearing minerals such as monazite-(Ce) have very low solubilities [7275] and will sequester some of the LREE from the fluids.

An uncertainty of the current model resides in the difficulty in determining the accuracy of the thermodynamic properties for REE carbonate, bicarbonate, and hydroxyl aqueous species at high temperature. These data rely on theoretical extrapolations of Haas et al. [36] from ambient condition experiments and have not yet been experimentally verified, except for certain hydroxyl species (see review by Migdisov et al. [39]). Experimental data from Wood et al. [37] for Nd hydroxyl species show that the theoretical predictions of Haas et al. [36] may significantly overestimate the stability of these complexes at elevated temperature. This is in contrast to the study by Pourtier et al. [38], who determined the solubility of monazite-(Nd) from 300 to 800°C at 2 MPa. These experiments better agree with the predictions by Haas et al. [36] at 300°C in the low pH range, whereas their experiments showed better agreement with Wood et al. [37] at high pH. Due to the limited availability of experimental data at high temperature, the present simulations were created using the extrapolations from Haas et al. [36], which permit us to keep an internally consistent dataset of the REE hydroxyl species. Simulations with high salinities and low pH are not likely to be affected by these data since the major REE ligands are predominantly chloride species based on experimentally derived thermodynamic data at hydrothermal conditions [39]. Our results show, however, that REE hydroxyl complexes could contribute significantly to the overall dissolved REE content within ore-forming fluids at near-neutral to high pH. It is therefore imperative to gather more data for these complexes in order to successfully understand and model these species at elevated temperature.

Lastly, thermodynamic data for Lu(OH)3(s) are not available, as it was omitted in the Navrotsky et al. [56] study. Cerium-, Dy-, Er-, and Tm(OH)3(s) thermodynamic data from the aforementioned study were derived from thermochemical data reported in Diakonov et al. (1998) [76]. These thermodynamic properties were calculated from available literature data, and unknown values were found using the correlations given in Bratsch and Lagowski [77]. This highlights another need for more experimental data for REE-bearing minerals and further verification of theoretical predictions.

6.2. Hydrothermal REE Transport and Comparison to Previous Modeling Studies

The REE commonly form complexes with the following ligands in crustal fluids: Cl, , F, OH, , and . The stability of these REE complexes from ambient temperature up to ~150°C can be predicted using Pearson’s rules [78, 79]. These rules are based on the hard/soft ligand theory and state that monovalent ligands, when complexed with REE, have the following order of stability: F > OH > > Cl > Br. Divalent ligands are predicted to have the following order of stability: > > . Based on ionic REE properties, the stability of REE complexes increases with atomic number from La to Lu, where the stability of Y would be between Ho and Er, and the stability of Sc would be after Lu [78, 79].

Experimental studies have shown that at higher temperature (above 150°C) the speciation behavior of the REE deviates considerably from these predictions. At these conditions, the REE chloride and fluoride complexes show a decrease in stability from La to Lu leading to a possible fractionation of the REE [35, 39]. Despite REE forming stronger complexes with fluoride than chloride, a numerical modeling study by Migdisov and Williams-Jones [71] concluded that chloride is the main transporting ligand in natural ore-forming fluids. These simulations were carried out between 200 and 400°C at 1 kbar, focusing on several subsystems, including REE-F-Cl-H-O. Their simulated predictions indicate that, at low to weakly acidic pH, the REE chloride species are dominant, and their predominance field increases to higher pH with increased temperature, leading to dissolved REE concentrations of >1 ppm. At near-neutral and alkaline pH, when the weak hydrofluoric acid (HF0) dissociates to release a proton (H+) and F, fluoride is available to complex with the REE. However, due to the presence of insoluble REE fluorides (REEF3) in their simulations, the REE concentrations were decreased considerably in the fluids. This led Migdisov and Williams-Jones [71] to interpret the fluoride ions as depositional ligands (i.e., leading to precipitation), as opposed to a transporting ligand represented by the chloride ions. Even though solid REE fluorides are not commonly observed in natural systems, their major point was that insoluble F-bearing minerals will control the availability of REE at high pH. These may include fluorocarbonates (bastnäsite and synchysite) and even REE-bearing fluorite.

Our study is in agreement with these previous simulations by Migdisov and Williams-Jones [71], as seen in the speciation diagrams (Figures 13). In our simulations, the REE chloride species were stable at low to near-neutral pH at 150 and 350°C, with increased stability at the higher temperature. Since our models are free of F but contain CO2, the mobility of the REE is not limited by the solubility of REE fluorides at high pH, but by the solubility of the REE hydroxides, which become stable at pH values ranging between ~3.75 and ~6, depending on the simulated conditions. This limits the mobility of the REE at near-neutral pH. However, significant addition of a base such as NaOH at high temperature will buffer the pH to values >7.5–8 to 10 and permit the dissolution of significant amounts of REE due to the increased stability of the species (Figures 13; models at 350°C). This REE hydroxyl species could therefore function as a main transporting ligand, especially in high temperature alkaline fluids that may be present in carbonatite systems [23].

The mobility of the REE will also depend on the solubility of the REE-bearing calcite considering a solid solution between the REE hydroxides and CaCO3 (Figures 4 and 5). In contrast to the speciation model (Model Ia–c), the simulated REE partitioning model (Model II) assumes a pH that is buffered uniquely by fluid-calcite equilibria (without additional NaOH). At these conditions, the predicted transporting ligands for the REE are the bicarbonate/carbonate complexes in a high CO2 and low salinity system, and the chloride complexes in a high CO2 and high salinity system (Figure 5). In addition, the REE are predicted to more soluble <200°C, which is closely related to the retrograde solubility behavior of calcite; this favors the mobilization of the REE. The increased stability of calcite at high temperature will lead to very low REE activities, and REE hydroxyl species seem therefore unimportant for REE mobilization. This implies that, in a calcite vein, lower temperature and the presence of chloride and carbonate/bicarbonate complexes will play an important role for transport of REE in hydrothermal systems.

6.3. REE Partitioning in Natural Calcite

Fluid reservoirs and REE mineral solubilities are commonly assumed to be the main limiting factors for the partitioning of REE between calcite-fluid in natural systems [27]. Chondrite-normalized REE profiles of natural calcite are generally LREE-enriched, due to their closer ionic radii to Ca (i.e., 1.03 Å for La3+ and 1.00 Å for Ca2+ in 6-fold coordination). In carbonatites such as at Bear Lodge, WY, calcite compositions display variations in their REE fractionation behavior between LREE and HREE, which has been ascribed to hydrothermal processes [9, 31], and this is also observed in other carbonatite complexes [32]. A compilation of REE concentrations in natural calcite from carbonatite complexes (Table 1) reveals that magmatic REE profiles are commonly LREE-enriched, whereas hydrothermal overprint may lead to enrichment in middle (M)REE and HREE, leading to characteristic REE profiles as shown in Figure 6.

It is presently difficult to quantify the exact processes responsible for these observed compositional variations of calcite in natural systems. Even though our simulations show the closed-system behavior of the REE solubility during calcite-fluid interaction while natural systems are expected to display dynamic open-system behavior, we can still begin to quantify some of the REE variations observed in natural calcite. The fluid composition plays an important role for the speciation and availability of the REE, in particular pH, salinity, and XCO2 (Figures 13). In addition, the solubility of REE minerals plays a significant role for the mobility of the REE, though simple REE compounds such as the REE hydroxides are not commonly formed in natural systems. One of the reasons for this is the very low solubility of REE minerals such as the fluorocarbonates bastnäsite-(Ce) and parisite-(Ce) [75]. However, some calcite veins are devoid of these minerals, and the solubility of the REE hydroxides is controlled by the retrograde solubility of calcite (Figures 4 and 5), which could explain the absence of these minerals in natural systems at such conditions.

Comparison of the simulated REE concentration using our partitioning model (Model II) with data of natural calcite compositions from Bear Lodge and other carbonatite complexes (Figure 6) indicates that lower temperatures and/or higher salinities will favor fractionation between LREE and HREE. In contrast, increasing temperature will tend to yield flatter REE chondrite-normalized profiles. Despite the simplified approach taken in the present study for simulating the partitioning of the REE, the model fits remarkably well with the REE profiles of natural hydrothermal calcite from Bear Lodge. Besides temperature and salinity, the fluid/rock ratio affects the overall enrichment of REE in calcite for a given initial REE concentration of the hydrothermal fluid (i.e., the model had a starting concentration of 100 ppm for each of the REE). In these simulations, the initial calcite to fluid ratio was set to 1000 g calcite per kg fluid to reproduce absolute REE values within the range observed in natural calcite. Using lower calcite/fluid ratios results in a shift towards higher chondrite-normalized profiles, but it does not affect the overall REE patterns observed in Figure 6. These findings could potentially be used in future studies to determine the calcite/fluid ratios in well characterized natural systems. Additionally, some variations can be seen between neighboring REE of the simulated chondrite-normalized calcite compositions, whereas the natural calcite compositions generally display smooth chondrite-normalized trends. The exact nature of these variations may be related to either the initial REE concentrations in the model or to the solid solution model assumed in this study. Determining these differences between the partitioning model and natural systems will require additional well-constrained experimental investigations.

7. Conclusions

To better understand the mechanisms of REE mobilization and mineralization in natural hydrothermal ore-forming systems requires a fundamental knowledge of the effects of fluid chemistry on metal transport and mineralization processes. To our knowledge, this study is the first attempt to model the partitioning of REE between calcite and fluids in hydrothermal ore deposits. The simulated REE partitioning model indicates that, for a system controlled by calcite-fluid equilibria (e.g., in a calcite vein), the main transporting ligands for the REE will be the bicarbonate/carbonate complexes in a high CO2 and low salinity system, and this shifts to the chloride complexes in a high CO2 and high salinity system. In addition, the REE are predicted to be more soluble at <200°C, which is related to the retrograde solubility behavior of calcite. In more alkaline fluids, REE hydroxyl complex may start to play an important role for the solubility and transport of the REE in high temperature fluids.

The current lack of experimental data for the REE hydroxyl and carbonate/bicarbonate complexes at elevated temperature leads to uncertainties in our prediction capabilities of the mobility of REE in hydrothermal fluids. Additional experimental work is needed to provide new REE partitioning data for calcite-fluid at hydrothermal conditions, which is part of our ongoing investigations. Further development of the simulated model may be extended to other types of ore deposits where hydrothermal calcite plays an important role to vector fluid-rock reactions and ore mineralization.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.


This study was supported by the Colorado School of Mines with a professional development grant to Dr. A. P. Gysi. Additionally, the authors are very grateful for the insightful comments provided by Dr. A. A. Migdisov.