The electroneutrality constraint in nonlocal models

In 1855, observations of water and nutrient transport through biological membranes prompted Adolf Fick to propose the phenomenological law, $J=−D∇C$⁠, which states that the flux of particles (J) from regions of high concentration to low concentration is proportional to the negative of the concentration gradient, $∇C$⁠.^1,2 The constant of proportionality D, which measures how quickly particles migrate, is the diffusion constant. Microscopically, Fick’s law arises from Brownian motion—the random zig-zag motion of particles at equilibrium due to thermal energy. The mean-squared displacement, $⟨x2(t)⟩$⁠, of an ensemble of such particles is proportional to time t; in one dimension, $⟨x2(t)⟩=2Dt$⁠, where D is the same diffusion coefficient.^1–4 This general relationship, $⟨x2(t)⟩=2Dt$⁠, holds for particles in homogeneous porous media [see Fig. 1(a)] saturated with a quiescent fluid.^5–11 In such porous materials, the correlation length of permeability fluctuations is small compared to the overall size of the domain and diffusion is predominantly Fickian.

Unlike normal or Fickian diffusion, anomalous diffusion is characterized by $⟨x2(t)⟩∼tα$⁠, where $α≠1$⁠.^12–14 Transport is sub-diffusive when $α<1$ and super-diffusive when $α>1$⁠. Anomalous diffusion is common in heterogenous porous media, such as naturally occurring aquifers and oil reservoirs, synthetic foams, and membranes. Due to heterogeneity in such systems, they exhibit variations in permeability starting from grain or pore size at the sub-micron level to the millimeter length scale. Sub-diffusive behavior, which is often observed in systems with a broad distribution of waiting times before random hops, is manifested when the trapping of a contaminant at stagnation points in groundwater flow leads to longer travel times than expected by normal diffusion.¹⁵ Another source of sub-diffusive transport in porous media is the adsorption and desorption of solutes onto surfaces.¹⁴ On the other hand, due to the presence of cracks or long channels in porous media that are comparable with the overall size of the domain, long-range correlations in the fluid velocity result in super-diffusion [Fig. 1(b)]. For such systems, Koch and Brady showed that $α$ is related to the covariance of the permeability field.^16,17

Anomalous diffusion arises from temporal or spatial nonlocality. Therefore, nonlocal formulations of mass conservation laws, which allow for non-Fickian diffusion, are potentially useful in the study of flow through heterogeneous porous media.^18,19 Similarly, nonlocal electrochemical models are of interest for investigating the behavior at the metal/electrolyte interface in problems such as corrosion.^20,21 There is evidence in the literature for anomalous dissolution of metals due to electrochemical corrosion,^22–24 and it has been reported that metal/electrolyte interfaces exhibit a nonlocal character.^25–28 There is increasing interest in developing models that can predict mass transport and electrochemical reaction kinetics in applications where such nonlocal behaviors exist.

Corrosion and groundwater transport are complex multiphysics systems that involve numerous physical processes such as networks of reactions, specific interactions, complex geometries, and moving boundaries that are often coupled with an underlying fluid or solid mechanics problem. In this paper, we isolate one physical process—electrochemical reactions in heterogeneous porous media—and ask the question “what is the best mathematical model (or governing equations) to describe charge transport in such systems?”

To do so, we consider a 1D nonlocal model for reaction and diffusion in multicomponent ionic solutions that describes the classic “liquid junction potential” problem in electrochemistry in which contact is initiated between two electrolytic solutions with different ionic compositions.^29,30 Two commonly studied variations include the (i) constrained-diffusion junction problem in which a porous membrane separates the two ionic solutions and (ii) free-diffusion junction problem in which there is no such membrane and direct contact is initiated.^31,32 Here, we consider the former, the constrained-diffusion liquid junction problem in which an electric field is induced across the membrane. This field adjusts the transport of charged species across the membrane until a stable concentration and potential profile is established. Before outlining the specific goals of the paper and developing the nonlocal formulation, we briefly review the classical (local) transport equations in multicomponent ionic solutions.

The conservation of mass in multicomponent ionic systems leads to the Nernst-Planck equation^30,32 that in the absence of convection, for an n-component system, takes the form

where C_j, J_j, and R_j represent the concentration, flux, and rate of generation of component $j=1,2,…,n$ in the solution, respectively. In dilute electrolytic solutions, the flux arises from concentration and potential gradients,

