A Non-Electroneutral Model for Complex Reaction-Diffusion Systems Incorporating Species Interactions

In this study we develop a general framework for describing reaction-diffusion processes in a multi-component electrolyte in which multiple reactions of different types may occur. Our motivation for this is the need to understand how the interactions between species and processes occurring in a complex electrochemical system. We use the framework to develop a modiﬁed Poisson-Nernst-Planck model which accounts for the excluded volume interaction (EVI) and incorporates both electrochemical and chemical reactions. Using this model, we investigate how the EVI inﬂuences the reactions and how the reactions inﬂuence each other in the contexts of the equilibrium state of a system and of a simple electrochemical device under load. Complex behavior quickly emerges even in relatively simple systems, and deviations from the predictions of ideal solution theory, together with how they may inﬂuence the behavior of more complex system, are discussed.

Electrochemical energy storage devices play a crucial role in the modern world, having enabled the development of a wide range of portable and mobile devices in a vast range of applications. They also have significant future potential in facilitating a shift away from environmentally damaging fossil fuels as our primary source of energy, through the electrification of transport and as load balancing for the variability suffered by most forms of renewable energy. However, modeling these systems can be challenging because the overall behavior is typically the emergent result of a large number of processes and interactions at the microscopic scale, making linking the microscopic behavior to the macroscopic performance complicated. Furthermore, individual processes may themselves be complex, so simplifications have to be made if we wish to understand the device behavior at macroscopic length and time-scales.
By way of example, the particular system in which we are interested is that of a lithium-sulfur (LiS) cell, a promising post-lithium-ion technology with both an expected practical energy density of 500-600 Wh kg −1 and a lower raw materials cost. 1,2 The overall discharge process of a LiS cell involves the reduction of solid phase S 8 to solid phase Li 2 S, according to the reaction While the overall process is bound by the dissolution of S 8 and the precipitation of Li 2 S, the intermediate steps occur between species in the solvent/electrolyte phase, involving the electrodissolution of lithium from the anode and a range of electrochemical and chemical reactions involving a number of ionic sulfur species at the cathode. While these types of process are not uncommon in traditional (i.e. non-intercalation) battery chemistries, the sheer number of species and intermediate elementary reaction steps involved, together with the integral role played by chemical reaction processes, make understanding the LiS mechanism complex. 3 The common approach to modeling LiS cells is similar to that taken for many electrochemical devices, with the reduction of the cell structure to a one-dimensional model in which the porous structures are homogenized, ideal solution theory is applied and electroneutrality of the electrolyte is assumed. [4][5][6][7][8][9] This approach has significantly improved our understanding of LiS behavior, but the underlying simplifications to some extent limit our ability to look at the system below the homogenized level. As a consequence, it is difficult to probe how z E-mail: geraint.minton@oxisenergy.com species interactions affect reaction processes, how the processes themselves interact or to develop an understanding of how the geometrical structure might affect them.
In particular, homogenized models do not explicitly account for the electric double layer (EDL), formed by the attraction of counter-ions to a charged surface and the repulsion of co-ions from it. This surface charge might be the result of an externally applied voltage, such as in a capacitor, or be generated internally, as in a battery, but it is almost always present in an electrochemical device, and therefore so is the EDL. The EDL therefore forms the interface between the electrode and the electrolyte, being where all surface interactions between the two (e.g., electrochemical reactions) occur.
At all but the lowest concentrations and voltages, the structure of the EDL differs significantly from that of the bulk electrolyte. These changes, driven by forces causing species migration and diffusion and complicated by species interactions, have measureable effects on the macroscopic properties of both equilibrium and dynamic systems without reactions, despite the fact it exists on a nanometer lengthscale, which can be orders of magnitude smaller than the overall system. [10][11][12] Following on from this, it is likely that accounting for species interactions will impact the behavior of systems in which reactions occur, such as the LiS cell, wherein there are a large number of different species interacting in a poorly understood manner. In order to describe this, a framework for building the model is required.
An array of methods exist for describing a reaction-diffusion system, ranging from continuum ideal solution theory through to molecular dynamics, but to describe an electrochemical device the need to describe macroscopic length scales and long time-scales places a constraint on the complexity that can be accounted for. For example, molecular dynamics may provide highly detailed information on the structure of an electrolyte, but it is limited to describing very short time-spans and a limited number of molecules.
In this present work, we outline an approach to building general reaction-diffusion models in which species interactions are incorporated, the EDL is described, and macroscopic transient behavior (e.g., voltage curves) can be estimated. This has been developed in the context of understanding the mechanism driving a LiS cell, but the framework itself is of general form. As such, it is adaptable to a range of electrochemical devices which may benefit from an understanding of how interactions between species or between reaction processes may affect the system behavior. This includes the double layer structure of supercapacitors or the behavior of pseudocapacitors, certain types of fuel cell (for example, the direct methanol fuel cell has a relatively E3277 complex reaction mechanism), the electrolyte component of intercalation batteries such as lithium-ion, and other electrochemical systems in which homogenization may oversimplify the electrolyte-surface interaction, such as the lithium-air battery. Given the generality of the framework, herein we only discuss its application in general terms, albeit tied to examples drawn from our experience with LiS.
The derivation stems from the free energy of the system and how it links to the structure, dynamics and reactions in the system, although to manage the complexities introduced by describing the EDL, we ultimately apply relatively simple models for the components of the system. This should provide a clear basis around which improvements can be made, since alterations to the free energy expression feed through to the rest of the model.
In our present version, the equations resulting from the theory take the form of a modified Poisson-Nernst-Planck (mPNP) model, in which species transport is described by a modified Nernst-Planck equation and the electric field is solved for using the Poisson equation. Surface electrochemical reactions are accounted for, as are chemical reactions in the bulk electrolyte or at the surfaces, where the latter may be coupled to a precipitation model. For the purposes of the reaction processes, we homogenize the surface in the current treatment, meaning that the formation of precipitate is implicit, rather than being explicitly described in the modeling domain.

