The concepts presented in this section applies to flat electrodes and to porous materials where the pore size is much bigger than the curvature of the surface (big pores). The presence of a charged surface in contact with an electrolyte will produce a re-arrangement of charges in the liquid phase. Ions with the same charge as the solid surface (co-ions) will be repelled while ions of opposite charge (counter-ions) will be attracted to the surface. Molecules of the solvent, water for example, can also be oriented in a particular direction due to the presence of the surface charge. The level of orientation will depend upon the local electric field that is determined by the surface charge or the applied potential. Classical models of the EDL were originally derived for flat solid surfaces (Helmholtz, 1853; Stern, 1924; Bockris et al., 1963; Bockris et al., 1998; Newman and Thomas-Alyea, 2004; IUPAC, 2019). Critical understanding of the structure of the EDL was achieved through the work of Grahame (Grahame, 1947) who experimentally determined the differential capacity of mercury in contact with aqueous solutions of sodium fluoride as a function of concentration. Several different models were proposed to reproduce Grahame’s experimental data (Bockris et al., 1963; Bockris et al., 1998).

Presently, there is general agreement on the model depicted in Figure 1. The region between the solid surface and the inner Helmholtz plane contains a layer of solvent molecules and/or some specifically (chemically) adsorbed ions. According to IUPAC, the inner Helmholtz plane (IHP) is the locus of the electrical centers of specifically adsorbed ions (IUPAC, 2019). There is a layer of counterions located, on average, a distance from the electrode equal to the hydrated radius of the counterion. The outer Helmholtz plane (OHP) is the locus of the electrical centers of these ions (IUPAC, 2019). In the diffuse layer the ions follow a Boltzmann distribution and there is a net charge given by the counterions concentration (Landau and Lifshitz, 1980).

In the previous sections we have established a qualitative picture of ion storage in charged porous materials. In this section we try to review the most common models for ion transport and storage reported in the literature. The different theoretical models to simulate transport phenomena in charged porous media can be classified in two main groups, phenomenological models based on general and/or empirical assumptions, and models based on rigorous mathematical methods of upscaling transport phenomena in charged porous media. Several upscaling techniques can be used to study transport in porous media, multi-scale asymptotic expansions, mathematical homogenization, spatial averaging, moment methods, stochastic-convective approaches, various Eulerian and Lagrangian perturbation schemes, projection operators, renormalization group techniques, space transformational methods, continuous time random walks, etc. (Cushman et al., 2002). In this work we will concentrate on the volume averaging technique as there are several relevant applications published in literature (Whitaker, 1999; del Rìo and Whitaker, 2002a; Schmuck and Bazant, 2014; Gabitto and Tsouris, 2015; Gabitto and Tsouris, 2016; Gabitto and Tsouris, 2017; Gabitto and Tsouris, 2018; Gabitto and Tsouris, 2019).

The article by Johnson and Newman (Johnson and Newman, 1971) presented for the first time a porous electrode model that considered the characteristic properties of capacitive ionic transport. The authors considered adsorption of ions in electrical double layers and concluded that in the absence of concentration gradients the system behaves as a network of resistances and capacitances. The authors also suggested the important concept that a porous system could show different behavior at different length scales; for example, transient behavior at the macroscale and steady-state behavior at the microscale (Johnson and Newman, 1971).

Biesheuvel and Bazant (Biesheuvel and Bazant, 2010) presented a model for capacitive charging and desalination by ideally polarizable porous electrodes. This model excluded effects of Faradaic reactions or specific adsorption of ions. The authors discussed the theory for the case of a dilute, monovalent, binary electrolyte using the Gouy-Chapman-Stern (GCS) model of the EDL. The model is strictly valid in the limit of “thin EDLs” and applies to big pores (macropores and mesopores >10 nm). The authors mentioned in their paper the use of a volume averaging technique, but they did not present derivations that support this claim. Therefore, the model should be considered as based upon qualitative concepts of the volume averaging methodology (Whitaker, 1999). Biesheuvel and Bazant (Biesheuvel and Bazant, 2010) presented the following equations: ∂ C i ∂ t ⏟ A c c u m u l a t i o n = ∇ → · { ∇ → C i + z i C i ∇ → ψ } ⏟ I n − O u t + R i ⏟ G e n e r a t i o n + a v J i , s a l t ⏟ I n − O u t   t o   E D L (10)