where $DjF$ and z_j represent the Fickian diffusivity and charge of component j, $ϕ$ is the potential, F is Faraday’s constant, R is the universal gas constant, and T is the absolute temperature.³³ To solve for the n + 1 unknowns (n concentrations and the potential $ϕ$⁠), the n equations in (1) have to be complemented by an additional relation. The electroneutrality condition provides this additional constraint. It is incorporated into the Nernst-Planck equation, usually in one of three different forms: (i) strict electroneutrality, (ii) charge conservation, and (iii) Poisson equation or Gauss law.

The most widely used relation in the literature for electrochemical modeling, due to its simplicity, is local electroneutrality. It asserts that the net charge density at any point is always zero, i.e.,

In this paper, we refer to this condition as strict electroneutrality (SE). Physically, SE contends that there is no charge separation anywhere in the solution. Empirically, SE is found to be a remarkably useful approximation which is valid in bulk solutions away from interfaces. Near electrodes and other charged surfaces, SE is violated in the thin double layer at length scales comparable with the Debye length. This length scale typically ranges from 1 to 10 nm and depends on the ionic strength of the solution.^31,34 SE provides a convenient shortcut and is widely used in electrochemical modeling; however, it is not a fundamental law. It is a poor approximation near interfaces at short length scales in weak ionic solutions. As noted by Planck over a hundred years ago,^31,35,36 the assumption of SE also presents a logical paradox: if there is no net charge anywhere in the domain, where does the electric field, $−∇ϕ$⁠, arise from?

Alternately, electroneutrality is sometimes modeled by asserting that the net current density due to moving ions at any point is zero

The relationship between Eq. (4) and SE can be made explicit by multiplying the $jth$ equation in Eq. (1) by z_j and summing up. From stoichiometry, chemical reactions do not produce net charge, hence $∑zjRj=0$⁠. Therefore,

If the net current density is zero (I = 0), then $∂Z/∂t=0$⁠. Here, Z(x, t) = Z(x, 0); the net charge density profile does not change with time. In this scenario, if the solution obeys SE initially, Z(x, 0) = 0, then it continues to obey SE, indefinitely. Conversely, SE (Z(x, t) = 0) implies $∇⋅I=0$⁠. If no external current is applied at the boundaries of the domain, it can be shown that this implies I = 0. Thus, SE and the “no net current” condition are equivalent under most common scenarios.³⁷ 

A less restrictive model for electroneutrality can be obtained from Eq. (5) by assuming that capacitive effects are negligible. Then, $∂Z/∂t=0$⁠, which yields the principle of charge conservation (CC),

This is more general than the no net current condition in that it allows for externally applied currents, which are important in certain processes such as corrosion. In contrast to SE which is an empirically motivated shortcut, CC can be derived from the Nernst-Planck equation. The principle of CC allows for charge separation and can be a reasonable choice for electrochemical modeling.³⁸ 

Besides SE and CC, a third method for complementing the Nernst-Planck equation is a Poisson equation which is sometimes referred to in the electrochemical modeling literature as Gauss’ law (GL). Given a spatial distribution of charged ions, GL allows us to solve for the resulting electric field. The Poisson equation for electrostatics is given by

where $𝜖$ is the permittivity of the medium. Like charge conservation, GL is a fundamental law of nature that can be derived from Maxwell’s equations. In Eq. (7), $𝜖/F≪1$ is small; for a uniform medium, this implies $Z≈0$ except at short length scales where gradients in $ϕ$ may be significant. We revisit this point later in the discussion. Note that it is possible to recast GL in terms of the electric field and current.^33,39 In principle, GL and CC are equivalent under most common scenarios. However, depending on circumstances, it may be preferable to invoke one over the other.

The primary goal of this paper is to determine how to best model nonlocal charge transport. To address this question, we develop a nonlocal Nernst-Planck model and impose electroneutrality using nonlocal counterparts of CC and GL. We study their relationship to the condition of SE by considering a variety of different initial scenarios. These scenarios allow us to consider the relative merits of the two methods. We perform a simple stability analysis and use both explicit and implicit solvers to numerically solve the model.

