Impact of Network Heterogenity on Nonlinear Electrokinetic Transport in Porous Media

TL;DR

This study uses numerical simulations to explore how heterogeneity in porous media affects nonlinear electrokinetic transport, revealing the dominance of electroconvection in certain conditions and highlighting the limitations of homogenized models.

Contribution

It introduces a detailed numerical framework to analyze the impact of pore network topology on electrokinetic phenomena, emphasizing the role of heterogeneity and flow loops.

Findings

Electroconvection dominates in large potential regimes for weak surface conduction.

Heterogeneous pore networks exhibit different transport behaviors compared to homogenized models.

The framework can inform design of porous media for water purification applications.

Abstract

We present a numerical study of nonlinear electrokinetic transport in porous media, focusing on the role of heterogeneity in a porous microstructure on ion concentration polarization and over-limiting current. For simplicity, the porous medium is modeled as a network of long, thin charged cylindrical pores, each governed by one-dimensional effective transport equations. For weak surface conduction, when sufficiently large potential is applied, we demonstrate that electrokinetic transport in a porous network can be dominated by electroconvection via internally induced flow loops, which is not properly captured by existing homogenized models. We systematically vary the topology and "accessivity" of the pore network and compare with simulations of traditional homogenized parallel-pore (capillary-bundle) models, in order to reveal the effects of regular and hierarchical connectivity. Our computational framework sheds light on the complex physics of nonlinear electrokinetic phenomena in microstructures and may be used to design porous media for applications, such as water desalination and purification by shock electrodialysis.

Figures

Tables

Table 1. Table 1 : Network dimensionless parameters

Dimensionless parameter	Description	Range
$p$	Pore blocking probability	[0, 0.5]
$C_v$	ratio of standard deviation to	[0, 0.667]
	mean of the pore size distribution
$C_s$	Network effective surface conduction	[0.004, 16]
${\lambda }_D$	Network effective Debye length	[0.002, 0.5]
$V/V_T$	Normalized applied potential	[0, 60]
$\alpha$	Accessivity	[0.02, 1]

Equations

$${\tilde{I}}_{\text{dl}}=I_{\text{lim}}(1-\exp (-V/V_T)),$$

$$I=I_{\text{dl}}+{\sigma }_{\text{OLC}}V/V_T,$$

$$S(x)\frac{\partial \overline{C^{-}}}{\partial t}+\frac{\partial }{\partial x}\{S(x)\overline{F^{-}}\}=0,$$

$$\overline{F^{-}}=\frac{\kappa }{2{\lambda }_D^2}\{\overline{C^{-}}\overline{g^p}+C_s\overline{g^{p-}}\}\frac{d{P}_0}{dx}+\{\kappa (\overline{C^{-}}\overline{g^e}+C_s\overline{g^{e-}})+\overline{C^{-}}\}\frac{d{\mu }^{+}}{dx}$$

$$+\{\frac{\kappa }{2{\lambda }_D^2}(\overline{C^{-}}\overline{g^c}+C_s\overline{g^{c-}})-2\bar{g}\}\frac{d{C}_0}{dx}.$$

$$h_p=\frac{\text{network total volume}}{\text{network total surface}}.$$

$$C_s=-\frac{\sigma }{ze{C}_{\text{ref}}{h}_p}=-2q_s^{\ast }{\lambda }_D^2,$$

$$I_{\text{dl}}={\tilde{I}}_{\text{dl}}/I_{\text{lim}}=\left(1-\exp (-\Delta {\mu }^{+}/2)\right).$$

$${\tilde{\sigma }}_{\text{olc,adv}}\sim q_s^{2/5}{C}_{\text{ref}}^{4/5}.$$

$$\frac{I}{I_{\text{lim}}}=1-{\exp }^{-\Delta V/V_T}+C_s\Delta V/V_T,$$

$${\sigma }_{\text{olc}}=\frac{dI/I_{\text{lim}}}{d\Delta V/V_T}=C_s.$$

Full text

\SectionNumbersOn\altaffiliation

Center for Turbulence Research, Stanford University, Stanford, California 94305, USA

\alsoaffiliationDepartment of Mechanical Engineering, Flow Physics and Computational Engineering, Stanford University, Stanford, California 94305, USA \altaffiliationCenter for Turbulence Research, Stanford University, Stanford, California 94305, USA

\phone+1650 7251325 \fax+1650 7253525

1 Introduction

Nonlinear electrokinetic phenomena, such as ion concentration polarization (ICP) 1, 2, electro-osmotic instability and electro-convection  3, 4, 5, 6, 7, 8, 9, 10, over-limiting current (OLC) 11, 12, 13, 14, 15, 16, deionization shock waves 17, 18, 19, 20, 21, current-induced membrane discharge 22, and induced-charge electro-osmosis 23, 24, have mostly been studied in simple geometries, defined by planar membranes or electrodes, isolated colloidal particles, or straight microfluidic channels. On the other hand, emerging applications, such as capacitive deionization 25 and ion separation 26, electro-osmotic pumping 27, resistive switching 28, electrokinetic fingering 29, shock electrodialysis 30, 31, 32, and shock electrodeposition 33, 34, naturally involve heterogeneous porous media, which leverage nanoscale electrokinetic phenomena at engineering scales more easily and effectively than is possible with microfluidic devices35, 36, 37, 38, 39. As such, there is a need to develop mathematical theories and simulation methods that go beyond classical “leaky membrane”  40, 17, 21, 41 and “porous electrode” 42, 43, 44 homogenized models to capture the complex physics of nonlinear electrokinetics in porous microstructures. In the absence of fluid flow, surface conduction (SC) has recently been shown to depend on the microstructure of porous electrodes 45, but to our knowledge, the coupled effects of SC and electrokinetic flows in complex porous networks have not yet been analyzed, despite their being implicated in experimental observations of OLC and deionization shock waves in porous media 31, 33, 34. Here, we leverage our recently developed computational framework for nonlinear electrokinetics in porous media 46, 47 to begin answering these fundamental questions related to OLC in charged porous media.

Diffusion-limited current is a general property of electrochemical systems 42, which separates distinct operating regimes in processes ranging from electrodialysis 11 to electrodeposition 48, 49. For steady conduction through a dilute, symmetric, binary bulk electrolyte from a reservoir to an ideally cation-selective membrane or electrode surface, the limiting current scales as $I_{\text{lim}}=\frac{2zeDC_{\text{ref}}}{L}A$, in which $C_{\text{ref}}$ is the reservoir concentration and $z,e$ and $D$ are ion valence, elementary charge, and ion diffusivity respectively. As discussed in the supplementary information, the dimensional diffusion-limited current ($\tilde{I}_{\text{dl}}$) in this prototypical case of ICP exhibits a diode-like dependence on the applied voltage ($V$) across electrolyte 1, 50, 12:

[TABLE]

where $V_{T}=\frac{k_{B}T}{\text{z}e}$ is the thermal voltage, $k_{B}$ is Boltzmann constant, and $T$ is absolute temperature. (Note that throughout this paper we distinguish dimensional variables from dimensionless values by using “~” on top of the symbols associated with dimensional variables.) The diffusion-limited current can only be exceeded (leading to OLC) if the preceding assumptions break down, by providing additional ions (via water splitting 11, 51 or membrane discharge 22), activating new transport mechanisms (such as electro-osmotic instability 9, 11, surface conduction or electro-osmotic flow 12, 13, 31, 16), or effectively shortening the diffusion length (by dendritic growth instability 49, 33, 34).

If electrolyte transport is confined laterally by charged surfaces in a microchannel or porous medium during the passage of OLC, then nonlinear dynamics such as deionization shocks 17, 18, 19, 20 may occur, which can be interpreted as shock waves in the concentration field. Deionization shocks can propagate in porous medium and trigger higher electric current 17, 31, 12, 40, 21 than $I_{\text{lim}}$. As first predicted by Dydek et al. 12, the current-voltage relation is nearly linear above the diffusion-limited regime across a range of experiments and theoretical models 31, 16, 40, 21, 15, 33, 34, 28:

[TABLE]

where $\sigma_{\text{OLC}}$ is normalized overlimiting conductance (scaled to $I_{\text{lim}}$), but the scaling of $\sigma_{\text{OLC}}$ with properties of the electrolyte and confining geometry has not yet been fully understood.

There have been a number of studies that aimed to understand the physical mechanisms that sustain OLC in electrokinetic systems. The basic mechanism of OLC sustained by electro-osmotic convection in a microchannels was first proposed by Yaroshchuk et al. 13, albeit using a Taylor-Aris model for hydrodynamic dispersion in the pressure-driven backflow that balances electro-osmotic forward flow, which only applies to very long microchannels. Dydek at al.  12 first predicted theoretically that OLC in a micro-channel can be dominated by one of three fundamental mechanisms – SC through electric double layers (EDLs), electro-osmosis (EO) along the sidewalls or electro-osmotic instability at the end surface, depending on the geometry and surface properties – and they proposed a modified scaling theory of OLC by EO convection. The predicted transition between the SC and EO mechanisms was later confirmed by Nam et al.  16 in microfluidic experiments, consistent with the predicted scalings of $\sigma_{\text{OLC}}$ with channel height in each regime. A few years earlier, Deng et al. 31 conducted the first experiments to investigate the contributing mechanisms to OLC in porous media. Their measurements revealed an approximate a power law dependence of $\sigma_{\text{OLC}}$ on reservoir concentration ($C_{\text{ref}}$) and pH-regulated surface charge $q_{s}$, which was consistent with theoretical prediction of Dydek et al. 12: $\sigma_{\text{OLC}}\sim q_{s}^{2/5}c_{0}^{4/5}$. This agreement was remarkable given that the theory of Dydek et al. was not derived for a porous material, but rather a single 2D micro-channel having only internal flow recirculation near the wall region. The general effect of surface convection leading to mixing during OLC in a microchannel was later studied by Rubinstein and Zaltzman 14, but again without considering any porous network effects. In order to rationalize the experimental results, Deng et al. 31 hypothesized that the scaling of Dydek must also hold at the network scale, based on eddies recirculating around loops in the pore network, rather than within individual pores.

While OLC and its mechanisms have been studied comprehensively in a micro-channel setting 12, 52, 13, 14, 15, 16 to the best of our knowledge there has not been any theoretical study of nonlinear electrokinetic phenomena in porous structures, going beyond macroscopic “leaky membrane” continuum models 17, 21, 40. More specifically, the impact of network heterogeneity and the induced flow on OLC and its connection with SC in porous structures have never been investigated through a rigorous numerical setting. In this article, we present numerical simulations of electrokinetic transport in porous structures that involve random networks of micro-scale and nano-scale pores, which appear to be the first computational investigations of nonlinear electrokinetic phenomena in porous networks.

We employ our previously developed computational framework 46, 47 to conduct the simulations of electrokinetic transport in massive random networks of pores. In particular, we examine the internally induced flows due to random pore coupling and connectivity loops and investigate how such flows can lead to enhancement of OLC in systems subject to ICP. As a measure of OLC, we report the I-V characteristics computed for a variety of porous networks and compare them to those obtained from corresponding homogenized models. As described in Section (2.3), homogenized models consist of decoupled straight pores and have the same total volume to surface area ratio (the same surface conduction) as the random porous networks. To develop comprehensive quantification we explore the impacts of network topology and pore size randomness, as well as electrolyte salinity on the overlimiting conductance. Lastly, we assess the applicability of the previously proposed theories 31, 12 based on the electroosmotic flow mechanism for OLC in porous materials.

2 Problem Setup

2.1 Computational Framework

In our previous two-paper series 46, 47, we presented comprehensive derivation and validation of our multi-scale model that enables pore-scale simulation of nonlinear electrokinetic phenomena in porous structures. We model a porous structure as a network of micro-scale and nano-scale pores, each is to be long and thin. We solve non-dimensionalized area-integrated transport equation for coin species:

[TABLE]

where $x$ coordinate is in pore’s axial direction, $S(x)$ is the local pore area cross-section, and $\overline{C^{-}}$ and $\overline{F^{-}}$ represent respectively the area-averaged concentration of coions and the area-averaged flux of coins. In our formulation, we write $\overline{F^{-}}$ in terms of local gradients of three driving potentials multiplied with pre-factors that capture the effects of local EDL and surface charge density:

[TABLE]

[TABLE]

We use the cross-sectional invariants of the three driving potentials, which are virtual total pressure, $P_{0}$, virtual counterion electro-chemical potential, $\mu^{+}$, and virtual ion concentration $C_{0}$. Detailed explanation of this choice is presented in our two-paper series 46, 47. $\bar{g}$, $\overline{g^{e}}$, $\overline{g^{c}}$, $\overline{g^{p-}}$, $\overline{g^{e-}}$, $\overline{g^{c-}}$ are the area-averaged coefficients that are tabulated with respect to the dimensionless EDL thickness, $\lambda^{*}$, and dimensionless surface charge density $q_{s}^{*}$. To close equation (3), we compute the fields of $P_{0}$ and $\mu^{+}$ by solving coupled conservation equations of fluid mass and electric charge. These equations involve additional area-averaged pre-factors which along with the aforementioned coefficients are tabulated with respect to $\lambda^{*}$ and $q_{s}^{*}$. Thus, when solving the discretized version of equation (3), we apply an interpolation procedure to extract the area-averaged coefficients from tabulated quantities based on the local values of $\lambda^{*}$ and $q_{s}^{*}$.

