Multicomponent Model for the Prediction of Nuclear Waste/Rare-Earth Extraction Processes.

We develop a minimal model for the prediction of solvent extraction. We consider a rare earth extraction system for which the solvent phase is similar to water-poor microemulsions. All physical molecular quantities used in the calculation can be measured separately. The model takes into account competition complexation, mixing entropy of complexed species, differences of salt concentrations between the two phases, and the surfactant nature of extractant molecules. We consider the practical case where rare earths are extracted from iron nitrates in the presence of acids with a common neutral complexing extractant. The solvent wetting of the reverse aggregates is taken into account via the spontaneous packing. All the water-in-oil reverse aggregates are supposed to be spherical on average. The minimal model captures several features observed in practice: reverse aggregates with different water and extractant content coexist dynamically with monomeric extractant molecules at and above a critical aggregate concentration (CAC). The CAC decreases upon the addition of electrolytes in the aqueous phase. The free energy of transfer of an ion to the organic phase is lower than the driving complexation. The commonly observed log-log relation used to determine the apparent stoichiometry of complexation is valid as a guideline but should be used with care. The results point to the fact that stoichiometry, as well as the probabilities of a particular aggregate, is dependent on the composition of the entire system, namely the extractant and the target solutes' concentrations. Moreover, the experimentally observed dependence of the extraction efficiency on branching of the extractant chains in a given solvent can be quantified. The evolution of the distribution coefficient of particular rare earth, acid, or other different metallic cations can be studied as a function of initial extractant concentration through the whole region that is typically used by chemical engineers. For every chemical species involved in the calculation, the model is able to predict the exact equilibrium concentration in both the aqueous and the solvent phases at a given thermodynamic temperature.

Introduction

In the context of selective recovery of rare earth elements (REEs) or removal of lanthanides and actinides in nuclear waste processing, liquid–liquid extraction basically represents the first and only choice for the development of efficient large-scale processes.[1,2] These two overlapping fields, REEs recycling and nuclear waste processing, in fact constitute a major branch of hydrometallurgy. The importance of this branch is tightly linked to and influenced by industrial and the economic growth worldwide.[2−5]

REEs have a major role in sustaining a green, low-carbon economy. Their numerous applications include permanent magnets, lamp phosphors, batteries for hybrid cars, etc. However, use of REEs also puts importance on their recovery from production scraps and end-of-life products.[1,4−6] Extraction of lanthanides and actinides constitutes a challenge in the field of nuclear energy, as they represent the second stage of purification of spent fuel in fission reactors. In the French Alternative Energies and Atomic Energy Commission (CEA), the implemented process involves the separation of plutonium and uranium in the first stage (the PUREX process) and co-extraction of the remaining lanthanides and actinides in the second stage (the DIAMEX process).[7,8]

An efficient liquid–liquid extraction process requires a particular extractant molecule dissolved in a solvent, which in turn forms a sort of weak complex with the target metal cation.[9] In the case of DIAMEX process, the extractant is DMDOHEMA (N,N′-dimethyl-N,N′-dioctyl-2-(2-hexyloxyethyl)malonamide).[7] The extractant molecules tend to self-assemble into reverse aggregates of various compositions. The compositions depend on the type of the extractant, the extracted solute, the temperature, pH, the salt concentration in the aqueous phase, and the type of the organic solvent.[10−13]

Over the years, various liquid–liquid extraction processes for different systems have been optimized and augmented on the industrial scale. Along with the well-understood experimental methods, a few thermodynamic models have been proposed.[14,15] Still, most of the models are based on the principle that all possible equilibria are established and then fitting over the experimental data is employed.[16] As a result, a set of apparent equilibrium constants and adjusted parameters is obtained.[1,17,18] Often there are more adjusted parameters than observable quantities. Even though satisfactory fits are usually obtained, the question arises whether such models could, in fact, be generalized for different extractants. Moreover, water molecules and organic solvent are neglected within the law of mass action.[19] By neglecting a solvent’s influence, the obtained apparent ion transfer constants do not depend on the branching of the extractant’s alkyl chains, which is completely opposite to experimental observations. Note also that the majority of the models are made to reproduce results of simple laboratory-size systems (which take into account only a few components).

In order to acquire insight into the forces that influence the aggregation process and the overall extraction, we propose a model derived from statistical thermodynamics coupled with ideas and models of molecular self-assembly of the extractant molecules. Our goal is to propose a new methodology with a minimal number of adjusted parameters. Moreover, the parameters themselves ought to reflect the molecular nature of reverse aggregate constituents. Within this paper, we aim to show how our model captures some of the relevant properties of systems formulated from an efficient extraction. An example will be made using the practical case of multicomponent extraction for acidic media. The system in study is composed of heptane solvent containing DMDOHEMA extractant and an aqueous phase containing three extracted solutes, namely, nitric acid, europium(III) nitrate, and iron(III) nitrate. We will also show how the choice of geometrical parameters, widely used to describe curved interfaces, influences the microscopic picture of reverse micelles. Finally, a quantitative description of the reversible extraction formulations will be provided.

Theory

Free Energy of the Reverse Micelle

The model system consists of two phases in contact, namely, the aqueous solution containing multiple ion species and the organic phase containing monomeric extractant molecules (DMDOHEMA) and self-assembled water-in-oil reverse micelles (i.e., the aggregates). We keep this historical name even though we always have a minimum of three components, namely, the water, the extractant, and the solvent. All species are in thermodynamic equilibrium. The aqueous solution represents the brine, whereas the organic phase is the solvent phase in the DIAMEX process. A model system is applicable to any hydrometallurgical process which utilizes an uncharged extractant. It must be noted that hydrophilicity in terms of partitioning of the DMDOHEMA extractant in aqueous phase is not taken into account. It was show that, globally, the effect is small for uncharged chelating and solvating extractants.[20,21] The aggregates are considered as spheres with two distinct parts. The outer part is assumed to be a layer of extractant molecule chains with the average length lchain. The inner part, or the core of the aggregate, consists of extractant polar head groups immersed in the droplet of an aqueous solution of ions that takes up the remaining volume, Vcore.

