Articles | Volume 16, issue 16
Development and technical paper
24 Aug 2023
Development and technical paper |  | 24 Aug 2023

Validating the Nernst–Planck transport model under reaction-driven flow conditions using RetroPy v1.0

Po-Wei Huang, Bernd Flemisch, Chao-Zhong Qin, Martin O. Saar, and Anozie Ebigbo

Reactive transport processes in natural environments often involve many ionic species. The diffusivities of ionic species vary. Since assigning different diffusivities in the advection–diffusion equation leads to charge imbalance, a single diffusivity is usually used for all species. In this work, we apply the Nernst–Planck equation, which resolves unequal diffusivities of the species in an electroneutral manner, to model reactive transport. To demonstrate the advantages of the Nernst–Planck model, we compare the simulation results of transport under reaction-driven flow conditions using the Nernst–Planck model with those of the commonly used single-diffusivity model. All simulations are also compared to well-defined experiments on the scale of centimeters. Our results show that the Nernst–Planck model is valid and particularly relevant for modeling reactive transport processes with an intricate interplay among diffusion, reaction, electromigration, and density-driven convection.

1 Introduction

Reactive transport processes are fundamental in large-scale subsurface applications such as geological hydrogen storage (Aftab et al.2022; Hassanpouryouzband et al.2022), geothermal energy extraction (Randolph and Saar2011; Fleming et al.2020; Ezekiel et al.2022), biogeochemical subsurface treatments (Carrera et al.2022), contaminant remediation (Frizon et al.2003; Reddy and Cameselle2009), recovery of minerals (Chang et al.2021), and carbon sequestration and storage (Luhmann et al.2013; Tutolo et al.2014; Pogge von Strandmann et al.2019; Grimm Lima et al.2020; Ling et al.2021). On a smaller scale, reactive transport governs the processes of steel corrosion in reinforced concrete (Yu et al.2020), cementation processes (Samson and Marchand2007; Cochepin et al.2008), dendrite growth in microelectronics (Illés et al.2022), and electrochemical reactors (Rivera et al.2021). In all aforementioned settings, the reacting fluid consists of multiple ionic species and may be subject to an external electric field. In some applications, the electromigration of ionic species can play a significant role in the overall transport process. A typical model that accounts for electromigration is the Nernst–Planck equation.

Modeling electrochemical processes using the Nernst–Planck equation is usually coupled with the Poisson equation for electrical potential, also known as the Poisson–Nernst–Planck (PNP) model. The PNP model has been applied in various electrochemical applications (Jasielec2021). In geoscientific applications, benchmarks of the Nernst–Planck model have been proposed and verified considering charged clay surfaces (Tournassat and Steefel2021). Charged surfaces affect the water permeability of a porous material (in contrast to air permeability, which is not affected by electric charges) (Revil and Leroy2004; Kwon et al.2004a, b; Cheng and Milsch2020). The PNP model is utilized at the pore scale to determine such water permeabilities (Priya et al.2021). Pore-scale experiments with external charge effects have also been conducted to investigate reactive transport processes (Sprocati et al.2019; Sprocati and Rolle2022; Rolle et al.2022; Izumoto et al.2022). Rasouli et al. (2015) verified numerical software using electromigration benchmarks (Lichtner1994; Glaus et al.2013) and multicomponent electromigration experiments considering dispersion effects (Rolle et al.2013). Finally, multicomponent transport experiments in porous media confirm the validity of the Nernst–Planck model (Muniruzzaman et al.2014; Muniruzzaman and Rolle2015, 2017; Rolle et al.2018).

The above benchmark and validation studies focused on reactions and the mixing of fluids, which are subject to flow introduced by fluid pressure or electric fields. In this work, we propose another set of benchmark problems for numerical models, based on chemically driven convection experiments (Avnir and Kagan1984; Eckert and Grahn1999; Eckert et al.2004), which focus on how reaction and mixing of fluids influence fluid flow. The concept of reaction-driven flow is as follows: consider reaction A+BC, where A and B are the reactants and C is the product. If A, B, and C have different material properties, for example, viscosity, diffusivity, and/or density, the reaction can induce fluid flow and hydrodynamic instabilities in such systems (De Wit2016, 2020). We are interested in aqueous systems, where ionic species exist, such that the Nernst–Planck model is applicable.

Zalts et al. (2008) performed experiments using an acid and a base to study chemically driven convection. Related experiments and theories of the underlying mechanism have emerged: unequal species diffusivities (Citri et al.1990; Almarcha et al.2010, 2011; Trevelyan et al.2011; Lemaigre et al.2013; Kim2019; Jotkar et al.2021), the same species diffusivity but with varying viscosity (Hejazi and Azaiez2013), and concentration-dependent species diffusivities (Bratsun et al.2015, 2017, 2021, 2022). Chemically driven convection has also been investigated using a weak acid and a weak base (Cherezov et al.2018). Since the Nernst–Planck model is expected to be able to account for unequal species diffusivities and has the characteristics of concentration-dependent diffusive fluxes, we examine the validity of the Nernst–Planck model by comparing our Nernst–Planck simulations with chemically driven experiments conducted by others.

In addition, we present simulations of convective dissolution of CO2, a process which is important in the underground storage of CO2 (Singh et al.2019). Thomas et al. (2016) performed experiments of convective dissolution of CO2 into alkaline brines, including LiOH and NaOH, and observed that different cations lead to varying onset times and fingering patterns. Thomas et al. (2016) proposed that the inert cation, which does not participate in the acid–base reaction between the dissolved CO2 and the brine, plays a major role in the development of convective instabilities. Since dissolved CO2 forms charged and uncharged states (CO2(aq), HCO3-, CO32-) when the pH changes as a consequence of neutralization reactions, it is natural to apply the Nernst–Planck model in such cases.

In the following sections, we first establish the fundamental equations of the conservation laws and the accompanying numerical methods. We then provide details concerning the two reaction-driven flow experiments and a CO2 dissolution experiment – all conducted by others. We compare the differences between the simulations of the single-diffusivity transport model and the Nernst–Planck model and elucidate the necessity of considering electromigration effects caused by differing diffusivities during reactive transport modeling.

2 Materials and methods

2.1 Fundamental equations for reactive transport processes

In this section, we delineate the fundamental equations describing the reactive transport of aqueous species.

2.1.1 Mass balance of the fluid components

We consider a single-phase fluid composed of N components, which can be ionic species, dissolved gas, or minerals. The single-phase fluid occupies a certain volume Vf, and the ith component of the fluid has the mass mi. The mass concentration of the ith component is defined as

(1) ρ i = m i V f .

We use the term mass concentration to clarify that ρi is not the density of the component. The fluid density, ρ, can be obtained by summing up the mass concentrations of the components, i.e.,

(2) ρ = i = 1 N ρ i = i = 1 N m i V f .

In general, the fluid density is a function of fluid pressure, p; temperature, T; and fluid composition (i.e., the mass concentrations, ρi). The equation of state of the fluid density can be written as

(3) ρ = ρ ( p , T , ρ i ; i = 1 , 2 , , N ) .

The mass balance of each fluid component is described by the continuity equation,

(4) ρ i t + ( ρ i v i ) = Q i ,

where ρivi is the mass flux and Qi is the source or sink term introduced by chemical reactions. Since chemical reactions preserve fluid mass, we have iQi=0. We express the mass flux of the ith component as

(5) ρ i v i = ρ i u + ρ i u EP , i + J diff , i ,

where u is the barycentric velocity, uEP,i is the electrophoretic velocity, and Jdiff,i is the diffusive flux. In the following sections, we define the models of the diffusive flux, the electrophoretic flux, and the barycentric flux.

2.1.2 Modeling of the diffusive flux

We use standard Fickian diffusion to model the diffusive flux,

(6) J diff , i = - D i M i C i i = 1 , 2 , , N - 1 ,

where D is the diffusivity, M is the molar mass, and C is the molar concentration. The contribution of the diffusing solvent can be neglected, considering ρi/ρN1. Such an assumption leads to a particular choice of the Nth component, which contributes the most mass to the fluid. In our experiments of interest, the Nth component is H2O(l). Fickian diffusion utilizes the diagonal diffusivity matrix assumption, which is valid for many electrolytes (Dreyer et al.2013). A diagonal diffusivity matrix means that short-range interactions between ions are negligible, which restricts our model to dilute solutions (Kontturi et al.2008).

2.1.3 Modeling of the barycentric flux

The type of experiments of interest for this work is performed in a Hele-Shaw cell, for which fluid flow can be described by Darcy's law,

(7) u = - k η p - ρ g ,

where k is the permeability, η is the dynamic viscosity of the fluid, and g is the gravitational acceleration. When combined with the incompressibility assumption,

(8) ( ρ u ) = 0 ,

u and p can be obtained. The scalar permeability of a Hele-Shaw cell, with a gap width of w, is given by

(9) k = w 2 12 .

2.1.4 Modeling of the electrophoretic flux

When the components in the fluids are electrically charged, for example, the ions Na+ or Cl, we have to consider the effect of Coulombic forces among the ions. Such Coulombic forces are typically modeled by incorporating the electric field E into the mass fluxes of the fluid components (Boudreau et al.2004). The electrophoretic flux is described by

(10) ρ i u EP , i = z i ρ i D i F R T E ,

where zi is the charge number of the fluid component, F is the Faraday constant, R is the ideal gas constant, and T is the fluid temperature. We require the charge conservation equation to define the electric field in a multicomponent fluid. To obtain the charge conservation equation of the ionic species, it is useful to use the flux defined in molar concentration. Using ρi=MiCi, we have the mass conservation equation (without considering the source or sink terms introduced by reactions)

(11) C i t = - ( C i u + C i u EP , i + J diff , i / M i ) .

Following the derivations of Kirby (2010), the charge density is defined as

(12) ρ E = F i = 1 N - 1 z i C i .

The charge conservation equation is the sum of Eq. (11) over all ionic species multiplied by the charge number and the Faraday constant,

(13) ρ E t = - F i = 1 N - 1 z i C i u + z i C i u EP , i + z i J diff , i / M i .

We assume the aqueous solution is locally electroneutral,

(14) i = 1 N - 1 z i C i = 0 ,

which is an approximation considering the Debye length (typically on the order of nanometers) vanishes at the length scale of the considered system (Dickinson et al.2011). The electroneutrality assumption eliminates the barycentric flux term in the charge conservation equation. We consider Fickian diffusive flux, Eq. (6), and the electrophoretic flux, Eq. (10), so that the charge conservation is given by

(15) ρ E t = - i = 1 N - 1 D i C i ( z i F ) 2 R T E - i = 1 N - 1 D i z i F C i .

The multiplier of the electric field, E, is the electrical conductivity (C2m-3kg-1s),

(16) σ = i = 1 N - 1 D i C i ( z i F ) 2 R T .

When the diffusivities of all species are equal, Di=D, the diffusive flux in the charge conservation equation is

(17) D i = 1 N - 1 z i F C i = D F i = 1 N - 1 z i C i = 0 .

Then, the charge conservation equation yields Laplace's equation of the electric potential, i.e., (σϕ)=0, where E=-ϕ. However, the diffusivity of the chemical species varies with molecular size and can differ by up to an order of magnitude. Such variations in molecular diffusivity lead to diffusive fluxes that produce an electric field. The model that combines mass conservation, Eq. (11), and Poisson's equation of electrostatics, is known as the Poisson–Nernst–Planck (PNP) model (Pamukcu2009).

The PNP model aims at resolving both the electric potential and the molar concentrations subject to the boundary conditions of the electric potential. The mathematical properties of the PNP model, with no external flow (Filipek et al.2017), coupled with Darcy flow (Herz et al.2012; Ignatova and Shu2022), and coupled with Navier–Stokes flow (Lee2021; Constantin et al.2021, 2022) are still being explored. Rigorous upscaling and homogenization of the PNP–Stokes system has been performed (Ray et al.2012a; Kovtunenko and Zubkova2021), and the upscaled formulations are examined by numerical modeling (Frank et al.2011; Ray et al.2012b). Numerical methods and solution techniques of the PNP equations are an ongoing research field (Flavell et al.2014; Song et al.2018; Shen and Xu2021; Liu and Maimaitiyiming2021; Yan et al.2021; Liu et al.2022; Zhang et al.2022). Models similar to the PNP model also arise in modeling multicomponent interdiffusion of solids (Bożek et al.2019), where the experimental comparison and the development of numerical methods are studied by Sapa et al. (2020) and Bożek et al. (2021).

In this work, we enforce the null current condition

(18) i = 1 N - 1 D i C i ( z i F ) 2 R T E - i = 1 N - 1 D i z i F C i = 0 ,

to ensure that no charge accumulates at the continuum scale of interest,

(19) ρ E t = 0 .