Equation (4) involves a number of dimensionless parameters. $\kappa=\frac{\varepsilon}{\mu D}V_{T}^{2}$ is electrohydrodynamic coupling constant. $\lambda_{D}(x)=\frac{1}{h_{p}}\sqrt{\frac{\varepsilon k_{B}T}{2\tilde{C}_{\text{ref}}\text{z}^{2}e^{2}}}$ is the ratio of reference Debye length to local effective pore size that is defined as the ratio of pore volume to pore surface area. $C_{s}=2|q_{s}^{*}|\lambda_{D}^{2}$ represents the excess couterions needed to neutralize surface charge at each local cross-section. $C_{s}$ indicates the significance of surface conduction effect.

2.2 Generation of Random Networks of Pores

Although our model does not have any constraint on the topology of a porous medium, we here focus our study on the networks that are constructed on a square background lattice, where pores are coupled horizontally or vertically to each other. We assume all pores have circular cross sections and all have the same length and surface charge density. To enforce the randomness of the network we consider three manipulations: first, we sample pores’ diameters from a log-normal distribution for a given mean and standard deviation. Second, we block pores randomly using a uniform random distribution. Pores with probability values below a certain quantity, $p$, are eliminated from the network. Third, we consider hierarchical structures by connecting series of lattices, in which each lattice has a scale twice finer than the previous lattice. In this strategy each lattice is constructed by the two previously mentioned manipulations. The final network is constructed by selecting the pore diameters based on the superposition of all scales.

In the present study, we consider two types of pore networks with varying degrees of parallel versus series connectivity. The first type has random pore sizes and random pore blocking with only a single-scale lattice, as shown in Figure (1(a)), and the second is constructed by all three manipulations mentioned above, as shown in Figure (1(b)). For the rest of this article, we call the former, “non-hierarchical network” and denote the latter by “hierarchical network”. For each network topology, we investigate how induced conduction due to network random coupling could contribute to OLC and the propagation of deionization shocks in overlimiting regime.

The two types of networks correspond to different ranges of “accessivity”, $\alpha$, a recently proposed parameter that quantifies the degree of parallel versus series connectivity in a finite porous network 53. The classical picture of a “capillary bundle” of parallel straight pores corresponds to the limit $\alpha=1$, since every pore is directly connected to the surface, like resistors in parallel. As we shall demonstrate, our hierarchical networks have intermediate accessivity, which is associated with greatly enhanced electrokinetic convection. In these networks internal pores are connected to the membrane by a chain of one or two pores in series, but there are also loops within the medium, which can support electrokinetic eddy formation. In contrast, our non-hierarchical regular networks have low accessivity, $\alpha\to 0$, since the number of accessible pores connected to the membrane surface scales with the (small) surface to volume ratio, measured in units of pore length, and pathways through the system involve many changes in pore diameter. The internal pores are thus mostly connected in like resistors in series, although again flow loops are possible, which are the focus of this study.

Although accessivity can be defined and measured experimentally by fitting capillary pressure versus saturation curves during immiscible fluid displacement (as in cyclic mercury porosimetry) to an invasion percolation theory, the measured accessivity can often be approximated by the fraction of pores that are directly connected to the external surface 53, or more generally to a well connected network of larger pores that transmit forces and flows to the representative volume of interest. The latter concept, measuring the fraction of micropores exposed to the external environment via well-connected macropores, has also been used recently to fit hysteretic moisture sorption isotherms (moisture content versus humidity during dry/wetting cycles) in a variety of hierarchical porous media, such as cement paste, porous glasses, carbon black, and dental enamel 54, 55, 56.

Motivated by these studies, we approximate the accessivity as follows for our different types of networks. For non-hierarchical porous networks, we define the effective accessivity as the ratio of the number of pores connected to the membrane surface (source of the largest electrokinetic forcing entering the network) to the total number of pores. This approximation procedure results in accessivities near zero for non-hierarchical lattices and equal to one for the homogenized capillary-bundle model. For hierarchical networks, due to the existence of multiple layers of pores at different scales, we utilize a slightly different method to estimate the relevant accessivity. Since transport in hierarchical networks is dominated by well-connected larger pores, we need to use a definition that captures network heterogeneity at that scale. To this end, we eliminate micro-scale pores, which are located primarily in the finest scale of the hierarchy. Due to higher hydrodynamic resistance, these pores have lower contribution to ion transport towards the membrane and are neglected, and the accessivity is computed from the network of larger pores using the same definition as we used for non-hierarchical networks. In other words, we estimate the accessivity of a hierarchical network by the accessivity of the macropores only, and below we show that this measure of network connectivity is indeed correlated with enhanced electro-osmotic convection.

2.3 Model Problem and Homogenized Model

Figure (2(a)) shows the schematic of the model problem consisting of a random network of pores confined by an impermeable cation-selective membrane. The membrane imposes zero net flux of anions ($\bar{F}^{-}(x=L_{x})=0$) and zero net flow ($\bar{u}(x=L_{x})=0$). The membrane is connected to a big reservoir that contains a symmetric electrolyte with reference concentration of $C_{\text{ref}}$. The electric potential is set to zero in this reservoir. Since the membrane and the reservoir are in equilibrium, the electrochemical potential of cations is nearly constant across the membrane-reservoir interface, as well as across the entire membrane. Thus, the electrochemical potential is fixed and equal to zero at the membrane interface with the random network ($\mu^{+}(x=L_{x})=0$). At x=0 there is a well-mixed reservoir with the same electrolyte concentration as the bottom reservoir. The reservoir pressure is set to atmospheric pressure ($P(x=0)=0$) and the dimensionless applied potential is $V/V_{T}$, which is equivalent to $\mu^{+}(x=0)=V/V_{T}$ at this boundary. Throughout this paper, we report the I-V characteristics with respect to the gradient of electrochemical potential of cations across the random network ( $\Delta\mu^{+}=V/V_{T}$).

We report an effective $C_{s}$ for a random network, which represents the significance of the surface conduction effect in ion transport and is defined based on network effective pore size, which can be computed as follows:

[TABLE]

The effective $C_{s}$ of a network is then defined as:

[TABLE]

