# Direct numerical simulation of electrokinetic transport phenomena in fluids:
variational multi-scale stabilization and octree-based mesh refinement111Here,
by Direct numerical simulation, we mean that we resolve the physics of
interest at appropriate time and length scales

Sungu Kim  Kumar Saurabh  Makrand A. Khanwale  Ali Mani  Robbyn K. Anand
Baskar Ganapathysubramanian<EMAIL_ADDRESS>Department of Mechanical
Engineering, Iowa State University, Iowa, USA 50011  Department of Chemistry,
Iowa State University, Iowa, USA 50011  Center for Turbulence Research,
Stanford University, Stanford, California, USA 94305

###### Abstract

Computational modeling of charged species transport has enabled the analysis,
design, and optimization of a diverse array of electrochemical and
electrokinetic devices. These systems are represented by the Poisson-Nernst-
Planck (PNP) equations coupled with the Navier-Stokes (NS) equation. Direct
numerical simulation (DNS) to accurately capture the spatio-temporal variation
of ion concentration and current flux remains challenging due to the (a) small
critical dimension of the diffuse charge layer (DCL), (b) stiff coupling due
to fast charge relaxation times, large advective effects, and steep gradients
close to boundaries, and (c) complex geometries exhibited by electrochemical
devices.

In the current study, we address these challenges by presenting a direct
numerical simulation framework that incorporates (a) a variational multiscale
(VMS) treatment, (b) a block-iterative strategy in conjunction with semi-
implicit (for NS) and implicit (for PNP) time integrators, and (c) octree
based adaptive mesh refinement. The VMS formulation provides numerical
stabilization critical for capturing the electro-convective flows often
observed in engineered devices. The block-iterative strategy decouples the
difficulty of non-linear coupling between the NS and PNP equations and allows
the use of tailored numerical schemes separately for NS and PNP equations. The
carefully designed second-order, hybrid implicit methods circumvent the harsh
timestep requirements of explicit time steppers, thus enabling simulations
over longer time horizons. Finally, the octree-based meshing allows efficient
and targeted spatial resolution of the DCL. These features are incorporated
into a massively parallel computational framework, enabling the simulation of
realistic engineering electrochemical devices. The numerical framework is
illustrated using several challenging canonical examples.

###### keywords:

electrokinetics , electrohydrodynamics , Navier-Stokes Poisson Nernst-Planck ,
octrees , variational multiscale approach

††journal: Journal of Computational Physics

##  1 Introduction

A diverse set of energy conversion, energy storage, manufacturing, and
healthcare processes involve the transport and interaction of charged species
in aqueous environments. These interactions span multiple time and length
scales, often involving complex geometries [1, 2, 3, 4, 5, 6]. Examples of
such processes include electrolysis for hydrogen production, electrodialysis
for water desalination, gas diffusion electrodes for chemical conversion,
carbon capture, as well as several other processes that have a direct impact
on sustainability and climate resilience [7, 8, 9, 10, 11, 12].

Quantitative analyses of transport processes in devices are performed using
continuum mechanics-derived partial differential equations (PDEs),
specifically the Navier-Stokes (NS) governing the flow and the Poisson-Nernst-
Plank (PNP) equations governing chemical species transport and electrostatic
interactions. Detailed, high-fidelity simulations- i.e., direct numerical
simulations- are critical for comprehensive understanding, design, and
quantitative prediction of system-level responses of electrokinetic systems.
The coupling between the NS and PNP equation is responsible for producing a
variety of very interesting phenomena [7, 8, 9, 10, 11, 12, 13, 14, 15]. For
instance, even in systems with small hydrodynamic velocities, i.e., low
Reynolds (Re) numbers, a strong applied electric field triggers strongly non-
linear and multi-scale electrokinetic instabilities [16]. However, the
complexities arising from the tight coupling of multiscale and multiphysics
phenomena in these processes make simulating these systems very challenging. A
comprehensive simulation framework can significantly accelerate the
exploration and understanding of complex electrochemical processes and enable
the design and optimization of societally relevant applications like more
efficient water electrolysis and cheaper and lightweight dialysis machines.

One can broadly identify three challenges to performing direct numerical
simulations of NS-PNP equations in complex geometries: First, electrokinetic
systems involve coupled mass transport, fluid flow, and electrodynamics.
Additionally, the coupled non-linear NS-PNP system involves a variety of
stiffness. Therefore, such coupling and stiffnesses in the NS-PNP system call
for numerical methods that can robustly and efficiently the multiscale
physics. Second, such systems often exhibit thin regions with a high electric
field and high gradients in ion concentration in the diffused charge layer.
Accurately resolving these near-boundary gradients is essential, as they
determine not only the bulk behavior but also impact the boundary fluxes,
which are often the key quantities of interest from an experimental point of
view [17, 18, 19, 20]. The challenge of accurately resolving these regions is
further exacerbated by the complex geometries exhibited by various
electrochemical devices (packed beds, dendritic surfaces). This calls for
approaches that are endowed with the ability to spatially resolve complex
geometries in a scalable fashion. Third, simulating practical systems of
engineering interest requires performing simulations across a long time
horizon. This calls for time discretizations that are second-order (or higher)
and have good stabilizing properties, thus offering the ability to
successfully capture long time horizon simulations affordably [21].

Resolving these challenges requires careful design of spatial stabilization
and temporal discretization as well as appropriate strategies of (de)coupled
solution of the equations and adaptive mesh refinement. To resolve those
challenges, Druzgalski et al. [16] performed pioneering direct numerical
simulations of electroconvective instabilities. They used a non-uniform mesh
in the wall-normal direction to resolve the DCLs for small Debye lengths and a
semi-implicit time scheme. The framework was then pushed for more complex
surfaces with patterned electrochemical properties in [22, 23]. One of the
main challenges in simulating these systems is that if one resolves the thin
DCLs with a fine mesh resolution, one often requires very small timesteps for
stability reasons. The problem of resolving thin DCLs is exacerbated when the
problem length scale is larger, i.e., small non-dimensional Debye lengths,
requiring finer meshes and smaller timesteps. Druzgalski et al. [16] resolved
the challenge of small timesteps using implicit time-stepping. Specifically,
the framework in [16] utilized Fourier methods to efficiently handle linear
equations arising from stiffnesses in NS and PNP independently. However, such
strategies limit simulations to uniform meshes and periodic boundary
conditions in the transverse directions, preventing their extendability to
complex geometries. Additionally, structured meshes enforced fine resolution
in the wall-tangential direction, even in regions away from the walls, where
no physical small-scale structure exists.

There have been computational studies using commercial software to analyze the
non-linear behavior of electrokinetics [24, 25, 26]. It is important to note
that most commercial softwares are unsuitable for such electrochemical
systems, either because the numerical methods are not designed for these
phenomena or are not scalable enough to accurately resolve the spatial and
temporal scales required. Karatay et al. [27] showed that with a carefully
designed numerical method, their code was 20$\times$ faster than the
commercial alternative (COMSOL in this case) on a single processor. This
advantage amplified to about 160$\times$ speed up when both the codes were
deployed in parallel but on a single node. Therefore, careful design of
numerical methods and efficient parallel implementation are essential for
affordable high-fidelity simulations of electrochemical systems. While there
has been some progress in developing flexible frameworks [28, 29], numerical
challenges involved in solving the PNP equations still remain.

In this work, we seek to resolve all three challenges detailed above for the
direct numerical simulation of electrokinetic phenomena. We build upon our
previous work on a finite elements based framework for solving the NS and PNP
equations [30]. In contrast to [13], we deploy a novel framework that utilizes
a Variational multiscale method-based finite element approach and fast octree-
based meshing. Finite element approaches have been successfully applied for
high-resolution simulations of NS (see [31] for an overview) and PNP [32, 33,
34, 35, 36, 37, 38, 39, 40] equations. They are well suited for the NS-PNP
system of equations due to (a) the use of a variational formulation that
allows the weakening of derivatives, (b) the natural incorporation of
heterogeneous and mixed boundary conditions often seen in these systems, and
(c) ability to construct rigorous a posteriori error estimates for mesh
adaptivity that enable substantial computational gains.

The main contributions of this paper are as follows:

  1. 1.

We decouple the solution procedure of NS and PNP equations in a block-solve
strategy. The PNP equations are solved using a fully implicit Newton method-
based non-linear solver, ensuring the choice of large time steps even for
cases with fast charge relaxation times and very fine mesh resolution used for
thin DCLs. On the other hand, the Navier-Stokes equations are discretized
using a semi-implicit linear discretization, which also ensures stability for
larger time steps. In conjunction with this strategy, we utilize a second-
order time discretization for both NS and PNP to ensure accuracy across long
time horizons. Along with carefully chosen preconditioners, this block-
iterative strategy allows robust and efficient DNS computations.

  2. 2.

We use a variational multi-scale (VMS) based stabilization for the pressure
that allows for a pressure-coupled solver in conjunction with continuous
Galerkin finite elements. Additionally, the VMS method models fine-scale
velocities and ion concentrations not resolved by the grid size. This allows
us to use adaptive mesh discretization without compromising the accuracy of
the field variables.

  3. 3.

We use a massively parallel octree-based meshing framework to resolve the thin
DCLs for small Debye lengths and mesh complex geometries in 3D. We use Dendro-
KT [41], a highly scalable parallel octree library, to generate, manage and
traverse adaptive octree-based meshes in distributed memory. Octrees are
widely used in computational sciences [42, 43, 44, 45, 46, 47, 48, 49, 50],
due to its conceptual simplicity and ability to scale across a large number of
processors. We leverage the highly optimized and scalable implementations of
general finite element kernels in Dendrite-KT [51, 48, 52, 47, 53, 50, 54]
built on top of the scalable Dendro-KT framework. In this work, we use the
octree-based refinement to achieve highly localized refinement. We would like
to note that the framework also supports adaptive mesh refinement for
extensions to velocity/vorticity-based refinement or refinement based on
interfaces when the framework is extended to handle multiphase fluid flow.

  4. 4.

We illustrate the framework with a thorough series of test cases with
increasing complexity, including electro-convective instabilities and complex
geometries.

The rest of the paper is structured as follows: Section 2 details the
equations in non-dimensional form. Section 3 lays out the spatial and temporal
discretization and associated properties. Section 4 describes the
implementation and algorithmic details of the framework. Section 5 reports on
the numerical results, starting from convergence tests using manufactured
solutions to a full 3D simulation of a bio-diagnostic device. We conclude in
section 6.

##  2 Governing equations

###  2.1 Dimensional form:

Consider a bounded domain $\Omega\subset\mathbb{R}^{d}$ (for $d=2,3$) and the
time interval $[0,T]$. The domain contains a fluid (usually aqueous) with $N$
ionic species, and we seek to model the behavior of the $N$ species
concentration under the imposition of electric fields. We consider the case of
strong electrolytes with no reactions. The concentration distribution of the
$N$ species is represented using Nernst-Plank transport equations and the
Poisson equation for the electric potential, together making up the Poisson-
Nernst-Plank (PNP) equations. The pressure and fluid velocity are represented
using the Navier Stokes (NS) equations. There is a two-way coupling between
the species concentration and the underlying fluid motion, given by an
advection term in the PNP equation and an electrostatic body force term in the
NS equation.

We consider $N>1$ number of charged species with subscript $s$ indicating the
species index. The species flux, $j^{s,*}_{i}$, which is driven by diffusion,
migration, and convection is written as:

$\displaystyle j^{s,*}_{i}=-D^{s}\frac{\partial c^{s,*}}{\partial
x^{*}_{i}}-D^{s}\frac{z^{s}F}{RT}c^{s,*}\;\frac{\partial\phi^{s,*}}{\partial
x^{*}_{i}}+v^{*}_{i}c^{s,*}\quad\quad\text{\text{for $s^{th}$ species,}}$ (1)

$D^{s}$ is diffusivity of $s^{th}$ species, $c^{s,*}$ is concentration of
$s^{th}$ species, $z^{s}$ is valance of $s^{th}$ species, $F$ is Faraday
constant 96,485.33 $\mathrm{C}\mathrm{/}\mathrm{mol}$, $R$ is gas constant
8.314 $\mathrm{J}\mathrm{/}\mathrm{mol}\,\mathrm{K}$, $T$ is temperature,
$\phi^{*}$ is electric potential, and $v_{i}^{*}$ is fluid velocity.

######  Remark 1.

Note that we use Einstein notation throughout this work; in this notation,
$v_{i}$ represents the $i^{\text{th}}$ component of the vector $\mathbf{v}$,
and any repeated index is implicitly summed over. Superscript $s$ is the
species index. To avoid confusion powers of the terms are represented by
writing the power outside the braces. For example square of $s^{th}$ species
of concentration is written as $\left(c^{s,*}\right)^{2}$. For non-dimensional
quantities, $*$ in the superscript will be dropped.

The time variation of concentration can be obtained by taking the divergence
of species flux,

$\displaystyle\frac{\partial c^{s,*}}{\partial t^{*}}=-\frac{\partial
j^{s,*}_{i}}{\partial x^{*}_{i}},\quad\quad\text{\text{for the $s^{th}$
species.}}$ (2)

Combining eqs. 1 and 2 produces the Nernst-Planck equation,

$\frac{\partial c^{s,*}}{\partial
t^{*}}+\frac{\partial\left(v^{*}_{i}c^{s,*}\right)}{\partial
x^{*}_{i}}=\frac{\partial}{\partial x^{*}_{i}}\left(D^{s}\frac{\partial
c^{s,*}}{\partial
x^{*}_{i}}+D^{s}\frac{z^{s}F}{RT}c^{s,*}\;\frac{\partial\phi^{s,*}}{\partial
x^{*}_{i}}\right).$ (3)

The second term on the left-hand side of Eqn. 3 represents the advective
transport of species, which is the coupling between the species concentration
and the flow field. The first term on the right-hand side of Eqn. 3 represents
the diffusion of species, while the second term represents the electric
migration of species.