Armed with insights from this important modeling problem, we can approach more complex experimental systems. We discuss implications for multiphysics problems such as corrosion-induced mechanical failure, where charge transport is only one of many coupled processes.^40,41 In a typical corrosion problem, a load bearing metallic solid is in contact with a salt solution, e.g., ship hulls in sea water. Electrochemical reactions between the ionic species and the metal result in the degradation of the solid due to corrosion which increases its susceptibility to failure by fatigue or fracture. Indeed, the development of peridynamics-based corrosion damage models,^20,21,40,41 where nonlocal mechanics models are tweaked to capture corrosion damage processes through suitable coupling with diffusive fields, has recently attracted considerable attention. Metal/electrolyte interfaces tend to exhibit nonlocal effects,^25–28 and there is evidence on the anomalous dissolution of metals due to electrochemical corrosion.^22–24 However, most peridynamics models published to-date for corrosion phenomena do not explicitly consider the evolution of corrosion electrochemistry in a nonlocal setting. We note that nonlocal diffusion behaviors have been investigated in different domains and with different constraints.^42,43 Recent efforts have also focused on interface problems in nonlocal diffusion⁴⁴ and coupling of nonlocal mechanics and diffusion⁴⁵ in the realm of peridynamics theory.

In the general corrosion-induced mechanical failure problem, nonlocal electrochemical reaction diffusion equations, such as those discussed in this paper, have to be systematically augmented with charge transport and momentum conservation equations.^46–49 Such studies are currently underway in our group and results will be reported in the future. The development of peridynamic models for treatment of the metal/electrolyte interface is a natural extension of classical electrochemical models. They can be applied to study phenomena such as corrosion, battery electrochemistry, and heterogeneous catalysis.

In this section, we describe a nonlocal reaction diffusion model in 1D. This is complemented by two different conditions for electroneutrality, viz. nonlocal forms of the charge conservation and Poisson equations. We assume there are three species, A, B, and the product P, dissolved in a homogeneous electrolytic solution, undergoing the reversible chemical reaction

Let C_j(x, t) denote the concentration of species $1≤j≤3$ at location x and at time t. In this paper, the 3 species are ordered as A, B, and P, respectively. Thus C₁ denotes the concentration of A, etc. The shorthand C(x) = {C₁(x), C₂(x), C₃(x)} is used to denote the set of concentrations at a particular location. The reaction rates are assumed to follow from the stoichiometry of the reaction. Thus, the rates of generation are given by

where k_f and k_b are the forward and backward reaction rate constants.

Here, we propose nonlocal counterparts of the Nernst-Planck equation, charge conservation, and Poisson equations.

The nonlocal reaction-diffusion model with homogeneous reactions may be written, for $1≤j≤3$⁠, as

where $δ$ is the nonlocal horizon, R_j is the rate of generation of j, and the flux density operator J_j encodes the transport of species j due to concentration and potential gradients. In Eq. (9), $Jj(x,x^,t)$ depends explicitly on the concentrations C_j(x, t) and $Cj(x^,t)$ and the potential difference $Δϕ(x,x^)=ϕ(x^)−ϕ(x)$ between the points $x^$ and x. Temporarily suppressing the dependence on time in the notation and letting $ΔCj(x,x^)=Cj(x^)−Cj(x)$ and $C¯j(x,x^)=(Cj(x)+Cj(x^))/2$⁠, the flux density may be written by suitably adapting the Nernst-Planck equation as

where z_j is the charge on species j, D_j is the effective diffusion coefficient, R is the universal gas constant, F is Faraday’s constant, and T is the absolute temperature. The kernel $K(x,x^)$ may be chosen to reflect Fickian or for anomalous diffusion,⁵⁰ 

As described in the Introduction, anomalous diffusion may be appropriate for transport in heterogeneous membranes and porous media. For instance, superdiffusion has been observed in corrosion-related transport,^22–24 in which the MSD $⟨x2(t)⟩∼tα$ with $α>1$⁠. s in the kernel is related to the $α$ via $s=1/α$⁠. The nonlocal formulation described in the paper can handle the full range of $s∈(0,1)$⁠. In this work, we set s = 1/4. One reason for choosing s = 1/4 is to explore this superdiffusive regime that cannot be easily accessed using conventional models. However, insofar as the goals of the paper are concerned (discrimination between CC and GL), the choice of s does not affect any of the conclusions drawn.

We render (9) and (10) dimensionless. In 1D, the dimensions of C_j and J_j are mol/m and mol/m²s, respectively. The units of the effective diffusivity, which is related to Fickian diffusivity via $Dj=DjF/δ2$⁠, are $s−1$⁠. Similarly, the units of the kernel K and the factor RT/F are $m−1$ and volts, respectively. Thus, we can non-dimensionalize (i) the concentration with a representative concentration C₀ so that $Cj*=Cj/C0$⁠, (ii) time and diffusivities with a representative effective diffusivity D₀ so that t^* = tD₀ and $Dj*=Dj/D0$⁠, (iii) lengths with the effective domain size L so that x^* = x/L and $δ*=δ/L$⁠, (iv) the potential with RT/F so that $ϕ*=Fϕ/RT$⁠, (v) the kernel with L so that K^* = KL, (vi) the reaction rates with C₀D₀, so that $Rj*=Rj/(C0D0)$⁠, and (vii) the flux operator by C₀D₀/L so that $Jj*=JjL/(C0D0)$⁠. The nondimensional reaction diffusion equations can then be written as

