Geometrically Induced Acceleration for Charging Dynamics of Electrical Double-Layers in a Nanopore with Sloped Walls

TL;DR

This paper shows how sloped walls in nanopores can speed up charging while improving energy storage in supercapacitors.

Contribution

The study introduces a new theoretical framework linking pore geometry to double-layer dynamics, enabling faster and more efficient supercapacitor designs.

Findings

Sloped pore walls induce additional ionic flux, accelerating charging and enhancing charge storage.

Theoretical predictions align closely with simulations while reducing computational cost by 5–6 orders of magnitude.

A modified circuit model captures geometric variations and can be integrated into pore-network models.

Abstract

Confinement strongly
influences electrochemical systems, where
structural control has enabled advances in nanofluidics, sensing,
and energy storage. In electric double-layer capacitors (EDLCs), or
supercapacitors, energy density is governed by the accessible surface
area of porous electrodes. Continuum models, built on first-principles
transport equations, have provided critical insight into electrolyte
dynamics under confinement but have largely focused on pores with
straight walls. In such geometries, a fundamental trade-off emerges:
wider pores charge faster but store less energy, while narrower pores
store more charge but charge slowly. Here, we apply perturbation analysis
to the Poisson–Nernst–Planck (PNP) equations for a single
pore of gradually varying radius, focusing on the small potential
and slender aspect ratio regime. Our analysis reveals that sloped
pore walls induce an additional ionic flux, enabling simultaneous
acceleration of charging and enhancement of charge storage. The theoretical
predictions closely agree with direct numerical simulations while
reducing computational cost by 5–6 orders of magnitude. We
further propose a modified effective circuit representation that captures
geometric variation along the pore and demonstrate how the framework
can be integrated into pore-network models. This work establishes
a scalable approach to link pore geometry with double-layer dynamics
and offers new design principles for optimizing supercapacitor performance.

Figures

A schematic of generalized nanopore, with red circles representing cations and blue circles representing anions. Although initially uncharged, at time τ = 0, a negative potential of Φ w is applied to the pore wall, driving cations into the pore and expelling anions. There exist three distinct regions: the reservoir, the static diffusion layer, and the pore. The reservoir acts as an unchanging repository of charged ions, allowing ions to flow in and out freely. The reservoir remains electroneutral, serving as the reference level for the potential. The static diffusion layer (SDL) serves as a bridge between the reservoir and the pore. While electrically neutral, its potential increases linearly as ions approach the pore entrance. The radius and length scale of this region are adjustable parameters, when formulated appropriately, define the dimensionless Biot numbera quantity characterizing the resistance of electrochemical flux from the reservoir region to the pore entrance. The pore region is where the EDL forms, and is the main focus of analysis. Importantly, the walls are ideally blocking, and the pore’s axial length is much larger than its radial dimension. The radius of the pore depends on the axial coordinate, given by α­(Z).

A schematic of an equivalent circuit diagram for the dimensional electrochemical potential of charge, μ w *, inside a converging pore. The reservoir region acts as the reference potential, with the SDL region bridging the reservoir and the pore regions. This SDL region is treated as a lumped resistor (R SDL), which is inversely related to the Biot number. Inside the pore, each axial point has an associated resistance (R p ) and capacitance (C p ), which are functions of radius, Debye length, diffusivity, and electric permittivity. Each axial junction has current i flowing through it, and the current di is being stored in the EDL. Driving this electrokinetic flow is an applied wall potential of μ w * = 2Φ w *.

Electrochemical potential μ̂(Z,τ), the electric potential Φ̂(Z,R,τ), and the charge density ρ̂(Z,R,τ) plotted at various timesearly (a,c,e) and intermediate times (b,d,f)against direct numerical simulations. (a,b) Represent plots of the electrochemical potential as a function of axial position, and include the SDL region when Z < 0. The radial electric potential (c,d) and charge density (e,f) profiles are plotted midway, (Z = 0.5), through the pore as a function of radial position. Schematics of all the pores are depicted on the right side of the figure: diverging, narrow, wide, and converging. The entrance-to-end ratio for the converging and diverging cases are 2:1 and 1:2, respectively, with corresponding equations α­(Z) = 2 – Z and α­(Z) = 1 + Z, respectively. The wide α­(Z) = 2 and narrow α­(Z) = 1 pores represent the widest and narrowest radii for the changing radius cases. τ denotes the nondimensional time, Z the nondimensional length along the pores, and R the nondimensional radial length.

Electrochemical potential μ̂(Z,τ) and the electric potential Φ̂(Z,R,τ) plotted at various timesearly (a,c) and intermediate times (b,d)against direct numerical simulations with larger changes in slope. (a,b) Represent plots of the electrochemical potential as a function of axial position, and include the SDL region when Z < 0. The radial electric potential (c,d) profiles are plotted midway, (Z = 0.5), as a function of radial position. Three converging pores and one diverging pore are tested. The converging pores share the same entrance radius of 20 nm, but differing in their end radii. The 20–10 nm case corresponds to the original converging pore discussed previously, while the 20–1 nm and 20–0.1 nm pores explore steeper slopes, testing the limits of the perturbation analysis. Finally, the diverging pore is also added with a radius starting of 2 nm and an ending radius of 20 nm.

Contour plots of charge density (ρ̂) throughout the pore at various times. (a–c) Depict the converging geometry, and (d–f) depict the diverging geometry, described by α­(Z) = 2 – Z and α­(Z) = 1 + Z, respectively. The charge density contours are shown at early (a,d), intermediate (b,e), and equilibrium (c,f) times. Lighter regions correspond to higher charge densities, while darker regions represent a more neutral state. τ denotes the nondimensional time, Z the nondimensional length along the pores, and R the nondimensional radius. The initial and system conditions are μ̂(Z,τ = 0) = 1, κ = 2, l s * = l p * and a s * = 4a p *.

Quiver plot of ∇ρ̂(Z,R,τ) and ∇Φ̂(Z,R,τ), representing the diffusive (light pink) and electromigrative (dark blue) fluxes, respectively. The flux profiles within the pores at early (a,d), intermediate (b,e), and equilibrium (c,f) times. The entrance-to-end ratio set at 2:1 for the converging pore and 1:2 for the diverging pore, defined by α­(Z) = 2 – Z and α­(Z) = 1 + Z, respectively. τ denotes the nondimensional time, Z the nondimensional axial coordinate, and R the nondimensional radial coordinate. The initial and system conditions are μ̂(Z,τ = 0) = 1, κ = 2, l s * = l p * and a s * = 4a p *.

Plot of charge fractionquantified as total instantaneous charge (Q) over total steady state charge (Q SS)as a function of time (a), where μ̂(Z,τ = 0) = 1, κ = 2, and a s * = 4a p *. The dotted line represents the τ1/2 or the time it takes a given pore to reach half of its equilibrium charge. For comparison, in plots (b,c), each τ1/2 is normalized by the half charge time of the narrow pore τ n . This ratio plotted on logarithmic scale versus κ for two different characteristic length scales of the SDL. These plots demonstrate the effect of entrance resistance; specifically, plot (b) exhibits a higher entrance resistance than plot (c).

Full text

Introduction