The electric potential, $\phi^{*}$ in the domain, depends on the concentration
distribution of the $N$ species and the imposed electric field. The potential
is the solution to the Poisson equation [55]

$\frac{\partial}{\partial
x^{*}_{i}}\left(-\varepsilon\frac{\partial\phi^{*}}{\partial
x^{*}_{i}}\right)=\rho_{e}^{*}.$ (4)

$\varepsilon$ is the electric permittivity of fluid, and $\rho_{e}$ is charge
density. Charge density can be expressed with species concentration.

$\rho_{e}^{*}\ =F\sum c^{s,*}z^{s}$ (5)

Then, equation (4) becomes (assuming that the permittivity, $\varepsilon$, is
constant)

$-\varepsilon\frac{\partial}{\partial
x^{*}_{i}}\left(\frac{\partial\phi^{*}}{\partial x^{*}_{i}}\right)=F\sum
c^{s,*}z^{s}.$ (6)

Finally, the Navier-Stokes equation is written as,

$\rho\frac{\partial v^{*}_{i}}{\partial
t^{*}}+\rho\frac{\partial\left(v^{*}_{j}v^{*}_{i}\right)}{\partial
x^{*}_{j}}+\frac{\partial p^{*}}{\partial
x^{*}_{i}}-\eta\frac{\partial}{\partial x^{*}_{j}}\left({\frac{\partial
v^{*}_{i}}{\partial x^{*}_{j}}}\right)-\rho_{e}^{*}E_{i}^{*}=0$ (7)

$\rho$ is the density of the fluid, $p$ is pressure, $\eta$ is dynamic
viscosity, and $E_{i}$ is $i^{th}$ component of the electric field. The last
term in the above equation is the body force due to the electric field
resulting in NS being coupled with the NS and the PNP. Using the relationship
between electric potential and electric field,

$-\frac{\partial\phi^{*}}{\partial x^{*}_{i}}=E^{*}_{i},$ (8)

and from the expression of charge density (equation (5)), equation (7) can be
written as

$\rho\frac{\partial v^{*}_{i}}{\partial
t^{*}}+\rho\frac{\partial\left(v^{*}_{j}v^{*}_{i}\right)}{\partial
x^{*}_{j}}+\frac{\partial p^{*}}{\partial
x^{*}_{i}}-\eta\frac{\partial}{\partial x^{*}_{j}}\left({\frac{\partial
v^{*}_{i}}{\partial x^{*}_{j}}}\right)+F\sum
c^{s,*}z^{s}\frac{\partial\phi^{*}}{\partial x^{*}_{i}}=0.$ (9)

######  Remark 2.

Note that some papers in the literature use elementary charge $e$ and
Boltzmann constant $k_{B}$ instead of Faraday constant $F$ and gas constant
$R$ [56]. The choice of these parameters is simply dependent on how to define
concentration. Here, we choose the Faraday constant, $F$, and gas constant $R$
since molarity $M$ was used for concentration.

We write all the governing equations together as,

$\displaystyle\begin{split}\text{Momentum Eqns:}&\quad\rho\frac{\partial
v^{*}_{i}}{\partial
t^{*}}+\rho\frac{\partial\left(v^{*}_{j}v^{*}_{i}\right)}{\partial
x^{*}_{j}}+\frac{\partial p^{*}}{\partial
x^{*}_{i}}-\eta\frac{\partial}{\partial x^{*}_{j}}\left({\frac{\partial
v^{*}_{i}}{\partial x^{*}_{j}}}\right)+F\left(\sum
c^{s,*}z^{s}\right)\frac{\partial\phi^{*}}{\partial x^{*}_{i}}=0,\end{split}$
(10) Solenoidality: $\displaystyle\quad\frac{\partial v^{*}_{i}}{\partial
x^{*}_{i}}=0,$ (11) Poisson:
$\displaystyle\quad-\varepsilon\frac{\partial}{\partial
x^{*}_{i}}\left(\frac{\partial\phi^{*}}{\partial x^{*}_{i}}\right)=F\sum
c^{s,*}z^{s},$ (12) Nernst-Planck: $\displaystyle\quad\frac{\partial
c^{s,*}}{\partial t^{*}}+v^{*}_{i}\frac{\partial c^{s,*}}{\partial
x^{*}_{i}}=\frac{\partial}{\partial x^{*}_{i}}\left(D^{s}\frac{\partial
c^{s,*}}{\partial
x^{*}_{i}}+D^{s}\frac{z^{s}F}{RT}c^{s,*}\;\frac{\partial\phi^{*}}{\partial
x^{*}_{i}}\right).$ (13)

Table 1 shows all the dimensional parameters involved.

Parameters | Units | Value/Formula
---|---|---
Species diffusivity, $D^{s}$ | $m^{2}/s$ |
Species valence, $z^{s}$ |  |
Faraday constant, $F$ | $C/mol$ | 96,485.33 $C/mol$
Gas constant, $R$ | $J/K\cdot mol$ | 8.314 $J/K\cdot mol$
Temperature, $T$ | $K$ |
Vacuum permittivity, $\varepsilon_{0}$ | $F/m,\>C/V\cdot m,\>N/V^{2}$ | 8.854$\times 10^{-12}F/m$
Relative permittivity, $\varepsilon_{r}$ |  |
Permittivity, $\varepsilon$ | $F/m,\>C/V\cdot m,\>N/V^{2}$ | $\varepsilon=\varepsilon_{r}\varepsilon_{0}$
Charge density, $\rho_{e}$ | $C$ |
Density, $\rho$ | $kg/m^{3}$ |
Viscosity, $\eta$ | $Pa\cdot s,\>cP$ | 8.90 $\times 10^{-4}P\cdot s$ for water at $25\text{\,}\mathrm{\SIUnitSymbolCelsius}$
Ionic strength of bulk, $I_{b}$ | $M,\>mol/m^{3}$ | $I_{b}=\frac{1}{2}\Sigma z_{i}^{2}c_{i}$
Thermal voltage, $V_{T}$ | $V$ | $V_{T}=RT/F$
Table 1: Dimensional parameters in the governing equations.

###  2.2 Non-dimensional form:

Each variable is scaled with its associated characteristic quantity to define
its non-dimensional counterpart.

$x_{i}=\frac{x_{i}^{*}}{L_{ch}},\quad t=\frac{t^{*}}{t_{ch}},\quad
v_{i}=\frac{v_{i}^{*}}{U_{ch}},\quad p=\frac{p^{*}}{p_{ch}},\quad
c_{i}=\frac{c_{i}^{*}}{c_{ch}},\quad\phi=\frac{\phi^{*}}{\phi_{ch}}.$ (14)

And for the derivative terms,

$\frac{\partial}{\partial x_{i}}=L_{ch}\frac{\partial}{\partial
x^{*}_{i}},\quad\frac{\partial}{\partial x_{i}}\left(\frac{\partial}{\partial
x_{i}}\right)=L_{ch}^{2}\frac{\partial}{\partial
x^{*}_{i}}\left(\frac{\partial}{\partial x^{*}_{i}}\right)$ (15)

To non-dimensionalize the governing equations, it is important to select
appropriate reference quantities in eq. 14. The reference concentration and
reference potential are selected as follows:

$c_{ch}=I_{b}=\frac{1}{2}\sum\left(z^{s}\right)^{2}c^{s}_{initial}\quad\text{and}\quad\phi_{ch}=V_{T}=\frac{RT}{F}$
(16)

$I_{b}$ is the ionic strength of bulk fluid at a reference time (initial), and
$V_{T}$ is the thermal voltage. The rest of the variables have multiple
options for the reference quantity. We choose the following reference
quantities: The characteristic velocity $U_{ch}$ is given by $D/L_{ch}$,
characteristic time $t_{ch}$ is given by $L_{ch}^{2}/D$ and characteristic
pressure given by $\eta D/L_{ch}^{2}$. Here, $D$ is the mean species
diffusivity ($D=1/N\sum^{N}_{s=1}D^{s}$), and $\eta$ is the dynamic viscosity
of the fluid medium. For the cases in this paper, we assume $D^{s}$ is the
same for all species; therefore, the non-dimensional $D^{s}$ equals 1.

With these characteristic quantities, the non-dimensional equations are given
as follows,

$\displaystyle\begin{split}\text{Momentum
Eqns:}&\quad\frac{1}{Sc}\left(\frac{\partial v_{i}}{\partial
t}+\frac{\partial\left(v_{j}v_{i}\right)}{\partial
x_{j}}\right)+\frac{\partial p}{\partial
x_{i}}-\frac{\partial^{2}v_{i}}{\partial
x_{j}^{2}}+\frac{\kappa}{2\Lambda^{2}}\left(\sum
c^{s}z^{s}\right)\frac{\partial\phi}{\partial x_{i}}=0,\end{split}$ (17)
Solenoidality: $\displaystyle\quad\frac{\partial v_{i}}{\partial x_{i}}=0,$
(18) Poisson: $\displaystyle\quad-2\Lambda^{2}\frac{\partial^{2}\phi}{\partial
x_{i}^{2}}=\sum c^{s}z^{s},$ (19) Nernst-Planck:
$\displaystyle\quad\frac{\partial c^{s}}{\partial t}+v_{i}\frac{\partial
c^{s}}{\partial x_{i}}=\frac{\partial}{\partial x_{i}}\left(\frac{\partial
c^{s}}{\partial x_{i}}+z^{s}c^{s}\;\frac{\partial\phi}{\partial
x_{i}}\right).$ (20)

Where the Schmidt number, $Sc=\frac{\eta}{\rho D}$, the electrohydrodynamic
(EHD) coupling constant, $\kappa=\frac{\varepsilon}{\eta
D}\left(\frac{RT}{F}\right)^{2}$, and non-dimensionalized Debye length,
$\Lambda=\lambda_{D}/L_{ch}$, where $\lambda_{D}$ is the Debye length,
$\lambda_{D}=\sqrt{\frac{1}{2}\frac{\varepsilon RT}{F^{2}I_{b}}}$.

##  3 Numerical method and its properties

We use a second-order backward difference formula (BDF) scheme as a time-
stepping strategy for eqs. 17 to 20. We solve momentum equations and Poisson
Nernst-Planck (PNP) equations in two separate blocks.

Let $\delta t$ be a time-step; let $t^{k}:=k\delta t$; let $n+1$ be the
current time-point, and $n$ be the previous time-point. We define an
intermediate velocity at time $n+1$ as follows:

$\widetilde{\mathbf{v}}^{k+1}:=2\mathbf{v}^{k}-\mathbf{v}^{k-1},$ (21)

We define the time-discretized weak form of the Navier-Stokes-Poisson-Nernst-
Plank (NS-PNP) equations are as follows:

Definition 1 Let $(\cdot,\cdot)$ be the standard $L^{2}$ inner product. We
state the time-discrete variational problem as follows: find
$\mathbf{v}^{k+1}(\mathbf{x})\in\mathbf{H}^{r}(\Omega)$,
$p^{k+1}(\mathbf{x})$, $\phi^{k+1}(\mathbf{x})$, $c^{s,k+1}(\mathbf{x})$
$\in{H}^{r}(\Omega)$ such that $\displaystyle\begin{split}\text{Momentum
Eqns:}&\quad\frac{1}{Sc}\left(w_{i},\,\frac{{\beta}_{0}v^{k+1}_{i}+{\beta}_{1}v^{k}_{i}+{\beta}_{2}v^{k-1}_{i}}{\delta
t}\right)_{\Omega}+\frac{1}{Sc}\left(w_{i},\,\,\widetilde{v}_{j}^{k+1}\,\frac{\partial
v_{i}^{k+1}}{\partial x_{j}}\right)_{\Omega}\\\
&\quad+\left(w_{i},\,\frac{\partial p^{k+1}}{\partial
x_{i}}\right)_{\Omega}+\left(\frac{\partial w_{i}}{\partial
x_{j}},\,{\frac{\partial v_{i}^{k+1}}{\partial
x_{j}}}\right)_{\Omega}+\frac{\kappa}{2\Lambda^{2}}\left(w_{i},\left(\sum
c^{s,k+1}z^{s}\right)\frac{\partial\phi^{k+1}}{\partial
x_{i}}\right)_{\Omega}\\\ &\quad-\left(w_{i},\,\hat{n}_{j}{\frac{\partial
v_{i}^{k+1}}{\partial x_{j}}}\right)_{\partial\Omega}=0,\end{split}$ (22)
Solenoidality: $\displaystyle\quad\quad\left(q,\,\frac{\partial
v^{k+1}_{i}}{\partial x_{i}}\right)_{\Omega}=0,$ (23) Poisson:
$\displaystyle\quad 2\Lambda^{2}\left(\frac{\partial q}{\partial
x_{i}},\frac{\partial\phi^{k+1}}{\partial
x_{i}}\right)_{\Omega}-2\Lambda^{2}\left(q,n_{i}\frac{\partial\phi^{k+1}}{\partial
x_{i}}\right)_{\partial\Omega}-\left(q,\sum c^{s,k+1}z^{s}\right)_{\Omega}=0,$
(24) $\displaystyle\begin{split}\text{Nernst-
Planck:}&\quad\left(q,\frac{{\beta}_{0}c^{s,k+1}+{\beta}_{1}c^{s,k}+{\beta}_{2}c^{s,k-1}}{\delta
t}\right)_{\Omega}+\left(q,\,\widetilde{v}_{i}^{k+1}\frac{\partial
c^{s,k+1}}{\partial x_{i}}\right)_{\Omega}\\\ &\quad+\left(\frac{\partial
q}{\partial x_{i}},\,\frac{\partial c^{s,k+1}}{\partial
x_{i}}+z^{s}c^{s,k+1}\;\frac{\partial\phi^{k+1}}{\partial
x_{i}}\right)_{\Omega}\\\ &\quad-\left(q,\,\hat{n}_{i}\left(\frac{\partial
c^{s,k+1}}{\partial x_{i}}+z^{s}c^{s,k+1}\;\frac{\partial\phi^{k+1}}{\partial
x_{i}}\right)\right)_{\partial\Omega}=0,\end{split}$ (25)
$\forall\mathbf{w}\in\mathbf{H}^{r}(\Omega)$, $\forall q\in H^{r}(\Omega)$,
given $\mathbf{v}^{k},\;\mathbf{v}^{k-1}\in\mathbf{H}^{r}(\Omega)$, and
$\phi^{k},\phi^{k-1},c^{s,k},c^{s,k}\in H^{r}(\Omega)$. $\hat{\mathbf{n}}$ is
the outward pointing normal to the boundary $\partial\Omega$

