One-and Two-Equation Models to Simulate Ion Transport in Charged Porous Electrodes

Energy storage in porous capacitor materials, capacitive deionization (CDI) for water desalination, capacitive energy generation, geophysical applications, and removal of heavy ions from wastewater streams are some examples of processes where understanding of ionic transport processes in charged porous media is very important. In this work, oneand two-equation models are derived to simulate ionic transport processes in heterogeneous porous media comprising two different pore sizes. It is based on a theory for capacitive charging by ideally polarizable porous electrodes without Faradaic reactions or specific adsorption of ions. A two-step volume averaging technique is used to derive the averaged transport equations for multi-ionic systems without any further assumptions, such as thin electrical double layers or Donnan equilibrium. A comparison between both models is presented. The effective transport parameters for isotropic porous media are calculated by solving the corresponding closure problems. An approximate analytical procedure is proposed to solve the closure problems. Numerical and theoretical calculations show that the approximate analytical procedure yields adequate solutions. A theoretical analysis shows that the value of interphase pseudo-transport coefficients determines which model to use.


Introduction
In systems comprising different phases, conservation equations for the properties under study are normally derived for each separate phase.This approach leads to the formulation of a multi-equation model.However, it is desirable to mathematically represent the system with a single set of conservation equations.This desirable approach is called one-equation modelling (Whitaker [1]).A one-equation model can be extracted from a multi-equation model based on the principle of local equilibrium (Whitaker,[2,3]; Kaviani, [4]; Quintard and Whitaker, [5,6]; del Rio and Whitaker, [7,8]; among others).
The goal of this work is to derive one-and two-equation models for capacitive processes in ideally polarizable porous electrodes with a binary pore size distribution, e.g., macropores with pore diameters bigger than 100 nm and micropores with pore diameters of the order of 2 nm.
It is generally accepted that porous electrodes for electrochemical applications should have a hierarchical pore size distribution including macropores (d p ≥ 50 nm) as ion reservoirs, mesopores (2 nm ≤ d p ≤ 50 nm) to facilitate ion transport, and micropores (d p ≤ 2 nm) for increased ion storage (Tsouris et al. [9], among others).
Johnson and Newman [10] presented a porous electrode model that considered, for the first time, many of the characteristic properties of capacitive ionic transport.The authors considered adsorption of ions in electrical double layers (EDLs) and concluded that in the absence of concentration gradients, the system behaves as a network of resistances and capacitances.The authors also suggested that a porous system could show different behaviors at different length scales; for example, transient behavior at the microscale and steady-state behavior at the macroscale.Simulating the capacitive deionization (CDI) process, Biesheuvel and Bazant [11] presented a model for capacitive charging and desalination by ideally polarizable porous electrodes.This model excludes effects of Faradaic reactions or specific adsorption of ions.The authors discussed the theory for the case of a dilute, binary electrolyte using the Gouy-Chapman-Stern (GCS) model of the EDL.Gabitto and Tsouris [12] used the physical model presented by Biesheuvel and Bazant [11] to study the capacitive deionization process in homogeneous porous electrodes by means of a volume averaging technique [1].The authors derived equations identical to the ones reported by Biesheuvel and Bazant [11] plus an alternative formulation of the same problem.Gabitto and Tsouris [12] also calculated the value of the transport coefficients in isotropic porous media by solving the appropriate closure problems.
The model of Biesheuvel and Bazant [11] has been extended to include Faradaic reactions and a dual-porosity (macropores and micropores) approach [13] which considers that the electrodes are composed of solid particles that are porous themselves.Biesheuvel et al. [13] used, for the first time, a modified Donnan (mD) approach for the micropores, valid for strongly overlapped double layers.Biesheuvel et al. [14] also presented a porous electrode theory for ionic mixtures, including Faradaic reactions.Recently, Biesheuvel et al. [15] improved the mD model by making the ionic attraction term dependent on total ion concentration in the carbon pores.The authors reported that the new mD model significantly improves predictions of the influence of salt concentration on the CDI performance.Gabitto and Tsouris [16] presented a model to simulate the CDI process in heterogeneous porous media comprising two different pore sizes.A two-step volume averaging technique was used to derive the averaged transport equations in the limit of thin electrical double layers.The authors derived a one-equation model based on the principle of local equilibrium (Whitaker,[2,3]; Quintard and Whitaker, [5,6]; del Rio and Whitaker, [7,8]; among others).Gabitto and Tsouris [16] derived the constraints determining the range of applications of the one-equation model and reported 'approximate' two-equation model equations.Gabitto and Tsouris [16] used the thin EDL assumption inside the micropores and electroneutrality in the solution occupying the bulk of the pores.This assumption imposes severe restrictions on the pore size and it can be considered unrealistic in many circumstances.Recently, Gabitto and Tsouris [17] derived the 'exact' version of the two-equation model for this problem.The authors found that the derived new equations coincide with the original ones proposed by Gabitto and Tsouris [16].
The state of the art in this field is given by the phenomenological mD model presented by Biesheuvel and co-workers.However, this model represents a limiting condition as it is only strictly applicable to micropores and it precludes transport inside the pores.In this work, we propose a general model which could be applied to pores of different sizes and could allow for transport inside the pores.We will derive one-equation and two-equation model equations describing ionic concentrations and electrostatic potential transport.We will also calculate effective transport coefficients for isotropic porous media assuming that there is not convective transport inside the macropores.Finally, we will theoretically compare the one-and two-equation models.