Confinement is common in electrochemical systems, and careful control of structure in these systems has enabled breakthroughs in nanofluidics, single-molecule sensing, and enhancements in energy storage devices. ?−? ? Focusing on energy storage devices, electric double-layer capacitors (EDLCs), or supercapacitors, offer high power density along with minimal cyclical degradation. ?,? EDLCs have diverse current and potential applicationsincluding flexible self-powered electronics, balancing consumption/generation within the electric grid, and powering electric vehicles ?−? ? ? ? but are limited by their energy densities. Advancements in room temperature ionic liquids (RTILs), and their corresponding theory, have yielded improved understanding with increased performance and safety with their extended voltage window, thermal stability, and nonflammable nature. ?−? ? ? ? Since the performance of EDLCs is proportional to surface area, porous carbon electrodes are employed to enhance energy storage.? This raises a key question: how do we optimize pore design inside electrodes? To this end, an understanding of electrolyte transport under charged confinement is crucial to obtain structur–property relationships and optimize performance.

Molecular dynamics (MD) simulations are commonly used to model electrolyte–electrode interactions under confinementoffering insights into atomic-scale phenomenasince pore length scales are often comparable to the ionic radius. Recently, MD studies have revealed insights into the dynamics of ionic liquids, ?,? optimal charging–discharging cycles, ?,? the importance of gradual desolvation in nanopores,? and benefits of ionophobic pore walls. ?,? Although MD simulations are able to capture atomic-level resolution inside nanopores, their computational expense typically limits their simulation size to a single nanopore and short time scales. Hence, parallel advances are required in continuum modeling to predict electrolyte transport in complex geometries and at long time scales.

At the continuum level, de Levie pioneered modeling EDL formation within porous electrodes with an equivalent circuit model, ?,? which remains influential today. In the past few years, there have been exciting developments in this area, such as the inclusion of nonlinear potentials,? Faradaic reactions,? surface conduction,? arbitrary Debye lengths,? diffusivity asymmetry, ?,? floating electrodes concentrated electrolytes,? impact of convection? and a network of pores. ?−? ? Notably, some of these approaches model a macrohomogeneous electrode, representing the electrode collective behavior rather than at the pore scale. These models can produce equivalent circuits similar to de Levie transmission line models; however, the underlying assumptions and starting points of these approaches are different. While these models provide deep insights into the charging of double-layer inside porous materials, many studies focus on cylindrical or slit-pore ?−? ? shapes with a constant radius or height, i.e., straight walls. One key result that these studies reveal is that for small to moderate potentials, wider pores charge faster at the expense of a lower charge density.? Here, we address the question: is it possible to use sloped walls to enhance the charging rates, while also increasing the density of charge stored? Our results reveal, rather surprisingly, that this is indeed possible due to the additional flux that the sloped walls create.

By building on our prior work, ?−? ?,? we employ perturbation analysis to develop a theoretical model to predict double-layer charging for an axisymmetric pore with gradually varying radius, under the long-slender pore assumption. The theoretical model reveals that sloped walls create an additional flux for all Debye lengths, which can be used to speed up (or slow down) the charging of double-layers. The model is in excellent agreement with direct numerical simulations and can successfully capture the spatiotemporal variations of charges and potential with a 5–6 order-of-magnitude reduction in computational cost. We also proposed a modified effective circuit that is able to capture geometrical variations along the pore. Finally, given that our theory is computationally friendly, it can easily be integrated with network models to simulate thousands of interconnected pores ?−? ? and provides a crucial step toward obtaining structure–property relationships in supercapacitors to optimize performance.

Methods

Numerical Solution for the Linearized Model

The theoretical framework developed provides the governing equations, boundary, and initial conditions necessary to solve the transient pore charging problem. Using the method of lines, the time-dependent partial differential equations are transformed into a system of ordinary differential equations (ODEs). A finite difference scheme is employed to discretize the spatial domain, and the resulting ODE system is then solved over the specified time domain. SciPy’s built in function odeint was used to solve this ODE system, which uses the Livermore Solver for Ordinary Differential Equations (LSODA). The resulting solution describes the time evolution of the electrochemical potential of charge throughout the pore. From the electrochemical potential, the electric potential and charge density can be found. Reworking eqs and ?, the governing expressions for these quantities as functions of radius can be derived. The final expressions for the average charge density, average potential, charge density, and potential are as follows

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$

Direct Numerical Simulation Details

Open-source software OpenFOAM was used to perform direct numerical simulations (DNS) and solve eq. A pore of length 10 μm and a cylindrical static diffusion layer (SDL) region of length 5 μm, were connected via a transition region with an arc to smooth out the geometrical transition and reduce the numerical errors. Four pore geometries were compared against DNS: converging, diverging, narrow and wide. For the narrow and wide cases, the pore radius was held constant at 10 and 20 nm, respectively. The converging pore starts from the radius of the wide pore and linearly decreases until it reaches the radius of the narrow nanopore. Conversely, the diverging pore starts with the radius of the narrow pore and linearly increases to the radius of the wider pore. The SDL radii were fixed to be 40 nm except for the narrow pore cases, in which case the radius was set to 20 nm. The boundary conditions were kept identical to our prior work of Gupta et al.? The vacuum permittivity is 8.85 × 10^–12^ F/m, electrical permittivity of water is ϵ = 7.1 × 10^–10^ F/m (assuming relative permittivity to be 80.2), diffusion coefficient D = 1.34 × 10^–9^ m^2^/s, and reservoir concentration c 0 = 0.94 mM were used. The simulations were performed using a 28-core workstation with a computational cost of 4 million seconds for each second elapsed in the simulation.

Results and Discussion

Problem Setup

We consider an idealized axisymmetric single pore where the pore radius changes along the axial direction; see Figure. We assume the electrolyte to be binary, monovalent with valences z ± = ± 1 and symmetric diffusivities such that D ± = D. The concentration of the cation is described by c +; conversely, the concentration of the anion is described by c –. We ignore surface redox reactions and assume the surface to be ideally blocking, i.e., ion fluxes through a surface are negligible.

A schematic of generalized nanopore, with red circles representing cations and blue circles representing anions. Although initially uncharged, at time τ = 0, a negative potential of Φ w is applied to the pore wall, driving cations into the pore and expelling anions. There exist three distinct regions: the reservoir, the static diffusion layer, and the pore. The reservoir acts as an unchanging repository of charged ions, allowing ions to flow in and out freely. The reservoir remains electroneutral, serving as the reference level for the potential. The static diffusion layer (SDL) serves as a bridge between the reservoir and the pore. While electrically neutral, its potential increases linearly as ions approach the pore entrance. The radius and length scale of this region are adjustable parameters, when formulated appropriately, define the dimensionless Biot numbera quantity characterizing the resistance of electrochemical flux from the reservoir region to the pore entrance. The pore region is where the EDL forms, and is the main focus of analysis. Importantly, the walls are ideally blocking, and the pore’s axial length is much larger than its radial dimension. The radius of the pore depends on the axial coordinate, given by α­(Z).

We assume there exist three distinct regions, see Figure: the reservoir, the static diffusion layer, and the pore. ?,?,?,?,?−? ? Far from the electrode surface, the reservoir consists of ions with concentration c ± = c 0 and potential ϕ = 0. This reference potentialchosen for mathematical convenienceis ultimately arbitrary and different system perspectives exist; for instance, a 1D two-electrodes system would have a point reference potential.? An electroneutral static diffusion layer (SDL) forms between the pore entrance and the reservoir. ?,?,?,?,? Although the pore regionwhere ions accumulate to form EDLsis the main focus of this analysis, the SDL is meant to capture transport limitations (entrance resistance) outside the pore. The SDL and pore region have associated length scales l s ^^ and l p ^^ respectively, as well as radii a s ^^ and a p ^^noting that () terms are dimensional quantities. We assume the pore length scale to be much greater than the radial length scale, such that a p ^^≪l p ^*^. We write the shape equation of the pore as