Here, z _i is the ionic charge, C i is the ion concentration, R i is the ionic volumetric generation in the bulk of the electrolyte, ψ is the dimensionless electrostatic potential ( ψ = φ   V T ), V _T is the thermal voltage (R T/F), ϕ is the dimensional potential, and J i ,   s a l t is the mass flow of the i-species into the EDL. In the absence of chemical reaction and for a monovalent, binary electrolyte, Eq. 10 can be written as: ∂ C s a l t ∂ t ⏟ A c c u m u l a t i o n = ∇ 2 C s a l t ⏟ I n − O u t + a v   J s a l t ⏟ I n − O u t   t o   E D L (11) Here, J s a l t is the mass flow of the ionic species into the EDL ( J s a l t = J + ,   s a l t + J − ,   s a l t ).

A charge balance for the same electrolyte leads to: 0 A c c u m u l a t i o n = − ∇ → · C s a l t ⏟ I n − O u t ∇ → ψ ) + a v   J c h a r g e ⏟ I n − O u t   t o   E D L (12) Here, I → e = − C s a l t   ∇ → ψ is the ionic current and J c h a r g e is the charge flow into the EDL ( J c h a r g e = J + ,   s a l t − J − ,   s a l t ).

The last terms in Eqs 11, 12 can be replaced by the rate of change of the excess salt concentration in the EDL ( w ) and the rate of change of excess charge density ( σ ), respectively by using the boundary conditions between the bulk of the solution and the EDL. J c h a r g e = ∂ w ∂ t a n d J s a l t = ∂ σ ∂ t (13)

Gabitto and Tsouris (Gabitto and Tsouris, 2015) studied the ionic transport process in homogeneous porous electrodes using a volume averaging method formulation (Whitaker, 1999). The method of volume averaging is a homogenization technique that is used to derive continuum equations for multiphase systems. Equations that are valid in a particular phase can be spatially smoothed to produce equations that are valid everywhere in the domain (Whitaker, 1999). The volume averaging methodology not only allows for the derivation of the continuum equations but provides a way to estimate the transport coefficients appearing in the derived equations by solving the necessary closure problems.

Based upon the same constraints, the equations derived by Gabitto and Tsouris (Gabitto and Tsouris, 2015) were identical to the ones reported by Biesheuvel and Bazant (Biesheuvel and Bazant, 2010). The authors also derived an alternative formulation of the same problem.

The Biesheuvel and Bazant (Biesheuvel and Bazant, 2010) model was extended to include Faradaic reactions and a dual-porosity (macropores and micropores) approach (Biesheuvel et al., 2011; Biesheuvel et al., 2012) which considers that the electrodes are made of porous solid particles. This model has become very popular, and it has been used by many different researchers in many different studies [(Biesheuvel et al., 2010; Biesheuvel and van der Wal, 2010; Biesheuvel et al., 2011; Biesheuvel et al., 2012; Biesheuvel et al., 2014; Porada et al., 2014; Biesheuvel and Dykstra, 2021), among others]. It has also been extended to other research areas as ionic flow in membranes (Biesheuvel et al., 2010; Biesheuvel and van der Wal, 2010) and flow capacitive deionization [(Porada et al., 2014; Nativ et al., 2017), among others]. Suss et al. (2015) reviewed critically the CDI technology discussing these issues. The authors also discussed equilibrium and kinetic modeling aspects of the different CDI processes.

Biesheuvel et al. (2011), Biesheuvel et al. (2012) proposed the following equations: ε M ∂ C i M ∂ t + ε m ∂ C i m ∂ t ⏟ A c c u m u l a t i o n = ∇ → · D i M ∇ → C i M + z i C i M ∇ → ψ ⏟ I n − O u t + R F , i ⏟ G e n e r a t i o n b y   r e d o x   r e a c t i o n   i n   m i c r o p o r e s (14)

Here, D i M is the effective diffusivity of the i-component in the macropores, ε m and ε M are the micro and macro porosities, respectively. Equation 14 has the shortcoming of including two dependent variables, the macro- and micropores concentrations ( C i M and C i m ); therefore, an equation relating both concentrations is needed. Biesheuvel et al. (2011) used the modification of the Donnan approach (mD) given by Eqs 1–5. Biesheuvel et al. (2012) also proposed a quadratic variation for the volumetric Stern layer capacity with either charge or potential drop. The authors reported that the quadratic dependance leads to better fit of their experimental data (Biesheuvel et al., 2012). C i S t = C i , o S t + A ∆ ψ S t e r n 2 (15)