######  Remark 3.

Note that the variational forms above are for general Hilbert spaces with
regularity $r$ (order of allowed derivatives in the space). However, for the
cases in this paper, we restrict ourselves to $r=1$, which corresponds to
bilinear basis functions.

The semi-implicit discretization for eq. 22 is useful as it naturally
linearizes the momentum equations. This discretization is particularly useful
for chaotic systems where the timesteps are naturally small to capture small-
scale flow structures. Solving this linearized semi-implicit system helps with
performance as opposed to the fully implicit discretization, which requires a
Newton iteration and multiple sub-iterations (each iteration being a linear
solve).

###  3.1 Spatial discretization and the variational multiscale approach

We use continuous Galerkin finite elements for spatial discretization of the
time-discrete system presented in section 3. However, imposing the d’Alembert
constraint (eq. 23) when solving $\mathbf{v},\,p$ together presents a
challenge. An equal order polynomial approximation of velocity and pressure
leads to pressure instability because of the violation of the discrete inf-sup
condition (Ladyzhenskaya-Babuska-Brezzi condition, e.g., see Volker [31, page
31]). There are two methods to circumvent this issue: 1) add pressure
stabilization to eq. 23, which converts the inf-sup stability restriction to a
V-elliptic coercivity problem [57], thereby curbing the artificial pressure
instabilities; 2) The second method is to solve a separate pressure Poisson
equation using the Helmholtz Hodge decomposition of solenoidal velocity
fields. For our specific problem of NS-PNP equations, we use the first method.
This method allows us to use the Variational multi-scale method (VMS) [58] to
stabilize pressure. In addition to providing pressure stabilization, VMS also
has the following useful features. It allows us to naturally perform Large
Eddy simulation type decomposition and appropriate fine scale modeling [59].
The VMS enrichment for modeling fine-scale velocities, pressures, and ion
concentrations allows grid coarsening away from the extremely fine resolutions
near the boundary needed for the EDLs. Note that providing targeted resolution
near the boundary for resolving EDLs is relatively simpler; however, resolving
emerging small-scale structures due to electrochemical instabilities for high
electric fields is extremely difficult without doing adaptive meshing
222Adaptive mesh refinement to capture small-scale instabilities also presents
its challenges, primarily of constructing aposteriori error estimates for
solutions to Navier-Stokes which is an open research field. Therefore, the
subgrid-scale modeling away from the fine resolutions of the EDLs is valuable
for accurately capturing electroconvective instabilities with relatively
coarser resolutions. However, in this manuscript, we focus on much finer
resolutions, where the primary contribution of VMS is to circumvent the inf-
sup condition and provide pressure stabilization.

The VMS approach uses a direct-sum decomposition of the function spaces as
follows. If $\mathbf{v}\in\mathbf{V}$, $p\in Q$, and $\phi\in Q$ then we
decompose these spaces as:

$\displaystyle\mathbf{V}=\mathbf{V}^{c}\oplus\mathbf{V}^{f}\qquad\text{and}\qquad
Q=Q^{c}\oplus Q^{f},$ (26)

where $\overline{\mathbf{V}}$ and $\overline{Q}$ are the cG(r) subspaces of
$\mathbf{V}$ and $Q$, respectively, and the primed versions are the
complements of the cG(r) subspaces in $\mathbf{V}$ and $Q$, respectively. We
decompose the velocity and pressure as follows:

$\mathbf{v}=\mathbf{v}^{c}+\mathbf{v}^{f},\quad\phi=\phi^{c}+\phi^{f},\quad
c^{s}=c^{s,c}+c^{s,f},\quad\text{and}\quad p=p^{c}+p^{f},$ (27)

where the coarse scale components are $\mathbf{v}^{c}\in\mathbf{V}^{c}$ and
$p^{c},\phi^{c},c^{s,c}\in Q^{c}$, and the fine scale components are
$\mathbf{v}^{f}\in\mathbf{V}^{f}$ and $p^{f},\phi^{f}\in Q^{f}$. We define a
projection operator, $\mathscr{P}:\mathbf{V}\rightarrow\mathbf{V}^{c}$, such
that $\mathbf{v}^{c}=\mathscr{P}\\{\mathbf{v}\\}$ and
$\mathbf{v}^{f}=\mathbf{v}-\mathscr{P}\\{\mathbf{v}\\}$. Similar operators
decompose the other variables, $p$, $\phi$, $c^{s}$.

Substituting these decompositions in the original variational form in section
3 yields:

$\displaystyle\begin{split}\text{Momentum
Eqns:}&\quad\frac{1}{Sc}\left(w_{i},\;\frac{{\beta}_{0}v^{c,k+1}_{i}+{\beta}_{1}v^{c,k}_{i}+{\beta}_{2}v^{c,k-1}_{i}}{\delta
t}\right)_{\Omega}\\\
&+\frac{1}{Sc}\left(w_{i},\widetilde{v}_{j}^{k+1}\frac{\partial
v_{i}^{c,k+1}}{\partial
x_{j}}\right)_{\Omega}-\frac{1}{Sc}\left(\frac{\partial w_{i}}{\partial
x_{j}},\widetilde{v}_{j}^{k+1}v_{i}^{f,k+1}\right)_{\Omega}\\\
&+\left(w_{i},\frac{\partial\left(p^{c,k+1}+p^{f,k+1}\right)}{\partial
x_{i}}\right)_{\Omega}+\left(\frac{\partial w_{i}}{\partial
x_{j}},\frac{\partial v_{i}^{c,k+1}}{\partial
x_{j}}\right)_{\Omega}+\boxed{\left(\frac{\partial w_{i}}{\partial
x_{j}},\frac{\partial v_{i}^{f,k+1}}{\partial x_{j}}\right)_{\Omega}}\\\
&+\frac{\kappa}{2\Lambda^{2}}\left(w_{i},\left(\sum\left(c^{s,c,k+1}+c^{s,f,k+1}\right)z^{s}\right)\frac{\partial\phi^{c,k+1}}{\partial
x_{i}}\right)_{\Omega}\\\
&\boxed{+\frac{\kappa}{2\Lambda^{2}}\left(w_{i},\left(\sum\left(c^{s,c,k+1}+c^{s,f,k+1}\right)z^{s}\right)\frac{\partial\phi^{f,k+1}}{\partial
x_{i}}\right)_{\Omega}}=0,\end{split}$ (28)
$\displaystyle\begin{split}\text{Solenoidality:}&\quad\left(q,\frac{\partial
v_{i}^{c,k+1}}{\partial x_{i}}\right)_{\Omega}+\left(q,\frac{\partial
v_{i}^{f,k+1}}{\partial x_{i}}\right)_{\Omega}=0,\end{split}$ (29)
$\displaystyle\begin{split}\text{Poisson:}&\quad
2\Lambda^{2}\left(\frac{\partial q}{\partial
x_{i}},\frac{\partial\left(\phi^{c,k+1}+\phi^{f,k+1}\right)}{\partial
x_{i}}\right)_{\Omega}+2\Lambda^{2}\left(q,n_{i}\frac{\partial\left(\phi^{c,k+1}+\phi^{f,k+1}\right)}{\partial
x_{i}}\right)_{\partial\Omega}\\\ &-\left(q,\sum
c^{s,c,k+1}z^{s}\right)_{\Omega}-\left(q,\sum
c^{s,f,k+1}z^{s}\right)_{\Omega}=0,\end{split}$ (30)
$\displaystyle\begin{split}\text{Nernst-
Planck:}&\quad\left(q,\frac{{\beta}_{0}c^{s,c,k+1}+{\beta}_{1}c^{s,c,k}+{\beta}_{2}c^{s,c,k-1}}{\delta
t}\right)_{\Omega}+\left(q,\,\widetilde{v}_{i}^{k}\frac{\partial\left(c^{s,c,k+1}+c^{s,f,k+1}\right)}{\partial
x_{i}}\right)_{\Omega}\\\ &\quad+\left(\frac{\partial q}{\partial
x_{i}},\,\frac{\partial\left(c^{s,c,k+1}+c^{s,f,k+1}\right)}{\partial
x_{i}}+z^{s}\left(c^{s,c,k+1}+c^{s,f,k+1}\right)\;\frac{\partial\left(\phi^{c,k+1}+\phi^{f,k+1}\right)}{\partial
x_{i}}\right)_{\Omega}\\\
&\quad-\left(q,\,\hat{n}_{i}\left(\frac{\partial\left(c^{s,c,k+1}+c^{s,f,k+1}\right)}{\partial
x_{i}}+z^{s}\left(c^{s,c,k+1}+c^{s,f,k+1}\right)\;\frac{\partial\left(\phi^{c,k+1}+\phi^{f,k+1}\right)}{\partial
x_{i}}\right)\right)_{\partial\Omega}=0,\end{split}$ (31)

where $\mathbf{w},\mathbf{v}^{c},\in\mathscr{P}\mathbf{H}^{r}(\Omega)$,
$p^{c},\;\phi^{c}\in\mathscr{P}H^{r}(\Omega),\mathbf{v}^{f}\in(\mathscr{I}-\mathscr{P})\mathbf{H}^{r}(\Omega)$,
$\phi^{f},c^{s,f},p^{f}\in(\mathscr{I}-\mathscr{P})H^{r}(\Omega)$, and
$c^{s,c},q\in\mathscr{P}H^{r}(\Omega)$. Here $\mathscr{I}$ is the identity
operator, and $\mathscr{P}$ is the projection operator.

We can further simplify eq. 31,

$\displaystyle\begin{split}&\left(q,\frac{{\beta}_{0}c^{s,c,k+1}+{\beta}_{1}c^{s,c,k}+{\beta}_{2}c^{s,c,k-1}}{\delta
t}\right)_{\Omega}+\left(q,\,\widetilde{v}_{i}^{k}\frac{\partial\left(c^{s,c,k+1}+c^{s,f,k+1}\right)}{\partial
x_{i}}\right)_{\Omega}\\\ &\quad+\left(\frac{\partial q}{\partial
x_{i}},\,\frac{\partial\left(c^{s,c,k+1}+c^{s,f,k+1}\right)}{\partial
x_{i}}+z^{s}\left(c^{s,c,k+1}+c^{s,f,k+1}\right)\;\frac{\partial\left(\phi^{c,k+1}+\phi^{f,k+1}\right)}{\partial
x_{i}}\right)_{\Omega}\\\
&\quad-\left(q,\,\hat{n}_{i}\left(\frac{\partial\left(c^{s,c,k+1}+c^{s,f,k+1}\right)}{\partial
x_{i}}+z^{s}\left(c^{s,c,k+1}+c^{s,f,k+1}\right)\;\frac{\partial\left(\phi^{c,k+1}+\phi^{f,k+1}\right)}{\partial
x_{i}}\right)\right)_{\partial\Omega}=0,\end{split}$ (32)
$\displaystyle\begin{split}\implies&\quad\left(q,\frac{{\beta}_{0}c^{s,c,k+1}+{\beta}_{1}c^{s,c,k}+{\beta}_{2}c^{s,c,k-1}}{\delta
t}\right)_{\Omega}+\left(q,\,\widetilde{v}_{i}^{k}\frac{\partial
c^{s,c}}{\partial
x_{i}}\right)_{\Omega}+\left(q,\,\widetilde{v}_{i}^{k}\frac{\partial
c^{s,f}}{\partial x_{i}}\right)_{\Omega}\\\ &\quad+\left(\frac{\partial
q}{\partial x_{i}},\,\frac{\partial c^{s,c,k+1}}{\partial
x_{i}}\right)_{\Omega}+\boxed{\left(\frac{\partial q}{\partial
x_{i}},\,\frac{\partial c^{s,f,k+1}}{\partial x_{i}}\right)_{\Omega}}\\\
&\quad+\left(\frac{\partial q}{\partial
x_{i}},\,z^{s}c^{s,c,k+1}\;\frac{\partial\phi^{c,k+1}}{\partial
x_{i}}\right)_{\Omega}+\left(\frac{\partial q}{\partial
x_{i}},\,z^{s}c^{s,f,k+1}\;\frac{\partial\phi^{c,k+1}}{\partial
x_{i}}\right)_{\Omega}\\\ &\quad+\boxed{\left(\frac{\partial q}{\partial
x_{i}},\,z^{s}c^{s,c,k+1}\;\frac{\partial\phi^{f,k+1}}{\partial
x_{i}}\right)_{\Omega}}+\boxed{\left(\frac{\partial q}{\partial
x_{i}},\,z^{s}c^{s,f,k+1}\;\frac{\partial\phi^{f,k+1}}{\partial
x_{i}}\right)_{\Omega}}\\\ &\quad-\left(q,\,n_{i}\frac{\partial
c^{s,c,k+1}}{\partial
x_{i}}\right)_{\partial\Omega}-\boxed{\left(q,\,n_{i}\frac{\partial
c^{s,f,k+1}}{\partial x_{i}}\right)_{\partial\Omega}}\\\
&\quad-\left(q,\,n_{i}z^{s}c^{s,c,k+1}\;\frac{\partial\phi^{c,k+1}}{\partial
x_{i}}\right)_{\partial\Omega}-\boxed{\left(q,\,n_{i}z^{s}c^{s,f,k+1}\;\frac{\partial\phi^{c,k+1}}{\partial
x_{i}}\right)_{\partial\Omega}}\\\
&\quad-\boxed{\left(q,\,n_{i}z^{s}c^{s,c,k+1}\;\frac{\partial\phi^{f,k+1}}{\partial
x_{i}}\right)_{\partial\Omega}}-\boxed{\left(q,\,n_{i}z^{s}c^{s,f,k+1}\;\frac{\partial\phi^{f,k+1}}{\partial
x_{i}}\right)_{\partial\Omega}}\\\ &\quad=0,\end{split}$ (33)