$$[eqn]$$

where α­(z) is the local radius, which changes along the z–direction. We assume that the length scales of the SDL and pore region are of the same order, hence l s ^^/l p ^*^ = O(1).

Initially, the concentration everywhere inside SDL and the pore is assumed to be the same as the concentration of the reservoir, c ± = c 0. To initiate charging, we set the pore surfaces to be a constant potential ϕ = Φ_ w _ ^*^ at t > 0, where t is time. To illustrate the physical process, we assume the potential to be negative. Consequently, cations will preferentially transport into the pore and anions will transport out of the pore, forming a double-layer. The double-layer will form initially at the pore mouth only. As time progresses, ions will continue to flux in and out of the pore, developing the double-layer further into the pore until equilibrium is reached. Our objective is to quantify this dynamic process and especially elucidate the effect of pore geometry.

Derivation

Outline

Due to the length of the derivation section, we provide a brief summary here and highlight the key assumptions of our work. Given the number of symbols introduced, a summary of the nomenclature is included below the Acknowledgments. Our goal is to capture the transient charging behavior within a single shape-changing pore. We use the Planck–Nernst–Poisson (PNP) equations to model ion fluxes within the pore. We apply a two-term perturbation expansion in the small potential and slender aspect ratio regimes to linearize the PNP equations. From this expansion, we derive the electrochemical potential of charge and find that it is radially constant. This quantity, which we will later define, allows us to recover both the potential and charge density profiles within the pore. After radial averaging, we recover a transient diffusion equation that depends only on the axial position. For the boundary conditions used to solve our PDE, we assume the walls to be ideally blocking and the pore to be connected to an electrically neutral regionthe static diffusion layerwhere a flux-matching condition links the two domains. This electroneutral region introduces additional parameters that account for transport inefficiency from the reservoir to pore region.

Governing

Equations

We first focus on the pore region and capture the SDL afterward. We assume that the ion transport inside the pore is described by the Nernst–Planck equation, which is given as

$$[eqn]$$

where the ion fluxes N ± are given as

$$[eqn]$$

where ϕ is the electrical potential, e is the fundamental charge, k _ B _ is the Boltzmann constant and T is the absolute temperature. The first term of the RHS is the diffusive flux and the second term is the electromigrative flux. We note that several other modifications of the Poisson–Nernst–Planck exist, including finite ion-size, ?,?,?,? dielectric decrement ?,? and ion–ion correlations, ?−? ? ? ? but are not considered for simplicity. While these effects will modify our results, we anticipate that the qualitative influence of geometry is model-independent; however, the quantitative impacts may vary, and a full analysis of the modified equations will be required to determine its magnitude. We also note that recently, convection has been argued to be important in the pore charging process,? but it is not included here. This is consistent with the linearized regime we investigate, as explained later. However, interested readers are referred to for the effect of convection in the nonlinear regime.?

To complete the model, we invoke the Poisson equation

$$[eqn]$$

where ϵ is the electrical permittivity of the electrolyte solution, which is treated as constant in this work, but ϵ can depend on ion concentration due to dielectric decrement effects. ?,?

Nondimensional Equations

For simplification, we nondimensionalize by the cation and anion reservoir concentration, c 0, into a charge density, $[eqn]$ , and a salt concentration, $[eqn]$ . The potential is normalized by the latent thermal voltage, $[eqn]$ . The time is nondimensionalized by the diffusion time $[eqn]$ . The Debye length is given as $[eqn]$ and nondimensional Debye ratio $[eqn]$ , where a p ^*^ is a reference pore radius. The two different length scales (radial and axial) are used to nondimensionalize the length, so $[eqn]$ and $[eqn]$ . We also define $[eqn]$ . Converting eq in terms of these nondimensional parameters

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

where $[eqn]$ and () denotes a dimensional length. We introduce the reference radius, a _ p _ ^^, to enable relative comparison between pores of different sizes. For example, if pore A is twice as large as pore B, then pore A’s nondimensional wall equation would correspond to 2α­(Z) given the same a _ p _ ^^. Although other normalization schemes can be used, we keep a _ p _ ^^ constant between all pores and simulations for simplicity.

Perturbation Expansion

To make progress in solving eq, we consider a two-term perturbation expansion for the parameters Φ_ w _ and δ, where both Φ_ w _ ≪ 1 and δ ≪ 1. The first parameter, Φ_ w _ corresponds to a small applied wall potential, and our prior work has shown that a linearized solution compares favorably to direct numerical simulation up to Φ_ w _ ≈ 4.? The second parameter, δ, represents the slender aspect ratio, defined as δ = a _ p _ ^*2^/l _ p _ ^*2^≪1. Beyond restricting the pore’s aspect ratio, δ imposes further constraints on pore wall variations, which are quantified later on. The regular perturbation expansions are then proposed as follows

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

here, the first subscript denotes the order of the solution in Φ_ w , and the second subscript denotes the order in δ, e.g. ρ_21 corresponds to O­(Φ_ w _ ^2^δ^1^). Although not immediately evident, except for S 00, all O­(Φ_ w _ ^0^δ^ n ^) terms vanish. In the absence of an applied electric field, there is no driving force for ions to move in or out of the pore, and thus O­(Φ_ w _ ^0^δ^ n ^) corrections are not included in the expansion. The inclusion of δ is vital, as it enables a lubrication-like approximation that accounts for variations in pore radius.

The leading-order problem is trivial to solve, as it corresponds to a binary electrolyte with no applied potential. Consequently, no charge density or potential gradient develops within the electrolyte, yielding ρ_00_(R, Z, τ) = Φ_00_(R, Z, τ) = 0. The ion concentration remains constant and equal to that of the reservoir concentration; hence, the leading-order salt concentration S 00(R, Z, τ) = 2.

A careful examination of the next order reveals that radial scaling plays a crucial role in determining which terms appear. In particular, due to geometric scaling in front of the gradient operator, eq has a term at O­(Φ_ w _δ^–1^), which reads

$$[eqn]$$

after integrating and employing symmetry at R = 0 becomes

$$[eqn]$$

Next, we find at O­(Φ_ w _ ^1^δ^0^), eq and ? become

$$[eqn]$$ $$[eqn]$$

We note that eq is neglected, as it is not required to determine the potential and charge density profiles. We also note that the additional factor of 2 for Φ_11_ and Φ_10_ is a result of the leading order salt concentration, S 00 = 2. As eqs and ? provide the lowest-order nontrivial solution, i.e., O­(Φ_ w _ ^1^δ^0^), we limit the analysis to them, as they capture the essential features of charging. However, we would like to eliminate the ρ_11_ + 2Φ_11_ from eq.

To do so, we use the no-flux condition at the wall, n·(**∇**ρ + S **∇**Φ) = 0. Mathematically from eq, n∝∇(r-α^*^(z)), which in dimensionless variables yields

$$[eqn]$$

Substituting eq in n·(**∇**ρ + S **∇**Φ) = 0 at R = α­(Z) and expanding to O­(Φ_ w _ ^1^δ^0^), we obtain

$$[eqn]$$