The free energy of an aggregate of particular composition, , can be written as follows:where Fchain is the free energy associated with the layer of extractant molecules (in short, the chain term) and Fcore is the free energy of the core of the aggregate, defined asFdroplet is the free energy of a droplet of aqueous electrolyte solution, and Fcomplex is the term which describes interactions between cations and extractant head groups. Fdroplet and Fcomplex will be discussed later in this section. Fchain has been taken into account in an already well-established approach through the relation[14]where Nl is the number of extractant molecules assembled into a reverse aggregate, κ* represents the generalized bending constant for one extractant molecule in the extractant film, p is the packing parameter of the extractant molecule, and p0 is the intrinsic spontaneous packing parameter for a certain type of extractant.[22−24] For a fixed chain length, the packing parameter can be written in an explicit form as[25−27]where the radius of the core isand Ni, Nw, and Nl are respectively the numbers of ions, water molecules, and extractant head groups that constitute the core of the aggregate, whereas Vm,i, Vm,w, and Vm,l are respectively the specific molar volumes of ions, water molecules, and extractant head groups. The molar volumes are taken from the literature.[28]

To calculate Fcore, first we need to calculate the free energy of a single droplet of aqueous solution, Fdroplet, immersed in a medium characterized with a low dielectric constant. The partition function Z̃ in a canonical ensemble for a single droplet of an aqueous solution of ions can be written as[29,30]where the index k sums over the total number of particles (ions and water molecules) in the droplet, h is Planck’s constant, r and p are respectively the position and momentum of each particle in the droplet, β is defined as β = 1/kBT, where kB is the Boltzmann constant and T is the thermodynamic temperature, m is the mass of the k-th particle, and V is the interacting potential among particles. The factorials in the denominator of eq account for the indistinguishability of particles. For the sake of simplicity, we consider all particles as spherical objects with no internal degrees of freedom (as sketched in Figure ). In such a formulation, the integral of eq over momenta of the particles giveswhere Λw and Λi are the effective de Broglie thermal wavelengths of water molecules and ions, respectively. Equation holds for large numbers of particles. In contrast to such conditions, our droplet is composed of typically up to 10 particles, with small variation depending on the composition of the polar core. Therefore, Z̃ needs to be corrected for such a small number of particles. If Stirling’s approximation is written as N! ≃ Ne–, thenwhere Z̃approx is the canonical partition function for a droplet of aqueous electrolyte solution, corrected for the small number of particles. Free energy in a canonical ensemble can be evaluated through the relation Fdroplet = −kBT ln Z̃approx, and the following expression is obtained:where Fdroplet denotes the free energy of the droplet made out of a small number of particles. After applying logarithm rules and sorting all the terms, we end up with the following expression:The first terms of the right side of eq represent a correction of the partition function for a small number of particles, whereas the last term, kBT ln Z̃, represents the free energy Felect of the equivalent system in the bulk, where Stirling’s approximation is applicable. Since we consider the condensed phase, namely, a liquid, we can equalize Felect ≃ Gelect, where the Gibbs energy of the electrolyte solution can be written aswhere μworg and μiorg are respectively chemical potentials of water molecules and ions confined in the core of the aggregate.[31] Within this paper we consider only ideal aqueous solutions both in the core of the aggregate and in the aqueous phase with both activity and osmotic coefficients equal to 1.[31] Furthermore, an approximation has been made that the standard chemical potentials of ions and water confined inside this droplet are the same as the ones in aqueous solution in contact, i.e., the same reference state is understood.[32] Therefore, we obtainandwhere μiorg, μworg, μi°, μw°, miorg, mi°, xiorg, and xworg are respectively the chemical potentials of ions and water in the core of the aggregate, the standard chemical potentials of ions and water in the core of the aggregate, the molal concentration of ions in the core, the molal concentration of ions at standard state, and the mole fraction of ions and water in the core. Equation is the consequence of eq when the Gibbs–Duhem relation is used.

Figure 1

Schematic representation of the extraction process. Various types of aggregates are present in the solvent, and their probability at equilibrium is determined by the composition of their cores. Considering the surfactant nature of the extractant, the interface is at least partially covered by the extractant molecules (not shown here).

In order to calculate the free energy of the core, Fcore, we still need to include the complexation free energy term, Fcomplex. The complexation energy per particle, E0,Cat, is defined as a favorable interaction that lowers the potential energy of the system, thus benefiting the extraction. It takes into account the first-sphere interactions between extractant head groups and solvated ions.[14] It typically represents the bond energy measured in the EXAFS measurements.[33,34] Now we can rewrite the corrected partition function Z̃approx in the following way:where NCat depicts the number of ions that are interacting with the extractant head groups (with appropriate energy E0,Cat). Since eβ is a constant term, it can be extracted from the integral. We obtain

As there are multiple head groups available to interact with the “complexing” particle, the actual number of microstates scales with the multiplication factor Ncomplex within the partition function. Ncomplex takes into account all possible configurations of interacting extractant head groups and complexing particles. If we assume that every particle interacts with two head groups, then Ncomplex can be written as[35,36]Note that the denominator in eq does not contain a factorial for indistinguishability of the complexed particle, since it has already been included in the expression for Z̃approx. The term 2 in the denominator accounts for the swap of two head groups to prevent double counting of the same configurations. Now the free energy of the core of the aggregate can be written aswhich givesTherefore, the complexation free energy term (recall eq ) is

We can conclude that Fcomplex consists of configuration entropy term −kBT ln Ncomplex and the internal energy term described as −NCatE0,Cat. Indeed, this was a desired outcome when we defined Fcomplex in eq and implemented the additional stabilizing potential E0,Cat in eq .