where, $q_{s}^{*}$ and $\lambda_{D}$ are defined based on the network effective $h_{p}$. Dydek et al. 12 demonstrated that for a long and narrow pore the I-V characteristics in units of limiting current and thermal voltage indicates a constant overlimiting conductance equal to pore $C_{s}$. Later, we report our numerical simulations of electrokinetic transport in random networks of pores and investigate whether there are other mechanisms contributing to overlimiting conductance in these networks.

The network effective $h_{p}$ is used to construct an equivalent homogenized (capillary bundle) model of a random network of pores with 100% accessivity ($\alpha=1$). As shown in Figure (2(b)), the homogenized model is composed of parallel replicas of a single straight pore with the same $h_{p}$ as the network effective $h_{p}$. Given the results of Dydek et al. 12 the adapted method of constructing homogenized model ensures that its $C_{s}$ matches the one for the random network. Therefore, in the absence of mechanisms other than surface conduction we expect that the ratio of overlimiting conductance to limiting current in the homogenized model matches that in the random network. To match the limiting currents, one may manipulate either the number of single pores or their length (or their combination), but since this is a linear operation, it does not provide useful new information. In our computations we simulated a single pore of length, $L_{x}$, representing the homogenized model but compared I-V responses after normalizing current by the limiting value. Similar to the porous network, the homogenized model is blocked by a cation-exchange membrane at $x=L_{x}$. Thus, there is no net flow from the reservoir towards the selective membrane in both systems.

We compare the steady state behaviors of the two systems and investigate how pores’ random sizes and coupling in a network can influence the system performance against the corresponding homogenized model.

2.4 Dimensionless Parameters

We explore the dependency of overlimiting conductance on the pore blocking probability ($p$), the coefficient of pore size variability $C_{v}$ defined as the ratio of standard deviation of the pore size distribution to its mean, the effective $C_{s}$, and the effective $\lambda_{D}=\frac{\lambda_{\text{ref}}}{h_{p}}$. The electrolyte electrohydrodynamic coupling constant, $\kappa$, is set to 0.5 to limit the investigation to typical aqueous electrolytes. Table (1) lists the important dimensionless parameters and the ranges considered in this study.

The first two parameters, p and $C_{v}$, control network geometry, and as a part of this study we explored their influence on system response. We also examine a wide range of $C_{s}$ and $\lambda_{D}$ since the sensitivity to geometry itself may be dependent on these dynamically important parameters. Experimentally, these parameters may be varied by changing the surface material/chemistry and reservoir concentration. As a measure of system response, we consider steady state OLC over a wide range of applied voltages. As we will show, in the overlimiting regime the I-V response can be approximated by a constant slope, and thus we mostly report overlimiting conductance.

Since for each set of network parameters ($p$, $C_{v}$, $C_{s}$, $\lambda_{D}$, and $\kappa$) we construct the network using random numbers, our numerical results are prone to statistical variations. Only in the limit of very large networks with many pores, this statistical variation becomes negligibly small. In our investigation, we considered larger domain using many pores, as shown in Figure (1). Additionally, for each value of $p$ and $C_{v}$ we repeat our simulations for a number of random distributions and report the ensemble averaged results for I-V characteristics. In a separate study, to mimic the experimental investigation of Deng et al. 31 we considered variation of dimensional $C_{\text{ref}}$ to investigate the scaling of overlimiting conductance with respect to electrolyte salinity.

2.5 Computation of Current

We compute the current of a network at the membrane interface by summing up the current in the pores that are directly connected to the membrane. We report I-V curves normalized by $I_{\text{lim}}$ and $V_{T}$. $I_{\text{lim}}$ is the maximum amount of charge that can be driven by diffusion mechanism at zero surface charge and zero advection limit. This current can be computed by solving the diffusion equation in the network with zero surface charge with the boundary condition of $C=C_{\text{ref}}$ at the reservoir and $C=0$ at the membrane. While here we use a mathematical trick to compute $I_{\text{lim}}$ by solving diffusion equation with Dirichlet boundary conditions, in practice $I_{\text{lim}}$ can be realized in the limit of high voltages when surface charge is zero and $\lambda_{D}$ is very small. As shown in Supplementary Information, when the applied voltage is finite, $I_{\text{dl}}$ in equation (1) can be written in terms of electrochemical potential of counterions as follows:

[TABLE]

In the settings considered in this study, current above $I_{\text{lim}}$ is expected due to the effects of surface conduction and internally induced flows. Below, we introduce the idea of current decomposition that we will employ to quantify the impact of surface conduction as well as internally induced flows on the I-V curves. We break down the total current ($I_{\text{tot}}$) into diffusion limited current ($I_{\text{dl}}$) plus components that are driven by surface conduction through EDLs ($I_{\text{sc}}$) and advection effects ($I_{\text{adv}}$). One should note that the advection-driven overlimiting current is composed of two mechanisms, which are shown schematically in Figure (3). There is an overlimiting current driven by recirculations within pores (Figure(3(a))), which are generated when pressure driven flow and electroosmotic flow balance each other inside a single pore, leading to zero net flow rate. These recirculation zones can exist even in isolated single pores such as those in the homogenized models 31, 12. The second transport mechanism via advection involves network-level flow loops (Figure(3(a))), which are provoked by the competition between different modes of flow when pores of different sizes are connected in a loop.

We compute different components of $I_{\text{tot}}$ via the following calculations. First, $I_{\text{tot}}$ is computed by solving the full system of equations without any constraint. Second, the diffusion limitted current, $I_{\text{dl}}$, is computed via equation (7) and the computation of $I_{\text{lim}}$ from the solution to the diffusion equation, as discussed above. Thirdly, $I_{\text{dl}}+I_{\text{sc}}$ is computed by turning off all modes of advection. This can be done effectively by setting the tabulated coefficients associated with advection fluxes to zero, and therefore allowing only diffusion and surface conduction modes of transport. $I_{\text{adv}}$ can be computed by subtracting $I_{\text{dl}}+I_{\text{sc}}$ from $I_{\text{tot}}$. Similarly, we denote the overlimiting conductance due to surface conduction and advection by $\sigma_{\text{olc},\text{sc}}$ and $\sigma_{\text{olc},\text{adv}}$ respectively. These conductances can be computed respectively from the slopes of $I_{\text{sc}}$ and $I_{\text{adv}}$ versus the applied voltage. Following Dydek et al. 12, we define a conductance ratio $\Gamma=\frac{\sigma_{\text{olc},\text{adv}}}{\sigma_{\text{olc},\text{sc}}}$ that indicates the relative importance of advection effect to surface conduction. As we will discuss, for electrokinetic systems with high surface conduction and/or regular structures, $\Gamma<1$, whereas porous structures with low surface conduction and/or highly random topology can result in $\Gamma\geq 1$.