Equation imposes an additional constraint on the pore geometry, that being α­(Z) must be smooth and $[eqn]$ (note that in the nonstretched coordinates, this would be O­(δ)). Sharp variations in wall geometry would violate the expansion at this order. Physically, this requirement can be understood by noting that steep geometric gradients would introduce local flux restrictions within the pore, breaking the assumption of slow axial variation along the pore. Equation can be used to eliminate ρ_11_ + 2Φ_11_ in eq. However, since eq is only valid at R = α­(Z), we radially integrate eq to write

$$[eqn]$$

where the first term on the right-hand side can integrated to write

$$[eqn]$$

Finally, substituting the boundary condition from eq in the above, we obtain the following equation only in O­(Φ_ w _ ^1^δ^0^) variables

$$[eqn]$$

Due to eq, $[eqn]$ is constant for all R values, and thus does not need to be specifically evaluated R = α­(Z).

Expressing Equations in

Terms of the Electrochemical Potential

Hereafter, we omit the subscript notation for O­(Φ_ w _ ^1^δ^0^) terms. For clarity, we define $[eqn]$ and $[eqn]$ . We also introduce $[eqn]$ , referred to as the electrochemical potential of charge; briefly restoring dimensions, we observe, $[eqn]$ . This quantity corresponds to the linearized electrochemical potential difference between the cation and the anion, $[eqn]$ ,? and this combined quantity will be shown to be a useful metric for describing the charging dynamics along the pore. Our goal is to express the Nernst–Planck equation at O­(Φ_ w _ ^1^δ^0^) in terms of μ̂.

We start by noting that eq can be transformed into μ̂ to be read as

$$[eqn]$$

or equivalently

$$[eqn]$$

as μ̂ is constant with respect to a given radial cross-section. Equation can also be expressed in terms of the average charge density $[eqn]$ and potential $[eqn]$ within a given radial cross-section

$$[eqn]$$

While these quantities are still unknown, ρ̅ and Φ̅ only depend on the axial position, as the averages are not functions of the radial coordinate.

First, we employ eq and integrate the Poisson equation in eq to eliminate Φ_10_ and find the relationship between ρ̅ and ρ̂. Second, we radially integrate eq to obtain a self-contained equation in μ̂(Z,τ) only.

Substituting $[eqn]$ from eq in eq and restoring the variables ρ̂, Φ̂ and μ̂, we get

$$[eqn]$$

which has solutions in terms of modified Bessel functions of order zero

$$[eqn]$$

We note that eq is implicitly time dependent, since the equation is linear in ρ̂ and a separation of variables holds. The time dependency will be introduced by the gradual screening of the double layers along the pore’s axis. In order for the solution to remain physical at R = 0, C 2 = 0. Next, C 1 can be written in terms of the average charge density $[eqn]$ by integrating eq from the center to the pore wall

$$[eqn]$$

or

$$[eqn]$$

where ρ̅(Z,τ) is still unknown. Calculating ρ̂ from eq at R = α and utilizing $[eqn]$ at R = α, eq yields

$$[eqn]$$

Next, we substitute $[eqn]$ and $[eqn]$ in eq to write

$$[eqn]$$

Using eq, μ̂ is factored out of the integral on the right-hand side of the above equation. Finally, integrating and using $[eqn]$ on the left-hand side yields

$$[eqn]$$

or equivalently

$$[eqn]$$

Or rewriting ρ̅ in terms of μ̂ using eq, we obtain the crucial result of this article, i.e.

$$[eqn]$$

This result showcases a critical feature of the variable pore radius, namely, the transport of the electrochemical potentials of charge by cross-sectional area variations. The term $[eqn]$ is akin to a diffusion with a variable diffusion coefficient. Upon expanding, eq can be rewritten in the useful form

$$[eqn]$$

where $[eqn]$ . The impact of change in cross-sectional area is captured by $[eqn]$ , which conveys the axial component of wall-normal ion fluxes in the formation of EDLs on a slanted surface and can be interpreted as a “pseudo-advective” term. When $[eqn]$ (a diverging profile), it corresponds to a negative velocity and transports counterions away from the pore end, hindering charging. In contrast, when $[eqn]$ , this mechanism accelerates charging. A comparison of the mechanisms of axial diffusion/migration and a tangential component of EDL formation can be estimated by the magnitude of the “pseudo-velocity coefficient” $[eqn]$ .

Boundary

Conditions

In order to solve eq, two boundary conditionsone at the pore entrance and one at the pore endalong with an initial condition are required. The initial condition is straightforward, i.e., at τ = 0 there is only an applied potential from the wall and no charge, $[eqn]$ . The pore end has a no-flux condition, or $[eqn]$ . The pore entrance condition requires a more careful analysis since it is connected to reservoir via a SDL.

The inclusion of the SDL serves to capture transport limitations through an entrance boundary condition, while also matching the linear potential profiles from DNS. The region outside the pore is modeled as an electroneutral cylinder, in which a flux matching condition is imposed between the SDL and the pore. The geometric length scales of this cylinder are treated as adjustable physical parameters. Although the true ion flux in the external electrolyte is more complex than this representation, this SDL model provides an approximation of the entrance resistance in the linear regime. Inside the SDL, since there are no charged surfaces in the radial direction, the governing equations only vary in the axial direction. Additionally, it is common in continuum simulations to assume that regions outside the pore behave linearly. ?,? Furthermore, due to small potentials, we can argue that S 00 = 2. Therefore, to maintain a constant flux of charge, the electric field must be uniform or equivalently the potential profile linear. As the electric and electrochemical potentials can be related by eq, the governing equation in this region is

$$[eqn]$$

In our prior work,? we show how the continuity of electrochemical potential of charge is required across an interface. This criterion can be understood from the gradient of the electrochemical potential of charge providing an estimate of charge flux. Across an interface, if there is a difference between electrochemical potentials, there would be a spurious flux of charge, which thus necessitates the electrochemical potential of charge to be continuous. Mathematically, $[eqn]$ at the pore-SDL interface and $[eqn]$ at the SDL-reservoir interface. Applying these boundary conditions, we write

$$[eqn]$$

Conserving the charge flux across the SDL-pore interface, we write

$$[eqn]$$

The fluxes are governed by the Nernst–Planck equations and can be expressed in terms of the electrochemical potential. Notably, the scaling for the SDL is different than the pore region, requiring a description of the left flux, (N left), to be expressed in terms of the ratios of length scales

$$[eqn]$$

where Bi $[eqn]$ is the Biot number. Examining the limiting behavior of the Biot number offers insight into its physical significance. As Bi → ∞, the SDL becomes vanishingly thin, effectively placing the reservoir in direct contact with the pore entrance. In this limit, the SDL imposes no resistance to ion transport from the reservoir into the pore. Conversely, as Bi → 0, the SDL becomes infinitely thick, acting as a region of high resistance. In this case, only a negligible flux can pass through the SDL to maintain electroneutrality, effectively isolating the pore from the reservoir and resulting in significantly slower charging dynamics. With the governing equation, boundary conditions, and initial condition defined, the model is complete and can be numerically simulated. Our model’s reduced-order enables rapid exploration of different geometries, offering a significant computational reduction over traditional direct numerical simulations.

Equivalent Circuit Model

Based on eq, we seek to propose an equivalent circuit that incorporates the charging dynamics for an arbitrarily shaped pore; see Figure. This circuit is designed to capture the charging behavior of the electrochemical potential in agreement with the governing equations, which are redimensionalized and rearranged into per-unit-length resistances (R _ p _) and capacitance (C _ p _). Our goal is to develop a circuit analogue that aids in the understanding of how pore radius influences charging and how the system behaves in the overlapping and thin EDL limits. Unlike typical circuit models that focus on electrical potential, this formulation centers on the electrochemical potential of charge, μ̂, which is fundamentally distinct. We emphasize that the combined quantity $[eqn]$ , and its corresponding gradient, captures the two fundamental driving mechanisms of the system: electromigration and diffusion. By basing our formulation on μ̂, the equivalent circuit remains valid across a broader range of pore sizes and double-layer thicknesses, providing a natural representation of charging behavior from our governing equations.