This result is a simplification of the PNP model. The null current condition has been utilized in modeling multicomponent ionic transport (Lichtner1985; Cappellen and Gaillard1996; Giambalvo et al.2002; Muniruzzaman and Rolle2019; Cogorno et al.2022; López-Vizcaíno et al.2022). The combined assumptions of local electroneutrality and null current are referred to as strict electroneutrality by Lees et al. (2017). Furthermore, if there are no sources of the electric field (Tabrizinejadas et al.2021), then the electric field can be represented by

(20) E = R T F j = 1 N - 1 D j z j C j k = 1 N - 1 ( z k ) 2 D k C k ,

where we use j and k as summation indices. Hence, the molar electrophoretic flux of individual species can be formulated as

(21) C i u EP , i = z i D i C i F R T E = z i D i C i j = 1 N - 1 z j D j C j k = 1 N - 1 ( z k ) 2 D k C k .

There are many physical interpretations of such a flux. For example, one is that the dissolved species is advecting due to the electric field. Another point of view can be obtained by combining Eq. (21) with the Fickian diffusion flux,


which corresponds to the “natural generalization” of Fickian diffusion, denoted by Onsager (1945). Following Boudreau et al. (2004), the cross-coupling diffusivities are given by

(24) D i j = δ i j D i - z i z j D i D j C i k = 1 N - 1 ( z k ) 2 D k C k = δ i j D i - z i D i C i z j D j k = 1 N - 1 ( z k ) 2 D k C k ,

where Einstein's summation convention is not applied. Hence, locally electroneutral diffusion is equivalent to nonlinear multicomponent diffusion with a diffusivity tensor being a rational function of molar concentrations of the charged species (Rubinstein1990). Similar forms of the cross-coupling diffusivities can also be found in Lichtner (1985).

Combining Eqs. (11) and (22) yields

(25) C i t = - - D i C i + z i D i C i j = 1 N - 1 z j D j C j k = 1 N - 1 ( z k ) 2 D k C k + u C i ,

which is the Nernst–Planck equation for single-phase multicomponent systems. If all species have the same diffusivity, Di=D, and the aqueous solution is electroneutral, Eq. (14) is used. Then, the single-phase multicomponent Nernst–Planck equation reduces to

(26) C i t = - - D C i + u C i ,

which we refer to as the single-diffusivity model. Both the Nernst–Planck model and the single-diffusivity model satisfy the electroneutrality condition, Eq. (14), and the zero-charge accumulation assumptions, Eq. (19). However, of the two models, only the Nernst–Planck model is capable of handling the physics of varying species diffusivities. To communicate the effectiveness of the Nernst–Planck equation in reactive transport scenarios, we compare the modeling results with the single-diffusivity model, which is a common approach. We set the diffusivity of the single-diffusivity model as 3×10-3 mm2 s−1, which is within the common values for species diffusivities in water at 25 C.

2.1.5 Chemical equilibrium using Gibbs energy minimization

Another aspect of mass balance amongst the fluid components, which we now denote as the chemical species, is the reaction that changes the amount of the chemical species while conserving the mass balance of chemical elements. The reaction source or sink term in the mass balance of each fluid component, Qi, is implicitly formulated as the Gibbs energy minimization problem,

(27) min n G = i = 1 N n i μ i subject to  mass balance of chemical elements n 0 ,

where G is the Gibbs energy and ni is the molar amount of the ith chemical species. The chemical potential is defined as

(28) μ i = μ i + R T log a i ,

where μi is the standard chemical potential of the ith species and ai is the activity of the ith species. For the standard chemical potential, we employ the SUPCRT07 database (Johnson et al.1992). To evaluate the activities, we use the Helgeson–Kirkham–Flowers (HKF) activity model (Helgeson and Kirkham1974a, b, 1976; Helgeson et al.1981), which is based on the extended Debye–Hückel activity model. For the activity model of CO2(aq), we use the Drummond model (Drummond1981). The combination of the SUPCRT database and the HKF activity model is utilized in many geochemical modeling packages, such as EQ3/6 (Wolery1992), CHILLER (Reed1998), ChemApp (Petersen and Hack2007), Perplex (Connolly2009), GEM-Selektor (Kulik et al.2012), PHREEQC (Parkhurst and Appelo2013), Geochemist's Workbench (Bethke2007), and Reaktoro (Leal2022) (Miron et al.2019). General reviews of the activity models can be found in Marini (2007) and Hingerl (2012).

We use Reaktoro to solve the chemical equilibrium problem. The numerical methods and implementations of the Gibbs energy minimization problem are detailed in Leal et al. (2014, 2016, 2017). The equilibrium calculations are performed at every time step and for every discretized volume of the simulation domain, and calculations of similar chemical states can be cached for efficiency. Efficient lookup tables and machine learning approaches of chemical-equilibrium calculations are in active research (Huang et al.2018; De Lucia and Kühn2021; De Lucia et al.2021; Savino et al.2022; Laloy and Jacques2022), and significant speedups in solving reactive transport problems have been achieved (Kyas et al.2022; Bordeaux-Rego et al.2022). In particular, Reaktoro employs on-demand machine learning and physics-based interpolation, which is able to conserve mass and control interpolation accuracy (Leal et al.2020). Surrogate models of the full reactive transport process using machine learning have also been developed (Sprocati and Rolle2021). However, to avoid approximation errors, machine learning approaches and surrogate models are not employed in this work.

2.2 Numerical schemes

This section elaborates on the numerical schemes we use for solving the fundamental equations of the reactive transport processes.

2.2.1 The barycentric flux

We utilize the mixed finite-element formulation to solve for fluid velocity and fluid pressure under the constraints of mass balance, Eq. (8), and momentum balance, Eq. (7). To avoid spurious solutions, the finite-element spaces have to satisfy the Ladyzhenskaya–Babuška–Brezzi (LBB) condition. We refer the reader to other works for further discussions regarding the LBB condition in the context of fluid flow problems (Donea and Huerta2003), geodynamics (Thieulot and Bangerth2022), and elasticity and electromagnetic applications (Boffi et al.2013). We select LBB-compatible elements, namely the zeroth-order Raviart–Thomas (RT) element for the barycentric flux and the piecewise constant element for the pressure (Raviart and Thomas1977). Another popular LBB-compatible element is the Brezzi–Douglas–Marini element (Brezzi et al.1985). The mixed finite-element formulation can be written in matrix form as

(29) A B T B 0 u ^ p ^ = f 0 ,

where u^ and p^ are the solution vectors for the barycentric velocity and the fluid pressure. Such a formulation leads to a saddle point problem. Following Benzi et al. (2005)'s review of saddle point problems, one can apply Uzawa's method (Uzawa1958) to solve such problems. Given an initial condition of p^n at step n, Uzawa's method solves u^ and p^ iteratively,


where ω>0 is a relaxation parameter. Normally, Uzawa's method converges slowly and requires many iterations, meaning that improving Uzawa's method is an ongoing research effort (Bacuta2006). We use one of the improved methods, the augmented Lagrangian Uzawa's method, introduced by Fortin and Glowinski (1983):

(32) A + r B T B B T B 0 u ^ p ^ = f 0 ,

where r>0 is a tuning parameter. We solve the matrix problem iteratively,


When ω=r, we can choose sufficiently large r to accelerate the convergence of Uzawa's method. However, the matrix A+rBTB becomes more ill-conditioned as r increases, which may lead to a large number of iterations to solve when applying an iterative method (Fortin and Glowinski1983). Hence, we use the MUMPS (Amestoy et al.2001, 2019) direct solver to update the velocity, Eq. (33). The degrees of freedom of u^ in our simulations range from 17 000 to 140 000.

2.2.2 Transport of fluid components

In the previous section, we selected an LBB-compatible velocity and pressure pair. The basis function of pressure is a piecewise constant, which can be described as the average pressure of a certain cell volume in a given mesh. Therefore, it is straightforward to define the molar concentrations as piecewise constants. We use the finite-volume method to discretize the transport of the fluid components, Eq. (25):

(35) Ω q i C i t d Ω = - Γ q i ( - D i C i + z i D i C i j = 1 N - 1 z j D j C j k = 1 N - 1 ( z k ) 2 D k C k + u C i ) n d Γ ,

where the fluxes are yet to be defined. Considering no-flux boundary conditions of the fluid components, we use the two-point flux approximation scheme (TPFA) for the diffusive flux,

(36) Γ q i D i C i n d Γ = - Γ int [ q i ] D i [ C i ] h d Γ ,

where Γint denotes the interior cell interfaces, [•] is the jump operator that evaluates the difference in a function across a common interface of two cells, and h is the distance between the cell centers. We assumed that the diffusivity of each species is constant. For diffusivities varying in space, Tournassat et al. (2020) discussed and compared various averaging approaches of the diffusivities at the interfaces between cells. For the advective flux, we utilize the full upwinding method,

(37) - Γ q i C i u n d Γ = - Γ int [ q i ] C i up u n d Γ ,

where Ciup is the molar concentration in the upwind direction. For the Nernst–Planck term, we follow Tournassat et al. (2020)'s approach:

(38) - Γ q i z i D i C i j = 1 N - 1 z j D j C j k = 1 N - 1 ( z k ) 2 D k C k n d Γ = Γ int [ q i ] z i D i { C i } j = 1 N - 1 z j D j [ C j ] k = 1 N - 1 ( z k ) 2 D k { C k } h d Γ ,

where {} is the averaging operator. We use the arithmetic average to represent the molar concentration at the cell interface, which is suitable when all charged components are mobile (Gimmi and Alt-Epping2018) or when there are no membrane effects or large concentration gradients (Tournassat et al.2020). Unified Form Language (UFL; Alnæs2012) is utilized for defining the variational forms in our code implementation.

With regard to the time-stepping schemes, we use an explicit scheme for the upwind advection term and the Crank–Nicolson scheme for the diffusion and Nernst–Planck terms. The time step size Δt is determined in an adaptive manner, where the minimum and maximum time step sizes are prescribed for each simulation. To avoid negative concentration values, we perform a transformation of variables Ci=exp (ci), which leads to the logarithmically transformed Nernst–Planck equations (Kirby2010). The transport equations are assembled using FEniCS (Alnæs et al.2015) and solved using PETSc (Balay et al.2022).

2.2.3 Coupling flow, transport, and chemical equilibrium

Of many coupling approaches in reactive transport modeling (Steefel and MacQuarrie2018; Carrayrou et al.2004; Abd and Abushaikha2021), we apply the sequential non-iterative approach (SNIA) to couple flow, transport, and chemical equilibrium. Figure 1 illustrates the simulation procedure, where the fluid flow, the species transport, and the equilibrium steps are solved consecutively.

Due to the simplicity and effectiveness in applying the SNIA coupling between the transport solvers and the chemical equilibrium codes, this coupling method is widely adopted in the following software: OpenGeoSys (Shao et al.2009; Naumov et al.2022), poreReact (coupling of OpenFOAM (Weller et al.1998) and Reaktoro) (Oliveira et al.2019), CSMP++GEM (Yapparova et al.2019), FEniCS–Reaktoro (Damiani et al.2020), Osures (Moortgat et al.2020), FEniCS-based Hydro-Mechanical-Chemical solver (Kadeethum et al.2021), PorousFlow (based on the MOOSE Framework; Permann et al.2020) (Wilkins et al.2021), IC-FERST-REACT (Yekta et al.2021), COMSOL and PHREEQC (Jyoti and Haese2021), GeoChemFoam (coupling of OpenFOAM and PHREEQC) (Maes and Menke2021), coupling of Reaktoro and Firedrake (Rathgeber et al.2016) (Kyas et al.2022), and P3D-BRNS (Golparvar et al.2022). For reviews of reactive transport codes and the underlying coupling approaches, we refer the reader to the publications by Gamazo et al. (2015), Damiani et al. (2020).

Figure 1Flowchart of the simulation procedures that couple flow, transport, and chemical equilibrium using the sequential non-iterative approach (SNIA).


2.3 Selected experiments for model evaluation

We select three experiments, performed (by others) in a Hele-Shaw cell, to evaluate the effectiveness of the Nernst–Planck model in reactive transport situations. The diffusivities of the aqueous species are listed in Table 1. The density and dynamic viscosity of the relevant aqueous solutions are listed in Table 2. Our models assume constant dynamic viscosity, set as the average viscosity of the aqueous solutions. All experiments are performed in a lab environment, and we assume a constant temperature of 25 C and a constant background pressure of 100 kPa. No-flow boundary conditions are prescribed on all sides of the simulation domain. Table 3 specifies the chemical compositions of the aqueous solutions, which defines the initial conditions of the mass transport equations.

Table 1Diffusivities of aqueous species in 25 C water at infinite dilution.

The diffusivities of aqueous species other than HCO3-, CO32-, and CO2(aq) are taken from Vanýsek (2021). The diffusivities of HCO3-, CO32-, and CO2(aq) are calculated using Zeebe (2011).

Download Print Version | Download XLSX