Before proceeding to the result discussion, we would like to point out that our primary objective of this study is to characterize porous materials and examine the impact of network parameters on the character of the I-V curve in the nonlinear regime. Thus, all the I-V curves in our study are normalized by network $I_{\text{lim}}$. We also document the impact of the geometry change on $I_{\text{lim}}$ itself in Supplementary Information. $I_{\text{lim}}$ indicates the ohmic conductance of the material (according to equation (7)) as $I_{\text{lim}}$ is proportional to $dI/dV$ in the limit of $V=0$. It is an important factor when it comes to the design and the selection of porous materials, and has been examined extensively in the context of percolation and the quantification of the impedance of random geometries subject to Laplace equation and Dirichlet boundary conditions 57, 58.Therefore, by documenting $I_{\text{lim}}$ we report solely the impact of geometric manipulation on the ohmic resistance of the network. As noted before, we report ensemble averaged results over 2 or 3 simulations for each set of network parameters. Our results indicate that there is a statistical variation of 1% to 8% in $I_{\text{lim}}$. For an accurate quantification more number of ensembles are required to further reduce the statistical error. Nevertheless, in the present analysis the error is sufficiently small to allow capture of the trends intended to be revealed by this study.

3 Results

3.1 Electrokinetic Transport in Non-hierarchical Structures

Figure (4(a)) depicts one example of a non-hierarchical network with low accessivity, $\alpha=0.03$. This network consists of 936 pores and 633 nodes, and the pore diameters are sampled from a log-normal distribution with mean diameter $d=600nm$ and standard deviation to mean ratio of $C_{v}=0.33$. The pore blocking probability is set to $p=0.2$. The random network is connected to a large reservoir that contain KCl solution with $C_{\text{ref}}=50mM$. We assume that all pores in the network have the same surface charge density equal to $q_{s}=-10mC/m^{2}$. This setting results in the network effective $h_{p}$ to be $167.06nm$, which corresponds to the network effective $\lambda_{D}=0.008$ and $C_{s}=0.012$. In this regime, each individual pore is thin enough to be in the SC dominated regime 12, so any significant effects of EO convection must involve loops in the porous network.

We used our computational model to study the electrokinetic transport in this network for a wide range of applied voltages. We here present the steady state response of the system for $V/V_{T}=50$. The concentration fields in the network and the corresponding homogenized model are presented in Figure (5). All pores in the homogenized model have the same $h_{p}$ equal to 167.06 nm and exhibit SC dominated OLC 12, 40, 15, 16. In the network the random coupling of pores with different pore diameters has caused the deionization shocks to propagate relatively heterogeneously away from the cation-selective membrane, however, the impact of network randomness does not seem to be prominent in this setting.

Figure (5(a)) depicts the induced flow field in the random network. Note that the macroscopic flow field is zero in both random network and the homogenized model. One can see that the internally induced flow exists only near the membrane region, where there is a dominant backward pressure-driven flow, which balances the induced local electroosmotic flows 13, 12, 31. The high electroosmotic flow in the near membrane region results from the higher electric field in this zone. As depicted in Figure (5), this zone has lower concentration due to ICP, and thus requires higher electric field to sustain the same current. The length scale of the flow loops (some loops are marked by red circles on Figure (5(a))) are at the scale of pore’s length or the length scale of the depletion zone. The induced flow loops have slightly helped the propagation of the deionization shock and made the shock front look nonuniform. Yet, the average impact on conductance in this example is not significant.

Figure(5(b)) depicts the I-V curves of the random network and the homogenized model, compared with the diffusion-limited curve or an uncharge material. All I-V curves are normalized by $I_{\text{lim}}$ and $V_{T}$. For high voltages the random network as well as the homogenized model exhibit overlimiting currents with constant conductances. The random network shows higher overlimiting current with slightly higher conductance, which results from the network random topology and the induced advection effect. The impact is not significant though, as a large portion of the domain involves zero flow and the induced flow loops are only active in a small region close to the membrane. The Péclet number associated with the flow loops generated in the network for $V/V_{T}=50$ is $Pe\leq 1$, which demonstrates that the advection effect is not significant in this configuration. To further investigate the contribution of advection to the overlimiting current of the system, we utilized the current decomposition method described in Section (2.5). The I-V characteristics for this setting (green curve) is close to the one obtained for the homogenized model. Thus, it confirms that the difference between the random network and the homogenized model is primarily due to the induced advection. However, the primary mode of overlimiting transport is surface conduction. We examined a wide variety of non-hierarchical networks using different combinations of $C_{v}$ and $p$ and the results are presented in the Supplementary Information. Our results demonstrate that despite a finite advection effect induced in all the non-hierarchical networks, surface conduction is the dominant mechanism of OLC with $\Gamma<1$ for all non-hierarchical networks of thin (locally SC dominated) charged pores.

3.2 Electrokinetic Transport in Hierarchical Structures

We here present one sample simulation of a hierarchical network composed of four scales having moderate accessivity, indicating a combination of both parallel and series network connections between pores of different sizes. Using the methodology described in Section (2.4), we obtained an intermediate accessivity is $\alpha\approx 0.114$ for this particular network.The mean size of the coarsest “macropores” is $2\mu m$, and finer “micropores” have proportionally smaller mean pore diameters at each level of recursion. The pore variability $C_{v}$ and pore blocking probability ($p$) are set equal to 0.3 and 0.2 respectively for all scales. The final network is depicted in Figure (6(a)). It consists of 557 pores and 340 nodes. The top reservoir contains KCl solution with concentration of $C_{\text{ref}}=50mM$ and the applied potential is $V/V_{T}=60$ at this boundary. The considered network has a uniformly distributed surface charge density equal to $q_{s}=-13mC/m^{2}$. This setting results in the effective $h_{p}=612.3nm$, which corresponds to the effective $\lambda_{D}=0.002$, and the effective $C_{s}=0.004$. Alongside the random network, we considered the equivalent homogenized model with the same $h_{p}$, and $C_{s}$. Similar to Section (3.1), we compare the steady state responses of the two systems including their concentration fields, flow fields and I-V characteristics.