*A schematic of an equivalent circuit diagram for the dimensional electrochemical potential of charge, μ w

*, inside a converging pore. The reservoir region acts as the reference potential, with the SDL region bridging the reservoir and the pore regions. This SDL region is treated as a lumped resistor (R SDL), which is inversely related to the Biot number. Inside the pore, each axial point has an associated resistance (R

p ) and capacitance (C

p ), which are functions of radius, Debye length, diffusivity, and electric permittivity. Each axial junction has current i flowing through it, and the current di is being stored in the EDL. Driving this electrokinetic flow is an applied wall potential of μ w

• = 2Φ w

.

In the limit of thin EDLs in a pore, de Levie proposed an RC transmission line model ?,? that remains widely used today.? This model yields a Warburg element in an impedance plot, which is commonly argued to be a signature of electrolyte diffusion in a porous material.? Several modifications to the transmission line model have been proposed, which include nonlinear effects,? Faradaic reactions,? and surface conduction.? Additionally, these corrections have focused primarily on the thin EDL limit.? More recently, studies have started to focus on arbitrary Debye lengths ?,?,?,? and asymmetric diffusivities. ?,? In general, the analyses in the literature focus on straight cylindrical or slit-pores and do not consider a change in pore radius or height. To the best of our knowledge, the only works that consider pore shape and its effect on the equivalent circuit are the works of Black and Andreas,? Keiser et al.? and Cooper et al.? Black and Andreas? employ an ad-hoc approach to the equivalent circuit and are not grounded in first principles. While our results bear significant resemblance to those of Keiser et al.? and Cooper et al.,? this agreement holds only in the thin EDL limit. In this regime, the capacitance per unit length scales with the pore’s surface area, and the resistance varies inversely with cross-sectional area, consistent with Keiser et al. However, unlike previous models, our approach extends to a broad range of EDL thicknesses, capturing the charging dynamics of both moderate and overlapping double-layers; thus, representing a significant advancement beyond earlier work.

Starting from eq, we redimensionalize our governing equation

$$[eqn]$$

where $[eqn]$ is the dimensional electrochemical potential of charge, λ is the Debye length, α*­(z) is the dimensional radius function, and I _ n _ are modified Bessel functions of the first kind of order n. As in de Levie’s analysis, we find it convenient to define capacitance per unit length as C _ p _ and the resistance per unit length as R _ p . Consider a node at an arbitrary axial position z in the circuit of Figure. We apply a current balance to this differential node, yielding di _ z _ = i| z _ – i|_ z+dz _. Using Ohm’s law, the incoming and outgoing currents into this node can be expressed as

$$[eqn]$$ $$[eqn]$$

Any current accumulated in this region must be stored as a capacitive current in the EDL, at location z. From the definition of capacitance

$$[eqn]$$

where $[eqn]$ is the applied dimensional wall potential.

We now seek appropriate expressions for R _ p _ and C _ p _, such that the resulting equivalent circuit formulation is consistent with the governing equation and limiting behavior; hence, before providing the mathematical definitions, we discuss the physical intuition of the result. Since the total ionic concentration remains constant, the conductivity of the electrolyte also stays constant. Therefore, R _ p _ should be inversely proportional to the local cross-sectional area. In addition, the local capacitance per unit length, C _ p _, should approach different values in the thin and overlapping double-layer limits. In the thin double-layer limit of α*/λ → ∞, the capacitance per unit surface area should approach the Debye–Huckel limit of $[eqn]$ , which implies that C _ p _ = πϵα*/λ (note that this is a factor of 1/2 different than the de Levie model because $[eqn]$ contains a factor of 2Φ* and hence capacitance has to be adjusted to maintain equal charge). For the limit of $[eqn]$ , physically, no potential is screened and the charge varies only axially. For $[eqn]$ , redimensionalizing eq, $[eqn]$ reveals $[eqn]$ . With these limits in mind, we argue that R _ p _ and C _ p _ are given as

$$[eqn]$$ $$[eqn]$$

Substituting these definitions for resistance and capacitance into the current balance yields

$$[eqn]$$

which is identical to eq. Finally, we note that resistance per unit length in the SDL, R _ p,SDL_, will simply be scaled by the area of the SDL, or $[eqn]$ . While we propose this interpretation of the equivalent circuit; alternative formulations are possible. The present representation should be viewed primarily as a conceptual tool to interpret the charging dynamics from a transmission-line perspective, rather than as a unique or definitive model.

Pore Descriptions

We focus on the transient charging behavior inside conical and cylindrical geometries, although our framework can be readily incorporated into other geometries. To describe the conical geometries, we employ

$$[eqn]$$

where m < 0 corresponds to a converging pore and m > 0 to a diverging pore. For the converging geometry, b = 2 and m = – 1, whereas for the diverging geometry, we used b = 1 and m = 1. The values of b and m are chosen such that the converging and diverging pores represent the same geometry but with the entrance and ends flipped. Since there are no reactions and μ̂→0 at equilibrium, the charge density and potential profiles of both the converging and diverging pores are identical, but mirrored. For the wide pore, we choose b = 2 and m = 0, and for the narrow geometry, b = 1 and m = 0. The nondimensional geometries remain fixed throughout the paper, and all discussion refers to these geometries, unless otherwise noted. To elucidate geometric effects on the transient charging behavior, we analyze charge density distributions, compare electromigrative and diffusive fluxes, and vary entrance resistance. For model validation, we compare our reduced-order framework to a fully resolved Planck–Nernst–Poisson simulation, demonstrating excellent agreement to DNS while providing major computational savings. It is important to note that κ, or the ratio of reference length scale (a _ p _ ^*^) with the Debye length (λ), is held constant for all the simulations at κ = 2, unless otherwise noted. While some quantitative values change, we show in the Supporting Information that the trends hold for κ = 0.1 and κ = 10 as well; see Figures S1, S2.

After solving eq for the four geometries, we obtain μ̂(Z,τ), which we substitute in eq to obtain ρ̂(R,Z,τ) and Φ̂(R,Z,τ). The calculations were performed for Φ_ w _ = 0.5 for all simulations, except the comparison to the direct numerical simulation (DNS), where at 10 mV, Φ_ w _ ≈ 0.389. At τ = 0, the pore is uncharged but set to be at a potential of Φ_ w _. Hence, $[eqn]$ for the current discussion. A value of $[eqn]$ corresponds to an uncharged region of the pore, whereas values approaching μ̂→0 indicate a charged region. At equilibrium, μ̂(Z,τ) = 0. Although the time is nondimensionalized by the pore diffusion time, τ = Dt/l _ p _ ^*2^, this scaling does not reflect the pore charging time scale as other factors affect charging times more: pore geometry, entrance resistance, and the Debye length ratio κ. For instance, in the thin double-layer limit, the formation of the EDL occurs on a much shorter time scale; these trends are illustrated in the κ = 10 case found in Supplementary Figure S2.

