ARTICLE pubs.acs.org/ac
Propagating Concentration Polarization and Ionic Current Rectification in a NanochannelNanofunnel Device Dzmitry Hlushkou,† John M. Perry,‡ Stephen C. Jacobson,*,‡ and Ulrich Tallarek*,† † ‡
Department of Chemistry, Philipps-Universit€at Marburg, Hans-Meerwein-Strasse, 35032 Marburg, Germany Department of Chemistry, Indiana University, Bloomington, Indiana 47405-7102, United States ABSTRACT: We study ionic current rectification observed in a nanofluidic device with a nanofunnel positioned between two straight nanochannels. Ion transport is simulated by resolving the coupled three-dimensional NernstPlanck, Poisson, and Navier Stokes equations. In the modeled system, the electric double layer extends into the channel, and consequently, the funnel tip exhibits charge-selective properties, which results in the formation of enriched and depleted concentration polarization (CP) zones within the nanofunnel in the high- and low-conductance states, respectively. This scenario is similar to the one observed for ion transport through a charged conical nanopore connecting two macroscopic reservoirs. However, the presence of the adjacent straight nanochannels allows the CP zones to propagate out of the funnel into the adjoining channels. The condition for propagation of the CP zones is determined by several parameters, including the electroosmotic flow velocity. We demonstrate that in the high-conductance regime the modeled system is characterized by increased ionic concentrations in the entire cathodic nanochannel, whereas in the low-conductance state the depleted CP zone does not propagate out of the funnel and remains localized. The required three-dimensional modeling scheme is implemented on a parallel computational platform, is general as well as numerically efficient, and will be useful in the study of more advanced nanofluidic device designs for tailoring ionic current rectification.
R
ecently, nanochannels and nanopores have received strong attention because of the intrinsic properties they exhibit. Many systems, such as synthetic membranes and microfabricated channels, already have diverse and outstanding applications, e.g., for analyte sensing, trapping, preconcentration, and separation processes. The simple geometries applied in such systems are commonly symmetric. At the same time, systems with broken geometrical symmetry, such as nanopipets, nanofunnels, nanonozzles, and conical nanopores, possess several unique and fascinating features. For instance, the transport properties of asymmetric channels and the corresponding phenomenological coefficients (e.g., hydraulic and electroosmotic permeabilities, ionic electric conductivity, diffusion coefficient, and mean residence time) are dependent on the direction of macroscopic mass and charge transport. This kind of anisotropy results in a number of intriguing phenomena, such as asymmetric diffusive flux,14 valveless hydraulic pumping,5,6 asymmetric electroosmotic flow (EOF),7,8 and ionic current rectification (ICR).915 Particularly, the latter two phenomena are respectively encountered when the EOF velocity and ionic current through a channel depend on the polarity of the applied bias. Though ICR is interesting from both fundamental and application points of view,16 there is still controversy about the origins of this effect. ICR occurs in systems with asymmetric interaction between the charged channel walls and ions of an electrolyte solution filling the channel. In 1997, Wei et al.9 proposed a model to explain the asymmetric behavior of the currentvoltage (IV) characteristics detected with opposite polarities of voltage r 2011 American Chemical Society
applied to a quartz nanopipet electrode system. According to that model, the asymmetry of the IV characteristics originates from differences in limiting transport rates of electrolyte ions flowing into and out of the pipet segment whose diameter is comparable with the electric double layer (EDL) thickness. The model includes parameters (e.g., access radii and selectivity parameter) that are very difficult to approach experimentally, and the system asymmetry is described with a heuristic factor relating the cationic and anionic components of the flux and not accounting for the electric properties of the system, such as the surface charge density and EDL thickness. Later, Woermann1719 developed another phenomenological approach to explain the data obtained by measurements of ionic current passing through single conical pores in poly(ethylene terephthalate) membranes. That approach attributes asymmetric branches of an IV curve to different transference numbers for the negative and positive bias. Woermann18 proposed as well a semiquantitative expression for the IV relationship which predicts two straight lines with different slopes for positive and negative applied voltages. The experimental IV data, however, demonstrate commonly more complex, nonlinear behavior. A probable explanation of that discrepancy is that Woermann’s model assumes a negligible distortion from the equilibrium state by an externally applied voltage, which is frequently not the case. To explain nonohmic IV characteristics of Received: September 21, 2011 Accepted: November 23, 2011 Published: November 23, 2011 267
dx.doi.org/10.1021/ac202501v | Anal. Chem. 2012, 84, 267–274
Analytical Chemistry systems rectifying ionic current, Siwy and co-workers11,12,20 proposed a model based on the concept of the electrostatic ratchet/trap mechanism for counterions of a background electrolyte. This model assumes that the electric potential profile within a channel can be represented as a superposition of the potential profiles due to a surface charge at the walls and an externally applied electric field. For an asymmetric pore, the former is characterized by a nonuniform shape and can be associated with an asymmetric potential well (electrostatic trap). Then imposition of an external electric field results in either a reduction (or even disappearance) of the potential barrier for counterions or an increase of the potential barrier, which depends on the applied bias polarity. An expression to describe a relationship between the ionic current and the local potential distribution along a conical nanopore was proposed.11,20 A notable difference with Siwy’s model is that the system is assumed to be far from equilibrium (in contrast to the models of Wei et al. and Woermann); the system is driven from equilibrium by the externally applied electric field as well as by thermal fluctuations. Along with the aforementioned theoretical approaches, a number of numerical models were recently applied to simulate ICR in conical nanopores. The first numerical model allowing a quantitative description of ionic currents through a conical nanopore was developed by Cervera et al.21 The model is based on a solution of the Poisson, NernstPlanck, and continuity equations in spherical coordinates relative to the cross-sectional averaged electric potential and species concentrations and assumes that the ionic fluxes have only a radial component. Output and input boundary conditions at the tip and base of the pore were determined from the Donnan equilibrium relations, and the only adjustable parameter of that model was the surface charge density at the pore walls. Later, several groups reported results of their numerical studies of the effects of the pore geometry, surface charge density, electrolyte ionic strength, and EOF and pressure-driven flow on ICR in conical pores.2229 All these models are based on a numerical analysis of the coupled partial differential equations governing mass and charge transport (the NernstPlanck, Poisson, and, sometimes, NavierStokes equations). The most relevant conclusion from the results of the numerical studies is that local concentrations of co-ions and counterions within a pore are different for positive and negative applied voltages, which results in different conductances in the system for opposite biases. Here, all aforementioned theoretical and numerical approaches either do not account for mass charge transport in the regions adjoining the conical pore at all (i.e., the boundary conditions are assumed to be imposed directly at the edges of the pore) or include these regions only as macroscopic reservoirs (i.e., their dimensions are much larger than the EDL thickness and the pore tip size) to account for eventual entrance effects. At the same time, the latter studies indicate that a difference in ionic species concentrations corresponding to positive and negative applied voltages is observed as well outside the conical pore.23,25,2729 In 2009, Jung et al.30 published their experimental study of ICR in a symmetric nanopore fabricated in series with two micropores of different diameters. The authors also proposed a qualitative explanation of the observed ICR, based on the asymmetry of concentration polarization (CP) developed at the nano/micro interfaces of the system they studied. CP is the electrochemical term for effects related to the formation of concentration gradients (depletion and enrichment zones) in the electrolyte solution adjacent to an interface between two domains with different ion
ARTICLE
transport characteristics upon the passage of a direct electric current.31,32 Because CP changes the local conductivity and electric field distribution, asymmetric CPs developing under opposite applied voltages result in different ionic currents through the same system. This explanation is, to the best of our knowledge, the first proposal to explain ICR in a geometrically asymmetric configuration by accounting for phenomena occurring within a symmetric pore with adjacent asymmetric regions.30 Recently, we reported an experimental study of ion transport in a nanofluidic system composed of two straight nanochannels separated by a nanofunnel that rectified ionic current similar to conical pores.33 The planar format of the device allowed simultaneous electric and optical characterization which showed that an anionic fluorescent probe was enriched and depleted not only within the funnel, but in the adjoining straight channels as well, even though these channels were several times longer than the funnel. To understand and in the future tailor ICR behavior of these more complex nanochannelnanofunnel systems, which so far have not been characterized through numerical simulations, we present a general, three-dimensional modeling approach and a set of first simulation results on ion transport in a nanofluidic system with a channelfunnelchannel configuration. These results show that depleted and enriched CP zones not only are developed within the nanofunnel, but also can propagate outward and occupy the adjacent nanochannels. This phenomenon is an essential component of device functionality and can be reasoned by the theory of propagating CP.34,35 This theory proposes the criterion determining whether the CP zones remain stationary in the vicinity of a charge-selective domain or propagate outward. That criterion is defined by several system parameters, in particular the EOF velocity. The difference in CP dynamics observed with opposite applied voltages can be explained by different EOF rates generated in the system under positive and negative bias of the same magnitude.
’ EXPERIMENTAL SECTION AND MODELING APPROACH As a model system we used a nanofluidic device that was fabricated by casting high-modulus poly(dimethylsiloxane) (h-PDMS) onto an SU-8 master written by e-beam lithography.33 The device was composed of a nanofunnel and two straight nanochannels that connect the tip and base of the funnel to two microchannels, which were spaced 85 μm apart. Each straight nanochannel was ∼40 μm long, 1.05 μm wide, and 120 nm high, whereas the funnel was 80 nm wide and 120 nm high at the tip and 1.05 μm wide and 120 nm high at the base. Three different designs were realized by using different lengths of the funnel, resulting in 5°, 10°, and 20° taper angles. A more detailed description of the fabrication technique can be found elsewhere.33 For our simulations we chose the geometry of the 5° funnel as it exhibited the highest ICR ratio during previous experiments. For the sake of illustration, the three-dimensional atomic force microscope image of the 5° funnel is shown in Figure 1a. In parts b and c of Figure 1, the layout of the integrated microchannel nanochannel device and a scanning electron microscope image of the 5° funnel are shown, respectively. Figure 1d also shows a 2D projection of the 5° funnel system with critical dimensions used in the simulations. The height (third dimension) of the system is constant and equal to 120 nm. The lengths of the straight channels shown in Figure 1d are comparable to the funnel length. 268
dx.doi.org/10.1021/ac202501v |Anal. Chem. 2012, 84, 267–274
Analytical Chemistry
ARTICLE
flow velocity field to the pressure field and the body force due to the electrostatic potential (Lorentz force). Equations 13 have been complemented with the corresponding boundary conditions formulated at the solidliquid interface (eq 4) and at the inlet/outlet of the system (eq 5). n 3 Ji ¼ 0,
n 3 ∇ϕ ¼ σ=ε,
∇|| ci, in ¼ ∇|| ci, out ¼ 0,
u¼0
ϕin ϕout ¼ A,
ð4Þ pin pout ¼ B ð5Þ
I ¼
Our modeling approach is based on a numerical resolution of the equations describing the coupled advectivediffusive and electromigration masscharge transport in the system: zi F Di ci ∇ϕ þ ci u RT
ε∇ 3 ð∇ϕÞ ¼ F
∑i zi ci
1 u∇ 3 u ¼ ð ∇p þ η∇2 u F F
ð1Þ ð2Þ
∑i zi ci ∇ϕÞ
S
F
∑i zi Ji 3 n dS
ð6Þ
The above mathematical model was realized as a numerical algorithm based on discrete three-dimensional spatiotemporal schemes optimized for parallel computations. Particularly, the kinetic D3Q19 latticeBoltzmann equation method (LBM)36,37 and the D3Q19 lattice-based approaches developed by Warren38 and by Capuani et al.39 were employed to resolve numerically the NavierStokes, Poisson, and NernstPlanck problems, respectively. The D3Q19 lattice is a projection of the four-dimensional face-centered hypercubic lattice onto three-dimensional space, where each lattice node is connected with 18 nearest and nextnearest neighbors. Such lattice connectivity is indispensable to achieve isotropic and Galilean-invariant transport behavior of the modeled system. The body force was introduced into the LBM by the method of calculating the equilibrium distribution function with an altered velocity.40,41 The no-slip boundary condition at the solidliquid interface for the NavierStokes equation was realized by the halfway bounce-back rule,42 whereas conventional Neumann boundary conditions assuming constant surface charge density and zero normal flux at the solidliquid interface were implemented for the Poisson and NernstPlanck problems, respectively. All numerical schemes were realized as parallel codes in C language using the message passing interface (MPI) standard. The developed computational model was implemented at a supercomputer (SuperMIG) of the LeibnizRechenzentrum der Bayerischen Akademie der Wissenschaften (Garching, Germany). A typical simulation of masscharge transport in the funnel system required ∼30 h at 1024 processor cores and around 70 GB of system memory. In all numerical schemes, space and time steps of Δx = 2.5 nm and Δt = 2 1011 s were used. To guarantee that the above spatiotemporal resolution is sufficient to perform accurate modeling of masscharge transport, we compared the fluid flow velocity and electric potential fields computed for geometrically similar, but physically smaller systems with various resolutions. For example, by using a still finer resolution of Δx = 1.25 nm and Δt = 0.5 1011 s, the computed flow velocity and potential fields were found to differ only slightly and locally from those obtained with Δx = 2.5 nm and Δt = 2 1011 s, by less than 4%. During the simulations, we assumed that the nanofluidic device is filled with a 1 mM potassium chloride electrolyte solution (cbulk)
Figure 1. (a) Atomic force microscope image of the nanofunnel with a 5° taper angle. (b) Layout of the integrated microchannelnanochannel device. (c) Scanning electron microscope image of the 5° funnel. (d) Schematic of the modeled nanofluidic channelfunnelchannel configuration with a 5° taper angle. The height of the modeled device is 120 nm.
Ji ¼ Di ∇ci
Z
)
)
where n is the unit normal vector directed from the liquid toward the solid phase, σ is the surface charge density, r ci,in and r ci,out are the inlet and outlet axial components of the ith species concentration gradient, respectively, and A and B are constant values corresponding to the voltage and pressure differences applied to the system, respectively. The resulting ionic current was determined by integration of the flux over a cross-sectional area S of the system:
ð3Þ
In eqs 13, Ji, Di, ci, and zi are the ionic flux, diffusion coefficient, concentration (mol/m3), and algebraic valence of species i, respectively, ϕ and u are the local electric potential and fluid flow velocity, ε, F, and η are the permittivity, density, and viscosity of the fluid, respectively, p is the pressure, and F, R, and T denote the Faraday constant, molar gas constant, and absolute temperature, respectively. The NernstPlanck equation (eq 1) describes local mass and charge transport due to diffusion, electromigration, and advection. The Poisson equation (eq 2) establishes a relationship between the local electric potential and the species concentration distributions. Finally, the generalized NavierStokes equation (eq 3) relates the fluid 269
dx.doi.org/10.1021/ac202501v |Anal. Chem. 2012, 84, 267–274
Analytical Chemistry
ARTICLE
potential distribution obtained for the funnel base shown in Figure 2b. Due to the lateral symmetry of the potential surface, we present only half of the cross-section to keep the same scale as in Figure 2a. In Figure 2c, the one-dimensional potential profile extracted from the potential surface in Figure 2b is shown. This profile was obtained for the direction parallel to the device height (120 nm) in the center of the funnel base. The corresponding surface potential determined at that location is 56.36 mV, which is very close to the value of the surface potential, ζ = 56.77 mV, calculated by the Grahame equation:43 " # 2RT σ sinh1 ζ¼ ð7Þ F ð8εε0 RTcbulk Þ1=2 We conclude that the employed numerical approach is able to reproduce accurately the EDL formation in an electrolyte solution due to the presence of surface charges. Another relevant conclusion is that both transverse dimensions must be accounted for by an adequate model of the system. Although the straight nanochannels and the funnel base are rather wide (1.05 μm) compared to the EDL thickness (10 nm), the device height (120 nm) is comparable to λD and a significant portion of any cross-section of the investigated system is occupied by the EDLs developed at the substrate and the cover (Figure 2). Hence, an accurate description of the modeled system cannot be obtained with the assumption that the system is homogeneous along the third dimension. Thus, an adequate model must be three-dimensional to reproduce all relevant features of the system. Next we studied transport properties of the modeled system, in particular, its ability to rectify ionic current. This feature was characterized by the rectification ratio determined as the ratio of the resulting ionic currents through the system for biases with the same magnitude but opposite sign. We designated these two regimes as the high-conductance (HG) and low-conductance (LG) regimes. According to the system geometry presented in the schematic in Figure 1d, the counter electrode was assumed to be grounded, and the HG and LG regimes were realized by setting the potential of the working electrode jwork = (0.4 V to be, correspondingly, negative and positive. The above absolute value of jwork results in an average applied electric field strength across the system of ∼11.8 kV/m, which is close to that used for the experimental setup when a 1 V bias was applied to the ∼85 μm long system.33 The simulations of ion transport under the aforementioned conditions result in electric currents of 73.0 and 50.7 pA for the HG and LG regimes, respectively, providing a rectification ratio of 1.44. The corresponding experimental values are 105 pA, 50 pA, and 1.9. Eventual causes for the disagreement between simulated and experimental data include differences in the surface charge density, electrolyte solution, physical parameters, and system geometry. However, our primary goal during the present work was to reveal the origin of ICR in the system we studied, but not to reproduce the exact experimental values. Therefore, we did not intend to fit the model parameters to get a perfect quantitative agreement between simulations and experimental results. Then we analyzed the simulated distributions of the species concentrations in the system. A conventional explanation of ICR in conical pores is the formation of enriched and depleted regions within the pore.9,12,14,23 According to this mechanism, ICR is a result of the selectivity and, as a consequence, a large transference number for counterions at the pore tip. In turn, rectification is a consequence of an excess of mobile counterions in that region. If
Figure 2. Simulated cross-sectional equilibrium potential distributions for the (a) funnel tip (width 80 nm, height 120 nm) and (b) funnel base (width 1050 nm, height 120 nm). For the base cross-section only half of the whole area is shown. (c) Distribution of electric potential for the direction parallel to the device height in the center of the funnel base. The surface charge density is 5 mC/m2, and a 1 mM 1:1 electrolyte solution is assumed.
characterized by η = 103 Pa 3 s, ε = 7.0 1010 F/m, F = 103 kg/m3, and T = 293 K. The diffusion coefficients for potassium and chloride ions were assumed to be 2 109 m2/s.28 The surface charge density at the solidliquid interface is equal to 5 mC/m2.
’ RESULTS AND DISCUSSION Before presenting and discussing simulation results for mass and charge transport in the modeled system, we consider and analyze equilibrium electric potential and species concentration distributions obtained with the conditions ϕin = ϕout and pin = pout, i.e., when the distributions are governed exclusively by the surface charge at the solidliquid interface and the system geometry. In parts a and b of Figure 2, we present the simulated electric potential distributions for cross sections corresponding to the funnel tip and funnel base, respectively. For a 1 mM 1:1 electrolyte solution, the Debye screening length characterizing the EDL thickness is λD ≈ 10 nm. Therefore, the EDLs formed at the walls in the tip region extend significantly into the nanofunnel as compared to its local transverse dimensions of 80 nm 120 nm (Figure 2a). At the same time, although the width of the funnel increases gradually from 80 nm (tip) to 1050 nm (base), the funnel height remains constant at 120 nm, resulting in a still significant extension (relative to the system height) of the EDLs associated with the substrate and cover surfaces into the central region of the system cross-section. The degree of EDL extension is illustrated by the cross-sectional 270
dx.doi.org/10.1021/ac202501v |Anal. Chem. 2012, 84, 267–274
Analytical Chemistry
ARTICLE
Figure 3. Simulated normalized concentrations of cations and anions along the centerline of the system with a 5° funnel (filled with 1 mM potassium chloride solution) for the HG and LG regimes (surface charge density 5 mC/m2, applied electric field strength (11.765 kV/m). Dashed lines indicate the location of the funnel tip and base.
Figure 4. Bright-field image of the 5° funnel and nanochannels connecting two microchannels (a) and fluorescence images of disodium fluorescein transported through the device in the HG regime (b) and LG regime (c). The inset in panel a is a schematic of the funnel with the rectangle highlighting its position along the 85 μm long nanochannel.
region providing ionic flux imbalance due to EDLs that occupy a significant fraction of the tip cross-section. Consequently, depletion and enrichment zones develop, respectively, in the anodic and cathodic regions adjacent to the funnel tip. As a result, extended depletion and enrichment zones develop in the counter electrode channel in the HG and LG regimes, respectively (cf. Figure 3). Another relevant observation is that the enriched CP zone which develops in the HG regime within the funnel is extended over the whole working electrode channel. Concentrations of both cations and anions in that compartment are ∼1.12cbulk. The above results can be compared with those obtained by monitoring the fluorescence intensity of an anionic probe (disodium fluorescein) presented in Figure 4.33 In the HG regime (Figure 4b), fluorescein accumulates in the whole working electrode channel, whereas the counter electrode channel is characterized by a lower concentration compared to when no bias is applied. The difference in the ionic species concentrations leads immediately to a difference in the local conductivity and total system conductance. In Figure 5, we present the local conductivity of the electrolyte solution along the centerline of the system, normalized by the conductivity corresponding to cbulk. Similar to the data presented in Figure 3 for ionic concentration distributions, the whole working electrode channel is characterized by an increased conductivity in the HG regime. Thus, the geometrical asymmetry of the nanofunnel results in different ionic species concentrations in the HG and LG regimes not only within the funnel, but also in the adjoining straight nanochannels. The above finding can be explained by propagating CP zones. Recently, a two-paper series concerning the theoretical and experimental study of ion transport in variableheight micro- and nanochannels and the dynamics of CP was published.34,35 Two distinct physical regimes were identified: (1) CP without propagation where the polarization effect stays local to the micronanochannel interface and (2) CP with propagation in which enrichment and depletion shocks propagate outward. The existence of each regime was found to depend on two nondimensional system parameters: a nanochannel Dukhin
the potential of the pore base is lower than that of the pore tip, electromigration flux of co-ions is directed from the base of the pore toward its tip. By approaching the pore tip, co-ions are rejected because the co-ion transference number becomes small. In contrast, counterions migrating toward the tip from the electrolyte solution in the external reservoir are free to traverse the tip. As a result, a region of locally high ionic concentration builds up inside the pore near its tip. The electrolyte conductivity in this region increases, and the system enters the HG regime. If the potential of the pore base is higher than that of the tip, electric transference causes a decrease of ionic concentration within the pore and the system enters the LG regime. Summarizing the above illustrative explanation, the ionic concentration within the pore in the HG regime must be higher than that in the LG regime. Our simulations indicate that the above explanation is generally applicable as well to explain ICR in the funnel systems. The results of our simulations were obtained as three-dimensional data structures, and we used their complete content for further analysis. However, three-dimensional data structures are difficult to represent clearly. Therefore, for visualization purposes, we used reduced-dimensional extractions from the complete data structures, which permits characterization of the basic properties of the system. In Figure 3 simulated concentration profiles along the centerline of the studied nanochannelnanofunnel system are shown. Dotted lines indicate the location of the funnel tip and funnel base. We normalized local concentrations with cbulk = 1 mM. The data presented in this figure indicate that the local species concentrations within the funnel in the HG regime are higher than those in the LG regime. In the HG regime, these concentrations are higher than cbulk, whereas in the LG regime they are lower than cbulk. Additionally, depletion and enrichment zones develop in the region of the counter electrode channel adjacent to the funnel tip in the HG and LG regimes, respectively. This picture resembles CP developing around a charge-selective domain, e.g., an ion-exchange membrane44,45 or an ion-permselective particle46 due to a difference in the ion transport numbers. This similarity can be explained by the fact that the funnel tip is a 271
dx.doi.org/10.1021/ac202501v |Anal. Chem. 2012, 84, 267–274
Analytical Chemistry
ARTICLE
Figure 5. Simulated normalized conductivity of the electrolyte solution along the centerline of the system with a 5° funnel (filled with 1 mM potassium chloride solution) for the HG and LG regimes (surface charge density 5 mC/m2, applied electric field strength ( 11.765 kV/m). Dashed lines indicate the location of the funnel tip and base.
Figure 6. (a) Schematic of a linear microchannelnanochannel microchannel system with developed CP zones and its equivalent electric circuit composed of three variable resistors representing the resistance of the anodic and cathodic regions and of the charge-selective domain (nanochannel). (b) Schematic of a linear nanochannelnanofunnel nanochannel system and its equivalent electric circuit.
number (ratio of the EDL conductivity to the bulk conductivity) and electrophoretic mobility of the co-ions divided by the EOF mobility. Thus, a nanochannel (or any other charge-selective spatial domain) can have significant long-range effects in nanofluidic systems (up to ∼cm long). Moreover, the steady-state concentrations of co-ions and counterions in the propagating enrichment and depletion zones are as well governed by the aforementioned system parameters; in particular, they depend on the EOF velocity. In a subsequent analysis,47 the data from 56 sets of experiments reported in 19 publications were analyzed according to that approach and indicated that the model is a good predictor of the above two different physical regimes. Hence, an equivalent circuit for the electric conductance of microchannelnanochannel systems can be presented in the form of three series-connected variable resistors. For geometrically symmetric systems (Figure 6a) their resistances do not depend on the polarity of the applied voltage; they are determined exclusively by the inherent system parameters, e.g., system geometry, surface charge density, electrolyte ionic strength, and ion mobilities. For systems with asymmetry (Figure 6b), these values are additionally affected by the polarity of the externally applied voltage. The differentiation between the two physical regimes of nonpropagating and propagating CP zones in a microchannel nanochannel system was proposed according to the following condition (assuming a 1:1 electrolyte solution):34,35 cbulk
Fh < maxðν, 2ν 1Þ σ
We applied a similar approach to evaluate the HG and LG regimes in the modeled system. Here, the left-hand side (lhs) of eq 8 is determined only by quantities characterizing the physicochemical and geometric system properties, which are not affected by the masscharge transport characteristics. This parameter combination is equivalent to the inverse Dukhin number. By contrast, the right-hand side (rhs) of eq 8 includes the values of parameters that quantitatively characterize mass and charge transport. The width of the funnel tip (80 nm) was used as the characteristic length (h) describing the transverse dimension of a charge-selective domain, resulting in 1.544 as a value for the lhs of eq 8. Then we used the simulation results to calculate the values of uep and ueof (cf. eq 9) averaged over the tip crosssectional area for the HG and LG regimes. Our simulation results gave a different EOF rate (or superficial flow velocity) for these two regimes. The ratio of the average EOF velocities for jwork = (0.4 V was 1.29. However, the higher EOF rate was observed in the LG regime, in contrast to the ICR. Similar behavior was observed with simulations of ICR in conical nanopores reported by White and Bund23 and in experiments carried out with pyramidal-pore mica membranes;8 the ionic current and cathodic EOF are rectified in opposite manners. The theoretical explanation of this phenomenon is beyond the scope of the current paper, and we refer the interested reader to a related publication.48 Here we want to point out that the difference in the EOF velocity fields originates from a difference in induced back-pressure flows. Figure 7 shows the distributions of the axial component of the EOF velocity in the middle plane of the system (i.e., in the plane at channel half-height) for the HG and LG regimes. To facilitate a comparison of the two velocity fields, we assume that the positive velocity direction coincides with that of the average EOF for the corresponding regime (in fact, these directions are opposite for the HG and LG regimes). In the HG regime, recirculation flow domains (the blue-colored domains in the upper panel of Figure 7) develop in the region of the counter electrode channel, which is adjacent to the funnel tip. This region is characterized by an abrupt reduction in the channel cross-section that results in an
ð8Þ
with ν being the dimensionless mobility of the co-ions, defined as uep ν¼ ð9Þ ueof where h is the nanochannel height (or smaller lateral dimension), uep is the electrophoretic velocity of co-ions, and ueof is the EOF velocity in the nanochannel. The enrichment and depletion zones propagate outward from a charge-selective domain if the inequality expressed by eq 8 is satisfied; otherwise, the CP zones are stationary and remain localized with respect to that domain. 272
dx.doi.org/10.1021/ac202501v |Anal. Chem. 2012, 84, 267–274
Analytical Chemistry
ARTICLE
positioned between two straight nanochannels with the same height (120 nm). Thus, a significant portion of any cross-section of the investigated system is occupied by electric double layers, which requires a three-dimensional modeling approach. The results of the simulations indicate that the ICR phenomenon in the modeled system can be explained by a mechanism similar to that observed in conical nanopores; i.e., the HG regime is characterized by increased concentrations of ionic species within the geometrically asymmetric domain (funnel), whereas species concentrations for the LG regime are reduced. However, the presence of the adjoining straight nanochannels permits the developed enriched and depleted CP zones to be extended over significant distances beyond the funnel segment. The funnel tip is a region characterized by charge-selective properties, and the development of an enriched/depleted CP zone at one side of the tip is accompanied by the formation of a depleted/enriched zone at the other side. Moreover, the HG regime is characterized by propagating CP zones, whereas CP zones remain local in the LG regime. The presented numerical simulation approach has identified propagating CP as an essential component during ICR with the investigated nanochannelnanofunnelnanochannel system. The approach is also numerically efficient, which will be helpful in the study of nanofluidic ICR devices under general conditions when a three-dimensional analysis of physical phenomena that govern flow and transport needs to be performed, i.e., when the inherent device geometry and/or the evolving (electro)hydrodynamics require the full three-dimensional information to be held on, e.g., species concentrations, electric potential, and flow velocity.
Figure 7. Distribution of the axial flow velocity in the middle plane of the system with a 5° funnel (filled with 1 mM potassium chloride solution) simulated for the HG and LG regimes (surface charge density 5 mC/m2, applied electric field strength ( 11.765 kV/m). The positive velocity direction is assumed to coincide with the average EOF for each regime.
induced pressure gradient due to the geometrical constriction and in an increase of the inertial force contribution to the local hydrodynamics of the modeled system. In contrast, no remarkable flow recirculation develops in the LG regime. Asymmetric EOF patterns in systems with a converging and diverging microchannel were observed as well experimentally by Wang et al.7 The formation of asymmetric flow patterns in systems with an asymmetric configuration can be explained by different conditions for the mutual transformation of kinetic energy associated with the fluid flow velocity and potential energy associated with the induced pressure gradient.49 Returning to our analysis, the maxima of the rhs of eq 8 for the HG and LG regimes calculated from the simulation results are 5.18 and 0.48, respectively. According to eq 8, the HG regime should be characterized by propagating CP zones, whereas for the LG regime the CP zones must be stationary and localized around the funnel. The above regime differentiation expressed by eqs 8 and 9 was developed for systems where a charge-selective domain is assumed to be axially uniform. The capability of the funnel to differentiate transport of co-ions and counterions decreases gradually from its tip toward the base region. Moreover, the increasing width of the funnel exceeds its constant height at some distance from the tip. In that region the chargeselective properties arise mainly due to a restriction along the third direction (channel height 120 nm). Therefore, one can expect that the inverse Dukhin number based on the tip width (80 nm) is lower than that determined by some average transverse characteristic length of the funnel, which should be less than 120 nm. Nevertheless, the maxima of the rhs in eq 8 for the HG and LG regimes differ by 1 order, and this indicates that two different regimes can be well distinguished with the presented approach.
’ AUTHOR INFORMATION Corresponding Author
*E-mail: tallarek@staff.uni-marburg.de (U.T.); jacobson@ indiana.edu (S.C.J.).
’ ACKNOWLEDGMENT This work was supported in part by the National Science Foundation (Grants CHE-0750295 and CHE-0832651). Numerical simulations were run at the Leibniz-Rechenzentrum der Bayerischen Akademie der Wissenschaften (Garching, Germany), supported by Project HLRB pr26wo. ’ REFERENCES (1) Siwy, Z.; Kosi nska, I. D.; Fuli nski, A.; Martin, C. R. Phys. Rev. Lett. 2005, 94, 048102. (2) Berezhkovskii, A. M.; Pustovoit, M. A.; Bezrukov, S. M. J. Chem. Phys. 2007, 126, 134706. (3) Berezhkovskii, A. M.; Pustovoit, M. A.; Bezrukov, S. M. Phys. Rev. E 2009, 80, 020904(R). (4) Berezhkovskii, A. M.; Pustovoit, M. A.; Bezrukov, S. M. Chem. Phys. 2010, 375, 523–528. (5) Olsson, A.; Stemme, G.; Stemme, E. J. Micromech. Microeng. 1999, 9, 34–44. (6) Hwang, J. J.; Tseng, F. G.; Pan, C. Int. J. Multiphase Flow 2005, 31, 548–570. (7) Wang, S.; Hu, X.; Lee, L. J. Lab Chip 2008, 8, 573–581. (8) Jin, P.; Mukaibo, H.; Horne, L. P.; Bishop, G. W.; Martin, C. R. J. Am. Chem. Soc. 2010, 132, 2118–2119. (9) Wei, C.; Bard, A. J.; Feldberg, S. W. Anal. Chem. 1997, 69, 4627–4633.
’ CONCLUSIONS In this paper we report the first three-dimensional numerical simulations of coupled mass and charge transport in a nanofluidic funnel system. The simulated system is composed of a nanofunnel 273
dx.doi.org/10.1021/ac202501v |Anal. Chem. 2012, 84, 267–274
Analytical Chemistry
ARTICLE
(10) Siwy, Z.; Heins, E.; Harrell, C. C.; Kohli, P.; Martin, C. R. J. Am. Chem. Soc. 2004, 126, 10850–10851. (11) Siwy, Z.; Fuli nski, A. Am. J. Phys. 2004, 72, 567–574. (12) Siwy, Z. S. Adv. Funct. Mater. 2006, 16, 735–746. (13) Cervera, J.; Schiedt, B.; Neumann, R.; Mafe, S.; Ramírez, P. J. Chem. Phys. 2006, 124, 104706. (14) Cervera, J.; Alcaraz, A.; Schiedt, B.; Neumann, R.; Ramírez, P. J. Phys. Chem. C 2007, 111, 12265–12273. (15) Kovarik, M. L.; Zhou, K.; Jacobson, S. C. J. Phys. Chem. B 2009, 113, 15960–15966. (16) Cervera, J.; Ramírez, P.; Mafe, S.; Stroeve, P. Electrochim. Acta 2011, 56, 4504–4511. (17) Woermann, D. Nucl. Instrum. Methods Phys. Res., B 2002, 194, 458–462. (18) Woermann, D. Phys. Chem. Chem. Phys. 2003, 5, 1853–1858. (19) Woermann, D. Phys. Chem. Chem. Phys. 2004, 6, 3130–3132. (20) Siwy, Z.; Fuli nski, A. Phys. Rev. Lett. 2002, 89, 198103. (21) Cervera, J.; Schiedt, B.; Ramírez, P. Europhys. Lett. 2005, 71, 35–41. (22) Liu, Q.; Wang, Y.; Guo, W.; Ji, H.; Xue, J.; Quyang, Q. Phys. Rev. E 2007, 75, 051201. (23) White, H. S.; Bund, A. Langmuir 2008, 24, 2212–2218. (24) Ramírez, P.; Apel, P. Y.; Cervera, J.; Mafe, S. Nanotechnology 2008, 19, 315707. (25) Qian, S.; Joo, S. W.; Cheney, M. A.; Hou, W. J. Colloid Interface Sci. 2009, 329, 376–383. (26) Ai, Y.; Zhang, M.; Joo, S. W.; Cheney, M. A.; Qian, S. J. Phys. Chem. C 2010, 114, 3883–3890. (27) Momotenko, D.; Cortes-Salazar, F.; Josserand, J.; Liu, S.; Shao, Y.; Girault, H. H. Phys. Chem. Chem. Phys. 2011, 13, 5430–5440. (28) Kubeil, C.; Bund, A. J. Phys. Chem. C 2011, 115, 7866–7873. (29) Wan, W.-J.; Holden, D. A.; White, H. S. J. Am. Chem. Soc. 2011, 133, 13300–13303. (30) Jung, J.-Y.; Joshi, P.; Petrossian, L.; Thornton, T. J.; Posner, J. D. Anal. Chem. 2009, 81, 3128–3133. (31) H€oltzel, A.; Tallarek, U. J. Sep. Sci. 2007, 30, 1398–1419. (32) Kim, S. J.; Song, Y.-A.; Han, J. Chem. Soc. Rev. 2010, 39, 912–922. (33) Perry, J. M.; Zhou, K.; Harms, Z. D.; Jacobson, S. C. ACS Nano 2010, 4, 3897–3902. (34) Mani, A.; Zangle, T. A.; Santiago, J. G. Langmuir 2009, 25, 3898–3908. (35) Zangle, T. A.; Mani, A.; Santiago, J. G. Langmuir 2009, 25, 3909–3916. (36) Chen, S.; Doolen, G. D. Annu. Rev. Fluid Mech. 1998, 30, 329–364. (37) Succi, S. The Lattice Boltzmann Equation: For Fluid Dynamics and Beyond; Oxford University Press: Oxford, U.K., 2001. (38) Warren, P. B. Int. J. Mod. Phys. C 1997, 8, 889–898. (39) Capuani, F.; Pagonabarraga, I.; Frenkel, D. J. Chem. Phys. 2004, 121, 973–986. (40) Shan, X.; Doolen, G. J. Stat. Phys. 1995, 81, 379–393. (41) Martys, N. S.; Chen, H. Phys. Rev. E 1996, 53, 743–750. (42) Gallivan, M. A.; Noble, D. R.; Georgiadis, J. C.; Buckius, R. O. Int. J. Numer. Methods Fluids 1997, 25, 249–263. (43) Grahame, D. C. J. Phys. Chem. 1953, 21, 1054–1060. (44) Hlushkou, D.; Dhopeshwarkar, R.; Crooks, R. M.; Tallarek, U. Lab Chip 2008, 8, 1153–1162. (45) Kim, P.; Kim, S. J.; Han, J.; Suh, K. Y. Nano Lett. 2010, 10, 16–23. (46) Ehlert, S.; Hlushkou, D.; Tallarek, U. Microfluid. Nanofluid. 2008, 4, 471–487. (47) Zangle, T. A.; Mani, A.; Santiago, J. G. Chem. Soc. Rev. 2010, 39, 1014–1035. (48) Dorrepaal, J. M. J. Eng. Math. 1993, 27, 343–356. (49) Olsson, A.; Stemme, G.; Stemme, E. Sens. Actuators, A 2000, 84, 165–175. 274
dx.doi.org/10.1021/ac202501v |Anal. Chem. 2012, 84, 267–274