Figure (6) indicates the steady state concentration fields for the random network and the homogenized model, when a potential gradient of $V/V_{T}=60$ is applied. As opposed to the homogenized model, where the concentration field is 1D, the random network involves a heterogeneous pattern for the propagation of deionization shocks leading to a 2D and complex concentration field. Moreover, the depletion zone is located close to the membrane interface in the homogenized model, while there is a dominant deionization shock in the middle of the random network, which has propagated almost all the way to the top boundary. As shown in Figure (7(a)), this shock has propagated in the direction of the pressure-driven flow away from the membrane and has completely depleted the bulk concentration in the middle vertical pores. Additionally, the rightmost pores in the network involve a dominant electroosmotic flow towards the membrane, which impedes the propagation of the deionization shock.

Figure (7(a)) and (7(b)) depict respectively the steady state flow field and the electric current field in the hierarchical network. As discussed before, although the macroscopic mean flow rate of the system is zero, there are internally induced flow loops resulting from the competition between electroosmotic flow and the backward pressure-driven flow. The internal flow loops would not occur if the pore sizes were uniform over the whole network, and having pore size variability is a key factor to provoke the competition of different modes of advection. As shown in Figure (7(a)), there are two major flow loops located on the right section of the network, which is the region containing pores with less hydraulic resistance than the bottom section of the medium. The strong backward pressure-driven flow in the middle vertical pore has been balanced by the dominant electroosmotic flow on the right side of the network. The length scale associated with the largest loop is equal to the length scale of the system in the x direction. We obtained the Pe$=uL/D=\frac{\Sigma_{i=1}^{n}A_{i}u_{i}}{\Sigma_{i=1}^{n}A_{i}}\frac{L}{D}$ (where n is the number of pores forming the loop) associated with this loop to be 8.1, which demonstrates that the induced advection is remarkably dominant over the molecular diffusion in this porous structure.

Considering both diagrams in Figure (7), one can see that the majority of the charge transport occurs in the rightmost vertical pores, where the fast flow is towards the membrane. The induced flow transports salt from the reservoir enhancing local electric conductance. Moreover, the current is almost zero in the middle vertical pathway as the backward pressure-driven flow has pushed the deionization shock towards the reservoir, leading to a substantial increase in the electric resistance. Since the system operates at low $C_{s}$, the charge transport due to surface conduction effect is small.

The I-V characteristics of the hierarchical network and the homogenized model were obtained over a wide range of applied potential, as shown in Figure (8). The homogenized model which is driven primarily by surface conduction has hardly moved above the limiting current, given that the system surface conduction is small. However, the hierarchical network has exhibited a significantly different response. For $V/V_{T}<20$, the two systems show similar I-V responses. For intermediate voltages ($5<V/V_{T}<20$), the induced flow field is not strong enough and the system is mainly controlled by surface conduction. However, after $V/V_{T}=20$ the overlimiting conductance increases considerably in the hierarchical random network compared to the homogenized model. This voltage can be regarded as a transition point from which the induced flow loops dominates the weak surface conduction effect and allows a substantial charge transport towards the membrane. Our results demonstrate that after the transition point, the Péclet numbers associated with the flow loops of the system become higher than 1.

In Figure (8), the green I-V curve denoted by $I_{dl}+I_{sc}$ belongs to the simulation of the hierarchical network where the advection terms are turned off. The I-V curve decomposition demonstrates that the surface conduction effect is well represented by the homogenized model, and it reconfirms the fact that all of the significant excess current in the full model must be due to the advection effect by the flow through the closed loops connecting multiple pores. This important advection mechanism is missing in the homogenized capillary bundle model.

We used our model to identify the dependency of overlimiting transport in hierarchical networks to $C_{v}$, $p$ and reported our investigation in Supplementary Information. Our results indicate that in all hierarchical networks there is a remarkable contribution of advection to overlimiting transport, and it is the dominant transport mechanism in the range of $C_{s}$ considered in this study. We consider hierarchical structures and aim at finding scaling laws relating the overlimiting conductance to system salinity, and examined the validity of theoretical scaling laws reported in the literature 31, 12 in high and low surface conduction regime.

3.3 Dependency on Bulk Salinity

Following the experiments of Deng et al. 31, we here present our investigation on the dependency of the I-V curve on system salinity in both low surface conduction and high surface conduction regime. Measurements of Deng et al. 31 matched the power law dependence of $\sigma_{olc}$ on $C_{\text{ref}}$ which was derived by Dydek et al. considering a single micro-channel12:

[TABLE]

Our analysis presented in Section (S2.1) in Supplementary Information, as well as the scaling estimates based on the theory of Dydek et al 12, suggests that a single-pore circulation mechanism cannot result in significant $\sigma_{\text{olc,adv}}$ in the regimes considered in the experiment. Given the results presented in earlier sections, only hierarchical structures may sustain significant $\sigma_{olc,adv}$. Therefore, in this section we investigate the scaling of the overlimiting conductance with respect to reservoir concentration ($C_{ref}$) for only hierarchical porous networks.

In our study we maintained $\kappa=0.5$, $p=0.2$, $C_{v}=33$, and $C_{s}$ fixed as we vary $C_{\text{ref}}$. This is equivalent to only changing $\lambda_{D}$. We intentionally chose to maintain fixed $C_{s}$ (instead of fixed $q_{s}$) since change of $C_{s}$ can transition the system from being dominated by surface conduction ($\sigma_{\text{olc}}\sim\sigma_{\text{olc,sc}}$) to that dominated by advection ($\sigma_{\text{olc}}\sim\sigma_{\text{olc,adv}}$). As we shall see, different scalings with respect to $C_{\text{ref}}$ are obtained in high and low surface conduction regimes. As constant $C_{s}$ implies $\sigma\sim C_{\text{ref}}$, the scaling proposed by Dydek et al. 12 in equation (8) suggests a 6/5 scaling of $\sigma_{\text{olc,adv}}$ with respect to $C_{\text{ref}}$.

Note that to eliminate the effect of statistical variations, we require to keep the network geometry fixed, as commonly done in experiments. As will be discussed, in high $C_{s}$ regime, the overlimiting conductance is primarily driven by surface effect ($\sigma_{\text{olc}}\simeq\sigma_{\text{olc, sc}}$). Conversely, in the limit of low surface conduction, the advection mechanism becomes dominant, and the charge transport is controlled primarily by the flow loops in the network ($\sigma_{\text{olc}}\simeq\sigma_{\text{olc, adv}}$).

3.3.1 High Surface Conduction Regime

We considered different $C_{s}$ values in the limit of high surface conduction regime, and computed the I-V characteristics over a range of $C_{ref}$. The hierarchical network introduced in Section (3.2) was considered for this study. To keep $C_{s}$ fixed, one needs to increase the surface charge of the network proportionally as $C_{ref}$ increases. Our results demonstrate that in the high surface conduction limit the contribution of surface conduction to charge transport always dominates advection mechanism, leading to $\Gamma<1$.