Theory
We consider a general electrolyte system in which ions of type i are treated as charged hard spheres with valence z i and diameter d i that are immersed in a continuum solvent with relative permittivity r . The ions are allowed to react and undergo chemical transformations. To develop an approximate dynamical model for this system, we assume that its evolution can be derived from an underlying equilibrium free energy functional, which is dependent on the concentration distribution of species throughout the system. The dynamics of the species concentrations are taken to be related to gradients of the free energy functional, subject to local conservation constraints, so that the free energy monotonically approaches a minimum with time.
The resulting model is composed of four main components: an equilibrium free energy functional, a species transport model, a reaction model and a precipitation model. While forms of each part of the model have appeared separately in the literature, we present them here together within a coherent formalism in order to present a consistent theory and to understand directions in which each aspect of the model can be improved and how this might be achieved.
In the next section, we develop an approximate free energy functional to describe the equilibrium electrolyte system. This accounts for their relative formation free energies, in addition to their mutual interactions, and we demonstrate how this functional can be used to describe chemical reaction equilibria. Next, we present a phenomenological approach to extending this equilibrium theory to describe the transient behavior in non-equilibrium systems, including species transport and reactions. In the final part of the section we outline the precipitation model that we couple to the system. Approximate free energy functional.-We first develop an approximate expression for the free energy functional for the electrolyte system. The total Helmholtz free energy F of the system is written as the sum of two contributions The first is the free energy of a reference system F ref , which describes the contributions from entropy and non-electrostatic interactions, and the second is the contribution due to electrostatic interactions F el . In principle, these two contributions are closely coupled together, however, we make the approximation that their individual effects can be linearly added together. A number of approaches have been developed to describe the reference system, ranging from a simple ideal gas description through to more sophisticated density functional theories. 13,14 In order to man-age the overall model complexity, herein we choose to use a simple model for the electrolyte structure by working within the local density approximation (LDA). Within this scheme, F ref is a function only of the local concentrations c i (r) of each species in the system: [ 3 ] where f ref is the Helmholtz free energy density of a uniform reference fluid. While the use of the LDA is known to cause errors in the prediction of the detail of the EDL structure 15 and also to break down at high electrode potentials, 16 the resulting model structure is much more amenable to describing the long time-scale transient behavior of a system. It is also worth noting that ideal solution theory, which leads to the Poisson-Boltzmann theory and underpins most electrochemical modeling to date, is of the LDA form, so nothing has been lost by making this choice, however, it does represent an opportunity for future model development.
The free energy density f ref is typically separated into an ideal (or entropic) component and a residual component, which describes the non-electrostatic interactions within the system. We account here for the excluded volume interaction (EVI) of the species using a form of the van der Waals equation of state for mixed hard spheres, 11,17 [4] where μ i is the standard state chemical potential of particle i (defined as an ideal solution of non-interacting species c i at a concentration of 1 M at temperature 298 K and pressure 1 bar), k B is the Boltzmann constant, and T is the temperature. The first term is the formation free energy of each component in the system. The second term represents the ideal entropic contribution to the free energy. The final term in Equation 4 is the residual contribution to the free energy due to the EVI, in whichv ii is defined as where v ii is the excluded volume per particle between particles i and i . This is itself defined in terms of the solvated diameters d i of the particles in question: This includes an empirical modification of 3 3 2π ≈ 0.78, which is used as to adjust the free energy density in this model to more closely resemble that of the lattice model used by Bikerman. 18 The lattice model is inapplicable for particles of different sizes, but does not suffer as badly from the aforementioned breakdown of the LDA as the more accurate Boublik-Mansoori-Carnahan-Starling-Leyland model, 19,20 which is known to overestimate the influence of the EVI at high electrode potentials within the LDA. 16 The energy of the electrostatic interactions is given by [7] where φ(r) is the electrostatic potential, (r) is the fixed charge density (e.g., the electrode surface charge) and Q(r) = i z i e 0 c i (r) is the mobile (ionic) charge density, in which e 0 is the elementary charge and z i the valence of species i. In general, the addition of electrostatic interactions will alter the manner in which the particles in the system organize themselves, and, consequently, it will alter the manner in which non-electrostatic interactions contribute to the free energy. However, in this work, we neglect this coupling.