Model Description
The physical model developed in this work is based on the one proposed by Biesheuvel et al. [13] and used by Gabitto and Tsouris [16] to develop a one-equation model formulation for salt concentration and electrostatic potential in the limit of thin EDLs.The porous electrode is divided into pore space filled with quasi-neutral electrolyte and a solid matrix.The solid matrix is itself porous, and these smaller pores are also filled with electrolyte.This representation is schematically shown in Figure 1.The solid electrode presents two length scales.The macroscale is represented by a representative elementary volume (REV M ) of radius R M , comprising big pores (γ-phase) and porous particles (κ-phase).The porous particles are represented by an REV of radius R m comprising small pores (α-phase) and solid (β-phase).The 'small' and 'big' pores are related by the constraint, l α << l γ .Based on this constraint, we can say that ion adsorption will occur mostly in the small pores inside the solid particles.This assumption is based on the fact that the external area of the particles will be negligible compared to the internal area.pores (α-phase) and solid (β-phase).The 'small' and 'big' pores are related by the constraint, lα << lγ.
Based on this constraint, we can say that ion adsorption will occur mostly in the small pores inside the solid particles.This assumption is based on the fact that the external area of the particles will be negligible compared to the internal area.

Electrostatic Potential
Starting from the point formulation for the Poisson equation, Gabitto and Tsouris [18] derived the relevant equations corresponding to this case: where zi is the ionic charge, av is the specific surface area (Aαβ/Vm), Vm is the volume of the small scale REV, Aαβ is the interface α-β area, VT is the thermal voltage (R T/F), R is the gas constant, F is the Faraday constant, T is the temperature, εα is the volume fraction , respectively, given by: is the area average of the solid charge (σ) given by: The effective permittivity tensor ( eff  ) can be expressed as:

Electrostatic Potential
Starting from the point formulation for the Poisson equation, Gabitto and Tsouris [18] derived the relevant equations corresponding to this case: where z i is the ionic charge, a v is the specific surface area (A αβ /V m ), V m is the volume of the small scale REV, A αβ is the interface α-β area, V T is the thermal voltage (R T/F), R is the gas constant, F is the Faraday constant, T is the temperature, ε α is the volume fraction (ε α = V α /V m ), φ α , and c i α are the intrinsic averages for the dimensionless potential (ψ = φ/V T ) and the ionic concentrations (c i ), respectively, given by: σ αβ is the area average of the solid charge (σ) given by: The effective permittivity tensor (η ) can be expressed as: Here, η is the electrolyte permittivity, n αβ is the normal vector pointing from the α-phase into the β-phase, and g is the closure variable used to define the potential deviations as: A nomenclature is included in the Supplementary Materials.