Biesheuvel et al. (2014) improved the mD model by making the ionic attraction term dependent on total ion concentration in the carbon pores. The new model significantly improved the agreement of calculated results and experimental data for the operation of CDI processes. Gabitto and Tsouris (Gabitto and Tsouris, 2016) presented a model to simulate the CDI process in heterogeneous porous media comprising macro- and micropores. The authors used a two-step volume averaging technique to derive the averaged transport equations at both length scales. The authors derived a one-equation model, valid for both phases, based on the principle of local equilibrium [(del Rìo and Whitaker, 2002a; Whitaker, 1986; Whitaker, 1991; Quintard and Whitaker, 1995), among others] and approximate models for each phase (Gabitto and Tsouris, 2016). Gabitto and Tsouris (Gabitto and Tsouris, 2018) derived using a more precise theoretical procedure the two-equation model equations for the dual-porosity porous medium. The authors found out that these equations coincided with the equations originally proposed by them in their previous work (Gabitto and Tsouris, 2016).

Gabitto and Tsouris (Gabitto and Tsouris, 2018) derived the following equations for the liquid (γ-phase) and the solid (κ-phase) in the macropores: ε γ ∂ c i , γ γ ∂ t = ∇ →   •   ε γ D ¯ ¯ i , γ γ • ∇ → c i , γ γ + ∇ → • ε γ D ¯ ¯ i , γ κ • ∇ → c i , κ κ + z i ∇ → • U ¯ ¯ i , γ γ ε γ c i , γ γ ∇ → ψ γ γ + z i ∇ → • U ¯ ¯ i , γ κ ε γ c i , γ γ ∇ → ψ κ κ − ε γ v γ → γ • ∇ → c i , γ γ + 1 V M ∫ A γ κ n → γ κ • D i ∇ → c i , γ + z i   c i , γ   ∇ → ψ γ d A i n γ − p h a s e (16) ε ∝   ε κ ∂ c i , κ κ ∂ t = ∇ → • ε ∝ ε κ   D ¯ ¯ i , κ γ   •   ∇ → c i , γ γ + ∇ →   •   ε ∝ ε κ   D ¯ ¯ i , κ κ   •   ∇ → c i , κ κ + z i   ∇ →   •   U ¯ ¯ i , κ κ   ε κ c i , κ κ   ∇ → ψ κ κ + z i   ∇ →   •   U ¯ ¯ i , κ γ   ε κ c i , κ κ   ∇ → ψ γ γ + 1 V M ∫ A κ γ n → κ γ •   ε α   D ¯ ¯ i , e f f   •   ∇ →   c i , κ + z i   D i   ∇ →   •   c i , κ   ∇ → ψ κ   d A i n κ − p h a s e (17) Here, ε ∝   ε κ is the void fraction in the micropores ( V α / V M ), V α is the volume of micropores, V γ is the volume of macropores, V M is the macroscale representative elementary volume (REV), A γ κ is the macropores interphase liquid-solid area, ε γ is the void fraction in the macropores ( V γ / V M ) , ε κ is the solid-liquid volume ratio in the macropores, A = i , j k are the various transport coefficients that control the flow of the i-species, c i , j j and ψ j j are the intrinsic j-phase averages for the i-species in the j-phase given by: c i , j j = 1 V j   ∫ V j c i , j   d V = ε j c i (18) ψ j = 1 V j   ∫ V j ψ   d V = ε j ψ (19)

Here, c i and ψ , are the phase averages of the species concentrations and the electrostatic potential. The last integral term in the right-hand-side of Eqs 16, 17 represents the liquid-solid interfacial transport of the average properties. Gabitto and Tsouris (Gabitto and Tsouris, 2018) proved that the integral terms are related by: 1 V M ∫ A γ κ n → γ κ   •   D i ( ∇ →   c i , γ + z i   c i , γ   ∇ → ψ γ )   d A = − 1 V M ∫ A κ γ n → κ γ •   ε α D = i , e f f •   ∇ →   c i , κ + z i   D i   ∇ →   •   c i , κ   ∇ → ψ κ   d A , (20)

The two-equation model derived by Gabitto and Tsouris (Gabitto and Tsouris, 2018) can be compared to the phenomenological formulation of Biesheuvel et al. (2011), Biesheuvel et al. (2012) by using the following assumptions:

• There are no transport processes inside the micropores. This assumption is justified on the grounds of the micropores sizes. The characteristic time for a transport process is proportional to the square of the micropores characteristic length, for example, t c = l m 2 /   D i .

Introducing all the assumptions listed above into Eqs 16, 17 leads to: ε γ ∂ c i , γ γ ∂ t + ε ∝   ε κ ∂ c i , κ κ ∂ t = ∇ →   •   ε γ   D ¯ ¯ i , γ γ   •   ∇ → c i , γ γ + z i   ∇ →   • U ¯ ¯ i , γ γ   ε γ c i , γ γ   ∇ → ψ γ γ − ε γ v γ → γ   •   ∇ → c i , γ γ i n γ − p h a s e (21)

Equation 21 is like Eq. 14 presented by Biesheuvel et al. (2011), Biesheuvel et al. (2012). However, these equations are not identical as the time variation of the average concentration in the solid (κ-phase) is not equal to the average concentration in the micropores ( c i , α α ) as it appears in Eq. 14. Both concentrations are related through the definition of the point concentration in the κ-phase, but not to the average concentration (Gabitto and Tsouris, 2018). c i , κ = ε α 〈 c i , α 〉 α (22)

The two kinds of representative volumes (REVs) are shown in Figure 5.

The small REV allows calculation of averaged concentrations inside the solid particles (R_m). Using the Donnan model we can assume that the ion concentrations are constant inside the solid particles ( c i , α α ). However, the average at the macroscale associates the average value of the property in a macro representative elementary volume (R_M) to every point in the domain. Even if the properties are constant in every particle the property values will change for particles in contact with electrolyte at different positions within the R_M. So, the average of the properties at the macroscale will be different from the averages inside the particles.

It is important to stress that not only the three assumptions listed above must hold for Eq. 21 to be valid. Equation 16 and Equation 17 were derived by Gabitto and Tsouris (Gabitto and Tsouris, 2018) using also other constraints that are listed in reference (Gabitto and Tsouris, 2018). These constraints need also to be satisfied for Eq. 21 to be valid. Gabitto and Tsouris (Gabitto and Tsouris, 2018) also derived expressions to calculate the values of the transport coefficients appearing in Eq. 21 using the closure procedure in the volume average technique.

In recent years, there has been an exponential growth in the number of articles discussing different cell configurations and mechanisms to improve the CDI and related processes. Some researchers review critically the literature to summarize recent developments [(Suss et al., 2015; Tang et al., 2019a), among others]. New cell structures comprise the use of ion exchange membranes (MCDI) (Biesheuvel et al., 2010; Biesheuvel and van der Wal, 2010; Galama et al., 2016; Hassanvand et al., 2017; Tang et al., 2019b), flow through electrodes (FTE CDI) (Suss et al., 2012; Guyes et al., 2017), fluidization electrodes (Doornbusch et al., 2016), flow-electrode capacitive deionization (FCDI) without (Hatzell et al., 2014), and with membranes (Jeon et al., 2013; Porada et al., 2014; Choo et al., 2017; Rommerskirchen et al., 2017; Ma et al., 2018).

There has been increased interest in the use of ion exchange membranes in combination with porous electrodes. Several cell configurations have been proposed for water desalination using ion exchange membranes placed in front of porous electrodes allowing the passage of counterions while restricting the flow of co-ions. The process received the name membrane capacitive deionization (MCDI). Biesheuvel et al. (2010) proposed a theory for MCDI including the electrostatic EDLs inside the micropores of porous particles, but also incorporating the role of the macropores, i.e., the void space between the particles. In MCDI the macro porosity becomes a salt storage reservoir that increases the ion storage capacity of the electrodes. The authors used modified Donnan model to simulate ion transport between macro- and micropores. The cell configuration studied by the authors is shown in Figure 6. There is a central water channel, and each electrode is separated from the water channel by an ionic membrane. Horizontal (x-axis) mass flows were only considered by using an average flow while vertical (y-axis) were considered by dividing the domain into k-ideally stirred volumes without diffusion mass transport. This assumption is justified as the height of electrodes, membranes, and water channel is much higher than the respective thicknesses. Convective transport in the y-direction inside the water channel is also higher than the corresponding diffusive transport.