We set the boxed terms to zero (see next paragraph, and also, Table. 2). These
assumptions reduce the Nernst-Plank equations to the following,

$\displaystyle\begin{split}\text{Nernst-
Planck:}&\quad\left(q,\frac{{\beta}_{0}c^{s,c,k+1}+{\beta}_{1}c^{s,c,k}+{\beta}_{2}c^{s,c,k-1}}{\delta
t}\right)_{\Omega}+\left(q,\,\widetilde{v}_{i}^{k}\frac{\partial
c^{s,c,k+1}}{\partial
x_{i}}\right)_{\Omega}+\left(q,\,\widetilde{v}_{i}^{k}\frac{\partial
c^{s,f,k+1}}{\partial x_{i}}\right)_{\Omega}\\\ &\quad+\left(\frac{\partial
q}{\partial x_{i}},\,\frac{\partial c^{s,c,k+1}}{\partial
x_{i}}\right)_{\Omega}+\left(\frac{\partial q}{\partial
x_{i}},\,z^{s}c^{s,c,k+1}\;\frac{\partial\phi^{c,k+1}}{\partial
x_{i}}\right)_{\Omega}+\left(\frac{\partial q}{\partial
x_{i}},\,z^{s}c^{s,f,k+1}\;\frac{\partial\phi^{c,k+1}}{\partial
x_{i}}\right)_{\Omega}\\\ &\quad-\left(q,\,n_{i}\frac{\partial
c^{s,c,k+1}}{\partial
x_{i}}\right)_{\partial\Omega}-\left(q,\,n_{i}z^{s}c^{s,c,k+1}\;\frac{\partial\phi^{c,k+1}}{\partial
x_{i}}\right)_{\partial\Omega}=0\end{split}$ (34)
$\displaystyle\begin{split}\implies&\quad\left(q,\frac{{\beta}_{0}c^{s,c,k+1}+{\beta}_{1}c^{s,c,k}+{\beta}_{2}c^{s,c,k-1}}{\delta
t}\right)_{\Omega}+\left(q,\,\widetilde{v}_{i}^{k}\frac{\partial
c^{s,c,k+1}}{\partial x_{i}}\right)_{\Omega}-\left(\frac{\partial q}{\partial
x_{i}},\,\widetilde{v}_{i}^{k}c^{s,f,k+1}\right)_{\Omega}\\\
&\quad+\left(\frac{\partial q}{\partial x_{i}},\,\frac{\partial
c^{s,c,k+1}}{\partial x_{i}}\right)_{\Omega}+\left(\frac{\partial q}{\partial
x_{i}},\,z^{s}c^{s,c,k+1}\;\frac{\partial\phi^{c,k+1}}{\partial
x_{i}}\right)_{\Omega}+\left(\frac{\partial q}{\partial
x_{i}},\,z^{s}c^{s,f,k+1}\;\frac{\partial\phi^{c,k+1}}{\partial
x_{i}}\right)_{\Omega}\\\
&\quad+\boxed{\left(q,\,n_{i}\widetilde{v}_{i}^{k}c^{s,f,k+1}\right)_{\partial\Omega}}-\left(q,\,n_{i}\frac{\partial
c^{s,c,k+1}}{\partial
x_{i}}\right)_{\partial\Omega}-\left(q,\,n_{i}z^{s}c^{s,c,k+1}\;\frac{\partial\phi^{c,k+1}}{\partial
x_{i}}\right)_{\partial\Omega}=0.\end{split}$ (35)

Table 2 lists the various boxed terms in the coupled set of NS-PNP equations
that are set to zero, along with the rationale for setting them to zero.

Equation |  Terms |  Rationale for assumption
---|---|---
Momentum equations |  $\left(\frac{\partial w_{i}}{\partial x_{j}},\frac{\partial v_{i}^{f,k+1}}{\partial x_{j}}\right)_{\Omega}$ |  Projection: Orthogonality conditions from the projection decomposition, as $w_{i}\in\mathscr{P}\mathbf{H}^{r}(\Omega)$ and $v_{i}^{f,k+1}\in(\mathscr{I}-\mathscr{P})\mathbf{H}^{r}(\Omega)$ [60, 61, 58]
Momentum equations |  $\frac{\kappa}{2\Lambda^{2}}\left(w_{i},\left(\sum\left(c^{s,c,k+1}+c^{s,f,k+1}\right)z^{s}\right)\frac{\partial\phi^{f,k+1}}{\partial x_{i}}\right)_{\Omega}$ |  $\phi^{f,k+1}$ assumed to be zero
Nernst-Planck equations |  $\left(q,\,n_{i}\widetilde{v}_{i}^{k}c^{s,f,k+1}\right)_{\partial\Omega}$ |  Depends on the boundary condition, could be zero by either a Dirichlet condition of velocity or a Dirichlet condition of concentration
Nernst-Planck equations |  $\left(\frac{\partial q}{\partial x_{i}},\,\frac{\partial c^{s,f,k+1}}{\partial x_{i}}\right)_{\Omega}$ |  Projection: Orthogonality conditions from the projection decomposition, as $q\in\mathscr{P}H^{r}(\Omega)$ and $c^{s,f,k+1}\in(\mathscr{I}-\mathscr{P})H^{r}(\Omega)$ [60, 61, 58]
Nernst-Planck equations |  $\left(\frac{\partial q}{\partial x_{i}},\,z^{s}c^{s,c,k+1}\;\frac{\partial\phi^{f,k+1}}{\partial x_{i}}\right)_{\Omega}$ |  $\phi^{f,k+1}$ assumed to be zero
Nernst-Planck equations |  $\left(\frac{\partial q}{\partial x_{i}},\,z^{s}c^{s,f,k+1}\;\frac{\partial\phi^{f,k+1}}{\partial x_{i}}\right)_{\Omega}$ |  $\phi^{f,k+1}$ assumed to be zero
Nernst-Planck equations |  $\left(q,\,n_{i}\frac{\partial c^{s,f,k+1}}{\partial x_{i}}\right)_{\partial\Omega}$ |  flux of $c^{s,f,k+1}$ assumed to be zero on boundaries
Nernst-Planck equations |  $\left(q,\,n_{i}z^{s}c^{s,f,k+1}\;\frac{\partial\phi^{c,k+1}}{\partial x_{i}}\right)_{\partial\Omega}$ |  $c^{s,f,k+1}$ assumed to be zero at the boundary
Nernst-Planck equations |  $\left(q,\,n_{i}z^{s}c^{s,c,k+1}\;\frac{\partial\phi^{f,k+1}}{\partial x_{i}}\right)_{\partial\Omega}$ |  $\phi^{f,k+1}$ assumed to be zero
Nernst-Planck equations |  $\left(q,\,n_{i}z^{s}c^{s,f,k+1}\;\frac{\partial\phi^{f,k+1}}{\partial x_{i}}\right)_{\partial\Omega}$ |  $c^{s,f,k+1}$ assumed to be zero at the boundary
Poisson equations |  $2\Lambda^{2}\left(\frac{\partial q}{\partial x_{i}},\frac{\partial\phi^{f,k+1}}{\partial x_{i}}\right)_{\Omega}$ |  $\phi^{f,k+1}$ assumed to be zero
Poisson equations |  $2\Lambda^{2}\left(q,n_{i}\frac{\partial\phi^{f,k+1}}{\partial x_{i}}\right)_{\partial\Omega}$ |  $\phi^{f,k+1}$ assumed to be zero
Table 2: Fine scale terms assumed to be zero and along with the rationale.

With all the simplifications listed in Table 2 we can now write the fully
discrete variational form for the NS-PNP equations.

Definition 2 Let $(\cdot,\cdot)$ be the standard $L^{2}$ inner product. We
state the time-discrete variational problem as follows: find
$\mathbf{v}^{c,k+1},\in\mathscr{P}\mathbf{H}^{r}(\Omega)$, and
$p^{c,k+1},\;c^{s,c,k+1}(\mathbf{x}),\phi^{c,k+1}\in\mathscr{P}H^{r}(\Omega)$
such that $\displaystyle\begin{split}\text{Momentum
Eqns:}&\quad\frac{1}{Sc}\left(w_{i},\;\frac{{\beta}_{0}v^{c,k+1}_{i}+{\beta}_{1}v^{c,k}_{i}+{\beta}_{2}v^{c,k-1}_{i}}{\delta
t}\right)_{\Omega}\\\
&+\frac{1}{Sc}\left(w_{i},\widetilde{v}_{j}^{k+1}\frac{\partial
v_{i}^{c,k+1}}{\partial
x_{j}}\right)_{\Omega}-\frac{1}{Sc}\left(\frac{\partial w_{i}}{\partial
x_{j}},\widetilde{v}_{j}^{k+1}v_{i}^{f,k+1}\right)_{\Omega}\\\
&+\left(w_{i},\frac{\partial\left(p^{c,k+1}+p^{f,k+1}\right)}{\partial
x_{i}}\right)_{\Omega}+\left(\frac{\partial w_{i}}{\partial
x_{j}},\frac{\partial v_{i}^{c,k+1}}{\partial x_{j}}\right)_{\Omega}\\\
&+\frac{\kappa}{2\Lambda^{2}}\left(w_{i},\left(\sum
c^{s,k+1}z^{s}\right)\frac{\partial\phi^{c,k+1}}{\partial
x_{i}}\right)_{\Omega}=0,\end{split}$ (36) Solenoidality:
$\displaystyle\quad\quad\left(q,\frac{\partial v_{i}^{c,k+1}}{\partial
x_{i}}\right)_{\Omega}+\left(\frac{\partial q}{\partial
x_{i}},v_{i}^{f,k+1}\right)_{\Omega}=0,$ (37)
$\displaystyle\begin{split}\text{Poisson:}&\quad
2\Lambda^{2}\left(\frac{\partial q}{\partial
x_{i}},\frac{\partial\phi^{c,k+1}}{\partial
x_{i}}\right)_{\Omega}-2\Lambda^{2}\left(q,n_{i}\frac{\partial\phi^{c,k+1}}{\partial
x_{i}}\right)_{\partial\Omega}\\\ &\quad-\left(q,\sum
c^{s,c,k+1}z^{s}\right)_{\Omega}-\left(q,\sum
c^{s,f,k+1}z^{s}\right)_{\Omega}=0,\end{split}$ (38)
$\displaystyle\begin{split}\text{Nernst-
Planck:}&\quad\left(q,\frac{{\beta}_{0}c^{s,c,k+1}+{\beta}_{1}c^{s,c,k}+{\beta}_{2}c^{s,c,k-1}}{\delta
t}\right)_{\Omega}\\\ &\quad+\left(q,\,\widetilde{v}_{i}^{k}\frac{\partial
c^{s,c,k+1}}{\partial x_{i}}\right)_{\Omega}-\left(\frac{\partial q}{\partial
x_{i}},\,\widetilde{v}_{i}^{k}c^{s,f,k+1}\right)_{\Omega}\\\
&\quad+\left(\frac{\partial q}{\partial
x_{i}},\,z^{s}c^{s,c,k+1}\;\frac{\partial\phi^{c,k+1}}{\partial
x_{i}}\right)_{\Omega}+\left(\frac{\partial q}{\partial
x_{i}},\,z^{s}c^{s,f,k+1}\;\frac{\partial\phi^{c,k+1}}{\partial
x_{i}}\right)_{\Omega}\\\
&\quad-\left(q,\,n_{i}z^{s}c^{s,c,k+1}\;\frac{\partial\phi^{c,k+1}}{\partial
x_{i}}\right)_{\partial\Omega}\\\ &\quad+\left(\frac{\partial q}{\partial
x_{i}},\,\frac{\partial c^{s,c,k+1}}{\partial
x_{i}}\right)_{\Omega}-\left(q,\,n_{i}\frac{\partial c^{s,c,k+1}}{\partial
x_{i}}\right)_{\partial\Omega}=0.\end{split}$ (39)
$\forall\mathbf{w}\in\mathscr{P}\mathbf{H}^{r}(\Omega)$,
$q\in\mathscr{P}H^{r}(\Omega)$ , and
$\mathbf{v}^{f,k+1}\in(\mathscr{I}-\mathscr{P})\mathbf{H}^{r}(\Omega)$,
$c^{s,f,k+1},p^{f,k+1}\in(\mathscr{I}-\mathscr{P})H^{r}(\Omega)$ , given
$\mathbf{v}^{k},\;\mathbf{v}^{k-1}\in\mathbf{H}^{r}(\Omega)$, and
$\phi^{k},\phi^{k-1},c^{s,k},c^{s,k}\in H^{r}(\Omega)$.

Note that the above variational problem is not closed, as the fine-scale
velocities, pressure, and ion concentrations are unknown. To close the system,
we use the residual-based approximation proposed in Bazilevs et al. [60] for
the fine-scale components. Defining $\mathcal{R}_{m}$, $\mathcal{R}_{sol}$,
and $\mathcal{R}_{c^{s}}$ as the residuals of the momentum, solenoidality and
Nernst-Plank equations, the fine-scale components are:

$v_{i}^{f}=-\tau_{m}\mathcal{R}_{m}\left(v_{i}^{c},p^{c},c^{s,c}\right),\qquad
p^{f}=-\tau_{sol}\mathcal{R}_{sol}(v_{i}^{c}),\qquad\text{and}\qquad
c^{s,f}=-\tau_{c^{s}}\mathcal{R}_{c^{s}}\left(c^{s,c},v_{i}^{c},\phi^{c}\right),$
(40)

where,

$\displaystyle\begin{split}\tau_{m}&=\left(\frac{4}{\Delta
t^{2}}+v_{i}^{c}G_{ij}v_{j}^{c}+C_{I}\left(\frac{1}{Re}\right)^{2}G_{ij}G_{ij}\right)^{-1/2},\end{split}$
(41)
$\displaystyle\begin{split}\tau_{sol}&=\frac{1}{tr(G_{ij})\tau_{m}},\end{split}$
(42) $\displaystyle\tau_{c^{s}}$ $\displaystyle=\left(\frac{4}{\Delta
t^{2}}+v_{i}^{c}G_{ij}v_{j}^{c}+C_{I}G_{ij}G_{ij}\right)^{-1/2}.$ (43)

Here, we set $C_{I}$ and $C_{\phi}$ for all our simulations to 6, and
$G_{ij}=\sum_{k=1}^{3}\partial\xi_{k}/\partial x_{i}\partial\xi_{k}/\partial
x_{j}$ is a mesh-based tensor that accounts for the inverse mapping
($\partial\xi_{k}/\partial x_{i}$) between the parametric and the physical
domain of the element. For octrees with an equal aspect ratio ($\Delta
x=\Delta y=\Delta z=h$), such a transformation simplifies $G_{ij}$ to
$(2/h)^{2}$, where $h$ is the element’s size. Subsequently, the final fully-
discrete variational form that is solved can be written as follows.

Definition 3 Let $(\cdot,\cdot)$ be the standard $L^{2}$ inner product. We
state the time-discrete variational problem as follows: find
$\mathbf{v}^{c,k+1},\in\mathscr{P}\mathbf{H}^{r}(\Omega)$, and
$p^{c,k+1},\;c^{s,c,k+1}(\mathbf{x}),\phi^{c,k+1}\in\mathscr{P}H^{r}(\Omega)$
such that $\displaystyle\begin{split}\text{Momentum
Eqns:}&\quad\frac{1}{Sc}\left(w_{i},\;\frac{{\beta}_{0}v^{c,k+1}_{i}+{\beta}_{1}v^{c,k}_{i}+{\beta}_{2}v^{c,k-1}_{i}}{\delta
t}\right)_{\Omega}\\\
&+\frac{1}{Sc}\left(w_{i},\widetilde{v}_{j}^{k+1}\frac{\partial
v_{i}^{c,k+1}}{\partial
x_{j}}\right)_{\Omega}+\frac{1}{Sc}\left(\frac{\partial w_{i}}{\partial
x_{j}},\widetilde{v}_{j}^{k+1}\tau_{m}\mathcal{R}_{m}\left(v_{i}^{c,k+1},p^{c,k+1},c^{s,c,k+1}\right)\right)_{\Omega}\\\
&+\left(w_{i},\frac{\partial p^{c,k+1}}{\partial
x_{i}}\right)_{\Omega}+\left(\frac{\partial w_{i}}{\partial
x_{i}},\tau_{sol}\mathcal{R}_{sol}\left(v_{i}^{c,k+1}\right)\right)_{\Omega}\\\
&+\left(\frac{\partial w_{i}}{\partial x_{j}},\frac{\partial
v_{i}^{c,k+1}}{\partial x_{j}}\right)_{\Omega}\\\
&+\frac{\kappa}{2\Lambda^{2}}\left(w_{i},\left(\sum
c^{s,k+1}z^{s}\right)\frac{\partial\phi^{c,k+1}}{\partial
x_{i}}\right)_{\Omega}=0,\end{split}$ (44) Solenoidality:
$\displaystyle\quad\quad\left(q,\frac{\partial v_{i}^{c,k+1}}{\partial
x_{i}}\right)_{\Omega}-\left(\frac{\partial q}{\partial
x_{i}}q,\tau_{m}\mathcal{R}_{m}\left(v_{i}^{c,k+1},p^{c,k+1},c^{s,c,k+1}\right)\right)_{\Omega}=0,$
(45) $\displaystyle\begin{split}\text{Poisson:}&\quad
2\Lambda^{2}\left(\frac{\partial q}{\partial
x_{i}},\frac{\partial\phi^{c,k+1}}{\partial
x_{i}}\right)_{\Omega}-2\Lambda^{2}\left(q,n_{i}\frac{\partial\phi^{c,k+1}}{\partial
x_{i}}\right)_{\partial\Omega}\\\ &\quad-\left(q,\sum
z^{s}c^{s,c,k+1}\right)_{\Omega}+\left(q,\sum
z^{s}\tau_{c^{s}}\mathcal{R}_{c^{s}}\left(v_{i}^{c,k+1},\phi^{c,k+1}\right)\right)_{\Omega}=0,\end{split}$
(46) $\displaystyle\begin{split}\text{Nernst-
Planck:}&\quad\left(q,\frac{{\beta}_{0}c^{s,c,k+1}+{\beta}_{1}c^{s,c,k}+{\beta}_{2}c^{s,c,k-1}}{\delta
t}\right)_{\Omega}\\\ &\quad+\left(q,\,\widetilde{v}_{i}^{k}\frac{\partial
c^{s,c,k+1}}{\partial x_{i}}\right)_{\Omega}+\left(\frac{\partial q}{\partial
x_{i}},\,\widetilde{v}_{i}^{k}\tau_{c^{s}}\mathcal{R}_{c^{s}}\left(v_{i}^{c,k+1},\phi^{c,k+1}\right)\right)_{\Omega}\\\
&\quad+\left(\frac{\partial q}{\partial
x_{i}},\,z^{s}c^{s,c,k+1}\;\frac{\partial\phi^{c,k+1}}{\partial
x_{i}}\right)_{\Omega}-\left(\frac{\partial q}{\partial
x_{i}},\,z^{s}\tau_{c^{s}}\mathcal{R}_{c^{s}}\left(v_{i}^{c,k+1},\phi^{c,k+1}\right)\;\frac{\partial\phi^{c,k+1}}{\partial
x_{i}}\right)_{\Omega}\\\
&\quad-\left(q,\,n_{i}z^{s}c^{s,c,k+1}\;\frac{\partial\phi^{c,k+1}}{\partial
x_{i}}\right)_{\partial\Omega}\\\ &\quad+\left(\frac{\partial q}{\partial
x_{i}},\,\frac{\partial c^{s,c,k+1}}{\partial
x_{i}}\right)_{\Omega}-\left(q,\,n_{i}\frac{\partial c^{s,c,k+1}}{\partial
x_{i}}\right)_{\partial\Omega}=0.\end{split}$ (47)
$\forall\mathbf{w}\in\mathscr{P}\mathbf{H}^{r}(\Omega)$,
$q\in\mathscr{P}H^{r}(\Omega)$ , and
$\mathbf{v}^{f,k+1}\in(\mathscr{I}-\mathscr{P})\mathbf{H}^{r}(\Omega)$,
$c^{s,f,k+1},p^{f,k+1}\in(\mathscr{I}-\mathscr{P})H^{r}(\Omega)$ , given
$\mathbf{v}^{k},\;\mathbf{v}^{k-1}\in\mathbf{H}^{r}(\Omega)$, and
$\phi^{k},\phi^{k-1},c^{s,k},c^{s,k-1}\in H^{r}(\Omega)$.

##  4 Solution strategy and design of the numerical framework

###  4.1 Solution strategy

It is important to note that because we are using a block solution method, the
fully discrete Navier-Stokes equations, eqs. 44 to 45, and the Poisson Nernst
Planck equations, eqs. 46 to 47, are solved as two different sub-problems.
Note that the NS (eqs. 44 to 45) is a linear system, whereas, PNP (eqs. 46 to
47) is a non-linear system. We use the Newton method to reduce the non-linear
system (eqs. 44 to 45) to an iteration of linear problems. Symbolically, we
can write this nonlinear algebraic system as

$T_{i}\left(U^{k}_{1},U^{k}_{2},\dots,U^{k}_{n}\right)=0,$ (48)

where $\mathbf{U}^{k}$ is a finite-dimensional vector that contains all of the
degrees of freedom at the discrete-time level $t^{k}$. In the system of eqs.
44 and 45, $\mathbf{U}^{k}$ contains ion concentrations and electric potential
at time level $k$. Then, the Newton iteration can be posed as,

$\displaystyle A_{ij}^{m,k}\;\delta
U_{j}^{m,k}=-T_{i}\left(U_{1}^{m,k},U_{2}^{m,k},\dots,U_{n}^{m,k}\right),\quad\quad
A_{ij}^{m,k}:=\frac{\partial}{\partial
U_{j}}F_{i}\left(U_{1}^{m,k},U_{2}^{m,k},\dots,U_{n}^{m,k}\right),$ (49)

where $U_{j}^{m,k}$ is the vector containing all the degrees of freedom at the
$k^{\text{th}}$ time step and at the $m^{\text{th}}$ Newton iteration. $\delta
U_{j}^{m,k}$ is the “variation” vector that will be used to update the current
Newton iteration guess:

$\displaystyle U_{j}^{m+1,k}=U_{j}^{m,k}+\delta U_{j}^{m,k}.$ (50)

and $A_{ij}^{m,k}$ is the Jacobian matrix computed analytically by calculating
the variations (partial derivatives) of the operators with respect to the
degrees of freedom. The iterative procedure begins with an initial guess,
which is set to the solution from the previous time step:

$U_{i}^{0,k}=U_{i}^{k-1},$ (51)

and ends once we reach the desired tolerance:

$\|\delta U_{j}^{m,k}\|\leq\text{TOL}.$ (52)

The block-iteration strategy for solving the system is given in Fig. 1. The
flow chart in fig. 1 shows that the linear semi-implicit Navier Stokes system
is solved first from the available data from the previous timestep.
Subsequently, the Poisson Nernst-Planck system is solved from the data from
the previous timestep and the first block of Navier-Stokes. The block
iteration is repeated until a desired convergence is reached.

Note that the momentum block is linear due to the semi-implicit
discretization. The number of iterations of the non-linear Poisson-Nernst
Planck (PNP) block depends on the timestep. We generally observe that the non-
linear PNP block converges in two to three Newton iterations (each iteration
being a linear solve).

We choose the timestep for accuracy reasons close to the advective CFL of
momentum equations. The implicit treatment of higher order operators and the
coupling operator allows us to do this (Note that this timestep is still
larger compared to the case if the operators were treated explicitly). Then,
both PNP and momentum blocks are generally iterated twice. In some cases, like
the chaotic electroconvective instability (section 5.2), where the timestep
($1\text{\times}{10}^{-4}$) is chosen to be of the order of smallest element
size (0.0009765625). The timestep is small enough that the block iteration
converges in just one block of momentum and PNP. The accuracy of quantities of
interest (QOI) from our solutions is compared in Fig. 7, and it shows very
good agreement with high-resolution simulations from Druzgalski et al. [16].
We have added this explanation to the manuscript.

For chaotic cases where a small timestep is naturally chosen to capture the
flow structures, the PNP system relaxes much faster than the momentum system.
Therefore, when a small timestep is chosen for momentum (of the order of the
smallest mesh size), the convergence of PNP is even faster, and the block
iteration converges in one iteration. Numerical experiments reveal that two
block iterations produce tight convergence for all examples, with one block
iteration often providing good results for most examples.

Figure 1: Flowchart for the block iteration technique as described in section
4.1 .

###  4.2 Software design

The software is designed in a layer-by-layer fashion. The numerical simulation
presented in this work is implemented on top of massively parallel framework
Dendrite-KT. Dendrite-KT uses Dendro-KT for adaptive octree mesh generation
and Petsc for linear algebra solver by utilizing KSP (Krylov Subspace) context
for solving a linear system of equations and SNES context for Newton
iteration. Dendro-KT allows allows for targeted, dynamic mesh refinement based
on either the location or aposteriori estimates.

The organization of the overall framework is illustrated in Fig. 2. The yellow
layers represent the computational framework. The modeling framework - shown
as three blocks on top of the computational framework, represents the fully
discrete scheme presented in section 3. The non-linear solve for the PNP
system (eqs. 46 to 47) is solved by coupling the data structures on the octree
meshing framework with Petsc using the SNES context. Subsequently, the inner
linear solves required for the PNP system and the NS system (eqs. 44 to 45)
are solved using the KSP construct. This allows us to exploit the plethora of
preconditioner and linear iterative solvers in Petsc. The final solution
vectors are saved in the parallel VTK format, which makes visualization using
ParaView straightforward.

Dendro-KT meshing framework: An “octant“ represents the fundamental building
block of the octree mesh. At the root of the octree, an octant is a cube that
encompasses the full domain. In order to construct the octree, all $p$ compute
processes start at the root node. To avoid communication during the refinement
stage, we perform redundant computations on all $p$ processes. Starting from
the root node, all processes refine (similar to the sequential code) until at
least $\mathcal{O}(p)$ octants are produced. Then using a weighted space-
filling-curve (SFC) based partitioning, we partition the octants across all
processes. Further, proceeding in a top-down fashion, an octant is refined
based on user-specified criteria. The refinement criteria is a user-specified
function that takes the coordinates of the octant and the associated values,
and returns true or false. Since the refinement happens in an element-local
fashion, this step is embarrassingly parallel. In distributed memory, the
initial top-down tree construction (which describes the mesh) also enables an
efficient partitioning of the domain across an arbitrary number of processes.