with

Because the current density is $I*=∑zjJj*$⁠, the nonlocal counterpart of the charge conservation (CC) relation may be written in dimensionless form as

Equation (14), like Eq. (6), is valid only when capacitive effects are negligible and the spatial distribution of charge density is time invariant.

The nonlocal version of Gauss’ law (GL), Eq. (7), can be written in dimensionless form as

where the dimensionless permittivity, $𝜖*=𝜖(LRT)/(F2C0)$⁠.

Equations (12) and (13) define the nonlocal Nernst-Planck relations that must be complemented by a nonlocal electroneutrality condition either via the CC [Eq. (14)] or the GL [Eq. (15)] relations.

The geometry for the problem corresponds to the liquid junction potential problem in electrochemistry for which we assume that two dissimilar electrolytic solutions are brought into contact via a permeable thin film membrane. The spatial domain for the computation is the thin film between the two bulk solutions that extends over $x*∈[−δ*,1+δ*]$⁠. “Boundary” conditions corresponding to the two different bulk solutions are imposed over the intervals $[−δ*,0)$ on the left and $(1,1+δ*]$ on the right.

In the base case, $Δx*=0.05$⁠, the horizon $δ*=3Δx*=0.15$⁠, and $Δt*=0.05$⁠. The reaction rates are $kf*=kf/C0D0=0.1$ and $kb*=kb/C0D0=0.1$⁠, whereas the dimensionless permittivity $𝜖*=1000$⁠. The index, charge, and dimensionless diffusion coefficients of the 3 dissolved species are listed in Table I.

Dirichlet conditions are specified at the left and right boundary regions. In $[−δ*,0)$⁠, we set the bulk solution concentrations of A and B equal to 0.3 and 0.6, respectively, so that the strict electroneutrality condition, $∑zjCj*=0$⁠, is satisfied. For the same reason, on the right boundary $(1,1+δ*]$⁠, we set the concentrations of A and B equal to 0.5 and 1.0, respectively. In both the left and right boundary regions, the concentration of P is set to $C3=(kf*/kb*)C1*C2*2$⁠, where C₁ and C₂ are the concentrations of A and B, respectively. This ensures that the concentrations of the three species are in chemical equilibrium in both boundary regions. Furthermore, the potential $ϕ*=0$ in both boundary regions.

We consider four different sets of initial conditions; these cases are summarized in Fig. 2. The four initial conditions are chosen not to mimic real systems but to help discriminate between charge using conservation (CC) and Gauss’ law (GL). The goal is to consider a parsimonious set of initial conditions from which general inferences regarding the differences between CC and GL can be drawn.

In case 1, we impose a quadratic profile for A on the domain [0, 1], $C1*(x*,0)=0.55+0.2x*−(x*−0.5)2$⁠, which is consistent with the boundary region conditions. For B, we set $C2*=2C1*$ so that the SE condition is satisfied and the net charge is zero. Similarly, we set $C3*=(kf*/kb*)C1*C2*2$ so that the solution is in local chemical equilibrium throughout the domain.

In case 2, the A and B concentration profiles are identical to those in case 1. Therefore, the solution obeys SE. However, the concentration of P is varied linearly across the solution domain. Thus $C3*$ is not in local chemical equilibrium with A and B.

Case 3 is somewhat artificial. Here, we impose a quadratic profile for A as in case 1; however, we set a linear profile for B, $C2*(x*,0)=0.6+0.4x*$⁠. Because $C2*≠2C1*$⁠, there is excess positive charge on the domain. We also set $kf*=kb*=0$ so that chemical reactions are suppressed, i.e., $Rj*=0$⁠. Thus, in this case, the component P does not participate in the calculation so we set $C3*(x*,0)=0$⁠.

Finally, for case 4, we set the concentration of all the species in the solution domain equal to zero. This corresponds to a physical situation in which two electrolytic solutions in chemical and charge equilibrium are first brought into contact through a pristine permeable membrane. Note that SE and chemical equilibrium are “satisfied” inside the solution domain.

To summarize, the initial condition in cases 1 and 4 obey both SE and chemical equilibrium. Case 2 satisfies SE but not chemical equilibrium, while case 3 violates SE and has no chemical reaction.