When all terms from the expressions for Fcore and Fchain are summed, the free energy of the aggregate of particular composition is obtained (eq ). The whole expression can be found in Appendix A. In fact, is by definition the standard chemical potential, , of the reverse micelle in the particular organic solvent (recall that the partition function was written for a single droplet of an aqueous solution with complexation).[37] Therefore, we can write

Special Case of Pure Water Extraction

The derivation so far was concerned with aggregates containing an electrolyte in the core. When only water molecules are present in the core, i.e., when we consider a pure phase, the corrected partition function, Z̃pure,w, of such a system reduces towhere Vpure,w is the potential between water molecules. The free energy of the pure water droplet is thenwhich is equal to

In the absence of all terms describing the complexation, Fcore is equal to Fdroplet, whereas Fchain is calculated in the same manner as for the case with the ions present. It is obvious that addition of an electrolyte imposes much greater complexity in the system.

Global Equilibrium

Within this model we consider equilibria between an aggregate in the organic phase (solvent) and its constituents, namely, extractant molecules, and solvent, ions, and water molecules in the aqueous phase. We can write the law of mass action aswhere W, Ii, L̅, and are respectively symbols for the water, ions in an aqueous phase, and extractant and aggregate in a solvent. Note that since DMDOHEMA is a neutral extractant, the cation is always transferred from aqueous phase to solvent along with an appropriate number of nitrate anions (NO3–) to balance the charge. Salt molecules considered within this work are HNO3, Eu(NO3)3, and Fe(NO3)3. The chemical potentials of species involved in the chemical reaction described by eq can be written aswhere , μl, μl°, , cl, and c° are respectively the chemical potentials, the standard chemical potentials, the equilibrium molar concentration, and the concentrations at standard state of the aggregates and the extractant in a solvent. μiaq, miaq, μwaq, xiaq, and xwaq are respectively the chemical potentials and the molal concentration of ions, the chemical potential of water, and the mole fractions of ions and water in the aqueous phase. Equation is the consequence of eq when the the Gibbs–Duhem relation is used. To complete the calculation we need to write the law of mass action (eq ) in terms of chemical potentials of all involved species:which is equal toAt this point, it is convenient to definewhere is the reduced standard chemical potential of the aggregate, obtained by subtracting the chemical potentials of ions and water confined in the core from . still contains all the other terms, namely, the chain, the complexation, and terms for the correction for a small number of particles. After inserting eq into eq , the standard chemical potentials of ions and water cancel out, and the following expression is obtained:After some rearrangement of eq , multiplying with β, and exponentiation of the whole expression, we obtainwhere D is defined as

To satisfy the law of mass action, the polynomial of degree l (eq ) needs to be solved for any composition of the aggregate. The equilibrium concentration of extractant in the system, cl, is the root of the polynomial of degree l. Furthermore, the system needs to be solved in such a way that the following conditions are satisfied:where cltot and nitot are respectively the total or the “initial” molar concentration of extractant and mole number of each ion species in the system. Both quantities are inputs of the model. The sum goes over all possible configurations of the aggregate.

To conclude this section it is worth noting that, instead of concentrations at equilibrium, aggregates will be described through their probabilities, since normalized quantities are easier to discuss. The equilibrium aggregate probability, P(w,l,i), is defined aswhere is the equilibrium concentration of the aggregate with a particular composition. The sum in the denominator goes over all possible aggregate concentrations.

Results and Discussion

Input for the Model

In order to perform the calculations, we need a certain set of measurable quantities as inputs for the model. We require the molar volumes of water, acid, salts, and extractant head groups (Vm,w, Vm,i, Vm,l).[28] The length of the extractant molecule chains averaged of all conformations (i.e., the average length), lchain, is also needed to perform the calculations. lchain can be assessed through a combination of small-angle neutron scattering (SANS) and small-angle X-ray scattering (SAXS) pattern fits.[38]lchain of DMDOHEMA molecule chains can also be calculated with molecular dynamics simulations. A recent study showed that results of simulations with explicit n-heptane solvent are in agreement with the experiments.[27]lchain used for our model equals 9.6 Å, and it corresponds to typically 80% of fully stretched chain length. The length of the extractant molecule chains is considered not to change for different compositions of the core of the aggregate.

Within this study, the solvent is not included explicitly in the law of mass action, but it is still indirectly taken into account through the lchain value used for the calculation. The influence of penetrating and non-penetrating solvents on the overall extraction process can therefore easily be included in this model.[39]

Along with measurable quantities, we need to specify the system in terms of initial concentrations of the species, namely, extractant molar concentration, cltot, and ion molal concentrations, mitot.

Beside measurable quantities and initial composition of the system, our model requires a set of parameters, namely, the standard chemical potential of monomeric extractant μl°, the spontaneous packing parameter p0, the rigidity constant κ*, and the complexation parameter E0,Cat, for solutes (except water) that can be extracted into the solvent phase. Among these values, only μl° is accessible by experiments. p0 and κ* can generally be assessed by fitting a three-component phase diagram.[40] The procedure for adjusting the parameters and more detailed discussion about its properties are presented in the next section.

In calculations, the upper limit of both water and extractant molecules was set to Nw = Nl = 10. Our model is self-consistent, so values of Nw and Nl higher than 10 are not required (see next section). Moreover, these intervals represent what is typically observed in the vapor pressure osmometry measurements and in the fits of the SAXS/SANS patterns.[18,41] The complexation term Fcomplex (eq ) is added under the condition that every extracted acid or salt molecule requires at least four extractant molecules to assemble the aggregate; therefore, we impose the rule Nl,min = 4NCat.

It is also important to emphasize that our model is made entirely for the case of spherical micelles, meaning that quantitative interpretation is possible for systems up to cl,initial = 0.605 mol dm–3 of DMDOHEMA (where worm-like micelles usually do not exist) and up to 0.5 mol dm–3 Eu(NO3)3 concentration (before the experimentally observed formation of the third phase occurs).[42]

Model Parameters