In the membrane Biesheuvel et al. (2010) assumed linear ionic mass flows. The total ionic flow through the membrane ( J i o n ,   m e m ) is calculated using: J i o n , m e m y = − D m e m δ m e m ∆ C T m e m y + w X ∆ ψ m e m y (23)

Here, D m e m is the effective diffusivity of the salt inside the membrane, δ m e m is the membrane thickness, C T m e m is the total concentrations of ions ( C c a t i o n m e m + C a n i o n m e m ), w is the sign of the membrane chemical charge, X is the magnitude of the membrane chemical charge, and ∆ ψ m e m is the potential drop calculated across the membrane. On both solution/membrane interphases Biesheuvel et al. (2010) used a Donnan equilibrium equation to calculate the Donnan potentials between solution and membrane. ∆ ψ D o n n a n = ψ m e m − ψ s o l = sin ⁡ h − 1 w   X 2   C s a l t s o l (24) Here, C s a l t s o l is the salt concentration in the solution just outside the membrane.

The ionic current flowing through the membrane is given by: I i o n , m e m y = − D m e m δ m e m 〈 C T m e m 〉 y + w X ∆ ψ m e m y (25) Here, C s a l t s o l is the salt concentration in the solution averaged at both ends of the membrane ( C s a l t s o l = C a n o d e   s a l t s o l + C c a t h o d e   s a l t s o l / 2 ).

In the water channel the authors used the following equation: ∂ C j C h ∂ t = − 2 J i o n δ c h − 2 k τ c h C j c h − C j − 1 c h (26) Here, C j c h   is the salt concentration in the j-ideally stirred volume, δ c h is the water channel thickness, and τ c h is the residence time in the water channel measured in the y-direction.

In Eq. 26 the y-direction concentration changes by the convective transport in the y-direction, last term in the right-hand-side, and by the average ionic flow in the x-direction.

The flow in the electrodes was simulated by using a simplified version of Biesheuvel et al. (2014) model. For example, ionic mass flows are calculated by: ε m   ∂ C m , c a t i o n , j e l e c ∂ t + ∂ C m , a n i o n , j e l e c ∂ t + 2   ε M   ∂ C , M , j e l e c ∂ t = −   J i o n δ e l e c − k   τ e l e c   C j e l e c   − C j − 1 e l e c   (27) Here, δ e l e c is the electrode thickness, and τ e l e c is the residence time in the electrode measured in the y-direction, C , M , j e l e c is the salt concentration in the macropores in the j-ideally stirred volume, and C m , a n i o n , j e l e c and C m , a n i o n , j e l e c , are the salt concentrations in the micropores of the j-ideally stirred volume for the anion and cation, respectively.

Galama et al. (2016) studied theoretically and experimentally the membrane potential generated by different ionic concentration in both sides of an ion exchange membrane. The authors reported an elegant derivation of the Teorell-Meyer-Sievers (TMS) equation The authors extended the theory considering the presence of stagnant layers on both sides of the membrane and differences in ionic mobility between co-ions and counterions inside the membrane.

Tang et al. (2019b) studied experimentally and theoretically saltwater desalination by an MCDI technique. The salt removal capacity of carbon electrodes is limited due to the high ionic concentration of seawater. The authors proposed to overcome this limitation by applying over-potentials, close to the Faradaic limit, and by reversing cell polarity in the ion-desorption step. The proposed technique (OP-MCDI-RP) was found to be acceptable for seawater desalination (Tang et al., 2019b). Tang et al. (2019b) used a cell like the one depicted in Figure 6. The authors simulated the process by a combined approach. They simulated the operation of the membrane by using Biesheuvel et al. (2010) membrane model, Eqs 23, 25. They simulated the flow inside the water channel by 1-D mass balance in the y-direction using the assumption of average x-direction ionic flows. The operation of the porous electrodes was simulated using the one-equation model presented by Gabitto and Tsouris (Gabitto and Tsouris, 2018). ∇ _ • η _ _ * • ∇ _ ψ * = − ε κ a v V T σ α β κ − ε F V T ∑ i z i c i * (28) ε ∂ c i * ∂ t = ∇ _ • D _ _ * • ∇ _ c i * + z i ∇ _ • U _ _ * • c i * ∇ _ ψ * (29) Here, η _ _ * is the global effective permittivity, D _ _ * is the global effective diffusivity, U _ _ * is the global effective mobility, ε is the global liquid phase void fraction, a * is the one-equation model average of the a variable, and σ σ 〉 σ α β κ is the average α−β interphase EDL charge density.

