Modeling Preparative Chromatographic Separation of Heavy Rare Earth Elements and Optimization of Thulium Purification

Rare Earth Elements are in growing demand globally. This paper presents a case study of applied mathematical modeling and multi objective optimization to optimize the separation of heavy Rare Earth Elements, Terbium-Lutetium, by means of preparative solid phase extraction chromatography, which means that an extraction ligand, HDEHP, is immobilized on a C18 silica phase, and nitric acid is used as an eluent. An ICP-MS was used for online detection of the Rare Earths. A methodology for calibration and optimization is presented, and applied to an industrially relevant mixture. Results show that Thulium is produced at 99% purity, with a productivity of 0.2 0.5 kg Tm per m3 stationary phase and second, with Yields from 74% to 99%.


Introduction
The Rare Earth Elements (REEs) are Scandium and Yttrium and the members of the Lanthanide group in the periodic table.Many of these metals are used in batteries, lasers, capacitors, superconductors [1], which make the purification process demanding through high purity requirements in some of these applications.Current industrial separation methods include liquid-liquid extraction, selective oxidation/reduction, and ion exchange chromatography [2]- [5].Extraction chromatography is commonly used for analytical purposes [1] [6]- [9].
The REEs are of growing economic importance, with China limiting the export quotas [10] [11], which creates a demand for alternative sources, and alternative purification processes, as the currently dominating liquidliquid extraction uses high quantities of organic compounds for the separation.Solid phase extraction or extraction chromatography uses smaller amounts of organic chemicals than liquid-liquid extraction.This is why solid phase extraction could be one of the alternative purification processes.This work is a continuation of the work presented in [12]- [14] where a mixture of the elements from Neodymium to Gadolinium was used in modeling and optimization of a chromatographic step, and also in experimental optimization.In the present work a mixture of all REEs except Scandium was separated on a modified column that operates in a solid phase extraction mode, meaning that the extractant is immobilized as a ligand in the column.To find the optimal operating conditions for the separation, a Dispersive-Equilibrium model was constructed, based on Langmuirian type binding with mobile phase modulator.The parameters for the model were calibrated from isocratic and gradient experiments.The optimal operating points were found by applying an optimization algorithm on the calibrated model and performance functions.In addition to this one of optimal operating points was verified experimentally.

Theory
To design and optimize the chromatographic system, a performance function is required.This in turn depends on a mathematical formulation of the physical chromatography system.This mathematical formulation requires retention parameters, mass transfer parameters and column capacity to be calibrated.To reduce the complexity of the calibration, only the retention parameters are calibrated with a less computationally expensive model, i.e. the retention volume model.These retention parameters are then used as a starting guess in the more advanced model of the chromatography column, which takes dispersion and mass transfer effects into account as well as the retention data.For the optimization, Productivity and Yield were selected as performance functions, subject to a Purity constraint.

The Retention Volume Model
The retention volume model is an ideal model, which can be seen as a simplified version of the film mass transfer model, presented in 2.2, that includes the isocratic modifier concentration effect on the retention volume: where V R is the retention volume, V R0 is the retention volume under non-binding conditions, ε C is the column porosity, ε P is the particle porosity, K D is the exclusion factor and H is the Henry constant [15]: where H 0 is the Henry constant when the modifier concentration is unity, s modifier is the modifier concentration (mol/m 3 ), ν is the stochiometric coefficient.This expression is valid only when the loaded amount is much smaller than the total adsorption capacity.