We use two methods to solve the Nernst-Planck equations numerically: an explicit and an implicit scheme. We describe the explicit scheme first because it is easier to articulate. We discretize the time derivative in (12) via

We divide the spatial domain into a uniform grid with spacing $Δx*$ and evaluate the integral over the flux density operator at a grid point $xi*$ using the Riemann sum,⁵¹ 

To propagate the solution from t^* to $t*+Δt*$⁠, we use Euler’s method.⁵¹ In the explicit method, concentration and potential values from the previous time step t^* are used on the RHS in the nonlocal Nernst-Planck equation. This allows us to determine all the concentrations at the new time $C(t*+Δt*)$ explicitly. We then solve the appropriate electroneutrality equation to find the updated value of the potential $ϕ*(t*+Δt*)$⁠.

For charge conservation, Eq. (14), we use a similar Riemann sum to resolve the integral, which yields

Equation (18) leads to one linear equation at each grid point $xi*$ in the computational domain. We can solve the resulting linear system to obtain the profile for the potential. Letting $ϕi*=ϕ*(xi*)$ and rearranging terms, we obtain the banded linear system

The inner summation is over the n = 3 chemical species and the outer summation is over the nonlocal horizon around $xi*$⁠. Here,

For Gauss’ law, Eq. (15), we get a similar banded system when we discretize

The linear system (19) or (22) can be solved to obtain $ϕ*$ at the internal nodes over the entire domain. To summarize, the explicit method of propagating the solution involves repeatedly the following:

• Solving for the current concentrations $C*(x*,t*+Δt*)$ of all species based on the concentrations C^*(x^*, t^*) and potential $ϕ*(x*,t*)$ from the previous time step using (16) and (17) in (12).

• Solving for the potential $ϕ*(x*,t*+Δt*)$ by solving the linear system corresponding to the CC (19) or GL (22) electroneutrality conditions.

In the implicit scheme, instead of using species concentrations and potential from the previous time in the first step, we use their current values. This involves setting up a nonlinear system of equations in which we guess $C*(x*,t*+Δt*)$⁠, solve the linear system for $ϕ*(x*,t*+Δt*)$⁠, and refine our guess until Eq. (12) is satisfied. We use the Powell hybrid method as implemented in MINPACK to solve the nonlinear system with a tolerance of 10⁻⁵.⁵² A SciPy (v0.16) interface to the method is available through the optimize module in Python.^53–55 

Starting from the initial conditions depicted in Fig. 2, we carry out simulations up to a final time of $tf*=50$⁠. All the cases studied evolved to their corresponding steady states (all time derivatives, e.g., $∂C*/∂t*≈0$⁠) over this time frame. We expect the rates of approach to a steady state for the CC and GL electroneutrality conditions to be different as they are controlled by different parameters D^* (CC) and $𝜖*$ (GL). Nevertheless, the steady state profiles do not depend on these parameter choices and allow us to compare the equilibrium states reached by the two models.

We use a Riemann sum to numerically evaluate integrals for which the error is proportional to the grid spacing $Δx*$⁠. To propagate the solution in time, we use Euler’s method for which the global error is proportional to $Δt*$⁠.⁵¹ The base case choice of $Δx*=0.05$ and $Δt*=0.05$ was found to be adequate; decreasing $Δx*$ and $Δt*$ further did not produce any observable improvement in the steady state solution.

First we discuss results of simulations using the explicit scheme for each of the four different initial conditions. This is followed by a discussion of the relative merits of the two electroneutrality constraints: charge conservation and Gauss’ law. Finally, the stability of the explicit method and comparison with the implicit scheme are reported in the  Appendix.

Figure 3 depicts the initial and final profiles of $Cj*$ and $ϕ*$ using the two different criteria for electroneutrality. Let us compare $ϕ*$ computed from the initial concentration profiles by CC and GL. Because strict electroneutrality is imposed initially, there is no net charge, Z^*(x^*, 0) = 0. In the case of GL [Fig. 3(c)], this leads to $ϕ*(x*,0)=0$⁠, due to the boundary conditions. Note, however, that the initial potential profile in Fig. 3(a) corresponding to CC is slightly negative in the interior regions of the domain. Although SE is satisfied, the diffusivities of A and B are unequal $D1*≠D2*$⁠. Thus, $∑zjDj*Cj*≠0$ which results in a small, but nonzero $ϕ*(x*,0)$ as $Gjk≠0$ in Eq. (19). If SE is satisfied and the diffusivities of the charged species are equal, then CC leads to $ϕ*(x*,0)=0$⁠. In contrast, diffusivities do not enter GL directly.

