Steady Electrodiffusion in Hydrogel-colloid Composites: Macroscale Properties from Microscale Electrokinetics

A rigorous microscale electrokinetic model for hydrogel-colloid composites is adopted to compute macroscale profiles of electrolyte concentration, electrostatic potential, and hydrostatic pressure across membranes that separate electrolytes with different concentrations. The membranes are uncharged poly-meric hydrogels in which charged spherical colloidal particles are immobilized and randomly dispersed with a low solid volume fraction. Bulk membrane characteristics and performance are calculated from a continuum microscale electrokinetic model (Hill 2006b, c). The computations undertaken in this paper quantify the streaming and membrane potentials. For the membrane potential, increasing the volume fraction of negatively charged inclusions decreases the differential electrostatic potential across the membrane under conditions where there is zero convective flow and zero electrical current. With low electrolyte concentration and highly charged nanoparticles, the membrane potential is very sensitive to the particle volume fraction. Accordingly, the membrane potential – and changes brought about by the inclusion size, charge and concentration – could be a useful experimental diagnostic to complement more recent applications of the microscale electrokinetic model for electrical microrheology and electroacoustics (Hill and Ostoja-Starzewski 2008, Wang and Hill 2008).


INTRODUCTION
Hydrogels are an important class of soft matter that have gained widespread application in drug delivery (Qiu and Park 2001, Lin and Netters 2006, Peppas et al. 2000), tissue engineering (Khademhosseini and Langer 2007, Barndl et al. 2007, Drury andMooney 2003), advanced materials (Peppas et al. 2006, Eddington and Beebe 2004, Chaterji et al. 2007), and molecular separations (Wang et al. 1993, Kim andPark 1998).Novel characteristics can be achieved by immobilizing organic and inorganic colloidal particulates in the polymer skeleton.For example, embedding gold or gold-coated silica nanoparticles into a thermally responsive hydrogel induces light-wavelength-sensitive swelling to achieve optically active microfluidic flow control (Sershen et al. 2005).In biosensing, immobilizing silica nanoparticles in polyacrylamide hydrogels and applying an electric field increase the otherwise diffusion-limited flux of uncharged macromolecules across the composite membrane (Matos et al. 2006).Other applications include delivering growth factors for bone regeneration (Chung et al. 2007), improving the contrast of ultrasound imaging for early tumor detection (Liu et al. 2006, Dayton andFerrara 2002), and absorbing infrared energy for certain cancer treatment (Loo et al. 2005).Note also that polystyrene nanoparticles have been dispersed in neutral polyacrylamide hydrogels to increase the storage modulus and produce mechanoelectrical effects for artificial tactile perception and psycho-sensorial materials (Thévenot et al. 2007).Hill (2006c, b) developed an electrokinetic transport model to quantify how imposed gradients of electrostatic potential, ion concentration, and pressure perturb an equilibrium state where each colloidal inclusion in the hydrogel-colloid composite is enveloped by a diffuse layer of counterions.Such an equilibrium is widely acknowledged to be reasonably well described by the non-linear Poisson-Boltzmann equation (Verwey and Overbeek 1948), which itself is a special case of a much more general electrokinetic transport (non-equilibrium) model (e.g., Overbeek 1943, Booth 1950).Similar methodologies have been adopted for ordered and random consolidated porous media with the immobile charge uniformly distributed on the solid matrix, often under conditions where the Debye length is smaller than the characteristic pore size (e.g., Gupta et al. 2007, Wang andChen 2007).
Hill calculated perturbations to the equilibrium state of a single charged sphere immobilized in an uncharged porous medium (polymer hydrogel), and averaged the governing microscale equations to derive macroscale equations for dilute random dispersions, with transport coefficients derived from the microscale analysis of a single sphere.In principle, this theory could be advanced to handle higher particle concentrations by adopting cell models to account for particle interactions (e.g., Ahualli et al. 2006).Accordingly, the present macroscale equations are appropriate only when the inclusion volume fraction is sufficiently small, which, nevertheless, may often be the case.It should also be noted that the polymer skeleton is uncharged, so the immobile charge in the composite arises solely from the colloidal inclusions.Nevertheless, the calculations are valid for all practical values of the electrolyte concentration (Debye length), and particle size and surface charge density, which together determine the ζ -potential.
The hydrogel phase is modeled as a porous (Brinkman) medium with a very low volume fraction, so it hinders fluid flow due to hydrodynamic (Darcy) drag and viscous stresses, but does not hinder ion diffusion and electromigration.Predictions of the electroosmotic pumping capacity (termed the incremental pore mobility) and its relation to chemical and physical characteristics of the hydrogel and inclusions were recently found to compare well with experiments involving silica nanoparticles embedded in polyacrylamide (Hill 2007, Matos et al. 2006).Noteworthy is that the theoretical interpretation suggests that particle-particle and particle-polymer interactions (in the experiments) increase the effective hydrodynamic permeability of the hydrogel (Hill 2007).
The present model does not address elastic deformation, which is likely to be significant in situations where flow is driven by a sufficiently large applied pressure gradient (e.g., during streaming potential measurements), or the composite has to accommodate large gradients in the bulk concentration of electrolyte.However, recent advances of this model quantify how the electrolyte and inclusions couple (mechani- cally, electrically, and hydrodynamically) with the elastic polymer network (Hill 2006a, Hill and Ostoja-Starzewski 2008, Wang and Hill 2008), thereby facilitating calculations of the particle displacements and velocity when hydrogel-colloid composites are subjected to steady and oscillatory electric fields.Such experiments are relevant to the fields of microrheolgy (Cicuta and Donald 2007, Furst 2005, MacKintosh and Schmidt 1999) and electroacoustics (Hunter 1998, O'Brien 1990, 1988), which are widely used to probe the microstructure of complex fluids, including colloidal dispersions, polymer solutions and gels.
This paper establishes a systematic methodology for calculating macroscale transport properties of hydrogel-colloid composites based on earlier microscale analysis in the literature (Hill 2006c, b).These macroscale calculations are expected to complement experimental measurements of ion fluxes and membrane potentials.Accordingly, the paper begins by considering the continuum equations for a simple electrolyte in the absence of inclusions and hydrogel.This analysis of the well-known Poisson equation and ion-conservation equations is then used to rationalize an approximation (often referred to as an assumption) of bulk electroneutrality.The accompanying exact solution of this simple model provides a useful limiting case for interpreting subsequent numerical solutions of the non-linear coupled macroscale equations for hydrogel-colloid composite membranes: namely Hill's averaged ion-conservation equations, and fluid mass and momentum conservation equations.Several striking aspects of the results are physically interpreted using previous knowledge of the underlying microscale electrokinetic phenomena.

ELECTRODIFFUSION IN THE ABSENCE OF POLYMER AND INCLUSIONS
Here we consider electrodiffusion in an electrolyte without polymer or inclusions.Accordingly, the ions migrate by convection, diffusion and electromigration.If the fluid velocity is independently specified, then the equations required to determine the electrostatic potential ψ and N ion concentrations n j (z j is the valence and D j the diffusion coefficient with j = 1 . . .N ) under steady conditions are the well-known Poisson and ion-conservation equations: If the average fluid velocity u is specified, then this is a closed system of N + 1 equations involving E = −∇ψ and n j as unknowns ( o is the vacuum permittivity, kT is the thermal energy, and e is the fundamental charge).Note that electrical neutrality is not imposed explicitly, since it should emerge naturally from the boundary conditions.For a one-dimensional problem (0 ≤ x ≤ L) with uniform dielectric constant s , the equations above are (3) where subscripts x denote differentiation, and j j = const.is the flux of the jth species.
An Acad Bras Cienc (2010) 82 (1) Note that the electrical current density is so under conditions where i = 0 with a z-z electrolyte, j 1 = j 2 = j.Alternatively, consider solving If the average fluid velocity u is specified, then a closed system of N + 1 equations must be solved with ψ and n j as unknowns.However, while these equations could be solved with prescribed boundary conditions (at each end of the domain), the electrical current cannot be set to zero a priori.Rather, the electrostatic potential difference across the domain that yields zero electrical current must be established.Even so, a numerical solution may still be challenging because, as demonstrated below, the left-hand side of the Poisson equation above is extremely small compared to the right-hand side when the characteristic length scale is larger than the Debye length (typically between 1 and 100 nm).

ELECTRONEUTRALITY
Consider scaling the Poisson equation using a macroscopic length scale l c , with electrostatic potential scale ψ c = kT /e and ion-concentration (charge density) scale n c .Accordingly, the left-hand side is O[(κl c ) −2 ] with respect to the right-hand side, with κ −1 ∼ [ o s kT /(e 2 n c )] 1/2 the Debye length.Since κ −1 is of nanometer scale for charge densities in aqueous electrolytes, on macroscopic scales (κl c ) −2 1.Therefore, in the 'outer' region where l c κ −1 is indeed the appropriate length scale, the leading-order approximation of the Poisson equation becomes a statement of local electroneutrality, so the governing equations above become 0 = N j=1 n j z j e, (8) This important result was first established formally by MacGillivray (1968) using matched asymptotic expansions, and it is now widely adopted in standard texts.The following example establishes an explicit formula to be adopted in the following sections that address the much more difficult problem where charged colloidal inclusions are immobilized in an uncharged hydrogel skeleton.For a z-z electrolyte, the bulk electroneutrality approximation requires n 1 = n 2 = n, so the ionconservation equations are An Acad Bras Cienc (2010) 82 (1) or Eliminating E gives which is easily solved for n(x) given a constant value of u.
When the current density i = 0, j 1 = j 2 = j and Note that with zero convective flux (u = 0), and where γ = D 1 /D 2 and z 1 = −z 2 .Moreover, the potential difference across the membrane (membrane potential) is Let us examine the validity of the electroneutrality approximation by evaluating the left-hand side of the Poisson equation, which, recall, was neglected in reaching Eqn.(20): This charge density is equivalent to a molar concentration of ions, each with charge It is easily verified that this concentration is generally extremely small unless the characteristic length scale is O(nm).For example, for KNO 3 , γ = D + /D − ≈ 73.5/71.46≈ 1.029.Therefore, with n(0) = 1 mmol l −1 and n(L) = 100 mmol l −1 , Eqn. ( 22) gives a maximum at x = 0 of approximately 2 × 10 −6 mmol l −1 when L = 1 mm.In general, therefore, and consistent with MacGillivray (1968), the electroneutrality approximation is reasonable when the characteristic length scale is greater than the Debye length.

ELECTRODIFFUSION IN HYDROGEL-COLLOID COMPOSITE MEMBRANES
The averaged equations derived by Hill (2006b) are written below for unidirectional transport of a z-z electrolyte in hydrogel composites.Note that the electroneutrality approximation requires n 1 = n 2 = n , and the fluid conservation (continuity) equation requires u = const., so the remaining ion conservation equations and momentum equation are and These are more compactly written where Note that the asymptotic coefficients: provides an algebraic relationship between the two bulk ion fluxes.For the example below, however, there is zero electrical current, so j 1 = j 2 = j = const.It is helpful to write the macroscale equations above in a form that explicitly relates the unknown gradients of n , ψ and p to the (co-linear) constants j 1 , j 2 and u , i.e.,    where The components of matrix B depend on n (x), so this non-linear system of equations must be solved for y(x) = [ n (x), ψ (x), p (x)], given a constant b = ( j 1 , j 2 , u ).This is easily achieved by numerically integrating1 dy dx In general, such a computation is rather intensive, since, in addition to the bulk ion concentration, the asymptotic coefficients vary with κa and ζ e/(kT ) (with fixed /a), both of which are non-linear functions of n (x).However, the asymptotic coefficients can be computed beforehand, and subsequently interpolated from sufficiently refined tables.Representative values (and several other quantities derived from them) are provided in Tables I, II and III for NaCl electrolyte with /a = 0.1 (a = 10 nm and = 1 nm).For the macroscale computations presented below, the asymptotic coefficients were available for many more values of κa and ζ e/(kT ).The reader is referred to Hill (2006a, c) for details of these microscale calculations and the physical significance of the other quantities provided2 .Note that it is straightforward to permit the inclusions to have a constant charge, in which case the ζ -potential varies with position according to the particle size and surface charge density, and the bulk electrolyte concentration (Verwey andOverbeek 1948, Russel et al. 1989).
The hydrogels for the calculations presented below have a Brinkman screening length = 1 nm, which is representative of polyacrylamide (Hill 2007).The inclusion radius a = 10 nm, and the electrolyte (NaCl) is moderately asymmetric with γ = D 1 /D 2 ≈ 1.33/2.33.Note that the dimensionless parameter /a = 0.1, whereas previous microscale calculations reported in the literature for NaCl have /a ≈ 0.01

79
(a = 100 nm and ≈ 1 nm) (Hill 2006c, a).This difference affects all quantities derived from the microscale analysis that depend on convective flow3 .

STREAMING POTENTIAL
A fundamental characteristic of interest that has not been explicitly considered for hydrogel-colloid composites is the streaming potential.This is the differential potential ψ that prevails when n = i = 0 with u = 0 and p = 0. Accordingly, when φ 1, the streaming potential is where, recall, γ = D 1 /D 2 .This equation is obtained by multiplying Eqn. ( 23) for each species by z j e and setting their sum (electrical current) to zero.Note that terms involving the product φ ∇ψ are O(φ) smaller than the electromigrative term for the pure electrolyte, and a consistent approximation of Eqn. ( 24) is u ≈ −( 2 /η) p /L, since φ 1. Writing Eqn. ( 39) in terms of the dimensionless asymptotic coefficients provided in Table 3 above gives where ) with u * = o s (kT /e) 2 /(ηa).Note that n = I (bulk ionic strength).Equation ( 40) reveals that very large pressure differentials (∼ 10 5 Pa) are necessary to produce even small streaming potentials (∼ 1 mV).In general, therefore, the streaming potential is not a practically viable means of probing the microstructure.