Ionic Concentrations Equation
We again follow the work by Gabitto and Tsouris [18] deriving the following volume averaged equation: Here, the effective diffusivity (D i,e f f ) and effective mobility (U i,e f f ) tensors are given by, The f vector field required for calculation of the effective diffusivity tensor is calculated by solving the boundary value Problem 1, see Appendix A.
Gabitto and Tsouris [18] showed that the value of the effective diffusivity (D i,e f f ) is a function of the geometry and the charge at the electrolyte-solid interface, while the effective mobility (U i,e f f ) is only a function of the geometry.Therefore, the effective diffusivity and the effective mobility are not equal as is the case in the point equations.
The value of the salt concentration in the β-phase is equal to zero; therefore, there is no need to calculate a one-equation model.The assumption of the solid phase being an ideal conductor was used; therefore, the electrostatic potential in the β-phase is constant and an average potential in the system comprising the α-and β-phases can be calculated using: (10) According to the definition given by Equation (10) we find that: In this section, we will develop the spatially smoothed equations associated with the large scale REV M shown in Figure 1.In order to simplify the nomenclature in Equation (1), we will replace c αi α with c ki and ψ α with ψ κ using Equation (11).Both c ki and ψ κ represent the salt concentration and potential in the porous particles.
The boundary conditions in the interphase between the porous particles, κ-phase, and the γ-phase have received extensive attention in literature (Quintard Whitaker,[19]; Ulson de Souza and Whitaker, [20]; Valdès-Parada et al. [21]; Gabitto and Tsouris, [16], among others).Following this approach, we can obtain the following equations: The κ-Region We apply the phase average in the κ-phase to Equation ( 12) plus the spatial averaging theorem.After algebraic manipulations, we obtain: In the derivation of Equation ( 16) we have used Gray's decomposition [22] for the potential: The γ-Region We apply the phase average in the γ-phase plus algebraic manipulations to Equation ( 15) to obtain: In the derivation of Equation ( 18), we have used Gray's decomposition [22] for the potential: The final formulation either of the one-equation model or the two-equation model requires constitutive equations for the electrostatic potential ( ψ i ).This issue will be addressed in the following sections.

One-Equation Model for Electrostatic Potential
Equations ( 16) and ( 18) are the starting point for the one-and two-equation model derivations for the electrostatic potential.In order to derive a one-equation model, we assume that the mass and potential equilibrium principle holds.This assumption leads to the following identities: Colloids Interfaces 2018, 2, 4 6 of 20 We add Equations ( 16) and ( 18) and use Equations ( 20) and ( 21) to obtain: Equation ( 22) is the one-equation model for the electrostatic potential.Here, σ αβ κ is the intrinsic average of the charge average on the interfacial area α-β within the macroscale REV M .