As the system evolves, both GL and CC electroneutrality constraints lead to nearly identical profiles for the concentrations and the potential. SE is satisfied and the potential remains or approaches zero throughout the domain. Thus, for a system that initially obeys SE and is in chemical equilibrium, both electroneutrality criteria yield comparable results.

In this case, the concentration of P is not in chemical equilibrium with A and B. Because the initial concentration profiles of the charged species are the same as in case 1 and P is uncharged, the initial potential profiles in Fig. 4 are identical to those in Fig. 3.

As the system evolves, SE is essentially satisfied for both the CC and GL conditions. Indeed, the similarity of Figs. 4(b) and 4(d) allows us to make a fairly general claim. If SE is initially satisfied throughout the domain, then it does not matter whether we choose CC or GL as our condition for electroneutrality. There is one important exception to this claim which is discussed through case 4. Otherwise, both electroneutrality conditions yield nearly identical concentration profiles as the system evolves to a steady state. The next two cases highlight situations in which the results for CC and GL diverge.

Figure 5 depicts the initial and final profiles of $Cj*$ and $ϕ*$ for case 3. The concentration of P is zero throughout the calculation, and the chemical reaction does not enter the picture. Initially, there is a nonzero net charge distribution due to an excess of positive A ions. When the CC electroneutrality condition is invoked [Fig. 5(b)], the net charge profile persists. On the other hand, when the GL condition is imposed [Fig. 5(d)] both Z^* and $ϕ*$ decay to zero.

For CC, this case is artificial because that condition relies on the assumption that capacitive effects are negligible: the medium is unable to store charge. The initial condition is thus incongruent with this assumption. The response of the CC condition, while nonintuitive, follows the constraints placed on it. The net charge density does not change with time, i.e., $∂Z*/∂t*=0$⁠, and the medium is unable to discharge previously accumulated charge. This is reflected in Figs. 5(a) and 5(b): the net charge profile is time invariant.

It is interesting to note that the initial profiles of $ϕ*$ resulting from the CC and GL conditions are of opposite signs. CC yields a $ϕ*$ profile that “justifies” the accumulation of net positive charge. The trough in $ϕ*$ centered around x^* = 0.5 describes the accumulation of positively charged ions. Because SE is not satisfied at t^* = 0 and the net charge density is time invariant, CC violates SE throughout the simulation. This is in direct contrast to cases 1 and 2 for which SE is always satisfied.

When the GL condition is used, the potential $ϕ*$ reflects the accumulation of positively charged ions. The resulting electric field, which corresponds to the negative gradient of $ϕ*$⁠, pushes A towards the boundary. This current eventually results in linear steady state concentration profiles. As the net charge approaches zero, the system starts to obey SE on the domain.

Figure 6 presents an initial condition in which the concentration of the dissolved species in the domain is zero. In this scenario, we cannot use the CC condition even though SE is satisfied. This is due to the structure of the linear system that arises from the CC electroneutrality condition. When all the concentrations at any location are zero, the terms H_jk in Eq. (19) corresponding to that node become zero. The linear system is underspecified, does not have full rank, and hence cannot be solved. This problem is encountered even when the sum of the concentrations is small but non-zero; it is manifested as an ill-conditioned linear system. On the other hand, this scenario does not present any problems for GL. The linear system corresponding to GL remains well specified even when concentrations are zero.

Cases 3 and 4 show how GL may be used in situations where CC fails for structural or numerical reasons. In Case 3, we admittedly used an artificial initial condition for which CC behaved as expected. On the other hand, case 4 represented a perfectly natural initial setting. It highlighted an important limitation of CC: low or zero concentration at any point on the solution domain. Low concentrations are relevant in corrosion problems in which electrolyte solutions gradually infiltrate metal electrodes and initiate degradation.

GL has another significant advantage over CC. Both CC and GL lead to banded linear systems, $Mϕ*=r$⁠. The bandwidth of the systems is identical and depends only on the horizon $δ*$⁠. However, the elements in the matrix M are constant for GL; for CC they depend on the concentrations and hence need to be updated at each time step. For GL, this structure can be exploited to store or factorize the matrix M only once, making it much cheaper to evaluate. Thus, a case can be made in favor of GL over CC on grounds of both efficiency and versatility.