Tang et al. (2019b) reported very good agreement between experimental data and simulation results using the capacity of the Stern layer ( C S t e r n ) as an optimized parameter. Similar agreement has also been reported by other researchers [(Biesheuvel et al., 2010; Biesheuvel et al., 2014; Rommerskirchen et al., 2017), among others]. The authors concluded that at the high potentials used in the CDI/MCDI process the capacity of the EDL ( C E D L ) is approximately equal to the capacity of the Stern layer as (Tang et al., 2019b): 1 C E D L = 1 C S t e r n + 1 C D L = 1 C S t e r n + 1 C D H   cosh ∆ ψ D l 2 ≅ 1 C S t e r n (30) Here, C D H is the Debye-Hückel capacity, C D L is the diffuse layer capacity, and ( ∆ ψ D L ) is the potential drop in the diffuse layer.

Besides the use of membranes to improve the CDI process other configurations have been tried. Suss et al. (2012) proposed a cell configuration where the brackish feed flows through the electrodes in the direction of the applied electric field (FTE CDI). The authors used an electrode made up of carbon aerogel monoliths. They found that this configuration exhibits low flow resistance and high capacitance. They modelled the mass and charge balances using equations like the 1D heat equation, but the mass balance contains a sink term dependent on potential that couples both equations. Suss et al. (2012) concluded that the desalination process operates with a characteristic timescale equal to the cell’s RC constant.

Guyes et al. (2017) presented a 1-D model of water desalination by an FTE CDI cell based on the modified Donnan electric double layer theory and simplified boundary conditions. The authors reported good agreement for a comparison between experimental data measured using a laboratory FTE CDI cell and their simulation results. This agreement was obtained by simulating the experimental data using the capacity of the Stern layer, the micropore volume, and the attraction energy term as fitting parameters. Guyes et al. (2017) formulated mass balances for each ion species using Eq. 14 without a redox reaction and considering convective flow. In the case of a monovalent, binary electrolyte the total mass balance is calculated by adding the cation and anion differential equations while the net charge equation is derived by subtracting the anion equation from the cation one, to give: ∂ C e f f ∂ t = ε M ∂ C M ∂ t + ε m 2 ∂ C i o n s m ∂ t = − v   ∂ C M ∂ x + D M   ∂ 2 C M ∂ x 2 (31) ε m ∂ σ i o n i c ∂ t = 2   D M   ∂ ∂ x C M ∂ ψ ∂ x   (32) where C M is the salt concentration in the macropores, C i o n s m is the total ionic concentration in the micropores, and σ i o n i c is the net ion charge in the micropores ( C + m − C − m ). The authors used the modified Donnan model, Eq. 1, to relate the concentrations in the macro and micropores.

Rommerskirchen et al. (2017) reported a constant voltage a process model for a continuous, steady state FCDI process. The process involves a single module device operating in single-pass model. The system comprises, a central water channel plus two flow-through electrodes each one separated from the water channel by an ion exchange membrane. The electrodes are made of flowing porous conductive porous particles. The membranes are considered ideally permselective so, only allow flow of the counterions into the flow electrodes. Reduction of flow through the membrane can occur due to mass transfer limitations in the water channel-membrane interphase [concentration polarization (Ma et al., 2018)]. A limiting current density occurs when the ion concentration on the membrane surface is reduced to zero. The authors allowed for this phenomenon by calculating the limiting flow in a stagnant boundary layer on the water channel side. This limiting flow is compared to the unlimited flow through the membrane and the smaller value is selected.

The electrodes are considered as homogeneous carbon suspensions with constant particle concentration and void fraction. Charge neutrality is assumed outside the particles, macropores volume, while there is net ionic charge inside the porous particles, micropores volume, balanced by the induced charge on the particles’ solid matrix. The storage of ions in the micropores of the carbon flowing particles was modeled using the modified Donnan model of Biesheuvel et al. (2011), Biesheuvel et al. (2012). Inside the porous particles the potential drop is divided between a Stern layer region without ionic charge and a region with potential drop given by the modified Donnan model. The potential drop through the various parts of the cell was considered by relating the cell potential to the sum of all ohmic potential losses. V c e l l = ∆ ψ S t e r n A n o d e + ∆ ψ D A n o d e + ∆ ψ S t e r n C a t h o d e + ∆ ψ D C a t h o d e + ∆ ψ L o s s T o t a l + ∆ ψ T r a n s T o t a l (33)

