Balance Sheets vs. Poisson–Nernst–Planck: An Interactive Walk-through

Ankur Gupta

An interactive companion to a paper led by Grace Origer with Ritu R. Raj, Nathan Jarvey, P. N. Romero Zavala and Wilson A. Smith.

Heads up: this post is a work in progress. Text, widgets, and figures are still being polished — expect rough edges, occasional typos, and the odd misbehaving slider. If you find an issue, please let me know.

TL;DR: The "balance sheet" approach in Ch. 4 of Bard, Faulkner & White (BFW) is the first thing every electrochemistry student is taught about ion transport. It gets you a number, but it gets it for the wrong reason: 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 available at the journal and the source codes (Python and MATLAB) reproducing each figure are linked below.

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 a tidy bookkeeping table that splits the current into a "diffusive" part near the electrodes and an "electromigrative" part everywhere else. That table is the balance sheet. 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.

The balance sheet is appealing because it requires only arithmetic. But the price of the simplicity is a set of assumptions that are never really justified. 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 the paper it is based on — is not to dunk on the balance sheet (which is a perfectly fine first pass), but to show where it breaks, and to give you a hands-on tool to see for yourself how a first-principles calculation differs.

You don't need to take our word for any of it. Drag the sliders, watch the curves move, and compare.

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 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.

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:

Assume no concentration gradients in the bulk. Then the current there must be carried entirely by electromigration.
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$.
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.

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.
The diffusive flux is discontinuous between the bulk (zero) and the electrode (non-zero). Each species is not conserved across the cell.
The electromigrative flux is taken to be constant right up to the electrode — despite the very same approach assigning a concentration gradient there.
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.
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, in one screen

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. Putting these together for a 1D cell:

$$\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, obtained by combining the species balances with $j = \mathrm{const} = -\nu$.

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.

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.

Current density $\mathcal{J}$ 0.17

low current ← capped at $\sim 0.95\,\mathcal{J}_\mathrm{lim}$

Supporting electrolyte ratio $R$ 0

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

One thing the BSA cannot do is tell you when its own answer becomes nonsense. Push the current slider above to about $\mathcal{J}\approx 0.6$. The Cu(II) concentration profile gets steeper and steeper at the cathode, until it crosses through 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:

The slider here is capped at about $0.95\,\mathcal{J}_\mathrm{lim}$. We picked the cap so the in-browser Newton solver is stable; in the paper we go right up to the limit and beyond, which is where the BSA picture really falls apart.
The BSA dashed lines for the electromigrative flux do not change as you change $\mathcal{J}$ — in normalised flux units, BSA is a current-independent prediction. The PNP solid lines, in contrast, develop large gradients near the electrodes as you approach the limit. The discrepancy near the cathode is the BSA's blind spot.

What a supporting electrolyte does (and why people use one)

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:

The reactive species' electromigrative fluxes collapse to zero. Look at the Cu(II), Cu(I), and Cl$^-$ traces in the electromigration panel: they almost flatline. The inert Na$^+$ and ClO$_4^-$ now carry essentially all of the electromigration. Their large concentrations dominate the denominator of $\partial\tilde\phi/\partial\tilde{x}$, screening the field and shutting down electromigration of the trace reactive species.
The diffusive flux of Cu(II) and Cu(I) does the work. What was previously a mix of diffusion and electromigration becomes essentially pure diffusion for the reactive species. This is the regime where the textbook "diffusion-only" picture finally becomes a good approximation — and notice how, in this limit, the BSA's predictions for the reactive ions get much closer to the PNP curves.

So the BSA is not wrong everywhere. It is essentially right when you have a strong supporting electrolyte and you are looking at the trace species. The issue is that BFW presents the analysis without that caveat, in cells where the assumption is not satisfied.

What about early times? Before the surface reactions have had a chance to perturb the bulk, the concentration profiles really are nearly flat — just like BSA assumes. Maybe 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.

Time $\tau$ 0.000

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 really is nearly flat — you can see it in the concentration plot. 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.

Read this carefully: the BSA is an early-time approximation, even though it is presented (in BFW and elsewhere) as a general analysis. The transient solution makes that explicit.

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.

Time $\tau$ 0.05

just after start fronts about to meet

H⁺ Cl⁻ Salt $\tilde{C}_\mathrm{S}$ BSA

