Solving the Nernst‐Planck Equation in Heterogeneous Porous Media With Finite Volume Methods: Averaging Approaches at Interfaces

Molecular diffusion of dissolved species is a fundamental mass transport process affecting many environmental and technical processes. Whereas diffusive transport of single tracers can be described by Fick's law, a multicomponent approach based on the Nernst‐Planck equation is required for charge‐coupled transport of ions. The numerical solution of the Nernst‐Planck equation requires special attention with regard to properties that are required at interfaces of numerical cells when using a finite difference or finite volume method. Weighted arithmetic and harmonic averages are used in most codes that can solve the Nernst‐Planck equation. This way of averaging is correct for diffusion coefficients but inappropriate for solute concentrations at interfaces. This averaging approach leads to charge balance problems and thus to numerical instabilities near interfaces separating grid volumes with contrasting properties. We argue that a logarithmic‐differential average should be used. Here this result is generalized, and it is demonstrated that it generally leads to improved numerical stability and accuracy of concentrations computed near material interfaces. It is particularly relevant when modeling semipermeable clay membranes or membranes used in water treatment processes.


Introduction
Diffusion of aqueous species in geological or engineered media is a fundamental mass transport process. It is especially important for low permeability geological materials containing significant amount of clay minerals such as clayey shales, engineered materials such as clay barriers, or concrete structures. Their low permeability and diffusion properties make them ideal for waste confinement applications or technological materials such as filtration membranes used for water treatment. The characterization of diffusion processes is also essential for our ability to understand various hydrogeochemical observations such as isotopic fractionation coupled to transport processes (Bourg et al., 2010;Bourg & Sposito, 2007;La Bolle et al., 2008;La Bolle & Fogg, 2001;Peeters et al., 2002;Rolle et al., 2010), the dynamics of gas-water exchanges (Haghighi et al., 2013;Tokunaga et al., 2017), or the dynamics of contaminant accumulation and release in and from rocks and sediments having very heterogeneous pore structures (Bone et al., 2017;Chapman & Parker, 2005;Gouze et al., 2008;Hadley & Newell, 2014;Liu et al., 2006Liu et al., , 2011Robinet et al., 2012;Zachara et al., 2016). Ultimately, diffusion is the fundamental process that generates mixing of dissolved species and enables reactive fronts to appear between aqueous solutions having contrasted chemical compositions (de Anna et al., 2011(de Anna et al., , 2013Le Borgne et al., 2011. Diffusion processes are the result of random motion of dissolved species subject to thermal agitation and for which no interactions between the dissolved species are considered. Diffusion processes are commonly simulated with Fick's laws. However, ions are charged species, and their individual diffusion coefficients in solution are dependent on their charge, mass, and radius. As a consequence of the electroneutrality condition in aqueous environments, ions are affected by electrochemical migration effects, and multicomponent diffusion processes are thus better represented by the more general Nernst-Planck formulation rather than by the limiting Fick's laws. The importance of electrostatic interactions among charged species in the modeling of multicomponent diffusion processes was early emphasized to explain vertical profiles of ion concentrations in the pore water of marine sediments, that is, systems in which diffusion is the dominant transport process (Ben-Yaakov, 1972;Boudreau et al., 2004;Felmy & Weare, 1991a, 1991bGiambalvo et al., 2002;Lasaga, 1979). Later, the importance of multicomponent diffusion in our understanding of mixing processes in porous media has been demonstrated even for systems whose mass transport is dominated by advective flow (Chiogna et al., 2011;Muniruzzaman et al., 2014;Muniruzzaman & Rolle, 2015, 2017Rasouli et al., 2015;Rolle et al., 2018). In the field of reactive transport modeling, the use of multicomponent diffusion models is hindered by two factors: The first one is the scarcity of codes that are able to handle the Nernst-Planck formulation for the resolution of diffusive fluxes ; the second one is the computational cost associated with the use of the Nernst-Planck formulation rather than Fick's laws. In the last decade, the use of Nernst-Planck equation instead of Fick's laws has been shown to be essential to understand the apparent anomalous diffusion behavior of systems in which the diffusion of charged species is affected by the electrostatic properties of the surfaces present on the solid phases . Most of the related studies concerned the properties of clay and concrete materials, which are investigated with regard to their confinement properties for radionuclides or other toxic solutes (Appelo & Wersin, 2007;Appelo et al., 2008Appelo et al., , 2010Appelo, 2017;Alt-Epping et al., 2015Bourg & Tournassat, 2015;Gvirtzman & Gorelick, 1991;Glaus et al., 2013Glaus et al., , 2015Tinnacher et al., 2016;, 2019a, 2019b. However, the use of reactive transport models using the Nernst-Planck formulation can be foreseen to be increasingly important for the modeling of other types of materials and related applications including microbial electrochemical cells or membrane filtration technologies (Marcus et al., 2010).
The numerical solution of the Nernst-Planck equation in a reactive transport model using a finite difference/volume method is subject to a range of difficulties when applied to spatially heterogeneous media . In this study, we address the problem of the definition of averaged properties at the interface between porous domains having contrasting properties. This work should facilitate a rigorous implementation of the Nernst-Planck equation in reactive transport codes.

Governing Equations
In absence of an external electric potential, the electrochemical potential μ i (J mol −1 ) of an ion i can be expressed as where T is the temperature (K), R is the gas constant (8.314 J·K −1 ·mol −1 ), F is the Faraday constant (96,485 J·V −1 ·mol −1 ), ψ an (internal) electrical potential (V), m o is the standard state molality (1 mol kg −1 ), μ o i is the standard (electro)chemical potential of species i (J mol −1 ), a i is its chemical activity, z i is its charge number (dimensionless), m i is its molality (mol kg −1 ), and γ i is its activity coefficient (dimensionless). The diffusive flux J i,s (mol·m −2 ·s −1 ) of an ion in a solution is given by the Nernst-Planck equation: where u i is the mobility of species i (mol·m 2 ·s −1 ·J −1 ) and c i is its molarity (mol m −3 ), which can be expanded as: where ρ solv is the density of the solvent (in kg solvent m −3 solution ). The mobility u i refers here to the average velocity of a species in solution acted upon by a unit force, independent of the origin of the force (Steefel & Maher, 2009). The diffusion coefficient D i (m 2 s −1 ) of the species i is proportional to its mobility according to the Nernst-Einstein equation: In a porous medium, the diffusion coefficient of the species i is usually described as a function of the porosity ϕ, of the tortuosity factor τ i of the medium, which can be specific to each species, and of the self-diffusion coefficient of the species in solution D i,s (Shackelford, 1991): 10.1029/2019WR026832

Water Resources Research
TOURNASSAT ET AL.
The diffusive flux in a porous medium, J i,p , can thus be written as follows: In one dimension, for the sake of simplicity, equation 6 becomes As an additional simplifying condition, the value of the solvent density is often considered constant and equal to 1,000 kg solvent m −3 solution . It follows In the absence of an external electric field, there is no electrical current and so The combination of equations 7 and 9 provides an expression for the gradient of the diffusion potential Consequently, it is possible to express the Nernst-Planck equation with known parameters only, that is, concentrations, diffusion coefficients, and activity coefficients: The Nernst-Planck equation for the diffusion of individual charged species in a porous medium contains thus two main contributions: 1. a contribution related to the gradient of activity, − D i;e c i ∂ln ciγ i ð Þ ∂x , 2. and a contribution related to the diffusion potential þz i D i;e c i , which arises from the different mobilities of the diffusing species and the zero electrical current condition.
The contribution related to the gradient of activity can be itself split in two contributions: 1. a contribution related to the gradient of concentration − D i;e ∂ci ∂x that corresponds to the Fickian diffusion contribution 2. and a contribution related to the gradient of activity coefficient − D i;e c i ∂ln γ i ∂x . If the diffusive transport processes take place in the presence of a spatially homogeneous background electrolyte composition, the contribution of the activity coefficient gradient can be omitted and equation 11 is simplified to In addition, if all species have the same diffusion coefficient D e , equation 12 is simplified into 10.1029/2019WR026832

Water Resources Research
Because of the electroneutrality condition in solution ∑ j z j ∂cj ∂x is equal to 0, and equation 13 reduces then to the Fickian diffusion equation:

Problem
In the finite difference/volume numerical resolution scheme that is common to most of the reactive transport codes , the properties of the media, that is, porosity, tortuosity, and local concentrations, are defined at the center for each grid cell and apply to the whole of each grid cell. The flux terms, in contrast, have to be evaluated at the interface between two cells. The activity or concentration gradients between two adjacent cells can be evaluated directly for this purpose. However, equations 11, 12, and 14 contain several terms that must be averaged over two adjacent cells. After discretization, with consideration of activity gradients equation 11 becomes For the case where activity coefficient gradients are not considered, equations 12 and 14 become, respectively, where X represents an average value of the parameter X at the interface between two grid cells. Reminding that ABC, the average of A × B × C, is not equal to A × B × C, the product of the average values, in general, the discretization method on a grid makes it necessary to define proper averaging methods for the mean values present in equations 15, 16, and 17.
Most reactive transport codes handle only Fickian diffusion (equation 17), but some can handle the Nernst-Planck equation which includes the diffusion potential term (equations 15 and 16; e.g., Flotran, Crunchflow, MIN3P, and PHREEQC) . Among them, only PHREEQC resolves the dependence of the flux to the activity coefficient gradient (Appelo, 2017;Appelo et al., 2010;Appelo & Wersin, 2007). In the Fickian approximation, it is only necessary to define a proper evaluation of D e . Otherwise it is necessary to define the averaging method for D i;e c i and D i;e . In the following, rigorous averaging methods are derived for all of these terms, and the influence of the averaging methods on the computed diffusive flux is investigated.

Fickian Approximation and Average Value of D e at Interface
In the case where Fick's diffusion equation applies, the flux J x i;p;1→2 from Grid Cell 1 to Grid Cell 2 can be written as follows: where subscripts 1 and 2 indicate that the values are referred to Cell 1 and Cell 2, respectively. Δx 1 and Δx 2 are the lengths of Grid Cells 1 and 2, respectively. It is also possible to define J x i;p;1→int and J x i;p;int→2 , the flux from the center of Cell 1 to the interface and from the interface to the center of Cell 2, where the subscript "int" indicates that the values are referred to the interface between the two cells. The properties within each cell are homogeneous, and it follows Under stationary conditions, and thus, At steady state, the value of the effective diffusion coefficient at the interface D i;e is thus the weighted harmonic mean of D i,e,2 and D i,e,1 .

Average Concentration to be Used in the Nernst-Planck Equation at Interface
Gimmi and Alt-Epping (2018) explored this problem in the specific case of a Donnan membrane system in which a reservoir of electroneutral solution (subscript 1) was considered to be at equilibrium with another reservoir (subscript 2) that contained fixed charges. The solution in Reservoir 2 was not electroneutral, and its charge compensated the fixed charges. The fixed charges were simulated using immobile species (D e = 0), and the system was modeled with the Nernst-Planck equation. The system was considered to be at equilibrium when the net ion fluxes (including diffusion and electrochemical migration) were equal to 0 for each of the species. Because of the presence of the fixed charges in the Reservoir 2, solute species concentrations were not equal in Reservoirs 1 and 2 at equilibrium (zero flux condition). In these conditions, they were able to show analytically and numerically that the average concentration at the interface, c i,int , should be calculated for all mobile species according to It is possible to generalize this result to any diffusion problem in transient nonequilibrium conditions. The activity gradient terms in the Nernst-Planck equation can be expanded into and the concentration gradient term must respect the following mathematical equality: 10.1029/2019WR026832

Water Resources Research
TOURNASSAT ET AL.
It follows after discretization on a grid and thus, between two Cells 1 and 2 In a medium with spatially homogeneous properties (constant D i,e value) or generally when the average D e is independent of the concentrations c i (as is the case for equation 23 as long as the local D i,e are independent of the pore water chemistry), equation 28 becomes Equation 29 is identical to equation 24, but it was obtained for a more general case, that is, without requiring equilibrium or steady-state conditions and without the presence of immobile solute species. One must note that the terms related to the activity coefficient gradient cancel in equation 25; thus, the condition of the absence of activity coefficient gradients is not necessary to apply in equation 28 or 29.
This result shows that the simplifications made from equations 11 to 12 with the equality ∂lny ∂x ¼ 1 y ∂y ∂x might result in reduced accuracy of the results obtained after spatial discretization on a grid. In spatially heterogeneous media, and without any assumptions about equilibrium or steady-state conditions, equation 28

Evaluation of Alternative Averaging Methods on the Computation of Diffusive Fluxes
Reactive transport codes use different types of averaging methods to evaluate the diffusion parameters at interfaces between cells (Tournassat & Steefel, 2019a), and the influence of averaging schemes that are different from the correct one, which is given by equation 31 combined with equation 23, should then be evaluated. Two simple model systems were set up to illustrate these differences.
The first system was made by two reservoirs separated by a membrane. Na + and Cl − concentrations were set to 0.1 mol L −1 in the left reservoir and in the membrane, while the right reservoir contained a solution of 1

10.1029/2019WR026832
Water Resources Research mol L −1 Cl − , 2 mol L −1 Na + , and 1 mol L −1 of a large monovalent anionic molecule for which the membrane was impermeable. To this end, a tortuosity factor of 0 was specifically assigned to this species in the membrane. Consequently, all species were able to diffuse through the membrane except the large anionic molecule. The tortuosity factor of the reservoirs and membrane were set otherwise to 1 for all species. Self-diffusion coefficients (D 0 ) were set to 1.3 × 10 −9 , 2.1 × 10 −9 , and 10 −9 m 2 s −1 for Na + , Cl − , and the large anionic species. The length of the two reservoirs (porosity of 1) was set to 5 mm, and the thickness of the membrane (porosity of 0.1) was set to 200 μm (Figure 1). Each of the reservoir domains was discretized into 25 grid cells. The second system differed from the first one by the absence of the membrane between the two reservoirs, by the size of the grid cells in the left reservoir (100 μm for a total length of 2.5 mm), and by the presence of different tortuosity factors for the different species in the two reservoirs: 0.5 for all species in the left reservoir and 1, 0.7, and 0.2 for Na + , Cl − , and the large anionic species, respectively, in the right reservoir ( Figure 2). The charge of the large anionic molecule was also set to −2, and its concentration was decreased to 0.5 mol L −1 . Three different averaging methods were tested: the reference method given by equation 31 combined with equation 23 and two alternative methods described in Table 1. The alternative Method 1 lumped together the effective diffusion coefficient and the concentration before harmonic averaging at the interface, while the alternative Method 2 computed the harmonic average of D i,e and multiplied it with the weighted arithmetic average concentration at the interface. The diffusion calculations were run using the code 3Diff with an explicit, forward in time and central in space, numerical resolution scheme. This code and its resolution scheme have been benchmarked successfully with CrunchClay and PHREEQC using the arithmetic average method (alternative Method 2) (Tournassat & Steefel, 2019a).
System 1 is representative of a semipermeable membrane system for which a Donnan equilibrium is expected after equilibration. Indeed, the reference model predicted the correct concentrations in the right reservoir corresponding to the Donnan equilibrium (Figure 1, left), a result that was consistent with previous findings of Gimmi and Alt-Epping (2018), who showed the importance of using a logarithmic average for the

10.1029/2019WR026832
Water Resources Research computation of the concentration at the interface between two cells when solving the Nernst-Planck equation in the presence of immobile species. The alternative Method 2 also made it possible to predict the correct concentration but only far from the membrane-reservoir interfaces. Next to this interface, charge balance problems occurred, and electroneutrality was not achieved on both sides of the membrane. This problem illustrates the need to compute correctly the average concentrations in the interfacial terms of the Nernst-Planck equation. The alternative Method 1 resulted in large deviations from electroneutrality, which ultimately led to large concentration oscillations in the numerical solution of the transport equation (Figure 1, right). In System 2, which is very heterogeneous, the reference and

Water Resources Research
the alternative Method 2 led to similar results, while the alternative Method 1 resulted in large concentration oscillations after~300 s of simulated time (Figure 2). The alternative Method 2 (arithmetic averaging) is the method commonly used in reactive transport modeling codes. Our simple intercomparison exercise pointed out the adequacy of the arithmetic averaging method for problems, in which membrane behavior and/or large electrolyte concentration gradients are not present.

Conclusions
In the present study, the proper numerical method was defined to average the concentrations of dissolved species and the porous media properties at the interface between two grid cells in order to solve the Nernst-Planck equation with a finite difference/volume method. The computation of the weighted arithmetic average (alternative Method 2) has been historically the averaging procedure that is used in most reactive transport codes that can solve the Nernst-Planck equations. Our results emphasize the necessity to change this averaging method to one based on a logarithmic-differential average, that is, the reference method demonstrated in the present study and proposed previously by Gimmi and Alt-Epping (2018). The resulting improvement in the numerical stability and in the accuracy of concentration prediction is especially necessary to model semipermeable membrane properties such as those used in water treatment processes.