CONCENTRATION-GRADIENT DRIVEN FLUX WITH ZERO ELECTRICAL CURRENT AND FLOW
In the following example, the flux is specified as j = −D e [ n (L) − n (0)] [Eqn (15)], with n (0) = 0.01 mmol l −1 , n (L) = 10 mmol l −1 and L = 500 μm, so ln n ≈ 6.9.Recall, this flux prevails in the absence charged inclusions when u = 0.The electrical current density i = 0, so j 1 = j 2 .Note that the numerically exact solutions do not yield n (x = L) = 10 mmol l −1 , because the charged inclusions change the macroscale fluxes.Also, while the calculations could be performed with p = 0 and, hence, u = 0, for simplicity the computations were undertaken with u = 0, so p = 0.The asymptotic coefficients are determined by the local bulk electrolyte concentration.Accordingly, from the local values of κa and ζ e/kT , asymptotic coefficients are obtained by interpolating refined versions of Tables I, II and III.Details of the microscale calculations, and a discussion of various quantities derived from them, are given by Hill (2006a, c).
When ln n 1, it is reasonable to approximate the asymptotic coefficients as constants based on an approximately constant bulk concentration n , so where, recall, γ = D 1 /D 2 , i = u = 0, and φ 1.Again, this equation is obtained by multiplying Eqn. ( 23) for each species by z j e and setting their sum (electrical current) to zero; note also that n /L ≈ ( n /L) ln n .Writing Eqn.(41) in terms of the dimensionless asymptotic coefficients provided in Tables I and II gives where with u * = o s (kT /e) 2 /(ηa).Note that n = I (bulk ionic strength), and that the increment is the same as the more general expression of Hill (2006b).Furthermore, a careful inspection of Tables I and II reveals that the dimensionless asymptotic coefficients generally yield φ E ∼ 1 only when κa < 1 and |ζ | > kT /e.Note that microscale theory does not account for particle interactions, so φ * = φ[1 + (κa) −1 ] 3 should be small.This is restrictive on the particle volume fraction φ = c(4/3)πa 3 when κa 1 (c is the particle number density).For the specific example introduced above with ln n ≈ 6.9, the macroscale calculations were undertaken with φ = 0.64(0, 0.1, 0.2, 0.4, 0.8)/[1 + (κa) −1 ] 3 , with κa < 1 evaluated at x = 0, which is generally the position where the ionic strength is lowest (with j > 0) and, hence, where κa is smallest.With n (0) = 0.01 mmol l −1 , φ ≈ (0, 0.0534, 0.1069, 0.2137, 0.4275) × 10 −3 .Note also that calculations were performed with a constant particle surface charge density σ = 1 μC cm −2 , with the ζ -potential varying according to the (semi-empirical) formula (Russel et al. 1989) When φ = 0, the bulk concentration varies approximately linearly across the membrane, as given by Eqn. ( 15) when Pe → 0. As the inclusion volume fraction increases, the specified flux j is achieved with practically the same (almost uniform) electrolyte concentration gradient, even though the membrane is macroscopically inhomogeneous due to the varying ζ -potential.As seen in Figure 1, the high surface charge density and low bulk electrolyte concentration produce a high ζ -potential at x = 0.With higher electrolyte concentrations, the ζ -potentials are low (with fixed surface charge) and the diffuse double layers thin; accordingly, the inclusions behave as impenetrable uncharged spheres and the membrane potential tends to its value for pure electrolyte, and the effective diffusion coefficients tend to their Maxwell values (see below).These and other limiting cases were thoroughly discussed by Hill (2006b, c).
As expected from the very small particle volume fractions (φ ∼ 10 −4 ), the charge on the inclusions has a negligible influence on transport that occurs predominantly by gradient diffusion alone.However, as shown in the top panel of Figure 2, the inclusions have a significant impact on the electrostatic potential.Increasing the particle concentration decreases the membrane potential ψ = ψ (x = L).At the highest particle concentration of only φ ≈ 4 × 10 −4 , ψ is about 20 mV lower than in the absence of inclusions.From Eqn. (20), ψ ≈ 24.6 mV.Therefore, the increment E ≡ [ ψ (x = L; φ = 0)/ ψ (x = L; φ) − 1]/φ ∼ −10 3 , which is in good agreement with expectations for uniform membranes with small macroscale gradients (Hill 2006c, Fig. 7).
As identified by Hill (2006c), the negatively charged inclusions here with Na + counterion reduce the effective asymmetry of the electrolyte by increasing (decreasing) the effective diffusion coefficient of the otherwise less (more) mobile Na + (Cl − ) ion.This diminishes the accompanying electric field required to maintain equal bulk fluxes (zero electrical current).The bottom panel of Figure 2 shows the electrostatic potential for membranes with the same particle volume fractions and electrolyte, but with constant particle ζ -potential (ζ = −4kT /e ≈ −100 mV) rather than constant surface charge.In both cases, the gradient of electrostatic potential is negative in a thin region at x = 0.In the absence of an accompanying pressure gradient, the large, positive electric field would drive electroosmotic flow in the direction of the bulk concentration gradient (from left to right).However, because these calculations have been performed with zero convective flow (u = 0), the pressure, which is plotted in Figure 3, varies in a very similar manner to the electrostatic potential.
A useful measure of the overall influence of the inclusions on transport is the bulk diffusion coefficient This is plotted in Figure 4 scaled with the diffusion coefficient D e given by Eqn.(17) for electrodiffusion in the absence of inclusions (φ = 0).The fact that D * ≈ 1 when x ≈ L reflects the underlying (almost) linear bulk ion concentration profile.However, near x = 0, where ζ e/(kT ) is large and κa is small, D * is a complicated and rapidly changing function of position and particle concentration.Note that the effective diffusion coefficient of a tracer in a dilute random array of impenetrable spherical inclusions is (Maxwell 1873, Batchelor andO'Brien 1977) 4 so with the largest value of φ ≈ 0.000428, Eqn.(47) gives D * /D ≈ 0.9994.However, the top panel of Figure 4 indicates that D * /D e ≈ 1.0030 when φ = 0.000428, so, overall, the charge on the inclusions produces a higher effective diffusivity than expected from Eqn. (47).47), but consistent with expectations for uniform composites with small macroscale gradients (Hill 2006b).For example, the incremental contributions to the individual ion diffusion coefficients in Table II ( B 1 for Na + and B 2 for Cl + ) (see also Hill 2006b, Fig. 6) reveal that the inclusions increase (decrease) the effective diffusion coefficient of the otherwise less (more) mobile Na + (Cl − ) ion.Accordingly, the negatively charged inclusions decrease the asymmetry of the NaCl electrolyte, which explains the overall decrease in the membrane potential shown in Figure 1.It is also evident from Table II

SUMMARY
Numerical solutions of the coupled macroscale ion conservation equations and fluid mass and momentum conservation equations for electrodiffusion across hydrogel-colloid composite membranes were presented for the first time.These computations were undertaken with the approximation of local macroscale electroneutrality, using bulk transport properties established in earlier work for microscale electrokinetic transport processes.The principal example presented was characterized by a low bulk electrolyte concentration on one side of the membrane, with a very low solid volume fraction of small, highly charged colloidal inclusions.The influence of these inclusions on the bulk fluxes was negligible in absolute terms, but the incremental contribution of the particles was significantly larger than expected from Maxwell's classical theory for impenetrable particles in a permeable continuum.Significant from a practical perspective was that very low concentrations of inclusions have a measurable impact on the membrane diffusion potential.Such changes reflect the size, charge, and concentration of the inclusions.Therefore, experimental measurements of the membrane potential may be promising for testing the microscale theory, which underlies several other interesting and technologically important microscale electrokinetic problems, including those of electrical microrheology and electroacoustics.

ACKNOWLEDGMENTS
RJH gratefully acknowledges support from the Natural Sciences and Engineering Research Council of Canada (NSERC) and the Canada Research Chairs program.
An Acad Bras Cienc (2010) 82 (1) in general, functions of the local scaled particle radius κa, scaled ζ -potential ζ e/(kT ), and scaled Brinkman screening length κ (or /a) for a given electrolyte [η is the fluid viscosity and 2 is the hydrodynamic (Darcy) permeability of the polymer skeleton].Again, with the local macroscale electroneutrality approximation, the bulk current density i = z 1 e j 1 + z 2 e j 2 (35) An Acad Bras Cienc (2010) 82 (1) ELECTRODIFFUSION IN HYDROGEL-COLLOIDAL COMPOSITES 75 is independent of ψ and p , the solution provides y ≡ y(x = L) − y(x = 0) in the parameter space comprising b and y(x = 0).In practical terms, this means the flux j = | j |, pressure differential p , and electrostatic potential differential ψ can be calculated (implicitly) as a function of the bulk convective flow u = | u |, and bulk ion concentrations on each side of the composite membrane: e.g., with n = n (x = L) − n (x = 0) and n (x = 0) as two independent scalars.

TABLE II Scaled (dimensionless) asymptotic coefficients for bulk diffusion of NaCl in a Brinkman medium with charged spherical inclusions (see Hill 2006b, for details
): /a = 0.1 (a = 10 nm, = 1 nm)

TABLE III Scaled (dimensionless) asymptotic coefficients for bulk convec- tion of NaCl in a Brinkman medium with charged spherical in- clusions (see Hill 2006c, for details)
that when ζ e/(kT ) is large, the inclusions enhance the effective diffusivity of Na + more than Cl − , i.e., | B 1 / B 2 | > 1.This explains the overall tendency for D * (L)/D e to be greater than one.Clearly, lower ζ -potentials will yield D * (L)/D e less than one, as expected from Eqn. (47) when ζ → 0.