Ankur Gupta ← Back to home

Rethinking How We Teach Ion Transport: Balance Sheets vs. Poisson–Nernst–Planck

Grace Origer, Ritu R. Raj and Ankur Gupta

An interactive companion to our paper: G. Origer, R. R. Raj, N. Jarvey, P. N. Romero Zavala, W. A. Smith and A. Gupta, “Rethinking balance sheets: A Poisson-Nernst-Planck based approach for modeling concentration and flux profiles inside an electrochemical cell,” Chemical Engineering Science, 2026, 124600, doi:10.1016/j.ces.2026.124600.

TL;DR: The "balance sheet" approach in Ch. 4 of Bard, Faulkner & White (BFW) is the first thing many electrochemistry students are taught about ion transport. It gives you a convenient approach to build a model, but the reasoning does not always hold up: it assumes flat concentration profiles in the bulk, then introduces diffusion only at the electrodes, an inconsistency that breaks species conservation. The Poisson–Nernst–Planck (PNP) equations resolve this by predicting the full concentration and flux profiles from first principles. The widgets below let you solve the PNP equations in your browser and overlay the balance sheet predictions, so you can see exactly where the textbook story works, where it doesn't, and why. The full paper is published in Chemical Engineering Science, and the source codes (Python and MATLAB) reproducing each figure are provided with the supplementary material of the publication.

Why this matters

If you have taken a chemical engineering course on electrochemical systems, chances are you have seen something that looks like the figure below: a one-dimensional cell with a cathode on one side, an anode on the other, and, alongside it, a tidy bookkeeping table that splits the current into a "diffusive" part near the electrodes and an "electromigrative" part everywhere else. That table (below) is the balance sheet, and the recipe behind it is the balance sheet analysis (BSA). It is the entry point that BFW gives in Ch. 4 to a beginner trying to make sense of mass transfer in an electrochemical cell.

This post has a personal origin. Some of us took a course in which the balance sheet was taught exactly this way, and came away genuinely confused about where the diffusive contribution goes in the bulk and why the numbers stop balancing there. That confusion turned out not to be a student failing: when we later put the same questions to an expert in the field, they were puzzled by the construction too. That is what prompted us to work the examples out properly.

The balance sheet for the copper cell at $k=6$ electrons per unit time, in atomic transport rates.
SpeciesElectromigration
(bulk)
Diffusion
(near electrode)
Total
(at electrode)
Cu(II)−1−5−6
Cu(I)−1+7+6
Cl⁻+3−30

The balance sheet is appealing because it requires only arithmetic. But the cost of that simplicity is a set of assumptions that are never really justified, and some of them turn out to be wrong even for the very examples BFW uses to illustrate the method. The goal of this post, and of the paper it is based on, is to ask whether the balance sheet analysis (BSA) holds up when you solve the transport equations from first principles, and to give you a hands-on tool to see for yourself where the two pictures agree and where they part ways.

The cell

We focus on the same cell that BFW does: a one-dimensional copper redox cell. The electrolyte contains three ions — copper(II), copper(I), and chloride, in a background electrolyte of ammonia. At the left electrode (the cathode, $\tilde{x}=-1$) Cu(II) is reduced to Cu(I); at the right electrode (the anode, $\tilde{x}=+1$) the reverse happens. Chloride is inert (it does not react at either electrode), but it is required to keep the solution electroneutral.

Cathode x̃ = −1 Anode x̃ = +1 Cu(II) + e⁻ → Cu(I) Cu(I) → Cu(II) + e⁻ 2+ + 2+ + Cu(II) Cu(I) Cl⁻
A 1D copper redox cell. Cu(II) is consumed at the cathode and produced at the anode (so its flux at both walls is fixed by the surface reaction). Cu(I) does the reverse. Cl⁻ is inert. Throughout the cell, the three species' concentrations must satisfy electroneutrality: $2\tilde{C}_\mathrm{Cu(II)} + \tilde{C}_\mathrm{Cu(I)} - \tilde{C}_\mathrm{Cl^-} = 0$.

Following BFW we make one simplifying assumption that keeps the algebra clean: the limiting equivalent conductances of the three ions are equal. With $\lambda_g = (F^2/RT)|z_g|D_g$, this fixes the diffusivity ratios:

$$\mathscr{D}_\mathrm{Cu(II)} = \tfrac{1}{2}, \quad \mathscr{D}_\mathrm{Cu(I)} = 1, \quad \mathscr{D}_\mathrm{Cl^-} = 1.$$
Dimensionless diffusivities scaled by a reference $D$. Cu(II) is divalent, so its diffusivity is half that of the monovalents under the equal-conductance assumption.

The single dimensionless control knob in this problem is the imposed current density:

$$\mathcal{J} = \frac{\nu \ell}{F D C_0},$$
where $\nu = ke/A_\text{electrode}$ is the charge flux at the electrode, $\ell$ is the half-cell width, $D$ the reference diffusivity, and $C_0$ the initial Cu(II) concentration. Throughout, time is non-dimensionalised by $\tau = tD/\ell^2$ and length by $\tilde{x} = x/\ell$.

The balance sheet, briefly

The recipe BFW follows is short:

  1. Assume no concentration gradients in the bulk. Then the current there must be carried entirely by electromigration.
  2. Split the current using transference numbers. Each species' contribution to the bulk electromigrative current is $i_{m,g} = t_g\, i$, where $t_g = z_g^2 D_g C_g / \sum_g z_g^2 D_g C_g$.
  3. At the electrodes, force the imposed reaction rate. Whatever current the electromigration didn't account for is assigned to a "diffusive" boundary contribution.

For our cell at $k=6$ electrons per unit time, the recipe spits out the famous BFW table: an electromigrative atomic transport rate of $\{-1,-1,+3\}$ for $\{$Cu(II), Cu(I), Cl$^-\}$ everywhere in the bulk, plus a diffusive rate of $\{-5,+7,-3\}$ confined to a thin region near each electrode.

This works, in the sense that the numbers add up to the imposed reaction rate at the wall. But it has problems. We highlight five.

  1. Where does "no bulk gradient" come from? No mixing, no specific time scale, and no analysis of the actual transport equations is invoked to justify it.
  2. The diffusive flux is discontinuous between the bulk (zero) and the electrode (non-zero). That discontinuity is precisely what breaks species conservation: a flux that jumps has a non-zero divergence, so at steady state the species cannot be conserved across the cell.
  3. The electromigrative flux is taken to be constant right up to the electrode — despite the very same approach assigning a concentration gradient there.
  4. There is no notion of a limiting current. The recipe gives a number for any imposed $k$, even values that would require a negative concentration somewhere.
  5. There are no concentration profiles at all. So there is no way to check, from the recipe, whether the assumptions it relies on are even approximately satisfied.

The PNP equations resolve all five.

The PNP equations

The Poisson–Nernst–Planck (PNP) framework is just three building blocks: a species balance for each ion, a flux law (Nernst–Planck) that includes both diffusion and electromigration, and a constraint on the electric field (Poisson). For a thin-double-layer system, the regime relevant here, Poisson reduces to electroneutrality.

Start from the statement that each species is conserved:

$$\frac{\partial \tilde{C}_g}{\partial \tau} = -\frac{\partial \tilde{N}_g}{\partial \tilde{x}},$$
Conservation of species $g$: concentration changes only because flux enters or leaves.

and note that the flux $\tilde{N}_g$ has exactly two contributions:

$$\tilde{N}_g = \underbrace{-\mathscr{D}_g \frac{\partial \tilde{C}_g}{\partial \tilde{x}}}_{\textstyle \text{diffusion}} \;\underbrace{-\, z_g \mathscr{D}_g \tilde{C}_g \frac{\partial \tilde\phi}{\partial \tilde{x}}}_{\textstyle \text{electromigration}}.$$
The Nernst–Planck flux. Every panel in the widgets below is one of these two terms, or their sum. This is the split the balance sheet is trying to estimate by arithmetic.

Putting the two together gives the equation we actually solve:

$$\frac{\partial \tilde{C}_g}{\partial \tau} = \mathscr{D}_g \frac{\partial^2 \tilde{C}_g}{\partial \tilde{x}^2} + z_g \mathscr{D}_g \frac{\partial}{\partial \tilde{x}}\!\left(\tilde{C}_g\, \frac{\partial \tilde\phi}{\partial \tilde{x}}\right).$$
Species balance with diffusive and electromigrative fluxes.
$$\sum_g z_g \tilde{C}_g = 0, \qquad \frac{\partial \tilde\phi}{\partial \tilde{x}} \;=\; -\frac{-\mathcal{J} + \sum_g z_g \mathscr{D}_g\, \partial \tilde{C}_g/\partial \tilde{x}}{\sum_g z_g^2 \mathscr{D}_g \tilde{C}_g}.$$
Electroneutrality plus an explicit expression for the potential gradient.

That second expression is worth a sentence, because it is doing a lot of work. We never solve a separate equation for $\tilde\phi$. Instead, multiply each species balance by $z_g$ and sum over species: electroneutrality makes the left-hand side vanish, which forces the current $\tilde{j} = \sum_g z_g \tilde{N}_g$ to be independent of position. Evaluating it at either electrode fixes its value at the imposed current, and solving that statement for $\partial\tilde\phi/\partial\tilde{x}$ gives the expression above. In other words, the potential gradient is whatever it has to be, at each point, to keep the current uniform.

An aside on a common misconception: electroneutrality does not mean a constant electric field. It is a simplification of Poisson's equation in the thin-double-layer limit, but $\partial\tilde\phi/\partial\tilde{x}$ is generally a function of position because the concentrations $\tilde{C}_g(x)$ vary.

The boundary conditions impose the surface reaction rates: $\tilde{N}_\mathrm{Cu(II)} = -\mathcal{J}$, $\tilde{N}_\mathrm{Cu(I)} = +\mathcal{J}$, $\tilde{N}_\mathrm{Cl^-} = 0$ at both walls. Initial conditions are uniform: $\tilde{C}_\mathrm{Cu(II)}=\tilde{C}_\mathrm{Cu(I)}=1$, $\tilde{C}_\mathrm{Cl^-}=3$. That's the entire problem statement.

Interactive #1 · Steady-state explorer

The widget below solves the steady-state PNP equations in your browser. Drag the current slider $\mathcal{J}$ and watch four things happen: (i) the concentration profiles develop slopes; (ii) the electromigrative flux $\tilde{N}_{m,g}$ becomes spatially varying; (iii) the diffusive flux $\tilde{N}_{d,g}$ is non-zero throughout the cell, not just at the electrodes; and (iv) the total flux $\tilde{N}_g$ stays beautifully flat — as it must at steady state, because each species is conserved.

The dashed lines are the balance sheet predictions. Where they agree with PNP and where they don't is the whole story.

Steady-state PNP solver, with BSA overlay
low current ← capped at $0.95\times\mathcal{J}_\mathrm{lim}$ = 0.66
none $100\times$ excess
Cu(II) Cu(I) Cl⁻ PNP (solid) BSA (dashed)
Top-left: concentration profiles (PNP only — BSA makes no prediction here). Top-right: electromigrative flux. Bottom-left: diffusive flux. Bottom-right: total flux. All fluxes are normalised by $|\mathcal{J}|$. Move the current slider toward its right edge and watch the Cu(II) concentration plunge to zero at the cathode — that's the limiting current. Toggle the supporting-electrolyte slider to add NaClO$_4$ at up to $100\times$ excess; the inert ions take over the electromigration.

Pushing toward the limiting current

Without concentration profiles, one thing the BSA cannot do is tell you when its own answer stops being physically attainable. Push the current slider above to about $\mathcal{J}\approx 0.6$. Above this current the Cu(II) concentration profile gets steeper and steeper at the cathode, until it is driven to zero. At that point the surface reaction is starving for reactant: ions cannot be supplied to the wall fast enough to sustain the imposed current. This is the limiting current, and for this cell (with no supporting electrolyte) it sits at $\mathcal{J}_\mathrm{lim}\approx 0.7$.

Two things to notice:

Benefits of a supporting electrolyte

In real electrochemistry, you often add a large excess of an inert salt (a "supporting electrolyte" — here, NaClO$_4$). Drag the $R$ slider to $100$. Two things happen:

So the BSA is not wrong everywhere. It comes closest to being right when you have a strong supporting electrolyte and you are looking at the trace reactive species. The difficulty is that BFW presents the balance sheet as a general recipe, and applies it to cells (including the copper cell above, which has no supporting electrolyte) where that condition is not met.

Interactive #2 · The transient story