This part of the paper is dedicated to adjusting the model parameters and to the study of their influence on the properties of extraction systems. The model parameters (κ*,p0, μl°, and E0,Cat) were adjusted in such a manner that the three crucial conditions were satisfied. The first two conditions are, in fact, the experimentally observed properties of reverse aggregates and extraction systems. The first condition deals with the composition of the aggregates in terms of aggregation number and water content, whereas the second condition ensures that the calculated critical aggregate concentration (CAC) is in accordance with the experimentally observed one. The third condition has more to do with the numerical nature of the calculation. The method needs to be self-consistent, and all results ought to be invariant to the choice of an upper limit of Nl, Nw, and Ni used in calculations.

By satisfying these three crucial conditions, we end up with a small domain of possible sets of parameters, which means that the predicting power of the model is greatly enhanced. The prediction of an extraction process for various species is then a consequence of satisfying those conditions.

Note that the parametrization has been done step-by-step, which means that κ*, p0, and μl° are adjusted for the extraction of water and are therefore considered fixed in later fitting of E0,Cat for each solute. This is the only way to preserve reproducibility.

Before discussing the importance of a proper value of p0 used in the calculations, we need to study the dependence of the packing parameter p on aggregate compositions within the framework the concept used (eq ). We have already stated that the average chain length, lchain, is considered as a constant whatever the composition of the polar core of the aggregate. With such an approximation, the packing parameter p can be calculated using eq . Figure a shows the map of calculated p as a function of the composition of the core of the aggregate when one salt molecule, namely Eu(NO3)3, is present inside. Note again that these values are a consequence of lchain = 9.6 Å (heptane solvent). The different solvent would imply use of a different lchain value, thus changing the map of calculated p. An important feature that can be seen in Figure a is that p asymptotically approaches a value p ≈ 2.6 for the high aggregation numbers, Nl, and the water molecule content, Nw. Note that if we were to calculate p for a very large number of extractant and water molecules (when Nl, Nw → +∞), p would approach a value of 1 (eq ), which means that, by huge swelling of the aggregate, we would end up in the lamellar phase (plane-like structure).[43]

Figure 2

(a) Calculated packing parameter p as a function of the composition of the core of the aggregate (i.e., map of packing parameter for the reverse micelle). (b) Squared difference between calculated parameter p and spontaneous packing parameter p0 = 3.5 as a function of composition of the core of the aggregate. In both figures the core contains one salt molecule, namely Eu(NO3)3, Nl depicts the number of extractant molecules (the aggregation number), and Nw depicts the number of water molecules present in the core (the water content).

Since Fchain is a function of p (eq ), obviously it is important to use a proper value of spontaneous packing parameter p0 for the calculation. p0 describes the position of the “chain energy valley”, where Fchain is close to 0 or sufficiently small to allow the formation of the aggregates. Consequently, the calculated equilibrium aggregate probabilities will be affected by the choice of p0.

Equation represents a simple harmonic approximation of the potential of mean force.[27] Therefore, it is convenient to study the squared differences between calculated p (each p corresponds to a particular composition of the core of an aggregate) and p0, as plotted in Figure b. Figure b shows (p – p0)2 calculated for p0 = 3.5. The calculations have also been made with p0 = 2.5 and 3. The results of these calculations are presented in Figures S1 and S2 in the Supporting Information. It is worth mentioning that, even though the potential well is not perfectly symmetrical with respect to p0, the approach is still quite suitable for the description of Fchain and provides acceptable results.[44−46]

When p0 = 3.5, the preferred compositions correspond to four extractant molecules, and the number of water molecules varies between 1 and 7, depending on the type of solute present inside the core (Figure ). This outcome of the model is in agreement with experimental reports and theoretical studies.[11,18,27,33,37,41,47,48] For p0 = 2.5, the valley of low Fchain will correspond to high numbers of extractant and water molecules. A similar but less pronounced effect is obtained when p0 = 3. Favored aggregation numbers are then 5 and 6, but the water content is still very high. The compositions of such aggregates do not correspond to experimental observations. The equilibrium aggregate probabilities calculated for p0 = 2.5 and 3 are presented in Figures S3 and S4 in the Supporting Information.

Figure 3

Calculated equilibrium aggregate probabilities as a function of the composition of the core of the aggregate. Nl depicts number of extractant, whereas Nw depicts number of water molecules present in the core. Scaled results of the upper left figure are shown in its inset. The model parameters are p0 = 3.5, κ* = 16 kBT per extractant molecule, μl° = 2.5 kJ/mol, = 5 kBT, , and per complexed ion. The system in study is as follows: cl,initial = 0.605 mol dm–3, = 3 mol kg–1, and = = 0.05 mol kg–1.

Furthermore, there is another limitation for use of p0. For p0 = 2.5 and 3, the method is not self-consistent. This means that the range of Nl and Nw affects the prediction of overall extraction (for all solutes present in the system). The result of the calculation ought to be invariant to the upper limit of Nl and Nw. That is another constraint which may as well be a crucial condition when deciding what p0 value to take for the calculation. The case of p0 = 3.5 gives a self-consistent calculation, where large Nl and Nw do not contribute to the result of the calculation and can, therefore, be neglected.

Another geometrical parameter in our model is the generalized bending constant, κ*. Therefore, we studied the influence of κ* on the “valley” of low Fchain and the equilibrium aggregate probabilities. Previously reported values of adjusted κ* for reverse aggregates were 2.5 kBT per DMDOHEMA molecule.[25,26] Recently, a molecular dynamics study done by our group provided a value of 16 kBT per extractant molecule in heptane solvent.[27] This value points to very high curvature toward water (in the reverse micelles) and was attributed to the strong interactions of Eu3+ and DMDOHEMA molecules. As a starting point, we used this value for calculations and then changed it by ±10 kBT per extractant molecule. The results are provided in Figures S5–S7 in the Supporting Information. κ* influences the width of the chain energy valley in such a way that an increase of κ* increases the gradient of Fchain plane. Consequently, higher values of κ* (in our calculation 26 kBT per extractant molecule, Figure S7) allow assembly of aggregates with a smaller number of water molecules. For κ* = 6 kBT, the valley of low Fchain is rather wide, so dilution of the core of the aggregate is highly favorable. As a consequence, aggregates of very different stoichiometry coexist. Most such aggregates are unrealistic compared to the experiments.

