REVIEW 4 major objections 5 minor 10 references
Modeling Concentration Profiles in Electrolytes by Solving 3-D Poisson-Nernst-Planck Equations via Finite Difference Method
T0 review · 4 major / 5 minor · reviewed 2026-08-10 · deepseek-v4-flash
Pith's one-line read The paper derives explicit finite-difference stability conditions for the 3-D Poisson-Nernst-Planck equations and shows that, under those conditions, explicit time-stepping of the coupled system is practically infeasible for realistic…
desk verdict The paper's stability criteria don't follow from its own analysis, and without them the central claim of a validated explicit 3D PNP solver collapses. read the letter →
The pith
A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.
The reading
What carries the argument
The machinery is the explicit finite-difference discretization of the 3-D Nernst-Planck and Poisson equations on a uniform cubic grid, together with Fourier-mode (von Neumann) stability analysis. The central objects are the amplification factors for the diffusion part, migration part, and Poisson part, whose combination yields the safe-step inequalities. Then the Poisson equation is removed from the time-stepping loop and solved as an elliptic problem by Gauss-Seidel iteration, which sidesteps the stricter Poisson stability constraint. Boundary conditions are handled with no-flux Neumann conditions on lateral and top faces and the Butler-Volmer electrode reaction at the bottom surface.
What would settle it
Run an explicit 3-D PNP code with parameters satisfying Eqs. (30) and (36) but with a steep potential gradient near the electrode, and check whether oscillations or divergence appear before the nominal stable time horizon; if they do, the additive combination of amplification factors is disproven.
Extended reading notes
Core claim
The central discovery is a closed-form stability criterion for explicit finite-difference solution of the 3-D PNP system. By inserting Fourier modes into the discretized Nernst-Planck equation and treating the electric potential as a constant coefficient in the migration term, the paper obtains an amplification factor whose magnitude must stay below one; this yields Eq. (30), linking diffusion coefficients, time step, grid spacing, and migration strengths. A second Fourier analysis of the discretized Poisson update gives Eq. (36), a constraint on the grid spacing alone. Since the coupled-system amplification factor is taken as the sum of the two individual factors, the combined bound is stricter than either alone. The paper concludes that explicit time-stepping of the fully coupled PNP equations is generally not viable for realistic electrochemical systems, and that solving Poisson as an elliptic equation with Gauss-Seidel is the practical path; with that choice it produces 3-D concentration profiles for single-ion and 1:1 electrolyte systems under finite-layer and Butler-Volmer boundary conditions.
Load-bearing premise
The proof assumes the electric potential can be treated as a fixed coefficient while Fourier-analyzing the migration term, and that the overall amplification factor is the sum of the Nernst-Planck and Poisson factors; if potential gradients vary sharply in space or time, those assumptions break and the safe-step inequalities no longer guarantee stability.
Editorial extensions
If this is right
- Anyone using an explicit finite-difference scheme for 3-D PNP can pre-check the time step and grid against Eqs. (30) and (36) and avoid runs that diverge.
- Realistic battery parameters, such as high concentration, strong electric fields, and large domains, fall outside the stable window, which explains why fully explicit PNP coupling tends to fail in practice.
- Treating the Poisson equation as an elliptic problem rather than time-stepping it is not just a convenience; under the paper's analysis it is the difference between feasible and infeasible simulation.
- The parameter studies give practical guidance: roughly 50 to 80 time steps with moderate grid resolution balance stability and accuracy for the test cases.
- The framework can serve as a benchmark or training-data generator for machine-learning models of ion transport, as the paper itself suggests.
Reading between the lines
- Because the stability derivation freezes the potential inside the migration term, the true stability limit in steep concentration or potential gradients is likely stricter than Eq. (30); a practical rule of thumb would be to stay well inside the predicted window.
- The same amplification-factor construction could be adapted to implicit or semi-implicit time stepping, where the Poisson constraint would likely relax; this is a natural next test the paper does not run.
- The reported instability at longer times, even inside the nominal stable parameter set, suggests the linear stability analysis misses nonlinear feedback between concentration and potential; checking against a highly resolved reference solution would quantify how much the practical window shrinks.
Editorial analysis
A structured set of objections, weighed in public.
Referee Report
Summary. The manuscript proposes an explicit finite-difference scheme for solving the three-dimensional Poisson-Nernst-Planck (PNP) equations with Butler-Volmer boundary conditions. It derives stability inequalities for the Nernst-Planck diffusion and migration terms and for the Poisson equation, and uses these to argue that explicit PNP coupling is severely restricted, leading to practical guidance on time steps and grid sizes. The paper includes simulations of single-ion and two-ion (1:1 electrolyte) systems and discusses the influence of grid resolution and time-step count. Its central deliverable is the set of stability criteria in Eqs. (30), (36), and (43).
Significance. If the stability criteria were correct, they would provide quantitative guidance for choosing time steps and grid sizes in explicit FDM solutions of 3D PNP equations, a practically relevant problem in electrochemical modeling. The paper also demonstrates a step-by-step discretization and applies Butler-Volmer boundary conditions, which could be useful as a starting point for similar implementations. However, the central derivation is built on an invalid Fourier analysis of the migration term and an incorrect additive combination of amplification factors for the coupled system. Because the paper's main conclusions and practical recommendations rest on these steps, the contribution as presented is not reliable. The manuscript also contains internal inconsistencies in the boundary conditions and important presentation gaps, so the claimed 'precise simulation' capability is not substantiated.
major comments (4)
- [§2.1.4, Eqs. (20)–(26)] The stability analysis of the migration term treats the electric potential Φ as a constant while substituting a Fourier mode for the concentration. For a non-constant Φ, the finite-difference migration operator is not translation-invariant, so Fourier modes do not diagonalize the update. If Φ is treated as a constant multiplier, the migration term contributes an imaginary amplification factor i μ_E sin(kΔx), not the real expression written in Eq. (26). Consequently, the real-valued derivation of cases 1 and 2 leading to inequalities (28)–(30) does not govern the actual discretized scheme.
- [§2.1.5.2, Eq. (43)] The total amplification factor for the coupled PNP system is defined additively as λ_total = λ_N−P + λ_Φ. The coupled update is a composition of the Nernst-Planck and Poisson updates; even for scalar factors, |λ1 + λ2| ≤ 1 is neither necessary nor sufficient for stability of the product/composition. Thus the claim that the coupled stability condition is 'more restrictive' and the conclusion that 'it is more impossible to simulate successfully a real electrochemical system' do not follow from the analysis presented.
- [§2.1.4.2, Eqs. (35)–(36)] The Poisson equation is elliptic and does not involve time stepping, yet an 'amplification factor' and a stability condition (36) are derived for it and then used as a constraint on the spatial grid. No derivation shows how this condition relates to the temporal stability of the coupled explicit scheme; absent such a link, Eq. (36) does not support the time-step recommendations made later in Section 2.2.1.
- [§2.1.6, Eqs. (44)–(45)] The boundary conditions for the concentration are internally contradictory: Eq. (44) imposes zero Neumann conditions at the lateral and top boundaries, while Eq. (45) imposes zero Dirichlet conditions (c = 0) at the same lateral boundaries. Both cannot hold simultaneously, and the boundary treatment used in the simulations is therefore not well-defined. Additionally, Eq. (46) mixes the exchange current density J0, concentration gradient, potential gradient, and mobility ν in a single expression whose dimensional consistency is not evident as printed.
minor comments (5)
- [§2.1.5.1, Eq. (40)] Equation (40) is garbled and missing operator symbols, making the Gauss-Seidel update formula hard to parse; it should be rewritten with clear placement of Φ and ρ terms.
- [§2.2.1] The phrase 'the time step size cannot be less than 50' is ambiguous and dimensionally incomplete; a time step size should be expressed with units (e.g., seconds), and this recommendation is not supported by the preceding stability analysis.
- [Figures 3–6] The figures lack axis labels and color scales, and the text does not specify which spatial slice or time instant is shown. This makes it difficult to verify the claimed concentration profiles.
- [§2.1.3.3, Eq. (12)] Equation (12) contains apparent typographical errors, such as a misplaced superscript and the notation 'c_c' and 'cα' mixed with 'cα' superscripts, which should be corrected for clarity.
- [References] The references are inconsistently formatted (some with page ranges, some without), and the citations for the Nernst-Planck and Poisson equations [9,10] are to computational studies rather than to the original equation formulations; the authors should provide canonical references.
Circularity Check
No circularity: the stability conditions are derived from the paper's own discretization, with no fitted parameters, no self-citation chain, and no target quantity defined in terms of its inputs.
full rationale
I walked the derivation chain and found no circular step that satisfies the required standard of exhibiting a specific reduction. The stability inequalities in Eqs. (28)-(30), (35)-(36), and (43) are obtained by substituting Fourier modes into the paper's own explicit finite-difference discretizations of the Nernst-Planck and Poisson equations. The coefficients entering those inequalities are the discretization parameters (diffusion coefficients, grid spacings, time step, charge, Faraday constant, RT), not quantities fitted to the concentration profiles that the paper later claims to predict. The initial concentration, diffusion coefficients, applied potential, exchange current density, and transfer coefficient are all prescribed outright, and the Butler-Volmer boundary condition is imposed rather than inferred from the simulation output. There are no references to prior work by the present authors, so there is no load-bearing self-citation chain and no imported uniqueness theorem. The questionable step noted in the reader's take—treating the potential as a fixed multiplier in the Fourier analysis and adding rather than composing the amplification factors in Eq. (43)—is a mathematical-validity or correctness concern (the migration term may be imaginary and the coupled operator is a product of updates), not circularity, because this step does not assume the very stability bound it derives. Thus the appropriate circularity score is 0.
Assumptions & free parameters
free parameters (4)
- D_z =
1e-2 m^2/s
- J_0 =
1.0 mA/cm^2
- beta =
0.5
- N_t =
50 to 80 (recommended)
assumptions (6)
- standard math Central finite difference approximations for first and second derivatives.
- domain assumption Electrolyte is homogeneous with a uniform initial concentration of 1 M.
- domain assumption Ion-ion interactions beyond the mean-field Poisson coupling are neglected.
- domain assumption No ion flux through the lateral and top boundaries.
- ad hoc to paper The electric potential is treated as a constant coefficient in the Fourier stability analysis of the migration term.
- ad hoc to paper The total amplification factor is the sum of the N-P and Poisson amplification factors.
Cite this review
Pith. "Pith review of Modeling Concentration Profiles in Electrolytes by Solving 3-D Poisson-Nernst-Planck Equations via Finite Difference Method." pith.science (2026). https://pith.science/paper/FGDP7YE6
@misc{pith2026250105917,
author = {Pith},
title = {Pith review of: Modeling Concentration Profiles in Electrolytes by Solving 3-D Poisson-Nernst-Planck Equations via Finite Difference Method},
year = {2026},
howpublished = {\url{https://pith.science/paper/FGDP7YE6}},
note = {Machine review of arXiv:2501.05917}
}
read the original abstract
The Poisson-Nernst-Planck (PNP) equations are fundamental for modeling ion transport in electrochemical systems, capturing the intricate interplay of concentration gradients, electric fields, and ion fluxes essential for applications such as energy storage devices and other electrochemical devices. This study introduces a refined numerical framework employing the finite difference method to solve the 3-D PNP equations, enabling precise simulation of ion concentration distributions under realistic boundary conditions and applied electric fields. By rigorously addressing stability criteria and integrating advanced boundary constraints, including the Butler-Volmer equation for surface reactions, the model provides comprehensive insights into ion dynamics, particularly near electrode surfaces where electric field and reaction effects dominate. This framework significantly enhances traditional PNP modeling by accommodating varied boundary conditions, diffusion anisotropy, and complex electrochemical environments, offering a robust tool for investigating electrochemical processes and guiding the design of advanced electrochemical systems.
Figures
Reference graph
Works this paper leans on
- [1]
-
[2]
T.L. Horng, T.C. Lin, C. Liu, B. Eisenberg, PNP equations with steric effects: a model of ion flow through channels, J Phys Chem B, 116 (2012) 11422-11441
work page 2012
-
[3]
G. Mitscha -Baude, A. Buttinger -Kreuzhuber, G. Tulzer, C. Heitzinger, Adaptive and iterative methods for simulations of nanopores with the PNP–Stokes equations, Journal of Computational Physics, 338 (2017) 452-476
work page 2017
-
[4]
J.S. Newman, C.W. Tobias, Theoretical Analysis of Current Distribution in Porous Electrodes, Journal of The Electrochemical Society, 109 (1962) 1183
work page 1962
-
[5]
J. Ding, Z. Wang, S. Zhou, Positivity preserving finite difference methods for Poisson –Nernst– Planck equations with steric interactions: Application to slit -shaped nanopore conductance, Journal of Computational Physics, 397 (2019) 108864
work page 2019
-
[6]
L.F. Souza, M.T. Mendonça, M.A.F. Medeiros, The advantages of using high -order finite differences schemes in laminar -turbulent transition studies, International Journal for Numerical Methods in Fluids, 48 (2005) 565-582
work page 2005
-
[7]
H. Liu, Z. Wang, A f ree energy satisfying finite difference method for Poisson –Nernst–Planck equations, Journal of Computational Physics, 268 (2014) 363-376
work page 2014
-
[8]
Z. Xu, M. Ma, P. Liu, Self -energy-modified Poisson -Nernst-Planck equations: WKB approximation and finite-difference approaches, Phys Rev E, 90 (2014) 013307
work page 2014
Show all 10 references
-
[9]
Dolatabadi, A
R. Dolatabadi, A. Mohammadi, M. Baghani, A computational simulation of electromembrane extraction based on Poisson - Nernst - Planck equations, Anal Chim Acta, 1158 (2021) 338414
2021
-
[10]
E.Y.K. Ng, S.T. Tan, Study of EDL effect on 3‐D developing flow in microchannel with Poisson– Boltzmann and Nernst–Planck models, International Journal for Numerical Methods in Engineering, 71 (2006) 818-836
2006
Reviewed August 10, 2026 · model on record in the stance chip above.
Discussion (0). Continue with ORCID to comment.