The Dispersive-Equilibrium Model
The model including mass transfer and non-linear effects is a Dispersive-Equilibrium model using Langmuir type binding with mobile phase modifier, consisting of a dispersive-equilibrium description of the mobile phase with all mass transfer effects lumped in the dispersion, and Langmuir equilibrium adsorption for the solid phase [16].For component i, its liquid-phase concentration in the column is described by the following equation: where c i is the concentration of component i in the mobile phase (mol/m 3 ), D AX is the apparent dispersion coefficient (m 2 /s) see Equations ( 6)- (8).x is the axial coordinate along the column (m), v app is the apparent velocity, which is the superficial velocity divided by the total void fraction (m/s), ε T is the void fraction of the packed bed (m 3 mobile phase/m 3 column), q i is the concentration of adsorbed species i (mol/m 3 column) and t is the time (s).The column equation is subject to two boundary conditions, one at the inlet and one at the outlet.A Robin boundary condition describes the concentration at the inlet: where c inlet,i is the concentration of component i at the inlet.The concentration of species i at x = 0 may be lower than the inlet concentration due to dispersion.The other boundary condition is at the outlet of the column, at x = L, and is a von Neumann condition: The concentration of species i adsorbed on the stationary phase is described by the following equation, which is a mobile phase modified version of the Langmuir adsorption description [21]: where H 0 is the Henry constant, s is the modifier concentration, ν is the stochiometric constant, q max,i is the maximum adsorption capacity for component i.This equation must at equilibrium satisfy the adsorption equilibrium isotherm, which means that the parenthesis must be 0 at equilibrium.The flow rate dependence of the apparent dispersion, including mass transfer resistance and kinetics, is approximated with a second degree polynomial.This comes from the Knox equation [17]: where h is the reduced plate height, v app is the apparent velocity, d p is the particle diameter.A, B and C are constants coming from flow structure and non equilibrium effects.The apparent axial dispersion can then be described as: where D AX is the apparent dispersion coefficient.Combining equations 7 and 8 give: Estimating the apparent dispersion with a second degree polynomial gives a small error in the second term, however this is acceptable with current model precision.

Optimization
The target component for the optimization is Thulium, as producing pure Thulium would produce pure Ytterbium also, and then both of the major components have been purified, see Figure 4 and Figure 5.The rest of the REEs are considered impurities.The optimal performance of the separation system depends on the technical and economic constraints and incentives that are placed upon this system.Common performance attributes used are Productivity, Yield, Purity and Specific Productivity.In this work, Productivity, presented in Equation (10) and Yield, presented in Equation ( 11) are used as performance functions, while the Purity is used as a constraint.The productivity can be calculated from: where t cut,1 and t cut,2 are the times between which the target component can be collected at required purity (s), F is the flow rate (m 3 /s), c outlet is the outlet concentration of the component (mol/m 3 ), t cycle is the total cycle time (s) and V col is the column volume (m 3 ).The yield can be described as: In this case t load is the time at the end of the loading(s), c feed is the concentration in the feed (mol/m 3 ).When two competing objectives, in this case productivity and yield are competing, the result will be a Pareto front, which means that for a given optimal point on the front, increasing one objective must decrease the other [18]- [20].

Materials
The experiments were performed on an Agilent 1260 Bio-Inert HPLC system, consisting of two Agilent 1260 Bio-Inert Quaternary pumps, one Agilent Bio-Inert 1260 Autosampler and one Agilent 1260 Bio-Inert UV/Vis detector connected to an Agilent 7700 ICP-MS for detection of the rare earths.The column preparation was performed using a GE Healthcare ÄKTA Purifier 100.The individual rare earths were acquired from Merck at 1 g/L in nitric acid.The concentrated nitric acid of HPLC grade was acquired from Merck.The columns used were Kromasil C18, 100 Å, 16µm columns with immobilized HDEHP acting as ligand.HDEHP was acquired from Merck.Methanol was acquired from Merck.To be able to do higher flow rate experiments, an extra isocratic pump was added to the system after the column, and used to feed the ICP-MS with a constant flow of 1 ml/min, the rest going to waste or collection.The heavy mixture of rare earths used in the Thulium optimization was mixed to mimic an industrially relevant REE-stream, with a composition as described in Table 1.

Column Preparation
The column was modified by passing a mixture of 22 mM HDEHP in 65% (v/v) Methanol-water, using the conductivity meter of the ÄKTA to determine the breakthrough volume, and thereby the ligand concentration of 345 mol/m 3 column volume.As the ligand forms dimers, the concentration of sites was 172.5 mol/m 3 column volume.No overloading experiments were performed, as the column capacity had been decided during the preparation of the column.