Here, V c e l l is the voltage applied between anode and cathode, ∆ ψ S t e r n i + ∆ ψ D i , represent the potential drops in the porous particles in the i-electrode, ∆ ψ L o s s T o t a l gives the ohmic losses by flow through the ion exchange membranes and the flow channel, and ∆ ψ T r a n s T o t a l is the potential gradient for ionic transport through the flow electrodes and the water channel.

Rommerskirchen et al. (2017) reported good agreement between experimental data and simulation results calculated using the Stern layer capacity as the only fitting parameter. The authors also studied an operating process comprising 2 cells one operating in charging mode and the other in discharging mode. They reported a strong influence of flow rate on the desalination performance.

Tang et al. (2019c), Tang et al. (2020) studied experimentally and theoretically the influence of applied potential and electrolyte concentration on the FCDI and electrodialysis (ED) desalination processes. The model presented by the authors relied upon the same general framework presented by Rommerskirchen et al. (2017) but the operation of the flow electrodes was based upon the “hybrid” model for flow electrodes by Gabitto and Tsouris (Gabitto and Tsouris, 2019), and the membrane operation was modelled after Galama et al. (2016).

Gabitto and Tsouris (Gabitto and Tsouris, 2019) simulated the operation of slurry carbon electrodes by combining a phenomenological approach for the point equations with a two-steps volume averaging method. This ‘hybrid approach’ led to the formulation of the equations describing the operation of flow electrodes. The authors used the macroscale REV shown in Figure 6. The γ-phase and κ-phase represent the liquid and solid particles, respectively. ε γ ∂ c γ γ ∂ t = ∇ →   •   ε γ D ¯ ¯ e f f •   ∇ → c γ γ − ε γ   v → γ •   ∇ → c γ γ − ε κ   v → κ •   ∇ → c κ κ − ε κ   ∂ c κ ∂ t (34) 0 = ε γ   ∇ →   • U ¯ ¯ e f f • c γ γ   ∇ → ψ γ γ − ε κ   v → κ •   ∇ → ρ κ κ − ε κ ∂ ρ κ ∂ t (35) Equation 34 and Eq. 35 are valid in the γ-phase. In Eq. 34 the authors defined a total mass parameter, c j = c + , j j + c − , j j 2 ; c i , j j is the specific volume averaged concentration of the i-species in the j-phase; ψ γ γ is the specific volume averaged potential in the γ-phase; ε γ and ε κ are the volume fractions of both phases in the REV; D ¯ ¯ e f f and U ¯ ¯ e f f , are the effective diffusivity and effective mobility tensors for the salt concentration in the γ-phase; respectively. In Eq. 35 ρ_j is a charge parameter, ρ j = c + , j j − c − , j j 2 , where j is the phase index, j = γ, κ. In Eqs 34, 35 the assumption of electro-neutrality in the macropores ( c + , γ γ ∼   c − , γ γ ) has been used.

The system of Eqs 34, 35 requires calculation of the derivatives of c _κ and ρ _κ . Gabitto and Tsouris (Gabitto and Tsouris, 2019) proposed that: c κ = ε α   〈 c i , γ 〉 γ   cosh ∆ ψ D o n n a n , (36) ρ κ = − ε α   〈 c i , γ 〉 γ   sinh ∆ ψ D o n n a n . (37)

The Donnan potential ( ∆ ψ D o n n a n ) is calculated using Rubin et al. (2016): ψ β − 〈 ψ γ 〉 γ = ψ β − 〈 ψ α 〉 α + 〈 ψ α 〉 α − 〈 ψ γ 〉 γ = ∆ ψ S t e r n + ∆ ψ D o n n a n (38) ∆ ψ S t e r n = 2 F 〈 c i , γ 〉 γ   sinh ∆ ψ D o n n a n a v m   C S t e r n V T (39) ψ β − 〈 ψ γ 〉 γ = ∆ ψ D o n n a n + 2 F 〈 c i , γ 〉 γ ⁡ sinh ∆ ψ D o n n a n C S t e r n a V m V T (40)

Here, C _Stern is the surface capacitance of the Stern Layer, a V m is the specific area of the pores in the solid particles (A _αβ / V _m ), ψ _β is the solid particle potential, and V _m is the REV inside the porous particles. The Donnan potential used in Eqs 36, 37 is calculated by solving Eq. 40.