The key distinguishing factor for Dendro-KT, as compared to other octree
libraries [43, 62, 63], is the ability to carve out any arbitrarily shaped
object based on a simple Active \- InActive classification of a point
$\in\mathbb{R}^{d}$, leading to the generation of incomplete octree. An
element is said to be Active, if all the nodes of the element falls within the
active region of the domain $\Omega$ and InActive if all the nodes falls
within the inactive region. If some of the nodes fall within the active region
and some within the inactive, the element is marked Intercepted. The basic
idea of Dendro-KT is to retain the Active and Intercepted elements and remove
the InActive elements. Such a construction leads to generation of octree that
is incomplete. An octree is called incomplete if siblings of at least one
octant (at any given level) are missing. During the top-down construction, if
any given octant falls into the InActive region of the domain, then the given
octant is not further considered a candidate for refinement. This introduces a
very efficient way to generate octrees by pruning the tree at a coarser level.
Additionally, incomplete octrees provide saving in the total number of degrees
of freedom to be solved, in comparison to traditional complete octrees. This
translates to additional savings in time-to-solve as well as conditioning of
the matrix. We refer interested readers to our previous work [64, 52] for
further details of incomplete octree construction and scalability performance.

FEM operations: Dendro-KT supports both matrix and matrix-free computations.
To prevent indirect memory access, Dendro-KT does not store any elemental to
global map data structure but relies on performing top-down and bottom-up
traversals of the mesh tree. The top-down and bottom-up traversals can be
thought of as scatter and gather operations, respectively.

The top-down phase selectively copies nodes from coarser to finer levels until
the leaf level is reached. We create buckets for all child subtrees. Looping
through the nodes, a node is copied into a bucket if the node is incident on
the child subtree corresponding to that bucket. A node that is incident on
multiple child subtrees will be duplicated. By recursing on each child subtree
and its corresponding bucket of incident nodes, we eventually reach the leaf
level. Once the traversal reaches a leaf octant, the elemental nodes have been
copied into a contiguous array. The elemental vector is computed directly
without the use of an element-to-node map. The result is stored in a
contiguous output buffer the same size as the local elemental input vector.

After all child subtrees have been traversed, the bottom-up phase returns
results from a finer to a coarser level. The parent subtree nodes are once
again bucketed to child subtrees, but instead of the parent values being
copied, the values of nodes from each child are accumulated into a parent
output array. For any node that is incident on multiple child subtrees, the
values from all node instances are summed to a single value. The global vector
is assembled after the bottom-up phase executes at the root of the octree. A
detailed description of this process is provided in our previous work,
specifically Ishii et al. [51], Saurabh et al. [64], Khanwale et al. [50]

Scalability: We deploy the numerical methods on a scalable framework named
Dendrite-KT, which combines the octree-based meshing with continuous Galerkin
finite elements. We leverage the highly optimized and scalable implementations
of general finite element kernels in Dendrite-KT built on top of the scalable
Dendro-KT framework. We have shown massive scalability of Dendrite-KT for
various multiphysics phenomena in [51, 48, 52, 47, 53, 50, 54]. For example,
Saurabh et al. [54] recently showed the scalability of the multiphase flow
framework called Proteus based on Dendrite-KT till $\sim 114000$ MPI processes
on TACC Frontera.

Figure 2: Design schematic of the parallel framework for solving NS-PNP.

##  5 Results

###  5.1 Method of manufactured solutions for convergence analysis

Figure 3: Convergence analysis using the method of manufactured solutions. (a)
spatial, and (b) temporal convergence of the NS-PNP solver. The dashed line is
a reference line of slope 2.

We use the method of manufactured solutions to assess the temporal and spatial
convergence of our proposed scheme. In this work, we restrict ourselves to a
linear basis function ($r=1$). We construct the appropriate forcing function
to satisfy the analytical solution of the form:

$\begin{split}u&=\cos(2t)\sin(2\pi x)\cos(2\pi y)\\\ v&=-\cos(2t)\cos(2\pi
x)\sin(2\pi y)\\\ p&=\cos(2t)\sin(2\pi x)\cos(2\pi y)\\\
\phi&=-\cos(2t)\cos(2\pi x)\sin(2\pi y)\\\ c^{+}&=\cos(2t)\cos(2\pi
x)\sin(2\pi y)\\\ c^{-}&=\cos(2t)\sin(2\pi x)\cos(2\pi y)\end{split}$ (53)

We solve the NS-PNP equations on a square domain $[0,1]^{2}$, with Dirichlet
conditions enforced on all boundaries. We note that the manufactured solution
for the fluid velocity in eq. 53 is divergence-free. Here, we use a uniform
mesh with varying refinement levels for the convergence study. Fig. 3 shows
both the spatial and temporal discretization convergence results. For the
spatial convergence analysis, we fixed the timestep $\Delta t$ at 0.0157 and
varied the spatial resolution. Next, to study the temporal convergence
behavior, we chose the spatial resolution to be level 9 ($\Delta x=\Delta
y=1/2^{9}$). The $L_{2}$ error is reported at the final time of $t=\pi$. We
observe an expected second-order spatial ( Fig. 3(a)) and temporal convergence
( Fig. 3(a)). At an error of about $\mathcal{O}(10^{-7})$, we see the
flattening out of the error curves, especially in temporal convergence
studies. This flattening can be attributed to the errors from the spatial
resolution limit starting to dominate the $L_{2}$ error. These results
indicate that the implementation exhibits the expected second-order spatial
and temporal convergence behavior.

###  5.2 Electro-convective instability

We next demonstrate electrokinetic transport near an ion-selective surface,
which is a good benchmarking problem because 1) the problem can be found in a
wide range of applications, 2) the length scale of the EDL near the surface
requires an extremely small mesh size, which makes the problem a good test
case for the proposed octree-based mesh, and 3) The chaotic nature of the flow
triggers structures over a wide range of scales in all field quantities,
necessitating numerical stabilization from VMS. We compare our results against
the benchmark results of [16].

The simulation was performed for the transport of monovalent binary
electrolytes in a rectangular domain whose dimension is 8 $\times$ 1\. The top
boundary is a reservoir at which the potential and the concentration of both
species are fixed (Dirichlet). The ion-selective surface is located at the
bottom boundary, where the potential is grounded. The cation concentration was
$c^{+}=2$, and no anion flux was allowed across the ion-selective surface. A
periodic boundary condition is applied to the side walls. Fig. 4 shows the
schematic of the domain, and Table 3 lists the boundary conditions for the
different variables. We consider two cases, one when a moderate potential
difference of $\Delta\phi=20$ is imposed across the domain and another when a
large potential difference of $\Delta\phi=120$ is imposed. We anticipate
significant electro-convective instabilities to occur in the second case.

Figure 4: Schematic diagram (not on scale) for the electro-convective instability.  | Variable | Type | Value
---|---|---|---
NS | Top Wall | u | Dirichlet | 0
v | Dirichlet | 0
Bottom Wall | u | Dirichlet | 0
v | Dirichlet | 0
Side walls | Periodic
PNP | Top Wall | $\phi$ | Dirichlet | 20 or 120
$c_{+}$ | Dirichlet | 1
$c_{-}$ | Dirichlet | 1
Bottom Wall | $\phi$ | Dirichlet | 0
$c_{+}$ | Dirichlet | 2
$c_{-}$ | Zero flux | 0
Side walls | Periodic
Table 3: Boundary condition for the electro-convective instabilty $\mathbf{\Delta\phi}$ | Mesh size | # of MPI tasks | # of timestep | Final time | Wall time
---|---|---|---|---|---
$20$ |  $\sim$ 1M (8.3 M equivalent) |  224 CPU processes (4 nodes on TACC Frontera) | 1000 | 10 | 5 hours
$120$ | 12000 | 1.2 | 4 days
Table 4: Mesh details and time taken for the electro-convective instability.

Similar to the previous example, finer refinement was used near the ion-
selective surface to resolve steep gradients in the boundary layer. The bottom
part of the domain, extending from y = 0 to y = 0.1, consists of elements
refined to level 13 ($\Delta x=\Delta y=8/2^{13}$), while the rest of the
region consists of elements refined to level 10 ($\Delta x=\Delta
y=8/2^{10}$). The setup resulted in a domain with approximately 1M cells. In
contrast, an equivalent uniform mesh at the finest resolution would have
$\sim~{}10$M cells. The initial condition with fully developed EDL was
generated by solving only the PNP equations, then randomly perturbing the
resultant concentrations locally by 1% to initiate the instability before
solving the fully coupled NS-PNP equations. The timestep was $10^{-2}$ for the
case with potential difference, $\Delta\phi=20$, and a time step of $10^{-4}$
for the case of potential difference $\Delta\phi=120$. We ran the case of
$\Delta\phi=20$ until about 10 seconds when it reached steady-state
(approximately 1000 timesteps). For the case of $\Delta\phi=120$ we ran the
simulation for 1.2 seconds, which encompasses the time to reach statistical
stationarity and the time used to collect data (approximately 12000
timesteps). Table 4 briefs the mesh details along with the number of MPI tasks
and the wall time elapsed for both the cases.

We first compare the temporal evolution of the charge density ($c^{+}-c^{-}$)
and velocity fields under the imposition of the potential drops. This is shown
in Fig. 5 for the $\Delta\phi=20$ case and Fig. 6 for the $\Delta\phi=120$
case. Once a potential drop is applied, the ion concentration profile near the
ion-selective surface forms an electric double layer (EDL) and diffusion layer
(DL). Depending on the level of the potential drop, the extended space charge
layer (ESC) may form between EDL and DL [16]. When $\Delta\phi=20$, the ESC is
not formed, and the magnitude of convection is insufficient to trigger
instabilities. A series of counter-rotating vortices are observed in Fig. 5.
On the other hand, an ESC was observed in the initial stage of
$\Delta\phi=120$. The interplay between ESC and convection triggered
electrohydrodynamic instability, as reported in Druzgalski et al. [16]. As the
simulation proceeds, the development of electrohydrodynamic instability
becomes even more evident. We note that the average velocity at
$\Delta\phi=120$ was greater than that of $\Delta\phi=20$ case by two orders
of magnitude.

We next extract quantitative properties for comparison. The current measured
on the boundary is an often used property, as it is the easiest property to
measure experimentally. Fig. 7 (a) shows the temporal variation of the
magnitude of the net species flux, $|j_{T}^{\pm}|$, through the top boundary,
calculated as

$j_{T}^{\pm}=-\frac{1}{L_{T}}\sum_{l=1}^{N_{T}}\left[\left(q,\,\hat{n}_{T}\left(\frac{\partial
c^{+,t}}{\partial x_{i}}+z^{+}c^{+,t}\;\frac{\partial\phi^{t}}{\partial
x_{i}}\right)\right)_{\Gamma_{l}}-\left(q,\,\hat{n}_{T}\left(\frac{\partial
c^{-,t}}{\partial x_{i}}+z^{-}c^{-,t}\;\frac{\partial\phi^{t}}{\partial
x_{i}}\right)\right)_{\Gamma_{l}}\right],$ (54)

where ${L_{T}}$, $N_{T}$, $\hat{n}_{T}$, and
$\Gamma=\partial{\Omega_{\Gamma}}$ are the length of the top boundary, the
number of the top boundary elements, the normal vector to the top boundary,
and the top boundary, respectively. We see that, in the absence of electro-
convective instabilities, the current (for the $\Delta\phi=20$ case) was
uniform over time, as expected in the absence of ESC and not enough convection
to perturb established EDL. In contrast, significant current fluctuations were
observed for the $\Delta\phi=120$ case due to the instability.

We finally plot $x-$averaged field values for the potential (Fig. 7 (b)),
concentration (Fig. 7 (c)) and charge density (Fig. 7 (d)). Note the
logarithmic axis in these figures to ensure that the variation close to the
ion-selective surface is emphasized. These figures also allow detailed
comparison with previous DNS simulations of the same case by Druzgalski et al.
[16], which shows excellent agreement. Note that our numerical scheme allows
three orders of magnitude larger time steps as compared to the benchmarking
problem by Druzgalski et al. [16], with identical results.

|  |
---|---|---
t = 0.01 |  |
t = 0.1 |  |
t = 1.0 |  |
t = 5.0 |  |
Figure 5: Contour plots of charge density and velocity magnitudes at different times for $\Delta\phi=$ 20\.  |  |
---|---|---
t = 0.001 |  |
t = 0.005 |  |
t = 0.01 |  |
t = 0.05 |  |
t = 0.1 |  |
Figure 6: Contour plots of charge density and velocity magnitudes at different
times for $\Delta\phi=$ 120\.

(a) Variation of current with time

(b) Variation of potential gradient with $y$

(c) Variation of concentration with $y$

(d) Variation of charge density with $y$

Figure 7:  Electro-osmotic instability: Variation of current with time (in s)
for different $\Delta\phi$ is shown in Fig. 7(a). Time and X averaged
variation of potential gradient ($\partial_{y}\phi$) (Fig. 7(b)),
concentration (Fig. 7(c)) and charge density (Fig. 7(d)) with wall normal
distance $y$ for $\Delta\phi=$ 20 and 120. The dotted markers in the plot
corresponds to the previous DNS results by Druzgalski et al. [16]

###  5.3 Packed bed ion concentration polarization: Complex geometries, small
EDLs, and upstream vortex formation

Our final example illustrates a simulation of ion concentration polarization
(ICP) – a phenomenon used in a wide array of electrochemical unit operations –
in a realistic, state-of-art engineering device. The device consists of a
packed bed of metallic beads that act as a set of electrodes embedded within a
microfluidic device. This device was recently studied experimentally by
Berzina et al. [65] to circumvent two challenges faced by conventional ICP
occurring at a planar electrode, detailed in the next paragraph. The device
configuration offers a challenging canonical simulation example due to (a) the
complex geometries of the electrodes involved, (b) the steep gradients
(especially near these complex geometries), and small EDL’s produced by the
operating conditions, (c) the ensuing channel scale vortical structures that
have to be accurately resolved, and (d) the long time horizon over which the
current measurements stabilize.