We have plotted the dimensional overlimiting conductance ($\tilde{\sigma}_{olc}$) versus $C_{\text{ref}}$ for three values of $C_{s}$ in Figure (9(a)). This result reveals that for high surface conduction regime, the overlimiting conductance is linearly proportional to $C_{ref}$. Similarly, we observed the linear relation between the overlimiting conductance and $C_{s}$, which is shown in Figure (9(b)). Using the homogenized model in the limit of thin double layers, one can use the area-averaged transport equation in high surface conduction regime to demonstrate the linear scaling theoretically as follows 12:

[TABLE]

where $\Delta V$ refers to the applied electric potential difference between the reservoir and the membrane. In the limit of high $\Delta V$, the overlimiting conductance is

[TABLE]

Note that although we compute $\sigma_{\text{olc}}$ with respect to cation electrochemical potential difference (we use $V/V_{T}=\Delta\mu^{+}$ instead of $\Delta V/V_{T}$ as the input), the relation in (10) still holds for our numerical results. This is because most of the electric potential gradient occurs in the depleted region, where $\bar{C}^{+}$ is almost constant, $\bar{C}^{+}\simeq C_{s}$, and the cation electrochemical potential difference is primarily due to electric potential gradient, and hence, $\Delta\mu^{+}\simeq\Delta V/V_{T}$.

Given $I_{\text{lim}}$ is the outcome of a solution to the linear Laplace equation, it must be proportional to $C_{\text{ref}}$ regardless of the network geometry. Therefore, from equation (10), it is expected that $\tilde{\sigma}_{\text{olc}}$ which is equal to $\sigma_{\text{olc}}\times I_{\text{lim}}/V_{T}$, grows proportional to $C_{\text{ref}}$ for a fixed $C_{s}$. Additionally, when $C_{\text{ref}}$ is fixed, this expression predicts the proportionality of $\sigma_{\text{olc}}$ with $C_{s}$, as observed in Figure (9(b)). These observations further confirm the fact that in the high $C_{s}$ limit advection plays negligible role in manipulating the I-V response.

3.3.2 Low Surface Conduction Regime

We here employed our model to simulate the electrokinetic transport in hierarchical networks for different values of $C_{\text{ref}}$ while fixing $C_{s}$ to 0.004. In Section (3.2) we showed that for a low value of $C_{s}$ advection could be the dominant mechanism for overlimiting transport in this configuration. We considered a wide range of $C_{\text{ref}}$ and investigated whether there is a power law relating the overlimiting conductance to $C_{\text{ref}}$, as explored in 31. To maintain a small $C_{s}=0.004$, the surface charge density of the system is proportionally decreased as we decreased $C_{\text{ref}}$.

We present our results for two hierarchical porous networks shown in Figure (10). Both porous structures are generated using the same input parameters, $C_{v}=0.33$, and $p=0.2$, $C_{s}=0.004$, and $\kappa=0.5$. For the given input parameters, we sampled different sets of random variables for blocking pores as well as the pore diameters. The effective $h_{p}$ in these structures are very close and equal to 612 nm.

When the reservoir concentration is $C_{\text{ref}}=50mM$, the internally induced flow fields at $V/V_{T}=50$ are shown in Figure (10) demonstrating that the strength and the length scales of the induced flow loops are different due to the randomness in the network geometries. We have compared the I-V response of the two porous systems obtained for $C_{\text{ref}}=50mM$ in Figure (11). The majority of the charge transport is driven by advection mechanism in both systems as the surface conduction effect is small. Yet, the advection effect is more dominant in Network 1 leading to the transport of higher amount of charge through this network.

Figure (12) depicts the dimensional overlimiting conductance versus $C_{\text{ref}}$ in the log-scale for two hierarchical networks. For low to moderate values of $C_{\text{ref}}$ we obtained the 1.4 scaling of $\tilde{\sigma}_{olc}$ with respect to $C_{\text{ref}}$ for Network 1, whereas the scaling associated with Network 2 turns out to be 1.2. Furthermore, the overlimiting conductance deviates to lower values for $C_{\text{ref}}$ larger than 50 mM, and the scaling does not hold.

Here we sought a power law scaling following the experimental work by Deng et al. 31. However, the plots shown in Figure (12) suggest that power laws are not the best fit to the observed trends. Furthermore, even the fitted power depends on the specific geometry, and thus is not universal. Interestingly, for one of the networks the fitted power law is equal to 1.2, which matches the theoretical prediction of Deng et al. 12, 31. In this work, however, demonstrate that a power law with exponent 1.2 (or any other value) cannot be universal for all types of porous networks. Since the glass frit of used by Deng et al. consists of fused and sintered silica particles, the porous network was likely not very hierarchical, so it may also be that additional mechanisms affected the results, such as charge regulation of the surfaces, local pH variations, etc.

Next, we provide deeper investigation on the trends of $\sigma_{\text{olc,adv}}$ obtained from our simulations. Before doing so, we provide some context leading to a plausible explanation for the observed behavior. By analogy to the scaling of $I_{\text{lim}}\sim zeDC_{\text{ref}}/L_{x}$ where $L_{x}$ is the length scale of the network, consider $I_{\text{olc}}\sim zeD_{\text{eff}}C_{\text{ref}}/L_{x}$ for overlimiting regime, where the effective diffusivity $D_{\text{eff}}$ incorporates the additional dispersion due to induced flow loops. In the overlimiting regime with dominant advection mechanism, $D_{\text{eff}}$ can be viewed as eddy diffusivity which scales as $U_{\text{eddy}}l_{\text{eddy}}$. As demonstrated in Section (3.2), in the overlimiting regime the length scale of the induced flow loop is comparable to the network length scale ($l_{\text{eddy}}\sim L_{x}$). Using maximum magnitude of velocity ($U_{\text{max}}$) for the scaling of $U_{\text{eddy}}$ yields $I_{olc}\sim zeD_{\text{eff}}C_{\text{ref}}L_{x}\sim U_{\text{max}}C_{\text{ref}}$. We investigated the relation between $U_{\text{max}}$ and $C_{\text{ref}}$ (Figure (13)) and observed the power fits with the exponents 0.4 and 0.2 respectively for Networks 1 and 2. This yields $I_{\text{olc}}\sim C_{\text{ref}}^{1.4}$ for Network 1 and $I_{\text{olc}}\sim C_{\text{ref}}^{1.2}$ for Network 2, which are close to the fitted results of $\sigma_{\text{olc}}$ obtained for low and moderate values of $C_{\text{ref}}$. This observation could imply that the flow loops with high Pe number are the key controlling factor of the network behavior in the low surface conduction regime.