Table 2Dynamic viscosity and density of aqueous solutions and water at 25 C.

Linearly interpolated using a Nishikata et al. (1981), b Sipos et al. (2000), c Zaytsev and Aseyev (1992), and d Åkerlöf and Teare (1938). The viscosity and density of water are taken from Lemmon and Harvey (2021).

Download Print Version | Download XLSX

Table 3The chemical compositions of the aqueous solutions. The chemical compositions are in units of mol L−1. Dashes denote “not applicable” because the chemical species are not considered in the associated simulation cases.

Download Print Version | Download XLSX

2.3.1 Chemically driven convection of acid–base systems

Almarcha et al. (2010) experimentally showed that simple acid–base reactions, such as A+BC, can lead to hydrodynamic instabilities due to differences in diffusivities and densities of the chemical species. For the first case in our numerical modeling study, we follow the experimental setup described in Almarcha et al. (2010). The acid (A) is 1 M HCl and the base (B) is 1 M NaOH, where M is the molar concentration (mol L−1). The experiment was performed by Almarcha et al. (2010) in a Hele-Shaw cell, which was 3.1 cm wide, 5 cm high, and had a 0.5 mm gap width. To prevent Rayleigh–Taylor instability, the less dense acid was placed on top of the base. We consider the following species as the main fluid components in this case: H+, OH, Na+, Cl, and H2O(l). Please refer to Fig. 2 for the setup. Table 3 lists the chemical composition of 1 M HCl and 1 M NaOH. The prescribed minimum time step size is 3×10-2 s, and the maximum time step size is 2.0 s.

For the second modeling case, we choose a similar experiment with a strikingly different instability process. In their chemically driven convection experiments, Bratsun et al. (2017) reported a shockwave-like structure of the acid–base interface. In further studies, this phenomenon has subsequently been reproduced using many different acid–base pairs (Mizev et al.2021; Bratsun et al.2021). For our numerical comparison study, we use the acid–base pair from Bratsun et al. (2021), i.e., HNO3NaOH. The acid consists of 1.5 M HNO3 and is placed on top of 1.4 M NaOH. The Hele-Shaw cell of this experiment has a width of 2.5 cm, a height of 9.0 cm, and a gap width of 1.2 mm. We consider the following species as the main fluid components in this case: H+, OH, Na+, NO3-, and H2O(l). The setup is shown in Fig. 2, and Table 3 lists the initial chemical composition of 1.5 M HNO3 and 1.4 M NaOH. The prescribed minimum time step size is 1×10-3 s, and the maximum time step size is 1.0 s. We perform the simulation in a 2D geometry with the same width, height, and species initial conditions as in the experiment and compare the simulation results with the time series of interferometry snapshots taken during the experiments.

2.3.2 Convective dissolution of CO2 in reactive alkaline solutions

Thomas et al. (2016) performed experiments of convective dissolution of CO2 in reactive alkaline solutions. The experiment setup consisted of a Hele-Shaw cell saturated with alkaline solutions, where gaseous CO2 flows over the top of the Hele-Shaw cell. Eventually, gaseous CO2 dissolves into the alkaline solution, leading to density changes in the solution and, as a result, the formation of convection cells. Furthermore, Thomas et al. (2018) studied the convective dissolution of CO2 in various molar concentrations of NaCl solutions using the aforementioned experimental setup. Likewise, Mahmoodpour et al. (2019) studied the effect of various salt solutions on the onset of convective instabilities by conducting experiments within porous media. Similar experimental studies in a sandbox have been performed using dense miscible fluids (Neufeld et al.2010; Agartan et al.2015; Guo et al.2021; Wang et al.2021; Tsinober et al.2022). Although such experiments do capture the effects of density-driven flow, electromigration effects are not involved. Other experimental studies of convective dissolution of gaseous CO2 in aqueous solutions have been conducted in pressure–volume–temperature (PVT) cells (Yang and Gu2006), Hele-Shaw cells (Kneafsey and Pruess2010; Class et al.2020; Zhang et al.2020; Jiang et al.2020; Teng et al.2021), and 3D granular porous media (Mahmoodpour et al.2020; Brouzet et al.2022). Amarasinghe et al. (2020) reviewed CO2 dissolution experiments and performed experiments of convective dissolution of CO2 in a Hele-Shaw cell and porous media under supercritical CO2 conditions.

Numerical studies of CO2 convective dissolution have been performed considering the reactive transport of a single component (Azin et al.2013), a single component with mineral reactions (Babaei and Islam2018; Shafabakhsh et al.2021), a single component in two fluid phases with a resolved gas–fluid boundary (Martinez and Hesse2016), multiple components with mineral reactions (Zhang et al.2011), and multiple components in two fluid phases with mineral reactions (Audigane et al.2007; Soltanian et al.2019; Sin and Corvisier2019). However, the Nernst–Planck model was not invoked in any of these works.

Using lattice Boltzmann methods with consideration of electrostatic forces, Fu et al. (2020, 2021) investigated how ionic species of various salt solutions with varying concentrations affect the onset time and mass transfer rate of CO2 convective dissolution. So far, these are the only works we found that account for electromigration effects when modeling the convective dissolution of CO2.

In the third case of our numerical study, we model the experiment of Thomas et al. (2016) and investigate the reactive transport of CO2 in a Hele-Shaw cell including the electromigration of species using the Nernst–Planck model. Instead of comparing our results to the experiments at the full scale (165 mm by 210 mm) (Thomas et al.2016), a simulation domain with a width of 100 mm, a height of 40 mm, and a gap width of 0.5 mm is utilized; see the right panel of Fig. 2. The smaller simulation domain was chosen, as the relevant processes remained within this smaller domain. On the top boundary of the simulation domain, we assume that the gaseous CO2 and the alkaline brine are always in equilibrium. The alkaline brine is 0.01 M LiOH. We consider the following species as the main fluid components in the case of convective dissolution of CO2: H+, OH, Li+, CO2(aq), HCO3-, CO32-, and H2O(l), and the initial conditions are presented in Table 3. The prescribed minimum and maximum time step sizes are 5×10-3 and 1.5 s, respectively. In the results section, we present the development of the fingers and the movement of the “shock wave” induced by density, reaction, and electromigration effects.

Figure 2The selected experimental setups. On the left, the chemically driven convection experiments with 1 M HCl–1 M NaOH (Almarcha et al.2010) are shown. In the middle, the chemically driven convection setup with 1.5 M HNO3–1.4 M NaOH (Mizev et al.2021) is shown. On the right, a smaller setup of CO2 dissolution with 0.01 M LiOH solutions compared to the original setup (165 mm by 210 mm) designed by Thomas et al. (2016) is shown From left to right, the gap widths of the Hele-Shaw cells are 0.5 mm, 1.2 mm, and 0.5 mm, respectively. All Hele-Shaw cells are closed, meaning that no fluid leaves the domain.


3 Results

This section compares the simulations and the experimental results of the chemically driven convection problems and the convective dissolution of CO2 into a reactive alkaline solution. All simulations are performed using the Nernst–Planck and the single-diffusivity models.

3.1 Comparing the chemically driven convection experiments with the single-diffusivity model and the Nernst–Planck model

We present the results of the chemically driven convection between 1 M HCl and 1 M NaOH solutions. Figure 3 compares the digital interferometry images obtained from experiments with the simulated fluid densities. The interferometry images have a size of 13 mm by 23.075 mm, from which we select a sub-region of the same size as the simulations. This region is indicated by the red box in Fig. 3. Note that the interferometry images show the density contrast of the fluids (induced by changes in refractive index) rather than the density itself. Hence, the simulated density plots cannot be compared with the interferometry images pixel by pixel. Nonetheless, in the left panel of Fig. 3, only the simulations with the Nernst–Planck model replicate the experimentally observed fingering effect with a similar number of fingers. At the simulation time of 60 s, a low-density layer (shown in yellow) emerges in the Nernst–Planck model, an effect which cannot be observed in the single-diffusivity model. At 115 s, the low-density layer develops instabilities, causing density-driven flow. However, the time it takes for fingering to commence is longer in the simulations than in the experiments. This delay can be attributed to the finite representation by the numerical scheme of the diffusion and electromigration fluxes at the initially sharp interface between the acid and the base, which could be improved by adaptive mesh refinement or a higher mesh resolution.

Figure 4 compares the images from the physical experiment with our numerical simulations that utilize the Nernst–Planck model with a lag time of 40 s. We further present the filled contours of the Na+ molar concentration and the molar concentration of Na+ and Cl combined. Na+ is transported upwards, closely following the low-density fingers. The rightmost plots of Fig. 4 show that the low-density region is caused by lower concentrations of both Na+ and Cl. The dark red region indicates lower salinity and corresponds to the low-density yellow and turquoise region. More interestingly, a stable region of higher salinity (shown in blue) forms close to the initial contact line of the acid and the base.

Figure 3Comparing the simulations using the Nernst–Planck model and those using the single-diffusivity model with the experiments of chemically driven convection in HClNaOH solutions. The column on the right shows the full simulation domain (31 mm by 50 mm), and the red box indicates the comparison window. For the four columns on the left, the time increases from left to right, and the simulation time corresponds to the time shown in the experiment snapshots. The experimental figures are based on digital interferometry, and we adapted the figures from Almarcha et al. (2010).

Figure 4Comparing the simulation results using the Nernst–Planck model with the experimental results of chemically driven convection in HClNaOH solutions. The leftmost column shows the experimental figures, adapted from Almarcha et al. (2010). The filled contour plots present the simulation results, showing the fluid density, molar concentration of Na+, and the sum of the molar concentrations of Na+ and Cl with a lag time of 40 s.

Figure 5 shows a comparison of the experiments of 1.5 M HNO3 and 1.4 M NaOH solutions and the simulations. Both the Nernst–Planck model and the single-diffusivity model reproduce the shockwave-like structures seen in the experimental interferometry images, i.e., the rapid downward movement of the acid–base interface. In the interferometry images at the 3 s snapshot, fingers appear, which the simulations cannot reproduce at the 3 s simulation time. As discussed before, during the HClNaOH experiments the mismatch is due to the poor approximation of diffusive fluxes at the initially sharp interface. For the later snapshots, 150 and 700 s, we observe the dark lines becoming more sparse in the interferometry images. This indicates that the fluid densities become more homogeneous, and the simulations replicate such mixing physics.

In Fig. 6, we show the simulations of the molar concentrations of Na+ and the streamlines of the fluid velocity field at 25, 50, 75, 100, and 125 s. One can clearly see that the Na+ plumes follow the streamlines. The maximum velocity in the y direction of the Nernst–Planck simulations and the single-diffusivity simulations during the five presented snapshots are 4.086 and 4.184 mm s−1, respectively. Both models replicate the experiments, and no particular difference can be observed. In such experimental settings, where density-driven convection is dominant, a distinctive description of the diffusive fluxes is less crucial.

Figure 5Comparing the simulations using the Nernst–Planck model and those using the single-diffusivity model with the experiments of chemically driven convection in HNO3NaOH solutions. The column on the right shows the full simulation domain (25 mm by 90 mm), and the red box indicates the comparison window. In the left columns, the time increases from left to right, where the snapshots correspond to 3, 150, and 700 s. The figures from the experiments are based on digital interferometry, and they were adapted from Mizev et al. (2021).

Figure 6Comparing the simulations using the Nernst–Planck model and those using the single-diffusivity model in modeling chemically driven convection of HNO3NaOH solutions. The top and the bottom rows show the simulations employing the Nernst–Planck model and the simulations using one single diffusivity, respectively. The time increases from left to right, and the snapshots correspond to 25, 50, 75, 100, and 125 s.


3.2 Simulation of convective dissolution of CO2 in reactive alkaline solutions

Figure 7 compares the schlieren images of convective dissolution of CO2 in 0.01 M LiOH solution from the experiments with the density plots using the Nernst–Planck and the single-diffusivity models. The onset time of the Nernst–Planck model and the single-diffusivity model are approximately 600 and 1400 s, respectively. Therefore, the densities of the two simulation results are shown for the onset time plus the experimental snapshot, 780 s. At these time steps, the convective fingers in both models reach similar depths of ∼17 mm, as shown in the schlieren image from the experiment. The simulation result of the Nernst–Planck model is comparable to the experiment of CO2 dissolving into the LiOH brine and shows a broader range of finger sizes compared to those of the single-diffusivity model. In contrast, the convective fingers in the single-diffusivity model resemble the fingers seen in the schlieren image, where CO2 dissolves in H2O and the electromigration effects are weaker compared to the LiOH solution.