Chromatography Method
The chromatographic experiments performed, consisted of four consecutive steps; loading, isocratic elution, Cleaning in Place (CIP), and reequilibration.
The injected amount during loading was 10 µl of a mixture containing 0.0357 g/L of each of the lanthanides, plus yttrium.This was well within the linear range of the isotherm, so any competitive effects could be neglected.After loading, a prepared eluent of a predetermined concentration was passed through the column for 120 minutes, this was due to a limitation of the ICP-MS software, after elution a CIP at 7M nitric acid was performed.After the CIP the system was reequilibrated with water.All experiments were performed with three replicates.

Validation Experiment Chromatography Method
For the validation experiment, the experimental setup was modified so that it no longer had the post column pump.This decreased the post column dispersion.The validation experiment was adjusted somewhat from the very aggressive points on the Pareto front, to a longer gradient and with a lower flow rate.The validation experiment consisted of five consecutive steps: loading, wash, gradient elution, CIP and reequilibration.The injected amount during loading was 200 µl heavy mixture, diluted to 33.6 g/L with deionized water, then a 4.5 ml wash with deionized water, followed by a 17 ml gradient from 0 to 7 M nitric acid.Then 4 ml CIP at 7 M nitric acid, followed by 10 ml reequilibration with deionized water.

Simulation Method
The chromatography model was discretized in the axial direction, to turn the partial differential equation into a system of Ordinary Differential Equations (ODE) and Algebraic equations, Differential Algebraic Equation (DAE).The discretization used a 2 point backward finite difference scheme for the flow, and a 3 point central finite difference scheme for the dispersion.The assembled DAE system was solved with MATLAB ODE/DAE solver ode 15 s, which solves stiff DAEs.Discretization and Simulation is handled by the Preparative Chromatography Simulator, PCS, developed at the Department of Chemical Engineering, Lund University [21].In order to reduce the size of the computational problem for the optimization, all components lighter than Gadolinium were lumped together and added to the mass of Gadolinium, this is done, as in the worst case perspective of separating Thulium from the heavy mixture, they would act as Gadolinium, where they would pass more readily through the column in reality.

Calibration Method
First the Henry constant and the stochiometric coefficient were calculated from the slope and intercept of the logarithm of the modified retention volume, versus the logarithm of the modifier concentration.These were then used for the inverse method, i.e. minimizing the least squares error for the normalized detector response and the normalized simulated chromatograms, where the mass transfer effects were fitted, using a Nelder-Mead simplex algorithm.The maximum binding capacity of the columns was decided during preparation of the columns, and based on the HDEHP concentration in the backbone.

Optimization Method
The calibrated model was used with the performance functions presented in 2.3, and used with a modified Dif-ferential Evolution (DE) algorithm and parallel computations in a cluster environment [13].For the optimization, four decision variables were chosen, load volume, gradient length, CIP length and flow rate.

Results
The modified columns had some known characteristics, such as the porosities, ε c = 0.6, ε p = 0.4, and the maximum binding capacities, which were decided during the impregnation of the column, and calculated from equation 12: where λ is 172.5 (mol/m 3 column), and ν i is the stochiometric constant for component i.

Model Calibration
The Henry constants and characteristic charges, presented in Table 2, were calculated from the retention volume experiments presented in Figure 1.There is a wide span in retention data for the replicates, as can be seen in Figure 1.This is probably a result of non-isocratic and inexact mixing of the eluent as a result of the rotating inlet valve combined with the non-ideal mixing properties of water and nitric acid.

The Flow Rate Experiments
To calibrate the mass transfer effects, flow rate experiments were performed at 1, 1.25, 1.5, 1.75, 2, 2.5 and 3 ml/min, the normalized detector response was used against a normalized simulated response from the PCS, and parameters fitted by a least squares approach.The experimental and model results are presented in Figure 2, with model parameters for the axial dispersion as in Equation ( 9), with the terms: A = 1.05 × 10 −8 , B = 1.00 × 10 −8 and C = 0.9936 × 10 −8 .