E3278
Journal of The Electrochemical Society, 164 (11) E3276-E3290 (2017) The complete form of the Helmholtz free energy functional which we use is the combination of Eqs. 3, 4 and 7: [8] From knowledge of the free energy functional, it is possible to determine all equilibrium thermodynamic and structural properties of the system. The two quantities in which we are interested are the electrostatic and electrochemical potentials, as these together dictate the electrolyte structure. The electrochemical potential of a species is equal to the change in the free energy with respect to its concentration: If the species volumes are all zero, the electrochemical potential of an ideal solution is recovered. Specifically accounting for the contributions of additional interaction energies in the Helmholtz free energy will lead to additional terms in the electrochemical potentials, from where they will follow through to the rest of the model.
Chemical reactions equilibria..-We now consider the situation where species interconversion can take place. This can occur either chemically, electrochemically or even physically (e.g., precipitation, which converts a dissolved component in solution to a solid), shown pictorially in Figure 2. These reactions can be written in the general form: [10] where the subscripts k and l represent reactant and product species, respectively: b k j is the number of moles of reactant species B k consumed by reaction j, and b l j is the number of moles of product species B l produced by the reaction. The stoichiometric coefficient ν k j of a reactant species k in reaction j is equal to −b k j , while the stoichiometric coefficient ν l j of a product species l in reaction j is equal to b l j . One example is an electron transfer reaction Ox + n e e − Re [11] where Ox and Re represent a redox couple, and n e is the number of electrons e − transferred to Ox during the reduction process. Another example is the precipitation of a species A from solution: [12] where A (d) is the component in solution and A (s) is the same component in the solid phase. The overall amount that reaction j proceeds in a forward or reverse direction is quantified by its extent of reaction j . If the system is initially charged with n 0 i molecules of type i, then the number of molecules n i in the system when the reaction has moved forward by j is Conditions for equilibrium.-The equilibrium structure and thermodynamic properties of the system can be determined by minimizing the free energy functional, while maintaining any physical constraints on the system. We consider a closed system, i.e. one in which no mass enters or exits the system. In this case, the only manner in which the initial number of molecules of type i can change is through chemical reaction. Consequently, n i , the total number of particles of type i in the system, is directly related to the extents of all possible reactions in the system: For non-uniform systems, these relations impose constraints on the species concentration profiles: [15] where ξ j (r) is a local extent density of reaction j at position r in the bulk, andξ j (r) is a local extent density of reaction j at position r on the surface of the system.
Minimising the free energy in Equation 8 subject to these constraints leads to the conditions required for the system to be at equilibrium. Using the method of Lagrange multipliers, where we introduce the Lagrange multipliers λ i to maintain the constraints and find the minimum of the functional we find the relations δF δc i (r) − λ i = 0, [17] δF δφ(r) = 0, [18] δF δξ j (r) + i λ i ν i j = 0, [19] δF At equilibrium, the system must satisfy Equations 17-20. Equation 17, combined with Equation 9, shows that the Lagrangian multipliers can be identified with the electrochemical potentials of each species λ i = μ i (r). [21] This indicates that the electrochemical potential of every species must be uniform throughout the system (even if the concentration profiles are not) at equilibrium. From Equation 18, we find that the shape of the electric field can be determined from the principle that it always adjusts itself to minimize the free energy. This leads directly to the Poisson equation: −∇ · [ 0 r ∇φ(r)] = i z i e 0 c i (r) + (r). [22] Equation 19 gives the condition for chemical reaction equilibria. If we define a local affinity A j (r) of reaction j Journal of The Electrochemical Society, 164 (11) E3276-E3290 (2017) E3279 then this condition for reaction equilibria is given by A j (r) = 0. [24] From this relation, we see that the affinity characterizes the deviation of a reaction from equilibrium.
The dynamical model.-To describe the transient behavior of a system that is out of equilibrium, we assume that it evolves in a manner that tries to decrease its free energy, subject to physical constraints (e.g., conservation of mass, etc.). We expect that the rate of change of the system will be related to the gradients of the free energy. In our description of the equilibrium system, the species concentration distributions are considered independent variables and the electrostatic potential is assumed to react instantaneously to changes in the species concentrations such that it always minimizes the free energy. It is therefore quasi-time-independent and described by the Poisson equation. Consequently, the rate of change of F with time is given by ∂t . [25] The rate of change of the concentration profiles is restricted by the physical constraint of the local conservation of species, which for species i is given by where J i (r, t) is the local molecular flux of species i and R j is the rate of reaction j. In addition, for a closed system, this partial differential equation is subject to the boundary condition wheren is a unit vector pointing in an outward normal direction from the system surface, andR j (r, t) is the rate of reaction j at position r on the surface of the system. Models for J i (r, t), R j (r, t) andR j (r, t) are required to complete the theory. Substituting Equation 26 for the time derivative of the species concentrations, Equation 25 becomes Applying the divergence theorem and substituting Equation 23 into the result yields In the following, we will develop phenomenological expressions for the species flux and reaction rates to guarantee the decrease of the free energy with time.
Species flux.-Motivated by the standard expression for the flux J i of species i as given by Maxwell-Stefan diffusion, we write: = −D i c i (r, t)∇βμ i (r, t) [31] This form guarantees that the free energy will always decrease or remain constant with time. We note that this phenomenological expression for the flux can also be motivated from more fundamental statistical mechanical arguments. The one-body density (and hence the free energy) is a unique function of a time-dependent external field, 21 and under the assumption that two-particle correlations are identical in equilibrium and non-equilibrium fluids, the dynamics of a species i in the system can be approximated by 22 where ω c i = k B T /(m i D i ) is the collision frequency of particle i, m i is its mass and D i its self-diffusion coefficient. The high collision frequency of particles in the electrolyte means that the second order time-derivative is negligible.
Combining Eqs. 9, 26 and 30, we arrive at a modified Nernst-Planck (mNP) equation describing the species transport in the cell: In the case of an ideal solution μ res i , the residual component of the electrochemical potential which accounts for all non-electrostatic non-ideality, is zero and so we recover the standard Nernst-Planck equation.
Reaction kinetics.-From physical considerations, the affinity and the reaction rate always have the same sign, which implies that the reaction terms in Equation 29 are also always negative or zero. From this result, we can see that the free energy of the model system spontaneously decreases to its minimum value (the equilibrium state) with time.
The local affinity quantifies how far a particular reaction is from equilibrium at a point in space as well as the direction the system must move to reach equilibrium, under which condition its value is zero. The affinity is related to the ratio of the forward and reverse elementary rates, R f j and R r j respectively, by the relationship 23-25 = exp(βA j (r)), [34] and, furthermore, the two elementary reaction rates are related to the overall reaction rate: where the connection between the extent of the reaction and its rate is also indicated. Combining the previous equations we can write the overall rate in terms of the affinity and one of the elementary rates: The same set of relations applies to the surface reaction ratesR j (r), and we henceforth mean the use of the variable R to imply both the surface and bulk reaction rates.
The reaction expressions.-In order to apply Equation 36 we require an explicit expression for the affinity, for which expressions for the electrochemical potentials of the reactants are needed. For species in the solvent phase these are given by Equation 9, but in the case of an electrochemical reaction we must also specify the electrochemical potential of an electron in the electrode phase. This is related to the Fermi energy E f of the surface: where φ el is the electrostatic potential of the electrode relative to its uncharged state. Since the electrode exists in the same system as an electrolyte, its electrostatic potential in the uncharged state must be equal to the electrostatic potential of the bulk electrolyte, which is also in an uncharged state. For this reason the potential of the electrode can also be read as the potential difference between the electrode surface and the bulk electrolyte, which is effectively any point in the electrolyte outside the EDL, where the electric field is zero. For the general reaction scheme of Equation 10, the affinity is given by the combination of Equations 9, 23 and 37: where the spatial dependence of the affinity, concentration, residual chemical potential and electrostatic potential have been removed for clarity, and the term G RXN,j = i ν i j μ i − n e, j E f is the standard Gibbs free energy change of reaction j. Substituting this expression into Equation 36, the rate expression becomes The exponential involving the Gibbs free energy is related to the standard state equilibrium constant of the reaction K j , which is the ratio of the standard state forward and reverse rate constants, k f, j and k r, j : Despite substituting in this expression, the forward rate is still unknown and so the rate cannot yet be determined. Making the assumption that the forward rate is a function only of the reactant species properties, together with a proportion of the term involving the electrostatic potential, the forward rate can be approximated as where γ is the transfer coefficient commonly found in the Butler-Volmer equation. This gives the final expression for the general reaction rate: This very general expression simplifies for each of the reaction types considered in the model, each of which is described briefly below.
Electrochemical reactions.-We assume that electrochemical reactions occur as n-electron single-step transfer processes in which a molecule of species Ox is reduced to a molecule Re or as the corresponding oxidation process, as stated in Equation 11. Since this reaction type requires a source of electrons, it is only able to occur at an electrode surface, meaning that the reactant species must be close enough to the surface for electron transfer to occur. Under the assumption that the edge of the reactant molecule (more specifically for this model, the edge of the molecule's solvation shell) must be in contact with the surface for this process to take place, only molecules whose centers are at the Stern layer are able to react (see Figure 1). For simplicity, we assume that the Stern layer is identical for all species in the system, the consequence of which is that the spatial dependence of the variables c i , μ res i and φ reduce simply to their values at the Stern layer. This mirrors its use as a fitting parameter in equilibrium continuum double layer modeling.
For the electrostatic potential, the terms corresponding to the potential difference between the surface and the Stern layer are replaced with the Stern layer potential difference φ S = φ el − φ(s). Under these assumptions, the rate expression reduces to which is of the form of the generalized Frumkin-Butler-Volmer kinetics model, 26 but includes the effect on the rate of the residual chemical potential.
Within the theory, the specific values of the forward and reverse rate constants are determined from the equilibrium rate constant. The standard state Gibbs free energy of the reaction is linked to the standard state reduction potential according to the relationship  Substituting this into Equation 40, we get We assume that the rate constants can be split under the same assumptions by which we split the reaction rate expression, while noting that k f, j must be larger than k r, j when the reduction potential is positive, to write their individual forms as Chemical reactions.-No electron transfer occurs in a purely chemical reaction, so n e, j = 0 and the rate equation reduces to the following expression: Within the model, the activity a i of a species i is defined as so the rate expression can be seen to reduce to the standard equation for the rate of a chemical reaction. In terms of the properties of the chemical reactions we are clear here that for the sake of simplicity, we assume that the system is isothermal and the heat of reaction is not accounted for. Including this would be a significant improvement to the model, particularly for LiS cells which exhibit interesting thermal behavior on cycling. 27 The precipitation model.-To model the formation of precipitates, we assume a simple hemispherical growth model with a fixed seed point density per unit surface area ρ i . Rather than explicitly model the precipitate species, we also homogenize the surface, ultimately using the precipitate coverage to modify the various surface reaction rates. A unit precipitate of species i, having radius r i , is depicted in Figure 3.
From geometric considerations, we can determine the total electrode area covered by all precipitate species as well as the surface area of precipitate i, both per unit area of the electrode surface, θ S and θ P i respectively: In these equations, the radius of the precipitate is not a constant, because the precipitation reaction alters the volume of the precipitate.
Assuming that precipitation occurs evenly over the entire surface of the hemisphere and is dependent on the species concentrations at the electrode surface reaction plane, the following relationship defines the rate of change of the radius with time: where m i is the molar volume of precipitate species i, and R θ i j = θ P i R j is the homogenized (effective) surface precipitation reaction rate.
The rate of the reaction is calculated by Equation 48. However, since the precipitate is a solid phase, its activity has unit value by definition.
The homogenization of the precipitate on the electrode surface introduces a limitation in terms of the maximum amount of solid phase which can be held per unit area of surface before the overlap of hemispheres causes the available surface area to become unphysically negative. The limit depends on the specifics of the molar volume of the solid phase species and the seed point density, but is of the order 7 g m −2 . This is comparable to, but smaller than, the amount found in a real LiS cell, which is of the order of 20 g m −2 or more, depending on the method of cathode production. The two situations are not directly comparable, however, because in a real cell the sulfur is in a 3D porous network while here we only have a flat surface. This will lead to differences in the electrochemical performance of the model cell compared to a real cell, because the highest solid phase loadings in the model correspond to there being very little electrochemically active surface area, inevitably altering the electrochemical reaction rates. This is less of an issue for a real cell, where the porous structure provides a very high active surface area. While there would be some advantage in building a 3D model for the description of precipitation, such a model would likely be prohibitively complex to solve with current techniques and so we make the assumption that lower dimensional models are capable of providing sufficient insight into the system performance.
Boundary conditions.-The modeling domain extends as far as the Stern layer/reaction plane at the surfaces of the model system. As mentioned, for simplicity we assume a single value for the Stern layer width for all species throughout the domain. For the mNP equations, the boundary condition is defined by the normal surface fluxes of the species, which are defined in terms of the effective reaction rates where R θ i j depends on the reaction type: In order to solve the Poisson equation, a reference potential is required in the system, which we take to be the anode surface, defined as zero volts. However, the since the modeling domain extends only as far as the Stern layer, we require the potential at this point to use as the boundary potential. To determine this, we note that its value is related to the potential gradient perpendicular to the surface and the potential of the electrode φ el : where φ S is the electrostatic potential at the Stern layer, located a distance s from the surface. The electric potential gradient is defined by the structure of the electrolyte, allowing the above equation to be used as an implicit boundary condition for the system. For all other surfaces in the system, the potential gradient is defined in terms of the surface charge density: The actual value of (r) is not directly calculable, but its rate of change is a function of the current density I flowing in the electrode

Condition
Explanation Defined initial volume of the unit precipitate and the charge generated by the electrochemical reactions What this equation implies is that the surface charge density, and therefore also the electric field in the electrolyte, and therefore also the voltage, is a function of the reaction rate and the current drawn, two facts which we know from real cells.
Initial conditions.-Finally in terms of model development, the initial state of the system must be defined. One of the benefits of this type of model is that the initial state is simple to define: the electrolyte is initially homogeneous, there is no potential difference between the electrodes and there is a set amount of precipitate on the electrode surfaces, specified by the volume of the unit precipitate, from which both the unit precipitate radius and the total quantity of precipitate per unit area of electrode can be calculated. The conditions are summarized in Table I.
Physically, these initial conditions are similar to those in a real cell, if it were possible to instantaneously fill a cell with electrolyte, or otherwise to prevent any reactions occurring until the cell was filled. A real cell has no initial voltage until electrochemical reactions spontaneously charge the electrodes or an external current is applied, meaning that there is an initial charging process before the cells can be discharged. The charging process is driven by the fact that the system is not initially in equilibrium, so spontaneous reaction-diffusion processes occur until equilibrium is reached. For this type of model, the fact that the electric potential can initially be assumed as zero decouples the mPNP equations, making defining the initial state simple. This is where homogenized cell models differ: for these, the initial state is defined after the initial charging of the electrodes, so the PNP equations are coupled and all species concentrations have to be such that the reaction Nernst potentials are all equal. 4

Results
At this stage, we wish to try and understand how the components of the system interact, rather than how the geometry affects the behavior. For this reason we reduce the geometrical complexity by considering only a simple slit-pore structure. In doing this, we know from the symmetry of system that its properties are invariant in the plane parallel to the surfaces, meaning we only have to solve for the variations in the perpendicular (z) direction. The description of the system is therefore reduced to a simpler one-dimensional problem.
The modeling domain extends between the Stern layers located at the two boundaries, positioned at z = 0 and z = L, as shown in Figure 4. The actual pore surfaces are located at z = −s and z = L +s, although because s L the surface separation is henceforth referred to as L.
We draw the example systems from two processes which play an important role in the behavior of a LiS cell: the precipitation/ dissolution of a species into the solvent and an electrochemical reaction. The first of these represents both the first and last step of the charge and discharge processes, and is potentially responsible for a number of features in the voltage profile, including the dip in the voltage between the two plateaus in the discharge curve 4 and the shape of the flat second plateau. The second process is fundamental to current flow in a LiS cell and many other electrochemical devices.
Although we are motivated by the processes occurring in a LiS cell, we are interested in a more general sense in how the behavior of the processes alters when species interactions are introduced or when multiple processes occur simultaneously. For this reason we do not phrase the test cases in the context of the species present in a LiS cell specifically, but will make reference to how they may impact this particular system. By using a general system, we are able to consider the trends in behavior with the variation in what would otherwise be reaction-or species-specific properties, for example. It also allows us to bypass a common problem with complex electrochemical systems, which is that many parameters are unknown and unmeasureable, because their values cannot be measured within the system of interest.
We assume the temperature to be 298 K and the relative permittivity to be 78, the value for water. Unless otherwise stated, the electrode separation is 40 μm, which is in the range of the typical separation of the electrodes in a LiS or Li-ion cell. In terms of the species properties, we assume that the diffusion coefficients of all species are 10 −9 m 2 s −1 , a typical order of magnitude estimate for dissolved species, and we fix the transfer coefficient in all electrochemical reactions to 0.5.
The model was solved numerically using the COMSOL Multiphysics software package version 5.2a, which uses the finite element method.
Precipitation/dissolution.-The simplest process which the model describes is a dissolution reaction, whereby a solid phase dissolves into a solvent. In the example of a LiS cell, this process occurs as both the first and last steps of the reaction mechanism, and is implicated in causing some of the features of the discharge curves. Although it is conceptually straightforward, the inclusion of the EVI alters some of the dissolution behavior relative to the ideal solution case, and this may lead to changes in a more complex system.
We assume that a general precipitate species A (s) is located on the model boundary at z = 0 in the slit-pore structure and that it participates in the following reaction: We define a baseline case in which the domain length is 40 μm, the initial solute concentration 0.5 M, and the hard-sphere radius of the solvated solute particle 0.32 nm. The properties of the precipitate are listed in Table II, together with the reaction details. The only process which can occur in this system is for the precipitate to either dissolve in to the solvent or for the solute to precipitate onto the surface. As illustrated in Figure 5, it is the former which occurs: with time, the radius of the unit precipitate decreases and the average concentration of the solute increases. Eventually, both equilibrate at new values, the specific values of which are defined by the rate constants and the EVI.
Since the solid phase always has unit activity, the ratio of the forward and backwards rate constants determines the equilibrium solute  activity, i.e. they determine the solubility of the solute in the solvent. In this case, the ratio is 1, and so the equilibrium activity must have unit value. Figure 5 shows the activity at z = 0 and z = L, with both converging to unity as expected. The lag between the two arises because the species entering the solvent at z = 0 need to diffuse to across the domain before they affect the activity at the z = L boundary, but the fact that they ultimately attain the same activity indicates that there are no net concentration (activity) gradients at equilibrium.
The extent to which a solvent is able to dissolve a fixed amount of precipitate depends on the number of free solvent molecules in the system, which is affected by both the the domain length and the concentration of solute, as shown in Figures 6 and 7.
Increasing the domain length provides more free molecules of solvent per molecule of dissolved solid, so the activity of the solute phase increases less rapidly as the solid phase dissolves. A solvent in a longer domain can therefore dissolve more precipitate before becoming saturated. In the model system, for domain lengths L > 65 μm all of the precipitate can be accommodated in the solvent, and the precipitate dissolves completely. The equilibrium activity also drops at this point because the amount of dissolved solute remains constant while the number of free solvent molecules grows, making the final solution more dilute.
Increasing the initial concentration of the solute means reducing the initial number of free solvent molecules, which reduces the capac- ity of the solvent to dissolve the solid. Effectively, the solute already in the solvent has a higher activity, and, therefore, the system is closer to its saturated state. When the concentration is low enough, all of the precipitate is able to dissolve into the solvent, as can be seen from the final radius of the precipitate being zero for c 0 A < 0.4 M in Figure 7. One difference when the initial concentration of the solute is varied is that the initial solvent may be supersaturated, in which case it precipitates out onto the surface. The point at which this occurs is marked, and increases to the concentration above this point lead to the radius of the unit precipitate at equilibrium being larger than its initial value.
While the behavior of the preceding cases is identical to what would be predicted by a model based on ideal solution theory, in such a model the precipitation reaction would not show any dependence upon the concentrations of other species in the system. To examine this effect, we now consider the addition of a completely dissociating salt BC to the solvent; the solution will then contain a mixture of species A (d) as well as the monovalent ions B − and C + . All species have a radius of 0.32 nm.
In Figure 8, we show how the presence of the ionic species changes the equilibrium of the precipitation reaction for a system with an initial A (d) concentration of 0.5 M and a domain length of 40 μm. Adding ions to the system increases the activity of the solute molecules, meaning that the dissolution process ends sooner, so the solute concentration is lower and the precipitate radius larger at equilibrium. As with increasing the initial concentration of solute, there is a threshold in   Table III. the ion concentration above which the activity of the solute becomes greater than one, at which point the ionic species cause the solute to precipitate out of the solution, as marked in the figure. Regardless of the ion concentration, however, the solute activity at equilibrium always has unit value, as it must, unless all of the precipitate dissolves while the system is still able to hold more solute.
The example system we have considered so far is simple, but in a more complex reaction mechanism this effect could begin to play a role in the overall system behavior. In a model which treats the solution as ideal, a given reaction is only indirectly affected by any species not directly involved in the reaction, and only then if they are actually connected as part of a reaction chain. By accounting for species interactions, however, all species in a system will affect all reactions, whether or not they directly participate in them.
The addition of ions to the system alters the equilibrium state, but the interaction between the dissolving solute molecules and the ions can also lead to a diffusiophoresis-like process 28 occurring in the system, causing the development of a transient voltage during the dissolution process. The dissolution of the solute creates an activity gradient in in the solute species. As well as driving the diffusion of the solute across the domain, this also induces an activity gradient in the ionic species because of the EVI. This causes the ions to become temporarily displaced from the surface (or to move toward it if precipitation occurs). As the precipitation process equilibrates and diffusion drives the solute activity to become spatially constant, the force causing the displacement of the ions is removed and they too return to a homogeneous configuration. If the ions have different sizes, the larger ions are displaced more than the smaller ones, causing a local charge density to develop in the electrolyte and a voltage to develop across the system. This effect is illustrated in Figure 9.
The figure shows results for two cases of three different systems. The difference between the cases is the concentration: cases A have an ion (salt) concentration of 1 M, while in cases B it is 0.1 M. Thus we can see from Figure 8 that precipitation occurs in cases A and dissolution in cases B. The case numbers relate to the changes in the domain length and cation size, as summarized in Table III. In all cases the neutral species radius is 0.32 nm and that of the anion is 0.3 nm.  The direction of the precipitation reaction alters the sign of the potential: if dissolution occurs, the potential is negative, while for precipitation it is positive. As the cation is larger than the anion, the negative activity gradient created during dissolution causes the cation to be displaced more than the anion, making the region near z = 0 negatively charged and the region further away positively charged. Conversely, the positive activity gradient caused by the precipitation process means that the cations drift closer to the surface than the anions, and so positive charge develops there.
Since the relative sizes of the species determines how strongly they respond to the activity gradient, the charge separation increases with the difference in ion size. Additionally, the length of the domain affects the quantity of species which has to precipitate or dissolve before equilibrium is reached, which affects the duration of the transient voltage. Finally, the inset shows a sharp change in the gradient of the transient in the low concentration/long domain case, which is associated with the loss of all precipitate from the surface. The sudden change occurs because the loss of precipitate means that the solute near the surface stops being replenished as it diffuses away. Because of this, the solute concentration drops rapidly, decreasing the activity gradient and allowing the ionic species to drift closer together, decreasing the local charge density and thereby the potential. Regardless of the system, the voltage is only temporary, because the activity gradient in the solute will always disappear as the reaction reaches equilibrium.
By building a dynamical model for a charged electrolyte system which accounts for excluded volume interactions and coupling this to a simple description of precipitation, we are able to probe how these components interact. While the results are not necessarily unexpected, they do begin to indicate the limitations of applying ideal solution theory to electrochemical systems, especially as they become more complex.
In the case of a LiS cell, the formation of precipitate plays an important role in the system behavior: the dissolution of a solid phase is required at the start of the discharge process and the precipitation of a solid occurs at the end. As will be discussed shortly, these processes can impact the electrochemical behavior of a model cell, but the state of the cell, in terms of the concentration of the electrolyte, can be seen to affect the precipitation process itself.
One particular problem in LiS cells is the formation of Li 2 S, which on the one hand is thought to be responsible for the transition dip between the two parts of the discharge curve and to contribute to the flatness of the second plateau, but on the other hand passivates the cathode surface, leading to problems with recharging the cell. By better understanding how the electrolyte environment affects the precipitation process, it may be possible to modify the behavior of the precipitation to mitigate some of these problems.
One electrochemical reaction.-The second simplest case described by the model is that of an electrochemical reaction occurring at the surface. This is a fundamental process for many electrochemical systems, and many of the elementary reaction steps which take place as part of the LiS mechanism are of this type. To develop an understanding of how this process behaves within the model, and how species interactions influence its behavior, we assume that a single electrochemical reaction is able to take place on the surface at z = 0.
We define a mixed electrolyte formed by the addition to a solvent of the completely dissociating salts AC and BC, to give a mixture of ions A − , B − and C + . Furthermore, we assume the presence of a neutral species A to also be dissolved in the solvent. The species A and A − are those which participate in the following electrochemical reaction: In all examples studied in this section, the initial concentrations of A and A − (the active species) are equal, and we require c 0 C + = c 0 A − + c 0 B − to ensure overall electroneutrality in the initial state. Unless otherwise stated, the species radii are all 0.3 nm.
In the first study, shown in Figure 10, we examine how the cell potential evolves toward equilibrium at four different initial concentrations of A, each with 1 M of supporting electrolyte. Initially, the reaction is out of equilibrium, and A is reduced to A − , causing the electrode to gain a positive charge. This occurs at a faster rate at higher concentrations, both because there are more molecules at the surface able to react and because a larger concentration of active species A can support a larger activity gradient, facilitating quicker mass transport. Although modeled as a single value, the electrochemical potential of a particular species is a distribution around an average value. The consumption of the reactant lowers the average electrochemical potential of the particles, while the formation of product increases its average electrochemical potential. At the same time, the change in the electrode charge alters the energy required by the molecular species to react. As the reaction proceeds, the shrinking product energy and growing electron energy tend to make the forward reaction slow while the reverse process speeds up, because the electrochemical potential of the product species increases. The reaction therefore reaches a dynamic equilibrium when both rates become equal.
Because the system is closed, if the concentration of active species is lower, there are a smaller number of particles with enough electrochemical potential to react, so their consumption has a larger effect on the average electrochemical potential of the species. Similarly, the formation of the product has a larger effect on increasing its average electrochemical potential when the product concentration is low. As a result, the reaction reaches equilibrium with a smaller number of molecules having reacted, causing the charge on the electrode, and therefore its potential, to be lower when equilibrium is reached.
When the reaction equilibrates, the difference in the electrochemical potentials of the reactant and product must equal the electrochemical potential of the electrons in the electrode. Combined with the knowledge that the electrochemical potentials of the species must be spatially constant at equilibrium, this means that the difference in the average electrochemical potential of reactant and product species far from the surface (which is simply their activity) must also equal the electrochemical potential of the electrons in the electrode. This is the essence of the Nernst equation, linking the species activities far from the electrode to the deviation from the standard state reduction potential of the reaction. For a closed system, however, the Nernst potential is not a constant: in the initial state of the example system, the activities of the reactant and product are equal (their concentrations and ion sizes are the same) and so the Nernst potential equals the standard state reduction potential. However, as the reaction proceeds, the Nernst potential decreases because of the relative changes in the species activities. By using the species activities at the z = L boundary, we approximate the instantaneous Nernst potential, shown also in Figure 10. In all cases, the Nernst potential decreases as the reaction proceeds, although the decrease only really becomes significant at sub-millimolar concentrations, as shown more clearly in Figure 11. As well as illustrating how the Nernst potential changes with the initial species concentration, it also indicates that the Nernst potential varies with the supporting electrolyte concentration, tending to be lower as this increases. Since the Nernst potential depends upon the active species activities, which depend on the species interactions, this is unsurprising, but does indicate a further deviation from ideal solution theory which may become important in concentrated reaction systems.
The effect of the supporting electrolyte concentration and the EVI is explored further in Figure 12. Regardless of the relative sizes of the active species, the supporting electrolyte acts to suppress the equilibrium potential even at low concentrations, causing an initially rapid decrease followed by a roughly linear region if the ions have the same size. This type of behavior has been observed in the reduction of ferrocine in a supporting ionic liquid 29 and has previously been attributed to ion-pair formation in the electrolyte as well as to changes in the solvation energy of the reactive species. 30 The change in reduction potential arises from the change in the species activities. Although this change is the same for all species (assuming the species all have the same size), a growth in the EVI amplifies differences in the active species chemical potentials caused by changes in their concentrations due to the reaction (see Equation 49). Because of this, the reaction ends sooner than it would in an ideal system. Figure 13 shows how the activity of the reactant species decreases for low supporting electrolyte concentrations, while the activity of the product increases. However, as the overall supporting electrolyte concentration increases, the reactant species equilibrium activity does not continue to decrease. As the variation in the potential enters the linear region, the equilibrium activity begins to grow, driven by the increased EVI contribution that results from the high supporting electrolyte concentration.
Where the model begins to deviate from the previous studies is when there is an asymmetry in the active species radii. For both of the active species concentrations, a second curve in Figure 12 indicates the predicted behavior if the anion radius is increased to 0.34 nm. This causes its electrochemical potential to grow more quickly with the increased supporting electrolyte concentration, limiting how far the reaction can proceed before the reactant and product activities become equal. This limits the amount of charge transferred to the electrode causing the equilibrium potential to reduce.
One interesting feature of this model is the prediction that asymmetries in the sizes of the ions should lead to non-linear behavior at very high supporting ion concentrations. For these conditions, the EVI makes the activity of the larger anionic species (in this case) much larger than that of the neutral species, so the reaction does not proceed as far before it reaches equilibrium, meaning that less charge is transferred to the electrode and the resulting potential is lower. Since the growth of the EVI energy is highly non-linear with the concentration, the extent of its effect on electrode charging grows rapidly at high concentrations. This effect does not contradict the Nernst potential, because the differing responses of the two species to the EVI alters the species activities, in turn altering the Nernst potential.
We finally consider the effect of the domain length on the reaction equilibrium. Shown in Figure 14 is the variation of the equilibrium potential with L at a selection of initial active species and supporting electrolyte concentrations, under the assumption that all species are equally sized. The general trend is for the equilibrium potential to be severely limited at short domain lengths but to converge on the standard state Nernst potential as the domain length increases. The reason for this is similar to that for effect of the domain length on the precipitation reaction: the reaction requires a finite amount of reactant to be consumed to generate the requisite electrode potential to put the reaction in equilibrium. As the domain length grows, the reaction removes an ever-decreasing fraction of the reactant from the electrolyte and increases the product concentration by an ever-decreasing fraction of its initial value. The equilibrium potential therefore converges on the Nernst potential for the initial state. It was seen in Figure 11 that the initial active species concentrations alter the Nernst potential, but the extent to which this is true depends on the size of the domain, or the actual number of molecules of the active species. Essentially, increasing the domain length lessens the impact of a low active species concentration. This is only true when the active species are symmetric, however, as shown in Figure 15. By assuming that the anionic reaction product is larger than the neutral reactant, the equilibrium potential is shown to be suppressed at all domain lengths. Furthermore, the trend for the equilibrium potential to grow with the initial active species concentration is broken in longer domains at larger supporting electrolyte concentrations.
At lower domain lengths, the changes to the quantities of the species in the system seem to dominate the equilibrium potential, leading to larger decreases in its value with the initial concentration. At longer domain lengths, when the changes to the concentration are relatively much smaller, the EVI energy begins to dominate, causing the potential to be lower even though the initial solute concentration is higher, contrary to what was seen in Figure 11. The extent of the reduction in the potential still depends on the concentration of the supporting electrolyte, however, being smaller for all domain lengths at lower concentrations, in agreement with the data shown in Figure 12.
Coupled precipitation-electrochemical reactions.-The previous sections show that the species interactions and properties of the system interact to alter the equilibrium state and how the system reaches that state, marked by deviations from the expected behavior of an ideal system. To continue this theme, we now briefly discuss the effect on the equilibrium state of coupling reaction types within the model.
To do this, we define a system in which a precipitation reaction is coupled to an electrochemical step via the neutral species A according to reactions 60 and 61. The precipitate properties are as listed in Table II and the equilibrium rate constant for the electrochemical reaction is K = 10 −8 .
To define the initial state, we again assume the dissolution of the salts AC and BC to form an electrolyte with 0.1 M A + ions, 0.5 M B − ions and 0.6 M C + ions, and that the 0.1 M the solute phase A (d) is also present. The domain length is set to 40 μm and all species have a radius of 0.32 nm.
In Figure 16 we show the time-evolution of the potential from the system's initial state until equilibrium is reached, together with the instantaneous Nernst potential and the radius of the precipitate species. The case with k f = k r = 0 is almost identical to the system shown in Figure 14, except that the presence of the precipitate on the surface slows the electrochemical reaction rate. By comparison, enabling the chemical reaction can be observed to alter both the equilibrium potential and the reaction rate, depending on the solubility of the solute.
When the ratio k f /k r is low (blue lines), the solubility of the solute species A is lower than its current activity, causing it to precipitate at the surface. This consumption lowers the species' surface activity, lowering the rate of the electrochemical reaction and thereby causing the electrode potential to rise more slowly. Moving onto the equilibrium state, we know that the electrochemical potentials of all species must be spatially invariant and that the reactions must be in equilibrium. Because the chemical process defines the activity of the neutral species, which in turn affects the Nernst potential of the electrochem-ical reaction, the dissolution reaction also affects the equilibrium potential of the cell. In this case, because the equilibrium activity of the solute is reduced compared to the case with no chemical reaction, the equilibrium potential is reduced.
Increasing the solubility of the solute (green lines) has the reverse effect: the dissolution of the precipitate increases the activity of A in the vicinity of the electrode, causing the electrochemical reaction to run faster and the electrode to charge quicker. Similarly, the cell moves to equilibrate at a larger potential than when there is no chemical reaction, because the solute activity is higher than it otherwise would be, pushing equilibrium toward to the formation of the A − anion and a larger electrode charge. There is a limit to the dissolution of the precipitate, however, and as it all is consumed from the surface there is a sharp drop in the electrode potential, followed by a more gradual decrease. As the precipitate runs out, the solute molecules diffusing out from the EDL are no longer replaced and so they are displaced by counter-ions being attracted to the surface. The sudden change in the surface activity reverses the electrochemical reaction, with electrons being added to the surface and the potential decreasing slightly. This continues until the electrochemical reaction reaches equilibrium, which now lies further toward the formation of the ionic species because of the increased quantity of the neutral reactant in the bulk.
By describing the evolution of the system from a well defined initial state, the model is able to estimate the equilibrium potential of an electrochemical system driven by spontaneous chemical and electrochemical processes while accounting for species interactions. The approach incorporates the influence of the electrical double layer and electrochemical reactions on each other. [31][32][33] This aspect of the model makes it increasingly useful when trying to describe complex electrochemical systems wherein the interplay of species interactions causes deviations from the standard state Nernst potential. This is reflected, for example, in the instability of homogenized LiS models, wherein small changes to the initial conditions can stop the model from functioning. 7 This is not to say that this model does not have its numerical complications. For example, the model predicts large changes in the gradients of the dependent variables in the EDL, which is a source of numerical error in the solution. 34 Because the surface flux boundary conditions are intimately linked to the solution at the surface, this error interacts with the boundary conditions contributing to instability in the model. However, these numerical difficulties can be managed through suitable choice of the mesh and solver, as opposed to requiring a brute force approach to determine a suitable initial state.

Cycling a cell with a coupled precipitation-electrochemical reaction.-
The ultimate goal of the model is to be able to understand how complex electrochemical systems behave under transient conditions. While we have so far seen that interactions between processes alter the way that the system develops its initial equilibrium potential, we now take a short look at how they can also lead to significant changes in, for example, a voltage charge-discharge profile of a model electrochemical cell.
To do this, we complete the model cell by placing a second electrode at the z = L boundary and allowing an electrochemical reaction to occur at its surface. In common with a LiS cell, this takes the form of an electrodissolution reaction, whereby an ionic species C + is released or absorbed at the electrode surface, releasing or consuming an electron in the process. Since the solid phase at this electrode is assumed to be electrically conducting, its formation does not affect the electrode's active surface area. Furthermore, as with the precipitation reaction, the addition or removal of the ionic species is assumed not to alter the separation between the electrodes. The reactions at the z = 0 boundary are the same as in the previous example. All reactions used in this section are summarized together with their properties in Table IV.
The addition of the second electrode means that a galvanostatic current can also be applied according to the boundary condition in Equation 57. A discharge current neutralizes the charge on each electrode, reducing the electrochemical potential of the electrons therein.

Reaction
Position This allows the electrochemical reaction to recommence, the rate of which depends on the activities of the active species and their transport properties.
The loss of charge from the electrode due to the current causes the electrolyte to shift from its equilibrium state, causing the reactiondiffusion processes to recommence in such a way as to replenish the electrode charge. Because the current continually acts to decrease the charge, if the reactions do not occur quickly enough to replenish it then electrode potential will drop, placing the system further from equilibrium. However, increasing the deviation from equilibrium also causes the reaction rates to increase, until the electrode charge is replenished at the same rate as it is consumed. For a given external current, therefore, there will always be a loss in the cell potential, but this ensures that the internal current matches the external current, ensuring the law of conservation of current is observed. The extent to which the system has to move from its equilibrium position in order to maintain the current depends on the properties of the reactiondiffusion process, but is the essence of the overpotential or polarization of a discharge curve Since the parameter space for discussing this hypothetical system is quite large for this system, we limit this section to two data sets which give an indication of how strongly the processes are able to affect the discharge profile of this model cell. In both data sets, the system comprises four types of mobile particles: A, A − , B − and C + , each with a radius of 0.32 nm and involved in the same reaction processes. The initial concentrations, reaction properties and cycling conditions for the two data sets are listed in Table V.
For the first set we compare the effect of the quantity of precipitate which is able to dissolve on the discharge profile. Three cases are considered, with differences in the initial quantity of precipitate. The differences between cases 2 and 3 are listed in Table V; case 1 is the same as case 2 but with the precipitation reaction disabled, meaning that the effect of the precipitate on the electrochemically active surface area is accounted for, but that the dissolution itself is disabled.
In all cases the discharge/charge curves are monotonic, with the second and subsequent (not shown) profiles featuring an initial rapid change in the voltage, a more gradual middle section and another rapid change at the end. The first discharge looks like the second half of a full discharge, which it essentially is: because some of the reactants are consumed to give the cell its initial equilibrium voltage, the amount of active material available for the first discharge is smaller than in  Figure 17. Effect of precipitation reaction and precipitate quantity on the cycling voltage (upper plot) and available electrochemically active surface area (lower plot). In case 1, precipitate is present but the reaction is disabled; cases 2 and 3 both have the reaction enabled but with differing initial quantities of precipitate (see Table V).
subsequent discharges. There is also a symmetry between the charge and discharge profiles. This shape of voltage curve is characteristic of when the discharge process is limited by the electrochemical step. The finite quantity of neutral species A decreases as the cell discharges, and the corresponding decrease of the activity in the domain and at the surface causes the diffusion and reaction processes to slow. To counter the loss of the active species, the system moves further from equilibrium in order to maintain the internal current, corresponding to a continual reduction in the cell potential. It should be noted that the conversion of electrochemical reactant to product on discharge also causes the Nernst potential of the reaction to decrease with time. Because of this, not all of the voltage change observed is due to the increased deviation from an equilibrium state -some is attributable to the equilibrium state itself changing.
The effect of the precipitation reaction on the system is most clearly visible in the duration of the discharges which, because the current is constant, translates to the capacity of the model cell. The cell capacity is the measure of how much charge can be passed from the electrodes into the electrolyte, which depends upon the availability of species to carry the charge. By providing a source of the active neutral species through the precipitation reaction, the amount of charge that the electrolyte can hold is increased, thereby increasing the cell capacity.
The lower plot of Figure 17 shows the fraction of the electrode surface available for electrochemical reactions, which is a function of the quantity of precipitate on the surface -if θ S = 1, there is no precipitate on the surface. The plot thus shows that cycling the cell cycles the amount of precipitate, even though the chemical reaction is not directly affected by the current. Furthermore, similar to the features seen in Figure 16, and for the same reasons, the termination or onset of the precipitation reaction alters the voltage profile, as shown for the final discharge of case 3 in the inset of Figure 17.
For the second set of data in this section, shown in Figure 18, we demonstrate that complex behavior can emerge from this relatively simple model, depending on the rates of the processes. The initial conditions are indicated in Table V, where the difference between the two cases is the rest period between consecutive charge and discharge phases. While the initial conditions are similar to those of the first data set, in this case the solubility of the neutral species is larger but the rate of the precipitation reaction is significantly lower. We also apply a larger current and compare the case in which there no rest between charge and discharge phases to that where there is a 30 s rest.  Figure 18. Restricting the rate of the precipitation reaction while increasing the current significantly alters the shape of the discharge curve (see Table V for system details). Furthermore, resting the cell between cycles, shown by cases 1 (no rest) and 2 (30 s rest) has an impact on the subsequent discharge).
Changing the chemical rate constants but removing the rest period has the effect of flattening the middle part of the discharge profile and breaking the symmetry between the charge and discharge. The chemical reaction introduces a bottleneck to the formation of species A, which then limits the rate of the electrochemical reaction. Because this is slowed, the electrochemical reaction consumes as much of species A as is required for the chemical reaction rate to increase until it replaces species A at a similar rate to that at which it is electrochemically consumed. In this way the voltage profile stabilizes, but at a lower voltage.
The effect on the model cell of resting it is to allow the system to fully equilibrate after each charge and discharge phase. This allows for some recovery in the voltage as the electrochemical reactions are no longer competing with the current to charge the electrodes, so they replenish the charge up to the limit of the reaction equilibrium. In terms of the chemical reaction, all of the precipitate has been lost during discharge, meaning that resting the cell at this point has little effect on the state of the cell, and so the charges are similar regardless of the rest period. On charge, the solute continues to precipitate onto the surface during the rest, which changes the composition of the electrolyte and the coverage of the electrochemically active surface area, leading to differences between the discharge profiles of the rested and unrested cases. In particular, following a rest, the voltage profile can be seen to be non-monotonic, with a sharp initial decrease in the voltage being followed by a humped voltage profile. This ties into the notion that precipitation is required in order to explain the shape of a LiS discharge curve.