Figure 8 shows the density plots of the simulation results using the Nernst–Planck model and the single-diffusivity model. In the top row, we show the onset of the convective instability. The single-diffusivity model develops instabilities roughly 800 s later than the Nernst–Planck model. The second panel from the top shows the convective fingers after their onset; the Nernst–Planck model yields twice the number of fingers compared to the single-diffusivity model. In these two panels, we observe that the Nernst–Planck model generates a low-density layer, colored in yellow, at the dissolution front, which causes a shorter onset time and more fingers due to a larger density contrast. In the two bottom panels, we show the evolution of the simulations at later times. There are differences in the fingering pattern, such as more secondary convective fingers being generated by the Nernst–Planck model. Moreover, the low-density layer in the Nernst–Planck model persists even at later times.

To further quantify the simulation results, Fig. 9 visualizes the pH plots and the molar concentrations of Li+, CO2(aq), and CO32-. We use a diverging color palette to present the acid regions (orange) and the base regions (purple). The molar concentrations of CO2(aq) and CO32- follow the pH color palette, where most CO2(aq) exists in the acidic region and most CO32- exists in the basic region. Such changes in the CO2 states lead to changes in the charge states z=0,1,2 for CO2(aq), HCO3-, and CO32-, respectively. In the Li+ plot of the Nernst–Planck simulations, we observe up to 30 % differences in the molar concentration compared to the initial condition (0.01 mol L−1). The red regions are where Li+ accumulates, and we can clearly see the interaction between Li+ and CO32-, whereby Li+ migrates towards the double-charged CO32-. Meanwhile, the single-diffusivity simulation shows almost no contrast in the filled contours other than due to numerical errors. The difference between the maximum and minimum molar concentrations of Li+ is 10-8 mol L−1. Using the Nernst–Planck model, we confirm the influence of the inert (spectator) ions in the reactive-transport situation discussed in Thomas et al. (2016).

Figure 7Comparing the simulations using the Nernst–Planck model and those using the single-diffusivity model with the experiments of convective dissolution of CO2 in 0.01 M LiOH solution. We plot the experimental figure of CO2 dissolving in pure water in the top-left to demonstrate the differences in fingering patterns compared with CO2 dissolving in LiOH solution. The filled contours show the fluid density. The experimental figures are adapted from Thomas et al. (2016).

Figure 8Comparing the simulations using the Nernst–Planck model and those using the single-diffusivity model, while numerically modeling the convective dissolution of CO2 in 0.01 M LiOH solution. The filled contours show the fluid density. Time increases from top to bottom and is indicated in the plots.


Figure 9Comparing the simulations using the Nernst–Planck model and those using the single-diffusivity model, while numerically modeling convective dissolution of CO2 in 0.01 M LiOH solution. The top row shows the pH values. The orange palette represents acidic regions (pH<7), and the purple palette represents basic regions (pH >7). The gray region shows the unperturbed part of the LiOH solution (pH >11.5). In the bottom row, we show the zoomed in view of the red boxes in the top panel. We plot the molar concentrations of Li+, CO2(aq), and CO32-. Notice that the molar concentrations of CO2(aq) and CO32- share the same color bar as the pH values. The vertical red line in the color bar indicates the minimum molar concentrations of CO2(aq) and CO32- shown in the plots.


4 Discussion

In this section, we discuss whether the Nernst–Planck model is valid, and its differences are compared to the commonly used single-diffusivity model.

4.1 The validity of the Nernst–Planck model under reactive flow conditions

Comparing to two chemically driven convection experiments conducted by Almarcha et al. (2010) and Mizev et al. (2021), we show that the Nernst–Planck numerical model is able to reproduce the convective fingering observed in the experiments, especially in the 1 M HCl–1 M NaOH case. The mechanism that leads to convection is known as the reactive diffusive layer convection (DLC) instability, which arises when a less dense solution of a fast-diffusing component overlies a denser solution of a slow-diffusing component (Lemaigre et al.2013). In our Nernst–Planck model, the fast-diffusing component is H+. Additionally, the convective plumes of the reactive DLC instability are only observed above the initial contact line between the acid and the base (Almarcha et al.2010; Lemaigre et al.2013). This is the asymmetric feature of the reactive DLC instability, which can also be observed in the Nernst–Planck simulations (Fig. 3). To directly quote Lemaigre et al. (2013), “the asymmetry of the DLC pattern is related to the chemical reaction which consumes the acid before it can accumulate in the lower layer and replaces it by a third species of intermediate contribution to density”. This third species is salt, a product of the acid–base neutralization reaction. Using the Nernst–Planck model, we not only replicate this stabilizing lower layer but also find that the accumulation of salt can exceed the salt concentration at chemical equilibrium (in this case, approximately 0.5 mol L−1), as shown in Fig. 4.

In the 1.5 M HNO3–1.4 M NaOH case, both the Nernst–Planck model and the single-diffusivity model reproduce the shockwave-like structure observed in the experiments performed by Mizev et al. (2021). The instability in this case is convection controlled (CC) (Bratsun et al.2017). The underlying mechanism of this CC instability has previously been attributed to unequal and concentration-dependent diffusivities (Bratsun et al.2017, 2021, 2022). However, using the single-diffusivity model, our simulations can reproduce the experimental observations. This indicates that modeling CC instabilities is less sensitive to the specific diffusivity model used, as the CC instability is mostly caused by the density-driven convective flow of the reactants. Nevertheless, there are no drawbacks (other than being computationally more intensive) in employing the Nernst–Planck equations to model the CC instability cases.

In the case of CO2 dissolving in 0.01 M LiOH solution, we observe a much later convective onset time in the numerical simulations than in the experiment by Thomas et al. (2016). This discrepancy is much larger than in the 1 M HCl–1 M NaOH case. It cannot only be explained by the ineffectiveness of the diffusive-flux approximations at the sharp interface of dissolving CO2(aq). Convective instability can be introduced in a variety of ways in laboratory experiments. For example, interfacial instabilities between the gas and the liquid can occur due to the gas flowing over the interface. Furthermore, if the injected CO2 is not saturated with water, evaporation of water can induce gravitational instability (Bringedal et al.2022). These localized physical processes at the gas–liquid boundary, which are not included in our numerical models, could be partly responsible for the later onset times we observe in our simulations.

Although the comparison between the experiment and the simulations of the two models in Fig. 7 can be ambiguous, i.e., the Nernst–Planck model is not the model that describes the schlieren images, we emphasize that the difference between the experiments with Li+ and those without it is apparent in the figure. When comparing the Nernst–Planck and the single-diffusivity models, the former exhibits an earlier fingering onset time. This is caused by electromigration effects, which generate the low-density layer shown in Fig. 8. The electromigration of Li+, an inert ion in this case, and subsequent effects on the overall transport processes shown in Fig. 9, is clearly a motivation for more studies on the electromigration effects of inert (spectator) ions.

To summarize, we use the experiments conducted by Almarcha et al. (2010), Thomas et al. (2016), Mizev et al. (2021) and our numerical simulation results to show that the Nernst–Planck model is valid in the considered reactive transport scenarios. Particularly in the 1 M HCl–1 M NaOH experiments, where the convective fingers are controlled by diffusion. In the convection-controlled regime, electromigration effects are less crucial, and both the Nernst–Planck model and the single-diffusivity model replicate the 1.5 M HNO3–1.4 M NaOH experiments well. In the case of CO2 dissolving in an alkaline solution, the use of the Nernst–Planck model uncovers intricate details of how inert ionic species can affect the overall reactive transport processes.

4.2 Electromigration and the relevance of inert ionic species

Using the simulations of chemically driven convection in 1 M HCl–1 M NaOH solutions, we elaborate the main differences between the Nernst–Planck model and the single-diffusivity model. The top panel of Fig. 10 shows the molar concentrations of salt, plotted against elevation in the Hele-Shaw cell. As expected, the single-diffusivity model shows the concentration profile of two diffusing Heaviside step functions of Na+ and Cl. We plot the sum of Na+ and Cl molar concentrations to demonstrate that in single-diffusivity models, Na+ and Cl exhibit an ideal behavior, known as equimolar counter-diffusion. The Nernst–Planck model shows unexpected results, such as the accumulation of Na+ and Cl close to the acid–base boundary. The accumulation can be attributed to a larger diffusivity of the acid components, H+ (9.311×10-3 mm2 s−1) and Cl (2.032×10-3 mm2 s−1), compared to the base components, OH (5.273×10-3 mm2 s−1) and Na+ (1.334×10-3 mm2 s−1), which leads to more Cl migrating into the base. One might get the impression that this observation can be modeled by employing multiple diffusivities in the single-diffusivity model. However, we emphasize that the Nernst–Planck model has to be utilized to maintain electroneutrality at the continuum scale. Closely inspecting the fluid density profile, the deficit in salt concentration leads to a lower density than the acid, which causes the reactive DLC instability. The single-diffusivity model is unable to depict the accumulation or deficit in the salt concentrations.

Figure 10Comparing the simulations using the Nernst–Planck model and those using the single-diffusivity model in modeling chemically driven convection in 1 M HCl–1 M NaOH solutions. The snapshot is given at t=60 s, when fingers have not yet appeared (see Fig. 3). We plot the molar concentrations of salt and the fluid density against elevation in the Hele-Shaw cell, focusing on the contact region between the acid and the base (25 mm).


Figure 11Comparing the simulations using the Nernst–Planck model and those using the single-diffusivity model in modeling convective dissolution of CO2 in 0.01 M LiOH solution. The snapshot is given at t=600 s, when fingers have not yet appeared (see Fig. 8). We plot the molar concentrations of Li+, OH, CO32-, HCO3-, and CO2(aq) and the fluid density against elevation in the Hele-Shaw cell, focusing on the regions close to the top boundary of the Hele-Shaw cell. The horizontal dashed gray lines in the bottom panel show the minimum density observed in the Nernst–Planck model at t=600 s.


Figure 12Simulation results using the Nernst–Planck model, showing the molar concentration of Li+ during the convective dissolution of CO2 in 0.01 M LiOH solution.


Moreover, the only reaction that occurs in this experiment is the neutralization reaction

(39) H + + OH - H 2 O ( l ) .

No reaction occurs for Na+ and Cl. Hence, they are inert (spectator) species. We elaborate on the transport of H+, given by

(40) C H + t = - - D H + H + C H + - D H + OH - C OH - - D H + Na + C Na + - D H + Cl - C Cl - + u C H + ,

where the generalized-diffusion interpretation of the Nernst–Planck equation is applied in Eq. (23) and Dij is defined by Eq. (24). In the Nernst–Planck model, the inert species contribute to the transport processes in the form of the cross-coupling diffusivities, DH+Na+ and DH+Cl-, and the concentration gradients, CNa+ and CCl-. It has also been shown that the cross-coupling diffusivities (without reactions) can lead to hydrodynamic instabilities (Budroni2015). Such tight coupling between ionic species does not exist in the single-diffusivity model. Hence, the Nernst–Planck model is more expressive and capable of modeling the intricacies of reactive transport.

To further demonstrate the electromigration effects of the inert species, Fig. 11 shows the species concentration profile of simulations of CO2 dissolution in 0.01 M LiOH solution. We plot the 600 s snapshot, which is the onset time of the Nernst–Planck model (see Fig. 8). Figure 11 focuses on the concentration profile close to the top boundary, which is the elevation at 40 mm. When CO2 dissolves into the LiOH solution, the following chemical reactions occur


where H2CO3(aq) is a reaction intermediate that we have not considered in our numerical models. These reactions represent the entire process of CO2(aq) acidifying aqueous solutions. If we add Eqs. (41)–(43), this yields a simplified reaction,

(44) CO 2 ( aq ) + OH - HCO 3 - .

In Fig. 11, at elevations above 38 mm, excessive amounts of CO2(aq) react with OH, producing HCO3-. However, due to the electroneutrality constraints, the concentration of HCO3- is limited by the concentrations of the cations, particularly Li+. When CO2(aq) migrates deeper into the more basic regions, i.e., elevations between 36–38 mm, another reaction happens,

(45) HCO 3 - H + + CO 3 2 - .

In this basic region, HCO3- donates H+ and forms CO32-. We can still observe Li+ constraining the amount of CO3- owing to electroneutrality. Another cation that can contribute to electroneutrality is H+. However, in this basic region, the amount of H+ is negligible. Although Li+ is inert and does not take part in the reactions, Eqs. (41)–(45), it is still relevant in the overall reactive transport process. The transport of Li+ in the top left of Fig. 11 is evidence of its involvement in the reactive transport process.

Focusing on the density plots in Fig. 11, not only does the Nernst–Planck model produce lower densities than the single-diffusivity model, the density gradient between 37–38 mm is also higher in the Nernst–Planck model. This can explain the earlier onset time of the Nernst–Planck model compared to the onset time of the single-diffusivity model. We end the discussion with Fig. 12, which presents Li+ migrating as CO2 is dissolving and reacting in the solution. The tight coupling between the equilibrium reactions, electromigration effects, and density-driven flow elicits further investigation (both numerically and experimentally).

5 Conclusions