Optimization
The multi-objective optimization was performed using Productivity and Yield, see Equations (10) and (11), as objective functions.The optimization output is a the Pareto front, presented in Figure 3, which gives the system Figure 2. Model response and experimental response for dysprosium at different flow rates.There is a reasonably good fit, and using a least squares approach to model calibration will result in wider peaks where there is variation in the experiments.performance for different operating points, comparing in this case Productivity and Yield, under a 99% Purity constraint for the Thulium.Figure 3 shows that the Pareto front is much denser in the high yield region, than it is in the high productivity region.This is to be expected, as high yields are relatively easy to achieve by keeping the loading volume small, and the gradient slope flat.The selected points in Figure 3 are presented in Figure 4, with their decision variables presented in Table 3.
All of the operating points are under overloaded conditions, as can be seen in Figure 4, with a decreasing degree of overloading as we move from high productivity in chromatogram 1 to high yield in chromatogram 4. The curvature of the Pareto front is a result of peak overlapping, so that in chromatogram 4, the Thulium peak is slightly contaminated on the back side by Ytterbium but base line separated at the front, then going to chromatogram three, where the Thulium peak is contaminated more by Ytterbium at the back side, and starting to see Erbium at the front side.Chromatogram 2 shows even more Ytterbium in the back side of the peak, and a bit of Erbium at the front side as well.In chromatogram 1, the Erbium and Ytterbium peaks almost meet underneath the Thulium peak.

Validation
As described in 3.2.2, the validation experiment was performed at a slightly less aggressive operating point than max (Prod * Yield), with a longer gradient and a flow rate of 1 ml/min.The chromatogram from the validation experiment, presented in Figure 5 shows that base line separation on the front side of the Thulium peak is possible, with a small overlap at the back side, while both Thulium and Ytterbium show overloaded conditions.This   chromatogram has a calculated Productivity of 0.20 kg Thulium per cubic meter of column volume and second, with a Yield above 95%.The reason that the peaks are a lot sharper in the validation chromatogram than the simulated peaks, is that the post column pump was removed, which leads to much less post column dispersion.

Conclusion
Using the inverse method to fit retention data and kinetics at the same time for a system with 15 components is very hard, for this reason, performing isocratic experiments to determine the isotherm parameters, and using the inverse method for the mass transfer effects proved to be an easier path to follow.The model predicts preparative loads under gradient elution well, and therefore is applicable in optimizations of the separation system, for a wide range of mixtures and performance criterion.The performance of the system showed a Productivity range for Thulium from 0.18 kg/(m 3 •s) to 0.45 kg/(m 3 •s) at Yields of 98.5% to 73.3%.

Figure 1 .
Figure 1.Retention data for the experiments, with regressed lines.There is a wide span in retention data for the replicates, as can be seen in the figure above.

Figure 3 .
Figure 3. Pareto front Productivity versus Yield, selected Pareto points correspond to chromatograms in Figure 4.

Figure 4 .
Figure 4. Selected chromatograms from the Pareto front.All chromatograms are overloaded, with decreasing degree from chromatogram 1 to chromatogram 4.

Figure 5 .
Figure 5. Experimental chromatogram.This slightly overloaded chromatogram shows base line separation at the front side of the Thulium peak, and a slight overlap on the back side of the peak.This means high purity and high yield.

Table 1 .
Composition of heavy mixture, 336 g/L, 14 M HNO 3 .Isocratic experiments at 16 compositions of the mobile phase were performed to get retention volumes for all components, gradient experiments at different flow rates were performed to get peak broadening from mass transfer and kinetics.

Table 2 .
The fitted stochiometric coefficients and Henry constants.

Table 3 .
The chromatograms from Figure 3 and Figure 4, with their decision variable values and Productivity and Yield values.