Modelling the ion-exchange equilibrium in nanoporous materials

Distribution of a two component electrolyte mixture between the model adsorbent and a bulk aqueous electrolyte solution was studied using the replica Ornstein-Zernike theory and the grand canonical Monte Carlo method. The electrolyte components were modelled to mimic the HCl/NaCl and HCl/CaCl_2 mixtures, respectively. The matrix, invaded by the primitive model electrolyte mixture, was formed from monovalent negatively charged spherical obstacles. The solution was treated as a continuous dielectric with the properties of pure water. Comparison of the pair distribution functions (obtained by the two methods) between the various ionic species indicated a good agreement between the replica Ornstein-Zernike results and machine calculations. Among thermodynamic properties, the mean activity coefficient of the invaded electrolyte components was calculated. Simple model for the ion-exchange resin was proposed. The selectivity calculations yielded qualitative agreement with the following experimental observations: (i) selectivity increases with the increasing capacity of the adsorbent (matrix concentration), (ii) the adsorbent is more selective for the ion having higher charge density if its fraction in mixture is smaller.


Introduction
A variety of phenomena of interest for basic science and technology occur when an electrolyte solution or a mixture of electrolytes is in equilibrium with a porous phase containing fixed charges. Comprehension of the mechanism governing the equilibrium distribution of ions between the two phases is therefore of considerable interest. Materials containing charged micropores, such as membranes, glasses and certain clays, have a property of partially rejecting an electrolyte being filtered through them. In addition, such materials may be ion-selective; the ions having different sizes or charges adsorb to a different amount. This makes a basis for the separation process used in water softening and purification, in treatment of radioactive wastes, and for analytical purposes as well. Ion-exchange is, therefore, not only a powerful tool in chemical analysis and in scientific research, but also of great importance for everyday life [1][2][3][4].
For a quantitative description of these phenomena, some knowledge of the spatial distribution of the counter-ion and co-ion species in the micropores is needed. To begin a mathematical description, a suitable model of the adsorbing material should first be assumed. The ion-exchange resins are quite often modelled as spherical holes or cylindrical channels. In such cases the calculations are, due to the known geometry, relatively simple (see, for example, [5]) and the Poisson-Boltzmann theory is most often accurate enough to make viable predictions of measurable properties. Alternatively, the Monte Carlo or Brownian Dynamics methods may be used for this purpose [6][7][8][9][10].
In contrast to the models with well defined (spherical, cylindrical, etc.) geometry described above, we are here interested in modelling the porous systems where distribution of charges in a system is not governed by a simple geometry. The porous material (matrix) is pictured as a set of charged obstacles. Distribution of the obstacles is obtained by a certain recipe (see, for example, [9,11]). Inhomogeneity that causes the field in which mobile ions are distributed is given on a molecular level. Such a model is a representative of a quenched-annealed systems which have been more extensively studied only in the last 25 years. The progress in this area of research is documented in several review papers [12][13][14][15].
Following these pioneering studies, Pizio and coworkers [33] succeeded in extending the replica theory to the systems with Coulomb interactions. In a series of papers [34][35][36][37][38][39] the quenched system was modelled as a set of charged obstacles containing positive and negative charges so that the matrix subsystem (the set of obstacles) was electroneutral. The other electroneutral subsystem, that is electrolyte solution, was assumed to occupy the void space in-between the obstacles, i.e., to invade the matrix. To calculate the thermodynamic properties and spatial distribution between the quenched and annealed particles the replica Ornstein-Zernike (ROZ) integral equations were applied. In addition, an expression was derived [38] for the mean activity coefficient valid within the replica hypernetted-chain approximation. The canonical and grand canonical Monte Carlo method was used to verify the theoretical approximations inherent to the replica theories. Very good agreement between the theory and simulations was obtained for both, thermodynamic parameters and for the pair distribution functions (see, for example, [37][38][39]). More recently, good agreement was confirmed by the independent Brownian Dynamics simulations based on the same model [9].
In the last years, the model for a quenched system was modified to allow studying the ion-exchange phenomena. The model proposed in references [11,40,41] differs from the one presented above in one important aspect. We assume again that the distribution of charged obstacles corresponds to an equilibrium distribution of ions in an electrolyte solution under conditions given by the temperature and the dielectric constant before quenching. However, in contrast to the previous studies, we assume that only anions are quenched and the cations are allowed to equilibrate (anneal). We also assume that the model system is in thermodynamic equilibrium with an external electrolyte solution, which penetrates into the matrix. A similar model, but for a somewhat simpler system, was previously studied by Pastore et al. [42]. The model was recently studied by the grand canonical Monte Carlo method [11,40] as well as the Brownian Dynamics approach [10]. The simulations indicate that the ROZ theory in the hypernetted-chain approximation, adapted to this model in [11], yields very good agreement with computer simulations.
In the present article, we extend the approach of the previous studies [11,40] to the investigation of the ion-selectivity effect in the matrix representing a charged nanoporous system. The system of quenched obstacles formed from the primitive model ions is negatively charged. Within this system (matrix) the annealed ions are distributed; an excess of cations is present to maintain the electroneutrality condition. In equilibrium, the concentration of the invading electrolyte is determined by the properties of the bulk solution of the same chemical composition. Such a system can represent a crude model of the ion-exchange resin. The distribution of ionic species within the matrix is calculated by the ROZ integral equations and in several cases also by the grand canonical Monte Carlo (GCMC) method. For a certain composition of the external electrolyte, a different equilibrium composition of the internal electrolyte solution is established, favouring the ionic species, which interact stronger with the matrix particles. Structural, as well as the thermodynamic properties are calculated, such as the mean activity or Donnan exclusion coefficients, varying the composition of an external electrolyte. The ion-exchange isotherms, reflecting the ion selectivity in HCl/NaCl and HCl/CaCl 2 mixtures, are presented as functions of the model parameters.

