|Home | About | Journals | Submit | Contact Us | Français|
A simple but powerful numerical simulation for analyzing the electrochemical behavior of ion-selective membranes and liquid junctions is presented. The computer modeling makes use of a finite-element procedure in the space and time domains, which can be easily processed (e. g., with MS Excel software) without the need for complex mathematical evaluations. It leads to convincing results on the dynamic evolution of concentration profiles, potentials, and fluxes in the studied systems. The treatment accounts for influences of convection, flow, or stirring in the sample solution that act on the boundary diffusion layer and it is even capable of including the effects of an electrolyte flow through the whole system. To minimize the number of arbitrary parameters, interfacial reactions are assumed to be near local equilibrium, and space-charge influences are considered via phase-boundary potential differences. The applicability of the computer simulation is demonstrated for different ion-selective membranes as well as for liquid junctions. The numerical results are in excellent agreement with experimental data.
For the majority of conventional applications, the response of polymeric membrane ion-selective electrodes (ISEs) can be described quantitatively on the basis of a steady-state approach [1-5]. Since the relevant ion-transfer reactions at the interface between the membrane and the aqueous sample are usually very fast, a model assuming a thermodynamically controlled phase-boundary potential is sufficient to explain the emf behavior as a function of the sample composition [4, 5]. Accordingly, a more fundamental analysis of space- and time-dependent processes within the membrane phase is not required for such cases. However, various new techniques and applications emerged recently, for which the understanding of time-dependent phenomena may be essential.
New polymeric membrane materials such as polyacrylates and polymethacrylates have found considerable interest in the past years. In such phases, the diffusion coefficients of dissolved species are by 2–3 orders of magnitude lower than in plasticized poly(vinyl chloride) (PVC), which is the most widely used membrane matrix . When freshly prepared PVC membranes are conditioned in a solution of the respective primary ions, a perfect steady state is reached after only 10–20 hours . In contrast, conditioning times of months to years would be needed in the case of the novel membranes with extremely low diffusion coefficients, which is far beyond the time scale of routine measurements. The response behavior of such ISEs may persistently depend on membrane-internal processes as a function of the contact time with different solutions. Hence, a detailed knowledge of concentration and potential profiles in the space- and time-domains would be helpful for a full understanding.
Very often, the response of ISEs to sample ions at submicromolar levels is also characterized by time-dependent phenomena. The actual behavior is here influenced by ion fluxes, coupled to interfacial ion-exchange processes, and thus depends on the previous history of the ISE membrane [3, 9]. For instance, the shape of calibration curves can be strongly varied just by changing the contact time with the different sample solutions . Again, an appropriate modeling of the underlying processes is required for a quantitative explanation.
In some cases, transmembrane fluxes of primary and interfering ions can even alter the composition of the membrane on its inner side and even the internal solution of the ISE, which results in a long-term variations of the potential at the inner membrane-solution phase boundary . This situation is encountered when the membrane is initially conditioned with interfering ions and then exposed to the primary ions, or vice versa. The effects are more serious when the volume of the internal solution is very small, as for a water film between the membrane and the metal surface of a solid-contact ISE . Similarly, the response of supported liquid-membrane electrodes is also sensitive to phase-boundary influences on both sides of the membrane .
More recently, potentiometric ISE measurements under controlled current conditions have been introduced in order to compensate for zero-current ion fluxes [13-15]. Novel techniques based on current pulses and transient potential recording have also been described for ISEs [16, 17]. It is evident that a current flow generally affects the concentration profiles of all mobile ionic species in the membrane. Unfortunately, explicit theoretical descriptions are not available for such complex cases. Instead of this, computer simulations may be the appropriate way to gain more insight into the response mechanisms and dynamic phenomena of ISEs.
A series of numerical simulations have been reported in the past for membranes of ion-selective electrodes [18-25], for liquid junctions of reference electrodes [22, 24] and for diffusion layers of other electrode systems [26-29]. In an early contribution, Buck’s group introduced digital simulations for investigating the electrical behavior of membranes based on liquid ion-exchangers [18, 19]. The work was later extended to ionophore membranes with mobile ion-exchanger sites . These studies made use of finite-difference (or finite-element) procedures and were based on the conventional Nernst-Planck formalism, assuming electroneutrality to be fulfilled throughout the system. Buck also pioneered a more general treatment that accounts for space-charge influences according to the Poisson equation . Recently, Lewenstam’s group presented a new version of this space-charge approach and reported different applications to membranes and liquid junctions [22, 24, 25].
There is no doubt that the Nernst-Planck-Poisson theory offers the most complete and universal description of membranes and related systems. On the other hand, the observable effects of space-charge regions in real experiments on ISEs are usually very small, apart from fast ac-measurements as applied for the study of complex impedance spectra [21, 25]. The reason is that space charge, as a rule, is restricted to the nanometer domain, in accordance with the respective Debye length . The resulting net charge concentrations are usually well below the range of the total ion concentrations, except for extremely thin membranes [6, 31]. A further limitation of space-charge approaches is their need for highly involved, implicit mathematical procedures, which may become a critical factor [22, 24, 25]. Even more seriously, none of the earlier attempts of simulating ISEs accounted for the presence of a stagnant solution film adhering to the electrode. It is well known that a diffusion layer arises between the bulk sample solution and the membrane surface, which is analogous and related to the formation of a hydrodynamic boundary layer [32-34]. Convection or flow in the sample solution finally determines the average thickness of this diffusion layer. Actually, ISE responses were reported to be very sensitive to the flow- or stirring-rate in the sample solution, especially at lower ion activities [13, 35]. For rotating ISEs, a similar dependence on the rotation frequency was observed [23, 36]. Although boundary-layer influences are widely manifest in real experiments on ISEs, they have so far not been considered in the respective digital simulations.
Here, we present a very simple but nevertheless powerful computer modeling of ion-selective membrane electrodes and liquid junctions. The numerical simulations make use of a finite-element procedure in the space and time domains, which can be easily processed with conventional software (e. g., MS Excel) without the need for complex mathematical evaluations. It leads to excellent results on the dynamic evolution of concentration profiles, potentials, and fluxes in the studied systems. The treatment accounts for influences of convection, flow, or stirring in the sample solution that act on the boundary diffusion layer and it is even capable of including the effects of an electrolyte flow through the whole system. For simplicity, space-charge influences are considered via interfacial potential differences. Also, all interfacial reactions are assumed to be in or near local equilibrium in order to minimize the number of parameters required for the calculations. The applicability of the computer simulation is demonstrated for different ion-selective membranes as well as for liquid junctions.
The unidimensional flux of ions Izi within a liquid phase (e.g., a membrane of thickness d) is usually described by Eq. (1), which considers contributions from diffusion (first term), migration (second term), and convection (third term) [1, 2, 32]:
where Ji, Di, Ci, and Zi represent the flux, the diffusion coefficient, the concentration, and the charge of the species, respectively, x is the distance (0 ≤ x ≤ d), t the time, ψ = Fϕ / RT a dimensionless potential function, ϕ the electrical potential, v the flow velocity in x-direction, F the Faraday constant, R the universal gas constant, and T the absolute temperature. Ion fluxes are also involved in the continuity equation, Eq. (2), and in the expression for the electrical current density j, Eq. (3):
where the sum includes all mobile ions in the system.
In the corresponding finite-element approach, the continuous liquid phase is replaced by N thin layers (elements) with individual but uniform properties. By this procedure, the differentials in the above equations can basically be approximated by the small differences found between neighboring elements. Hence, the relationships for ion transport in Eqs. (1) and (2) are transformed into:
where Ji, ν/ν+1 (t) is the flux from the ν-th to the (v+1)-th element, Ci, ν and Ci, ν+1 are the concentrations in the respective elements, δ = d / N is the thickness of each element, and Δt is a sufficiently small time interval. From Eqs. (4) and (5) we immediately get a result for the time-dependent concentration changes in the system:
where τ and ω are dimensionless expressions for the time and the flow velocity, respectively, and Do is an arbitrary diffusion coefficient, e. g., of a reference species. The potential differences entering in Eq. (6) are determined from Eqs. (3) and (4):
where the contributions on the right side represent the diffusion potential induced by concentration changes (first term), the streaming potential caused by a flow (second term), and the ohmic potential drop generated by a current (third term). Due to the electroneutrality condition, the streaming potential can usually be neglected even for cases with ω ≠ 0. An exception is found for a porous membrane with fixed ionic sites that is exposed to a flow [33, 37]. Here, it holds that Σ Zi Ci, ν ≠ 0 (for all ν) since the sum expressions in Eq. (9) consider only the mobile ions of the system. Nevertheless, a solvent flow may also have a minor, indirect influence in the case of site-free phases because it induces some changes in the concentration profiles of the solutes.
Equations (4), (6), and (9) offer the basis for a computer modeling of ion transport in liquid membranes and other solution phases. They permit it to simulate the evolution of concentration profiles, potentials, and ion fluxes in the space and time domains. The computations based on the finite-element method can be processed with minimal effort but very high quality. To this end, first, a given initial state has to be chosen for all elements of the system. Also, the complete set of experimental parameters should be specified, which includes the definition of the appropriate conditions for the boundary elements (see below). Then, the changes of concentrations, potentials, and fluxes with time (i. e., with increasing number of discrete steps Δτ) are evaluated using the calculation procedure for subsequent small time intervals. These time steps should be chosen sufficiently small in order to avoid instabilities or oscillations of the computed results. Finally, a steady state is established for the majority of the studied examples that are representative for different realistic systems. The whole numerical calculation is far from being involved and can easily be performed with conventional MS Excel software (Microsoft Corporation). Since the computation rate is on the order of 103 cycles per second with corresponding time steps of Δt ≥ 1 ms, the present virtual experiments are often carried out in a shorter time period than real experiments.
In the following, a series of computer simulations by the finite-element method are presented and discussed. The results include aqueous diffusion layers (liquid junctions), ion-selective electrode membranes in potentiometric experiments, and liquid ion-exchanger membranes under the influence of an electrical current. All calculations were performed for 25 °C.
In order to test the reported computer modeling and to assess its capabilities, the finite-element procedure was first applied to aqueous diffusion layers. Well-known examples are the liquid junctions formed between the electrolyte of a reference electrode and sample solutions [1, 2]. Liquid-junction potentials at steady state (zero current and zero flow) can be calculated exactly from the theories of Planck [38, 39] and Schlögl , or approximated very closely by Henderson’s formula [41, 42] (for a review, see [1, 2, 43]). Although liquid junctions are often applied with a controlled outward flow of electrolyte from the reference electrode, no such flow effects have so far been taken into consideration by theoretical treatments.
In the present finite-element approach, the aqueous diffusion layer was modeled by a total of only 21 elements, of which the elements with ν = 0 and ν = 20 were used for the respective boundaries. Accordingly, the distance between the central planes of the two boundary elements was 20 times the thickness δ of a single element. Since this distance is identical to the total thickness d of the diffusion layer, which was assumed to be 0.1 cm, one obtains δ = 0.005 cm. The boundary concentrations were set equal to the bulk concentrations of the adjoining solutions (reference electrolyte for ν = 0, and sample solution for ν = 20). As the initial state, the ion concentrations of the elements ν = 1 to ν = 19 were chosen to form linear profiles within the diffusion layer, which corresponds to Henderson’s assumption. Accordingly, the initially determined liquid-junction potentials,
were comparable to the values quoted from the Henderson equation. After applying the numerical procedure with time intervals of Δt = 0.02 s, the computer simulation finally approximated the steady-state behavior of the system. The final liquid-junction potentials were therefore in agreement with the respective Planck potentials (or the ones obtained from Schlögl’s extension). All calculations were performed for zero-current conditions. Values for the diffusion coefficients were quoted from equivalent ionic conductivities, λi = |Zi| Di F2/ RT [33, 44] (nearly identical data are given in ). Since these diffusion coefficients are slightly different from the ones used earlier [1, 43], the resulting potentials also show minor deviations.
Results on the numerical simulation of liquid junctions are shown in Fig. 1 and summarized in Table 1. For these computations the water flow in the system was assumed to be zero. Evidently, the liquid-junction potentials derived from the finite-element method are in excellent agreement with the values predicted either from Planck’s or from Henderson’s theory (Table 1). This is also true for the illustrated concentration profiles at steady state (Fig. 1), which exactly conform to Planck’s theory. It can be concluded that the present numerical procedure is excellently suited for the intended applications. As pointed out earlier , a very high quality of the computer simulation can be achieved even when the real system is modeled by only a relatively small number of finite elements.
The reported analysis by numerical procedures opens new possibilities for studying the exact influence of water flow in liquid junctions. Fig. 2 shows some results obtained for the system KCl HCl at flow rates of up to 0.004 cm s-1, which corresponds to a time of 25 s for crossing the liquid junction. The comparison with the zero-flow case given in Fig. 1c indicates that an increased flow leads to a pronounced polarization of the concentration gradients on one side of the liquid junction. For chloride ions, even a concentration maximum is observed (Fig. 2). It results from the action of oppositely directed driving forces, namely the water flow streaming from the left to the right side, and the electrical field acting in the other direction. In spite of these effects, the flow-induced change of the liquid-junction potential remains relatively low. Actually, the deviation from the zero-flow value is only <7 % (see Figs. 1 and and2).2). This agrees well with the practical experience that reference electrodes with electrolyte flow do not exhibit any abnormalities in their potential.
In an additional test, the present finite-element procedure was applied to liquid membranes with different ionic sites that are exposed to an electrical current. Early studies of this type were restricted to ion-exchanger systems with one kind of mobile ionic site or ligand [18, 19]. Later, valinomycin membranes with tetraphenylborate sites were also treated  when ionophore-based membranes with ionic additives became a focus of interest [45-47]. It has been recognized in the meantime that liquid membranes with the commonly used polymeric matrices incorporate nearly immobile or even fixed anionic sites as well [6, 30, 48]. Moreover, the addition of lipophilic salts that generate both cationic and anionic sites in the membrane phase is often recommended for certain ion-selective electrodes . Accordingly, a recent theoretical treatment focused on the steady-state electrochemical properties of membranes containing any number and any kind of ionic sites .
In the finite-element approach, the ion-transport system was again modeled by 21 elements, including the two boundaries. The thickness of each element was δ = 0.001 cm in the case of a 0.02 cm thick membrane. The species considered in the membrane were permeating cations M+, mobile cationic sites Q+, mobile anionic sites R-, and fixed anionic sites S-. The initial concentration of a given species was taken to be the same in all elements, which corresponds to a zero-current steady state. It was assumed that the mobile sites are confined to the membrane phase because of their high lipophilicity. Hence, the average total concentration of cations or anions remains constant at a level chosen here as 0.01 M. In addition, the condition of zero flux through the two membrane surfaces applies for these species. In classical theories of diffusion, such an impermeable surface can be treated by the method of reflection at a boundary  which assumes an identical, oppositely directed system existing on the other side of the plane (e. g., from ν = 0 to ν = -20). This implies that slightly modified boundary conditions hold instead of Eqs. (5) and (6) (with ω = 0):
where the index i denotes mobile ionic sites. For simplicity, the same diffusion coefficient Di = Do = 10-8 cm2 s-1 was inserted for all mobile species in the membrane. Analogous expressions were used for the second membrane surface.
The finite-element calculation with time intervals of Δt = 2 s then led to results for the time-dependent variation of all concentration profiles in the system. The membrane-internal potential difference at a given current density was obtained from Eq. (9). The two boundary potentials were described using the assumption of an equilibrium distribution of permeating cations between membrane surfaces and external solutions [1, 2, 15, 50, 52, 53]. Finally, the total membrane potential was given according to Eq. (13):
where a’m and a”m are the cation activities of the solutions contacting the boundaries at ν = 0 and ν = 20, respectively. For simplicity, these activities were set equal to the respective bulk values, and it was assumed that a’m = a”m. It should be noted that the above approach does not account for any effects of ion-pair formation. On the other hand, neutral-ionophore-based membranes are treated correctly as long as they contain an excess of free ionophores [1, 3, 52]. Extensions aiming at more complex systems are within the realm of possibility.
Results obtained for a membrane with permeating cations and anionic sites are given in Fig. 3. Both mobile and fixed sites were considered as membrane components at relative amounts of 75 % and 25 %, respectively. The dynamic evolution of concentration gradients under the influence of a current is illustrated. The current density was chosen as 4.488 μA cm-2, which results in a total membrane potential of -500 mV (and an internal potential difference of -438.21 mV) for the final steady state. Evidently, this current flow induces an enrichment of mobile sites on one side of the membrane, and a depletion on the other side [18-20, 50]. The numerical results obtained for the final state can directly be compared with exact data since a theoretical description of the steady state is available . This comparison leads to the highly surprising and remarkable fact that the numerical predictions perfectly correlate with the exact values. The relative deviations are found to be ≤0.001 % for the concentration profiles, the internal potential difference, as well as the current density. Apparently, the computer modeling of membranes is a highly successful and efficient method, even when the number of finite elements considered is relatively small.
Fig. 4 shows results for a membrane that also contains cationic mobile sites, in addition to permeating cations and anionic membrane components. The two cationic species were assumed to have the same initial concentration. After switching on a constant current, we observe characteristic changes in the concentration profiles, as illustrated in Fig. 4. Cationic and anionic mobile sites are evidently shifted in opposite directions, but the overall polarization of concentrations is much less pronounced than in the former case (Fig. 3). Actually, a considerably lower current density of 0.6073 μA cm-2 is here sufficient to yield a final membrane potential of -500 mV when using the parameters specified above. The numerical results obtained for the final steady state can again be compared with theoretical data reported earlier  and are found to be in excellent agreement. Further comments on the electrical properties of membranes with different ionic sites are found elsewhere .
The primary aim of the present computer modeling was to investigate the complex response behavior of ion-selective membrane electrodes (ISEs) in the diffusion-controlled domain. Deviations from Nernstian responses as well as variations in the apparent selectivities are often found when ISEs are exposed to samples of different ions, especially at lower activities [3, 10, 54]. These effects generally result from ion-exchange processes at the membrane surface and the concomitant counter-fluxes of primary and interfering ions. Since diffusion layers are formed that usually increase with time, the observed phenomena are highly involved. Accordingly, earlier theoretical approaches were restricted to steady-state [3, 9, 15, 53, 55] or pseudo-steady-state cases [1, 3]. In contrast, a computer simulation based on finite-element procedures opens new possibilities for analyzing the dynamic response characteristics of ISEs.
The model membrane studied here was again 0.02 cm thick and consisted of finite elements having a thickness of δ = 0.001 cm. The species considered in the membrane were two cations, I+ and J+, and any mobile or fixed anionic sites, X-. For simplicity, all mobile species in the membrane were assumed to have the same diffusion coefficient Do = 10-8 cm2 s-1. Accordingly, the total site (or cation) concentration, taken as Xtot = 0.01 M, became invariant within the membrane. The reason is that, in this case, no diffusion potentials were generated at zero current. This implies that the flux descriptions in Eqs. (4) and (6) could be applied in a considerably reduced form since ψν = ψν+1 holds (for ω = 0):
A finite-element approach was also used to account for the aqueous diffusion layer that is formed between the outer membrane surface and the bulk of the sample solution. This so-called Nernst layer is directly related to the hydrodynamically controlled stagnant layer [32, 33]. It was represented by elements being numbered from ν = 0 (membrane surface) to ν = -5 (bulk sample). Unless otherwise specified, the diffusion layer was assumed to have an average thickness of d’N = 0.005 cm, which corresponds to a thickness of δ’ = 0.001 cm for a single element. If the sample contains a sufficiently concentrated background electrolyte, the ion fluxes in the Nernst layer are basically given by pure Fickian diffusion. Hence, it was possible to make use of expressions similar to Eqs. (14) and (15):
where the quantities marked with (‘) refer to the aqueous diffusion layer, p is the permeability ratio and q the thickness ratio between aqueous and membrane elements. Analogous expressions were applied for the second cationic species in the sample boundary, and a common diffusion coefficient of D’o = 10-5 cm2 s-1 was used. Time steps of Δτ = 10-5 were chosen for the numerical evaluation of Eqs. (15) and (17), which corresponds to time intervals of Δt = 1 ms.
Among different possibilities of treating the ion transfer at the phase boundary, a description based on the continuity law in Eq. (20) and the distribution relationship in Eq. (21) turned out to be most appropriate and successful:
where a’i and a’j are the respective ion activities in the sample boundary, and Kij is the selectivity coefficient of the ISE for J+ relative to I+, taken here as 10-6. For simplicity, no diffusion layer was considered for the aqueous zone near the inner phase boundary of the ISE. Instead, the inner membrane surface was assumed to be in equilibrium with the bulk of the inner reference solution (activities a”i and a”j). Finally, the membrane potential of the ISE was obtained from Eqs. (13) and (21) as:
All virtual experiments performed in this work are based on the same set-up and follow a well-defined procedure. The ISE membrane is assumed to be freshly prepared (or perfectly preconditioned) with the ion J+ as the only cationic species. In common practice, this situation is encountered when ionic components are incorporated as membrane additives. The membrane is interposed between a sample solution and an inner reference solution that both contain the initial ion J+ at identical activities of 0.1 M. The experiment is started by adding the primary ion I+ on the sample side at a bulk activity of 0.1 M. After a given measuring period, which involves a fixed number of time steps, the activity of I+ is reduced by a factor of 10. The procedure is repeated until a minimal activity of 10-7 M is reached. Then, the measurements are continued, but this time with increasing I+ activities. The whole measuring cycle is carried out at a constant background activity of J+ (0.1 M) and under controlled stirring or flow conditions in the sample solution.
Fig. 5 displays potential responses of ISEs resulting from such virtual experiments. The numerical simulations were made for measuring times of 2.5, 10, 40 and 160 s per sample, respectively. All response curves show a distinct region with a super-Nernstian slope at lower activities. This phenomenon is the consequence of a siphoning-off of ions I+ from the sample boundary and their replacement by ions J+ leaving the membrane [1, 56]. Fig. 5a also documents a pronounced hysteresis, which indicates that the measuring periods were too short in this case. A comparison of the results in Figs. 5b-d leads to the conclusion that the analytical range with a linear response can be extended to lower activities by simply increasing the measuring times. This fact was clearly confirmed by real experiments on ISEs [10, 23, 36]
Fig. 6 illustrates the influence of another experimental parameter, namely the stirring or flow rate of the sample solution. Results are shown on the response characteristics of ISEs with aqueous diffusion layers of 100, 50, 25, and 12.5 μm thickness, respectively. These and all the following computer simulations were performed for fixed measuring periods of 10 s per sample. Surprisingly, the response curves in Figures 6 have very similar shapes as the former ones in Fig. 5. A reduction of the stagnant film on the electrode surface obviously has the same favorable effect as before an increase of the measuring time. This improvement can be achieved by an efficient sample stirring or flow, since the average thickness of the Nernst layer is usually inversely proportional to the square root of the stirring or flow rate [33, 34][. Another possibility is the application of rotating disk ISEs that were introduced recently [23, 36]. The results in Fig. 6 are in excellent agreement with response curves observed in real measurements on ISEs .
Fig. 7 shows response curves for ISE membranes that differ in the total concentration of anionic sites. In this case, the model system with a 50-μm thick aqueous diffusion layer was studied in numerical experiments. It becomes evident from Fig. 7 that the ion-exchange capacity of the ISE membrane is highly important with respect to super-Nernstian responses. In fact, a reduction of the site concentration (Figs. 7b and 7c) immediately reduces or even eliminates the activity range where the reference ISE (Fig. 7a) exhibits the typical non-ideal behavior. On the other hand, the addition of ionic membrane components is often a prerequisite for improving the ion selectivity and other specifications [45-47, 57]. Accordingly, a compromise has to be found when it is the aim to optimize the membrane composition for general applications.
Fig. 8 shows results for a series of ISEs for which the diffusion coefficient of ions in the membrane phase is different. The underlying viscosity changes of the membrane can easily be realized by just varying the polymer / plasticizer ratio [10, 58, 59]. Fig. 8 makes evident that the mobility of ions in the membrane is of similar importance as the ion-exchange capacity studied before in Fig. 7. Large deviations from a linear response are apparently found for ISE membranes of low viscosity, and vice versa. These findings well agree with observations made on ion-selective microelectrodes. Freshly prepared microelectrodes with pure liquid membrane phases often exhibit nonlinearities and even hystereses of the response curves. However, the performance of microelectrodes can be greatly improved by the addition of a certain percentage of PVC to the membrane liquid .
In conclusion, the present computer simulations allow it to perform virtual experiments on ISEs and to analyze the response behavior ab initio. The documented results are fully underscored by the experience gained from realistic measurements, whenever these are available.
This work was financially supported by the Swiss National Science Foundation, an internal research grant of the ETH Zurich, and the National Institutes of Health (grant R01-EB02189).