Subscriber access provided by Brown University Library

Article

Spatiotemporal pH dynamics in concentration polarization near ion-selective membranes Mathias Baekbo Andersen, David M. Rogers, Junyu Mai, Benjamin Schudel, Anson V Hatch, Susan B. Rempe, and Ali Mani Langmuir, Just Accepted Manuscript • DOI: 10.1021/la5014297 • Publication Date (Web): 03 Jun 2014 Downloaded from http://pubs.acs.org on June 23, 2014

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Langmuir is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

Spatiotemporal pH dynamics in concentration polarization near ion-selective membranes Mathias B. Andersen,† David M. Rogers,‡,§ Junyu Mai,¶ Benjamin Schudel,¶ Anson V. Hatch,¶ Susan B. Rempe,∗,‡ and Ali Mani∗,† Mechanical Engineering Department, Stanford University, Stanford, CA, 94305, Center for Biological and Material Sciences, Sandia National Laboratories, Albuquerque, NM, 87185, and Biological Sciences and Engineering Department, Sandia National Laboratories, Livermore CA 94551 E-mail: [email protected]; [email protected]

Abstract We present a detailed analysis of the transient pH dynamics for a weak, buffered electrolyte subject to voltage-driven transport through an ion-selective membrane. We show that pH fronts emanate from the concentration polarization zone next to the membrane, and that these propagating fronts change the pH in the system several units from its equilibrium value. The analysis is based on a 1D model using the unsteady Poisson–Nernst–Planck equations with non-equilibrium chemistry and without assumptions of electroneutrality or asymptotically thin electric double layers. Nonequilibrium chemical effects, especially for water splitting, are shown to be important ∗ To

whom correspondence should be addressed Engineering Department, Stanford University, Stanford, CA, 94305 ‡ Center for Biological and Material Sciences, Sandia National Laboratories, Albuquerque, NM, 87185 ¶ Biological Sciences and Engineering Department, Sandia National Laboratories, Livermore CA 94551 § Current address: Department of Chemistry, University of South Florida, Tampa, FL, 33620 † Mechanical

1

ACS Paragon Plus Environment

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

for the dynamical and spatiotemporal evolution of the pH fronts. Nonetheless, the model also shows that at steady state the assumption of chemical equilibrium can still lead to good approximations of the global pH distribution. Moreover, our model shows that transport of the hydronium ion in the extended space charge region is governed by a balance between electromigration and water self-ionization. Based on this observation, we present a simple model showing that the net flux of the hydronium ion is proportional to the length of the extended space charge region and the water self-ionization rate. To demonstrate these effects in practice, we have adopted the experiment of Mai et al. [ACS Nano 6, 10206 (2012)] as a model problem, and by including the full chemistry and transport we show that the present model can capture the experimentally observed pH fronts. Our model can, among other things, be used to predict and engineer pH dynamics, which can be essential to the performance of membrane-based systems for biochemical separation and analysis.

Introduction Micro- and nanofluidic lab-on-a-chip systems have attracted considerable attention in the past decade due to their favorable properties for small-scale biochemical analysis including sample separation, preconcentration, and detection. 1 Around a decade ago, it was realized that the interface region at an ion-selective element (such as a nanopore or a membrane) is effective for biochemical separation and analysis when subjected to an electrical bias. 2–5 However, even though the initial studies were promising, the fundamental understanding of basic physicochemical transport phenomena at such interfaces was lacking and is still not fully developed. This understanding is not only sought in the context of biochemical analysis, but is also essential in other applications relying on ion-selective interfaces such as electrodialysis and production of chemicals, 6 microfluidic pumps, 7 as well as soil remediation. 8 Depending on the geometric confinement and fluidity of the electrolyte several phenomena can play a role in the transport characteristics spanning over concentration polarization, 9

2

ACS Paragon Plus Environment

Page 2 of 36

Page 3 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

deionization shocks, 10,11 electrohydrodynamic instabilities, 12 and enhanced chemical reactions. 13 In this work, we focus on elucidating the role played by chemical interactions in the electrolyte within and outside of a membrane, and the consequences these have for biochemical transport. Our work is distinct from previous ones in a number of aspects, including non-equilibrium chemical effects, transport of pH, buffering weak acids/bases, consideration of non-electroneutrality effects, and numerical resolution of all electrochemical boundary layers. Hydrodynamic effects are assumed to be mostly negligible due to the presence of immobilizing polymers. ∆V |t>0

Figure 1: Schematic of the microfluidic membrane system. At t = 0 a potential difference ∆V is applied between the inlets (x1 and x6 ) and leads to concentration polarization around the membrane. Moreover, as predicted in this study, pH fronts emanate from the membrane and lead to variation in the pH of up to several units. Positive direction is from left to right. Figure Figure 1 shows a schematic of the prototypical microchannel-membrane system studied here. Such a system was recently used to elucidate aspects of chemical effects in experiments by Mai et al., 14 in which a ratiometric dye was used to make concentrationindependent measures of pH. These experiments allow for a much more rigorous comparison to theory due to the observations in both space and time, whereas previous studies have been restricted to only one of these dimensions and mostly the temporal one. 9,15–18 Here, we develop a model that captures the observed pH variations in both space and time thereby achieving heretofore unprecedented prediction capabilities. Some of the key features of the model are non-equilibrium chemical effects, transport of pH and directly capturing charge dynamics without assuming electroneutrality. The ability to predict and engineer local pH at nanoporous interfaces is crucial for the use 3

ACS Paragon Plus Environment

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

of such systems in the detection and analysis of biological tracers, which are often extremely sensitive to the local pH environment. Furthermore, the mesoscopic continuum model developed in the present work is envisioned to inform more detailed nanoscale models, such as molecular dynamics simulations of ionic transport into nanopores, 19 with local spatiotemporal pH, electric field, and ionic strength. Such a combined approach of multiscale modeling is a necessary predictive tool for the systematic development of microfluidic membrane systems for detection and analysis of biomolecules. Initial feasibility studies of microfluidic membrane systems for biochemical analysis and preconcentration have already been demonstrated. For instance, Hatch et al. 20 used a membrane for preconcentration of proteins prior to electrophoretic separation in a sieving gel and showed detection of proteins down to concentrations of 50 fM. It was also found that concentration polarization has a detrimental effect on the reproducibility, which may be understood by improved and more accurate models like the one pursued in the current work. Additionally, portable devices for the detection of biological toxins and pathogens was demonstrated by Meagher et al., 21 who showed the viability of integrating microfluidic membrane systems, similar to the one considered in this work, into a portable, self-contained device. General reviews of the various transport mechanisms at ion-selective membranes subject to high electrical forcing was recently given by Nikonenko et al. 13,22 In these reviews, overlimiting current, i.e., transport of charge at rates beyond diffusion limitation, is generally attributed to four distinct mechanisms: (i) water splitting, 23–25 (ii) current-induced convection, such as electro-osmotic instability, 26,27 (iii) co-ion leakage, and (iv) currentinduced membrane discharge. 28 Furthermore, the reviews point out that there is a strong need for additional understanding of the fundamental phenomena including how water splitting is proceeding at the membrane-solution interface and how weak acids and bases influence transport in the system. In the present work, we focus on such chemical interactions and assume that convection phenomena, such as the electro-osmotic instability, 27,29 can be ignored due to both the confinement effects by the small dimensions of the microchannel and

4

ACS Paragon Plus Environment

Page 4 of 36

Page 5 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