The model and theoretical methods
To model an ion-exchange phenomenon, a quenched-annealed system composed of two subsystems was considered. The quenched negatively charged ionic obstacles represented the matrix subsystem, while the annealed subsystem was a mixture of two primitive model electrolytes with the addition of cations needed to compensate the negative charge of the matrix. It was assumed that the matrix remained insensitive to the presence of the annealed fluid, i.e., its structure did not change upon the introduction of the annealed phase. The solvent was treated as a structureless continuum (McMillan-Mayer level of description) characterized by the dielectric constant of pure water under the conditions of the study. Throughout the paper, indices 0 and 1 designate the matrix and the annealed subsystems, respectively.
The matrix phase was prepared in the following way: the +1 : −1 restricted primitive model electrolyte, with the diameters of ions σ 0 + = σ 0 − = 4.25 Å, was quenched at a certain temperature T 0 , where the dielectric constant of the medium was ε 0 . The structure of the matrix was considered to represent one of the equilibrium structures, corresponding to a given set of conditions (concentration, temperature, relative permittivity), and governed by the Coulomb pair interaction potential is the Bjerrum length (e 0 denotes the unit charge, ε v the permittivity of vacuum, and k B the Boltzmann constant), r designates the distance between the centres of ions i and j , and β 0 = 1/(k B T 0 ). In the final stage of the matrix preparation, the cations were disregarded, so that the quenched phase consisted only of a spatially fixed arrangement of negatively charged (z 0 − = −1) particles. Conditions of the matrix preparation corresponded to λ B,0 = 7.14 Å (water solution at 298 K), and ion concentrations c 0 = 0.1, 0.5, 1.0, and 2.0 mol dm −3 .
The annealed subsystem was modelled as a mixture of two +1 : −1 electrolytes, or a mixture of a +1 : −1 and +2 : −1 electrolyte having a common anion. An additional number of univalent cations (common to cations of one of the electrolytes in the mixture), was present to neutralize the overall negative charge of the matrix. The diameters of the +1 ions were σ 1 + = 5.04 Å (model for hydrogen ions, H + ) or σ 1 + = 3.87 Å (model for sodium ions, Na + ), of the +2 ions σ 1 + = 7.03 Å (model for calcium ions, Ca 2+ ), and of the −1 ions σ 1 − = 3.63 Å (model for chloride ions, Cl − ). Models of HCl and NaCl, or HCl and CaCl 2 mixtures were, therefore, considered. Cations compensating the matrix charge were taken to be H + , mimicking in this way a H + -ion-exchange resins. The system was considered to thermally equilibrate under the conditions of T 1 and ε 1 , with the pair potential for a 0 − 1 interaction being equal to and for a 1 − 1 interaction is the Bjerrum length for the annealed solution, and β 1 = 1/(k B T 1 ). Although the conditions for the preparation of the quenched phase can differ from the conditions of observation, we assumed them to be equal in this study, i.e. λ B,1 = λ B,0 = 7.14 Å.