It can be concluded that higher rigidity of the extractant film reduces polydispersity in terms of water content. Recall that monodispersity in terms of aggregation number is already achieved through the choice of p0 and imposed rule Nl,min = 4NCat in the definition of the complexation term.

The choice of κ* also affects the overall extraction. By varying κ* from 26 to 6 kBT, the calculated distribution coefficient varied from 14.6 to 5.6, i.e., nearly a factor of 3. This counterintuitive variation of extraction efficiency with the branching of the chains is always observed in industrial applications but has never been predicted by any model of extraction to the best of our knowledge.[39] The variation of κ* was also reflected in the CAC value, which was decreased with the decrease of rigidity.

The choice of κ* = 16 kBT per extractant molecule is acceptable in both realistic aggregate compositions and the overall extraction prediction.

So far, we have introduced and validated the parameters p0 and κ*. In order to calculate the extraction of solutes, we still require standard chemical potential of extractant molecules, μl°, and complexation energy, E0,Cat, for the particular solute.

The value of μl° determines the transition energy between monomeric and aggregated states of the extractant. In our study, the aggregated state has a form of reverse micelles. It is accessible by experimental methods that can determine the mole fraction of the unbound extractant, e.g., NMR shift techniques, scattering extrapolated to zero micelle concentration, or derivatives analysis of liquid–liquid surface tension.[49] In the case of common extractants, the three techniques provided the same result.[50]

Increasing μl° lowers the transition energy between the two states (eqs –34), thus favoring the micellization and the extraction of solutes. A favored micellization is seen as a decrease in CAC and increase in distribution coefficients. A notable property is that the calculated equilibrium aggregate probabilities are invariant to the change of μl°. μl° was fitted accordingly to the experiments, and the value obtained in our study was 2.5 kJ/mol.[18]

In order to obtain E0,Cat for each solute, namely HNO3, Eu(NO3)3, and Fe(NO3)3, we made a fitting based on different known studies. First, we fitted the model to the data concerning the extraction of HNO3 alone.[11,42,51,52] This yielded = 5 kBT, which is a typical order of magnitude for the hydrogen bond formation. In order to obtain , we fitted the experimental data of the HNO3/Eu(NO3)3 system using .[13] Fitting resulted in = 15.6 kBT. The same procedure was followed for the HNO3/Fe(NO3)3 system, which in turn yielded = 13 kBT.[53] This type of fitting was proposed in order to “isolate” the complexation parameter of each particular cation (i.e., particular solute molecule). By doing this, of course, we neglected any type of interaction between different solutes in the organic phase, and also we forbid the existence of mixed-solute aggregates. Obviously, we made a very crude approximation, but this still makes a good starting point for a study of complex multicomponent systems. E0,Cat in our study can be associated with the extraction free energy, ΔG0, from previous studies. ΔG0 was defined as the difference between the free energy of an ion in an aqueous phase and that of an ion complexed by extractant molecules in solvent phase.[14] This definition allows a measurement of ΔG0 by a combination of calorimetry and EXAFS measurements (coupled with ab initio calculations).

In the remainder of this paper, we will show how this model captures a majority of the specific properties of extraction systems.

Equilibrium Aggregate Probabilities

With all the parameters determined and discussed, we have performed calculations in an attempt to quantitatively describe a properties of the extraction systems. All the following results are obtained from the calculation with p0 = 3.5, κ* = 16 kBT per extractant molecule, μl° = 2.5 kJ/mol, = 5 kBT, = 13 kBT, and = 15.6 kBT.

The first thing to study is the aggregate compositions. Figure shows calculated equilibrium probabilities of the aggregates for the practical system (the type of the extractant is used in DIAMEX process and aqueous phase represents a system for a REEs recovery via hydrometallurgy) composed of = 3 mol kg–1, and = = 0.05 mol kg–1 aqueous solution in contact with cl,initial = 0.605 mol dm–3 extractant in heptane solvent. The results are presented for the case where only one salt molecule is inside the core of the aggregates since the calculations showed that addition of second salt molecule pays the huge penalty in terms of chain energy thus making the probabilities negligible. The upper left graph (Figure ) shows the probabilities of aggregates filled only with water molecules. Compared to the other types of aggregates, their probability is more than 10 times smaller. This means that in the case of concentrated aqueous solutions in contact with solvent phase, the water extracted to organic phase originates almost entirely from aggregates containing acid or salt molecules. This result is quite understandable since it is known that upon the addition of acid or metal cations to the aqueous phase, the amount of co-extracted water increases.[13,42,52] The upper left and the two bottom graphs (also Figure ) show the probabilities for aggregates containing, in the same order, Eu(NO3)3, HNO3, and Fe(NO3)3 molecules. In contrast to the case when only water molecules are present inside the core, the highest probabilities in terms of Nw are at compositions from 4 to 5 water molecules per aggregate for HNO3, whereas from 6 to 8 water molecules for Eu(NO3)3 and Fe(NO3)3 molecules. This is a consequence of the differences in chemical potentials of ions and water between the aqueous and the organic (solvent) phases. The higher water content for Eu(NO3)3 and Fe(NO3)3 salts is due to the higher number of particles in the core of the aggregate. The dilution of the core is a highly favorable effect that stabilizes the aggregates core, but it works in opposite way with respect to the Fchain because the inclusion of additional particles causes the increase in core radius (an unfavorable swelling of the reverse micelle). Another important thing to add here is that the probability of particular aggregate is not only governed by the dominant term E0,Cat but also depends on the composition of the entire system in study, i.e., on the initial extractant, acid and ion concentrations, temperature, etc. For example, if there is a change of reservoir acid or metal cation concentrations, the distribution of aggregates, i.e., the probabilities at equilibrium, will be different. To emphasize this, we performed the calculations for different concentrations of acid and lanthanides. The results are presented in Figures S8 and S9 in the Supporting Information. This important property is also reflected in the aggregation threshold, the apparent stoichiometry, and the overall extraction efficiency, as we will show in the following part of the discussion. We emphasize again that all subsequent calculations were done with values of parameters p0 = 3.5, κ* = 16 kBT per extractant molecule, μl° = 2.5 kJ/mol, = 5 kBT, = 13 kBT, and = 15.6 kBT per complexed ion.