By comparing the single-diffusivity model and the Nernst–Planck model with reaction-driven flow experiments, we have demonstrated the importance and necessity under certain conditions of modeling the electromigration effects between charged species in aqueous environments using the Nernst–Planck model. Our results of simulating the reaction-driven flow experiments and convective dissolution of CO2 show that non-reacting species (Na+, Cl, Li+) influence the overall reactive transport process via electromigration effects, which couple the charged species. Such electromigration effects cannot be modeled using the single-diffusivity model, and further studies on quantifying the valid regimes of the single-diffusivity model are recommended. We make the following conclusions.

  • By comparing our numerical modeling results to previously published flow experiments, we have shown that the Nernst–Planck model is valid for modeling reactive transport processes. The processes in the experiments considered here are characterized by an intricate interplay between diffusion, reaction, electromigration, and density-driven convection.

  • Compared to the often-used single-diffusivity model, the Nernst–Planck model enables the numerical modeling of the electromigration of ionic species introduced by differing species diffusivities, resulting in more numerous and improved physical insights.

  • These reaction-driven flow experiments from literature can be further utilized for benchmarking reactive transport codes.

Code and data availability

The RetroPy code is available on GitHub (, last access: 10 August 2023), and version 1.0 is archived on Zenodo (; Huang2022). The RetroPy code is published under the GNU Lesser General Public License (LGPL). The code that produced the simulations is located in the example folder of RetroPy, where the chemical_convection folder contains the chemically driven convection examples and the CO2_convection folder contains the CO2 convective dissolution example. The data and scripts that recreate the figures are available on Zenodo (; Huang et al.2022a).

Video supplement

We provide the simulation results of Figs. 4, 6, and 12 in the mp4 format. They are located in (Huang et al.2022b).

Author contributions

PWH: Conceptualization, Methodology, Software, Formal analysis, Writing – Original Draft, Writing – Review & Editing, Visualization. BF: Writing – Review & Editing. CZQ: Writing – Review & Editing. MOS: Supervision, Writing – Review & Editing. AE: Conceptualization, Writing – Review & Editing, Project administration, Funding acquisition.

Competing interests

The contact author has declared that none of the authors has any competing interests.


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


We thank the Werner Siemens-Stiftung (Werner Siemens Foundation) for its support of the Geothermal Energy and Geofluids (, last access: 10 August 2023) Group at ETH Zurich, Switzerland.

We would like to thank our colleagues, Xiang-Zhao Kong and Isamu Naets, for all the helpful discussions. Furthermore, this work could not have been done without the visualization tool Matplotlib (Hunter2007). We utilized the color palettes provided by Harrower and Brewer (2003), Davis (2019), and Crameri (2021).

Financial support

This research has been supported by the Schweizerischer Nationalfonds zur Förderung der Wissenschaftlichen Forschung project entitled “Analysing spatial scaling effects in mineral reaction rates in porous media with a hybrid numerical model” (grant no. 175673).

Review statement

This paper was edited by Sylwester Arabas and reviewed by Lucjan Sapa and one anonymous referee.


Abd, A. S. and Abushaikha, A. S.: Reactive transport in porous media: a review of recent mathematical efforts in modeling geochemical reactions in petroleum subsurface reservoirs, SN Appl. Sci., 3, 401,, 2021. a

Aftab, A., Hassanpouryouzband, A., Xie, Q., Machuca, L. L., and Sarmadivaleh, M.: Toward a Fundamental Understanding of Geological Hydrogen Storage, Ind. Eng. Chem. Res., 61, 3233–3253,, 2022. a

Agartan, E., Trevisan, L., Cihan, A., Birkholzer, J., Zhou, Q., and Illangasekare, T. H.: Experimental study on effects of geologic heterogeneity in enhancing dissolution trapping of supercritical CO2, Water Resour. Res., 51, 1635–1648,, 2015. a

Åkerlöf, G. and Teare, J.: A Note on the Density of Aqueous Solutions of Hydrochloric Acid, J. Am. Chem. Soc., 60, 1226–1228,, 1938. a

Almarcha, C., Trevelyan, P. M. J., Grosfils, P., and De Wit, A.: Chemically Driven Hydrodynamic Instabilities, Phys. Rev. Lett., 104, 1–4,, 2010. a, b, c, d, e, f, g, h, i, j

Almarcha, C., R'Honi, Y., De Decker, Y., Trevelyan, P. M. J., Eckert, K., and De Wit, A.: Convective Mixing Induced by Acid–Base Reactions, J. Phys. Chem. B, 115, 9739–9744,, 2011. a

Alnæs, M. S.: UFL: a finite element form language, in: Automated Solution of Differential Equations by the Finite Element Method: The FEniCS Book, edited by Logg, A., Mardal, K.-A., and Wells, G. N., Springer Berlin Heidelberg, Berlin, Heidelberg, 303–338,, 2012. a

Alnæs, M. S., Blechta, J., Hake, J., Johansson, A., Kehlet, B., Logg, A., Richardson, C., Ring, J., Rognes, M. E., and Wells, G. N.: The FEniCS Project Version 1.5, Archive of Numerical Software, 3, 9–23,, 2015. a

Amarasinghe, W., Fjelde, I., Åge Rydland, J., and Guo, Y.: Effects of permeability on CO2 dissolution and convection at reservoir temperature and pressure conditions: A visualization study, Int. J. Greenh. Gas Con., 99, 103082,, 2020. a

Amestoy, P. R., Duff, I. S., L'Excellent, J.-Y., and Koster, J.: A Fully Asynchronous Multifrontal Solver Using Distributed Dynamic Scheduling, SIAM J. Matrix Anal. A., 23, 15–41,, 2001. a

Amestoy, P. R., Buttari, A., L'Excellent, J.-Y., and Mary, T.: Performance and Scalability of the Block Low-Rank Multifrontal Factorization on Multicore Architectures, ACM T. Math. Software, 45, 1–26,, 2019. a

Audigane, P., Gaus, I., Czernichowski-Lauriol, I., Pruess, K., and Xu, T.: Two-dimensional reactive transport modeling of CO2 injection in a saline aquifer at the Sleipner site, North Sea, Am. J. Sci., 307, 974–1008,, 2007. a

Avnir, D. and Kagan, M.: Spatial structures generated by chemical reactions at interfaces, Nature, 307, 717–720,, 1984. a

Azin, R., Raad, S. M. J., Osfouri, S., and Fatehi, R.: Onset of instability in CO2 sequestration into saline aquifer: scaling relationship and the effect of perturbed boundary, Heat Mass Transfer, 49, 1603–1612,, 2013. a

Babaei, M. and Islam, A.: Convective-Reactive CO2 Dissolution in Aquifers With Mass Transfer With Immobile Water, Water Resour. Res., 54, 9585–9604,, 2018. a

Bacuta, C.: A Unified Approach for Uzawa Algorithms, SIAM J. Numer. Anal., 44, 2633–2649,, 2006. a

Balay, S., Abhyankar, S., Adams, M. F., Benson, S., Brown, J., Brune, P., Buschelman, K., Constantinescu, E. M., Dalcin, L., Dener, A., Eijkhout, V., Gropp, W. D., Hapla, V., Isaac, T., Jolivet, P., Karpeev, D., Kaushik, D., Knepley, M. G., Kong, F., Kruger, S., May, D. A., McInnes, L. C., Mills, R. T., Mitchell, L., Munson, T., Roman, J. E., Rupp, K., Sanan, P., Sarich, J., Smith, B. F., Zampini, S., Zhang, H., Zhang, H., and Zhang, J.: PETSc Web page,
url (last access: 3 November 2022), 2022. a

Benzi, M., Golub, G. H., and Liesen, J.: Numerical solution of saddle point problems, Acta Numer., 14, 1–137,, 2005. a

Bethke, C. M.: Geochemical and Biogeochemical Reaction Modeling, Cambridge University Press, 2nd edn.,, 2007. a

Boffi, D., Brezzi, F., and Fortin, M.: Mixed Finite Element Methods and Applications, Springer, Berlin, Heidelberg,, 2013. a

Bordeaux-Rego, F., Sanaei, A., and Sepehrnoori, K.: Enhancement of Simulation CPU Time of Reactive-Transport Flow in Porous Media: Adaptive Tolerance and Mixing Zone-Based Approach, Transport Porous Med., 143, 127–150,, 2022. a

Boudreau, B. P., Meysman, F. J. R., and Middelburg, J. J.: Multicomponent ionic diffusion in porewaters: Coulombic effects revisited, Earth Planet. Sc. Lett., 22, 653–666,, 2004. a, b

Bożek, B., Sapa, L., and Danielewski, M.: Difference Methods to One and Multidimensional Interdiffusion Models with Vegard Rule, Math. Model. Anal., 24, 276–296,, 2019. a

Bożek, B., Sapa, L., Tkacz-Śmiech, K., Zajusz, M., and Danielewski, M.: Compendium About Multicomponent Interdiffusion in Two Dimensions, Metall. Mater. Trans. A, 52, 3221–3231,, 2021. a

Bratsun, D., Kostarev, K., Mizev, A., and Mosheva, E.: Concentration-dependent diffusion instability in reactive miscible fluids, Phys. Rev. E, 92, 011 003,, 2015. a

Bratsun, D., Mizev, A., Mosheva, E., and Kostarev, K.: Shock-wave-like structures induced by an exothermic neutralization reaction in miscible fluids, Phys. Rev. E, 96, 053 106,, 2017. a, b, c, d

Bratsun, D., Mizev, A., and Mosheva, E.: Extended classification of the buoyancy-driven flows induced by a neutralization reaction in miscible fluids. Part 2. Theoretical study, J. Fluid Mech., 916, A23,, 2021. a, b, c, d

Bratsun, D. A., Oschepkov, V. O., Mosheva, E. A., and Siraev, R. R.: The effect of concentration-dependent diffusion on double-diffusive instability, Phys. Fluids, 34, 034 112,, 2022. a, b

Brezzi, F., Douglad, J., and Marini, L. D.: Two families of mixed finite elements for second order elliptic problems, Numerische Mathematik, 47, 217–235,, 1985. a

Bringedal, C., Schollenberger, T., Pieters, G. J. M., van Duijn, C. J., and Helmig, R.: Evaporation-Driven Density Instabilities in Saturated Porous Media, Transport Porous Med., 143, 297–341,, 2022. a

Brouzet, C., Méheust, Y., and Meunier, P.: CO2 convective dissolution in a three-dimensional granular porous medium: An experimental study, Phys. Rev. Fluids, 7, 033 802,, 2022. a

Budroni, M. A.: Cross-diffusion-driven hydrodynamic instabilities in a double-layer system: General classification and nonlinear simulations, Phys. Rev. E, 92, 063007,, 2015. a

Cappellen, P. V. and Gaillard, J.-F.: Chapter 8. Biogeochemical Dynamics in Aquatic Sediments, in: Reactive Transport in Porous Media, edited by Lichtner, P. C., Steefel, C. I., and Oelkers, E. H., De Gruyter, 335–376,, 1996. a

Carrayrou, J., Mosé, R., and Behra, P.: Operator-splitting procedures for reactive transport and comparison of mass balance errors, J. Contam. Hydrol., 68, 239–268,, 2004. a

Carrera, J., Saaltink, M. W., Soler-Sagarra, J., Wang, J., and Valhondo, C.: Reactive Transport: A Review of Basic Concepts with Emphasis on Biochemical Processes, Energies, 15,, 2022. a

Chang, E., Brewer, A. W., Park, D. M., Jiao, Y., and Lammers, L. N.: Selective Biosorption of Valuable Rare Earth Elements Among Co-Occurring Lanthanides, Environ. Eng. Sci., 38, 154–164,, 2021. a

Cheng, C. and Milsch, H.: Permeability Variations in Illite-Bearing Sandstone: Effects of Temperature and NaCl Fluid Salinity, J. Geophys. Res.-Sol. Ea., 125, e2020JB020122,, 2020. a

Cherezov, I., Cardoso, S. S., and Kim, M. C.: Acceleration of convective dissolution by an instantaneous chemical reaction: A comparison of experimental and numerical results, Chem. Eng. Sci., 181, 298–310,, 2018. a

Citri, O., Kagan, M. L., Kosloff, R., and Avnir, D.: Evolution of chemically induced unstable density gradients near horizontal reactive interfaces, Langmuir, 6, 559–564,, 1990. a

Class, H., Weishaupt, K., and Trötschler, O.: Experimental and Simulation Study on Validating a Numerical Model for CO2 Density-Driven Dissolution in Water, Water, 12, 738,, 2020. a

Cochepin, B., Trotignon, L., Bildstein, O., Steefel, C., Lagneau, V., and Van der lee, J.: Approaches to modelling coupled flow and reaction in a 2D cementation experiment, Adv. Water Resour., 31, 1540–1551,, 2008. a