However, we glossed over an important limitation of Gauss’ law. In our calculations, we set the parameter $𝜖*=1000$⁠. For a medium like water, the nonlocal permittivity is, in fact, $𝜖∼10−10/δ$ C/V m. If we assume $δ∼10μ$m, then $𝜖∼10−5$ C/V m². The nondimensional permittivity $𝜖*=𝜖LRT/F2C0$⁠; if we use reasonable scales for $L∼0.1$ m, RT/F = 25 meV, and $FC0∼10$⁠, we obtain a range of $𝜖*=10−7–10−10$⁠. This is much smaller than the value of $𝜖*=1000$ that we used in our calculations. Physically, the small value of $𝜖*$ indicates that even a modest separation of charge (⁠$Z*≠0$⁠) would induce extremely large electric fields that would quickly neutralize any excess charge. Indeed, for liquid junctions, it has been shown that the relevant time scale in which much of this action occurs is of the order $τ=λD2/DF$⁠, where $λD≈10−7$ cm is a typical Debye length, and D^F is the Fickian diffusivity, which is of the order of 10⁻⁵ cm²/s.³⁶ Thus, $τ∼10−9$ s, which means that the net charge density rapidly approaches zero for time scales associated with experimental observations.

Why then did we not use a small $𝜖*$ in our numerical calculations based on GL? For scenarios in which SE is strictly satisfied, the value of $𝜖*$ is irrelevant because the Poisson equation reduces to the Laplace equation. In cases 1, 2, and 4, even though SE is satisfied and $Z*≈0$ is small, the linear system corresponding to GL can become susceptible to round-off error because $Z*/𝜖*$⁠, which is a ratio of two small quantities, can be noisy. In these cases, where SE is essentially satisfied, using an artificially large value of $𝜖*$ helps denoise or regularize the solution, without significant loss in accuracy.

In case 3, where SE is clearly violated, the true dynamics of the system are much faster than those that arise from using an artificially larger $𝜖*$⁠. Again, this is because a smaller $𝜖*$ creates stronger electric fields that neutralize charge disparities rapidly. In particular, the nonzero Z^* profile approaches zero and satisfies SE more quickly when the true $𝜖*$ is used. Once SE is satisfied and $ϕ*≈0$⁠, the dependence on $𝜖*$ becomes less important.

The potential $ϕ̃*$ inferred using a large $𝜖̃*$ differs from the true potential $ϕ*$ corresponding to the true $𝜖*$ by a multiplicative constant, i.e., $ϕ*=(𝜖̃*/𝜖*)ϕ̃*$⁠. Therefore one can circumvent the noise issue by rescaling $ϕ̃*$ to the true potential. However, this might imply that the potential gradient term in the Nernst-Planck equation dominates the concentration gradient term when $(𝜖̃*/𝜖*)$ is large. This may require a small time step and/or an implicit scheme for time propagation. An adaptive integration scheme that uses a small time step when SE is violated and transitions to a larger time step as $Z*→0$ could also be considered.

The ability to describe physical processes involving anomalous diffusion and heterogeneous transport of species is partially responsible for the increasing popularity of nonlocal models. Science and engineering applications with such processes include the following: charge transport through porous media, corrosion of metals, and ion/energy transport in biological membranes and cells. For many of these multiphysics problems, the appropriate prescription of electroneutrality is important.

In this paper, we addressed this problem of prescribing the electroneutrality constraint in nonlocal models that involve diffusion and reaction of charged species. First, we presented a simple nonlocal model corresponding to the Nernst-Planck equations for reaction and diffusion in multicomponent ionic systems. To close the set of equations, we proposed nonlocal counterparts of the charge conservation (CC) and the Gauss law (GL). Subsequently, we applied this model to describe the 1D liquid junction problem in which two electrolytic solutions are brought into contact via a permeable membrane. Using four different initial conditions for the concentration profiles in the permeable membrane, we investigated the relative merits of CC and GL as closure constraints that account for electroneutrality in different scenarios. The resulting systems of equations were solved using both explicit and implicit methods. For the time step $Δt*=0.05$⁠, which is much smaller than the critical time step for the onset of instability $Δtmax*≈0.7$⁠, the results of the explicit and implicit methods agreed with each other.

When strong electroneutrality is satisfied in the initial condition, results of CC and GL essentially agreed with each other. In these cases, using an artificially large value of the permittivity constant $𝜖*$ helps us to denoise and smoothen the linear system used to determine the potential, without significantly perturbing the evolution of the concentration profiles. Differences between CC and GL were observed when SE was not initially satisfied or when the concentration of all the charged species dropped to zero at any point in the domain. In the former case, the CC equation assumes that there are no capacitive effects and hence does not allow any built-up charge to discharge. The latter case poses a bigger problem, especially in problems like corrosion modeling where the concentrations of the infiltrating electrolyte is necessarily small. In these cases, the Gauss law performs as expected. In addition to its versatility, it is also cheaper to solve for the potential using GL because the matrix in the linear system to be solved at every time step does not change with time. Ideas on treating small values of $𝜖*$ are also presented. Insights developed from this study could be helpful in the development of nonlocal (peridynamics-type) models for phenomena like corrosion damage and charge transport across interfaces or membranes.