There has been a growing interest in ICP due to its ability to enrich and
separate charged species. The key idea of ICP is to manipulate target analytes
by leveraging a high electric field in the ion depletion zone (IDZ) created by
selective charge transport. Previous studies have successfully demonstrated a
variety of applications, including seawater desalination [8], hemodialysis
[66], analyte preconcentration [67], and manipulation of neutral species
following their partition into ionic micelles [68], all using planar ICP.
However, two challenges preclude the scaling up of ICP from enabling
commercial deployment. First, an IDZ formed at a planar membrane or electrode
may not extend across the full channel cross-section under the flow rate
employed for focusing, thereby allowing the analyte to “leak” past the IDZ.
Second, within the IDZ, large fluid vortices lead to mixing, which decreases
the efficiency of analyte enrichment and worsens with increased channel
dimensions. A promising approach has been to move away from a planar electrode
(along the bottom of the channel) to using engineered, porous electrodes
across the channel cross-section. Recently, Berzina et al. [65] introduced 3-D
packed bed electrodes that successfully control unwanted instabilities and
extend high electric field area for the enrichment of bioanalytes before
detection. In a follow-up study by Berzina et al. [69], an additional bed of
packed bioconjugated beads were incorporated for DNA capture. We deploy our
computational framework on this configuration to explore the impact of the
packed bed electrodes on the current generation, concentration distribution,
electric field extent, and vorticity patterns.

The domain consists of a cuboidal geometry with packed bed electrodes – made
up of spherical beads – downstream of the channel inlet. Fluid is pushed from
left to right, and a potential difference is maintained. The length, width,
and height of the simulation domain are $32\times 1\times 1$. We only simulate
a part of the width considering symmetry. On the floor of the packed bead bed,
a planar electrode is located. The geometry around the packed beds is shown in
Fig. 9 \- Fig. 11.

A challenging aspect in such simulations is the need to construct an
adaptively refined mesh of this complex geometry. Here, the octree-based
adaptive meshing feature in our framework allows us to efficiently construct a
high-quality mesh of the system. Fig. 8 illustrates the mesh refinement near
the packed bed on a vertical plane. The refinement level ranges from 9 (bulk)
to 13 (electrodes), which corresponds to the element sizes ranging from
$32/2^{9}\sim 0.0625$ and $32/2^{13}\sim 0.0039$, which captures physics in a
wide range of length scales from the EDL around the ”rough” electrode surface
to the fluid flow at bulk. The non-dimensional Debye length, in this case, is
$0.0036$. The details of boundary conditions are listed in Table 5. As before,
the initial conditions were obtained by solving only the PNP equations until
the EDL was fully developed. Then, the fully coupled NS-PNP equations were
solved.

| Type | Value
---|---|---
NS |  Packed Bed Electrode, Top & Bottom Walls (no - slip) | u | Dirichlet | 0
v | Dirichlet | 0
w | Dirichlet | 0
p | Neumann | 0
Inlet | u | Dirichlet | 0.371
v | Dirichlet | 0
w | Dirichlet | 0
p | Neumann | 0
Outlet | u | Neumann | 0
v | Neumann | 0
w | Neumann | 0
p | Dirchlet | 0
Side Walls | Periodic
PNP |  Packed Bed & Planar Electrodes | $\phi$ | Dirichlet | 0
$c_{+}$ | Dirichlet | 0
$c_{-}$ | Zero Flux | 0
Inlet | $\phi$ | Dirichlet | 277.1
$c^{+}$ | Dirichlet | 1
$c^{-}$ | Dirichlet | 1
Outlet | $\phi$ | Neumann | 0
$c^{+}$ | Neumann | 0
$c^{-}$ | Neumann | 0
Top & Bottom Walls (except electrode) | $\phi$ | Neumann | 0
$c^{+}$ | Zero Flux | 0
$c^{-}$ | Zero Flux | 0
Side Walls | Periodic
Table 5: Boundary condition for 3-D electrode packed bed ICP

Table 6 briefs the mesh information and the total time taken. It is worth
noting the importance of adaptivity. The ability to selectively refine only
the regions of interest resulted mesh comprised of around 1.6 M elements and a
total wall time of approximately two days. In the absence of any adaptivity, a
resultant mesh with a refinement level of 13 would result in approximately 536
million elements. Such a drastic increase in the number of elements would
render the problem expensive in a reasonable amount of time.

Mesh size | Number of MPI task | Wall time
---|---|---
1.6 M (536 M equivalent) |  544 CPU cores (8 KNL nodes on Stampede2) | 2 days
Table 6: Mesh details and time taken for the Packed bed ion concentration
polarization.

(a)

(b)

(c)

(d)

(e)

Figure 8: Mesh refinement of pack bed ion concentration polarization.

(a)

(b)

Figure 9: The magnitude of the electric field (a) only with the plate
electrode, and (b) with added 3-D electrode bed.

(a)

(b)

Figure 10: The magnified view of the electric field near the packed bed from
the side (a) and top (b)

(a)

(b)

Figure 11: The flow streamline and magnitude around the electrodes. (a) only
with the plate electrode, and (b) with added 3-D electrode bed.

(a)

(b)

Figure 12: The ionic strength (average concentration in this case) (a) only
with the plate electrode, and (b) with added 3-D electrode bed.

We compare the spatial variations in the electric field, charge densities, and
flow fields between two cases: a microchannel with a packed bead electrode as
described above versus a microchannel with only a planar electrode. We first
compare the magnitude and extent of the electric field in Fig. 9. We see that
the presence of the microbeads causes the local enhancement of the electric
field to extend throughout the channel cross-section rather than being
concentrated only in the bottom region. We note that these observations are
qualitatively consistent with the ones reported in the experimental study of
Berzina et al. [65]. Figure. 10 shows a magnified view of the magnitude of the
electric field near the bead bed.

Next, we compare the flow fields between the two cases in Fig. 11. Note the
set of vortices that span the channel width, with each vortex pair nearly
spanning the channel height. This is in contrast to a single large vortex in
the planar electrode case. The increased region over which the electric field
acts induces a larger velocity field upstream of the packed bed. However,
experimental evidence suggests that this does not cause any electroconvective
instabilities. The absence of electroconvective instabilities is confirmed by
visualizing the ionic strength in Fig. 12.

##  6 Conclusions and future work

Direct numerical simulations of the coupled NS-PNP equations are challenging
due to the 1) wide ranges of time and length scales exhibited by the
multiscale physics, 2) The chaotic nature of the flow triggers structures over
a wide range of scales in all field quantities, and 3) geometric complexities
that escalate the computational cost. In this work, we presented a finite
element-based numerical framework that addresses those challenges by
incorporating 1) Variational multiscale pressure stabilization, 2) hybridized
semi-implicit and fully implicit time integration schemes, and 3) highly
parallelized adaptive meshing.

We demonstrate the utility of the framework using a series of numerical
examples. The temporal and spatial convergence tests were performed to verify
the order of accuracy of the numerical framework. Next, we performed
simulations capturing electrohydrodynamic instabilities near an ion-selective
surface. These simulations showed excellent agreement with previously reported
direct numerical simulations. Moreover, the proposed numerical scheme enabled
using time steps three orders of magnitude larger than the benchmarking study,
allowing significant savings in computational costs. Finally, we deployed the
framework for simulating a practical microfluidics application for analyte
preconcentration, which involves a complex geometry. The adaptive meshing and
numerical approach are able to capture the impact of the complex geometry of
3-D electrodes. The simulation confirmed the effect of the 3-D electrodes on
suppressing unwanted vortices and extending the high electric field area, and
this phenomenon was also observed in the experimental work.

In summary, our proposed numerical framework addresses the challenges involved
in detailed simulations of electrokinetic phenomena for practical
applications. The numerical framework offers robustness for simulations with
numerical instabilities, the ability to take large time steps to reach long-
time horizon simulations, and adaptive octree-based meshing to resolve
geometric complexities, all packaged within a massively scalable software
stack.

##  7 Acknowledgements

The authors acknowledge XSEDE grant number TG-CTS110007 for computing time on
TACC Stampede2. The authors also acknowledge computing allocation through a DD
award on TACC Frontera and Iowa State University computing resources. BG, KS,
SK, and MAK were funded in part by NSF grants 1935255, and 1855902.

##  Appendix A Details of solver selection for the numerical experiments

For the cases presented in sections 5.2 and 5.3 we use the BiCGStab linear
solver (a Krylov space solver) with additive Schwarz-based preconditioning.
For better reproduction, the command line options we provide petsc are given
below which include some commands used for printing some norms as well. The
petsc options for the Navier-Stokes equations are:

⬇

solver_options_ns = {

ksp_atol = 1e-10

ksp_rtol = 1e-08

ksp_stol = 1e-14

ksp_max_it = 2000

ksp_type = "bcgs"

pc_type = "asm"

# residual monitoring

ksp_monitor = ""

ksp_converged_reason = ""

};

and for Poisson-Nernst-Planck equations, we use the Petsc SNES for the Newton
iteration. We used the BiCGStab linear solver (a Krylov space solver) with
additive Schwarz-based preconditioning for the inner linear solves of the
Newton iteration.

⬇

solver_options_pnp = {

ksp_atol = 1e-8

ksp_rtol = 1e-8

ksp_stol = 1e-8

#multigrid

ksp_type = "fgmres"

pc_type = "gamg"

pc_gamg_asm_use_agg = True

mg_levels_ksp_type = "gmres"

mg_levels_pc_type = "sor"

mg_levels_ksp_max_it = 10

};

We used default tolerance for SNES convergence.