Cogorno, J., Stolze, L., Muniruzzaman, M., and Rolle, M.: Dimensionality effects on multicomponent ionic transport and surface complexation in porous media, Geochim. Cosmochim. Ac., 318, 230–246,, 2022. a

Connolly, J. A. D.: The geodynamic equation of state: What and how, Geochem. Geophy. Geosy., 10, Q10014,, 2009. a

Constantin, P., Ignatova, M., and Lee, F.-N.: Interior Electroneutrality in Nernst–Planck–Navier–Stokes Systems, Arch. Ration. Mech. An., 242, 1091–1118,, 2021. a

Constantin, P., Ignatova, M., and Lee, F.-N.: Existence and stability of nonequilibrium steady states of Nernst–Planck–Navier–Stokes systems, Physica D, 442, 133536,, 2022. a

Crameri, F.: Scientific colour maps, Zenodo [data set],, 2021. a

Damiani, L. H., Kosakowski, G., Glaus, M. A., and Churakov, S. V.: A framework for reactive transport modeling using FEniCS–Reaktoro: governing equations and benchmarking results, Comput. Geosci., 24, 1071–1085,, 2020. a, b

Davis, M.: Palettable: Color palettes for Python, (last access: 25 November 2022), 2019. a

De Lucia, M. and Kühn, M.: DecTree v1.0 – chemistry speedup in reactive transport simulations: purely data-driven and physics-based surrogates, Geosci. Model Dev., 14, 4713–4730,, 2021. a

De Lucia, M., Kühn, M., Lindemann, A., Lübke, M., and Schnor, B.: POET (v0.1): speedup of many-core parallel reactive transport simulations with fast DHT lookups, Geosci. Model Dev., 14, 7391–7409,, 2021. a

De Wit, A.: Chemo-hydrodynamic patterns in porous media, Philos. T. Roy. Soc. A, 374, 20150419,, 2016. a

De Wit, A.: Chemo-Hydrodynamic Patterns and Instabilities, Annu. Rev. Fluid Mech., 52, 531–555,, 2020. a

Dickinson, E. J. F., Limon-Peterson, J. G., and Compton, R. G.: The electroneutrality assumption in electrochemistry, J. Solid State Electr., 15, 1335–1345,, 2011. a

Donea, J. and Huerta, A.: Finite Element Methods for Flow Problems, John Wiley & Sons,, 2003. a

Dreyer, W., Guhlke, C., and Müller, R.: Overcoming the shortcomings of the Nernst–Planck model, Phys. Chem. Chem. Phys., 15, 7075–7086,, 2013. a

Drummond, S.: Boiling and mixing of hydrothermal fluids: chemical effects on mineral precipitation, PhD thesis, Pennsylvania State University, (last access: 21 November 2022), 1981. a

Eckert, K. and Grahn, A.: Plume and Finger Regimes Driven by an Exothermic Interfacial Reaction, Phys. Rev. Lett., 82, 4436–4439,, 1999. a

Eckert, K., Acker, M., and Shi, Y.: Chemical pattern formation driven by a neutralization reaction. I. Mechanism and basic features, Phys. Fluids, 16, 385–399,, 2004. a

Ezekiel, J., Adams, B. M., Saar, M. O., and Ebigbo, A.: Numerical analysis and optimization of the performance of CO2-Plume Geothermal (CPG) production wells and implications for electric power generation, Geothermics, 98, 102270,, 2022. a

Filipek, R., Kalita, P., Sapa, L., and Szyszkiewicz, K.: On local weak solutions to Nernst–Planck–Poisson system, Appl. Anal., 96, 2316–2332,, 2017. a

Flavell, A., Machen, M., Eisenberg, B., Kabre, J., Liu, C., and Li, X.: A conservative finite difference scheme for Poisson–Nernst–Planck equations, J. Comput. Phys., 13, 235–249,, 2014. a

Fleming, M. R., Adams, B. M., Kuehn, T. H., Bielicki, J. M., and Saar, M. O.: Increased Power Generation due to Exothermic Water Exsolution in CO2 Plume Geothermal (CPG) Power Plants, Geothermics, 88, 101865,, 2020. a

Fortin, M. and Glowinski, R.: Augmented Lagrangian Methods in Quadratic Programming, in: Augmented Lagrangian Methods: Applications to the Numerical Solution of Boundary-Value Problems, edited by Fortin, M. and Glowinski, R., vol. 15 of Studies in Mathematics and Its Applications, chap. 1, 1–46, Elsevier,, 1983. a, b

Frank, F., Ray, N., and Knabner, P.: Numerical investigation of homogenized Stokes–Nernst–Planck–Poisson systems, Computing and Visualization in Science, 14, 385–400,, 2011. a

Frizon, F., Lorente, S., Ollivier, J., and Thouvenot, P.: Transport model for the nuclear decontamination of cementitious materials, Comp. Mater. Sci., 27, 507–516,, 2003. a

Fu, B., Zhang, R., Liu, J., Cui, L., Zhu, X., and Hao, D.: Simulation of CO2 Rayleigh Convection in Aqueous Solutions of NaCl, KCl, MgCl2 and CaCl2 using Lattice Boltzmann Method, Int. J. Greenh. Gas Con., 98, 103066,, 2020. a

Fu, B., Zhang, R., Xiao, R., Cui, L., Liu, J., Zhu, X., and Hao, D.: Simulation of interfacial mass transfer process accompanied by Rayleigh convection in NaCl solution, Int. J. Greenh. Gas Con., 106, 103281,, 2021. a

Gamazo, P., Slooten, L. J., Carrera, J., Saaltink, M. W., Bea, S., and Soler, J.: PROOST: object-oriented approach to multiphase reactive transport modeling in porous media, J. Hydroinform., 18, 310–328,, 2015. a

Giambalvo, E. R., Steefel, C. I., Fisher, A. T., Rosenberg, N. D., and Wheat, C. G.: Effect of fluid-sediment reaction on hydrothermal fluxes of major elements, eastern flank of the Juan de Fuca Ridge, Geochim. Cosmochim. Ac., 66, 1739–1757,, 2002. a

Gimmi, T. and Alt-Epping, P.: Simulating Donnan equilibria based on the Nernst-Planck equation, Geochim. Cosmochim. Ac., 232, 1–13,, 2018. a

Glaus, M. A., Birgersson, M., Karnland, O., and Van Loon, L. R.: Seeming Steady-State Uphill Diffusion of 22Na+ in Compacted Montmorillonite, Environ. Sci. Technol., 47, 11522–11527,, 2013. a

Golparvar, A., Kästner, M., and Thullner, M.: P3D-BRNS v1.0.0: A Three-dimensional, Multiphase, Multicomponent, Pore-scale Reactive Transport Modelling Package for Simulating Biogeochemical Processes in Subsurface Environments, Geosci. Model Dev. Discuss. [preprint],, in review, 2022. a

Grimm Lima, M., Schädle, P., Green, C. P., Vogler, D., Saar, M. O., and Kong, X.-Z.: Permeability Impairment and Salt Precipitation Patterns During CO2 Injection Into Single Natural Brine-Filled Fractures, Water Resour. Res., 56, e2020WR027213,, 2020. a

Guo, R., Sun, H., Zhao, Q., Li, Z., Liu, Y., and Chen, C.: A Novel Experimental Study on Density-Driven Instability and Convective Dissolution in Porous Media, Geophys. Res. Lett., 48, e2021GL095619,, 2021. a

Harrower, M. and Brewer, C. A.: An Online Tool for Selecting Colour Schemes for Maps, The Cartographic Journal, 40, 27–37,, 2003. a

Hassanpouryouzband, A., Adie, K., Cowen, T., Thaysen, E. M., Heinemann, N., Butler, I. B., Wilkinson, M., and Edlmann, K.: Geological Hydrogen Storage: Geochemical Reactivity of Hydrogen with Sandstone Reservoirs, ACS Energy Letters, 7, 2203–2210,, 2022. a

Hejazi, S. H. and Azaiez, J.: Nonlinear simulation of transverse flow interactions with chemically driven convective mixing in porous media, Water Resour. Res., 49, 4607–4618,, 2013. a

Helgeson, H. C. and Kirkham, D. H.: Theoretical prediction of the thermodynamic behavior of aqueous electrolytes at high pressures and temperatures; I, Summary of the thermodynamic/electrostatic properties of the solvent, Am. J. Sci., 274, 1089–1198,, 1974a. a

Helgeson, H. C. and Kirkham, D. H.: Theoretical prediction of the thermodynamic behavior of aqueous electrolytes at high pressures and temperatures; II, Debye-Huckel parameters for activity coefficients and relative partial molal properties, Am. J. Sci., 274, 1199–1261,, 1974b. a

Helgeson, H. C. and Kirkham, D. H.: Theoretical prediction of the thermodynamic properties of aqueous electrolytes at high pressures and temperatures. III. Equation of state for aqueous species at infinite dilution, Am. J. Sci., 276, 97–240,, 1976. a

Helgeson, H. C., Kirkham, D. H., and Flowers, G. C.: Theoretical prediction of the thermodynamic behavior of aqueous electrolytes by high pressures and temperatures; IV, Calculation of activity coefficients, osmotic coefficients, and apparent molal and standard and relative partial molal properties to 600 C and 5 KB, Am. J. Sci., 281, 1249–1516,, 1981. a

Herz, M., Ray, N., and Knabner, P.: Existence and uniqueness of a global weak solution of a Darcy-Nernst-Planck-Poisson system, GAMM-Mitteilungen, 35, 191–208,, 2012. a

Hingerl, F. F.: Geothermal electrolyte solutions: thermodynamic model and computational fitting framework development, PhD thesis, ETH Zurich,, 2012. a

Huang, P.-W.: pwhuang/RetroPy: First release of RetroPy (v1.0), Zenodo [code],, 2022. a

Huang, Y., Shao, H., Wieland, E., Kolditz, O., and Kosakowski, G.: A new approach to coupled two-phase reactive transport simulation for long-term degradation of concrete, Constr. Build. Mater., 190, 805–829,, 2018. a

Huang, P.-W., Flemisch, B., Qin, C.-Z., Saar, M. O., and Ebigbo, A.: Data supplement for: Validating the Nernst–Planck transport model under reaction-driven flow conditions using RetroPy v1.0 (v1.0.0), Zenodo [data set],, 2022a. a

Huang, P.-W., Flemisch, B., Qin, C.-Z., Saar, M. O., and Ebigbo, A.: Video supplement for: Validating the Nernst–Planck transport model under reaction-driven flow conditions using RetroPy v1.0, ETH Zürich [video],, 2022b. a

Hunter, J. D.: Matplotlib: A 2D graphics environment, Comput. Sci. Eng., 9, 90–95,, 2007. a

Ignatova, M. and Shu, J.: Global Smooth Solutions of the Nernst–Planck–Darcy System, J. Math. Fluid Mech., 24, 26,, 2022. a

Illés, B., Medgyes, B., Dušek, K., Bušek, D., Skwarek, A., and Géczy, A.: Numerical simulation of electrochemical migration of Cu based on the Nernst-Plank equation, Int. J. Heat Mass Tran., 184, 122268,, 2022. a

Izumoto, S., Huisman, J. A., Zimmermann, E., Heyman, J., Gomez, F., Tabuteau, H., Laniel, R., Vereecken, H., Méheust, Y., and Le Borgne, T.: Pore-Scale Mechanisms for Spectral Induced Polarization of Calcite Precipitation Inferred from Geo-Electrical Millifluidics, Environ. Sci. Technol., 56, 4998–5008,, 2022. a

Jasielec, J. J.: Electrodiffusion Phenomena in Neuroscience and the Nernst–Planck–Poisson Equations, Electrochem., 2, 197–215,, 2021. a

Jiang, L., Wang, S., Liu, D., Zhang, W., Lu, G., Liu, Y., and Zhao, J.: Change in Convection Mixing Properties with Salinity and Temperature: CO2 Storage Application, Polymers, 12, 2084,, 2020. a

Johnson, J. W., Oelkers, E. H., and Helgeson, H. C.: SUPCRT92: A software package for calculating the standard molal thermodynamic properties of minerals, gases, aqueous species, and reactions from 1 to 5000 bar and 0 to 1000 C, Comput. Geosci., 18, 899–947,, 1992. a

Jotkar, M., De Wit, A., and Rongy, L.: Control of chemically driven convective dissolution by differential diffusion effects, Phys. Rev. Fluids, 6, 053504,, 2021. a

Jyoti, A. and Haese, R. R.: Validation of a multicomponent reactive-transport model at pore scale based on the coupling of COMSOL and PhreeqC, Comput. Geosci., 156, 104870,, 2021. a

Kadeethum, T., Lee, S., Ballarin, F., Choo, J., and Nick, H.: A locally conservative mixed finite element framework for coupled hydro-mechanical–chemical processes in heterogeneous porous media, Comput. Geosci., 152, 104774,, 2021. a