Left: salt concentration $\tilde{C}_\mathrm{S}(\tilde{x})$ — the depletion fronts grow inward from each electrode. Middle: H$^+$ flux (electromigrative + diffusive + total). Right: Cl$^-$ flux. The BSA dashed lines are the textbook bulk + boundary-spike pattern in $k=10$ units. The faster H$^+$ carries most of the current via electromigration; Cl$^-$ does it via diffusion. The similarity solution is valid only while the depletion fronts have not yet met — so $\tau \lesssim 0.15$ here.

Take-aways for students and instructors

1. The BSA is right about the destination, often wrong about the journey.
At the electrodes, the BSA enforces the imposed reaction rate — that's a boundary condition you can't get wrong. But its claim of zero-gradient bulk and constant electromigration is generically violated, and the diffusive contribution does not vanish in the interior of a steady-state cell.

2. There is a regime where BSA works: trace reactive species in a sea of supporting electrolyte.
Drag the $R$ slider to $100$ in the steady-state widget and look at the Cu(II)/Cu(I) electromigrative fluxes. They flatten. The BSA assumption that the concentration is essentially uniform 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.

3. Early times are another regime where BSA is approximately right.
Hit play on the transient widget and pause around $\tau \approx 0.05$. The bulk is nearly flat. The BSA "bulk + boundary spike" decomposition is a one-time-slice picture that captures this moment.

4. PNP gives you something BSA cannot: the limiting current.
Crank the current slider toward $\mathcal{J}\approx 0.65$. The Cu(II) profile dives toward zero at the cathode. Past that point the imposed current is unphysical; the cell cannot supply Cu(II) fast enough. The BSA never tells you this, because it has no concentration field to look at.

If you teach this material, the punchline we hope to leave you with is: the BSA is a fine entry-level mnemonic, but it should be presented alongside its assumptions and the conditions under which they are or aren't satisfied. If you are learning this material, hopefully the widgets above let you see for yourself where the textbook story is loyal to the physics and where it isn't.

Under the hood

All three solvers run entirely in your browser — no server, no Python kernel.

Steady state. We discretise the integral-deviation form of the steady PNP equations on a uniform grid (Eq. 28 in the paper), then use Newton's method with continuation in $\mathcal{J}$ to march from the trivial solution at $\mathcal{J}=0$ to the requested current. The Jacobian is computed by forward differences and the linear system is solved with dense LU + partial pivoting. For $\sim 30$ grid points and $3$–$5$ ions, a slider drag completes in under $\sim 200$ ms.
Transient. Explicit finite-difference time stepping (FTCS) of the species balance, with $\partial\tilde\phi/\partial\tilde{x}$ recomputed at every step from the instantaneous concentration field. To keep the slider responsive, we precompute a frame buffer at startup so the playback is just animation. Stable for $\Delta\tau \lesssim h^2/(2\max\mathscr{D}_g)$.
H-cell similarity solution. Closed form. We approximate $\mathrm{erfc}$ with a degree-7 Chebyshev expansion accurate to $\sim 10^{-6}$. The widget evaluates the analytical formula at every $\tilde x$.

The original Python (steady state and h-cell) and MATLAB (transient) solvers used to produce all the figures in the paper are released alongside it. Those are the recommended tools if you want to extend the analysis — the in-browser solvers here are designed for pedagogy, not production.

A note on the slider caps. The current slider in widget #1 is capped at about $0.95\,\mathcal{J}_\mathrm{lim} \approx 0.65$ because the Newton solver becomes finicky near the limit (the Jacobian becomes ill-conditioned as $\tilde{C}_\mathrm{Cu(II)}\to 0$ at the cathode). The h-cell time slider in widget #3 is capped at $\tau = 0.15$ because the underlying similarity solution is only valid while the depletion fronts from the two electrodes haven't yet met. Past these caps the in-browser solvers misbehave; the full paper handles both cases properly.

The paper. Origer, Raj, Jarvey, Romero Zavala, Smith, Gupta, "Rethinking Balance Sheets: A Poisson-Nernst-Planck Based Approach for Modeling Concentration and Flux Profiles Inside an Electrochemical Cell," Chemical Engineering Journal, 2026 (under revision). Source codes for all figures are bundled with the supplementary material.

Acknowledgements. Major credit to Grace Origer who led the work and to Ritu R. Raj for substantial contributions to the analysis and figures. Thanks also to Nathan Jarvey, P. N. Romero Zavala, and Wilson A. Smith for collaboration on the manuscript. The interactive widgets here were built by AG with Claude as a coding assistant.