Ionic Concentration Equations
In this section, we will develop the spatially smoothed ionic concentration equations associated with the large scale REV M , which is shown in Figure 1.In order to simplify the nomenclature in Equation ( 1), we will replace c αi α with c ki , which represents the species i concentration in the electrolyte located inside the porous particles.The boundary conditions in the interphase between the porous particles, κ-phase, and the electrolyte, γ-phase, have been used extensively ( [16,[19][20][21], among others).Following this approach results in the following equations: The κ-Region We start from Equation ( 27) and, after complex algebraic manipulations, we obtain: In the derivation of Equation (31), we used Gray's decomposition [22] for the ionic concentrations in the κ-phase: A key step in the derivation is the calculation of the average of the product of ionic concentration and the gradient of the potential ( c κi ∇ψ κ ).Algebraic manipulations lead to: Here, c κi ∇ ψ κ is the ionic concentration dispersion produced by the electric field.In order to close Equation (31), we need constitutive equations for the ionic concentrations ( c κi ) and potential ( ψ κ ) deviations.This issue will be addressed in the following sections.

The γ-Region
Application of the phase average in the γ-phase plus algebraic manipulations to Equation (30) leads to: In the derivation of Equation (34), we have used Gray's decomposition for the ionic concentrations in the γ-phase: The average of the hydrodynamic dispersion term ( v s c γi ) has been calculated following Carbonell and Whitaker [24]: The average of ionic concentration dispersion in the γ-phase is given by:

One-Equation Model for Ionic Concentrations
We introduce the assumption of mass and potential equilibrium, Equations ( 20) and ( 21), into Equations ( 31) and (34), leading to: Equation ( 38) needs constitutive expressions for the potential and ionic concentration deviations.These expressions will be derived in the next section.

Closure for Ionic Concentrations One-Equation Model
Following the work of Gabitto and Tsouris [18] we propose: The closure variables ( f γi , f κi ) are calculated by solving boundary value Problem 3, see Appendix A.
Introduction of Equations ( 39) and (40) into Equation (38) leads to the final formulation for the ionic concentration using the one-equation model: Here, the effective global diffusivity (D * i,di f f ), the effective dispersion (D * i,disp ), and the effective mobility (U * i ) are given by: A total effective dispersion tensor can be obtained by combining the global diffusivity tensor (D * i,e f f ) with the effective dispersion tensor (D * disp ), Equations ( 22) and (45) form the foundation of a one-equation model to simulate ionic transport inside porous carbon electrodes.The validity of the one-equation model approach rests upon constraints similar to the ones presented by Gabitto and Tsouris [16].These constraints are based upon the idea of achieving thermodynamic equilibrium and, therefore, high transport rates between the two phases at the macroscale.This issue was discussed further by Gabitto and Tsouris [17].

Two-Equation Model for the Electrostatic Potential
Following Quintard and Whitaker [5], we propose the following expressions for the electrostatic potential deviations in the κ-and γ-phases: The parameters g i , h i , and s i are the closure variables for the electrostatic potential [5].The terms, ∇ c κ κ , ∇ c γ γ , and ( c κ κ − c γ γ ) are considered sources in the boundary value problems used to calculate the six closure variables.In two-equation models, the sources are normally considered constant [1].Quintard and Whitaker [5] conducted a detailed analysis of the closure problem and found that, due to the presence of higher order terms, one cannot a priori consider these source terms constant.However, we will assume that either the source terms are constant or the variation of these terms does not appreciably change the results.
(a) The κ-Region Introducing Equations ( 46) and (47) into Equation (31) leads to: Equation ( 48) can be rewritten as: Here, a vM is the effective area in the macroscale, σ αβ κ is the volume average of the interfacial charge σ αβ , and the transport tensors (η κκ and η κγ ) are given by: The interfacial transport coefficient k is defined as: Introducing Equations ( 46) and (47) into Equation (34) leads to: Colloids Interfaces 2018, 2, 4 10 of 20 The transport coefficients (η γκ and η γγ ) are calculated by: The closure parameters appearing in Equation ( 53) are calculated solving the appropriate closure problems, see appendices.

Two-Equation Model for the Ionic Concentrations
Following Quintard and Whitaker [5], we propose the following expressions for deviations in the ionic concentrations in the κ-and γ-phases: The parameters v γi , w γi , v γi , w γi , t κi, and t γi are the closure variables [5].

(a) The κ-Region
We start the derivation from Equation (31).Replacing c κi by Equation ( 27), after algebraic manipulations, leads to: The transport parameters are calculated from: The interfacial transport coefficients k i are given by: We start the derivation from Equation (34).Replacing c κi and c γi by Equations ( 54) and (55), after algebraic manipulations, leads to: The transport parameters are calculated from: (62) 2.3.9.Closure for the Ionic Concentrations Two-Equation Model Calculation of the closure variables (v κi , v γi , w κi , and w γi ) is required to evaluate the transport parameters (D ij and U ij ).This is achieved by solving boundary value Problems 4, 5, and 6, see Appendix A.
The relationship among the one-and two-equation transport parameters is given by: In the case where the value of the macroscale velocity (v s ) is zero, it can be proved, following Quintard and Whitaker [5], that the three problems for calculation of the ionic concentration closure variables will be equal to the ones derived for potential.The values of the transport parameters are given by equations similar to the ones reported by Gabitto and Tsouris [17] for the limiting case of thin EDLs.This symmetry is broken by the presence of the convective term.The transport parameters (D i,κκ , D i,γκ ) related to the closure parameters (v γi , v κi ) are still calculated by the same analytical expressions (see results section below), while the transport parameters (D i,γγ , D i,κγ ), which depend upon the (w γi , w κi ) closure parameters, can only be calculated numerically.Quintard and Whitaker [6] provide examples of these calculations.

Introduction
In this section, we will estimate the values of the macroscopic effective diffusivity tensors (D κκ , D κγ , D γκ , and D γγ ) plus the pseudo-mass-transport coefficients k and k i .In order to simplify the math, we will use the assumptions of no convective flow inside the macropores (ν s = 0) and that all ionic bulk diffusivities are equal (D i = constant).These assumptions allow us to drop the i-subscripts in the effective diffusivity tensors.We have to solve closure Problems 3, 4, and 5 for the salt concentration, and then use the appropriate integral equations to calculate the values of the macroscopic diffusivity tensors.In order to do so, we will solve the closure problems in isotropic spatially periodic unit cells.The use of these cells will reduce the problem to calculation of only one component of the effective diffusivity tensors (D ij,xx = i•D ij •i).

Calculation of One-Equation Model Transport Parameters
Gabitto and Tsouris ( [16,18]) solved numerically the closure problem in Chang's unit cells (Chang [25,26]), see Figure 2. A spatially periodic porous medium, where the particles are either cylinders of infinite length or spheres, can be represented by a square unit cell containing a circular particle in the case of infinite-length cylinders, or a cubic cell containing a body-centered spherical particle.Chang [25,26] proved that in the approximate (circular) unit cell shown in Figure 2, the periodic boundary condition on the cell boundary (r = 1) can be replaced by a zero Dirichlet boundary condition.Ochoa [27] examined the issue in detail numerically and concluded that reasonably good agreement exists between the numerical solutions in the real unit cell and the analytical ones in Chang's unit cells.An example of the calculations to solve the closure problems is shown in Appendix B.
Gabitto and Tsouris [17] imposed a Dirichlet boundary condition at r = 0 instead of the normal symmetry condition.Obviously, this choice leads to an approximate solution; however, the authors argued that, in order to calculate the value of the transport coefficients, only the values of the closure parameters on the A κγ are relevant.If the approximate solution gives a good approximation of the real solution at that location, the fact that it may be inaccurate somewhere else is not important.Similar arguments have been made concerning the use of the Dirichlet boundary conditions at r = 1.Gabitto and Tsouris [17] showed that this approximation leads to acceptable values of the global effective conductivity tensor.We used this concept to solve analytically the boundary value Problem 3 for the one-equation model following the procedure reported by Ochoa-Tapia et al. [28].The solution is presented in Appendix B.
(Chang [25,26]), see Figure 2. A spatially periodic porous medium, where the particles are either cylinders of infinite length or spheres, can be represented by a square unit cell containing a circular particle in the case of infinite-length cylinders, or a cubic cell containing a body-centered spherical particle.Chang [25,26] proved that in the approximate (circular) unit cell shown in Figure 2, the periodic boundary condition on the cell boundary (r = 1) can be replaced by a zero Dirichlet boundary condition.Ochoa [27] examined the issue in detail numerically and concluded that reasonably good agreement exists between the numerical solutions in the real unit cell and the analytical ones in Chang's unit cells.An example of the calculations to solve the closure problems is shown in Appendix B.
Gabitto and Tsouris [17] imposed a Dirichlet boundary condition at r = 0 instead of the normal symmetry condition.Obviously, this choice leads to an approximate solution; however, the authors argued that, in order to calculate the value of the transport coefficients, only the values of the closure parameters on the Aκγ are relevant.If the approximate solution gives a good approximation of the real solution at that location, the fact that it may be inaccurate somewhere else is not important.Similar arguments have been made concerning the use of the Dirichlet boundary conditions at r = 1.Gabitto and Tsouris [17] showed that this approximation leads to acceptable values of the global effective conductivity tensor.We used this concept to solve analytically the boundary value Problem 3 for the one-equation model following the procedure reported by Ochoa-Tapia et al. [28].The solution is presented in Appendix B.   After complex algebraic manipulations, we obtain: Here, δ is the ratio of the solid and liquid phase diffusivities (δ = ε α D eff,XX /D γ ).
In Equation (65) we have assumed that all ionic diffusivities are identical.This assumption allows us to drop the i-subscript.Equation ( 65) is identical to Maxwell's equation [29] for heat transfer in infinite cylinders.This result is similar to the numerical calculations reported by Gabitto and Tsouris [18] and the analytical solution found by Gabitto and Tsouris [17] using the thin layer approximation for the same system.

Calculation of Two-Equation Model Transport Parameters
In this case we start our calculations by inserting the solution proposed in Equation (B7) into the boundary value Problem 4. The procedure followed by Gabitto and Tsouris [17], taken from Ochoa-Tapia et al. [28] leads to the following equations for the effective diffusivity transport coefficients, see Appendix B.
Here, D κκ,XX = i•D κκ •i is the only non-zero term in the effective diffusivity tensor for isotropic porous media.
In Equation (70) we have defined a Sherwood dimensionless number using the unit cell characteristic length r 2 and the porous media effective area (a vM ) as geometric parameters.

Comparison of One-and Two-Equation Model Transport Parameters
Equations ( 63) and (64) give the relationship among the one-equation and two-equation transport parameters; for example, we can derive an algebraic expression for D * by using the equations for the D ij transport coefficients derived in the previous section: Here, D * is the global diffusivity tensor component calculated using the two-equation model formulation.Simple inspection shows that Equation ( 71) is identical to Equation (65), which was derived using the one-equation model formulation.This result confirms that the two-equation model reduces to the one-equation model when the local equilibrium conditions, Equations ( 20) and ( 21), hold.

Transport Parameters
The algebraic Equations (67) to (69) allow calculation of the effective diffusivity transport coefficients from the values of the γ-phase volume fraction (ε γ ) and the ratio of the effective diffusivities in both phases (δ).In Figures 3-5, we show typical variation in the effective diffusivities with the porous medium void fraction.In all these Figures, we plot the value of the isotropic effective diffusivity versus the ratio of the porous solid and liquid phases effective diffusivities (δ = ε α D e f f ,XX /D γ ) using the liquid phase void fraction (ε γ ) as a fixed parameter.A constant value for the microscale void fraction (ε α ) equal to 0.5 was used when needed.In CDI processes the δ ratio has values smaller than 1; however, in these figures we plot δ values from 0.01 to 100 for comparison with the heat transfer results from Quintard and Whitaker [5] who estimated numerically the values of the two-equation model effective heat conductivities in cell (a) depicted in Figure 2 for the case of infinite-length cylinders.In all these figures the lines represent the results calculated in this work.
In Figure 3, we plot the value of the diagonal porous solid diffusivity (D κκ ).It can be seen that the value of the effective diffusivity increases as the value of the δ parameter increases.The effective diffusivity also decreases as the liquid void fraction (ε γ ) increases.The values of the transport parameter calculated in this work agree very well with the numerical calculations of Quintard and Whitaker [5] for ε γ values higher than or equal to 0.5.For lower liquid phase void fractions, our calculated values are smaller than the ones reported in [5].
In Figure 4, we plot the value of the cross-term effective diffusivities, D κγ .The transport coefficient increases as the δ ratio increases.The cross-term effective diffusivities (D κγ and D γκ ) are related by Equation (63); therefore, while the κ-phase coefficient increases, the γ-phase coefficient decreases as ε γ increases.There is very good agreement between the values calculated in this work and the numerical calculations of Quintard and Whitaker [5].calculated values are smaller than the ones reported in [5].
In Figure 4, we plot the value of the cross-term effective diffusivities, Dκγ.The transport coefficient increases as the δ ratio increases.The cross-term effective diffusivities (Dκγ and Dγκ) are related by Equation (63); therefore, while the κ-phase coefficient increases, the γ-phase coefficient decreases as εγ increases.There is very good agreement between the values calculated in this work and the numerical calculations of Quintard and Whitaker [5].In Figure 5, we plot the variation in the γ-phase diagonal effective diffusivity (Dγγ) with the δ ratio.It is shown in the Figure that the γ-phase diagonal effective diffusivity increases as the liquid void fraction increases.The results depicted in Figure 5 show that the values of the effective diffusivity slightly increase as the δ parameter increases.We can also observe that, as the void fraction increases, the value of the transport parameter tends to the bulk solution diffusivity value (Dγ).In all cases there is remarkable agreement with the numerical results of Quintard and Whitaker [5].In Figure 5, we plot the variation in the γ-phase diagonal effective diffusivity (D γγ ) with the δ ratio.It is shown in the Figure that the γ-phase diagonal effective diffusivity increases as the liquid void fraction increases.The results depicted in Figure 5 show that the values of the effective diffusivity slightly increase as the δ parameter increases.We can also observe that, as the void fraction increases, the value of the transport parameter tends to the bulk solution diffusivity value (D γ ).In all cases there is remarkable agreement with the numerical results of Quintard and Whitaker [5].
ratio.It is shown in the Figure that the γ-phase diagonal effective diffusivity increases as the liquid void fraction increases.The results depicted in Figure 5 show that the values of the effective diffusivity slightly increase as the δ parameter increases.We can also observe that, as the void fraction increases, the value of the transport parameter tends to the bulk solution diffusivity value (Dγ).In all cases there is remarkable agreement with the numerical results of Quintard and Whitaker [5].In the comparison of the one-equation and two-equation models, the evaluation of the interphase transport parameters (avM k and avM ki) plays a crucial role.In Figure 6 In the comparison of the one-equation and two-equation models, the evaluation of the interphase transport parameters (a vM k and a vM k i ) plays a crucial role.In Figure 6, we study the variation in the Sh interphase transport parameter with void fraction and the δ ratio.The results in Figure 6 show that the interphase mass transfer coefficient increases as the δ ratio increases.The mass transfer coefficient values also decrease as the void fraction values increase.These results are in agreement with the ones reported by Quintard and Whitaker [5]; however, our calculated results are smaller than the ones reported by these authors.In the Figure , we corrected the values predicted by Equation (70) by a geometrical factor (π) to allow for the different unit cells, (a) and (b) in Figure 2, used in [5] and in our calculations.The values of the dimensionless mass transfer coefficients (Sh and Sh i ), calculated using values of the δ ratio commonly used in the CDI process, are O{1-10}.

Comparison of One-Equation and Two-Equation Models
Equations ( 56) and ( 60) can be made dimensionless as:

Comparison of One-Equation and Two-Equation Models
Equations ( 56) and (60) can be made dimensionless as: (74) Inspection of Equations ( 72) to (74) shows that the one-equation model will be acceptable for high values of the interphase mass transfer coefficients Sh kM and Sh kiM .The following constraint can be derived considering that the microscale Sh and Sh i values are O{1-10}: The constraint given by Equation ( 75) is consistent with the constraint derived by Gabitto and Tsouris [16] using the thin EDL approximation.Furthermore, it approximates the mixed length scale l k by the unit cell characteristic length (r 2 ) and the transport process characteristic length (L) by the porous electrode thickness (L e ).Gabitto and Tsouris [17] numerically verified that, as the mass transport coefficients increase, the solution calculated using the two-equation model tends to the solution predicted by the one-equation model.
An interesting conclusion can be drawn from the values of the Sh kM number estimated from Equation (75).At the microscale, the Sh values are O{1-10}, while at the macroscale, the Sh kM values are O{>10 3 }.Therefore, we can conclude that, at the microscale, the interphase mass transport process is transient, while at the macroscale, it is at steady-state.This conclusion agrees well with the observations made by Johnson and Newman [10] in their classic article on this subject.

Conclusions
The volume average method has been used to derive the average one-and two-equation model equations describing salt capture in dual-porosity porous electrodes by electrosorption.We have derived the complete form of the volume averaged equations starting from the point equations and the appropriate boundary conditions.The assumption of thin EDLs has been dropped.In the porous solid-particle phase, the concentration and electrostatic potential deviation fields were found to be functions of the corresponding averaged variables.The numerous effective diffusivities appearing in the transport equations were estimated by solving the appropriate closure problems in Chang's unit cells.The approximate analytical solution predicts Maxwell's [29] solution for infinite cylinders geometry.
We have shown that the one-and two-equation models converge to similar results as the interphase mass transfer coefficients (Sh M and Sh iM ) increase.This conclusion is in agreement with the constraints reported by Gabitto and Tsouris [16] in their work on the one-equation model for

Figure 1 .
Figure 1.Schematic representation of the model.

Figure 1 .
Figure 1.Schematic representation of the model.
calculated by solving boundary value Problem 2 in an appropriate local domain, see Appendix A.

Figure 2 .
Figure 2. Schematic of the spatially periodic (a) and Chang's approximate (b) unit cells.

Figure 2 .
Figure 2. Schematic of the spatially periodic (a) and Chang's approximate (b) unit cells.

Figure 4 .
Figure 4. Variation in the εa Dκγ effective diffusivity with the δ ratio for several void fractions (εγ).

Figure 4 .
Figure 4. Variation in the ε a D κγ effective diffusivity with the δ ratio for several void fractions (ε γ ).

Figure 5 .
Figure 5. Variation in the γ-phase diagonal effective diffusivity with the δ parameter.

Figure 5 .
Figure 5. Variation in the γ-phase diagonal effective diffusivity with the δ parameter.

21 Figure 6
Figure6show that the interphase mass transfer coefficient increases as the δ ratio increases.The mass transfer coefficient values also decrease as the void fraction values increase.These results are in agreement with the ones reported by Quintard and Whitaker[5]; however, our calculated results are smaller than the ones reported by these authors.In the Figure,we corrected the values predicted by Equation (70) by a geometrical factor (π) to allow for the different unit cells, (a) and (b) in Figure2, used in[5] and in our calculations.The values of the dimensionless mass transfer coefficients (Sh and Shi), calculated using values of the δ ratio commonly used in the CDI process, are O{1-10}.

Figure 6 .
Figure 6.Variation in the interphase migration mass transport parameter (Sh) with the δ parameter.

Figure 6 .
Figure 6.Variation in the interphase migration mass transport parameter (Sh) with the δ parameter.