For a fair comparison, we kept the length scales of the SDL constant across all simulations, except the DNS comparison, and appropriately adjusted Biot number based on the definition of Bi. We assume that l _ s _ ^^ = l _ p _ ^^ and a _ s _ ^^ = 4a* _ p _ ^*^. As a result, the Biot numbers are set at Bi = 4 for the wide entrances and Bi = 16 for narrow ones. Although the SDL is modeled as a cylinder with two characteristic length scales, other constructions are possible. This SDL formulation provides a controlled way to represent entrance resistance, while still preserving the effect of entrance size on the charging dynamics. By fixing the length scales in the SDL region, the size of the transport environment outside the pore remains consistent across cases with different entrance radii. Other constructionssuch as those that hold the entrance resistance (or Biot number) constantdo not necessarily share this feature. For example, if the Biot number is fixed, the resulting entrance boundary conditions become identical for both wide and narrow entrances, making a fair comparison difficult.

Direct Numerical Simulation Comparison

To validate the reduced-order model (ROM), we compared its predictions directly against DNS of the full nonlinear PNP equations. We find close agreement between the DNS and ROM, demonstrating that the ROM captures the influence of geometry on the essential features of pore charging, like the potential and charge distribution profiles. Moreover, additional DNS comparisons, Figure, demonstrate that the ROM retains accuracy even for larger wall slopes. For simulation parameters, we set the lengths of the cylindrical SDL to be, l _ s _ ^^ = 0.5l* _ p _ ^^+3max­(α­(Z)^^/2), and the radii to be a _ s _ ^^ = 2max­(α­(Z)^^). Therefore, the resulting Biot numbers are Bi ≈ 31.81 for the diverging pore and Bi ≈ 7.95 for all other pores, in Figure. The choice of dimensionless group in the DNS was largely driven by numerical convenience. In dimensional terms, the simulated pores are 1 μm in length and have varying radii ranging from 20 nm down to 0.1 nm, with a characteristic Debye length fixed at 10 nm. As we show later, we can recover the spatiotemporal variations of charge density, potential, and the electrochemical potential, giving us confidence to investigate the trends from our reduced-order model systematically.

Electrochemical potential μ̂(Z,τ), the electric potential Φ̂(Z,R,τ), and the charge density ρ̂(Z,R,τ) plotted at various timesearly (a,c,e) and intermediate times (b,d,f)against direct numerical simulations. (a,b) Represent plots of the electrochemical potential as a function of axial position, and include the SDL region when Z < 0. The radial electric potential (c,d) and charge density (e,f) profiles are plotted midway, (Z = 0.5), through the pore as a function of radial position. Schematics of all the pores are depicted on the right side of the figure: diverging, narrow, wide, and converging. The entrance-to-end ratio for the converging and diverging cases are 2:1 and 1:2, respectively, with corresponding equations α­(Z) = 2 – Z and α­(Z) = 1 + Z, respectively. The wide α­(Z) = 2 and narrow α­(Z) = 1 pores represent the widest and narrowest radii for the changing radius cases. τ denotes the nondimensional time, Z the nondimensional length along the pores, and R the nondimensional radial length.

To compare our linearized reduced-order model against DNS, we plot μ̂(Z,τ), Φ̂(Z = 0.5,R,τ), and ρ̂(Z = 0.5,R,τ) for the four different geometries at two different times. First, we note that we obtain excellent agreement between DNS and our proposed model, both for axial variation as well as the radial variation. This demonstrates that our model can capture the physics of EDL charging in the small potential limit for shape-changing pores without any loss of accuracy.

Early time dynamics are represented by τ = 0.13; see Figurea,c,e). One would expect that the dynamics of the converging and diverging pores would lie between the dynamics of the narrow and wide pores. However, the diverging pore with α­(Z) = 1 + Z charges slower throughout the majority of the pore compared to all other cases, and the converging pore α­(Z) = 2 – Z charges faster than all other cases. Although at early times, the diverging pore still maintains a slightly more developed EDL near the entrance, compared to the converging pore.

The majority of this entrance behavior at early times can be attributed to the pore mouth size, consequently affecting the Biot number. The fixed length scales of the SDL region result in a higher Biot number for narrow entrances. This enhances the entrance flux, leading to more rapid charging near the pore mouth, as observed in the diverging case. A high Biot number effectively reduces entrance resistance, enabling greater electrochemical flux into the system. Despite the higher initial entrance flux observed in the diverging case, the overall charging rate across the entire pore is lowera result influenced by several factors that will be discussed later.

During later times, as shown in Figureb,d,f, the trends observed at early times persist and become more substantial. During intermediate times, the difference in charging rates between the diverging–narrow pores and the converging–wide pores becomes more pronounced. The convergent pore exhibits more well-formed electric double-layers compared to the other cases, even near the entrance. Additionally, we observe that the diverging pore shows a greater separation in the electrochemical potential profile than the narrow pore.

The change in charging rate for later time is likely due to a multitude of factors, but we will highlight a few mechanisms here. From eq, it is evident that the cross-sectional area plays a significant role in determining the electrochemical behavior within the pore. In the converging geometry, the cross-sectional area decreases along the pore length, while in the diverging case, it increases. Consider two connected differential slices within the pore. Any current not accumulated by the EDL would be affected by the locally changing cross-sectional area, altering the electrochemical flux. Specifically, if the adjacent region has a smaller area, the electrochemical flux must increase to maintain the current, as seen in the converging pore. This results in higher local fluxes and consequently, more rapid charging. Conversely, in the diverging pore, the increasing cross-sectional area causes a reduction in local flux, slowing the charging process. Furthermore, as the pore widens or narrows, the total double-layer capacitance changes accordingly, increasing with a wider cross-section and decreasing with a narrower one.

Figurec,d shows the simulated radial potential profiles at Z = 0.5 at early and intermediate times. At this location, the converging and diverging pores share the same radius, and at equilibrium the potential and charge density profiles would be identical. For convenience and readability, the radial coordinate is normalized by the local wall distance, α­(Z). At early times, we observe that the potential profiles are relatively undeveloped, remaining close to the applied potential. Additionally, all the potential profiles monotonically increase from the pore center toward the wall, reaching the applied wall potential Φ_ w _ ≈ 0.389. The largest potential difference from the centerline to the wall occurs for the wide and converging geometries, consistent with their faster charging dynamics. In contrast, the diverging and narrow pores exhibit a smaller potential difference, indicative of slower charging. At intermediate times, these trends persist. Although the relative difference between individual profiles decreases, the absolute potential drop from the center to the wall continues to grow as the EDL becomes more developed.

Similar to the potential profiles, the corresponding charge density distributions are shown in Figuree,f. These profiles provide insight into the structure and charge state of the local EDL. In early times, we observe that the converging pore exhibits the highest charge density, then followed by the narrow pore. In contrast, the diverging and wide pores show lower charge densities through the radial cross section. Notably, the wide pore displays a more developed EDL near the wall than either the narrow or diverging pore. The narrow and diverging pores have relatively flat charge density profiles, while the wide and converging geometries display steeper, monotonically decreasing profiles. This behavior is indicative of more developed double-layers near the wallwhich follows from their electrochemical potential profiles depicted in Figurea,bsuggesting that these pores are closer to equilibrium conditions. These trends persist at intermediate times, with the profiles becoming more pronounced and steeper, similar to the potential profiles.