Predicting the Aggregation Threshold, i.e., the Critical Aggregate Concentration

We have already stated that one of the crucial properties of extraction system is measured CAC. We tried to mimic different systems, which are often found in the literature and to calculate corresponding CAC. Figure shows the calculated equilibrium concentrations of both monomeric and aggregated extractant molecules as a function of initial (or total) extractant concentration. Figure a shows the concentrations of the extractant for LiNO3 aqueous solution in contact with a solvent, which represents a typical background non-extracted salt system. As shown in the figure, CAC for such a system is evaluated to be 0.13 mol dm–3 initial extractant concentration, and only water molecules are present in the aggregates. Figures b–c show respectively the system of a HNO3 aqueous solution and a mixture of HNO3, Eu(NO3)3, and Fe(NO3)3 in contact with the solvent phase. The details of the concentrations in each system are given in the caption of Figure . Upon the addition of acid and cations in the system, CAC decreased to 0.07 and 0.06 M, respectively.

Figure 4

Calculated concentrations of monomeric and aggregated extractant as a function of the initial (or total) extractant concentration. Solvent phase in contact to (a) mLiNO = 3 mol kg–1 (nonextracted salt), (b) = 3 mol kg–1, and (c) = 3 mol kg–1 and = = 0.05 mol kg–1. The black dashed line shows the calculated critical aggregate concentration.

Figure demonstrates that the model captures a known decrease of CAC upon the addition of acid in the system. The results of the CAC as a function of initial HNO3 concentration in aqueous phase presented here are made for the system cl,initial = 0.6 mol dm–3, = = 0.05 mol kg–1. The enhancement of the micellization is due to the two factors contributing to the free energy of the aggregates. The first and more dominant factor is the increase of the solute concentration able to make a weak complex with the extractant molecule (defined by complexation energy ). The second factor is an increase of NO3– concentrations ratio between the aqueous phase and the core of the aggregate (eq ). The properties of the extracting systems shown in Figures and 5 correspond to experimental findings.[48,51]

Figure 5

Critical aggregate concentration as a function of HNO3 concentration. The system in study is cl,initial = 0.6 mol dm–3, = = 0.05 mol kg–1.

Prediction of the Extraction Process

It is often convenient to show the evolution of solute concentrations in the organic phase as a function of cl,initial. Therefore, in Figure we show the extraction curves of all species for a system = 3 mol kg–1, = = 0.05 mol kg–1. An important feature is the low extraction of HNO3 compared to the literature. This is a consequence of our approximation that only one type, acid or metal nitrate molecule, can occupy the aggregate. In future publications, the complexation term will be defined in a more realistic way, which will result in a proper HNO3 extraction curve. Still, different trends can be observed while inspecting Figure . A typical high water uptake and an increased metal nitrate extraction with an increase in cl,initial correlate extremely well with the experiments.[42] Also, it is worth to mention that before CAC the aggregation is controlled by water extraction because of a smaller penalty in chain energy while after CAC, the aggregation is entirely controlled by complexation of metal nitrates. The results showing the concentrations of the extracted solutes as a function of the initial extractant concentration for different acidity and initial salt concentrations are presented in Figure S10 in the Supporting Information. The results show that extraction efficiency is strongly dependent on the composition of the system, i.e., on the initial acid and salt concentrations.

Figure 6

Calculated equilibrium concentrations of all extracted solutes as a function of the initial extractant concentration. The system in study is = 3 mol kg–1, = = 0.05 mol kg–1.

The concentration of acid in the aqueous phase influences the extraction of cations.[42] In order to test our model we made a set of calculations for different HNO3 aqueous concentrations. The extraction efficiency is usually expressed aswhere DCat,i is a distribution coefficient of target solute between aqueous and solvent phases at equilibrium and cCat,iorg is the total concentration of solute in the organic phase.[54] If one wants to express DCat,i as a fraction of molar concentrations, a usual conversion of molality in aqueous phase is made with relation ciaq = ρaqmiaq/(1 + ∑miaqMi), where ρaq is the density of the aqueous solution.

DEu and DFe as a function of cl,initial have been calculated for different (different acidity of aqueous phase) and presented in Figure S8 in the Supporting Information. The results show the nonlinear increase of distribution coefficients with increasing cl,initial. The increase in causes an increase in DCat,i which means that the extraction is enhanced upon addition of HNO3 the system. This is again the consequence of an increase of NO3– concentrations ratio between the aqueous phase and the core of the aggregate.

