REE Separation for NdFeB Magnet Production

REE Separation for NdFeB Magnet Production#

This notebook designs the first separation section of a bastnasite refinery – the cut that takes praseodymium and neodymium away from the lanthanum that makes up most of the feed – and then prices it.

Background#

NdFeB magnets are critical for electric vehicle motors, wind turbine generators, hard drives and consumer electronics. The alloy needs:

  • Nd, the primary magnetic element,

  • Pr, usually left blended with it as didymium rather than separated out,

  • Dy (or Tb), a few percent, for coercivity at temperature.

Bastnasite supplies the first two and almost none of the third. That is a fact about the ore, not about the process, and section 1 makes it explicit rather than assuming a feed that would be convenient.

What this notebook is careful about#

Three things in here are easy to get wrong and are called out where they happen:

  1. One section makes one cut. An extract-scrub-strip circuit separates the feed at a single point on the extraction ladder. It cannot simultaneously reject the elements lighter than the target and the elements heavier than it. Section 4 measures the ceiling that imposes.

  2. PC88A’s correlation is valid over pH 0.1-2.5 and b = 3, so a tenth of a pH unit is a factor of two in every distribution coefficient. Every operating pH here sits inside that window, and section 3 shows what pH actually buys.

  3. The plant scale, the molar flows and the economics are one number. Section 1 fixes the capacity and everything downstream is derived from it.

Objectives#

  1. Model the La/(Pr+Nd) cut with PC88A

  2. Understand what extraction pH does, and does not, control

  3. Analyze sensitivity by automatic differentiation

  4. Estimate economics from the circuit’s own product, not from assumed yields

import jax
import jax.numpy as jnp
from jax import grad

jax.config.update("jax_enable_x64", True)

# Import difflow core
from difflow.streams import make_stream, get_flows

# Import REE plugin
from difflow_ree import (
    # Database
    get_element, list_ree_elements, get_extractant, ExtractantDatabase,
    # Equilibrium
    REEDistribution, get_distribution_coefficient,
    # Units
    REEExtractor, REEExtractorParams,
    # Flowsheets
    ExtractScrubStripCircuit, ExtractScrubStripParams,
    # Economics
    REEPricing, estimate_capex, capex_basis, estimate_opex, calculate_profit,
    ree_oxide_mass_flow,
)
from difflow_ree.provenance import explain

print("REE elements in the database:", len(list_ree_elements()))
print(list_ree_elements())

pc88a = get_extractant("PC88A")
print(f"\nPC88A pH validity range: {pc88a.valid_ph_range}")
print("Every operating pH below is inside it. The correlation is not clamped")
print("outside -- it warns and extrapolates -- and with b = 3 an extrapolation")
print("of one pH unit is three decades in D.")

REE elements in the database: 15
['La', 'Ce', 'Pr', 'Nd', 'Sm', 'Eu', 'Gd', 'Tb', 'Dy', 'Y', 'Ho', 'Er', 'Tm', 'Yb', 'Lu']

PC88A pH validity range: (0.1, 2.5)
Every operating pH below is inside it. The correlation is not clamped
outside -- it warns and extrapolates -- and with b = 3 an extrapolation
of one pH unit is three decades in D.

1. Feed and plant scale#

Everything downstream is derived from two numbers: the composition, which is a property of the ore, and the capacity, which is a choice. Fixing them here is what keeps the molar flows in section 4 and the dollars in section 6 talking about the same plant.

The composition is the Mountain Pass bastnasite REO distribution with cerium already removed – cerium is roughly half the contained REO and is taken out first by oxidation, which is cheaper than extracting it. The remainder renormalises to a feed that is about two-thirds lanthanum.

Read the Dy line before designing around it. Bastnasite is not a dysprosium ore. At 0.06 mol% there is no dysprosium product here at any recovery; the heavy rare earths come from ion-adsorption clays and from monazite/xenotime, not from this feed. Section 6 prices what the circuit actually makes.

# Mountain Pass bastnasite REO distribution with Ce removed, renormalised.
# Mole percent of contained rare earth (not oxide weight percent).
FEED_MOL_PCT = {
    "La": 65.44,
    "Pr": 8.62,
    "Nd": 23.84,   # the target, with Pr
    "Sm": 1.62,
    "Gd": 0.41,
    "Dy": 0.061,   # see the markdown above: this is not a Dy ore
}
ELEMENTS = tuple(FEED_MOL_PCT)

# Plant scale. These two numbers set every flow and every dollar below.
CAPACITY_TPY = 5000.0     # tonnes of contained REO per year
HOURS_PER_YEAR = 8000.0   # 91% availability

# Oxide mass per mole of metal. Prices are $/kg of oxide, so the oxide molar
# mass -- not the atomic weight -- is the conversion. And it is not one formula
# for all of them: the database carries Pr as Pr6O11 and Ce as CeO2 because that
# is how they are sold, so dividing oxide_mw by 2 is wrong for a third of the
# lanthanides. ree_oxide_mass_flow reads the metal-atom count off the formula.
KG_REO_PER_MOL = {e: float(ree_oxide_mass_flow({e: 1.0})) for e in ELEMENTS}
print("oxide basis:", {e: f"{get_element(e).oxide_formula}"
                       f" -> {KG_REO_PER_MOL[e]*1000:.1f} g/mol RE"
                       for e in ELEMENTS})

_total_pct = sum(FEED_MOL_PCT.values())
x_feed = {e: v / _total_pct for e, v in FEED_MOL_PCT.items()}
kg_per_mol_mixed = sum(x_feed[e] * KG_REO_PER_MOL[e] for e in ELEMENTS)

TOTAL_REE_MOL_S = CAPACITY_TPY * 1000.0 / kg_per_mol_mixed / (HOURS_PER_YEAR * 3600.0)
feed_composition = {e: x_feed[e] * TOTAL_REE_MOL_S for e in ELEMENTS}