nanopores, passivation of microchannel surface charge by covalent polyacrylamide grafting and the addition of viscosity-enhancing agents to the solution. 14 There has been less attention to how ionic transport in porous media couples with chemical interactions. The majority of previous studies have assumed either steady state, 30–32 electroneutrality, equilibrium chemistry, fixed pH values, binary electrolytes, or combinations thereof. 16,23,33–38 All of these assumptions are relaxed in the present work. A key point to emphasize in the context of chemical interactions is the presence of both homogeneous and heterogeneous reactions, i.e, reactions both in the bulk solution as well as at the surface groups on the pore wall. 39 Water splitting was one of the first mechanisms proposed to explain the ionic transport properties during strong electrical forcing at ion-selective membranes. Simons 40 considered non-equilibrium of the water splitting reaction and argued that the large fluxes of hydronium and hydroxyl ions could be attributed to catalyzing effects of the protonation and deprotonation reactions of water with the surface groups on the membrane. This catalysis hypothesis, which essentially corresponds to larger numerical values of the water splitting rate constants, has been supported in several other studies. 9,17,18,25,35,36,41 The catalysis effect would also provide an explanation for the difference in the water splitting between membranes constituting different surface groups per the observation that water splitting in many cases is faster at anion-selective membranes relative to cation-selective membranes. 9,17,18,41 A systematic review of the water splitting mechanism was given early on by Zabolotskii et al., 25 who stressed the importance of including water splitting in modeling of membrane systems in the overlimiting-current regime. In many studies, 9,17,18,41 temporal pH changes up to several units have been measured on both sides of anion- and cation-selective membranes, but without accounting for the spatial distribution, as done here. Measurements of the partial currents have also shown that the increased flux of H+ and OH− in general is too small to account for the observed overlimiting current. 9,41 Recently, bipolar membranes, essentially two membranes of opposite surface charge in close contact, have been used as effective water

5

ACS Paragon Plus Environment

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

splitters as a means for regulating the pH in microfluidic channels. 42,43 Previous work on systems related to the one studied here has been done in the group of Tallarek, where several studies of preconcentration in microfluidic channels have been performed. 44–46 However, where our focus is on the coupling of chemical effects with driftdiffusion, their focus has primarily been on the coupling of advection with drift-diffusion. The detailed model developed in the present work can be used to map out a reducedorder (or lumped parameter/short-cut) model, 43 which could be more feasible for engineering design and optimization studies. In particular, the present model has the capability to study membranes of lower acidity (higher pK value). Membrane acidity is another feature that can be leveraged in the engineering of biochemical separation and analysis. 47–50 Furthermore, the present model can be used to explore possible strategies to avoid unwanted concentration polarization and pH variations. These strategies include pulsed electrical forcing 51 and biasing the porous membrane with a gating voltage. 47 The paper is organized in the following format. In section Figure 1, we introduce the system geometry and definitions. The governing equations with boundary and initial conditions are developed. A short description of the numerical model is given. Then, in section Eq. (14b), we present the results and discussions. We start by giving a thorough motivation of the parameters used in the model. Then, we compare the theoretical predictions of pH dynamics to the experimental observations. Potential causes for discrepancies are discussed. After the benchmark, the model is used to investigate the detailed dynamics of the system including an analysis of the much-debated assumption of chemical equilibrium. Section Eq. (19) gives a summary and conclusions.

6

ACS Paragon Plus Environment

Page 6 of 36

Page 7 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

