Pith. sign in

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 →

arxiv 2501.05917 v1 pith:FGDP7YE6 submitted 2025-01-10 physics.chem-ph

classification physics.chem-ph
keywords Poisson-Nernst-Planckequationsfinitedifferencemethodconcentrationprofilestabilityanalysis3-DsimulationButler-Volmerboundaryconditionelectrochemicaltransportexplicittime-stepping
verification ladder T0 review T1 audit T2 compute T3 formal

The pith

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The reading

The paper tries to put 3-D Poisson-Nernst-Planck (PNP) simulation on firmer footing by giving a step-by-step finite-difference recipe and, crucially, the stability conditions that an explicit scheme must satisfy. It derives two inequalities, Eqs. (30) and (36), from Fourier analysis of the diffusion, migration, and Poisson updates, and argues that combining them makes explicit coupling of the full PNP system nearly impossible for realistic batteries or concentrated electrolytes. To get usable results, it treats the Poisson equation as an elliptic boundary-value problem solved by Gauss-Seidel iteration and applies Butler-Volmer kinetics at the electrode surface. If the stability analysis is right, modelers can check their time step and grid choice against closed-form limits before simulating, and will know why many existing 3-D PNP attempts diverge.

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.

Watch

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

Editorial extensions of the paper, not claims the author makes directly.

  • 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.
Share X Bluesky LinkedIn Reddit HN

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, and a circularity audit.

Referee Report

4 major / 5 minor

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)
  1. [§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. [§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.
  3. [§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.
  4. [§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)
  1. [§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.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.
  3. [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.
  4. [§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.
  5. [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

0 steps flagged · score 0.0 of 10

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 4 free parameters · 6 assumptions · 0 invented entities

The central claim rests on a standard explicit discretization, three domain assumptions about the electrolyte and boundaries, and two ad hoc mathematical premises: treating the potential as constant during migration stability analysis and adding amplification factors for the coupled system. No new physical entities are introduced and no parameters are fitted to experimental data; the free parameters are hand-selected simulation inputs.

free parameters (4)
  • D_z = 1e-2 m^2/s
    Hand-chosen diffusion coefficient in the z direction, 1000 times larger than D_x and D_y; strongly affects the stability limit and the concentration profiles.
  • J_0 = 1.0 mA/cm^2
    Hand-chosen exchange current density in the Butler-Volmer boundary condition; no experimental source or sensitivity study.
  • beta = 0.5
    Transfer coefficient in the Butler-Volmer equation; chosen within the typical 0.3 to 0.7 range, not measured.
  • N_t = 50 to 80 (recommended)
    Number of time steps recommended for stability; based on observed instability at N_t=100 and N_t=200 rather than on a derived criterion.
assumptions (6)
  • standard math Central finite difference approximations for first and second derivatives.
    Used throughout Sections 2.1.3 and 2.1.5; standard but supplies the discretization on which all results rest.
  • domain assumption Electrolyte is homogeneous with a uniform initial concentration of 1 M.
    Section 2.1.1 assumption (1) and Eq. (8); simplifies the model to a uniform starting state.
  • domain assumption Ion-ion interactions beyond the mean-field Poisson coupling are neglected.
    Section 2.1.1 assumption (5) and Section 2.2.1; limits validity to dilute electrolytes and excludes crowding effects.
  • domain assumption No ion flux through the lateral and top boundaries.
    Section 2.1.1 and Eq. (44); isolates the system but conflicts with the zero-concentration condition in Eq. (45).
  • ad hoc to paper The electric potential is treated as a constant coefficient in the Fourier stability analysis of the migration term.
    Eqs. (20)-(26); the analysis does not account for the feedback between concentration and potential.
  • ad hoc to paper The total amplification factor is the sum of the N-P and Poisson amplification factors.
    Eq. (43); no derivation is given for this sum, and it is not the standard way to analyze coupled systems.

how reviews work

0 comments
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

Figures reproduced from arXiv: 2501.05917 by the authors.

Figure 1
Figure 1. A general process of solving PDEs using the finite difference method. [PITH_FULL_IMAGE:figures/full_fig_p002_1.png] view at source ↗

Discussion (0). Continue with ORCID to comment.

Reference graph

Works this paper leans on

10 extracted references · 10 canonical work pages

  1. [1]

    Yan, M.Z

    D. Yan, M.Z. Bazant, P.M. Biesheuvel, M.C. Pugh, F.P. Dawson, Theory of linear sweep voltammetry with diffuse charge: Unsupported electrolytes, thin films, and leaky membranes, Phys Rev E, 95 (2017) 033303

  2. [2]

    Horng, T.C

    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

  3. [3]

    Mitscha -Baude, A

    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

  4. [4]

    Newman, C.W

    J.S. Newman, C.W. Tobias, Theoretical Analysis of Current Distribution in Porous Electrodes, Journal of The Electrochemical Society, 109 (1962) 1183

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

  6. [6]

    Souza, M.T

    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

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

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

Show all 10 references
  1. [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

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

Pith tools

Reviewed August 10, 2026 · model on record in the stance chip above.