# Aqueous flow set to make the feed 0.30 M in total rare earth, which is where a
# real chloride liquor sits. Water is 55.5 mol/L, so this is a molar proxy for
# the volumetric flow the circuit's O/A ratios act on.
FEED_MOLARITY = 0.30
FEED_L_PER_S = TOTAL_REE_MOL_S / FEED_MOLARITY
FEED_H2O_MOL_S = FEED_L_PER_S * 55.5

print(f"Plant: {CAPACITY_TPY:.0f} t REO/yr at {HOURS_PER_YEAR:.0f} h/yr")
print(f"Mean oxide mass: {kg_per_mol_mixed*1000:.1f} g REO per mol RE")
print(f"Total rare earth: {TOTAL_REE_MOL_S:.4f} mol/s "
      f"in {FEED_L_PER_S:.2f} L/s at {FEED_MOLARITY:.2f} M\n")

print(f"{'':>4} {'mol%':>7} {'mol/s':>10} {'t REO/yr':>10}  group")
print("-" * 48)
for elem in ELEMENTS:
    flow = feed_composition[elem]
    tpy = flow * KG_REO_PER_MOL[elem] * HOURS_PER_YEAR * 3600.0 / 1000.0
    print(f"{elem:>4} {x_feed[elem]*100:>7.2f} {flow:>10.5f} {tpy:>10.2f}  "
          f"{get_element(elem).group}")
print("-" * 48)
_check = sum(feed_composition[e] * KG_REO_PER_MOL[e] for e in ELEMENTS) \
    * HOURS_PER_YEAR * 3600.0 / 1000.0
print(f"{'':>4} {100.0:>7.2f} {TOTAL_REE_MOL_S:>10.5f} {_check:>10.2f}")

oxide basis: {'La': 'La2O3 -> 162.9 g/mol RE', 'Pr': 'Pr6O11 -> 170.2 g/mol RE', 'Nd': 'Nd2O3 -> 168.2 g/mol RE', 'Sm': 'Sm2O3 -> 174.4 g/mol RE', 'Gd': 'Gd2O3 -> 181.2 g/mol RE', 'Dy': 'Dy2O3 -> 186.5 g/mol RE'}
Plant: 5000 t REO/yr at 8000 h/yr
Mean oxide mass: 165.1 g REO per mol RE
Total rare earth: 1.0516 mol/s in 3.51 L/s at 0.30 M

        mol%      mol/s   t REO/yr  group
------------------------------------------------
  La   65.45    0.68826    3229.09  light
  Pr    8.62    0.09066     444.50  light
  Nd   23.84    0.25074    1214.89  light
  Sm    1.62    0.01704      85.56  middle
  Gd    0.41    0.00431      22.51  middle
  Dy    0.06    0.00064       3.45  heavy
------------------------------------------------
      100.00    1.05165    5000.00

2. Extractant selection: PC88A against D2EHPA#

The usual sentence is “PC88A is preferred for Nd/Pr because it has the higher separation factor.” Until #270 this database said the opposite, and this notebook said so at some length. It no longer does, and the reason the old answer was wrong is worth more than the new answer is.

The comparison has to be made carefully, because the two extractants do not work at the same acidity. PC88A’s validity window is pH 0.1-2.5 and D2EHPA’s, since the #270 refit, is 0.0-2.0; they overlap, but their working points are 0.71 pH units apart, so at any pH where one is doing its job the other is either saturated or asleep. Each is therefore evaluated at its own working point – the pH at which D(Nd) = 1, where a counter-current section on Nd would be designed – and the separation factors are compared there.

Both correlations are log10 D = a + b*pH + c*pH^2, and since #270 both have b pinned at 3, the stoichiometric value for RE3+ + 3(HA)2 <-> RE(HA2)3 + 3H+, with c = 0. A shared slope makes log10 beta_ij = a_i - a_j, so beta is a constant on each record and the whole of the selectivity lives in the intercepts. The cell below prints both columns flat down the pH sweep, which is the visible signature of that.

What changed the answer: D2EHPA’s old coefficients were HAND_TUNED and carried a staggered b (2.30 for La rising to 2.80 for Dy). A staggered slope makes beta rise with pH for no mechanistic reason, and at D2EHPA’s own working point the stagger printed a larger Nd/Pr beta than PC88A’s measured one. That was an artifact of the hand-tuning, and #270 removed it by refitting D2EHPA against named primary sources with the slope pinned. The textbook ordering came back on its own – which is a better argument for the refit than any of the fit statistics.

The rest of the case for PC88A is unchanged and does not rest on beta at all: it works at lower acidity and, much more importantly, it strips at ordinary acidity. D2EHPA holds the middle and heavy rare earths so strongly that stripping them needs 4-6 M acid, which is a reagent bill and a materials problem rather than a separation factor. Section 4 runs into the mild version of that same problem with PC88A and dysprosium.

dist_d2ehpa = REEDistribution(extractant="D2EHPA", elements=ELEMENTS)
dist_pc88a = REEDistribution(extractant="PC88A", elements=ELEMENTS)


def working_pH(extractant, element="Nd"):
    """pH at which D(element) = 1 -- where a section on that element sits.

    Solved from the correlation rather than searched, so no evaluation is made
    outside the extractant's validity window on the way to the answer.
    """
    c = get_extractant(extractant).ph_coefficients[element]
    if c.c == 0.0:
        root = -c.a / c.b
    else:
        disc = (c.b ** 2 - 4.0 * c.c * c.a) ** 0.5
        roots = [(-c.b + disc) / (2.0 * c.c), (-c.b - disc) / (2.0 * c.c)]
        lo, hi = get_extractant(extractant).valid_ph_range
        inside = [r for r in roots if lo <= r <= hi]
        root = inside[0] if inside else min(roots, key=abs)
    return root