Another important regime is early times. What happens before the surface reactions have had a chance to perturb the bulk? There the concentration profiles really are nearly flat, just as the BSA assumes. So perhaps the BSA is a snapshot of an early-time state?

It is, sort of. The widget below solves the time-dependent PNP equations starting from the uniform initial condition. Hit play and watch the profiles evolve. We track three things: the concentration field, the total flux $\tilde{N}_g(\tilde{x})$, and the current density $\tilde j(\tilde{x}) = \sum_g z_g \tilde{N}_g$. At every time, current is a constant in $\tilde{x}$ (electroneutrality enforces it), but the individual fluxes are only constant once steady state is reached.

Transient PNP — copper redox cell at $\mathcal{J} = 1/6$
Cu(II) Cu(I) Cl⁻ BSA reference
Left: $\tilde{C}_g(\tilde{x})$. Middle: total flux $\tilde{N}_g(\tilde{x})/|\mathcal{J}|$ with the BSA "bulk + boundary spike" pattern overlaid. Right: $\tilde{j}(\tilde{x})/|\mathcal{J}|$ — should be uniformly $-1$ at all times by electroneutrality. Useful checkpoints: $\tau \approx 0.05$ (bulk gradients haven't developed yet — close to BSA assumption), $\tau \approx 0.25$ (gradients propagating in), $\tau \approx 5$ (fully steady).

At $\tau \approx 0.05$, the bulk concentration profiles, and therefore the total flux profiles, are nearly flat. The BSA would then say "in the bulk, only electromigration matters," and at this single time slice it's a defensible statement. By $\tau \approx 0.25$ the gradients have invaded the bulk, and by $\tau \approx 5$ the system has settled into the steady-state profile from the previous widget.

The conclusion is worth stating plainly: the BSA is an early-time approximation (and even there it gets some of the details wrong), even though it is presented (in BFW and elsewhere) as a general analysis. The transient solution makes that explicit.

Interactive #3 · The hydrogen evolution cell

The third example in BFW is HCl electrolysis. Here both the cation (H$^+$) and the anion (Cl$^-$) react: H$^+$ is reduced to H$_2$ at the cathode, Cl$^-$ is oxidised to Cl$_2$ at the anode. The diffusivity ratio is large — H$^+$ diffuses about $4\times$ faster than Cl$^-$ — so the picture differs from the copper cell.

Because both species are electroactive, the salt concentration drops near both electrodes. For early times (before the depletion fronts from the two electrodes meet in the middle), there is a beautiful closed-form solution: a similarity-variable expression involving the complementary error function. We use it here as Interactive #3.

$$\theta(\eta) = \eta\,\mathrm{erfc}\!\left(\tfrac{\eta}{2}\right) - \frac{2}{\sqrt\pi}\,\exp\!\left(-\tfrac{\eta^2}{4}\right), \qquad \eta = \frac{1+\tilde x}{\sqrt{D_a \tau}}$$
Similarity variable solution for the salt concentration near the cathode. $D_a = 2 D_\mathrm{H} D_\mathrm{Cl}/(D_\mathrm{H}+D_\mathrm{Cl})$ is the ambipolar diffusivity. The full cell solution is the superposition of two such solutions, one anchored at each electrode.
H-cell similarity-variable solution at early times
just after start fronts about to meet
H⁺ Cl⁻ Salt $\tilde{C}_\mathrm{S}$ BSA
Left: salt concentration $\tilde{C}_\mathrm{S}(\tilde{x})$, showing the depletion fronts growing inward from each electrode. Middle: H$^+$ flux (electromigrative + diffusive). Right: Cl$^-$ flux (electromigrative + diffusive). The BSA dashed lines are the textbook bulk plus boundary-spike pattern in $k=10$ units. The faster H$^+$ carries most of the current by electromigration, while Cl$^-$ does it by diffusion. The similarity solution is valid only while the depletion fronts have not yet met, so $\tau \lesssim 0.15$ here.

Two features are worth pausing on, because they mirror what the first two widgets showed. First, look at the flat middle of the salt profile at small $\tau$: the interior has not heard from the electrodes yet, so the bulk really is gradient-free and the BSA numbers there are close to right. Slide $\tau$ up and watch that flat region shrink from both sides. Once the two fronts approach each other, the bulk the BSA relies on has effectively disappeared.

Second, compare the solid and dashed curves near the walls. The BSA puts a sharp step in the diffusive flux at the edge of its notional boundary layer and holds the electromigrative flux constant right through it. The similarity solution instead turns over smoothly, and the split between the two mechanisms changes continuously across the depletion zone. The asymmetry between the two panels is physical rather than numerical: $\mathscr{D}_\mathrm{H} = 4\mathscr{D}_\mathrm{Cl}$, so the slower Cl$^-$ needs a steeper gradient to carry the same reaction rate, which is why its front is sharper and its diffusive flux does relatively more of the work.

Take-aways for students and instructors

1. The BSA gets the boundary values right and the interior wrong.
At the electrodes the BSA enforces the imposed reaction rate, which is a boundary condition it cannot get wrong. Everything between the walls is where it slips. In the total-flux panel of Interactive #1 (bottom right) the PNP curves are flat, as conservation demands, while the BSA dashed lines jump at the edge of its notional boundary layer. In the diffusive-flux panel (bottom left) the PNP curves are non-zero all the way across the cell, whereas the BSA sets them to zero in the bulk.
2. There is a regime where the BSA comes close: trace reactive species in a sea of supporting electrolyte.
Drag the $R$ slider to $100$ in Interactive #1 and look at the Cu(II)/Cu(I) traces in the electromigrative-flux panel (top right). They flatten, and the BSA assumption of an essentially uniform concentration becomes a good approximation, because the inert ions buffer the field and the reactive species' contribution to it is negligible. This is also why supporting electrolytes are added in practice: they make the diffusive picture self-consistent. Even here, though, the agreement is not exact: the PNP solution still returns different electromigrative fluxes, because electroneutrality has to be satisfied at every point.
3. Early times are another regime where the BSA is approximately right, with some details missing.
Hit play on Interactive #2 and pause around $\tau \approx 0.05$: the concentration panel is nearly flat in the middle, and there the BSA numbers are close. Two things it still misses. Near the walls the PNP total flux turns over smoothly, while the BSA imposes an abrupt step at the edge of its boundary layer (visible in the middle panel, where the dashed lines jump and the solid ones do not). And the BSA has no time in it at all: it returns the same numbers at $\tau = 0.05$, $0.25$ and $5$, whereas the PNP profiles keep evolving until the gradients fill the cell. It is a snapshot presented as if it were the whole film.
4. PNP gives you something BSA cannot: the limiting current.
Crank the current slider in Interactive #1 toward its cap and watch the Cu(II) trace in the concentration panel (top left) dive toward zero at the cathode. Past that point the imposed current is unphysical: the cell cannot supply Cu(II) fast enough. Note too that the cap moves as you add supporting electrolyte, from $\mathcal{J}_\mathrm{lim}\approx 0.69$ at $R=0$ down to $\approx 0.5$ at large $R$, because the reactive ion loses its electromigrative assist and has to arrive by diffusion alone. The BSA never tells you any of this, because it has no concentration field to look at.

If you teach this material, we hope you can see that the BSA is a fine zero-level approximation, but that it should be presented alongside its assumptions and the conditions under which they are or are not satisfied. If you are learning this material, hopefully the widgets above help you understand the material better.

Under the hood

All three solvers run entirely in your browser. The steady state is a Newton solve on a discretised form of the steady PNP equations, the transient is explicit finite-difference time stepping, and the H-cell uses the closed-form similarity solution. Both sliders are capped where the in-browser solvers stop being reliable: the current at $0.95\times\mathcal{J}_\mathrm{lim}$ (which itself falls as you add supporting electrolyte), and the H-cell time at $\tau = 0.15$, where the depletion fronts meet and the similarity solution stops being valid. To go past either cap, use the Python and MATLAB solvers provided with the supplementary material of the publication, which is what produced every figure in it.

Acknowledgements. RR thanks the National Science Foundation (DGE-2040434) Graduate Research Fellowship for financial support. AG thanks the NSF CAREER award #2238412 for financial support. The authors thank Adam Holewinski for his useful input on the manuscript.

A note on AI. This blog post, including its text and the interactive widgets, was created with the help of AI. AI was not used to write the manuscript text; there, it was used only to help organise our code.