Another way of plotting these results is by employing so-called “log–log” plot, as it is traditionally made in the slope method to investigate the apparent stoichiometry of the system.[55,56] In this manner we have transformed the data from Figure S8 to decimal logarithms and presented them in Figure . Figure a shows again the extraction for different HNO3 concentrations. It can be noticed that depending on the region of cl,initial the slope of the extraction lines changes and typically three regimes are observed. Moreover, the trend in change of slope is dependent on the acidity of system, i.e., on HNO3 concentrations. Below typically = 2 mol kg–1 the calculations have shown a different behavior than for higher concentrations. In order to see the differences better, we have isolated the graphs for = 0.5 (black line) and 5 mol kg–1 (red line) and plotted them separately in Figure b. What is striking is the fact that slope changes substantially for = 0.5 mol kg–1, depending on cl,initial. A high initial slope of approximately 3.5 is followed by 1.8 and then 1.2. The case of = 5 mol kg–1 shows again the three regimes but in all regions slope is around 1.7. This observation points to the fact that high salt concentrations in aqueous phase tend to damp the fine-tuning influences on the apparent extraction stoichiometry. In practice, by considering a large concentration range the log–log plots are usually not straight lines. There is a deviation at both low and high extractant concentrations. When the central slope is used, the non-integer value is said to correspond to the average effective stoichiometry. If the slope at a given range of concentrations is not an integer, several different complexes are invoked. Furthermore, the complexation at low extractant concentration is higher. Obtaining a larger aggregation number is contrary to the Le Chatelier’s principle. This difficulty has been discussed in chemical engineering for different types of adducts, i.e., molecules which participate in the aggregate but are not complexed. All the extra parametrization which is necessary when the slope method is applied is no longer required using our general model. Without any extra parameters, the extraction and its intrinsic non-linearity are predicted.[57,58] This finding brings attention to the longstanding use of the slope method in determination of the stoichiometry for various hydrometallurgical processes. We wish to emphasize that, under a certain physical condition of the system, various regimes in stoichiometry can be “masked” by experimental error, thus leading to false simplicity in the understanding of the behavior of the system.[59] Besides the concentration of nitric acid, we have varied the concentrations of the lanthanide. As expected, the apparent stoichiometry is dependent on the initial concentration of target solutes. The results are presented in the Supporting Information.

Figure 7

(a) Decimal logarithm of Eu3+ distribution coefficient as a function of decimal logarithm of the initial extractant concentration. The system in study is cl,initial = 0.6 mol dm–3, = = 0.05 mol kg–1. The results are presented for various concentrations of nitric acid in aqueous phase. The equivalent curve for the extraction of Fe3+ is presented in the inset. (b) Slope method results for = 0.5 and 5 mol dm–3.

Figure shows the dependence of DEu from the mixture of salts on the HNO3 concentrations. This represents a part of the multicomponent phase diagram along the fraction of acid in the aqueous phase. Calculations are done for the extractant concentrations after the CAC, which ensures that we are in the regime where the aggregates are the dominant species in the solvent. The results show that, in fact, our model recovers typical Langmuir isotherms that have already been reported by both experiment and modeling.[13,51] For 0.1 mol dm–3, the extractant is entirely saturated and an additional increase of HNO3 concentration cannot enhance the extraction. In contrast, when the concentration of the extractant is 0.6 mol dm–3, there is a sufficient amount of monomers. Adding HNO3 in the aqueous phase pushes the equilibrium toward the creation of the additional aggregates containing Eu3+.

Figure 8

Eu3+ distribution coefficient as a function of nitric acid concentration in aqueous phase. The system in study is cl,initial = 0.6 mol dm–3, = = 0.05 mol kg–1.

Complexation Energy Study and Reversible Formulations

In hydrometallurgy, the extraction of a cation to the solvent phase is often identified as the complexation of a cation by the chelating agent (the extractant molecule). In fact, the affinity to form the complex is indeed the leading force for extraction of the cation, but it is only one of the terms in the global free energy of transfer.[14] In our model, complexation is counterbalanced by a few opposing forces such as energy cost for packing of extractant chains in ordered film, i.e., a steric hindrance, the differences in ion concentrations between the aqueous and the solvent phases, the differences in chemical potentials of water between two phases, etc. (see full expression for free energy of aggregate in Appendix A).

In order to clarify this misunderstanding, we plotted the negative value of the natural logarithm of the DEu as a function of the negative value of the complexation parameter . By declaring as a continuous variable, we mimic the strength (or affinity) of the extractant molecule to form a complex with Eu3+. Different types of extractants are characterized by a different .

The term −kBT ln DEu is often referred to as the apparent energy of the extraction.[26] In Figure , the assumption complexation = extraction (a green dashed line) is plotted for purpose of easier understanding of the context. Figures a−b show the −kBT ln DEu/kBT for a system = 0.1 mol dm–3 and = 1 mol dm–3, respectively. The extractant concentration was fixed to cl,initial = 0.6 mol dm–3. The two different initial Eu3+ concentrations represent the cases below and above the experimentally observed limiting organic concentration (LOC) of solutes. This view is quite useful when developing an efficient formulation for the extraction.

Figure 9

Negative value of natural logarithm of the distribution coefficient −ln DEu as a function of a negative value of the complexation energy parameter . The negative values of are taken for the purpose of visually easier reading of the saturation limit. (a) System: cl,initial = 0.6 mol dm–3, = 0.1 mol dm–3. (b) System: cl,initial = 0.6 mol dm–3, = 1 mol dm–3. (c) System: = 1 mol dm–3. The calculations are shown for the two initial extractant concentrations, namely cl,initial = 0.6 mol dm–3 and cl,initial = 0.3 mol dm–3.

Two immediate conclusions can be drawn. First, the extraction cannot be solely identified with the complexation energy, since solid and dashed lines are not collinear. Second, the choice of the extractant is also dependent on the concentration of the target metal cation and not only on its nature. When a concentration of metal cation is sufficiently low compared to the concentration of the extractant (Figure a), the can be very high (typically that of the ionic charged extractant) and saturation would still not occur. By saturation, in this context, we address the case where most of the extractant molecules are in aggregated form and the concentration of monomeric form is almost negligible. In Figure b for = 1 mol dm–3, the saturation is achieved already for the type of extractant described by typically = 10 kBT per complexed ion, thus showing an irreversible character of the formulation. This means that, in practical formulation, it would be sufficient to use a lower concentration of salt or to exchange DMDOHEMA extractant with some less efficient one.

An important feature of this model is the use of well-defined and justified parameters. The results presented in Figure provide a sort of justification of E0,Cat in general. E0,Cat has a proper value if and only if, for a given definition (recall Theory section, eq ), it provides a result which is in accordance with the experimental values. A value of = 15.6 kBT is in fact a good value for the description since it corresponds to approximately DEu = 11 (Figure a). A conclusion is that the system composed of = 0.1 mol dm–3 in 3 M nitric acid and DMDOHEMA extractant in a solvent represents a desirable reversible formulation.