pH_pc88a = working_pH("PC88A")
pH_d2ehpa = working_pH("D2EHPA")
print(f"Working point (D(Nd) = 1):  PC88A pH {pH_pc88a:.2f}   "
      f"D2EHPA pH {pH_d2ehpa:.2f}")
print(f"Validity windows:           PC88A {get_extractant('PC88A').valid_ph_range}"
      f"   D2EHPA {get_extractant('D2EHPA').valid_ph_range}")
_at_pc = float(get_distribution_coefficient("Nd", "D2EHPA", pH_pc88a))
_at_d2 = float(get_distribution_coefficient("Nd", "PC88A", pH_d2ehpa))
print(f"The windows overlap over pH 0.1-2.0, but the working points are")
print(f"{pH_pc88a - pH_d2ehpa:.2f} pH units apart, which is 3x that in log D. At PC88A's")
print(f"working point D2EHPA has D(Nd) = {_at_pc:.0f}; at D2EHPA's, PC88A has")
print(f"D(Nd) = {_at_d2:.4f}. A table at one shared pH compares a saturated")
print("extractant against an idle one.\n")

print("Distribution coefficients, each extractant at its own working point")
print("=" * 58)
print(f"{'Element':>8} {'D2EHPA':>14} {'PC88A':>14}")
print(f"{'':>8} {'(pH '+format(pH_d2ehpa,'.2f')+')':>14} "
      f"{'(pH '+format(pH_pc88a,'.2f')+')':>14}")
print("-" * 58)
for elem in ELEMENTS:
    D_d = float(get_distribution_coefficient(elem, "D2EHPA", pH_d2ehpa))
    D_p = float(get_distribution_coefficient(elem, "PC88A", pH_pc88a))
    print(f"{elem:>8} {D_d:>14.4g} {D_p:>14.4g}")

print("\nNd/Pr separation factor, and how it moves with pH inside each window:")
print(f"{'':>10} {'D2EHPA':>10}  {'':>10} {'PC88A':>10}")
print("-" * 46)
for pH_d, pH_p in ((0.5, 0.5), (1.0, 1.0), (1.5, 1.5), (2.0, 2.0)):
    b_d = float(dist_d2ehpa.get_separation_factor("Nd", "Pr", pH_d))
    b_p = float(dist_pc88a.get_separation_factor("Nd", "Pr", pH_p))
    print(f"  pH {pH_d:>4.1f} {b_d:>10.3f}    pH {pH_p:>4.1f} {b_p:>10.3f}")

beta_d = float(dist_d2ehpa.get_separation_factor("Nd", "Pr", pH_d2ehpa))
beta_p = float(dist_pc88a.get_separation_factor("Nd", "Pr", pH_pc88a))
print(f"\nAt their own working points: D2EHPA {beta_d:.3f}, PC88A {beta_p:.3f}.")
print("BOTH columns are constant, because b = 3 is shared on both records:")
print("beta is fixed by the intercepts alone. D2EHPA's column used to rise with")
print("pH, and that hand-tuned stagger is what used to make it 'win'.")
print("Provenance:")
for ext, elem in (("PC88A", "Nd"), ("PC88A", "Pr"), ("D2EHPA", "Nd"), ("D2EHPA", "Pr")):
    rec = explain("extractants", f"{ext}.ph_coefficients.{elem}.a")
    print(f"  {ext:>7} a({elem}) : {rec.cls}/{rec.source}")

Working point (D(Nd) = 1):  PC88A pH 1.03   D2EHPA pH 0.27
Validity windows:           PC88A (0.1, 2.5)   D2EHPA (0.0, 2.0)

The windows overlap over pH 0.1-2.0, but the working points are
0.76 pH units apart, which is 3x that in log D. At PC88A's
working point D2EHPA has D(Nd) = 192; at D2EHPA's, PC88A has
D(Nd) = 0.0052. A table at one shared pH compares a saturated
extractant against an idle one.

Distribution coefficients, each extractant at its own working point
==========================================================
 Element         D2EHPA          PC88A
              (pH 0.27)      (pH 1.03)
----------------------------------------------------------
      La        0.06457        0.09956
      Pr         0.3837         0.4667
      Nd              1              1
      Sm          6.243          5.171
      Gd          38.98          37.76
      Dy          243.3          486.5

Nd/Pr separation factor, and how it moves with pH inside each window:
               D2EHPA                  PC88A
----------------------------------------------
  pH  0.5      2.606    pH  0.5      2.143
  pH  1.0      2.606    pH  1.0      2.143
  pH  1.5      2.606    pH  1.5      2.143
  pH  2.0      2.606    pH  2.0      2.143

At their own working points: D2EHPA 2.606, PC88A 2.143.
BOTH columns are constant, because b = 3 is shared on both records:
beta is fixed by the intercepts alone. D2EHPA's column used to rise with
pH, and that hand-tuned stagger is what used to make it 'win'.
Provenance:
    PC88A a(Nd) : MEASURED/T21
    PC88A a(Pr) : DERIVED/CALC
   D2EHPA a(Nd) : DERIVED/CALC
   D2EHPA a(Pr) : DERIVED/CALC

3. What extraction pH actually controls#

optimal_pH_for_separation is in the API, and “which pH maximises beta” is the wrong question to ask of PC88A. With one shared slope, beta does not depend on pH at all. The search used to maximise that constant array and return whichever grid point argmax reached first, a different answer for a different pH_range. It now detects the flat scan, warns, and returns the pair’s cut point instead (below), the same pH for any range that contains it. The cell below shows both halves of that.