Conclusions
Starting from a description of the Helmholtz free energy of a system, we have outlined a framework for describing a general reactiondiffusion system. The species fluxes and the rates of the reactions are defined by gradients or differences in the species' electrochemical potentials, which follows from the form of the chosen form of the free energy functional, and we have shown that the model satisfies the constraint that the free energy is minimized with time by the processes occurring within it.
Working within the mean-field and local density approximations, we have developed this framework into a model describing the transient behavior of a system in which chemical and electrochemical processes are able to occur. We demonstrated this in a structurally simple slit-pore model, with which we explore the behavior of some example cases drawn from the reaction mechanism driving a LiS cell. While this particular system provides our motivation, the structure of the model can be applied to any reaction-diffusion system in which it might be desirable to understand the behavior of the system in more detail than can be provided by a homogenized electroneutral model.
We first applied this model to the dissolution of a surface precipitate into a solvent or electrolyte, showing how the inclusion of the species interaction term alters the precipitation/dissolution behavior for a given reaction. The reason for this is that the interactions alter the activities of the solute species, effectively altering their solubility. Because the LiS reaction mechanism is bounded by two such chemical processes, each occurring into a concentrated solution, the solute phases will inevitably be affected by these interactions, ultimately having an impact on the behavior of the cell. We also showed how the model predicts a small transient voltage during the dissolution process if the solvent also holds asymmetrically-sized ions.
Next, we used the model to describe the case of a single electrochemical reaction occurring at one of the domain boundaries. By considering changes in the initial active species concentration and domain length we observed deviations in the equilibrium potential of the reaction. These expected differences are caused by changes in the relative species activities caused by the consumption of reactant and formation of product in a closed system. The deviation of the equilibrium potential fits exactly with the change in the Nernst potential as the reaction progresses. Further deviations in the equilibrium potential are observed as a function of the supporting electrolyte concentration and the relative sizes of the electrochemically active species, becoming significant at high supporting electrolyte concentrations, even when the initial active species concentration is high. In the context of a LiS cell or many other electrochemical devices, wherein the electrolyte is highly concentrated, this variation indicates that the Nernst potential used in a homogenized cell model may be significantly inaccurate.
In the final two sections we considered how the equilibrium state and the dynamics of the system are affected by coupling a precipitation step to an electrochemical step. In terms of the equilibrium, the existence of the precipitation reaction was shown to both alter the rate at which the system equilibrates and the cell potential once equilibrium is reached. Furthermore, the shape of the voltage profile during the equilibration stage was observed to be strongly affected by the complete dissolution of the solid phase from the surface. Regardless of these effects, however, the electrochemical reaction still satisfies the Nernst potential at equilibrium, provided the influence of the solute phase on the species activities is accounted for. This result complicates the understanding of LiS behavior, Upon applying a current to the model cell, the presence of the solid phase precipitate was shown to increase the duration of the discharge, indicating an increase in the capacity of the model cell. The loss of the solid phase from the surface or the onset of precipitation gave rise to small features in the charge and discharge profiles. Additionally, severely restricting the rate of the chemical reaction while increasing the current led to a complete shift in the shape of the voltage profile. The existence of a chemical reaction bottleneck in the LiS reaction mechanism is thought to be responsible for the transition between the two plateaus seen in a LiS discharge profile. 4 That the emergent behavior of even a very small number of processes exhibits a reasonable degree of complexity suggests that gaining a true understanding of the mechanism may be a challenge, but the model provides a starting point for doing this.
In addition to allowing for the inclusion of species interactions, the formulation of the model provides a basic structure for studying the behavior of higher-dimensional systems, potentially providing the means to understand how reaction-diffusion processes may occur in confined geometries. As mentioned previously, it is known that the species interactions can lead to the exclusion of larger particles from the EDL which can be related to ion selectivity in porous structures. By also including the reaction processes, this model makes it possible to begin exploring how pore size and shape can alter a reaction mechanism.
A number of avenues are available for the general future improvement of the model within the framework. These may include the incorporation of additional species interactions or improved expressions for the reaction rates, but also extend to more fundamental developments. Examples are the inclusion of the heat of reaction for the chemical processes, entailing the coupling of a thermal model to the system, or a better description of the precipitation reaction. In the latter case this might include either a more realistic nucleation and growth model or, in the case of a two-or three-dimensional model, a physical description of the precipitate in the model domain.
The precipitation model in particular places a physical limit on the maximum amount of solid phase which can be present before the electrochemically active surface becomes blocked. This limits the one-dimensional model's capability to describe a porous electrode structure of a LiS or similar cell in which the areal precipitate loading may be larger than the model can describe.
Finally, while we are motivated by the desire to understand the LiS reaction mechanism, the model is generally applicable to a range of chemical and electrochemical systems, and may find use in helping to understand experimental data from a range of tests. One particular example may be electrochemical impedance spectroscopy data in weakly supported electrolytes, for which the equivalent circuit models used to fit the data may break down because of changes in the electrolyte resistance with concentration.