Lastly, we provide insights on why $U_{\text{max}}$ may increase with the electrolyte salinity when $C_{s}$ is held fixed. Using our tabulated area-averaged coefficients 46, we re-plotted $\bar{g}_{e}$ coefficient, representing the electroosmotic velocity, with respect to $1/\lambda_{D}^{2}\sim C_{\text{ref}}$ for different values of $C_{s}$ (Figure (13(c))). This diagram illustrates that for a fixed $C_{s}$ increasing $C_{\text{ref}}$ results in faster electroosmotic flow in a single channel subject to a fixed electric potential. In the thin EDL limit, this trend is expected given that channel zeta potential ($\zeta$) is increased as both $C_{\text{ref}}$ and $q_{s}$ are proportionally increased. However, moving from low surface conduction to high surface conduction regime can result in various scaling of electroosmotic flow even in a single pore.

The largest slope at any point on the curves for low $C_{s}$ indicates a 1/2 power law for the dependency of $U_{\text{max}}$ on $C_{\text{ref}}$, which is higher than the power law fits obtained for the two networks in Figures (13(a)) and (13(b)). The marked points on Figure (13(c)) indicate the regime examined in the considered random networks. Coincidentally, these points correspond to the maximum scaling slope. Therefore, the lower power exponent, obtained in the data fit from the measurement of $U_{\text{max}}$ can be attributed to the pore coupling effects and partial balance of the electroosmotic velocity by pressure-driven flows.

Aside from complexities arising from pore coupling and specific network geometry, figure (13(c)) suggests that even in the simplest geometries there is no universal power law for the dependence of advective effects on $C_{\text{ref}}$. As $C_{\text{ref}}$ varies, the transition from thin to thick EDL and low to high zeta potential can significantly shift the trends in the induced electroosmotic flow. This complexity as well as aforementioned sensitivity to geometric details of a network, brings significant challenges to the task of characterization of transport in porous media.

4 Correlation with Accessivity

We investigated the correlation of accessivity and $\sigma_{OLC}$ for different network topologies considered in this study. For non-hierarchical networks by computing the ratio of number of pores connected to membrane surface to the total number of pores, we obtained low accessivity $\alpha\simeq 0.03$. For hierarchical networks as discussed in previous sections, we utilized a slightly different method to compute accessivity and obtained moderate accessivity $\alpha\simeq 0.12$ for hierarchical lattices constructed for different network parameters $C_{v}$ and $p$.

Porous structures that exhibit strong EO eddies around network loops tend to have intermediate values of accessivity, far from the limits of predominantly parallel ($\alpha=1$) and series ($\alpha=0$) connections between pores of different sizes. Networks with intermediate values of accessivity and intricate combinations of parallel and series couplings also tend to be isotropic, with connected paths of accessible pores in at least two directions, thus allowing for many loops in the network. In contrast, the homogenized capillary-bundle model has the maximum accessivity ($\alpha=1$) for surfaces perpendicular to current direction (such as the bounding membrane), but the minimum accessivity ($\alpha=0$) for restricted to surfaces in the transverse directions, thus eliminating the possibility of any flow loops. Our hierarchical networks that exhibit massive increases in OLC from EO convection loops are isotropic with moderate values of accessivity to bounding surfaces, via connected pathways of macropores, in both directions.

Figure(14) depicts the variation of $\sigma_{OLC}$ with respect to accessivity for various networks that we investigated in this article. Our data have formed two main clusters and are colored based on the type of network topology. The results associated with non-hierarchical networks are depicted by blue points and the results obtained for hierarchical networks are shown in green. All data points associated with non-hierarchical networks are placed in the left bottom corner of the graph, where represents low values of $\sigma_{OLC}$ and low values accessivity. These networks cannot create the macroscopic flow loops that are required to sustain larger $\sigma_{OLC}$.

Hierarchical networks have resulted in intermediate values of accessivity ($\alpha\simeq 0.1$), but their $\sigma_{OLC}$ covers relatively a wide range. This implies that to have large convective over-limiting conductance, it is necessary to have combinations of parallel and series couplings in loops (intermediate accessivity), which corresponds to moderate values of accessivity. However, this is not the sufficient condition for having a large overlimiting conductance, and other parameters such as pore size distribution are required along with accessivity to precisely characterize nonlinear elektrokinetic phenomena in porous media.

5 Conclusions

In summary, we presented a comprehensive numerical study of electrokinetic transport in hierarchical and non-hierarchical networks of pores with random pore coupling and pore blocking. Our simulations reveal that pore random size and pores’ coupling are necessary conditions to provoke internal flow loops that can considerably enhance charge transport through the porous system. When the Péclet number associated with these flow loops becomes larger than 1, the induced advection can lead to significant increase in the overlimiting transport and system performance. We demonstrated that this effect is completely missed by the homogenized model, leading to the mis-prediction of the physics in the porous structures. Our results suggest that the nonlinear electrokinetic characteristics of porous structures can be exploited for designing efficient and compact water deionization devices, separation techniques, and high performance micro-pumps.

Using our computational framework, we investigated the dependency of the overlimiting conductance on different parameters including network accessivity, pore size variation, pore blocking probability, and electrolyte salinity. We observed that for non-hierarchical networks of low accessivity, having mostly series connections from the membrane through pores of different sizes, the role of advection in overlimiting transport is not significant, even for fairly low $C_{s}$ value of $\sim 0.01$. On the contrary, in hierarchical structures of moderate and isotropic accessivity, having many different-sized internal pores connected in parallel to the transport surfaces, the network advection mechanism can have a tremendous impact on overlimiting ionic transport. We demonstrated that there is no universal scaling on the dependence of $\sigma_{\text{olc, adv}}$ on $C_{\text{ref}}$, at least in finite networks considered here, since having a dominant advection regime in the low surface connection limit inherently depends on the network geometry. While this makes characterization a challenging task, our work suggests that there are many opportunities for geometrical optimization and robust design of porous materials by predicting the physical trends relating the system performance to network’s geometrical parameters.

6 Acknowledgments

The authors gratefully acknowledge the fruitful comments by Zongyu (Joey) Gu on the definition and characterization of accessivity in porous media. This work is supported by the National Science Foundation (NSF) under award no. 1553275. {suppinfo} Derivation of Diffusion-limited Current. Non-hierarchical Networks: Dependency on Network Parameters. Hierarchical Networks: Dependency on Network Parameters.