Kim, M. C.: Effect of the irreversible A+BC reaction on the onset and the growth of the buoyancy-driven instability in a porous medium: Asymptotic, linear, and nonlinear stability analyses, Phys. Rev. Fluids, 4, 073901,, 2019. a

Kirby, B. J.: Species and Charge Transport, in: Micro- and Nanoscale Fluid Mechanics: Transport in Microfluidic Devices, Cambridge University Press, 250–264,, 2010. a, b

Kneafsey, T. J. and Pruess, K.: Laboratory Flow Experiments for Visualizing Carbon Dioxide-Induced, Density-Driven Brine Convection, Transport Porous Med., 82, 123–139,, 2010. a

Kontturi, K., Murtomäki, L., andManzanares, A. J.: Ionic Transport Processes: in Electrochemistry and Membrane Science, Oxford University Press,, 2008. a

Kovtunenko, V. A. and Zubkova, A. V.: Existence and two-scale convergence of the generalised Poisson–Nernst–Planck problem with non-linear interface conditions, Eur. J. Appl. Math., 32, 683–710,, 2021. a

Kulik, D. A., Wagner, T., Dmytrieva, S. V., Kosakowski, G., Hingerl, F. F., Chudnenko, K. V., and Berner, U. R.: Theoretical prediction of the thermodynamic properties of aqueous electrolytes at high pressures and temperatures. III. Equation of state for aqueous species at infinite dilution, Computat. Geosci., 17, 1–24,, 2012. a

Kwon, O., Herbert, B. E., and Kronenberg, A. K.: Permeability of illite-bearing shale: 2. Influence of fluid chemistry on flow and functionally connected pores, J. Geophys. Res.-Sol. Ea., 109, B10206,, 2004a. a

Kwon, O., Kronenberg, A. K., Gangi, A. F., Johnson, B., and Herbert, B. E.: Permeability of illite-bearing shale: 1. Anisotropy and effects of clay content and loading, J. Geophys. Res.-Sol. Ea., 109, B10205,, 2004b. a

Kyas, S., Volpatto, D., Saar, M. O., and Leal, A. M. M.: Accelerated reactive transport simulations in heterogeneous porous media using Reaktoro and Firedrake, Computat. Geosci., 26, 295–327,, 2022. a, b

Laloy, E. and Jacques, D.: Speeding Up Reactive Transport Simulations in Cement Systems by Surrogate Geochemical Modeling: Deep Neural Networks and k-Nearest Neighbors, Transport Porous Med., 143, 433–462,, 2022. a

Leal, A. M., Blunt, M. J., and LaForce, T. C.: Efficient chemical equilibrium calculations for geochemical speciation and reactive transport modelling, Geochim. Cosmochim. Ac., 131, 301–322,, 2014. a

Leal, A. M., Kulik, D. A., and Kosakowski, G.: Computational methods for reactive transport modeling: A Gibbs energy minimization approach for multiphase equilibrium calculations, Adv. Water Resour., 88, 231–240,, 2016. a

Leal, A. M. M.: Reaktoro for Python and C++, (last access: 21 November 2022), 2022. a

Leal, A. M. M., Kulik, D. A., Smith, W. R., and Saar, M. O.: An overview of computational methods for chemical equilibrium and kinetic calculations for geochemical and reactive transport modeling, Pure Appl. Chem., 89, 597–643,, 2017. a

Leal, A. M. M., Kyas, S., Kulik, D. A., and Saar, M. O.: Accelerating Reactive Transport Modeling: On-Demand Machine Learning Algorithm for Chemical Equilibrium Calculations, Transport Porous Med., 133, 161–204,, 2020. a

Lee, F.-N.: Global Regularity for Nernst–Planck–Navier–Stokes Systems with mixed boundary conditions, Nonlinearity, 36, 255–286,, 2022. a

Lees, E., Rokkam, S., Shanbhag, S., and Gunzburger, M.: The electroneutrality constraint in nonlocal models, J. Chem. Phys., 147, 124102,, 2017. a

Lemaigre, L., Budroni, M. A., Riolfo, L. A., Grosfils, P., and De Wit, A.: Asymmetric Rayleigh-Taylor and double-diffusive fingers in reactive systems, Phys. Fluids, 25, 014103,, 2013. a, b, c, d

Lemmon, E. W. and Harvey, A. H.: Thermophysical Properties of Water and Steam, in: CRC Handbook of Chemistry and Physics, edited by: Rumble, J. R., CRC Press/Taylor & Francis, Boca Raton, FL, 102nd edn., 2021. a

Lichtner, P. C.: Continuum model for simultaneous chemical reactions and mass transport in hydrothermal systems, Geochim. Cosmochim. Ac., 49, 779–800,, 1985. a, b

Lichtner, P. C.: Principles and Practice of Reactive Transport Modeling, MRS Proceedings, 353, 117–130,, 1994. a

Ling, F. T., Plattenberger, D. A., Peters, C. A., and Clarens, A. F.: Sealing Porous Media through Calcium Silicate Reactions with CO2 to Enhance the Security of Geologic Carbon Sequestration, Environ. Eng. Sci., 38, 127–142,, 2021. a

Liu, H. and Maimaitiyiming, W.: Efficient, Positive, and Energy Stable Schemes for Multi-D Poisson–Nernst–Planck Systems, J. Sci. Comput., 148, 92,, 2021. a

Liu, H., Wang, Z., Yin, P., and Yu, H.: Positivity-preserving third order DG schemes for Poisson–Nernst–Planck equations, J. Comput. Phys., 452, 110777,, 2022. a

López-Vizcaíno, R., Cabrera, V., Sprocati, R., Muniruzzaman, M., Rolle, M., Navarro, V., and Yustres, Á.: A modeling approach for electrokinetic transport in double-porosity media, Electrochim. Acta, 431, 141139,, 2022. a

Luhmann, A. J., Kong, X.-Z., Tutolo, B. M., Ding, K., Saar, M. O., and Seyfried, W. E. J.: Permeability Reduction Produced by Grain Reorganization and Accumulation of Exsolved CO2 during Geologic Carbon Sequestration: A New CO2 Trapping Mechanism, Environ. Sci. Technol., 47, 242–251,, 2013. a

Maes, J. and Menke, H. P.: GeoChemFoam: Direct Modelling of Multiphase Reactive Transport in Real Pore Geometries with Equilibrium Reactions, Transport Porous Med., 139, 271–299,, 2021. a

Mahmoodpour, S., Rostami, B., Soltanian, M. R., and Amooie, M. A.: Effect of brine composition on the onset of convection during CO2 dissolution in brine, Comput. Geosci., 124, 1–13,, 2019. a

Mahmoodpour, S., Amooie, M. A., Rostami, B., and Bahrami, F.: Effect of gas impurity on the convective dissolution of CO2 in porous media, Energy, 199, 117397,, 2020. a

Marini, L.: Chapter 4 – The Aqueous Electrolyte Solution, in: Geological Sequestration of Carbon Dioxide, vol. 11 of Developments in Geochemistry, Elsevier, 53–77,, 2007. a

Martinez, M. J. and Hesse, M. A.: Two-phase convective CO2 dissolution in saline aquifers, Water Resour. Res., 52, 585–599,, 2016. a

Miron, G. D., Leal, A. M. M., and Yapparova, A.: Thermodynamic Properties of Aqueous Species Calculated Using the HKF Model: How Do Different Thermodynamic and Electrostatic Models for Solvent Water Affect Calculated Aqueous Properties?, Geofluids, 2019, 5750390,, 2019. a

Mizev, A., Mosheva, E., and Bratsun, D.: Extended classification of the buoyancy-driven flows induced by a neutralization reaction in miscible fluids. Part 1. Experimental study, J. Fluid Mech., 916, A22,, 2021. a, b, c, d, e, f

Moortgat, J., Li, M., Amooie, M. A., and Zhu, D.: A higher-order finite element reactive transport model for unstructured and fractured grids, Sci. Rep.-UK, 10, 15572,, 2020. a

Muniruzzaman, M. and Rolle, M.: Impact of multicomponent ionic transport on pH fronts propagation in saturated porous media, Water Resour. Res., 51, 6739–6755,, 2015. a

Muniruzzaman, M. and Rolle, M.: Experimental investigation of the impact of compound-specific dispersion and electrostatic interactions on transient transport and solute breakthrough, Water Resour. Res., 53, 1189–1209,, 2017. a

Muniruzzaman, M. and Rolle, M.: Multicomponent Ionic Transport Modeling in Physically and Electrostatically Heterogeneous Porous Media With PhreeqcRM Coupling for Geochemical Reactions, Water Resour. Res., 55, 11121–11143,, 2019. a

Muniruzzaman, M., Haberer, C. M., Grathwohl, P., and Rolle, M.: Multicomponent ionic dispersion during transport of electrolytes in heterogeneous porous media: Experiments and model-based interpretation, Geochim. Cosmochim. Ac., 141, 656–669,, 2014. a

Naumov, D. Y., Bilke, L., Fischer, T., Rink, K., Wang, W., Watanabe, N., Lu, R., Grunwald, N., Zill, F., Buchwald, J., Huang, Y., Bathmann, J., Chen, C., Chen, S., Meng, B., Shao, H., Kern, D., Yoshioka, K., Garibay Rodriguez, J., Miao, X., Parisio, F., Silbermann, C., Thiedau, J., Walther, M., Kaiser, S., Boog, J., Zheng, T., Meisel, T., and Ning, Z.: OpenGeoSys, Zenodo [data set],, 2022. a

Neufeld, J. A., Hesse, M. A., Riaz, A., Hallworth, M. A., Tchelepi, H. A., and Huppert, H. E.: Convective dissolution of carbon dioxide in saline aquifers, Geophys. Res. Lett., 37, L22404,, 2010. a

Nishikata, E., Ishii, T., and Ohta, T.: Viscosities of aqueous hydrochloric acid solutions, and densities and viscosities of aqueous hydroiodic acid solutions, J. Chem. Eng. Data, 26, 254–256,, 1981. a

Oliveira, T. D. S., Blunt, M. J., and Bijeljic, B.: Modelling of multispecies reactive transport on pore-space images, Adv. Water Resour., 127, 192–208,, 2019. a

Onsager, L.: Theories and Problems of Liquid Diffusion, Ann. NY. Acad. Sci., 46, 241–265,, 1945. a

Pamukcu, S.: Electrochemical Transport and Transformations, in: Electrochemical Remediation Technologies for Polluted Soils, Sediments and Groundwater, edited by: Reddy, K. R. and Cameselle, C., John Wiley & Sons, Ltd, chap. 2, 29–64,, 2009. a

Parkhurst, D. L. and Appelo, C. A. J.: Description of input and examples for PHREEQC version 3–A computer program for speciation, batch-reaction, one-dimensional transport, and inverse geochemical calculations, (last access: 10 August 2023), 2013. a

Permann, C. J., Gaston, D. R., Andrš, D., Carlsen, R. W., Kong, F., Lindsay, A. D., Miller, J. M., Peterson, J. W., Slaughter, A. E., Stogner, R. H., and Martineau, R. C.: MOOSE: Enabling massively parallel multiphysics simulation, SoftwareX, 11, 100430,, 2020. a

Petersen, S. and Hack, K.: The thermochemistry library ChemApp and its applications, Int. J. Mater. Res., 98, 935–945,, 2007. a

Pogge von Strandmann, P. A. E., Burton, K. W., Snæbjörnsdóttir, S. O., Sigfússon, B., Aradóttir, E. S., Gunnarsson, I., Alfredsson, H. A., Mesfin, K. G., Oelkers, E. H., and Gislason, S. R.: Rapid CO2 mineralisation into calcite at the CarbFix storage site quantified using calcium isotopes, Nat. Commun., 10, 1983,, 2019. a

Priya, P., Kuhlman, K. L., and Aluru, N. R.: Pore-Scale Modeling of Electrokinetics in Geomaterials, Transport Porous Med., 137, 651–666,, 2021. a

Randolph, J. B. and Saar, M. O.: Combining geothermal energy capture with geologic carbon dioxide sequestration, Geophys. Res. Lett., 38, L10401,, 2011. a

Rasouli, P., Steefel, C. I., Mayer, K. U., and Rolle, M.: Benchmarks for multicomponent diffusion and electrochemical migration, Computat. Geosci., 19, 523–533,, 2015. a

Rathgeber, F., Ham, D. A., Mitchell, L., Lange, M., Luporini, F., Mcrae, A. T. T., Bercea, G.-T., Markall, G. R., and Kelly, P. H. J.: Firedrake: Automating the Finite Element Method by Composing Abstractions, ACM Trans. Math. Softw., 43, 24,, 2016. a

Raviart, P. A. and Thomas, J. M.: A mixed finite element method for 2-nd order elliptic problems, in: Mathematical Aspects of Finite Element Methods, edited by Galligani, I. and Magenes, E., Springer Berlin Heidelberg, Berlin, Heidelberg, 292–315,, 1977. a