Table 1: Overview of the ionic species and their electrochemical properties. Throughout the text, we will interchangeably use index notation (i = 1, 2, ..., 8) and name-notation (cation=cat, anion=an, buffer=buf, water=w, membrane=mem), as well as denote concentration either by brackets [Czcat ] or by ci . Note, that as the unit of concentration is mM, the unit of the equilibrium constant and the reaction rate constant is mM as well. However, the + ] pH value is based on the concentration of hydronium with unit of M, i.e., pH = − log10 ( [H 1 M ). Mobile ions Fixed ions symbol description Salt ions Buffer ions Water ions Membrane ions i species index 1 2 3 4 5 6 7 8 zi valence zcat zan zbuf zbuf − 1 −1 +1 zmem zmem − 1 ci concentration [Czcat ] [Azan ] [BHzbuf ] [Bzbuf −1 ] [OH− ] [H+ ] [MHzmem ] [Mzmem −1 ] Ki equilibrium constant Kbuf Kw Kmem ki reaction rate constant kf,buf , kb,buf kf,w , kb,w kf,mem , kb,mem total concentration cbuf cmem

Mathematical model System geometry and definitions As shown in figure Figure 1, we employ a 1D model of the system, with the axial coordinate being x. Thus, the model assumes uniformity in the transverse directions and neglects fluid flow. An ion-selective membrane, which spans the entire cross section, is located between x3 and x4 . The inlets of the microchannel, between which a voltage ∆V is applied at time t = 0, are located at positions x1 and x6 . There is a large separation of length scales in the system, ranging from the small Debye screening length of O(1 nm) up to the microchannel length of ∼ O(1 cm), a ratio of 107 . That large range of length scales makes modeling of the system challenging since we do not resort to any asymptotic approximations of the thin boundary layers such as the electric double layer (EDL) and the extended space charge layer (ESC). 52 However, during the time of interest, which is of O(10 s), large parts of the microchannel remain unchanged. We can model these sections (from x1 to x2 and from x5 to x6 ) as Ohmic resistors of fixed concentration. This confinement of the variable regions close to the membrane is sketched in figure Figure 1, where the pH profile outside the membrane is uniform at equilibrium, but nonuniform between x2 and x5 during voltage bias. The model is verified by comparison to the experimental data from Mai et al., 14 in which phosphate buffered saline was used as electrolyte. Thus, as indicated in table Table 1, the 7

ACS Paragon Plus Environment

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 36

model is generalized to take into account multiple ionic species. Those ionic species are: a fully dissociated salt cation Czcat of valence zcat , a fully dissociated salt anion Azan of valence zan , a buffer consisting of a weak acid BHzbuf of valence zbuf and its conjugate base Bzbuf −1 of valence zbuf − 1, the hydroxyl ion OH− and hydronium ion H+ , and inside the membrane a fixed weak acid surface group MHzmem of valence zmem and its conjugate base surface group Mzmem −1 of valence zmem −1. We assume that the concentration of the dye that tracks the pH is sufficiently small to be neglected. The model accounts for three chemical reactions, namely the buffer reaction, the water splitting reaction, and the protonation/deprotonation reaction inside the membrane. These chemical reactions are not assumed to be in equilibrium, but are modeled directly using their respective forward kf,i and backward kb,i reaction rate constants, as well as their equilibrium dissociation constants Ki . As indicated in figure Figure 1, for the purpose of comparing to the experimental data, the ionic species are identified to be sodium 2− Na+ , chloride Cl− , and phosphate ions H2 PO− 4 and HPO4 . Also, since the membrane used

in Mai et al. 14 was negatively charged, the surface groups in the experiment are assumed to titrate between negative and neutral forms.

Governing equations We employ the standard continuum formulation of ionic mass conservation in the limit of an ideal dilute electrolyte, 53

∂ci ∂ji + = ri , ∂t ∂x

(1)

where ci , ji and ri are the concentration (with unit of mM = mol/m3 ), flux and source term, respectively, of species i. As we have neglected fluid flow, the flux is due solely to diffusion and electromigrative drift, and, continuing with the assumption of an ideal dilute electrolyte, given by the Nernst–Planck equation,

8

ACS Paragon Plus Environment

Page 9 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

ji = −αDi

!

∂ci zi ∂φ ci + , ∂x VT ∂x

(2)

in which Di , zi , VT = RT /F , R, F , and φ are the diffusion constant, ionic valence, thermal voltage, gas constant, Faraday constant and electrostatic potential, respectively. The membrane factor α is introduced to account for effects such as the porosity, tortuosity and constrictivity incurred by the membrane. α has a value of unity in the microchannel (x1 < x < x3 and x4 < x < x6 ) and a value below unity in the membrane (x3 < x < x4 ),     1,

in the microchannel,

α=   αmem ,

(3)

in the membrane.

Flux matching is always ensured at the microchannel-membrane interfaces (at x3 and x4 ). We emphasize our assumption of a dilute electrolyte, which corresponds to all activity coefficients being unity. Some non-ideal effects at high concentrations are therefore not captured, but we sacrifice this accuracy to allow a more direct probing of effects of interest, i.e., the importance of non-equilibrium chemistry and pH dynamics. The source term ri in (Eq. (1)) is due to chemical reactions between the ionic species. Our model accounts for three reactions of the protonation/deprotonation type: (i) a buffer reaction, (ii) the water splitting reaction, and (iii) the charge regulation reaction in the membrane, kf,buf

buf buf −1 BHz(aq) ⇋ Bz(aq) + H+ (aq) ,

kb,buf kf,w

+ H2 O(aq) ⇋ OH− (aq) + H(aq) , kb,w

mem MHz(s)

kf,mem

mem −1 + H+ ⇋ Mz(s) (aq) .

kb,mem

9

ACS Paragon Plus Environment

(4a) (4b) (4c)

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 36

Here, kf,i and kb,i are the forward and backward reaction rate constants, respectively, with their unit based on mM and not M. The phase of the species has been explicitly indicated in (Eq. (4)) to emphasize that, while the buffer and water splitting reactions occur exclusively between dissolved species, the charge regulation reaction occurs between the solid phase mem mem −1 and Mz(s) ) and the dissolved hydronium ion (H+ membrane species (MHz(s) (aq) ). Follow-

ing classical chemical kinetics, 54 the source term ri is the difference between the forward and backward reaction rates, found as the product of the rate constant by the concentrations (with r1 = r2 = 0 since the salt ions are completely dissociated), r3 = kb,buf [Bzbuf −1 ][H+ ] − kf,buf [BHzbuf ] 

(5a) 

= kb,buf [Bzbuf −1 ][H+ ] − Kbuf [BHzbuf ] , r5 = kf,w [H2 O] − kb,w [OH− ][H+ ]

(5b) (5c)





= kb,w Kw − [OH− ][H+ ] ,

(5d)

r7 = kb,mem [Mzmem −1 ][H+ ] − kf,mem [MHzmem ]

(5e) 



= kb,mem [Mzmem −1 ][H+ ] − Kmem [MHzmem ] ,

(5f)

where Ki = kf,i /kb,i is the equilibrium constant with unit based on mM and not M. In the special case of the water splitting reaction, Kw = kf,w [H2 O]/kb,w with the concentration of water [H2 O] being assumed constant. We note again, that an overview of the ionic species and their indices is given in table Table 1. The remaining source terms are expressed in terms of the ones in (Eq. (5)),

r4 = −r3 ,

(6a)

r6 = −r3 + r5 − r7 ,

(6b)

r8 = −r7 .

(6c)

The membrane ions (i = 7, 8) constitute a special case for which the modeling can be 10

ACS Paragon Plus Environment

Page 11 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

immediately reduced to a single variable in the following way. We add the conservation equation (Eq. (1)) whereby, by virtue of (Eq. (6c)), the source terms cancel. Additionally, since the membrane ions are fixed, their flux is identical to zero. The resulting equation states that the sum of the concentrations of the membrane ions is constant in time,

c7 + c8 = cmem ,

where cmem is the total concentration of fixed surface groups in the membrane. We make use of this condition to eliminate species i = 8 from the system of equations, whereby the source term for species i = 7 becomes

r7 = kb,mem [(cmem − c7 ) c6 − Kmem c7 ] .

(7)

Electrostatic interactions are governed by the Poisson equation,

8 X ∂ ∂φ zi ci , −ε = ρe = F ∂x ∂x i=1

!

(8)

where ε = εr ε0 is the dielectric permittivity in which εr and ε0 are the relative and vacuum permittivity, respectively, and ρe is the free charge density. Note that all ionic species are included in the sum of the charge density in (Eq. (8)), including the membrane ions that are responsible for the fixed charge in the membrane. It is this smeared-out, or volume-averaged, charge from the membrane ions that leads to membrane selectivity, thus concentration polarization, and other interesting effects to be analyzed in this work. The length scale over which electrostatic relaxation occurs is the Debye length

λD =

s

εRT , 2F 2 I

11

ACS Paragon Plus Environment

(9)

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 36

in which

I=

6 1X z 2 ci 2 i=1 i

(10)

is the ionic strength, essentially the effective concentration of the electrolyte. Note that our definition of the ionic strength only accounts for the mobile ions, i.e., it does not include the fixed ionic surface groups in the membrane. The Debye length is the characteristic width of the EDLs at the microchannel-membrane interfaces at x3 and x4 , and is around 1 nm for an electrolyte with an ionic strength of 100 mM.

Boundary and initial conditions Equilibrium ionic composition We assume that the system is fully equilibrated at t = 0 right before the voltage is applied. Thus, we find the initial ionic composition by solving the governing equations while neglecting unsteady and gradient terms, i.e., essentially solving ri = 0 subject to the electroneutrality condition, 6 X

zi ci = 0.

i=1

12

ACS Paragon Plus Environment

(11)

Page 13 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

Furthermore, ri = 0 corresponds to chemical equilibrium whereby it is possible to express through (Eq. (5)) all ionic concentrations in terms of the hydronium concentration, Kbuf , Kbuf + [H+ ] [H+ ] , [BHzbuf ] = cbuf Kbuf + [H+ ] Kw [OH− ] = + , [H ] Kmem , [Mzmem −1 ] = cmem Kmem + [H+ ] [H+ ] [MHzmem ] = cmem , Kmem + [H+ ] [Bzbuf −1 ] = cbuf

(12a) (12b) (12c) (12d) (12e)

where cbuf is the total concentration of buffer ions. Plugging these relations into the electroneutrality condition (Eq. (11)) yields a nonlinear algebraic equation for the hydronium concentration that can be solved numerically. The result is the equilibrium concentration of [H+ ] in free solution, which in combination with (Eq. (12)) gives the equilibrium concentration of all ionic species. These equilibrium concentrations are applied as fixed Dirichlet boundary conditions at x2 and x5 (see figure Figure 1) under the assumption that these positions are far enough from both the driving electrodes at x1 and x6 and the membrane interfaces at x3 and x4 so as not to experience any significant change in the concentrations during the time of interest. In practice, the widths L23 = x3 − x2 and L45 = x5 − x4 of the fully resolved parts of the microchannel are taken to be around 15 times the membrane width Lmem = x4 − x3 , while the widths L12 = x2 − x1 and L56 = x6 − x5 of the remaining parts of the microchannel is of order cm, i.e., tens of thousands times larger than Lmem . However, as we shall see, most of the voltage drop occurs in the computational region shortly after the voltage is applied.

13

ACS Paragon Plus Environment

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 36

Matching of electric field The electrostatic boundary condition is derived by assuming uniform conductivity from the driving electrodes at x1 and x6 to the positions x2 and x5 . This agrees with the assumption of negligible concentration variations. Furthermore, the overpotentials driving the electrode reactions are assumed to be at most of O(1 V), which is negligible compared to the typical applied voltages of O(100 V) and above. Thus, in accordance with Ohm’s law,

i = σE,

(13)

where i, σ and E = −∂x φ are the current density, conductivity and electric field. The electric field is uniform from x1 to x2 and from x5 to x6 . Additionally, since at x2 and x5 the current is Ohmic, continuity of the current density results in mixed Robin boundary conditions for the potential at x2 and x5 ,

(14a)



(14b)

φ|x2 − φ|x1 ∂φ , = L12 ∂x x2 φ|x6 − φ|x5 ∂φ = , L56 ∂x x5

in which the electrode potentials φ|x1 = ∆V and φ|x6 = 0. Initial condition with membrane The initial condition is found numerically as the steady-state solution to the governing equations (Eq. (1)) and (Eq. (8)) with the above-mentioned boundary conditions. In practice, to find the steady-state solution, an iterative procedure is employed in which the membrane concentration cmem is increased from a small fraction of its final value, and where the initial guess for each iteration is the solution from the previous step. In this manner, the initial condition for the unsteady calculations is ensured to be physically equilibrated and numerically consistent. The initial condition found in this manner produces a Donnan potential 14

ACS Paragon Plus Environment

Page 15 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

consistent with theory. 55

Numerical scheme We develop our own numerical scheme to solve the governing equations (Eq. (1)) and (Eq. (8)). The numerical integration is challenging in particular because the system is characterized by large separations in time and length scales, ranging from the small charge-relaxation time (~0.1 ns) to the diffusion time (~10 s), and from the small Debye length (~1 nm) to the microchannel length (~1 cm). Our numerical scheme uses a nonuniform grid that fully resolves the smallest features, such as the EDL and ESC boundary layers, and models the large microchannels through effective boundary conditions. Additionally, the scheme is based on an implicit time integration and the time step is adjusted appropriately throughout the simulation such that the transient dynamics is properly resolved. The scheme conserves mass at the level of discretization and is second-order accurate in both space and time. It is based on a so-called delta-form iteration to converge on the solution for each time step. The governing equations are solved in a dimensionless form, which can be seen in the supplemental material. Moreover, details of the discretization in delta-form are also found in the supplemental material. The scheme is implemented and solved in Matlab.

Results and discussion We model the experimental setup used in Mai et al. 14 and note that even though the data from this experiment is novel, it is also limited and there is uncertainty about the experimental value of some of the parameters in the model. A list of the parameters and their values is given in table Table 2. Thus, we do not aim for a conclusive determination of the model parameters by rigorously fitting to the experimental data. Rather, we will emphasize the fact that the model predicts qualitatively the same novel spatiotemporal pH dynamics as seen in the experiment. Additionally, the model provides detailed insight into the dynamics and the

15

ACS Paragon Plus Environment

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 36

Table 2: List of parameters. The specific choice of some parameter values are discussed in the main text. Note, that due to tradition in the literature, pK-values is computed based on an equilibrium constant K in units of molar M and not in mM, i.e., pK = − log10 ( 1KM ). Not listed in the table are: vacuum permittivity ε0 = 8.854 × 10−12 F m−1 , gas constant R = 8.314 J mol−1 K−1 , and Faraday constant F = 96485 C mol−1 . Symbol

Value

Unit

Na diffusivity Cl− diffusivity H2 PO− 4 diffusivity HPO2− 4 diffusivity OH− diffusivity

1.3 × 10−9

m2 s−1 56

D6

D1 D2 D3

Description +

Ref. Symbol Description zcat zan zbuf

2.0 × 10−9 m2 s−1 0.96 × 10−9 m2 s−1

56

0.76 × 10−9 m2 s−1

56

zmem

5.3 × 10−9

m2 s−1

56

pKbuf

H+ diffusivity

9.3 × 10−9

m2 s−1

56

pKw

ccat

Na+ concentration

160

mM

14

pKmem

can

Cl− concentration

150

mM

14

kb,buf

cbuf

Total concentration of buffer ions

6

mM

14

kb,w

cmem

Total volume concentration of surface groups in membrane Applied voltage Width of membrane Width of resolved parts of microchannel Width of remaining parts of microchannel

300

mM

100

V

14

αmem

5 × 10−5

m

14

εr

7.5 × 10−4

m

2 × 10−2

m

D4 D5

∆V Lmem L23 , L45

L12 , L56

56

kb,mem

T

Cation valence Anion valence Buffer top valence Membrane top valence Buffer reaction equilibrium constant Water splitting equilibrium constant Membrane reaction equilibrium constant Backward rate constant of buffer reaction Backward rate constant of water splitting reaction Backward rate constant of membrane reaction

1 −1 −1

Membrane factor Relative permittivity Temperature

0.2

14

16

Value Unit

ACS Paragon Plus Environment

Ref.

0 7.2

56

14

56

−2

107

mM−1 s−1

109

mM−1 s−1

107

mM−1 s−1

78 293

K

Page 17 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

electrochemical structure of the system, which is not observable in experiments. Moreover, the model illustrates to which degree the assumption of chemical equilibrium holds. Some of the reasoning behind the values in table Table 2 chosen for the model is as follows. Regarding the salt and buffer concentration, Mai et al. 14 used a commercial buffer specified as 150 mM sodium chloride (NaCl) with 10 mM of phosphate (NaH2 PO4 , Na2 HPO4 and Na3 PO4 ). For the phosphate component, we ascribe 10 mM of sodium and adjust cbuf until the pH value at equilibrium in the model match that in the experiment, which results in cbuf = 6 mM. We then follow the stated values and use can = 150 mM and ccat = 160 mM, with the latter resulting as the sum of the contributions from the sodium chloride and phosphate components. We note that for these concentrations the Debye length is λD ≈ 0.8 nm outside the membrane. The total concentration of surface groups in the membrane cmem is adjusted by hand until optimal agreement is found between model and experiment leading to cmem = 300 mM. This is a reasonable value considering that Mai et al. 14 used an acrylamide monomer and bisacrylamide cross-linker mixture with a concentration around 6×103 mM to polymerize the membrane. That concentration sets the upper limit for the concentration of surface groups in the membrane. Typically, only a small fraction of the amine groups on the acrylamide monomer hydrolyze into carboxylic groups, and the ratio of ∼ 1/20 here is consistent with typical numbers in the literature. We also note that for this membrane concentration the Debye length is λD ≈ 0.5 nm inside the membrane, which our model fully resolves. The nature of the chemical groups inside the membrane is unknown. Complementary experiments would be needed for an independent determination of pKmem . We use pKmem = −2 since this gives robust agreement between model and data. This low value of pKmem minimizes current-induced membrane discharge effects and renders the membrane fully ionized. A more rigorously determined value of pKmem would be based on additional experimental data currently not available. Rate constants can be difficult to find in the literature. For the surface protonation/deprotonation

17

ACS Paragon Plus Environment

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 36

reaction in the membrane, we note that Hoogerheide et al. 57 measured a backward reaction rate for silanol groups on silica of 7 × 108 mM−1 s−1 , which indicates that surface reactions can proceed fast. Based on this, we settle for kb,mem = 107 mM−1 s−1 , and note that the model shows little sensitivity to the actual value of kb,mem . For the buffer reaction, we use kb,buf = 107 M−1 s−1 based on estimates from literature. 58 Again, the model show little sensitivity to the exact value of kb,buf . For the ubiquitous water splitting, the rate constant in free = 1.4 × 108 mM−1 s−1 . 58 However, as explained in section Footfree solution is known, kb,w

note §, there is a large body of literature showing that the rate constant of water splitting is enhanced by order of magnitudes in concentration polarization. 9,13,15,18,25 We therefore use kb,w = 109 mM−1 s−1 as the nominal value and present the impact of variation of this constant over orders of magnitude. We emphasize that the simple uniform increase in kb,w employed here mimics the more detailed models quite naturally since the actual value of kb,w only becomes important in the concentration depletion zone. The membrane factor αmem accounts for the impact of porosity, tortuosity and constrictivity on the ionic transport in the membrane. It is difficult to determine a priori and it was not measured in the experiments. Thus, in accordance with the literature, 28,44,59 we use αmem = 0.2.

Comparison to experimental data: spatiotemporal pH fronts and current Figure Figure 2 shows model prediction of pH fronts around an ion-selective membrane in qualitative agreement with the experimental data. Both experiment and model show pH variations up to several units in magnitude emanating from the membrane over the course of tens of seconds. These pH modulations occur despite the relatively large phosphate buffer concentration. The immediate conclusion for biomolecule separation and analysis is clear, namely that pH around ion-selective membranes cannot be assumed to be uniform. There are several potential causes for the quantitative discrepancy between the model and 18

ACS Paragon Plus Environment

Page 19 of 36

pH

(a ) 10

(b ) pH

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

experiment

9 8 7 6 10 9 8 7 6

0s 10 s 20 s 30 s 40 s

model 12 10 8 6 4 2 −400−200 0 200 400

−400 −200

0 200 x [mu]

400

Figure 2: pH profiles across the membrane at different times: (a) Experimental data. (b) Model prediction. The left edge of the membrane is at x = 0 and its width is 50 µm. The inset shows the full range in pH predicted by the model. the data in figure Figure 2. For one, there is the experimental uncertainty, which manifests itself in various ways: (i) The accuracy of the ratiometric fluorescence microscopy is highest at a pH around 7 and saturates below 6 and above 9. Titration experiments in bulk solution show the ratiometric approach to be relatively insensitive to pH below 6 and the dye was depleted below its detection limit in the depletion zone. pH values outside the range between 6 and 9 from experimental measurements should be viewed as being uncertain, and, as seen in figure Figure 2, the pH moves significantly outside these limits. (ii) The system may not have been completely equilibrated before the start of the experiment, since at t = 0 (blue curve in figure Figure 2) the pH to the left of the membrane does not equal that to the right of the membrane, but is around 0.5 units higher. From this effect alone we could expect deviations between model and experiment of up to 0.5 pH units. Before moving on to a more quantitative analysis, we address some uncertainties in the model. As seen in figure Figure 2, the agreement in pH between model and experiment is less successful to the left of the membrane, where the propagation of the pH front is too fast in the model prediction. To the right of the membrane the model predicts very low pH values.

19

ACS Paragon Plus Environment

Langmuir

These discrepancies could be due to the omission of the two phosphate buffer reactions + H3 PO4 ⇋ H2 PO− 4 +H ,

pK = 2.1,

(15a)

3− + HPO2− 4 ⇋ PO4 + H ,

pK = 12.7.

(15b)

Given that the pK value of the phosphate reaction we do account for is 7.2, the neglect of (Eq. (15)) is less accurate for low and high pH values and especially critical for pH below 4 and above 10. Reactions (Eq. (15)) act to buffer, or slow down and minimize, the speed and magnitude of the pH fronts and, as is clear in figure Figure 2, this effect would improve the agreement between model and data. The reason that reactions (Eq. (15)) are not included in the model in its current implementation is primarily that the additional dependent variables would make the problem less numerically tractable. However, even with these limitations, the model reasonably captures the main trends of the data in figure Figure 2. (a ) x pH=7 .0 [µm]

400 300 200 100 (b) 0 i [kA/ m2 ]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 36

4 3 2 1 0 0

experiment k b,w model 6 4 2 0 0 0.05 0.1 0.15 0.2

10

20

experimental asymptote 30

40

50

t [s]

Figure 3: Comparison between the experiment (red circles) and the model employing different values of the water splitting backward reaction rate constant, kb,w /(mM−1 s−1 ) = 109 (blue curve), 1010 (cyan curve), 1011 (orange curve), and 1012 (purple curve). Kw is kept constant. (a) The position where the pH equals 7 versus time. (b) The current density versus time. The experimental asymptote is indicated (red dashed line). We proceed with a more quantitative comparison between model and experiment. As mentioned above, for both the experiment and the model, we have the largest confidence in intermediate pH values around 7. It therefore makes sense to use this pH value for a 20

ACS Paragon Plus Environment

Page 21 of 36

more quantitative comparison and we introduce xpH=7.0 as the position to the right of the membrane (x > x4 ) where the pH equals 7. Figure Figure 3(a) plots xpH=7.0 as a function of time and shows for kb,w = 109 mM−1 s−1 good agreement between the model and data, both of which exhibit a linear growth in time with similar slope, or speed of propagation of the pH front. The model prediction of xpH=7.0 is displaced slightly below the experimental one, but given the uncertainties discussed above this is not too surprising. The other model curves in figure Figure 3(a) pertain to increasing values of kb,w and will be addressed below. Measurements of the current provide a complementary observable for comparison between model and data. Figure Figure 3(b) plots the current as a function of time. The time resolution of the experimental current measurements is somewhat coarse with a measurement only at every 30 s. On the other hand, as seen in figure Figure 3(b) and its inset, the model predicts variations in the current on time scales down to 0.1 s. Given these few experimental data points in the time interval captured by the model, we also plot in figure Figure 3(b) the experimental asymptote in the limit of long times (the experimental current settles at this level after approximately 120 s). From Figure 3(b) we note that model and data agree on

10

0

−2

−0.01

0 50 x [µm]

(d)

(f )

0

t=0 s t = 0 .1 s t=1 s t= 10 s t= 40 s

−50

[O H − ] [mM ]

10

2

(e)

100

[H + ] [mM ]

10

(c )

[B 2− ] [mM ]

[C + ] [mM ]

(a ) 2 10 0 10 −2 10 −4 10 −6 10 −8 10 (b) 2 10 0 10 −2 10 −4 10 −6 10 −8 10

[B H − ] [mM ]

the magnitude of the current and its decreasing trend in time toward an asymptotic value.

[A − ] [mM ]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

−50

0 50 x [µm]

100

−50

0 50 x [µm]

100

Figure 4: Concentration versus position of all the ionic species at different times. The inset in panel (a) illustrates that the model fully resolves the rapid, but smooth and continuous, transition across the EDL boundary layers at the microchannel-membrane interfaces. 21

ACS Paragon Plus Environment

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 36

Concentration polarization It is well known that concentration polarization occurs around ion-selective membranes when subjected to an electric field. 53 In this phenomenon, concentration decreases on one side of the membrane and increases on the other. The concentration difference across the membrane grows both in time and with the strength of the applied electric field. For a suddenly applied strong electrical forcing, either at or above the limiting current, the Sand’s time,

τsand =

π zF c0 D 4 i 

2

,

(16)

is the characteristic time for the concentration to reach approximately zero on the depletion side of the membrane. We assume z = 1, D = 1.5 × 10−9 m2 s−1 and c0 = 150 mM, while the inset in figure Figure 3(b) gives the current density, i = 5 × 103 A m−2 . With these parameters, τsand ≈ 0.01 s, which is an order of magnitude less than the 0.1 s observed by the kink in the curve in the inset of figure Figure 3(b). However, this discrepancy is expected since the expression for Sand’s time in (Eq. (16)) is based on an ideally selective membrane. As the ratio of the electrolyte ionic strength to the membrane concentration is quite high, here I/cmem ≈ 150/300 = 1/2, the membrane is far from ideal. Thus, longer time is expected for salt depletion to develop and this is in accordance with the observed discrepancy. The concentration polarization phenomenon can be seen in figures Figure 4(a) and (b), which show the concentration of cations and anions at different times. These salt ion concentrations reach values down to 10−4 − 10−3 mM in the depletion zone to the left of the membrane (x < 0 µm), and they increase slightly to the right of the membrane (x > 50 µm). The transitions in the concentration fields across the EDL boundary layers (at x = 0 µm and 50 µm) appears discontinuous on the micrometer scale, but the inset in figure Figure 4(a) illustrates that it is in fact smooth and fully resolved by the model. While the salt ions behave as expected throughout the concentration polarization process, there is a different and more rich behavior of the chemically reactive ions. Figure Figure 4(c)

22

ACS Paragon Plus Environment

Page 23 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

through (f) indicate that the dynamics of the buffer and water ions are strongly influenced by chemistry. At the very initial times, t . 0.1 s (red curves in figure Figure 4), before enhanced water splitting sets in, the buffer ions follow classical concentration-polarization depletion/enrichment behavior. While the same is true for the hydroxyl ion, interestingly, the hydronium depletes to the right of the membrane and enriches to the left of the membrane, opposite to the behavior of all the other ions in the system. Then, at intermediate times 0.1 s . t . 10 s (green and black curves in figure Figure 4), the effect of water splitting is seen by propagating fronts of hydronium: an enrichment front going to the right and a depletion front going to the left. The change in the hydronium concentration drives the buffer reaction in (Eq. (4a)) toward protonation of the buffer to the right of x = 0, and vice versa to the left of x = 0. Finally, at longer times, t & 10 s (magenta curves in figure Figure 4), the propagating hydronium front moving right makes it through the membrane and continues in the microchannel to the right. Of course, these hydronium fronts are the same as the pH fronts in figure Figure 3, but instead of only observing the pH, we have through the model gained additional insight into the dynamics of all the ions in the system. As discussed above, water splitting leads to increased flux of water ions out of the depletion zone. As seen in figure Figure 4(f), the concentration of water ions reaches significant values up to at least 1 mM. Considering their uniquely high mobility, the water ions start to carry a considerable fraction of the current. A measure of the fraction of current carried by an ion is given by the transference number,

ti =

zi F ji . i

(17)

The transference numbers of the hydronium and hydroxyl ions are plotted in figure Figure 5. In accordance with the above observations, the front of hydronium propagating to the right causes an increasing fraction of the current in this part of the system to be carried by that ion. The same is true regarding the hydroxyl ion in the part of the system to the

23

ACS Paragon Plus Environment

(a )

10 10

Page 24 of 36

−1

−2

(b) 10−1 tH +

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

tOH −

Langmuir

t=0 s t = 0 .1 s t=1 s t= 10 s t= 40 s

−2

10

−200 −100

0 100 x [µm]

200

300

Figure 5: Transference number versus position at different times: (a) Hydroxyl transference number tOH− . (b) Hydronium transference number tH+ . left of the membrane. Moreover, the transference number shows that only around 1 % of the current is carried by the water ions.

Extended space charge We turn our focus to the concentration depletion zone next to the membrane on its left side. This is also where the extended space charge, or ESC, layer forms. Figure Figure 6 shows the pH, electric field, ionic strength and free charge density in the ESC layer. In fact, as seen in figure Figure 6(d), the ESC layer can be defined as the appearance of a considerable amount of free charge density, which extends much farther away from the membrane than the EDL. In this case, at t = 40 s, the ESC has an extension of approximately 15 µm, which is more than 104 times the EDL width of 0.8 nm. The large separation of the free charge from the membrane surface makes the ESC layer hydrodynamically unstable, which can lead to strong local advection and overlimiting current. 26,27,29 However, we neglect such hydrodynamic effects here. There are additional interesting phenomena in the ESC layer which are worth noting, especially in the context of biochemical separation and analysis. The state and transport properties of biomolecules primarily depend on pH, local electric field, and ionic strength. Thus, the observations made here can inform more detailed modeling, such as molecular dynamics simulations, 19 about the value of these fields in the critical microscopic region 24

ACS Paragon Plus Environment

Page 25 of 36

pH

(a ) 10 8 6 4 2

E [V /m ]

(b)

7

10

5

10

t=0 s t = 0 .1 s t=1 s t= 10 s t= 40 s

3

10

(c ) I [mM]

2

10

0

10

−2

10

(d) ρe /F [mM]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

−1

10

−2

10

−3

10 −40

−30

−20 −10 x [µm]

0

Figure 6: Dynamics of various fields in the ESC region: (a) pH. (b) Electric field E. (c) Ionic strength I. (d) Normalized free charge density ρe /F . where biomolecules enter the membrane. As seen in figure Figure 6(a), the pH is almost uniform and close to its original value in the ESC. However, the pH is much larger to the left of the ESC and much lower inside the membrane. Figure Figure 6(b) displays the huge electric field in the ESC, in this case up to at least 107 V m−1 . This field strength has to be contrasted to the average field between the electrodes of 2500 V m−1 , whereby it is clear that almost the entire applied voltage drops across the microscopic ESC region. These field strengths make the linearized Poisson–Boltzmann equation inapplicable. It is also interesting to compare these field strengths to those across biological membranes. The voltage drop across a cell (∼ 60 mV) or a muscle cell (∼ 80 mV) occurs over a 3-4 nm thick membrane. 60 That gives a field on the order of 107 V m−1 , which is comparable to those observed here. Finally, figure Figure 6(c) shows that the ionic strength reaches very low values in the ESC, in this case down to 10−2 mM from its equilibrium value of approximately 150 mM, i.e., a drop by a factor of more than 104 .

25

ACS Paragon Plus Environment

Langmuir

Non-equilibrium chemistry and sensitivity to the water splitting reaction rate constant Chemical equilibrium is often assumed when dealing with transport of ions. 16,33–38 By this assumption, it is often possible to reduce the number of dependent variables and eliminate numerically troublesome fast chemical time scales. In the present case, that assumption leads to the relations in (Eq. (12)). The assumption of chemical equilibrium is based on the chemical time scale τchem ∼ 1/(kb c0 ) being much faster than the transport time scales of diffusion τdiff ∼ L2 /D and electromigration τem ∼ LVT /(zDE). However, when considering hydronium in the ESC region, kb = 109 mM−1 s−1 (the rate enhanced by factor of 10), c0 = O(10−4 mM), E = O(107 V m−1 ) and L = O(1 µm) whereby τchem = O(10−5 s) and τem = O(10−7 s). Thus, by this simple estimate, transport by electromigration is faster than chemical reaction in the ESC region, and chemical equilibrium is doubtful. − 1

1

[OH − ][H + ] Kw

(a )

0 −1 −2

− 1

(b) [B 2− ][H + ] [B H − ]K buf

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 26 of 36

t=0 s t = 0 .1 s t=1 s t= 10 s t= 40 s

0

−1 −40

−30

−20 −10 x [µm]

0

Figure 7: Deviation from chemical equilibrium versus position at different times: (a) Water splitting reaction. (b) Buffer reaction. We can use our model to further investigate the assumption of chemical equilibrium. We quantify the deviation from equilibrium in the water splitting and buffer reactions by the

26

ACS Paragon Plus Environment

Page 27 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

departure from zero in a rearranged version of the relations in (Eq. (12)), [OH− ][H+ ] − 1 = 0, Kw [B2− ][H+ ] − 1 = 0. [HB− ]Kbuf

(18a) (18b)

That is, the left-hand-sides in (Eq. (18)) quantifies the deviation from chemical equilibrium. As shown in figure Figure 7, the left-hand-sides of (Eq. (18)) are not zero in the ESC region, and chemical equilibrium does not hold there. Moreover, the deviations approach minus unity, which can be explained in the following way. In and near the ESC layer, there is a local depletion of H+ and OH− , and a stronger local depletion of B2− relative to HB− . Thus, [OH− ][H+ ] and [B2− ][H+ ]/[HB− ] drop toward zero magnitude. When unity is subtracted, the result is close to −1 (see figure Figure 7). The deviation for the buffer ions extends further than that of the water ions due to the slower chemical time scale for the buffer reaction. The positive values of the deviation for the water ions (figure Figure 7(a)) is likely related to the coupling of H+ to the buffer reaction. To investigate the consequences of assuming chemical equilibrium of the water splitting reaction, formally corresponding to the limit kb,w → ∞, we consider in our model three additional values of kb,w , which are successively increased by an order of magnitude, while keeping Kw constant. It was shown in figure Figure 3(a) how these increases in kb,w lead to overprediction of the speed of the propagating pH front. This clearly illustrates that the assumption of chemical equilibrium is erroneous when studying dynamical phenomena in concentration polarization, like a propagating pH front as done here. Additionally, this is also a strong illustration of the power of the present study in accounting for both space and time, thereby gaining insight not obtainable when considering either of these two dimensions separately, as done in previous studies. Nevertheless, our model shows that the assumption of chemical equilibrium indeed can lead to good approximations of the global pH distribution in special cases such as, e.g., steady

27

ACS Paragon Plus Environment

Langmuir

(a ) pH x =2 5 µm

8 6 4 2 k b,w

(b) pH x =6 0 µm

8 6 4 2

inside membra ne

outside membra ne

k b,w −1

10

10

0

1

10 t [s]

Figure 8: The pH as a function of time at two different positions for increasing values of the water splitting backward reaction rate constant, kb,w /(mM−1 s−1 ) = 109 (blue curve), 1010 (cyan curve), 1011 (orange curve), and 1012 (purple curve): (a) Inside the membrane at x = 25 µm. (b) Outside the membrane at x = 60 µm. state. 28,30,31,38? Figure Figure 8 shows the evolution of the pH at two different positions in the system, inside the membrane (x = 25 µm) and outside the membrane (x = 60 µm), for the four different values of kb,w . While the transition from high to low pH proceeds faster with increasing kb,w , the steady-state pH profile appears relatively insensitive to kb,w . (a ) |∂ x j Hem+ | [mM/ s]

1

(b)

10 0 10 −1 10 −2 10

t=0 s t = 0 .1 s t=1 s t= 10 s t= 40 s

1

|r 5 | [mM/ s]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 28 of 36

10 0 10 −1 10 −2 10 −40

−30

−20 −10 x [µm]

0

Figure 9: Absolute magnitude of two dominating terms in the equation of mass conservation for hydronium: (a) The divergence of the electromigration flux. (b) The source term due to water splitting. Our model also allows us to gain insight into the mechanisms driving the water splitting in the ESC region. Figure Figure 9 shows the absolute magnitude of two of the terms in the governing equation for mass conservation of the hydronium ion, the divergence of the + em = ∂ (−αD electromigration flux ∂x jH + x H+ [H ]∂x φ/VT ) and the water-splitting source term

r5 = kb,w Kw (1 − [OH− ][H+ ]/Kw ). The remaining terms in the governing equation are not 28

ACS Paragon Plus Environment

Page 29 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

shown. As can be seen in figure Figure 9, the model nevertheless shows that the primary balance in the ESC region is between those two terms, ∂x jH+ ∼ r5 , where we have assumed em . This is in contrast to the traditional non-reacting models that assume a balance jH+ ≈ jH +

between diffusion and electromigration in the ESC. At the same time, figure Figure 7(a) shows that the water splitting term can be approximated by unity, 1 − [OH− ][H+ ]/Kw ≈ 1, whereby the balance simplifies to ∂x jH+ ∼ kb,w Kw . We integrate this balance from the left to the right edge of the ESC, while employing, as shown in figure Figure 5(b), that jH+ is approximately zero at the left edge of the ESC. Thus, the hydronium flux injection at the right edge of the ESC becomes

jH+ |x=0 ∼ kb,w Kw LESC ,

(19)

where LESC is the width of the ESC region. For t = 10 s, LESC ≈ 15 µm, and by (Eq. (19)), jH+ |x=0 = 1.5 × 10−4 mM m s−1 . Then, with a current of i = 103 A m−2 , this flux corresponds by (Eq. (17)) to a transference number of tH+ ≈ 0.015, in good agreement with the full simulation results in figure Figure 5. The approximation in (Eq. (19)) is very useful since the scaling of the length of the ESC is already known from the existing literature, 52 whereby (Eq. (19)) provides insightful a priori estimates of the hydronium ion flux.

Conclusion We have developed a 1D model capable of accurately integrating the unsteady Poisson– Nernst–Planck equations. The model is verified by comparison against experimental data by Mai et al. 14 of pH measurements around a nanoporous membrane in a microfluidic channel subjected to an external electric field. A major novelty of the present work is the consideration of both space and time for the study of pH dynamics in concentration polarization. This aspect is highlighted by the model being able to predict the propagating pH fronts observed 29

ACS Paragon Plus Environment

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 30 of 36

in the experiment. Both model and experiment clearly emphasize that pH cannot be assumed to be uniform, nor constant, during concentration polarization. This is an important point since this type of microfluidic-membrane system has gained considerable popularity as a means for biochemical processing and analysis in applications such as detection of biological warfare agents. 14,20,21,50 Since biomolecules are often very sensitive to pH, accounting for variations in pH with a model of the present type is important for interpreting observations. Furthermore, the model can be used to design and optimize microfluidic-membrane systems. We re-emphasize that this model and the dramatic shifts observed for pH warrants attention to possibility of similar phenomenon impacting broad-ranging devices where microchannels intersect with nanofeatures. Finally, there is potential in the present type of study to elucidate many of the still unanswered questions regarding transport in concentration polarization. The spatiotemporal analysis uniquely allows testing of aspects not possible through purely spatial or temporal approaches. For instance, the present study shows that the much debated assumption of chemical equilibrium is erroneous for dynamical phenomena in concentration polarization. Nonetheless, it is furthermore shown that chemical equilibrium can lead to a good approximation of the global pH distribution under special circumstances such as under steady state conditions reached after long times. Also, the present model gives insight that allows us to develop a simple scaling arguments. For instance, as shown in (Eq. (19)), we find that the hydronium ion flux scales with with the length of the extended space charge region, the rate constant, and the equilibrium constant.

Acknowledgement This project received support from the Defense Threat Reduction Agency-Joint Science and Technology Office for Chemical and Biological Defense (IAA number DTRA10027IA-3167). This work was also supported in part by Sandia’s LDRD program and the National Science 30

ACS Paragon Plus Environment

Page 31 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

Foundation under Grant No. PHYS-1066293 and the hospitality of the Aspen Center for Physics. Sandia National Laboratories is a multiprogram laboratory managed and operated by Sandia Corporation, a wholly owned subsidiary of Lockheed Martin Corporation, for the U.S. Department of Energy’s National Nuclear Security Administration under Contract DE-AC04-94AL8500.

Supporting Information Available A document detailing the non-dimensionalization and numerical discretization is given in the supporting information.

This material is available free of charge via the Internet at

http://pubs.acs.org/.

References [1] Jong, J. d.; Lammertink, R. G. H.; Wessling, M. Lab Chip 2006, 6, 1125–1139, 00197. [2] Khandurina, J.; Jacobson, S. C.; Waters, L. C.; Foote, R. S.; Ramsey, J. M. Anal. Chem. 1999, 71, 1815–1819, 00268. [3] Song, S.; Singh, A. K.; Kirby, B. J. Anal. Chem. 2004, 76, 4589–4592, 00127. [4] Wang, Y.-C.; Stevens, A. L.; Han, J. Anal. Chem. 2005, 77, 4293–4299, 00394. [5] Foote, R. S.; Khandurina, J.; Jacobson, S. C.; Ramsey, J. M. Anal. Chem. 2005, 77, 57–63, 00188. [6] Strathmann, H. Desalination 2010, 264, 268–288, 00094. [7] Suss, M. E.; Mani, A.; Zangle, T. A.; Santiago, J. G. Sens. Actuators, A 2011, 165, 310–315, 00017. [8] Probstein, R. F.; Hicks, R. E. Science 1993, 260, 498–503, 00454. [9] Krol, J. J.; Wessling, M.; Strathmann, H. J. Membr. Sci. 1999, 162, 145–154, 00146. 31

ACS Paragon Plus Environment

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

[10] Mani, A.; Zangle, T. A.; Santiago, J. G. Langmuir 2009, 25, 3898–3908, 00082. [11] Mani, A.; Bazant, M. Z. Phys. Rev. E 2011, 84, 061504, 00023. [12] Chang, H.-C.; Yossifon, G.; Demekhin, E. A. Annu. Rev. Fluid Mech. 2012, 44, 401–426, 00029. [13] Nikonenko, V. V.; Pismenskaya, N. D.; Belova, E. I.; Sistat, P.; Huguet, P.; Pourcelly, G.; Larchet, C. Adv. Colloid Interface Sci. 2010, 160, 101–123, 00065. [14] Mai, J.; Miller, H.; Hatch, A. V. ACS Nano 2012, 6, 10206–10215, 00001. [15] Simons, R. Electrochim. Acta 1984, 29, 151–158, 00209. [16] Ramírez, P.; Mafé, S.; Alcaraz, A.; Cervera, J. J. Phys. Chem. B 2003, 107, 13178– 13187, 00045. [17] Nikonenko, V.; Pis’menskaya, N.; Volodina, E. Russ. J. Electrochem. 2005, 41, 1205– 1210, 00023. [18] Zabolotskii, V.; Bugakov, V.; Sharafan, M.; Chermit, R. Russ. J. Electrochem. 2012, 48, 650–659, 00006. [19] Leung, K.; Rempe, S. B.; Lorenz, C. D. Phys. Rev. Lett. 2006, 96, 095504, 00043. [20] Hatch, A. V.; Herr, A. E.; Throckmorton, D. J.; Brennan, J. S.; Singh, A. K. Anal. Chem. 2006, 78, 4976–4984, 00141. [21] Meagher, R. J.; Hatch, A. V.; Renzi, R. F.; Singh, A. K. Lab Chip 2008, 8, 2046–2053, 00073. [22] Nikonenko, V. V.; Kovalenko, A. V.; Urtenov, M. K.; Pismenskaya, N. D.; Han, J.; Sistat, P.; Pourcelly, G. Desalination 2014, 00000. [23] Grossman, G. J. Phys. Chem. 1976, 80, 1616–1625, 00022. 32

ACS Paragon Plus Environment

Page 32 of 36

Page 33 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

[24] Simons, R. Desalination 1979, 28, 41–42, 00123. [25] Zabolotskii, V. I.; Shel’deshov, N. V.; Gnusin, N. P. Russ. Chem. Rev. 1988, 57, 801– 808, 00099. [26] Rubinstein, I.; Zaltzman, B. Phys. Rev. E. 2000, 62, 2238, 00227. [27] Druzgalski, C. L.; Andersen, M. B.; Mani, A. Phys. Fluids 2013, 25, 110804, 00006. [28] Andersen, M. B.; Soestbergen, M. v.; Mani, A.; Bruus, H.; Biesheuvel, P. M.; Bazant, M. Z. Phys. Rev. Lett. 2012, 109, 108301, 00036. [29] Davidson, S. M.; Andersen, M. B.; Mani, A. Phys. Rev. Lett. 2014, 112, 128302, 00000. [30] Lu, J.; Wang, Y.-X.; Zhu, J. Electrochim. Acta 2010, 55, 2673–2686, 00013. [31] Yeh, L.-H.; Zhang, M.; Qian, S. Anal. Chem. 2013, 85, 7527–7534, 00004. [32] Nielsen, C. P.; Bruus, H. Phys. Rev. E 2014, 89, 042405, 00000. [33] Nikonenko, V.; Lebedev, K.; Manzanares, J. A.; Pourcelly, G. Electrochimica Acta 2003, 48, 3639–3650, 00043. [34] Volgin, V.; Davydov, A. J. Membr. Sci. 2005, 259, 110–121, 00024. [35] Danielsson, C.-O.; Dahlkild, A.; Velin, A.; Behm, M. Electrochim. Acta 2009, 54, 2983– 2991, 00022. [36] Tanaka, Y. J. Membr. Sci. 2010, 350, 347–360, 00026. [37] Roelofs, S. H.; Soestbergen, M. v.; Odijk, M.; Eijkel, J. C. T.; Berg, A. v. d. Ionics 2014, 1–8, 00000. [38] Dykstra, J. E.; Biesheuvel, P. M.; Bruning, H.; Heijne, A. T. submitted to Phys. Rev. E 2014, 00000.

33

ACS Paragon Plus Environment

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

[39] Andersen, M.; Frey, J.; Pennathur, S.; Bruus, H. J. Colloid Interface Sci. 2011, 353, 301–310, 00028. [40] Simons, R. Electrochim. Acta 1985, 30, 275–282, 00132. [41] Choi, J.-H.; Moon, S.-H. J. Colloid Interface Sci. 2003, 265, 93–100, 00037. [42] Cheng, L.-J.; Chang, H.-C. Biomicrofluidics 2011, 5, 046502, 00014. [43] Stodollick, J.; Femmer, R.; Gloede, M.; Melin, T.; Wessling, M. J. Membr. Sci. 2014, 453, 275–281, 00001. [44] Dhopeshwarkar, R.; Crooks, R. M.; Hlushkou, D.; Tallarek, U. Anal. Chem. 2008, 80, 1039–1048, 00059. [45] Hlushkou, D.; Dhopeshwarkar, R.; Crooks, R. M.; Tallarek, U. Lab Chip 2008, 8, 1153– 1162, 00054. [46] Hlushkou, D.; Perdue, R. K.; Dhopeshwarkar, R.; Crooks, R. M.; Tallarek, U. Lab Chip 2009, 9, 1903, 00044. [47] Martin, C. R.; Nishizawa, M.; Jirage, K.; Kang, M.; Lee, S. B. Adv. Mater. 2001, 13, 1351–1362, 00132. [48] Chun, K.-Y.; Stroeve, P. Langmuir 2002, 18, 4653–4658, 00111. [49] Ali, M.; Ramírez, P.; Mafé, S.; Neumann, R.; Ensinger, W. ACS Nano 2009, 3, 603–608, 00095. [50] Sommer, G. J.; Mai, J.; Singh, A. K.; Hatch, A. V. Anal. Chem. 2011, 83, 3120–3125, 00007. [51] Malek, P.; Ortiz, J.; Richards, B.; Schäfer, A. J. Membr. Sci. 2013, 435, 99–109, 00002. [52] Yariv, E. Phys. Rev. E 2009, 80, 051201, 00012. 34

ACS Paragon Plus Environment

Page 34 of 36

Page 35 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

[53] Probstein, R. Physicochemical Hydrodynamics: An Introduction; John Wiley & Sons, 2005; 01988. [54] Bard, A.; Faulkner, L. Electrochemical methods: fundamentals and applications; Wiley New York, 1980; Vol. 2; 24329. [55] Galama, A. H.; Post, J. W.; Cohen Stuart, M. A.; Biesheuvel, P. M. J. Membr. Sci. 2013, 442, 131–139, 00008. [56] Haynes, W. M., ed. ed.; CRC Press/Taylor and Francis: Boca Raton, FL, 2010; 00004. [57] Hoogerheide, D. P.; Garaj, S.; Golovchenko, J. A. Phys. Rev. Lett. 2009, 102, 256804, 00053. [58] Eigen, M. Angew. Chem. Int. Ed. 1964, 3, 1–19, 01518. [59] Elattar, A.; Elmidaoui, A.; Pismenskaia, N.; Gavach, C.; Pourcelly, G. J. Membr. Sci. 1998, 143, 249–261, 00120. [60] Hille, B. Ion channels of excitable membranes; Sinauer Sunderland, MA, 2001; Vol. 507; 13774.

35

ACS Paragon Plus Environment

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Graphical TOC Entry ∆V |t>0

36

ACS Paragon Plus Environment

Page 36 of 36