A key strength of our reduced-order model is its accurate descriptions of the potential and charge density profiles. Figurec,d shows that the predicted potential profiles, given by simulating the full PNP equations, match very closely to those given by perturbation analysis with reduced computational cost. Similarly, Figuree,f demonstrates that the charge density recovered from the electrochemical potential agrees with directly simulating both ions using the PNP equations. Additionally, the constructed quantity μ̂the electrochemical potential of chargewhich is not explicitly included in the DNS, can be computed using eqs. As shown in Figurea,b, μ̂ demonstrates its effectiveness in characterizing the charging behavior, exhibiting strong agreement with the converted DNS data, and providing a useful charging metric.

To further validate our ROM and probe its limitations, we performed DNS of the nonlinear PNP equations for modified converging and diverging pores with larger changes in radii; see Figure. Due to numerical limitations, only one diverging pore geometry was tested. Simulation parameters were kept consistent with the other DNS cases, allowing for a systematic comparison with ROM predictions. Three converging pores were tested, all with entrance radii fixed at 20 nm and end radii of 10, 1, and 0.1 nm, representing 2×, 20×, and 200× decreases in radii. The 20 → 10 nm case corresponds to the original converging pore discussed previously. Because all converging pores share the same entrance size, they also share the same boundary condition, with Bi ≈ 7.95. The single diverging case tested had an entrance radius of 2 nm and a final radius of 20 nm; because its entrance size differs, the corresponding Biot number is Bi ≈ 795.

Electrochemical potential μ̂(Z,τ) and the electric potential Φ̂(Z,R,τ) plotted at various timesearly (a,c) and intermediate times (b,d)against direct numerical simulations with larger changes in slope. (a,b) Represent plots of the electrochemical potential as a function of axial position, and include the SDL region when Z < 0. The radial electric potential (c,d) profiles are plotted midway, (Z = 0.5), as a function of radial position. Three converging pores and one diverging pore are tested. The converging pores share the same entrance radius of 20 nm, but differing in their end radii. The 20–10 nm case corresponds to the original converging pore discussed previously, while the 20–1 nm and 20–0.1 nm pores explore steeper slopes, testing the limits of the perturbation analysis. Finally, the diverging pore is also added with a radius starting of 2 nm and an ending radius of 20 nm.

Despite these strong geometric variations, the ROM remains in close agreement with DNS, accurately capturing the evolution of both the electrochemical potential μ̂(Z,τ) and the electric potential Φ̂(Z,R,τ) across the pore. The ROM slightly underpredicts the charge state, with DNS showing the pores closer to equilibrium than the ROM predicts for the converging radii cases 20 → 1 nm and 20 → 0.1 nm. Additionally, the acceleration in the charging rate exhibits diminishing returns: the increase from 20 → 10 nm to 20 → 1 nm is larger than the subsequent increase from 20 → 1 nm to 20 → 0.1 nm, indicating a regime in which further geometric constrictions produce only marginal improvements. The steeper diverging case (2 → 20 nm) shows additional charging rate slowdown compared to the 10 → 20 nm diverging case (Figure) at the same nondimensionalized time. Additionally, the steeper diverging case exhibits a greater separation in charging rate relative to the original converging case, indicating that these widening pore openings substantially reduce the charging rate.

Overall, the DNS confirms that the ROM accurately captures the essential physics of pore charging and demonstrates that the electrochemical potential is a useful metric for quantifying charging dynamics in geometrically varying pores. While minor discrepancies exist between the ROM and DNS, particularly at later times, the close agreement at early times reinforces the model’s validity. Even with greater variations in radii, the ROM provides reliable qualitative insight into the influence of pore geometry on charging dynamics.

Geometry and

Charging

The formation of the EDL for the convergent and divergent geometries are depicted in Figure, where brighter contours indicate areas of high charge density. Figurea,d, which represent the early charging dynamics, show that the EDL penetrates deeper into the divergent pore compared to the convergent pore. For the intermediate charging dynamics, shown in Figureb,e, the EDL is more fully developed, as indicated by thicker contour sections in the convergent plot. The convergent pore exhibits a higher charge state with its thicker EDL. The equilibrium solutions are shown by Figurec,f. At this stage, the charge distributions are mirror images of each other, as the convergent and divergent pores are geometric reflections of each other. At equilibrium μ̂→0, resulting in time-independent charge and potential profiles from eqs

$$[eqn]$$ $$[eqn]$$

where α is a function axial position. Additionally, see supplemental Video S2 for animated visualization of charge density contours over time.

*Contour plots of charge density (ρ̂) throughout the pore at various times. (a–c) Depict the converging geometry, and (d–f) depict the diverging geometry, described by α­(Z) = 2 – Z and α­(Z) = 1 + Z, respectively. The charge density contours are shown at early (a,d), intermediate (b,e), and equilibrium (c,f) times. Lighter regions correspond to higher charge densities, while darker regions represent a more neutral state. τ denotes the nondimensional time, Z the nondimensional length along the pores, and R the nondimensional radius. The initial and system conditions are μ̂(Z,τ = 0) = 1, κ = 2, l

s

• = l

p

• and a

s

• = 4a

p

.

A useful and intuitive way to conceptualize this system is via the two different driving forces and their respective fluxes: the electromigrative (dark blue) and diffusive (light pink) fluxes. Figure illustrates a quiver plot of these two driving fluxes for κ = 2 at three representative time points: early, intermediate, and equilibrium. At the entrance, we find that electromigrative fluxes dominate over diffusion for the converging pore, whereas the electromigrative fluxes and diffusion are comparable for the diverging pore. This is expected because a smaller κα implies diffusion-dominated charging, whereas a larger κα implies that electromigration becomes the dominant transport mechanism.? Therefore, for κ = 0.1, both pores are dominated by diffusion, whereas for κ = 10, both pores are dominated by electromigration; see Supporting Information. We now emphasize why the charging of EDL speeds up for all κ values.

*Quiver plot of ∇ρ̂(Z,R,τ) and ∇Φ̂(Z,R,τ), representing the diffusive (light pink) and electromigrative (dark blue) fluxes, respectively. The flux profiles within the pores at early (a,d), intermediate (b,e), and equilibrium (c,f) times. The entrance-to-end ratio set at 2:1 for the converging pore and 1:2 for the diverging pore, defined by α­(Z) = 2 – Z and α­(Z) = 1 + Z, respectively. τ denotes the nondimensional time, Z the nondimensional axial coordinate, and R the nondimensional radial coordinate. The initial and system conditions are μ̂(Z,τ = 0) = 1, κ = 2, l

s

• = l

p

• and a

s

• = 4a

p

.

The electromigrative flux serves two purposes. First, it pushes the counterions toward the walls and the co-ions away from the walls, thus assisting with the formation of EDL. Second, electromigration also charges the pore axially, along with additional contributions from diffusion. Since the normal vector to the wall points toward the end of the pore for the converging scenario (Figurea–c), electromigration retains a favorable axial component, and the two purposes of electromigration assist each other. In contrast, for the diverging pore (Figured–f), the axial component of the electromigrative flux switches the direction, and the two purposes compete, slowing down the charging process. A similar trend is readily observable for κ = 10; see supplemental Video S1 and Supplementary Figure S2.

The discussion above highlights how a converging pore accelerates charging by using favorable electromigrative fluxes. However, for small κ values, electromigration ceases to exist; see Supplementary Figure S1. Even then, the converging pore charges faster. While the changing cross-sectional area and favorable electromigrative directions may initially appear to be two different effects, they are identical phenomena, since a varying cross-section necessarily affects the tangential wall component. Notably, the charging rate increases in converging geometries becomes more pronounced as the Debye length decreases.