The replica Ornstein-Zernike integral equation theory
A case of adsorption of a single electrolyte within the charged matrix was considered by our group several years ago [11]. Here, we extend the theoretical description to the study of mixtures of electrolytes with a common ion.
In order to define the spatial distribution of matrix ions we first need to obtain the structure for a +1 : −1 electrolyte. This is obtained by solving a set of Ornstein-Zernike (OZ) equations in the form h 00

23802-3
where h and c stand for the total and the direct correlation functions, respectively, and ⊗ denotes the convolution in r -space. Due to spherical symmetry, h 00 +− = h 00 −+ and c 00 +− = c 00 −+ . Equations (2.4) need to be renormalized before being solved numerically (see, for example, [43,44]). The hypernetted-chain (HNC) closure relation between h and c (c i j = h i j − ln h i j + 1 − β 0 U 00 i j ) was used to obtain the solution. Note that only h 00 −− is needed in the subsequent replica OZ equations (2.5). This function contains all the information of the matrix subsystem.
The replica Ornstein-Zernike equations for the quenched-annealed system (cf. [22,41]) read where indices i and j stand for cations or anions of the annealed electrolytes in the mixture (i.e., H + , Na + and Cl − , or H + , Ca 2+ , and Cl − ), while A and B denote only the cations of the annealed mixture (H + and Na + , or H + and Ca 2+ , respectively). We stress that the ρ 1 H + is the total density of H + ions in the system which belong to the HCl and to the cations compensating the negative charge of the matrix.
As in the case of matrix preparation (see expression (2.4)), equations (2.5)-(2.7) need to be renormalized prior to numerical solution (i.e. c and h need to be split into a short and long range part). The procedure is documented in [34,35]. The only difference from the published equations is that here κ 2 contains the sum over all the annealed species. Additionally, HNC closure condition was applied to correlate the h and c functions: The individual activity coefficient within the replica HNC formalism reads [36] ln γ 1 where c (s) denotes the short range part of the direct correlation function coming from the renormalization procedure, and dr = 4πr 2 dr . Again, A and B denote cations in the mixture (H + and Na + , or H + and Ca 2+ ).

