REE Extraction Basics: Using the difflow_ree Module#
This notebook demonstrates the difflow_ree plugin for rare earth element (REE) solvent extraction.
Background#
The difflow_ree module provides:
Database of 15 REE elements with properties
5 extractant systems (D2EHPA, PC88A, Cyanex272, TBP, naphthenic acid)
pH-dependent distribution models
Unit operations for extraction, scrubbing, stripping
Pre-built flowsheet templates
What You’ll Learn#
Access REE element and extractant databases
Calculate pH-dependent distribution coefficients
Simulate multi-stage extraction units
Perform sensitivity analysis and optimization
Use automatic differentiation for gradients
Read the provenance of every number before quoting it
A warning before any number in here is copied#
Since #270 all five extractant records have distribution coefficients fit to
named primary sources, and each record states its fit basis, its correction
arithmetic and its validity window beside the numbers. What is still
HAND_TUNED, and tagged as such, is the temperature_coefficients block on
D2EHPA, PC88A and Cyanex272 — every source behind the refit is isothermal, so
none of them says anything about dH. Only five of the fifteen element prices
carry a citation (USGS Mineral Commodity Summaries 2026); the rest are
estimates. Section 9 shows how to ask any field where it came from. Run the
gradients on all of it; put a dH-sensitive or price-sensitive number in a
paper only with its provenance beside it.
import jax
import jax.numpy as jnp
from jax import grad
# Enable 64-bit precision for numerical stability
jax.config.update("jax_enable_x64", True)
# Import difflow_ree components
from difflow_ree import (
# Database access
get_element,
get_extractant,
list_ree_elements,
list_extractants,
# Distribution model
REEDistribution,
# Unit operations
REEExtractor,
REEExtractorParams,
)
# Provenance: every database field can say where it came from
from difflow_ree.provenance import explain
# Import stream utilities
from difflow.streams import make_stream, get_flows
1. Exploring the REE Database#
The difflow_ree module includes a comprehensive database of REE properties.
# List all available REE elements
elements = list_ree_elements()
print(f"Available REE elements ({len(elements)}):")
print(" " + ", ".join(elements))
# Get detailed properties for specific elements
print("\nElement Properties:")
print("="*78)
print(f"{'Symbol':<8} {'Name':<15} {'Group':<10} {'Price ($/kg)':<15} {'Price source':<15}")
print("-"*78)
for symbol in ["La", "Nd", "Eu", "Dy"]:
elem = get_element(symbol)
prov = explain("elements", f"elements.{symbol}.price_usd_kg")
print(f"{elem.symbol:<8} {elem.name:<15} {elem.group:<10} "
f"${elem.price_usd_kg:<14.2f} {prov.cls + ' / ' + prov.source:<15}")
# Prices are oxide (REO) basis, 2025 annual averages where sourced. Five of the
# fifteen elements have a citation; Dy, Tb and the rest do not, and the numbers
# beside them are estimates. Do not read the dated citation on La as saying
# anything about Dy.
Available REE elements (15):
La, Ce, Pr, Nd, Sm, Eu, Gd, Tb, Dy, Y, Ho, Er, Tm, Yb, Lu
Element Properties:
==============================================================================
Symbol Name Group Price ($/kg) Price source
------------------------------------------------------------------------------
La Lanthanum light $1.00 REFERENCE / USGS26
Nd Neodymium light $69.00 REFERENCE / USGS26
Eu Europium middle $27.00 REFERENCE / USGS26
Dy Dysprosium heavy $450.00 ESTIMATED / EST
2. Exploring Extractant Database#
Industrial extractants for REE separation.
# List available extractants
extractants = list_extractants()
print(f"Available extractants ({len(extractants)}):")
for ext_name in extractants:
print(f" • {ext_name}")
# Get D2EHPA properties
d2ehpa = get_extractant("D2EHPA")
print(f"\nD2EHPA Properties:")
print(f" Full name: {d2ehpa.full_name}")
print(f" Formula: {d2ehpa.formula}")
print(f" MW: {d2ehpa.molecular_weight:.2f} g/mol")
print(f" pKa: {d2ehpa.pKa}")
print(f" Type: {d2ehpa.extractant_type}")
print(f" Typical concentration: {d2ehpa.typical_concentration} M")
print(f" Valid pH range: {d2ehpa.valid_ph_range}")
print(f" Cost: ${d2ehpa.cost_usd_kg}/kg")
Available extractants (5):
• D2EHPA
• PC88A
• Cyanex272
• TBP
• naphthenic_acid
D2EHPA Properties:
Full name: Di-2-ethylhexyl phosphoric acid
Formula: (C8H17O)2PO2H
MW: 322.43 g/mol
pKa: 3.24
Type: acidic_phosphoric
Typical concentration: 0.5 M
Valid pH range: (0.0, 2.0)
Cost: $8.0/kg
3. Distribution Coefficient Calculations#
Distribution coefficients (D) determine how REEs partition between aqueous and organic phases:
\[D = \frac{[\text{REE}]_{\text{organic}}}{[\text{REE}]_{\text{aqueous}}}\]
The model includes:
pH dependence: \(\log_{10}(D) = a + b \cdot pH + c \cdot pH^2\)
Temperature correction
Extractant concentration effects
# Create distribution model for D2EHPA
dist = REEDistribution(
extractant="D2EHPA",
elements=("La", "Nd", "Dy"),
concentration=0.5, # 0.5 M in organic phase
)
# Calculate D values at pH 1.0 -- inside D2EHPA's refitted window, [0, 2]
D_values = dist.get_D_all(pH=1.0, T=298.15)
print("Distribution Coefficients at pH 1.0, 25°C:")
print("="*50)
for elem, D in D_values.items():
print(f" D({elem}) = {float(D):8.4f}")
# Calculate separation factors
SF_Nd_La = float(D_values["Nd"] / D_values["La"])
SF_Dy_Nd = float(D_values["Dy"] / D_values["Nd"])
print(f"\nSeparation Factors:")
print(f" SF(Nd/La) = {SF_Nd_La:.2f}")
print(f" SF(Dy/Nd) = {SF_Dy_Nd:.2f}")
print("\n💡 Higher SF = easier separation")
Distribution Coefficients at pH 1.0, 25°C:
==================================================
D(La) = 10.1205
D(Nd) = 115.6112
D(Dy) = 13128.0433
Separation Factors:
SF(Nd/La) = 11.42
SF(Dy/Nd) = 113.55
💡 Higher SF = easier separation
4. pH Effect on Distribution#
D values increase exponentially with pH for acidic extractants.
# Effect of pH on distribution coefficients
pH_values = [0.0, 0.5, 1.0, 1.5, 2.0] # D2EHPA's fitted window, end to end
print("pH Effect on Distribution Coefficients (D2EHPA, 0.5 M):")
print("="*60)
print(f"{'pH':<6} {'D(La)':<12} {'D(Nd)':<12} {'D(Dy)':<12}")
print("-"*60)
for pH in pH_values:
D_vals = dist.get_D_all(pH=pH, T=298.15)
print(f"{pH:<6.1f} {float(D_vals['La']):<12.4f} {float(D_vals['Nd']):<12.4f} {float(D_vals['Dy']):<12.4f}")
# How steep is that, really? Read it off the table rather than asserting it.
print("\n📊 D rises by 10**b per pH unit, where b is the number of protons the")
print(" exchange releases. Measured across the rows above:")
for lo, hi in [(0.0, 1.0), (1.0, 2.0)]:
for elem in ["La", "Nd", "Dy"]:
ratio = float(dist.get_D(elem, hi) / dist.get_D(elem, lo))
print(f" D({elem}) x {ratio:8.1f} from pH {lo} to {hi}")
print("""
Cation exchange on a dimeric acidic extractant,
RE(3+) + 3 (HA)2 <-> RE(HA2)3 + 3 H+,
releases three protons, so mass action demands b = 3 exactly -- a factor of
1000 per pH unit, not 10. Every ratio above is 1000.0, for every element,
because the #270 refit pinned b = 3 on all four acidic records rather than
fitting a slope per element. It used to read ~300x here: D2EHPA's slope was
HAND_TUNED at b = 2.45, below the stoichiometric one, which is a
thermodynamic inconsistency dressed as a fit. One consequence shows up in
section 8 -- with a single shared slope, log10(beta) = a_i - a_j and the
separation factors are exactly pH-independent.""")
pH Effect on Distribution Coefficients (D2EHPA, 0.5 M):
============================================================
pH D(La) D(Nd) D(Dy)
------------------------------------------------------------
0.0 0.0101 0.1156 13.1280
0.5 0.3200 3.6559 415.1452
1.0 10.1205 115.6112 13128.0433
1.5 320.0369 3655.9479 415145.1800
2.0 10120.4541 115611.2242 13128043.2857
📊 D rises by 10**b per pH unit, where b is the number of protons the
exchange releases. Measured across the rows above:
D(La) x 1000.0 from pH 0.0 to 1.0
D(Nd) x 1000.0 from pH 0.0 to 1.0
D(Dy) x 1000.0 from pH 0.0 to 1.0
D(La) x 1000.0 from pH 1.0 to 2.0
D(Nd) x 1000.0 from pH 1.0 to 2.0
D(Dy) x 1000.0 from pH 1.0 to 2.0
Cation exchange on a dimeric acidic extractant,
RE(3+) + 3 (HA)2 <-> RE(HA2)3 + 3 H+,
releases three protons, so mass action demands b = 3 exactly -- a factor of
1000 per pH unit, not 10. Every ratio above is 1000.0, for every element,
because the #270 refit pinned b = 3 on all four acidic records rather than
fitting a slope per element. It used to read ~300x here: D2EHPA's slope was
HAND_TUNED at b = 2.45, below the stoichiometric one, which is a
thermodynamic inconsistency dressed as a fit. One consequence shows up in
section 8 -- with a single shared slope, log10(beta) = a_i - a_j and the
separation factors are exactly pH-independent.
5. Multi-Stage Extraction Unit#
Simulate a counter-current extraction cascade.
The organic solvent stream names its carriers explicitly — the extractant (D2EHPA) and the diluent (kerosene) — rather than a generic "Organic" species. The extractant molar flow is what sets the loading capacity of the organic phase (#191), and a solvent that names neither carrier is now rejected instead of being silently assigned an organic flow of 1.0 (#192).
# Create extraction unit parameters
params = REEExtractorParams(
n_stages=5,
extractant="D2EHPA",
elements=("La", "Nd", "Dy"),
# pH=None would take D2EHPA's own default -- the top of its window (#270).
# 0.3 is chosen instead because it is where this cascade is interesting:
# Dy is quantitatively extracted, Nd about two thirds, La barely at all.
pH=0.3,
)
# Create extractor
extractor = REEExtractor(params)
# Define feed streams
feed = make_stream(
flows={
"H2O": 10.0,
"La": 0.01,
"Nd": 0.02,
"Dy": 0.01,
},
T=298.15,
P=101325.0,
)
# Organic solvent stream.
#
# The stream must name the *extractant* and the *diluent* as species. A
# solvent whose carrier matches neither now raises instead of silently
# defaulting the organic flow to 1.0 (#192), and when loading is enabled the
# extractant molar flow is what sets the capacity of the organic phase
# (capacity = F_extractant / m, #191).
#
# 0.5 M D2EHPA in kerosene is roughly 10 mol% extractant (kerosene is
# ~0.75 g/mL and ~170 g/mol, so ~4.4 mol/L of diluent against 0.5 mol/L of
# extractant). The total organic flow is unchanged at 8.0 mol/s; it is just
# split 0.8 D2EHPA / 7.2 kerosene.
solvent = make_stream(
flows={
"D2EHPA": 0.8,
"kerosene": 7.2,
"La": 0.0,
"Nd": 0.0,
"Dy": 0.0,
},
T=298.15,
P=101325.0,
)
# Run extraction
raffinate, extract, info = extractor(feed, solvent)
# Analyze results
feed_flows = get_flows(feed)
raff_flows = get_flows(raffinate)
ext_flows = get_flows(extract)
print(f"Extraction Results (5 stages, pH {params.pH}):")
print("="*70)
print(f"{'Element':<10} {'Feed':<12} {'Raffinate':<12} {'Extract':<12} {'Recovery %':<12}")
print("-"*70)
for elem in ["La", "Nd", "Dy"]:
feed_val = float(feed_flows[elem])
raff_val = float(raff_flows[elem])
ext_val = float(ext_flows[elem])
recovery = (ext_val / feed_val) * 100
print(f"{elem:<10} {feed_val:<12.4f} {raff_val:<12.4f} {ext_val:<12.4f} {recovery:<12.1f}")
# Calculate purities
total_REE_extract = sum(float(ext_flows[e]) for e in ["La", "Nd", "Dy"])
nd_purity = float(ext_flows["Nd"]) / total_REE_extract * 100
print(f"\nExtract Composition:")
print(f" Nd purity (among REEs): {nd_purity:.1f}%")
Extraction Results (5 stages, pH 0.3):
======================================================================
Element Feed Raffinate Extract Recovery %
----------------------------------------------------------------------
La 0.0100 0.0094 0.0006 6.4
Nd 0.0200 0.0063 0.0137 68.3
Dy 0.0100 0.0000 0.0100 100.0
Extract Composition:
Nd purity (among REEs): 56.2%
6. Automatic Differentiation: Gradients#
Compute exact gradients for sensitivity analysis and optimization.
def nd_recovery(pH):
"""Nd recovery vs pH, by the same equation the REEExtractor solves."""
dist_pH = REEDistribution(
extractant="D2EHPA",
elements=("La", "Nd", "Dy"),
concentration=0.5,
)
D_val = dist_pH.get_D("Nd", pH, T=298.15)
S_F = 0.8 # solvent / feed ratio, as in section 5
N = 5.0 # stages, as in section 5
# Counter-current cascade -- Kremser:
#
# fraction extracted = (E**(N+1) - E) / (E**(N+1) - 1), E = D * S/F
#
# This cell used to compute 1 - 1/(1 + E)**N instead. That is the
# CROSS-CURRENT result -- N stages each contacted with FRESH solvent -- and
# it is not a small difference. Counter-current reuses one solvent stream,
# so every stage after the first meets organic that is already partly
# loaded.
#
# E = 1 is a removable singularity (the limit is N/(N+1)). This cell stays
# away from it; difflow_ree guards it with a jnp.where.
E = D_val * S_F
E_Np1 = E ** (N + 1.0)
recovery = (E_Np1 - E) / (E_Np1 - 1.0)
return recovery
# Same operating point as section 5.
pH_operating = 0.3
dR_dpH = grad(nd_recovery)(pH_operating)
R = float(nd_recovery(pH_operating))
# Where does E pass through 1? E = D(Nd) * S/F = 1 means
# log10 D(Nd) = -log10(S/F), and log10 D = a + 3 pH on this record, so
# the crossing is analytic.
coef = get_extractant("D2EHPA").ph_coefficients["Nd"]
pH_E1 = (-jnp.log10(0.8) - coef.a) / coef.b
print("Sensitivity Analysis:")
print("="*50)
print(f"Operating pH: {pH_operating}")
print(f"Extraction factor E = D(Nd)*S/F: {float(dist.get_D('Nd', pH_operating)) * 0.8:.4f}")
print(f"Nd recovery: {R*100:.2f}% (compare section 5's cascade above)")
print(f"\n\u2202(Nd recovery)/\u2202pH = {float(dR_dpH):.4f} per pH unit")
print(f"\n\U0001f4a1 Interpretation:")
print(f" \u2022 The derivative is local. 0.1 pH units up is worth about")
print(f" {float(dR_dpH)*0.1*100:.1f} points of recovery here, but the curve is an S,")
print(f" so do not extrapolate it a whole pH unit -- check by evaluating:")
for pH in [0.2, 0.3, 0.4, 0.5]:
print(f" pH {pH:.1f} recovery {float(nd_recovery(pH))*100:6.2f}%")
print(f" \u2022 pH control is critical: E passes through 1 at pH {float(pH_E1):.3f}, which is")
print(f" exactly where the cascade is most sensitive to it. Three protons per")
print(f" Nd(III) means one tenth of a pH unit is a factor of two in D.")
Sensitivity Analysis:
==================================================
Operating pH: 0.3
Extraction factor E = D(Nd)*S/F: 0.7347
Nd recovery: 68.52% (compare section 5's cascade above)
∂(Nd recovery)/∂pH = 3.5872 per pH unit
💡 Interpretation:
• The derivative is local. 0.1 pH units up is worth about
35.9 points of recovery here, but the curve is an S,
so do not extrapolate it a whole pH unit -- check by evaluating:
pH 0.2 recovery 36.66%
pH 0.3 recovery 68.52%
pH 0.4 recovery 94.78%
pH 0.5 recovery 99.69%
• pH control is critical: E passes through 1 at pH 0.345, which is
exactly where the cascade is most sensitive to it. Three protons per
Nd(III) means one tenth of a pH unit is a factor of two in D.
7. Optimization: Find Optimal pH#
Use gradient descent to maximize Nd purity.
def single_stage_nd_purity(pH):
"""Nd purity among all REEs after ONE equilibrium contact."""
dist_pH = REEDistribution(
extractant="D2EHPA",
elements=("La", "Nd", "Dy"),
concentration=0.5,
)
D_vals = dist_pH.get_D_all(pH, T=298.15)
feed_ratio = {"La": 0.01, "Nd": 0.02, "Dy": 0.01}
extract_amounts = {e: D_vals[e] * feed_ratio[e] for e in feed_ratio}
total = sum(extract_amounts.values())
return extract_amounts["Nd"] / (total + 1e-10)
pH_lo, pH_hi = get_extractant("D2EHPA").valid_ph_range
# First: look at the objective before optimizing it.
print("Single equilibrium contact, over D2EHPA's validity window:")
print("="*60)
print(f"{'pH':<8} {'Nd purity %':<16} {'d(purity)/dpH':<18}")
print("-"*60)
for pH in [0.0, 0.5, 1.0, 1.5, 2.0]:
print(f"{pH:<8.1f} {float(single_stage_nd_purity(pH))*100:<16.4f} "
f"{float(grad(single_stage_nd_purity)(pH)):<18.3e}")
print("""
The column is CONSTANT and the derivative is zero to machine precision. That is
not a numerical accident and it is not a bug: since the #270 refit every element
on this record shares one pH slope, b = 3, so every D moves by the same factor
per pH unit and their ratios -- which is all a single-stage purity is -- cannot
move at all. There is nothing here to optimize. Before #270 this same cell
climbed to the low-pH corner of the window, and the thing it was climbing was
a set of staggered, hand-tuned slopes with no thermodynamic basis.
pH still buys purity, but only through a CASCADE, where each element saturates
at a different rate. That objective is below, and it has a real interior
optimum.""")
# --- The cascade objective, which does depend on pH ------------------------
def cascade_nd_purity(pH):
"""Nd purity in the extract of the 5-stage counter-current cascade."""
dist_pH = REEDistribution(
extractant="D2EHPA",
elements=("La", "Nd", "Dy"),
concentration=0.5,
)
feed_ratio = {"La": 0.01, "Nd": 0.02, "Dy": 0.01}
S_F, N = 0.8, 5.0
extracted = {}
for elem, f in feed_ratio.items():
E = dist_pH.get_D(elem, pH, T=298.15) * S_F
E_Np1 = E ** (N + 1.0)
extracted[elem] = f * (E_Np1 - E) / (E_Np1 - 1.0)
total = sum(extracted.values())
return extracted["Nd"] / (total + 1e-30)
print("\n5-stage counter-current cascade:")
print("="*74)
print(f"{'pH':<8} {'Nd purity %':<14} {'d(purity)/dpH':<16} {'Nd recovery %':<14}")
print("-"*74)
for pH in [0.0, 0.2, 0.4, 0.6, 0.8, 1.0]:
print(f"{pH:<8.1f} {float(cascade_nd_purity(pH))*100:<14.2f} "
f"{float(grad(cascade_nd_purity)(pH)):<16.4f} "
f"{float(nd_recovery(pH))*100:<14.2f}")
# Projected gradient ascent on the cascade objective.
pH_opt = 0.1
learning_rate = 0.2
print("\nProjected gradient ascent:")
print("="*50)
print(f"{'Iter':<8} {'pH':<10} {'Nd Purity %':<15}")
print("-"*50)
for i in range(200):
gradient = grad(cascade_nd_purity)(pH_opt)
pH_opt = jnp.clip(pH_opt + learning_rate * gradient, pH_lo, pH_hi)
if (i + 1) % 40 == 0:
print(f"{i+1:<8} {float(pH_opt):<10.4f} "
f"{float(cascade_nd_purity(pH_opt))*100:<15.2f}")
final_purity = cascade_nd_purity(pH_opt)
final_grad = float(grad(cascade_nd_purity)(pH_opt))
at_bound = bool(abs(float(pH_opt) - pH_lo) < 1e-9 or abs(float(pH_opt) - pH_hi) < 1e-9)
print(f"\nConverged pH: {float(pH_opt):.4f} (window [{pH_lo}, {pH_hi}])")
print(f" Nd purity: {float(final_purity)*100:.2f}%")
print(f" Nd recovery: {float(nd_recovery(pH_opt))*100:.2f}%")
print(f" gradient at the solution: {final_grad:+.2e} per pH unit")
print(f" at a bound: {at_bound}")
print("""
This one IS stationary -- interior, with the gradient at zero rather than
pinned against a bound, and that is the check worth making every time. The
mechanism is worth naming too: Dy is quantitatively extracted across the whole
window, so raising pH adds La to the extract faster than it adds Nd, while
lowering it loses Nd before it loses Dy. The maximum sits where those two
losses balance.
It is still purity alone, and purity alone is never the whole objective. Read
the recovery column beside it: the optimum gives up a few points of Nd recovery
for its purity, and whether that trade is right is an economic question this
cell cannot answer. The industrial answer is not a cleverer single pH at all --
it is to separate La/Nd from Dy in one contactor and Nd from La in another,
which is what a cascade with scrub and strip sections is for (notebook 04).
""")
Single equilibrium contact, over D2EHPA's validity window:
============================================================
pH Nd purity % d(purity)/dpH
------------------------------------------------------------
0.0 1.7295 8.936e-11
0.5 1.7295 2.826e-12
1.0 1.7295 8.933e-14
1.5 1.7295 2.787e-15
2.0 1.7295 9.046e-17
The column is CONSTANT and the derivative is zero to machine precision. That is
not a numerical accident and it is not a bug: since the #270 refit every element
on this record shares one pH slope, b = 3, so every D moves by the same factor
per pH unit and their ratios -- which is all a single-stage purity is -- cannot
move at all. There is nothing here to optimize. Before #270 this same cell
climbed to the low-pH corner of the window, and the thing it was climbing was
a set of staggered, hand-tuned slopes with no thermodynamic basis.
pH still buys purity, but only through a CASCADE, where each element saturates
at a different rate. That objective is below, and it has a real interior
optimum.
5-stage counter-current cascade:
==========================================================================
pH Nd purity % d(purity)/dpH Nd recovery %
--------------------------------------------------------------------------
0.0 15.50 0.8976 9.25
0.2 41.53 1.5930 36.66
0.4 62.69 0.1302 94.78
0.6 57.11 -0.5240 99.99
0.8 50.19 -0.0531 100.00
1.0 50.00 -0.0001 100.00
Projected gradient ascent:
==================================================
Iter pH Nd Purity %
--------------------------------------------------
40 0.4188 62.80
80 0.4188 62.80
120 0.4188 62.80
160 0.4188 62.80
200 0.4188 62.80
Converged pH: 0.4188 (window [0.0, 2.0])
Nd purity: 62.80%
Nd recovery: 96.75%
gradient at the solution: +6.66e-16 per pH unit
at a bound: False
This one IS stationary -- interior, with the gradient at zero rather than
pinned against a bound, and that is the check worth making every time. The
mechanism is worth naming too: Dy is quantitatively extracted across the whole
window, so raising pH adds La to the extract faster than it adds Nd, while
lowering it loses Nd before it loses Dy. The maximum sits where those two
losses balance.
It is still purity alone, and purity alone is never the whole objective. Read
the recovery column beside it: the optimum gives up a few points of Nd recovery
for its purity, and whether that trade is right is an economic question this
cell cannot answer. The industrial answer is not a cleverer single pH at all --
it is to separate La/Nd from Dy in one contactor and Nd from La in another,
which is what a cascade with scrub and strip sections is for (notebook 04).
8. Comparing Different Extractants#
Compare D2EHPA, PC88A and Cyanex272 for Nd extraction — each inside its own validity window, because they no longer share one.
extractant_list = ["D2EHPA", "PC88A", "Cyanex272"]
# The three acidic extractants have no pH in common. Since the #270 refit
# each window is the span its own primary source measured: D2EHPA [0, 2],
# PC88A [0.1, 2.5], Cyanex 272 [1.5, 3.5]. Their intersection is empty, so
# there is no single pH at which the three can honestly be compared. This cell
# used to compare them all at pH 3.0, which extrapolated two of the three.
#
# Compare each at its OWN equal-split pH instead: the pH where D(Nd) = 1, where
# a 1:1 cascade puts half the Nd in each phase. That is the operating point the
# reagent is actually chosen around, and it is inside every window by
# construction.
def equal_split_pH(dist, lo, hi, elem="Nd"):
"""Bisect for the pH in [lo, hi] where D(elem) = 1."""
f = lambda p: float(jnp.log10(dist.get_D(elem, p, T=298.15)))
if f(lo) * f(hi) > 0:
return None
for _ in range(60):
mid = 0.5 * (lo + hi)
if f(lo) * f(mid) <= 0:
hi = mid
else:
lo = mid
return 0.5 * (lo + hi)
print("Each extractant at its own equal-split pH (D(Nd) = 1):")
print("="*86)
print(f"{'Extractant':<12} {'window':<14} {'conc M':<8} {'pH':<7} "
f"{'SF(Nd/La)':<11} {'SF(Dy/Nd)':<11} {'coefficients':<12}")
print("-"*86)
for ext_name in extractant_list:
ext = get_extractant(ext_name)
lo, hi = ext.valid_ph_range
dist_compare = REEDistribution(
extractant=ext_name,
elements=("La", "Nd", "Dy"),
concentration=ext.typical_concentration,
)
pH_eq = equal_split_pH(dist_compare, lo, hi)
D_vals = dist_compare.get_D_all(pH=pH_eq, T=298.15)
prov = explain("extractants", f"extractants.{ext_name}.ph_coefficients.Nd.a")
print(f"{ext_name:<12} [{lo}, {hi}]".ljust(27)
+ f"{ext.typical_concentration:<8.2f} {pH_eq:<7.3f} "
f"{float(D_vals['Nd']/D_vals['La']):<11.1f} "
f"{float(D_vals['Dy']/D_vals['Nd']):<11.1f} {prov.cls:<12}")
print("""
\U0001f4a1 What this table does and does not say:
* D2EHPA and PC88A are within 15% of each other on Nd/La, and Cyanex 272
beats both by a factor of five. Earlier versions of this notebook asserted
a clear PC88A advantage on light-REE selectivity; there is no such thing on
these records. PC88A's real industrial advantage over D2EHPA is not
selectivity, it is stripping: the weaker acid (pKa 4.10 vs 3.24) gives back
its loaded metal at far lower acid strength, which is why it displaced
D2EHPA for the middle and heavy separations (Nash 1993).
* Every separation factor in this table is pH-, temperature- and
concentration-INDEPENDENT: since #270 all four acidic records carry one
shared slope, b = 3, so log10(beta) = a_i - a_j exactly. Recompute any row
at any pH inside its own window and it will not move. That is also why the
choice of comparison pH above changes the D values but not this table.
* The equal-split pH column now runs the way acid strength says it should:
the strongest acid extracts at the lowest pH. It did not before #270,
when two of the three records' coefficients were HAND_TUNED. The last
column is the provenance of each `a`; ask for it before comparing absolute
D values across records.
* Cyanex 272 does have the largest Dy/Nd factor, and it is the reagent of
choice where a heavy/middle cut is what you need.""")
Each extractant at its own equal-split pH (D(Nd) = 1):
======================================================================================
Extractant window conc M pH SF(Nd/La) SF(Dy/Nd) coefficients
--------------------------------------------------------------------------------------
D2EHPA [0.0, 2.0] 0.50 0.312 11.4 113.6 DERIVED
PC88A [0.1, 2.5] 0.50 1.026 10.0 486.5 MEASURED
Cyanex272 [1.5, 3.5] 0.30 2.135 51.4 126.8 MEASURED
💡 What this table does and does not say:
* On these records D2EHPA has the BETTER light-REE selectivity, not PC88A.
Earlier versions of this notebook asserted the opposite. PC88A's real
industrial advantage over D2EHPA is not selectivity, it is stripping: the
weaker acid (pKa 4.10 vs 3.24) gives back its loaded metal at far lower
acid strength, which is why it displaced D2EHPA for the middle and heavy
separations (Nash 1993).
* Every separation factor in this table is pH-, temperature- and
concentration-INDEPENDENT: since #270 all four acidic records carry one
shared slope, b = 3, so log10(beta) = a_i - a_j exactly. Recompute any row
at any pH inside its own window and it will not move. That is also why the
choice of comparison pH above changes the D values but not this table.
* The equal-split pH column now runs the way acid strength says it should:
the strongest acid extracts at the lowest pH. It did not before #270,
when two of the three records' coefficients were HAND_TUNED. The last
column is the provenance of each `a`; ask for it before comparing absolute
D values across records.
* Cyanex 272 does have the largest Dy/Nd factor, and it is the reagent of
choice where a heavy/middle cut is what you need.
9. Where Did These Numbers Come From?#
Every field in the database carries a provenance class. Ask before you quote.
from difflow_ree.provenance import coverage
# The whole-database picture: how many fields of each class.
print("Provenance coverage across difflow_ree:")
for cls, n in sorted(coverage().items(), key=lambda kv: -kv[1]):
print(f" {cls:<12} {n:>4}")
# And any single field, in full.
print()
print(explain("extractants", "extractants.D2EHPA.ph_coefficients.Nd.b"))
Provenance coverage across difflow_ree:
CONVENTION 200
REFERENCE 188
MEASURED 125
DERIVED 50
HAND_TUNED 30
ESTIMATED 18
CONSTRUCTED 6
extractants:extractants.D2EHPA.ph_coefficients.Nd.b = 3.0
source Z1 [MEASURED]
citation Zhang, Jack; Zhao, Baodong; Schreiner, Bryan (2016). "Separation
Hydrometallurgy of Rare Earth Elements". Springer International
Publishing, Cham.
locus Eqs. (4.27)-(4.31), p. 124
note PINNED at 3.0, the stoichiometric slope this record declares
(protons_released: 3). Not fitted -- there is no pH series behind
this record to fit it to. Q1 Fig. 2.8 verifies slope 3 for 1 mol/L
P204 up to lg[H+] = -1. One slope for every element is what makes a
separation factor independent of pH, and it is also what lets
PPH63's separation factors be read directly as differences in `a`.
about Z1: Sec. 4.7. Table 4.36 is the ONLY retrieved full-series measured distribution
table for naphthenic acid, and is the sole source of every per-
element `a` coefficient on the naphthenic_acid record. It is
currently UNREPLICATED -- a second measured table is the highest-
value addition to this database. SINCE #270 IT IS ALSO THE ROUTE TO
CYANEX 272's SHAPE, as a secondary carrying ZL93: Sec. 4.3.1 (p.
138) states the conditions and Table 4.25 (p. 139) prints
separation factors against La for all fourteen stable lanthanides
plus Y. Two of its model statements are used throughout the acidic
records and are quoted on them: Eq. (4.11), free dimer =
C_nominal/2 - 3*[RE(HA2)3], and Eq. (4.12), the proton balance
[H+]_eq = [H+]_i + 3([RE3+]_i - x), which is what licenses
recovering an unreported equilibrium pH in an unbuffered batch
contact. Eqs. (4.27)-(4.31) (p. 124) are the dimer-basis cation-
exchange model itself, and p. 124 states that Y falls between Ho
and Er. Table 4.19 (p. 127, Chai 1998 = C98) is used ONLY as a
discordant cross-check on D2EHPA and is not the source of any
shipped number; Table 4.18 (p. 125) is rejected outright as non-
monotone. Text-extraction caveat: in this PDF `&` renders the minus
sign, and the Eu/Tm cell of Table 4.25 prints 8.0 for 82.0.
Summary#
This notebook demonstrated:
Database access — REE properties and extractant data
Distribution models — pH-dependent D values, and how steep they really are
Multi-stage extraction — counter-current cascades, by Kremser
Automatic differentiation — exact gradients for sensitivity
Optimization — and how to recognise a boundary solution when you get one
Extractant comparison — each inside its own validity window
Provenance — asking any field where it came from
Three things worth carrying out of here:
Counter-current is not cross-current.
(E**(N+1) - E)/(E**(N+1) - 1), not1 - 1/(1+E)**N; the second is roughly twice the first at the conditions in section 5.A converged-looking optimizer is not a converged optimizer. Section 7’s ascent stops at a bound with a non-zero gradient, and says so.
Only PC88A’s pH coefficients are measured, and only five of the fifteen prices are sourced. Section 9 is how you find that out for any field.
Next Steps:
See
21_custom_extractants.ipynbto learn how to define custom extractantsExplore pre-built flowsheet templates for complete separation trains
Add economic analysis to evaluate process profitability