We seek to quantify how the Debye length and entrance resistance interplays with charging rate; Figure demonstrates the influence of these two mechanism. Although the electrochemical potential of charge describes the “distance” from equilibrium, it is not analogous to total charge storage within the pore. Since the capacitance changes as a function of radiuslarger regions having higher capacitance and narrower regions with lower capacitancethe amount and distribution of charge are important considerations for the transient dynamics. For this, we introduce two new metrics Q(τ), the instantaneous total charge, and Q _ ss _, the equilibrium charge of the pore. Their ratio, Q/Q _ ss _, serves as a metric for a pore’s charge fraction. This charging metric provides a more comprehensive view of charging, as it represents a weighted average over all sections of the pore. These quantities can be simply derived from taking a volume integral of the average charge density at different points in time

$$[eqn]$$

*Plot of charge fractionquantified as total instantaneous charge (Q) over total steady state charge (Q SS)as a function of time (a), where μ̂(Z,τ = 0) = 1, κ = 2, and a

s

• = 4a

p

. The dotted line represents the τ1/2 or the time it takes a given pore to reach half of its equilibrium charge. For comparison, in plots (b,c), each τ1/2 is normalized by the half charge time of the narrow pore τ n . This ratio plotted on logarithmic scale versus κ for two different characteristic length scales of the SDL. These plots demonstrate the effect of entrance resistance; specifically, plot (b) exhibits a higher entrance resistance than plot (c).

From Figurea, the general trends in charging ratesconverging, wide, narrow, and then divergingcontinue to hold over time. Consistent with earlier observations, the narrow pore charges more rapidly at early times, but the wide pore eventually outpaces the narrow pore. The higher Biot number initially provides a boost to the narrow pore, but the wide pore’s lower resistance allows for faster charging at late times. Interestingly, their overall charging rates, at κ = 2, are relatively similar, unlike those of the converging and diverging pores. The converging pore demonstrates the most substantial improvements in charging rate compared to all other pore cases. In contrast, the diverging exhibits a significant decrease in charging rate. The charging rate decrease is even prominent at early times, where the higher Biot number benefits the narrow pore but not the diverging pore. This is because the total charge fraction includes contributions in the latter sections of the pore, where the narrow pore provides a constant resistance and capacitance throughout. In contrast, the diverging pore struggles to charge these regions throughout time. For all pores, the charging behavior is marked by a rapid initial charging phase, followed by a slow progression to equilibrium.

Important considerations are the effects of the Biot number and Debye length on charging dynamics. In Figureb,c, we compare charging rates with different length scales of the SDL, varying the relative dimensionless Debye ratio, κ, on a logarithmic scale. The SDL framework is a well-established technique for modeling the transition from reservoir to confined regions. ?,?,? Changing the characteristic length scale of the SDL affects the Biot number, which quantifies the resistance of electrochemical flux from the bulk region into the pore. We chose to change the SDL radial length scale, a _ s _ ^*^, which is related to Bi∝a _ s _ ^*2^ from eq. In the limit Bi → ∞, entrance resistance vanishes, whereas Bi → 0 corresponds to an infinitely resistive pore entrance. The influence of Biot number on the charging dynamics is nuanced, as it has a complex interplay with geometry and Debye length.

To evaluate charging rates, we introduce a new metric: the half-charge time, τ_1/2_, defined as the time required for a pore to reach 50% of its steady state charge. For clearer visualization and relative comparison, we normalize each pore’s τ_1/2_ by that of the narrow pore τ_ n _. In Figureb,c, we note two key limiting behaviors: the overlapping double-layer regime at low κ, and the thin double-layer regime at high κ, with a transition occurring around κ = O(1). We observe that high κ increase the charging rate generally, but the converging pore charges fastercompared to the narrow pore. Additionally, this effect is more prominent with lower entrance resistance, or equivalently, greater a _ s _ ^*^. This behavior stems from the reduced entrance resistance at high Biot, which permits greater ion flux into the porebenefiting geometries where transport is more limited by entrance effects, like the wide and converging pores. In contrast, diverging pores show minimal improvement, as their performance is primarily constrained by geometric inefficiencies rather than entrance resistance.

A surprising result emerges in the overlapping double-layer regime, where we found wide pores performing worse than narrow ones. While many studies highlight that wide pores tend to charge faster than narrow ones, ?,?,?,?,?

Figureb,c shows that entrance effects are an important consideration. While prior MD studies have identified nonmonotonic charging behavior in nanopores,? our work highlights a distinct continuum based mechanism of a similar phenomena. We attribute the reversal of the wide-pore charging trend primarily to the Biot number: by limiting the current entering the pore, the higher capacitance of the wider pore slows its charging rate. This effect manifests in the transition region from the thin to overlapping double-layer regime in Figureb,c, where charging rates are notably slower for the wide and converging cases for the higher resistance entrance, a _ s _ ^*^ = 4. While the entrance resistance strongly influences charging dynamics, it is not an easily tuned parameter, as it is set by several factors tied to the surrounding electrochemical environment. In practice, more controllable features lie within the electrode itselfmost notably the pore size and pore-size distribution. Because the entrance resistance in our model depends directly on the pore mouth radius, the entrance size provides a practical means to influence the Biot number. This highlights the importance of considering the interface between the electrode and the external electrolyte when evaluating charging behavior; in particular, determining the transport efficiency from the reservoir to the electrode, as it greatly affects the charging rate and which regime is limiting.

Conclusion

The result outlined in eq provides a differential equation that describes the charging of EDLs in a cylindrical pore with an axially varying pore radius. The reduced-order equation can recover all of the spatiotemporal details, while being 5–6 orders of magnitude faster in computation than the direct numerical simulations and provide mechanistic insight. Crucially, we find that a converging conical pore charges faster than both the wider and narrower pores for all Debye lengths. This acceleration occurs due to a boost to both electromigrative and diffusive fluxes as the cross-sectional area decreases. In the limit of thin EDL, where the analysis and impedance calculations of Keiser et al. hold,? a favorable electromigration flux direction further accelerates the charging process. Finally, this work highlights the importance of considering entrance resistance, as it can significantly impact electrode charging dynamics.

Looking forward, there exist numerous extensions and implications of this work. The most obvious extension is the expand this analysis to nonaxisymmetric pores and other shapes, and subsequently integrate these into a pore network model? with arbitrary geometries. This will allow one to explore how different shapes interact within complex, interconnected systemspotentially revealing structure–property relationships that do not emerge in the single pore case. For instance, Nguyen et al.? recently argued that electrode tortuosity factor may need to be adjusted; our work can allow such effects to be probed while capturing the effects of overlapping double-layers and pore shape changes. Additionally, investigating confinement effects in interconnected pore networks may enable the development of equivalent geometric representationssingle pore with an axially varying shape whose charging dynamics mimic those of the full network, allowing for more accurate and simplified models of real-world porous systems.

Beyond geometry, incorporating the effects for concentrated electrolytes, ?,? multi-ion systems,? diffusivity asymmetry, ?,? Faradaic reactions, ?,? and ionic liquids? opens several opportunities for inquiry. Beyond supercapacitors, our work can also advance understanding of electrolyte transport in membranes for separations? and biosensing applications.? Also, modeling of nanopores has shown potential in capacitive desalination and ion-selective nanopores, ?−? ? which shape could further manipulate the dynamics of these systems. Additionally, the response of conical geometries in an AC field has been shown to induce shape-dependent dynamicsparticularly relevant in the context of electrokinetic flows? where our theory could be employed to extend to arbitrary geometries.

Supplementary Material