## References

  * Rubinstein et al. [1997] I. Rubinstein, B. Zaltzman, O. Kedem,  Electric fields in and around ion-exchange membranes,  Journal of membrane science 125 (1997) 17–21.
  * Fleury et al. [1993] V. Fleury, J.-N. Chazalviel, M. Rosso,  Coupling of drift, diffusion, and electroconvection, in the vicinity of growing electrodeposits,  Physical Review E 48 (1993) 1279.
  * Fleury et al. [1994] V. Fleury, J. Kaufman, D. Hibbert,  Mechanism of a morphology transition in ramified electrochemical growth,  Nature 367 (1994) 435–438.
  * Huth et al. [1995] J. M. Huth, H. L. Swinney, W. D. McCormick, A. Kuhn, F. Argoul,  Role of convection in thin-layer electrodeposition,  Physical Review E 51 (1995) 3444.
  * Kwak et al. [2013] R. Kwak, K. M. Lim, J. Han, et al.,  Shear flow of an electrically charged fluid by ion concentration polarization: Scaling laws for electroconvective vortices,  Physical review letters 110 (2013) 114501.
  * Wessling et al. [2014] M. Wessling, L. G. Morcillo, S. Abdu,  Nanometer-thick lateral polyelectrolyte micropatterns induce macrosopic electro-osmotic chaotic fluid instabilities,  Scientific reports 4 (2014) 1–5.
  * Chi and Yu [2018] J. Chi, H. Yu,  Water electrolysis based on renewable energy for hydrogen production,  Chinese Journal of Catalysis 39 (2018) 390–394.
  * Kim et al. [2010] S. J. Kim, S. H. Ko, K. H. Kang, J. Han,  Direct seawater desalination by ion concentration polarization,  Nature nanotechnology 5 (2010) 297–301.
  * Kwak et al. [2016] R. Kwak, V. S. Pham, B. Kim, L. Chen, J. Han,  Enhanced salt removal by unipolar ion conduction in ion concentration polarization desalination,  Scientific reports 6 (2016) 1–11.
  * Nikonenko et al. [2014] V. V. Nikonenko, A. V. Kovalenko, M. K. Urtenov, N. D. Pismenskaya, J. Han, P. Sistat, G. Pourcelly,  Desalination at overlimiting currents: State-of-the-art and perspectives,  Desalination 342 (2014) 85–106.
  * Lees et al. [2022] E. W. Lees, B. A. Mowbray, F. G. Parlane, C. P. Berlinguette,  Gas diffusion electrodes and membranes for co2 reduction electrolysers,  Nature Reviews Materials 7 (2022) 55–64.
  * Weng et al. [2018] L.-C. Weng, A. T. Bell, A. Z. Weber,  Modeling gas-diffusion electrodes for co 2 reduction,  Physical Chemistry Chemical Physics 20 (2018) 16973–16984.
  * Kim et al. [2020] S. Kim, B. Ganapathysubramanian, R. K. Anand,  Concentration enrichment, separation, and cation exchange in nanoliter-scale water-in-oil droplets,  Journal of the American Chemical Society 142 (2020) 3196–3204.
  * Anand et al. [2022] R. K. Anand, S. Kim, B. Ganapathysubramanian, Concentration enrichment, separation and cation exchange in water-in-oil droplets, 2022\. US Patent App. 17/450,931.
  * Devarakonda et al. [2022] S. Devarakonda, S. Kim, B. Ganapathysubramanian, P. Shrotriya,  Designing asymmetrically modified nanochannel sensors using virtual eis,  Electrochimica Acta 403 (2022) 139694.
  * Druzgalski et al. [2013] C. Druzgalski, M. Andersen, A. Mani,  Direct numerical simulation of electroconvective instability and hydrodynamic chaos near an ion-selective surface,  Physics of Fluids 25 (2013) 110804.
  * Kim et al. [2017] J. Kim, I. Cho, H. Lee, S. J. Kim,  Ion concentration polarization by bifurcated current path,  Scientific reports 7 (2017) 1–12.
  * Dydek et al. [2011] E. V. Dydek, B. Zaltzman, I. Rubinstein, D. Deng, A. Mani, M. Z. Bazant,  Overlimiting current in a microchannel,  Physical review letters 107 (2011) 118301.
  * Han et al. [2014] J.-H. Han, E. Khoo, P. Bai, M. Z. Bazant,  Over-limiting current and control of dendritic growth by surface conduction in nanopores,  Scientific reports 4 (2014) 1–8.
  * Nam et al. [2015] S. Nam, I. Cho, J. Heo, G. Lim, M. Z. Bazant, D. J. Moon, G. Y. Sung, S. J. Kim,  Experimental verification of overlimiting current by surface conduction and electro-osmotic flow in microchannels,  Physical review letters 114 (2015) 114501.
  * Mani and Wang [2020] A. Mani, K. M. Wang,  Electroconvection near electrochemical interfaces: Experiments, modeling, and computation,  Annual Review of Fluid Mechanics 52 (2020) 509–529.
  * Davidson et al. [2014] S. M. Davidson, M. B. Andersen, A. Mani,  Chaotic induced-charge electro-osmosis,  Physical review letters 112 (2014) 128302.
  * Davidson et al. [2016] S. M. Davidson, M. Wessling, A. Mani,  On the dynamical regimes of pattern-accelerated electroconvection,  Scientific reports 6 (2016) 1–10.
  * Jia and Kim [2014] M. Jia, T. Kim,  Multiphysics simulation of ion concentration polarization induced by nanoporous membranes in dual channel devices,  Analytical chemistry 86 (2014) 7360–7367.
  * Gong et al. [2018] L. Gong, W. Ouyang, Z. Li, J. Han,  Direct numerical simulation of continuous lithium extraction from high mg2+/li+ ratio brines using microfluidic channels with ion concentration polarization,  Journal of membrane science 556 (2018) 34–41.
  * Seo et al. [2022] M. Seo, W. Kim, H. Lee, S. J. Kim,  Non-negligible effects of reinforcing structures inside ion exchange membrane on stabilization of electroconvective vortices,  Desalination 538 (2022) 115902.
  * Karatay et al. [2015] E. Karatay, C. L. Druzgalski, A. Mani,  Simulation of chaotic electrokinetic transport: Performance of commercial software versus custom-built direct numerical simulation codes,  Journal of colloid and interface science 446 (2015) 67–76.
  * Roberts et al. [2014] S. A. Roberts, V. E. Brunini, K. N. Long, A. M. Grillet,  A framework for three-dimensional mesoscale modeling of anisotropic swelling and mechanical deformation in lithium-ion electrodes,  Journal of The Electrochemical Society 161 (2014) F3052.
  * Roberts et al. [2016] S. A. Roberts, H. Mendoza, V. E. Brunini, B. L. Trembacki, D. R. Noble, A. M. Grillet,  Insights into lithium-ion battery degradation and safety mechanisms from mesoscale simulations using experimentally reconstructed mesostructures,  Journal of Electrochemical Energy Conversion and Storage 13 (2016).
  * Kim et al. [2022] S. Kim, M. A. Khanwale, R. K. Anand, B. Ganapathysubramanian,  Computational framework for resolving boundary layers in electrochemical systems using weak imposition of dirichlet boundary conditions,  Finite Elements in Analysis and Design 205 (2022) 103749.
  * Volker [2016] J. Volker, Finite Element Methods for Incompressible Flow Problems, volume 51 of Springer Series in Computational Mathematics Book, Springer, 2016.
  * Lim et al. [2007] J. Lim, J. Whitcomb, J. Boyd, J. Varghese,  Transient finite element analysis of electric double layer using nernst–planck–poisson equations with a modified stern layer,  Journal of Colloid and Interface Science 305 (2007) 159–174.
  * Lu et al. [2010] B. Lu, M. J. Holst, J. A. McCammon, Y. Zhou,  Poisson–nernst–planck equations for simulating biomolecular diffusion–reaction processes i: Finite element solutions,  Journal of computational physics 229 (2010) 6979–6994.
  * Sun et al. [2016] Y. Sun, P. Sun, B. Zheng, G. Lin,  Error analysis of finite element method for poisson–nernst–planck equations,  Journal of Computational and Applied Mathematics 301 (2016) 28–43.
  * Gao and He [2017] H. Gao, D. He,  Linearized conservative finite element methods for the nernst–planck–poisson equations,  Journal of Scientific Computing 72 (2017) 1269–1289.
  * He and Sun [2018] M. He, P. Sun,  Mixed finite element analysis for the poisson–nernst–planck/stokes coupling,  Journal of Computational and Applied Mathematics 341 (2018) 61–79.
  * Liu and Maimaitiyiming [2021] H. Liu, W. Maimaitiyiming,  Efficient, positive, and energy stable schemes for multi-d poisson–nernst–planck systems,  Journal of Scientific Computing 87 (2021) 1–36.
  * Liu et al. [2022] H. Liu, Z. Wang, P. Yin, H. Yu,  Positivity-preserving third order dg schemes for poisson–nernst–planck equations,  Journal of Computational Physics 452 (2022) 110777.
  * Dehghan et al. [2022] M. Dehghan, Z. Gharibi, R. Ruiz-Baier,  Optimal error estimates of coupled and divergence-free virtual element methods for the poisson–nernst–planck/navier–stokes equations,  arXiv preprint arXiv:2207.02455 (2022).
  * Chen et al. [2022] Q. Chen, X. Liu, Y. Lei, H. Zhu,  An electrokinetic preconcentration trapping pattern in electromembrane microfluidics,  Physics of Fluids 34 (2022) 092009.
  * Ishii et al. [2019] M. Ishii, M. Fernando, K. Saurabh, B. Khara, B. Ganapathysubramanian, H. Sundar,  Solving pdes in space-time: 4d tree-based adaptivity, mesh-free and matrix-free approaches,  in: SC’19: Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis, SC ’19, Association for Computing Machinery, New York, NY, USA, 2019, pp. 1–61. URL: https://doi.org/10.1145/3295500.3356198. doi:doi:10.1145/3295500.3356198.
  * Sundar et al. [2008] H. Sundar, R. Sampath, G. Biros,  Bottom-up construction and 2:1 balance refinement of linear octrees in parallel,  SIAM Journal on Scientific Computing 30 (2008) 2675–2708. doi:doi:10.1137/070681727.
  * Burstedde et al. [2011] C. Burstedde, L. C. Wilcox, O. Ghattas,  p4est: Scalable algorithms for parallel adaptive mesh refinement on forests of octrees,  SIAM Journal on Scientific Computing 33 (2011) 1103–1133. doi:doi:10.1137/100791634.
  * Fernando et al. [2019] M. Fernando, D. Neilsen, H. Lim, E. Hirschmann, H. Sundar,  Massively parallel simulations of binary black hole intermediate-mass-ratio inspirals,  SIAM Journal on Scientific Computing 41 (2019) C97–C138. URL: https://doi.org/10.1137/18M1196972. doi:doi:10.1137/18M1196972. arXiv:https://doi.org/10.1137/18M1196972.
  * Fernando et al. [2017] M. Fernando, D. Duplyakin, H. Sundar,  Machine and application aware partitioning for adaptive mesh refinement applications,  in: Proceedings of the 26th International Symposium on High-Performance Parallel and Distributed Computing, HPDC ’17, ACM, New York, NY, USA, 2017, pp. 231–242. URL: http://doi.acm.org/10.1145/3078597.3078610. doi:doi:10.1145/3078597.3078610.
  * Khanwale et al. [2020] M. A. Khanwale, A. D. Lofquist, H. Sundar, J. A. Rossmanith, B. Ganapathysubramanian,  Simulating two-phase flows with thermodynamically consistent energy stable cahn-hilliard navier-stokes equations on parallel adaptive octree based meshes,  Journal of Computational Physics 419 (2020) 109674.
  * Khanwale et al. [2022] M. A. Khanwale, K. Saurabh, M. Fernando, V. M. Calo, H. Sundar, J. A. Rossmanith, B. Ganapathysubramanian,  A fully-coupled framework for solving cahn-hilliard navier-stokes equations: Second-order, energy-stable numerical methods on adaptive octree based meshes,  Computer Physics Communications 280 (2022) 108501.
  * Saurabh et al. [2021] K. Saurabh, B. Gao, M. Fernando, S. Xu, M. A. Khanwale, B. Khara, M.-C. Hsu, A. Krishnamurthy, H. Sundar, B. Ganapathysubramanian,  Industrial scale large eddy simulations with adaptive octree meshes using immersogeometric analysis,  Computers & Mathematics with Applications 97 (2021) 28–44.
  * Xu et al. [2021] S. Xu, B. Gao, A. Lofquist, M. Fernando, M.-C. Hsu, H. Sundar, B. Ganapathysubramanian,  An octree-based immersogeometric approach for modeling inertial migration of particles in channels,  Computers & Fluids 214 (2021) 104764.
  * Khanwale et al. [2023] M. A. Khanwale, K. Saurabh, M. Ishii, H. Sundar, J. A. Rossmanith, B. Ganapathysubramanian,  A projection-based, semi-implicit time-stepping approach for the cahn-hilliard navier-stokes equations on adaptive octree meshes,  Journal of Computational Physics 475 (2023) 111874.
  * Ishii et al. [2019] M. Ishii, M. Fernando, K. Saurabh, B. Khara, B. Ganapathysubramanian, H. Sundar,  Solving pdes in space-time: 4d tree-based adaptivity, mesh-free and matrix-free approaches,  in: Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis, 2019, pp. 1–61.
  * Saurabh et al. [2022] K. Saurabh, M. Ishii, M. A. Khanwale, H. Sundar, B. Ganapathysubramanian,  Scalable adaptive algorithms for next-generation multiphase simulations,  arXiv preprint arXiv:2209.12130 (2022).
  * Khanwale et al. [2022] M. A. Khanwale, K. Saurabh, M. Ishii, H. Sundar, B. Ganapathysubramanian,  Breakup dynamics in primary jet atomization using mesh-and interface-refined cahn-hilliard navier-stokes,  arXiv e-prints (2022) arXiv–2209.
  * Saurabh et al. [2023] K. Saurabh, M. Ishii, M. A. Khanwale, H. Sundar, B. Ganapathysubramanian,  Scalable adaptive algorithms for next-generation multiphase flow simulations,  in: 2023 IEEE International Parallel and Distributed Processing Symposium (IPDPS), IEEE, 2023, pp. 590–601.
  * Kirby [2010] B. Kirby, Micro- and Nanoscale Fluid Mechanics: Transport in Microfluidic Devices, Cambridge University Press, 2010. URL: https://books.google.com/books?id=y7PB9f5zmU4C.
  * Druzgalski et al. [2013] C. L. Druzgalski, M. B. Andersen, A. Mani,  Direct numerical simulation of electroconvective instability and hydrodynamic chaos near an ion-selective surface,  Physics of Fluids 25 (2013) 110804–110804. doi:doi:10.1063/1.4818995.
  * Tezduyar et al. [1992] T. Tezduyar, S. Mittal, S. Ray, R. Shih,  Incompressible flow computations with stabilized bilinear and linear equal-order-interpolation velocity-pressure elements,  Computer Methods in Applied Mechanics and Engineering 95 (1992) 221–242.
  * Hughes et al. [2018] T. J. Hughes, G. Scovazzi, L. P. Franca,  Multiscale and stabilized methods,  Encyclopedia of Computational Mechanics Second Edition (2018) 1–64.
  * Hughes [1995] T. J. Hughes,  Multiscale phenomena: Green’s functions, the Dirichlet-to-Neumann formulation, subgrid scale models, bubbles and the origins of stabilized methods,  Computer Methods in Applied Mechanics and Engineering 127 (1995) 387–401.
  * Bazilevs et al. [2007] Y. Bazilevs, V. Calo, J. Cottrell, T. Hughes, A. Reali, G. Scovazzi,  Variational multiscale residual-based turbulence modeling for large eddy simulation of incompressible flows,  Computer Methods in Applied Mechanics and Engineering 197 (2007) 173–201.
  * Hughes et al. [2000] T. J. Hughes, L. Mazzei, K. E. Jansen,  Large eddy simulation and the variational multiscale method,  Computing and Visualization in Science 3 (2000) 47–59.
  * Zhang et al. [2019] W. Zhang, A. Almgren, V. Beckner, J. Bell, J. Blaschke, C. Chan, M. Day, B. Friesen, K. Gott, D. Graves, et al.,  Amrex: a framework for block-structured adaptive mesh refinement,  Journal of Open Source Software 4 (2019) 1370–1370.
  * Popinet [2003] S. Popinet,  Gerris: a tree-based adaptive solver for the incompressible euler equations in complex geometries,  Journal of computational physics 190 (2003) 572–600.
  * Saurabh et al. [2021] K. Saurabh, M. Ishii, M. Fernando, B. Gao, K. Tan, M.-C. Hsu, A. Krishnamurthy, H. Sundar, B. Ganapathysubramanian,  Scalable adaptive pde solvers in arbitrary domains,  in: Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis, 2021, pp. 1–15.
  * Berzina et al. [2022] B. Berzina, S. Kim, U. Peramune, K. Saurabh, B. Ganapathysubramanian, R. K. Anand,  Out-of-plane faradaic ion concentration polarization: stable focusing of charged analytes at a three-dimensional porous electrode,  Lab on a Chip (2022).
  * Berzina and Anand [2018] B. Berzina, R. K. Anand,  An electrokinetic separation route to source dialysate from excess fluid in blood,  Analytical chemistry 90 (2018) 3720–3726.
  * Ko et al. [2012] S. H. Ko, Y.-A. Song, S. J. Kim, M. Kim, J. Han, K. H. Kang,  Nanofluidic preconcentration device in a straight microchannel using ion concentration polarization,  Lab on a Chip 12 (2012) 4472–4482.
  * Berzina and Anand [2019] B. Berzina, R. K. Anand,  Continuous micellar electrokinetic focusing of neutral species driven by ion concentration polarization,  Lab on a Chip 19 (2019) 2233–2240.
  * Berzina et al. [2023] B. Berzina, U. Peramune, S. Kim, K. Saurabh, E. L. Claus, M. E. Strait, B. Ganapathysubramanian, R. K. Anand,  Electrokinetic enrichment and label-free electrochemical detection of nucleic acids by conduction of ions along the surface of bioconjugated beads,  ACS sensors 8 (2023) 1173–1182.