This work was funded by Naval Air Systems Command (NAVAIR) through a Small Business Technology Transfer (STTR) program, Contract No. N68335-15C-0032, awarded to Advanced Cooling Technologies, Inc. (ACT) in collaboration with Florida State University (FSU). The authors would like to thank the technical advice and support of Dr. Kishan Goel and Dr. Nam Phan, both from NAVAIR/Naval Air Systems Command, Patuxent River, MD.

In this appendix, we discus the stability of the explicit numerical method used in this work and compare it with the implicit scheme described earlier. For explicit schemes such as those used in the paper, we can perform a stability analysis to determine an upperbound for the time step, $Δt*$⁠, in terms of the grid spacing and other parameters. Figure 7 shows numerical experiments with different $Δt*$ for initial conditions corresponding to case 1 with the GL electroneutrality condition, using both implicit and explicit schemes.

In Sec. III, we used $Δx*=0.05$ and $Δt*=0.05$⁠, throughout. In Fig. 7, we reproduce two of the curves reported previously in Fig. 3(d) for the concentration of component A and the net charge Z^*. For the explicit scheme, as the time step is increased from $Δt*=0.05$⁠, the earliest signatures of instability in both $C1*$ and Z^* at a final time of $tf*=50.0$ are not clearly visible until $Δt*≈0.70$⁠. When we use a fully implicit scheme, with the same large time step $Δt*=0.70$⁠, the instability disappears, and the concentration and charge profiles essentially overlap with the $Δt*=0.05$ explicit or implicit scheme. Thus, from these numerical experiments, it appears that $Δt*=0.05$ is sufficiently smaller than the time step $Δt*≈0.70$⁠, at which the scheme becomes unstable.

We investigate the onset of instability more systematically, by performing a linear stability analysis on a simpler version of the problem. From previous results, we note that $ϕ*≈0$⁠, when the initial conditions obey SE (cases 1, 2, and 4). We also ignore the nonlinear reaction term and consider the following nonlocal diffusion problem for this analysis:

Let $Cjn$ represent $C(x,t)=C(jΔx,nΔt)$⁠. Then, the discretization used for the explicit scheme implies

We can perform a von Neumann stability analysis by setting $Cjn=ζneiθj$⁠, where $i=−1$⁠. Using $δ=mΔx$ and $xj=jΔx$⁠, in the equation above and simplifying, we get

The term in the summation

Each term in the sum over l in $F(θ;m,s)$ lies between 0 and 1. Thus, a conservative upperbound for $F(θ;m,s)$ is given by $Fmax=m$⁠. For stability, we require $|ζ|<1$⁠, which after simplification leads to $DΔtm2sFmax<1$ or

For m = 3 and s = 1/4 as used in this work, it can be shown that $Fmax≈1.2$⁠, which gives a more realistic upperbound. Using this value for $Fmax$ in Eq. (A5), it can be shown that the maximum value of $Δtmax*≈0.16$⁠, which is more than 3 times larger than the value of 0.05 used in Sec. III. An interesting feature of this analysis is that $Δtmax*$ depends on $Δx*$⁠, indirectly through m, i.e., for a fixed m, $Δtmax*$ is independent of $Δx*$⁠.

We can test this result numerically by suppressing the reaction and potential terms to make the calculations correspond to the assumptions made in the stability analysis. Using case 1 as the basis, we made two additional simplifications: (i) we set $kf*=0$ so that the reaction term (and species C) drops out and (ii) set $ϕ*=0$ throughout the calculation. Thus, we end up with only two components A and B, diffusing independently through a nonlocal diffusion equation. Because $ϕ*=0$ is held fixed, we do not need an additional electroneutrality constraint.

For this simplified problem, using $Δx*=0.05$ and m = 3, as in the base case, we found that $Δt*$ had to be increased to 0.75 before signs of instability were visible in the concentration profile of A, at $tf*=50.0$ (Fig. 8). For this value of the time step $Δt*=0.75$⁠, changing the spatial resolution 10-fold between $Δx*=0.01$ and $Δx*=0.10$ did not alter the time step required for the instability to manifest.