Ray, N., Muntean, A., and Knabner, P.: Rigorous homogenization of a Stokes–Nernst–Planck–Poisson system, J. Math. Anal. Appl., 390, 374–393,, 2012a. a

Ray, N., van Noorden, T., Frank, F., and Knabner, P.: Multiscale Modeling of Colloid and Fluid Dynamics in Porous Media Including an Evolving Microstructure, Transport Porous Med., 95, 669–696,, 2012b. a

Reddy, K. R. and Cameselle, C.: Electrochemical Remediation Technologies for Polluted Soils, Sediments and Groundwater, John Wiley & Sons,, 2009. a

Reed, M. H.: Calculation of Simultaneous Chemical Equilibria in Aqueous-Mineral-Gas Systems and its Application to Modeling Hydrothermal Processes, in: Techniques in Hydrothermal Ore Deposits Geology, Society of Economic Geologists,, 1998. a

Revil, A. and Leroy, P.: Constitutive equations for ionic transport in porous shales, J. Geophys. Res.-Sol. Ea., 109, B03208,, 2004. a

Rivera, F. F., Pérez, T., Castañeda, L. F., and Nava, J. L.: Mathematical modeling and simulation of electrochemical reactors: A critical review, Chem. Eng. Sci., 239, 116622,, 2021. a

Rolle, M., Muniruzzaman, M., Haberer, C. M., and Grathwohl, P.: Coulombic effects in advection-dominated transport of electrolytes in porous media: Multicomponent ionic dispersion, Geochim. Cosmochim. Ac., 120, 195–205,, 2013. a

Rolle, M., Sprocati, R., Masi, M., Jin, B., and Muniruzzaman, M.: Nernst-Planck-based Description of Transport, Coulombic Interactions, and Geochemical Reactions in Porous Media: Modeling Approach and Benchmark Experiments, Water Resour. Res., 54, 3176–3195,, 2018. a

Rolle, M., Albrecht, M., and Sprocati, R.: Impact of solute charge and diffusion coefficient on electromigration and mixing in porous media, J. Contam. Hydrol., 244, 103933,, 2022. a

Rubinstein, I.: Locally Electro-Neutral Electro-Diffusion without Electric Current, in: Electro-Diffusion of Ions, Society for Industrial and Applied Mathematics, chap. 3, 59–103,, 1990. a

Samson, E. and Marchand, J.: Modeling the transport of ions in unsaturated cement-based materials, Comput. Struct., 85, 1740–1756,, 2007. a

Sapa, L., Bożek, B., Tkacz-Śmiech, K., Zajusz, M., and Danielewski, M.: Interdiffusion in many dimensions: mathematical models, numerical simulations and experiment, Math. Mech. Solids, 25, 2178–2198,, 2020. a

Savino, M., Lévy-Leduc, C., Leconte, M., and Cochepin, B.: An active learning approach for improving the performance of equilibrium based chemical simulations, Computat. Geosci., 26, 365–380,, 2022. a

Shafabakhsh, P., Ataie-Ashtiani, B., Simmons, C. T., Younes, A., and Fahs, M.: Convective-reactive transport of dissolved CO2 in fractured-geological formations, Int. J. Greenh. Gas Con., 109, 103365,, 2021. a

Shao, H., Dmytrieva, S. V., Kolditz, O., Kulik, D. A., Pfingsten, W., and Kosakowski, G.: Modeling reactive transport in non-ideal aqueous–solid solution system, Appl. Geochem., 24, 1287–1300,, 2009. a

Shen, J. and Xu, J.: Unconditionally positivity preserving and energy dissipative schemes for Poisson–Nernst–Planck equations, Numerische Mathematik, 148, 671–697,, 2021. a

Sin, I. and Corvisier, J.: Multiphase Multicomponent Reactive Transport and Flow Modeling, Rev. Mineral. Geochem., 85, 143–195, (last access: 28 November 2022), 2019. a

Singh, M., Chaudhuri, A., Chu, S. P., Stauffer, P. H., and Pawar, R. J.: Analysis of evolving capillary transition, gravitational fingering, and dissolution trapping of CO2 in deep saline aquifers during continuous injection of supercritical CO2, Int. J. Greenh. Gas Con., 82, 281–297,, 2019. a

Sipos, P. M., Hefter, G., and May, P. M.: Viscosities and Densities of Highly Concentrated Aqueous MOH Solutions (M+ = Na+, K+, Li+, Cs+, (CH3)4N+) at 25.0 C, J. Chem. Eng. Data, 45, 613–617,, 2000. a

Soltanian, M. R., Hajirezaie, S., Hosseini, S. A., Dashtian, H., Amooie, M. A., Meyal, A., Ershadnia, R., Ampomah, W., Islam, A., and Zhang, X.: Multicomponent reactive transport of carbon dioxide in fluvial heterogeneous aquifers, J. Nat. Gas Sci. Eng., 65, 212–223,, 2019. a

Song, Z., Cao, X., and Huang, H.: Electroneutral models for dynamic Poisson-Nernst-Planck systems, Phys. Rev. E, 97, 012411,, 2018. a

Sprocati, R. and Rolle, M.: Integrating Process-Based Reactive Transport Modeling and Machine Learning for Electrokinetic Remediation of Contaminated Groundwater, Water Resour. Res., 57, e2021WR029959,, 2021. a

Sprocati, R. and Rolle, M.: On the interplay between electromigration and electroosmosis during electrokinetic transport in heterogeneous porous media, Water Res., 213, 118161,, 2022. a

Sprocati, R., Masi, M., Muniruzzaman, M., and Rolle, M.: Modeling electrokinetic transport and biogeochemical reactions in porous media: A multidimensional Nernst–Planck–Poisson approach with PHREEQC coupling, Adv. Water Resour., 127, 134–147,, 2019. a

Steefel, C. I. and MacQuarrie, K. T. B.: Chapter 2. Approaches to modeling of reactive transport in porous media, in: Reactive Transport in Porous Media, edited by Lichtner, P. C., Steefel, C. I., and Oelkers, E. H., De Gruyter, 83–130,, 2018. a

Tabrizinejadas, S., Carrayrou, J., Saaltink, M. W., Baalousha, H. M., and Fahs, M.: On the Validity of the Null Current Assumption for Modeling Sorptive Reactive Transport and Electro-Diffusion in Porous Media, Water, 13, 2221,, 2021. a

Teng, Y., Wang, P., Jiang, L., Liu, Y., and Wei, Y.: New Spectrophotometric Method for Quantitative Characterization of Density-Driven Convective Instability, Polymers, 13, 661,, 2021. a

Thieulot, C. and Bangerth, W.: On the choice of finite element for applications in geodynamics, Solid Earth, 13, 229–249,, 2022. a

Thomas, C., Loodts, V., Rongy, L., and De Wit, A.: Convective dissolution of CO2 in reactive alkaline solutions: Active role of spectator ions, Int. J. Greenh. Gas Con., 53, 230–242,, 2016. a, b, c, d, e, f, g, h, i, j

Thomas, C., Dehaeck, S., and De Wit, A.: Convective dissolution of CO2 in water and salt solutions, Int. J. Greenh. Gas Con., 72, 105–116,, 2018. a

Tournassat, C. and Steefel, C. I.: Modeling diffusion processes in the presence of a diffuse layer at charged mineral surfaces: a benchmark exercise, Computat. Geosci., 25, 1319–1336,, 2021. a

Tournassat, C., Steefel, C. I., and Gimmi, T.: Solving the Nernst-Planck Equation in Heterogeneous Porous Media With Finite Volume Methods: Averaging Approaches at Interfaces, Water Resour. Res., 56, e2019WR026832,, 2020. a, b, c

Trevelyan, P. M. J., Almarcha, C., and De Wit, A.: Buoyancy-driven instabilities of miscible two-layer stratifications in porous media and Hele-Shaw cells, J. Fluid Mech., 670, 38–65,, 2011. a

Tsinober, A., Rosenzweig, R., Class, H., Helmig, R., and Shavit, U.: The Role of Mixed Convection and Hydrodynamic Dispersion During CO2 Dissolution in Saline Aquifers: A Numerical Study, Water Resour. Res., 58, e2021WR030494,, 2022. a

Tutolo, B. M., Luhmann, A. J., Kong, X.-Z., Saar, M. O., and Seyfried, W. E. J.: Experimental Observation of Permeability Changes In Dolomite at CO2 Sequestration Conditions, Environ. Sci. Technol., 48, 2445–2452,, 2014. a

Uzawa, H.: Iterative methods for concave programming, in: Studies in linear and nonlinear progamming, edited by: Arrow, K. J., Hurwicz, L., and Uzawa, H., Standford University Press, 15, 1–46, 1958. a

Vanýsek, P.: Ionic Conductivity and Diffusion at Infinite Dilution, in: CRC Handbook of Chemistry and Physics, edited by: Rumble, J. R., CRC Press/Taylor & Francis, Boca Raton, FL, 102nd edn., 2021. a

Wang, S., Cheng, Z., Zhang, Y., Jiang, L., Liu, Y., and Song, Y.: Unstable Density-Driven Convection of CO2 in Homogeneous and Heterogeneous Porous Media With Implications for Deep Saline Aquifers, Water Resour. Res., 57, e2020WR028132,, 2021. a

Weller, H. G., Tabor, G., Jasak, H., and Fureby, C.: A tensorial approach to computational continuum mechanics using object-oriented techniques, Comput. Phys., 12, 620–631,, 1998. a

Wilkins, A., Green, C. P., and Ennis-King, J.: An open-source multiphysics simulation code for coupled problems in porous media, Comput. Geosci., 154, 104820,, 2021. a

Wolery, T. J.: EQ3/6, a software package for geochemical modeling of aqueous systems: Package overview and installation guide (Version 7.0),, 1992.  a

Yan, D., Pugh, M., and Dawson, F.: Adaptive time-stepping schemes for the solution of the Poisson-Nernst-Planck equations, Appl. Numer. Math., 163, 254–269,, 2021. a

Yang, C. and Gu, Y.: Accelerated Mass Transfer of CO2 in Reservoir Brine Due to Density-Driven Natural Convection at High Pressures and Elevated Temperatures, Ind. Eng. Chem. Res., 45, 2430–2436,, 2006. a

Yapparova, A., Miron, G. D., Kulik, D. A., Kosakowski, G., and Driesner, T.: An advanced reactive transport simulation scheme for hydrothermal systems modelling, Geothermics, 78, 138–153,, 2019. a

Yekta, A., Salinas, P., Hajirezaie, S., Amooie, M. A., Pain, C. C., Jackson, M. D., Jacquemyn, C., and Soltanian, M. R.: Reactive transport modeling in heterogeneous porous media with dynamic mesh optimization, Computat. Geosci., 25, 357–372,, 2021. a

Yu, Y., Gao, W., Castel, A., Liu, A., Feng, Y., Chen, X., and Mukherjee, A.: Modelling steel corrosion under concrete non-uniformity and structural defects, Cement Conrete Res., 135, 106109,, 2020. a

Zalts, A., El Hasi, C., Rubio, D., Ureña, A., and D'Onofrio, A.: Pattern formation driven by an acid-base neutralization reaction in aqueous media in a gravitational field, Phys. Rev. E, 77, 015304,, 2008. a

Zaytsev, I. D. and Aseyev, G. G.: Properties of Aqueous Solutions of Electrolytes, CRC Press, Boca Raton, FL, ISBN 9780849393143, 1992. a

Zeebe, R. E.: On the molecular diffusion coefficients of dissolved CO2, HCO3-, and CO32- and their dependence on isotopic mass, Geochim. Cosmochim. Ac., 75, 2483–2498,, 2011. a

Zhang, Q., Tu, B., Fang, Q., and Lu, B.: A structure-preserving finite element discretization for the time-dependent Nernst-Planck equation, Journal of Applied Mathematics and Computing, 68, 1545–1564,, 2022. a

Zhang, W., Li, Y., and Omambia, A. N.: Reactive transport modeling of effects of convective mixing on long-term CO2 geological storage in deep saline formations, Int. J. Greenh. Gas Con., 5, 241–256,, 2011. a

Zhang, Z., Fu, Q., Zhang, H., Yuan, X., and Yu, K.-T.: Experimental and Numerical Investigation on Interfacial Mass Transfer Mechanism for Rayleigh Convection in Hele-Shaw Cell, Ind. Eng. Chem. Res., 59, 10195–10209,, 2020. a

Short summary
Water in natural environments consists of many ions. Ions are electrically charged and exert electric forces on each other. We discuss whether the electric forces are relevant in describing mixing and reaction processes in natural environments. By comparing our computer simulations to lab experiments in literature, we show that the electric interactions between ions can play an essential role in mixing and reaction processes, in which case they should not be neglected in numerical modeling.