What pH does control is the cut point. An element is extracted when its Kremser extraction factor E = D * (O/A) exceeds one, and since log10 D = a + 3 pH, that happens above

\[\mathrm{pH}_{\mathrm{cut},i} = \frac{-a_i - \log_{10}(O/A)}{3}\]

So the ladder of intercepts becomes a ladder of pH thresholds, and choosing the extraction pH is choosing where on that ladder to cut the feed. Everything above the cut goes into the organic; everything below stays in the raffinate. That is the design decision, and beta only tells you how sharp the cut can be.

print("optimal_pH_for_separation on Nd/Pr, over two different ranges:")
for lo, hi in ((0.1, 2.5), (0.5, 2.0)):
    pH_opt, beta_opt = dist_pc88a.optimal_pH_for_separation(
        element1="Nd", element2="Pr", pH_range=(lo, hi))
    print(f"  range ({lo}, {hi}):  pH* = {pH_opt:.4f}   beta = {beta_opt:.9f}")

_scan = jnp.array([float(dist_pc88a.get_separation_factor("Nd", "Pr", float(p)))
                   for p in jnp.linspace(0.1, 2.5, 100)])
print(f"\nSpread of beta over that scan: {float(_scan.max() - _scan.min()):.2e}")
print("Same beta to nine figures: there is no maximum to find, and the")
print("intercepts already fixed beta at 10**(a_Nd - a_Pr). The pH returned is")
print("the Nd/Pr cut, where sqrt(D_Nd * D_Pr) * (O/A) = 1 (O/A = 1 by default),")
print("and the SeparationFactorFlatWarning above says so.")
print()

O_over_A = 1.2   # the solvent-to-feed ratio used in section 4
print(f"Cut points at O/A = {O_over_A}:  pH where E = D*(O/A) crosses 1")
print("=" * 52)
print(f"{'Element':>8} {'a':>10} {'pH_cut':>10}   at pH 1.15")
print("-" * 52)
coeffs = get_extractant("PC88A").ph_coefficients
for elem in ELEMENTS:
    a_i = coeffs[elem].a
    pH_cut = (-a_i - jnp.log10(O_over_A)) / 3.0
    side = "extracted" if 1.15 > pH_cut else "stays in raffinate"
    print(f"{elem:>8} {a_i:>10.4f} {float(pH_cut):>10.3f}   {side}")
print("-" * 52)
print("The gap between La at 1.333 and Pr at 1.110 is the whole design window:")
print("0.22 pH units wide, and the extraction pH has to land inside it.")

optimal_pH_for_separation on Nd/Pr, over two different ranges:

/var/folders/sj/cb6cr7_s45d67mqjmsfx8mt80000gn/T/ipykernel_47650/92341575.py:3: SeparationFactorFlatWarning: The Nd/Pr separation factor on 'PC88A' does not depend on pH (relative spread 2.3e-15 over [0.1, 2.5]): the returned pH 1.08 is the cut where their geometric-mean D * (O/A) = 1 at O/A 1, not an SF maximum.
  pH_opt, beta_opt = dist_pc88a.optimal_pH_for_separation(
/var/folders/sj/cb6cr7_s45d67mqjmsfx8mt80000gn/T/ipykernel_47650/92341575.py:3: SeparationFactorFlatWarning: The Nd/Pr separation factor on 'PC88A' does not depend on pH (relative spread 2.1e-16 over [0.5, 2]): the returned pH 1.08 is the cut where their geometric-mean D * (O/A) = 1 at O/A 1, not an SF maximum.
  pH_opt, beta_opt = dist_pc88a.optimal_pH_for_separation(

  range (0.1, 2.5):  pH* = 1.0810   beta = 2.142890601
  range (0.5, 2.0):  pH* = 1.0810   beta = 2.142890601

Spread of beta over that scan: 4.88e-15
Same beta to nine figures: there is no maximum to find, and the
intercepts already fixed beta at 10**(a_Nd - a_Pr). The pH returned is
the Nd/Pr cut, where sqrt(D_Nd * D_Pr) * (O/A) = 1 (O/A = 1 by default),
and the SeparationFactorFlatWarning above says so.

Cut points at O/A = 1.2:  pH where E = D*(O/A) crosses 1
====================================================
 Element          a     pH_cut   at pH 1.15
----------------------------------------------------
      La    -4.0795      1.333   stays in raffinate
      Pr    -3.4086      1.110   extracted
      Nd    -3.0776      0.999   extracted
      Sm    -2.3640      0.762   extracted
      Gd    -1.5006      0.474   extracted
      Dy    -0.3905      0.104   extracted
----------------------------------------------------
The gap between La at 1.333 and Pr at 1.110 is the whole design window:
0.22 pH units wide, and the extraction pH has to land inside it.

4. The La/(Pr+Nd) cut, and the ceiling on it#

Extraction at pH 1.15 puts the cut between La and Pr. Praseodymium and everything heavier loads onto the organic; lanthanum mostly does not. Scrubbing at pH 1.05 with a small aqueous stream sends the co-extracted lanthanum back – the scrub liquor is recycled to the feed in a real plant, so the lanthanum and the target metal it carries with it are not lost, they are recirculated. Stripping at pH 0.2 takes the product off the organic – most of it. The next cell reports the barren organic alongside the three aqueous exits, because with b = 3 the same acidity that releases neodymium (D of a few thousandths) is nowhere near enough for dysprosium (D still above one), which leaves on the recycled solvent instead of in the product.

The section cannot reject the heavies. Samarium, gadolinium and dysprosium all sit above the cut, they load more strongly than Nd, and scrubbing at a lower pH strips Nd before it strips them. So Sm + Gd + Dy report to the product no matter how the section is tuned, and they set a hard ceiling on Pr+Nd purity:

\[\text{ceiling} = \frac{x_{\mathrm{Pr}} + x_{\mathrm{Nd}}} {x_{\mathrm{Pr}} + x_{\mathrm{Nd}} + x_{\mathrm{Sm}} + x_{\mathrm{Gd}} + x_{\mathrm{Dy}}}\]

which the next cell evaluates. This is why an industrial refinery is a cascade: one section per cut, and the didymium section is fed by the section that took the heavies out, not by the raw liquor. What is designed here is a rougher, and it should be read as one.

EXTRACTION_PH = 1.15
SCRUBBING_PH = 1.05
STRIPPING_PH = 0.2
S_OVER_F = 1.2
SCRUB_RATIO = 0.3

params = ExtractScrubStripParams(
    extractant="PC88A",
    elements=ELEMENTS,
    target_elements=("Pr", "Nd"),   # didymium
    n_extraction_stages=10,
    n_scrubbing_stages=6,
    n_stripping_stages=5,
    extraction_pH=EXTRACTION_PH,
    scrubbing_pH=SCRUBBING_PH,
    stripping_pH=STRIPPING_PH,
    solvent_to_feed_ratio=S_OVER_F,
    scrub_to_solvent_ratio=SCRUB_RATIO,
    strip_to_solvent_ratio=0.4,
)

circuit = ExtractScrubStripCircuit(params)

feed_flows = {"H2O": FEED_H2O_MOL_S}
feed_flows.update(feed_composition)
feed = make_stream(flows=feed_flows, T=298.15, P=101325.0)

results = circuit(feed, T=298.15)

print("Three-section circuit")
print("=" * 56)
print(f"  Extraction: pH {params.extraction_pH}, "
      f"{params.n_extraction_stages} stages, O/A = {params.solvent_to_feed_ratio}")
print(f"  Scrubbing:  pH {params.scrubbing_pH}, "
      f"{params.n_scrubbing_stages} stages, scrub/O = {params.scrub_to_solvent_ratio}")
print(f"  Stripping:  pH {params.stripping_pH}, "
      f"{params.n_stripping_stages} stages")
print("All three are inside PC88A's validity range "
      f"{get_extractant('PC88A').valid_ph_range}.")

Three-section circuit
========================================================
  Extraction: pH 1.15, 10 stages, O/A = 1.2
  Scrubbing:  pH 1.05, 6 stages, scrub/O = 0.3
  Stripping:  pH 0.2, 5 stages
All three are inside PC88A's validity range (0.1, 2.5).

print("Where each element goes (% of its feed)")
print("=" * 76)
print(f"{'':>5} {'raffinate':>11} {'scrub liquor':>14} {'product':>10} "
      f"{'barren org.':>13} {'closure':>9}")
print("-" * 76)
raff = get_flows(results["raffinate"])
scrubl = get_flows(results["scrub_liquor"])
prod = get_flows(results["product"])
barren = get_flows(results["barren_organic"])
for elem in ELEMENTS:
    f_in = feed_composition[elem]
    r_ = float(raff.get(elem, 0.0)) / f_in
    s_ = float(scrubl.get(elem, 0.0)) / f_in
    p_ = float(prod.get(elem, 0.0)) / f_in
    b_ = float(barren.get(elem, 0.0)) / f_in
    print(f"{elem:>5} {r_*100:>10.1f}% {s_*100:>13.1f}% {p_*100:>9.2f}% "
          f"{b_*100:>12.2f}% {(r_+s_+p_+b_):>9.4f}")
print("-" * 76)
_D_strip = {e: float(get_distribution_coefficient(e, "PC88A", STRIPPING_PH))
            for e in ELEMENTS}
print("The barren-organic column is the one that has to be read. At the strip")
print(f"pH of {STRIPPING_PH}, D_Nd = {_D_strip['Nd']:.4f} so Nd comes off cleanly,"
      f" but D_Gd = {_D_strip['Gd']:.2f}")
print(f"and D_Dy = {_D_strip['Dy']:.2f}:")
print("PC88A holds dysprosium too strongly for a mild strip, so it leaves on")
print("the recycled organic and would accumulate there until it is bled off.")
print("A real circuit strips the heavies separately at much higher acidity.")

print("\nProduct composition (mol%)")
print("-" * 34)
for elem in ELEMENTS:
    print(f"  {elem:>3}: {results['product_purity'][elem]*100:>7.3f}%")

target_purity = results["target_purity"]
ceiling = ((x_feed["Pr"] + x_feed["Nd"])
           / sum(x_feed[e] for e in ("Pr", "Nd", "Sm", "Gd", "Dy")))
print(f"\nPr+Nd in product : {target_purity*100:.2f}%")
print(f"Ceiling set by the heavies above the cut : {ceiling*100:.2f}%")
print(f"Gap to the ceiling : {(ceiling - target_purity)*100:.2f} percentage points"
      "  (this is the La that got through)")

print("\nSingle-pass recovery to product:")
for elem, rec in results["target_recovery"].items():
    print(f"  {elem}: {float(rec)*100:.1f}%")
print("The balance is in the scrub liquor, not lost: recycling it to the feed")
print("is what turns a single-pass number into a circuit recovery, and it is")
print("why the scrub ratio is a purity knob rather than a yield penalty.")

Where each element goes (% of its feed)
============================================================================
        raffinate   scrub liquor    product   barren org.   closure
----------------------------------------------------------------------------
   La       58.0%          41.3%      0.73%         0.00%    1.0000
   Pr        0.1%          36.2%     63.74%         0.00%    1.0000
   Nd        0.0%          16.9%     83.07%         0.00%    1.0000
   Sm        0.0%           3.3%     96.73%         0.00%    1.0000
   Gd        0.0%           0.4%     98.31%         1.24%    1.0000
   Dy        0.0%           0.0%     16.45%        83.51%    1.0000
----------------------------------------------------------------------------
The barren-organic column is the one that has to be read. At the strip
pH of 0.2, D_Nd = 0.0033 so Nd comes off cleanly, but D_Gd = 0.13
and D_Dy = 1.62:
PC88A holds dysprosium too strongly for a mild strip, so it leaves on
the recycled organic and would accumulate there until it is bled off.
A real circuit strips the heavies separately at much higher acidity.

Product composition (mol%)
----------------------------------
   La:   1.730%
   Pr:  19.793%
   Nd:  71.343%
   Sm:   5.645%
   Gd:   1.452%
   Dy:   0.036%

Pr+Nd in product : 91.14%
Ceiling set by the heavies above the cut : 93.95%
Gap to the ceiling : 2.81 percentage points  (this is the La that got through)

Single-pass recovery to product:
  Pr: 63.7%
  Nd: 83.1%
The balance is in the scrub liquor, not lost: recycling it to the feed
is what turns a single-pass number into a circuit recovery, and it is
why the scrub ratio is a purity knob rather than a yield penalty.

5. Sensitivity by automatic differentiation#

A gradient is a local statement. With b = 3, d(log10 D)/d(pH) = 3, so recovery moves fast enough that a tangent line leaves the physically possible range within a tenth of a pH unit. The cell below reports the derivative and walks it out against the model it was taken from, which is the only honest way to quote one: the tangent and the truth agree at 0.02 pH units and have parted company by 0.10.

def nd_recovery_fn(pH, SF_ratio, n_stages=3):
    """Nd recovery from a three-element extractor, as a function of pH and O/A.

    n_stages is held fixed here; REEExtractorParams will accept a float for it
    (Kremser is E**(N+1)), but the point of this section is the two continuous
    levers.
    """
    params = REEExtractorParams(
        n_stages=n_stages,
        extractant="PC88A",
        elements=("La", "Nd", "Dy"),
        pH=pH,
    )
    extractor = REEExtractor(params)

    feed = make_stream(
        flows={"H2O": FEED_H2O_MOL_S, "La": feed_composition["La"],
               "Nd": feed_composition["Nd"], "Dy": feed_composition["Dy"]},
        T=298.15, P=101325.0,
    )
    # The solvent must name the extractant and the diluent as species, and it
    # must carry a real extractant flow: the loading capacity of the organic
    # phase is F_extractant / m (#191), so an extractant flow of 0.0 is zero
    # capacity and hence zero extraction. A stream whose carrier is neither
    # the extractant nor the diluent raises rather than silently defaulting
    # the organic flow to 1.0 (#192).
    #
    # Built the way ExtractScrubStripCircuit builds it, and at the same scale
    # as the real feed: diluent at F_aq * O/A, extractant at 0.5 M of that.
    # Get the scale wrong and the section runs out of loading capacity rather
    # than out of driving force, and the pH derivative changes SIGN -- more
    # acid then pulls La onto an organic that has no room for it and Nd comes
    # back off. That is a real regime, but it is not the one being measured.
    F_org = FEED_H2O_MOL_S * SF_ratio
    solvent = make_stream(
        flows={
            "PC88A": 0.5 * F_org,
            "kerosene": F_org,
            "La": 0.0, "Nd": 0.0, "Dy": 0.0,
        },
        T=298.15, P=101325.0,
    )

    _, extract, _ = extractor(feed, solvent)
    return get_flows(extract)["Nd"] / feed_composition["Nd"]

# pH 1.0 with 3 stages and O/A = 0.5 sits at partial recovery, where a
# derivative is informative. It is inside PC88A's window.
pH_base = 1.0
n_stages_base = 3
SF_base = 0.5

d_rec_d_pH = float(grad(nd_recovery_fn, argnums=0)(pH_base, SF_base, n_stages=n_stages_base))
d_rec_d_SF = float(grad(nd_recovery_fn, argnums=1)(pH_base, SF_base, n_stages=n_stages_base))
base_recovery = float(nd_recovery_fn(pH_base, SF_base, n_stages=n_stages_base))

print("Sensitivity of Nd recovery")
print("=" * 56)
print(f"Base case: pH = {pH_base}, N = {n_stages_base}, O/A = {SF_base}")
print(f"Recovery:  {base_recovery*100:.1f}%\n")
print(f"d(recovery)/d(pH)  = {d_rec_d_pH:.4f} per pH unit")
print(f"d(recovery)/d(O/A) = {d_rec_d_SF:.4f} per unit of O/A\n")

print("The pH tangent against the model it was taken from:")
print(f"{'delta pH':>9} {'linear':>10} {'actual':>10}")
for d_pH in (0.02, 0.05, 0.10, 0.20):
    linear = base_recovery + d_rec_d_pH * d_pH
    actual = float(nd_recovery_fn(pH_base + d_pH, SF_base, n_stages=n_stages_base))
    print(f"{d_pH:>9.2f} {linear*100:>9.1f}% {actual*100:>9.1f}%")
print("Recovery is bounded by 1 and the tangent is not, so the linear column")
print("has to fail; the table says where. Quote the derivative, not")
print("derivative x a finite step.\n")

# Stage count as a finite difference for comparison
rec_n3 = float(nd_recovery_fn(pH_base, SF_base, n_stages=3))
rec_n5 = float(nd_recovery_fn(pH_base, SF_base, n_stages=5))
print(f"Stages 3 -> 5: {rec_n3*100:.1f}% -> {rec_n5*100:.1f}% "
      f"({(rec_n5-rec_n3)/2*100:+.2f}% per stage)")

delta = 0.001
fd = (float(nd_recovery_fn(pH_base + delta, SF_base, n_stages=n_stages_base))
      - base_recovery) / delta
print(f"\nGradient check -- finite difference {fd:.4f}, autodiff {d_rec_d_pH:.4f}")

Sensitivity of Nd recovery
========================================================
Base case: pH = 1.0, N = 3, O/A = 0.5
Recovery:  55.7%

d(recovery)/d(pH)  = 2.8925 per pH unit
d(recovery)/d(O/A) = 0.8375 per unit of O/A

The pH tangent against the model it was taken from:
 delta pH     linear     actual
     0.02      61.5%      61.5%
     0.05      70.2%      70.1%
     0.10      84.6%      82.5%
     0.20     113.5%      96.0%
Recovery is bounded by 1 and the tangent is not, so the linear column
has to fail; the table says where. Quote the derivative, not
derivative x a finite step.

Stages 3 -> 5: 55.7% -> 60.1% (+2.18% per stage)

Gradient check -- finite difference 2.8941, autodiff 2.8925

6. Economics#

The capacity is the one fixed in section 1, and the production is what the circuit in section 4 actually made – not an assumed recovery times an assumed grade. That is the only reason the CAPEX, the OPEX and the revenue below refer to the same plant.

annual_capacity = CAPACITY_TPY  # section 1; the plant is sized on the feed

# Capital cost. `scope` sets the battery limits and is the single most
# consequential argument: the available anchors differ by more than a factor
# of ten at the same capacity because they enclose different amounts of plant.
# "separation_plant" is a standalone refinery -- cascade plus buildings,
# civils, utilities, effluent treatment, engineering and contingency.
basis = capex_basis("separation_plant")
print(f"Anchored to {basis['source']}: ${basis['capex_usd']/1e6:,.0f} M for "
      f"{basis['capacity_tpy']:,.0f} t/yr {basis['capacity_basis']} "
      f"({basis['year']})")
print(f"Method: {basis['derived_class']}, "
      f"{basis['derived_accuracy'][0]:+.0%} / {basis['derived_accuracy'][1]:+.0%}")

capex = estimate_capex(
    annual_ree_tonnes=annual_capacity,
    n_stages_extraction=params.n_extraction_stages,
    n_stages_scrubbing=params.n_scrubbing_stages,
    n_stages_stripping=params.n_stripping_stages,
    include_precipitation=True,
    year=2024,
    scope="separation_plant",
)

print("\nCapital Cost Estimate")
print("=" * 40)
for item, cost in capex.items():
    print(f"{item:20s}: ${cost/1e6:,.2f} M")

print(f"\nTotal CAPEX: ${capex['total']/1e6:,.2f} M")
print("\nOnly the TOTAL is anchored. The line items split it by conventional")
print("section shares -- the anchor published a total, not a breakdown.")

Anchored to AVALON_GEISMAR: $302 M for 10,000 t/yr separated REO (2012)
Method: AACE Class 5 (capacity-factored from one project), -50% / +100%

Capital Cost Estimate
========================================
mixer_settlers      : $55.90 M
tanks_vessels       : $15.24 M
pumps_piping        : $22.87 M
instrumentation     : $17.78 M
precipitation       : $15.24 M
installation        : $33.03 M
civil_structures    : $22.87 M
electrical_utilities: $17.78 M
effluent_treatment  : $10.16 M
engineering         : $20.33 M
contingency         : $38.11 M
total               : $269.31 M

Total CAPEX: $269.31 M

Only the TOTAL is anchored. The line items split it by conventional
section shares -- the anchor published a total, not a breakdown.

opex = estimate_opex(
    annual_ree_tonnes=annual_capacity,
    capex=capex["total"],
    extractant="PC88A",
)

print("Annual Operating Cost")
print("=" * 40)
for item, cost in opex.items():
    print(f"{item:20s}: ${cost/1e6:,.2f} M")

print(f"\nTotal OPEX: ${opex['total']/1e6:,.2f} M/year")

Annual Operating Cost
========================================
extractant          : $10.00 M
acid                : $5.00 M
base                : $2.50 M
precipitant         : $15.00 M
labor               : $1.68 M
utilities           : $2.11 M
maintenance         : $8.08 M
total               : $44.37 M

Total OPEX: $44.37 M/year

pricing = REEPricing()

# Prices are read from elements.yaml, not restated in the economics module --
# a second copy is the drift #268 removed from the Langmuir constants.
print("Oxide prices used, and where each came from:")
for elem in ELEMENTS:
    rec = explain("elements", f"{elem}.price_usd_kg")
    print(f"  {elem:>3}: ${pricing.get_price(elem, '99%', 'oxide'):>8.2f}/kg  "
          f"{rec.cls}/{rec.source}")

print("\nProduct stream, priced as the mixed oxide the section makes")
print("=" * 62)
print(f"{'':>5} {'t REO/yr':>11} {'$/kg':>9} {'$M/yr':>10}")
print("-" * 62)
total_revenue = 0.0
product_tpy_total = 0.0
for elem in ELEMENTS:
    mol_s = float(prod.get(elem, 0.0))
    tpy = mol_s * KG_REO_PER_MOL[elem] * HOURS_PER_YEAR * 3600.0 / 1000.0
    price = pricing.get_price(elem, "99%", "oxide")
    rev = tpy * 1000.0 * price
    total_revenue += rev
    product_tpy_total += tpy
    print(f"{elem:>5} {tpy:>11.2f} {price:>9.2f} {rev/1e6:>10.2f}")
print("-" * 62)
print(f"{'':>5} {product_tpy_total:>11.2f} {'':>9} {total_revenue/1e6:>10.2f}")

# The raffinate is a La stream. It is not free to dispose of, and at $1/kg it
# is not obviously worth purifying either -- but it is contained value, so it
# is worth seeing next to the product.
raff_tpy = sum(float(raff.get(e, 0.0)) * KG_REO_PER_MOL[e]
               for e in ELEMENTS) * HOURS_PER_YEAR * 3600.0 / 1000.0
raff_value = sum(float(raff.get(e, 0.0)) * KG_REO_PER_MOL[e]
                 * pricing.get_price(e, "99%", "oxide")
                 for e in ELEMENTS) * HOURS_PER_YEAR * 3600.0
print(f"\nRaffinate: {raff_tpy:.1f} t REO/yr worth ${raff_value/1e6:.2f} M/yr "
      "of contained oxide,")
print("almost all of it lanthanum. A cerium/lanthanum market in surplus is")
print("why a bastnasite refinery's economics rest on the didymium section.")

print(f"\nTotal Revenue: ${total_revenue/1e6:,.2f} M/year")
print(f"Every stream is priced at its pure-oxide value, but the section makes a"
      f" {target_purity*100:.0f}% Pr+Nd")
print("concentrate, not separated oxides. A 99.5% didymium product would price")
print("higher and needs the next section in the cascade to exist.")

Oxide prices used, and where each came from:
   La: $    1.00/kg  REFERENCE/USGS26
   Pr: $   69.00/kg  REFERENCE/USGS26
   Nd: $   69.00/kg  REFERENCE/USGS26
   Sm: $   15.00/kg  ESTIMATED/EST
   Gd: $   30.00/kg  REFERENCE/USGS26
   Dy: $  450.00/kg  ESTIMATED/EST

Product stream, priced as the mixed oxide the section makes
==============================================================
         t REO/yr      $/kg      $M/yr
--------------------------------------------------------------
   La       23.70      1.00       0.02
   Pr      283.32     69.00      19.55
   Nd     1009.23     69.00      69.64
   Sm       82.76     15.00       1.24
   Gd       22.13     30.00       0.66
   Dy        0.57    450.00       0.26
--------------------------------------------------------------
          1421.71                91.37

Raffinate: 1872.7 t REO/yr worth $1.89 M/yr of contained oxide,
almost all of it lanthanum. A cerium/lanthanum market in surplus is
why a bastnasite refinery's economics rest on the didymium section.

Total Revenue: $91.37 M/year
Every stream is priced at its pure-oxide value, but the section makes a 91% Pr+Nd
concentrate, not separated oxides. A 99.5% didymium product would price
higher and needs the next section in the cascade to exist.

profit = calculate_profit(
    revenue=total_revenue,
    opex=opex["total"],
    capex=capex["total"],
)

print("Profitability Analysis")
print("=" * 40)
print(f"Revenue:      ${profit['revenue']/1e6:>8,.2f} M/year")
print(f"OPEX:         ${profit['opex']/1e6:>8,.2f} M/year")
print(f"EBITDA:       ${profit['ebitda']/1e6:>8,.2f} M/year")
print(f"Depreciation: ${profit['depreciation']/1e6:>8,.2f} M/year")
print(f"Net Income:   ${profit['net_income']/1e6:>8,.2f} M/year")
print(f"\nPayback Period: {profit['payback_years']:.1f} years")
print(f"ROI: {profit['roi']*100:.1f}%")

print()
for _line in (
    "Read these as arithmetic on the inputs, not as an evaluation. The CAPEX",
    "total is anchored to a disclosed project cost, but the revenue above",
    "values every stream at the full oxide price -- a real offtake pays a",
    "fraction of contained value for an unfinished product, and the OPEX unit",
    "rates have no source at all. examples/10 works the same circuit with",
    "payability made explicit, and the sign of the answer changes.",
):
    print(_line)

Profitability Analysis
========================================
Revenue:      $   91.37 M/year
OPEX:         $   44.37 M/year
EBITDA:       $   47.00 M/year
Depreciation: $   26.93 M/year
Net Income:   $   15.05 M/year

Payback Period: 6.4 years
ROI: 5.6%

Read these as arithmetic on the inputs, not as an evaluation. The CAPEX
total is anchored to a disclosed project cost, but the revenue above
values every stream at the full oxide price -- a real offtake pays a
fraction of contained value for an unfinished product, and the OPEX unit
rates have no source at all. examples/10 works the same circuit with
payability made explicit, and the sign of the answer changes.

Summary#

What the notebook built: a 5,000 t REO/yr rougher that cuts a Ce-free bastnasite liquor between La and Pr with PC88A, and the economics of the concentrate it makes.

What is worth carrying out of it:

  1. The plant scale is one number. Composition and capacity were fixed in section 1 and every molar flow, tonne and dollar afterwards was derived from them. A notebook with three scales in it is three notebooks.

  2. beta is not a function of pH when b is shared, which is what stoichiometry says it should be. optimal_pH_for_separation therefore returns a cut point, not an SF maximum; the thing pH really sets is the cut point, pH_cut = (-a - log10(O/A))/3.

  3. One section, one cut. The heavies above the cut cannot be scrubbed out, so the Pr+Nd purity has a ceiling fixed by the feed. Section 4 measures both the ceiling and the distance to it. An industrial cascade is a sequence of these, not a single clever section.

  4. Bastnasite is not a Dy ore. The feed carries 0.06 mol% Dy and no design changes that – and the little there is does not even come off the solvent at a strip pH chosen for neodymium, so it recirculates rather than reporting to a product.

  5. A derivative is local. Section 5 quotes d(recovery)/d(pH) and then walks the tangent out until it leaves the unit interval, which with b = 3 takes about a tenth of a pH unit.

Next steps#

  • examples/10_bastnasite_separation.ipynb – the same feed through a full cascade, with payability made explicit

  • examples/24_custom_ree_elements.ipynb – adding an element, or an element’s coefficients for one extractant

  • examples/32_stochastic_ree_separation.ipynb – what the design looks like when the correlation coefficients are uncertain