Figure c shows the influence of total extractant concentration, cl,initial = 0.3 and 0.6 mol dm–3. The important feature is that saturation is achieved for the same value. The only difference is between the two is that formulation with the higher cl,initial can extract more of the target ion.

In hydrometallurgy, the supramolecular approach stipulates that the extraction free energy corresponds to the complexation of the cation by one or more chelating agents (synonym for the extractant molecules) associated with an entropy of mixing.[9] In the colloidal approach proposed here, a more general view is now possible. That is why we propose Figure to illustrate this.[14] Since the negative value of the natural logarithm of DEu represents the apparent free energy of the electrolyte extraction, we plot this quantity related to the efficiency on y-axis, while x-axis shows the chemical motor driving the transfer toward the phase containing the extractant. The latter is specific to each lanthanide/extractant couple. Figure a shows the case where the mole ratio of the extractant to lanthanide is a factor of 6, and we see that in these conditions, the usual supramolecular approximation holds. The two lines are separated by around 7 kBT, which corresponds to the sum of various contributions included within our model, namely, differences in ion concentrations between the core of the aggregate and the aqueous phase, the differences in chemical potentials of water between two phases, steric hindrance of extractant chains, etc. The opposite case is shown in Figure b, where now the amount of lanthanide cations is in excess relative to the amount of extractant molecules in the system. In this case, the curves exhibit the typical Langmuir isotherm’s behavior since, after the saturation of the extractants, the extraction efficiency no longer depends on the driving complexation energy. We show an order of magnitude for our practical case of DMDOHEMA/Eu3+ to which the complexation parameter has been attributed (via fitting procedure). Last but not the least, in this case, monomer concentrations are negligible, and this favors the danger of going in the three-phase triangles since the oil phase becomes unstable.[12]

Conclusion

In order to acquire insight into the forces that influence the aggregation process and to predict the overall extraction of solutes into a solvent phase, we proposed a minimal model for which the parameters are experimentally accessible. The minimal model was derived from statistical thermodynamics within a framework of molecular self-assembly of the extractant molecules. With this colloidal approach that goes beyond supramolecular chemistry considerations, the efficiency plots can be generated for any point of a Winsor II regime where the dominant aggregates are reverse spherical micelles.

Our model, in a global free energy difference approach, takes into account the dominant term called complexation free energy, which is well known in organometallic chemistry of supramolecular self-assembly. The complexation free energy is counterbalanced by weaker quenching terms associated with the packing of extractant chains, differences in ion concentrations between the two phases, water activity, etc. To the best of our knowledge, the model presented here uses quantitatively, for the first time, both generalized bending constant κ* and spontaneous packing parameter p0 larger than 2 for evaluating phase transfers.

The free energy associated with the film of extractant chains, Fchain, is a function of parameter p0. By adjusting p0, we showed that positioning of the valley of low Fchain greatly influences the composition of the aggregates at thermodynamic equilibrium. Realistic compositions were obtained for p0 = 3.5. Lower p0 favors unrealistically high aggregation numbers and water contents in the core of the aggregates.

Rigidity κ* adjusts the number of water molecules inside the core of the aggregate. A decrease in κ* increases polydispersity and causes an unrealistically high water uptake. A higher κ* decreases the number of water molecules inside the core, thus providing more realistic aggregate compositions, but it also quenches the extraction (and increases CAC). The value we have validated with our calculations is κ* = 16 kBT per extractant molecule.

Fitting of experimental CAC and extraction curves yielded the standard chemical potential μl° = 2.5 kJ/mol, and the complexation energies, namely, E0,HNO = 5 kBT, = 13 kBT, and = 15.6 kBT.

With the obtained parameters, we studied a practical system composed of HNO3, Eu(NO3)3, and Fe(NO3)3 aqueous solution in contact with a solvent containing DMDOHEMA extractant. The calculations showed that the most probable aggregates contain typically one salt molecule, four extractants, and from 4 to 8 water molecules inside the core (depending on the type of salt). Stable aggregates containing Eu(NO3)3 or Fe(NO3)3 are formed with an increased number of water molecules, since more ions in the core require a higher dilution in order to reach a stable form. The probabilities, and thus the concentrations, of the aggregates at equilibrium are dependent not only on the interaction of extractant and extracted solutes (i.e., on E0,Cat) but also on the composition of the entire system, e.g., the initial acid and salt concentrations, temperature, etc.

Our model predicts a decrease of CAC upon the addition of target salts in the aqueous phase. An increase of HNO3 concentration forces higher water co-extraction and also enhances the extraction of metal nitrates, namely Eu(NO3)3 and Fe(NO3)3.

The extractant concentration, especially above CAC for a particular system, plays a significant role in the extraction of target salts. The calculated distribution coefficient versus extractant concentration results shows nonlinear behavior, which is even more pronounced upon increasing the bulk HNO3 concentration. A slope method used to determine the apparent stoichiometry of complexation shows different trends that depend not only on extractant but also on the acid and lanthanide concentrations.

In the context of reversible and therefore desirable formulations for the extraction systems, we have performed calculations with varying . The results show that choice of extractant is dependent on, besides the nature of the target salt, also its total concentration in the aqueous phase. The calculations show that there is a threshold of after which the saturation of extractant is achieved (for a defined salt concentration), thus making an unfavorable formulation. Upon increasing the total extractant concentration, as expected, the extraction capacity of the used formulation increases, but the saturation threshold remains constant. Therefore, saturation in terms of is not a function of the total extractant concentration, cl,initial.

Our calculations also show a clear distinction between extraction of the solute and the complexation energy term.

In our next publication, extension to charged extractants will be considered. This implies ion exchange between oil and water phases, instead of extraction of a neutral salt molecule.[60] Such extension greatly broadens the applicability of the model to many commonly used industrial systems.[2]