Monte Carlo simulations
The matrix was prepared from a size symmetric (σ 0 The annealed electrolyte ions were then distributed within the matrix and the system was studied using the Monte Carlo method in the grand canonical ensemble. The methodology of this approach is well established and extensively described in several previous papers and is not, therefore, repeated here [11,40,41]. 168 The details of simulations are as follows: The number of matrix particles in the simulation box varied from 300 to 2500, depending on the concentration. The average number of fluid cation species distributed within the matrix varied from 500 to 3500. The periodic boundary conditions within the minimum image convention were applied. The ions within the matrix were first equilibrated over at least 10 7 Monte Carlo steps. After the equilibration, a production run of 2 · 10 8 attempted configurations was carried out to obtain the average concentration of the adsorbed electrolyte species. Note that thermodynamic properties were averaged over different annealed fluids, as well as over one to two different matrix configurations. To calculate the equilibrium concentration in the bulk electrolyte, the activity coefficients of the bulk electrolyte mixture should be known in advance. These values were obtained using the hypernettedchain (HNC) theory which has proved to be very successful in describing the properties of ionic fluids [45].

Results and discussion
The two theoretical methods described above were established to be complementary in the sense that they yield consistent results for structural and thermodynamic properties of electrolyte solutions adsorbed in nanoporous materials [11,37,38,40,41]. While the ROZ/HNC theory is computationally less demanding, and is, therefore, faster to use for systematic investigations of different effects on the system properties, the GCMC method is more convenient in studying the equilibrium distribution of ions between the bulk and matrix phase. In this paper we used both methods; the results are organized as follows. First, we tested the performance of the ROZ/HNC theory for the case studied here by comparing the structural (pair distribution functions) and thermodynamic properties (mean activity coefficients, γ ν + +ν − ± = γ ν + + γ ν − − ) obtained using the two methods. Next, the effect of the matrix concentration (e.g. adsorbent capacity) on the mean activity coefficient, which determines the chemical equilibrium in the system, was examined Table 1. The mean activity coefficients of electrolytes adsorbed in the matrix: γ ± , obtained from the

ROZ/HNC theory and GCMC simulation for different matrix concentrations. Concentrations of HCl and
NaCl (c HCl and c NaCl , respectively) and of the matrix (c 0 ) are given in mol dm −3 . I out = 0.5 mol dm −3 , λ B,0 = λ B,1 = 7.14 Å. The numerical error of both methods was estimated to be in the last digit given.

Comparison of the ROZ/HNC results with computer simulations
To test the accuracy of the ROZ/HNC theory for the model of interest, we compared the pair distribution functions and the mean activity coefficients with the results of computer simulations. Figures 1   and 2 show the pair distribution functions (++ and +−) of the annealed ions, while figure 3 shows the annealed ion-matrix distribution functions (+M and −M ). In all cases the concentration of the matrix was c 0 = 1.0 mol dm −3 . Figure 1 and the left panel of figure 3 apply to the mixture of HCl and NaCl, while figure 2 and the right panel of figure 3 describe the distributions in the mixture of HCl and CaCl 2 . Symbols denote the GCMC and the lines denote the ROZ results. Excellent agreement of the two methods is obtained for distribution functions of hydrogen, while the agreement is slightly worse, but still good, for the functions belonging to sodium and calcium ions. Note that the hydrogen ions are the ones compensating for the adsorbent charge and are, therefore, present in an excess to other ions. A similarly good agreement is obtained for the annealed ion-adsorbent pair distribution functions shown in figure 3.   c 0 = 1.0 mol dm −3 . Symbols represent GCMC data while lines show ROZ/HNC theory. λ B,0 = λ B,1 = 7.14 Å. The results for the mean activity coefficients, γ ± , in aqueous mixtures of HCl with NaCl, and HCl with CaCl 2 in the matrix are collected in tables 1 and 2, respectively. Note that in all cases, the particular electrolyte concentrations given, are in equilibrium with the bulk solution of ionic strength I out = 0.5 i z 2 i c out i = 0.5 mol dm −3 . The agreement between the ROZ/HNC and GCMC results is excellent for small matrix concentration, but it deteriorates for larger c 0 values. At c 0 = 1.0 mol dm −3 , the mean activity coefficients calculated within the ROZ/HNC theory are a bit too large, as already noticed in previous papers [11,36,38,40]. Nevertheless, the differences discussed here are within the estimated numerical errors of the two methods.

The effect of the adsorbent capacity on the mean activity coefficient of the invading electrolyte.
The ROZ/HNC method was used to investigate the effect of the adsorbent capacity, given here as the matrix concentration, on the mean activity coefficient of the invading electrolyte. Note that while the ionic strength of the annealed electrolyte mixture within the matrix (denoted by superscript "in") was kept constant, I in = 0.5 mol dm −3 , for all cases examined here, the composition given by   As the matrix concentration (capacity of the adsorbent) increases, the mean activity coefficients of all the species in the solutions increase from values smaller than unity to values higher than unity, as previously observed for a single electrolyte adsorbed in charged matrix [11,40]. The matrix charge repels the co-ions of the invading electrolyte mixture, which causes the exclusion of the electrolyte on the 172 23802-8 whole from the matrix. Concentration of a particular electrolyte is smaller in the matrix than in the corresponding equilibrium bulk solution. This can be clearly seen from the values of the Donnan exclusion coefficient, Γ i , for a particular salt i (HCl, NaCl, or CaCl 2 ). Γ i is defined as The results for the Donnan exclusion coefficients, calculated by the GCMC method (ionic strength of the bulk electrolyte was kept constant, I out = 0.5 mol dm −3 ), are shown in figures 6 (HCl/NaCl) and 7 (HCl/CaCl 2 ). As indicated by the behaviour of the mean activity coefficients (figures 4 and 5), the annealed electrolytes are excluded from the adsorbent (positive Γ values) in all cases. The exclusion is stronger for higher matrix concentration (i.e. higher adsorbent capacity, or higher charge), and depends on the composition of the invading electrolyte. The exclusion coefficient for all electrolytes in the mixture decreases with the increasing fraction of the HCl in the solution. While the Γ HCl approaches the value one as X in HCl → 0 (i.e. for no exclusion), the exclusion coefficient of the other electrolyte in the mixture approaches a somewhat lower value. Such a result has previously been obtained for electroneutral adsorbents [41]. The difference in the behaviour of the two electrolyte components is caused by the modelling of the system (see section 2), mimicking the H + -ion-exchange resin. The matrix charge is namely "neutralized" with hydrogen ions, which are accordingly always present in the annealed electrolyte mixture.

The ion-exchange isotherms
Due to a different degree of adsorption (exclusion), the concentration ratio of the competing ions in the adsorbent (cations) is generally different from that in the bulk solution. In other words, the adsorbent "selects" one species in preference to the other. The study of selectivity was one of the main goals of this work. Here we present the results in the form of ion-exchange isotherms, showing the mole fraction of one electrolyte component within the adsorbent, X in , as a function of the composition of the equilibrium external solution, X out . The ion-exchange isotherms are presented in figure 8 (HCl/NaCl mixture) and figure 9 (HCl/CaCl 2 mixture). Symbols show the grand canonical MC results at different matrix concentrations, and the connecting lines merely serve to guide the eye. The diagonal lines apply to a hypothetical case with no preference to the adsorbent. As noticed before [41,46], the selectivity increases with the increased capacity of the adsorbent (increased matrix concentration). In all cases, the mole fraction of the HCl solution in the adsorbent is smaller than in the external solution (left panels of figures 8 and 9), the opposite is true for the NaCl or CaCl 2 (right panels of figures 8 and 9), respectively. The selectivity strongly depends on the composition of the solution as well. For the component that is preferentially adsorbed in the matrix (NaCl and CaCl 2 in our case), the selectivity increases with an increasing fraction of the component in the solution, while the opposite is true for the other component. Similar trends were observed experimentally [1,5]. Due to a stronger electrostatic interaction between the matrix charges and the ion with higher charge density (Na + or Ca 2+ , respectively), the latter ions are excluded from the matrix to a smaller extent. 174

Conclusions
The model of ion-exchange resin was studied using the ROZ/HNC theory and Monte Carlo simulation in grand canonical ensemble. Theoretical results for the pair distribution functions were found to be in good agreement with computer simulation results; it is confirmed that ROZ theory represents a viable alternative to computer simulations. For all the electrolyte mixtures and matrix concentrations studied here, the electrolyte components are excluded from the adsorbent phase. This is most often explained as a consequence of electrostatic repulsion between matrix charges and co-ions from the electrolyte component. However, the ionic species are adsorbed to a different amount. Consistently with experimental observations, the ions with higher charge density (sodium and calcium ions in our case) adsorb to a greater extent.