A pore-scale model for electrokinetic in situ recovery of copper: the Influence of mineral occurrence, zeta potential, and electric potential

Electrokinetic in-situ recovery is an alternative to conventional mining, relying on the application of an electric potential to enhance the subsurface flow of ions. Understanding the pore-scale flow and ion transport under electric potential is essential for petrophysical properties estimation and flow behavior characterization. The governing physics of electrokinetic transport is electromigration and electroosmotic flow, which depend on the electric potential gradient, mineral occurrence, domain morphology, and electrolyte properties. Herein, mineral occurrence and its associated zeta potential are investigated for EK transport. The governing model includes three coupled equations: (1) Poisson equation, (2) Nernst--Planck equation, and (3) Navier--Stokes equation. These equations were solved using the lattice Boltzmann method within X-ray computed microtomography images. The proposed model is validated against COMSOL Multiphysics in a 2-dimensional microchannel in terms of fluid flow behavior when the electrical double layer is both resolvable and unresolvable. A more complex chalcopyrite-silica system is then obtained by micro-CT scanning to evaluate the model performance. The effects of mineral occurrence, zeta potential, and electric potential on the 3-dimensional chalcopyrite-silica system were evaluated. Although the positive zeta potential of chalcopyrite can induce a flow of ferric ion counter to the direction of electromigration, the net effect is dependent on the occurrence of chalcopyrite. However, the ion flux induced by electromigration was the dominant transport mechanism, whereas advection induced by electroosmosis made a lower contribution. Overall, a pore-scale EK model is proposed for direct simulation on pore-scale images. The proposed model can be coupled with other geochemical models for full physicochemical transport simulations.


Introduction
In situ recovery (ISR), also known as in situ leaching, is a technically and commercially feasible method for the extraction of minerals from high-grade ore bodies in a radically different way to conventional mining [1,2,3].ISR is the circulation of lixiviants through subsurface mineralized ore fractures to dissolve valuable minerals, especially metal minerals, without physically destroying ore formations [4,5,6,7].Compared with conventional mining, ISR has several advantages: (1) it eliminates mining costs, costs for removing ores and fragments to surface dumps, and costs of storage/disposal of tailings [2]; (2) it has lower noise, and greenhouse gas emissions [7]; and (3) it creates a safer working environment for mine workers.Starting in the early 1970s, ISR was developed and applied for uranium extraction from roll-front sandstone deposits, particularly in Kazakhstan and Uzbekistan [2,8,9,10].To date, ISR has been extensively applied in the production of other metals such as copper, gold, and lithium [11,2].Unlike sandstone deposits that contain many pores and fractures, the host formations for many metallic minerals are mainly low-permeability consolidated hard rocks.Therefore, it is difficult to inject lixiviant into hard rock formations by hydraulic force, which presents challenges for ISR and results in low recovery rates [7].
To resolve the issue of the limited permeability of hard rock formations, electrokinetic (EK) ISR (EK-ISR) has been proposed to enhance ion transportation within low-permeability formations via application of an external electric potential [12,13,14].EK is a mature technique that has been applied in many engineering fields, including soil remediation [15], wastewater treatment [16], and mine tailing remediation [17].For conventional ISR, flow is governed by hydraulic pressure gradients, whereby only a few preferential flow paths can be swept, which makes ISR highly unstable and unpredictable, as these preferential flow paths may host only a small fraction of the total metals [7,14].In EK-ISR, electric potential induces a more homogeneous flow through the heterogeneous ore.Therefore, EK-ISR is currently considered a promising method for the recovery of various metals, including gold and copper, from intact hard rocks [12,13,14].EK-ISR, However, has not yet been applied in any mining project, and thus requires further experimental and numerical work to better understand its efficacy for specific physico-chemical subsurface environments and rock surface properties.Until recently, several researchers experimentally and numerically studied EK transport in porous media and EK-ISR [18,19,20,21,14,22,23,24]. Especially, [18,22] studied the electrokinetic permanganate delivery in low-permeable porous media and found that electrokinetics significantly enhanced fluid permanganate delivery.[24] conducted lab-scale experiments on EK-ISR with different voltage and pressure field and [14] used COMSOL Multiphysics coupled with Phreeqc to perform a continuum-scale simulation of EK-ISR.However, to our best knowledge, the simulation and experimental EK-ISR study and the comparison between them are still lacking.This digital twin study can be potentially performed with micro-CT imaging and large-domain pore-scale direct simulation [25,26,27,28].The information extracted from the digital twin study can be used for the upscaling study.
The main transport mechanisms of EK-ISR include (1) electromigration, which involves the displacement of charged species towards an electrode of opposite charge, and (2) electroosmotic flow (EOF), which represents the net movement of fluid flow as a result of the excess charge adhered to the mineral surface [29].The electrical double layer (EDL) plays a key role in EOF.The EDL consists of a charged solid surface and a thin layer of counter ions in an aqueous solution.As counter ions in the EDL move towards the oppositely charged electrode, momentum is transferred to the surrounding fluid molecules, thereby inducing flow [29].The Helmholtz-Smoluchowsky equation (HS) is commonly used to determine EOF in porous media, and its application depends on the thickness of the EDL [29,30,31].Most studies in the literature assume a thin double layer, which means that the thickness of the EDL is considerably smaller than the pore size [30,31].The counter ions of the EDL screen the wall charge within a region that scales with Debye length [32].In addition, the thickness of the EDL varies with the electric potential at the particle surface.The zeta potential is a way to characterize the EDL based on the ionic concentration in the EDL and the pH as a result of the protonation/deprotonation reactions that occur at the particle surface [33,34,31,35,36].Other essential procedures and parameters for EK including geochemical reaction and local pH distribution as well as some larger-scale factors such as tortuosity and porosity were characterized for fluid-solid systems [37,38,39,40,31,20,41,42]. Geochemical reactions change ion composition and concentration and therefore, changes the thickness of EDL and zeta potential [43,39].The dissolution of the mineral during reaction will change the porous structure and flow pattern.These changes due to geochemical reaction influence the electroosmotic permeability and electromigration.Considering EK-ISR, the ionic concentration is high and its effect on the ion transport becomes non-trivial.Meanwhile, the local pH heterogeneity causes the heterogeneity of zeta potential and results in a nonlinear response of the electroosmotic velocity.With a significant change of pH, the zeta potential might be reversed and results in the reversed electroosmotic velocity [31].Tortuosity and porosity provide the morphological information of the porous ore and provides the bridge to study the micro-scale and macro-scale relationship for future upscaling [44,45,46].[46] compared the permeability estimated by EK and generated from flow experiments, the two permeability values are consistent after reconciling through the hydraulic tortuosity.[19,23] studied the effect of the heterogeneity porosity to the EK and found that the porous heterogeneities play an important role in EK and the coupling to hydraulic process.
Most of the metal ore formation for ISR is composed of a low porosity-permeability hard rock system, with only thin fractures presented.In this case, the thickness of EDL could be comparable to the fracture size.EDL can then be fully described and resolved.If these fractures are connected to the inlet, the velocity front for the lixiviant changes from plug-shape into parabolic-shape [31].For a negatively-charged mineral surface in the fracture, the flow velocity towards the cathode is lower near the mineral surface.Therefore, to understand EOF in EK-ISR, the thickness of the EDL and zeta potential must be characterized based on the chemical condition of the ore.EOF in microchannels has been extensively studied because of its significant applications in EK remediation [47,48,49,50].However, most studies are based on simplified pore-structure models [30,51,31,52].[51] developed a numerical method to simulate electroosmotic flow in a 2D microchannel, which consisted of the combination of nonlinear Poisson equation for the electric potential with lattice Boltzmann method (LBM) for fluid flow [51,53,54,55].Most EK studies using LBM are based on simultaneously solving the Poisson, Nernst-Planck, and Navier-Stokes equations in a self-consistent scheme.[51] studied the effect of electrically-driven and pressure-driven flows on flow velocity, electro-viscous effect, and electroosmosis in homogeneous microchannels.The LBM has also been used to study EOF in a heterogeneous pore structure reconstructed using random porous structures [31].They studied the effects of the heterogeneous zeta potential at the mineral-liquid interface at different pH values and reported that for a small electric potential strength, the effect of the electrical force on the distribution of pH causes a nonlinear response in the electroosmotic velocity.
Such studies are examples of the successful application of numerical schemes based on the LBM to analyzing EK flows.
However, these studies mainly focused on the fundamental physical perspective of EOF and used manually generated porous media as the study domain.At the same time, the effect of EOF on the flow behavior in complex porous media has not yet been investigated.Meanwhile, a more complex chalcopyrite-silica system is investigated and the efficacy of ion transport under various EK conditions and mineral distributions is evaluated.Specifically, the EK model based on LBPM was first validated against COMSOL Multiphysics in terms of EOF in cases where the EDL is both resolvable and unresolvable in a 2D microchannel.After validation, a chalcopyrite-silica system was created by mixing silica and chalcopyrite powders.
The synthetic ore was imaged using high-resolution X-ray micro-computed tomography (micro-CT), which allowed the visualization of the 3D structure and mineral distribution within the system and subsequent direct simulation with LBPM.Overall, the proposed model was designed to investigate the transport of fluid/ions to the target metals under electric potential at the pore-scale.For future work, considering the importance of surface potential in this study, a surface complexation model is planned to be introduced to characterize zeta potential at liquid-solid interfaces.The proposed model can be coupled with a geochemical solver such as PhreeqcRM [56,57,58] for the full characterization of physical-chemical process in the EK-ISR.

Numerical Methods
The governing model of the EK flow includes three coupled equations: Poisson equation for the electric potential, Nernst-Planck equation for ion transport driven by chemical and electric potentials, and Navier-Stokes equation for the flow of an electrolyte solution carrying ions.

Mathematical Models
The flow of the electrolyte solution is governed by the incompressible conservation of mass and Navier-Stokes equation: where u is the fluid velocity vector, ρ 0 is the fluid density, p is the fluid pressure, µ is the dynamic viscosity, and F is the body force, which in this study is primarily caused by an external electric potential.The Navier-Stokes equations were solved in pore spaces, whereas a standard nonslip boundary condition was applied to the solid spaces.In the case where the thickness of the EDL was much smaller than the characteristic length of the simulation, that is, below the resolution of an input image, an electroosmotic velocity boundary condition was introduced in the EDL which excluded the detailed flow field between the solid and slipping plane, and analytically calculated the velocity at the solid according to the local zeta potential of the solid surface.Herein, we adopted the commonly used HS equation where Ω denotes the fluid domain, ψ is the electric potential within the electrolyte, ζ is the local zeta potential of the solid surface, is the permittivity of the electrolyte solution, and ∇ T is the tangential part of the gradient operator, perpendicular to the solid surface orientation.When the electroosmotic velocity boundary condition is applied as the driving force of the flow, the electric body force in Eq.1 is set to zero, because the EDL is below resolution and the bulk fluid is considered electrically neutral [31].
Ion transport is governed by the Nernst-Planck equation, which incorporates electrochemical migration as an extra drift term into the mass flux: where C i is the concentration of the ith ion, z i is the ion algebraic valency, and V T = k B T /e is the thermal voltage, where k B is the Boltzmann constant and e is the electron charge.The ion mass flux is affected by three factors: the first and second terms on the left-hand side of Eq.3, which are the convection and electrochemical migration, respectively; the term on the right-hand side is the diffusion where D i is the diffusivity of the ith ion.
A non-flux boundary condition at the fluid-solid interface was applied to the ions: where n s is the unit normal vector of the solid surface, and The electric potential of the distribution of ions was solved by the Poisson equation: where ε 0 is the permittivity of vacuum, and ε r is the dielectric constant of the electrolyte solution.The net charge density ρ e (C/m 3 ) is related to the ion concentration as follows: where the sum runs over all ionic species and F is Faraday's constant given by F = eN A , where N A is Avogadro's number.The force on the body owing to the net charge density in the Navier-Stokes equation is given by: The fluid-solid boundary condition for the electric potential is typically specified in two forms: 1) the surface charge density σ e and 2) surface potential ψ s at the solid surface (when the EDL is unresolved, the surface potential is equivalent to the zeta potential).The former is a Neumann-type boundary condition given by whereas the latter is a Dirichlet-type boundary given by where ψ s is the user-specified electric potential of the solid surface.A detailed flow chart of how these equations were solved is shown in Fig. 1.
Figure 1: Flow chart for the proposed EK model.Debye length was calculated and compared with voxel length to check if the EDL is resolvable.

Lattice Boltzmann Methods
To solve the coupled transport and electric equations (mentioned above) in porous media, we adopted the commonly used LBM because of its inherent scalability of parallel computation and efficient handling of complex boundary conditions.Several coupled lattice Boltzmann (LB) frameworks dedicated to EK flow have been developed and studied over the past decade [51,53,31,54].Herein, we adopted the method proposed by [59], which was modified by [53,31], to incorporate an electroosmotic velocity boundary condition.
The LB method naturally suits parabolic partial differential equations.The Poisson equation, however, is elliptical, and thus an artificial time-dependent term is usually added so that the LBM yields a steady-state solution of the 'transient' Poisson equation of the following form: We deployed the D3Q7 lattice to solve the Poisson equation.The corresponding LB evolution equation for the distribution function h q of the electric potential ψ is given by with and the equilibrium distribution is where ξ q and ω q are the D3Q7 lattice velocity vector and weighting coefficient, respectively, with ω 0 = 1/4 and ω 1−6 = 1/8; and x is the spatial resolution of the simulation domain.It should be noted that the time in the LB Poisson equation is denoted as t; thus, it should be differentiated from the LB time t in the Navier-Stokes and Nernst-Planck equations to be covered later because only the steady-state solution of Eq.10 is of interest.Within each main evolution step ( Fig. 1 for flow chart), Eq.11 is executed iteratively until the standard mean squared error over certain amount of timestep, t , is smaller than the user specified tolerance: where N is the total number of fluid nodes, and err is the prescribed tolerance.
The LB relaxation time parameter τ ψ in Eq.11 is given by: where c 2 s is the LB speed of sound; for the D3Q7 lattice c 2 s = 1/4.
For the boundary condition, we adopted the formulation developed by [53], which is a completely localized scheme that is more suitable for complex porous media.For the Neumann-type boundary condition, where a surface charge density σ e is specified, after the LB collision, the normal streaming step is replaced by the following equation: where h q (x, t) is the post-collision distribution, and the index q indicates the direction opposite to q.The Dirichlet-type boundary condition, where the surface potential is specified, is given by: Incidentally, the electric field E is given by the gradient of the electric potential, that is, E = −∇ψ.According to [53], the gradient can be calculated locally as where the index α denotes Cartesian coordinates.
For ion transport, the LB evolution equation was also solved for the D3Q7 lattice and is given by for the ion distribution function g q .The equilibrium distribution function is given by where t Di is the time resolution; the method of determining t Di is covered further on; C i = q g q is the ion concentration of ith species; and u EP,i is the electrophoretic velocity of the ith ion in response to applied electric potential.The electrophoretic velocity is given by where D i,LB is the diffusivity of the ith ion species in the LB unit; it is related to the relaxation parameter τ Di as follows: For the Dirichlet-type boundary condition, that is, if the surface ion concentration C s is specified, the normal streaming step after LB collision is replaced by Here, g q denotes the post-collision distribution.For the non-flux boundary condition in Eq.4, it has been proved in [53] that it is equivalent to the half-way bounce-back boundary condition widely used in the LB method, which ensures no ion flux across the solid boundary: Regarding the LB Navier-Stokes solver for the electrolyte solution, because it has been extensively studied and used in numerous publications, the details of the formulation are not repeated here.Therefore, we implemented the formulation by [60,61,62], where a multi-relaxation LB method is deployed, and incorporated the slipping velocity boundary condition proposed by [63,31] into model cases in which the EDL is not resolved.
When solving a multi-physics problem where each transport equation has its own time scale and internal LB timestep, it is important to ensure that all of the coupled equations are synchronized in terms of a physical time scale.In other words, the relationship The time conversion factor t u for fluid flow is determined by the following relation: where ν phys and ν LB are the fluid kinematic viscosities in the physical and LB units, respectively.Notably, ν LB is linked to the Navier-Stokes LB relaxation time by ν LB = (τ u − 0.5)/3, where τ u is usually taken between 0.5 and 2 for numerical stability.
The time conversion factor t Di for ion transport is determined using the following relation: where D i,phys and D i,LB are the diffusivities of the ith ion in physical and LB units, respectively.Note that D i,LB is related to the LB relaxation time through Eq.22,where τ Di is set to 1.0 for numerical stability.
In summary, using the image resolution x, the input physical parameters (ν phys , D i,phys ) and user specified LB relaxation time (τ u and τ Di ), the time conversion factor for each solver was determined; this step was performed using a multi-physics controller in LBPM.The internal LB timestep of each solver was also subsequently determined based pressure-driven flow in a simple heterogeneous system.EK simulations were then performed on the chalcopyrite-silca system to evaluate the effect of the zeta potential and electric potential on the feasibility of EK transport for EK-ISR.
All simulation parameters are listed in Table S1 in Supplementary Materials Section 3. Simulations were performed on a local workstation with 64-core CPU, 24 GB of GPU memory and 256 GB of RAM.Simulations were computed on the GPU with a much faster computational speed than the CPU.

Micro-CT imaging and Image Processing
A chalcopyrite-silica system, which is a mixture of a Cu mineral (chalcopyrite) and gangue mineral (silica), was prepared to obtain a copper-rich porous system.This chalcopyrite-silica system was imaged using micro-CT to generate a digital 3D model for simulation.The chalcopyrite powder was obtained from Kremer Pigmente (Germany) with a particle size of approximately 80 µm.The mineral and elemental contents of the chalcopyrite powder were quantified using X-ray diffraction (XRD) analysis and X-ray fluorescence (XRF), and the results are shown in Table 1.The XRF results show that the major elements are Fe +3 and Cu +2 which are the main elements in chalcopyrite (CuFeS 2 ).It can be further confirmed from XRD that the powder contained over 72% chalcopyrite.To prepare the synthetic ore system, SiO 2 (2.5 g) and copper ore (40 mg) were mixed in NaCl solution (1.37 mL, 0.1 M) in a beaker (50 mL).The mixture was stirred for 5 min to ensure complete saturation of the system.The saturated mixture was moved to a container (height: 3 cm; inner diameter: 1cm; outer diameter: 1.2 cm) for micro-CT scanning, as shown in Fig. 2 (a).
Cotton cloth was used to cover the top of the chalcopyrite-silica system, ensuring minimal powder movement during the micro-CT scanning.The micro-CT image of the chalcopyrite-silica system is shown in Fig. 2  Identifying the mineral composition and distribution in the synthetic ore system was a pre-process for EK modeling using LBPM.Therefore, the first step was to perform multiphase segmentation of the synthetic ore system.The synthetic ore system was a mixture of silica, chalcopyrite, NaCl solution, and unsaturated voids.To segment these four phases,  we used trainable WEKA segmentation [64] to generate a training dataset of registered images that were then fed to a U-ResNet convolutional neural network(CNN) [65,66,67,68].The workflow is outlined as follows.
1.A 2D grayscale image was selected, and its pixels were manually clustered into the 4 phases as training data for the WEKA segmentation.
2. After training, the WEKA segmentation generated a 2D segmented slice of the input 2D grayscale image which works as ground truth for CNN training.
3. Once the CNN was trained, it could segment the entire 3D image of the synthetic ore system.The detailed U-ResNet architecture and training schedule can be found in Supplementary Materials Section 1.
Figs. 2 (c-d) show the 3D segmentation result of the 4-phase synthetic ore system, and Fig. 2 (e) shows the 3D distribution of silica and chalcopyrite, which are homogeneously distributed throughout the system.
3 Results and Discussion

Electrokinetic Model Validation
A detailed validation of the EK model built in LBPM is presented in this section.The EOF effect is validated against COMSOL Multiphysics in a 2D microchannel.The constant physical parameters used for validation are listed in Table S2.After validation, single phase flow behavior under only hydraulic pressure and electric potential is compared in the multi microchannels system.For the EOF validation in a 2D micro channel, as shown on the left of Fig. 3, two parallel plates are located at X = ± 0.5H, the domain is uniformly filled with a 1:1 electrolyte solution (NaCl).For the input micro channel to the LBPM and COMSOL, the pixel element is set to be 0.001 µm.Therefore, the total L and H of the channel are 0.1 µm and 0.05 µm, respectively.A meshed domain as the input for COMSOL is shown on the right of Fig. 3. Debye length plays an essential role in the flow behavior for electroosmosis, which is a function of ionic concentration.Fig. 4 demonstrates how the ionic concentration affects the Debye length and therefore influences the velocity profile.For the region where the EDL is resolvable, the velocity profile is parabolic-like.With an increase of ionic concentration, the EDL becomes unresolvable, and the velocity profile changes from parabolic-like to plug-like.
Figure 3: A visualization of the micro-channel and the meshed domain used for validation against COMSOL has 50 pixels in X direction and 100 pixels in Z direction.The physical size is the same as the input for LBPM.
Figure 4: The effect of ionic concentration on the EOF velocity with voxel size of 0.02 µm and a 0.02 V of electric potential.The Debye length is influenced by ionic concentration.When the ionic concentration is high enough, the Debye length becomes much smaller than the voxel size.The EDL is unresolvable and, therefore, the velocity profile appears plug-like.When the ionic concentration is small, the Debye length is comparable to the voxel size.Therefore, the EDL is resolvable and the velocity profile changed from plug-like to parabolic-like.
In the following validation, the EK model in LBPM is validated against COMSOL in both scenarios where the EDL is resovable and unresolvable.To meet the condition where the EDL is fully resolvable (M = H λ EDL >0.01) [31], the bulk ionic concentration is set to be 0.1mol/m 3 .The Debye length is determined as 0.0302µm and M = 1.656 (in the range of fully resolvable).The surfaces of the plates are charged with a constant charge density, σ.For EOF, an electric potential (E) is imposed along the channel in the z-direction.Coulomb's force (electromigration) acts on the electrolyte solution through the net charge within the EDL, and consequently the EOF occurs in the z-direction.Several EOFs under different σ and E at steady state are compared between LBPM and COMSOL.when the EDL is fully resolvable (M>1) [31].A visualization of the whole velocity domain is provided in Fig. 5 (b-c) to further demonstrate that the LBPM result achieves good agreement against the COMSOL result for the entire velocity domain.After validating the fully resolved EOF, the following validation is when the EDL is unresolvable.In this case, the thickness of the EDL is much smaller than the microchannel width so that the charge effect within the EDL is negligible and the domain remains electroneutral.The Helmholtz-Smoluchowski (HS) slip velocity is used to describe the slipping condition of the wall as a function of external electric potential and zeta potential.To meet the requirement that the EDL thickness is much smaller than H, the pixel size is set to be 1µm (H = 50µm), which is within a typical range of micro-CT image resolution.The ionic concentration is 0.1mol/m 3 , M is then calculated, which is 0.0006 ( 0.01).
The validation result is shown in Fig. 6.A piston-like velocity profile is obtained with the HS boundary condition.
The same shock front velocity (electroosmosis slip velocity) is obtained from LBPM and COMSOL under varied zeta potential and electric potential.

EK versus hydraulic pressure driven transport
To understand the difference between single-phase flow behavior under either solely EK or hydraulic pressure, a singlephase flow simulation was performed using the validated EK model in a 100 × 100 2D multi-microchannel system with a resolution of 4µm.The system contains only one type of solid phase with ζ = −0.02v.Three microchannels with sizes of 8, 20, and 40µm were established in the system.The hydraulic pressure ∆P was 3000kP a/m, and the electric field ψ was 2.5V /m.The simulation settings of the lixiviant and diffusion coefficient are listed in Table 2 and 3.For pressure-driven simulation, the external electric potential is set to be zero and only a pressure difference is applied.
The velocity profile results are shown in Fig. 7.With only ∆P, the flow was heterogeneous in the three microchannels, because with large channel size, the capillary entry pressure for fluid flow is lower than that of smaller channel size [69].Therefore, fluid prefers to enter the largest-sized channel owing to less hydraulic resistivity, as shown in Fig. 7(a).
Consequently, the flow velocities were greater in the larger channel.For ISR, only minerals that resided in the larger pores and fractures with preferential flow would be recovered efficiently, whereas minerals in the smaller pores would be transported less efficiently.On the other hand, a promising flow behavior was obtained when applying an electric potential, as shown in Fig. 7(b).The velocity in all three channels is homogeneous, and the velocity profile is sharp, indicating that all minerals in the pores or fractures can be efficiently exposed to lixiviant, independent of the pore size.Based on these results, it can be seen that flow driven by EK mechanism enables the minimization of the effect of pore-scale heterogeneity, thus promising a higher mineral recovery efficacy during ISR.Table 3: Settings for numerical simulation, which are the same as the previously reported experimental conditions (M represents molarity).At inlet and outlet, the numerical settings are the boundary conditions, and in the system, the numerical settings are the initial condition.

EK in a Chalcopyrite-Silica System
To apply EK-ISR, the migration of Fe 3+ through the ore system must be ensured for the dissolution of chalcopyrite into the fluid phase.Therefore, the discussion of Fe 3+ ion is focused on this section.Based on the validation results, the EK model in LBPM showed good agreement with COMSOL and was thus used for studying the EK transport in a complex chalcopyrite-silica system.To reduce the computational time for the simulation, a representative 256 cubic voxels subdomain was cropped from the original domain, as shown in Fig. S2 in the Supplemental Materials, Section 2. Representative elementary volume analysis of the chalcopyrite-silica system was also performed, and the results are provided in the Supplemental Materials.The typical chalcopyrite content in the sulfide ore is around 0.3-5 % [70,71,72].The mineral composition of the representative sub-domain that is used in the following simulation is 96.3% silica and 3.7% chalcopyrite, which is in the typical range.Moreover, to check the physical accuracy of the proposed EK model, and also considering the heterogeneous distribution of chalcopyrite in the small sub-section of the whole sample, two manually generated subdomains which have chalcopyrite concentration of 15.6% and 63.3 % were created and used to check the physical accuracy of the model and study the effect of zeta potential for different percentages of chalcopyrite.These subdomains were obtained by changing the silica grains to chalcopyrite grains.The material compositions for all simulated representative sub-domains are listed in Table 4.The following numerical studies were conducted on these sub-domains.

EK under Different Zeta Potentials
The simulation settings were set to match the EK-ISR laboratory-scale experimental settings in [14], as shown in Table 3.At inlet and outlet, the 0.2M HCl is set as boundary condition.In the system, 1 × 10 −7 M is set as initial condition.
During the simulation, H + is transported into the system under diffusion and electromigration.Chalcopyrite dissolution by the lixiviant (F eCl 3 ), which is the reaction with dissolved F e 3+ is described as follows: The oxidative dissolution of chalcopyrite occurs due to the reaction with F e 3+ ions present in the lixiviant.Under an electric field, F e 3+ moves from the source reservoir to the copper bearing reservoir, where Cu 2+ is leached.
The dissolved Cu 2+ and other cations are then transported to the target reservoir via electromigration.The reactive surface area plays a fundamental role on the mineral dissolution and reactive transport processes.The porous media properties such as tortuosity, porosity and permeability evolve as the mineral dissolution takes place.Such changes in the pore structure alter the magnitude of the velocity field and control the flow pattern.Moreover, mineral dissolution modifies the local pH in the system, affecting the zeta potential which has a great influence on the electroosmotic permeability.Therefore, understanding F e 3+ transport in the domain is non-trivial.The pH at the inlet and outlet was set to 0.7 (initial strongly acidic conditions).Other settings are listed in Table S1.The zeta potential ζ for silica and chalcopyrite is another important factor that controls the fluid/ion flow behavior.Moreover, transport is sensitive to the pH of the solution.For silica and chalcopyrite, experimental measurements of ζ at various pH values have been reported in [73,74].Under pH around 0.7, the ζ value for silica was approximately −0.005V .For chalcopyrite, ζ was approximately 0.002V , whereas ζ of chalcopyrite after the treatment with ferric chromium lignin sulfonate is approximately −0.03V .To further investigate the effect of ζ on EK, three simulations using three subdomains (3.7, 15.6, and 63.3% chalcopyrite) with ζ = 0.002V were compared.
With a high ionic concentration and large voxel size (5.4 µm), the EDL fell into the unresolvable regime; therefore, Eq.2 was used to compute the flow velocity.Using the HS equation, the velocity rapidly converged.The stopping criterion for the steady-state condition was the average velocity difference between the two time steps less than 10 −9 m/s.electroosmosis acted in the opposite direction to electromigration.However, when ζ for chalcopyrite was −0.03V , the fluid surrounding the chalcopyrite flowed at a higher velocity than that surrounding the silica, indicating that the effects of electroosmosis and electromigration were codirectional.The average velocities and electroosmotic permeabilities are listed in Table 5.The average velocity for ζ = −0.03Vwas slightly higher than that for ζ = 0.002V , owing to the EOF.The electroosmotic permeability governs the fluid flow under the electric potential similarly to how the hydraulic conductivity governs the flow under a hydraulic gradient.Therefore, a larger electroosmotic permeability was obtained for a more negative ζ.The difference is 7.28% because the amount of chalcopyrite in Domain1 is small compared to that in silica.
To investigate the effect of chalcopyrite occurrence, Domain2 and Domain3 were compared with Domain1 under the same conditions, where ζ = 0.002V .The results are presented in Table 6.Ion flux owing to advection, diffusion, and electromigration are also reported.With an increase in chalcopyrite fraction, the average velocity decreases significantly (up to 76.21% difference between Domain1 and Domain3) because the velocity near the chalcopyrite is negative owing to the positive zeta potential, leading to a decrease in the average velocity.For F e 3+ ion flux, with an increase in chalcopyrite fraction, the advective flux decreases because advection is related to fluid velocity.However, the diffusive flux remained almost constant for all the three domains, whereas the electrical flux slightly decreased.Comparing the ionic flux, the flux of advection due to EOF is 4 magnitudes less that electromigration and diffusion.Péclet number (P e) that calculates the ratio between advection transport and diffusion transport is also reported for all ions, as defined in Eq.
(1) and ( 2) in Supplementary materials.P e provides an indication of the dominant transport mechanism in the system independent of the geometry size of the domain.The advection term in P e is contributed by EOF and electromigration.
The P e result is listed in Table 1 in Supplementary materials.The P e shows that for all ions (except Fe 3+ ), the P e is around 1, while for Fe 3+ , the P e is around 3. The EOF and eletromigration induced velocity value is also reported.
Overall, the results reveal that for EK-ISR, where an external electric potential is applied, electromigration is the dominant ion transport mechanism.The diffusive flux is the second important transport mechanism in our simulation settings because the domain initially contained no F e 3+ , resulting in high diffusive flux.EOF contributes few to the ions transport.Percentage of difference 0.002% 0.01% Under an external electric potential, positively charged ions move towards the cathode, and negatively charged ions move towards the anode owing to electromigration.Fe 3+ , therefore, moves from the anode to the cathode, as shown in Timesteps 1-3 in Fig. 9 (a) and (b).The physical time corresponding to them is 3.6, 9.6, and 14.5sec, respectively.
which illustrates the transport of Fe 3+ ions under applied electric potential, .As Fe 3+ moves into the ore system, a larger region of chalcopyrite is exposed, which is essential for the reaction and recovery of Cu The velocity towards the cathode under a higher electric potential was greater than that under a lower electric potential at negative ζ for both silica and chalcopyrite.This can be explained using Eq. 2, where the velocity is proportional to the electric potential.The average velocity exhibited the same trend, as listed in Table 7. Notably, if ζ is positive for chalcopyrite, under a higher electric potential, the average velocity towards the cathode becomes smaller than that under a lower electric potential.This is because at positive ζ, electroosmosis leads to flow in the opposite direction (negative velocity) towards the anode at the locations surrounding chalcopyrite.
Additionally, the electroosmotic permeability is the same in both cases, which indicates that the electroosmotic permeability is independent of the external electric potential.This condition is valid for this case because a constant of 0.2M HCl was applied at the inlet and outlet, resulting in a homogeneous pH distribution in the system.However, this condition is not met when there is a variation in pH.For example, [31] found that electroosmotic permeability can change under heterogeneous pH conditions.Fig. 12 shows the relationship between F e 3+ and chalcopyrite at different timesteps.When ψ = 0.4V , at an early timestep, the majority of F e 3+ concentration near the chalcopyrite is in the range of 0 to 20 mol/m 3 , and only one voxel has an F e 3+ concentration of > 501mol/m 3 .At the late timestep, the number of voxels in the range 0 to 20 mol/m 3 decreases, whereas the number of voxels in other ranges increases, which is especially significant for the F e 3+ concentration of > 501mol/m 3 (from 1 to 2278).When ψ increases to 1.2V , as the simulation progresses, more voxels have F e 3+ concentration higher than 21 mol/m 3 , particularly with an increase in F e 3+ concentration greater  against COMSOL Multiphysics in terms of EOF in a 2D microchannel.The thickness of the EDL and its effect on the numerical model were discussed and validated.Specifically, at the thickness of EDL comparable to the domain size, the EDL can be fully resolved by the EK model, whereas at the thickness of EDL smaller than the domain size, the EDL is unresolvable by the model.Therefore, the HS equation was used to compute the slip boundary velocity.Good agreement was obtained between the EK model and COMSOL Multiphysics.
Subsequently, a chalcopyrite-silica system consisting of chalcopyrite and silica powder was prepared and imaged with micro-CT.The simulation with the chalcopyrite-silica system provides a more complex porous structure for characetrising fluid and ions flow under the EK condition.The simulation application studied the EK processes under various mineral occurrence, zeta potential, and electric potential.The results highlight the important influence of mineral occurrence, zeta potential, and electric potential to electroosmosis and electromigration for EK.The flexibility of the model opens the opportunity to coupling the surface complexion model for local pH characterization, such as 1-pk model [75] and triple layer model [76], coupling with a geochemical model, such as PhreeqcRM to fully capture the reactive transport under electric potential [58], and simulating large domain ore sample which is used in experiments.
Such model application can be applied for EK-ISR where the low permeability-porosity system is presented [14].The CNN model was trained for 100 epochs with an initial learning rate of 0.0001 and batch size of 16 using the Adam  For each full size slice segmentation, the approximate training time was 12 seconds.Therefore, 10 hours were used to 10 segment the whole 3D synthetic ore system.

Péclet number Calculation
Péclet number can be calculated as where µ is the average pore velocity, and L is the characteristic length, calcualted by .
where V is the bulk volume of the sample, and A is the surface area of pores.
Herein, we developed a EK model for simulation the physical transport of fluid and ions for EK using the lattice Boltzmann-Poisson method (LBPM) based on the fundamental principles of surface chemistry and EK transport.The governing model includes three coupled equations: (1) Poisson equation, (2) Nernst-Planck equation, and (3) Navier-Stokes equation.In this study, we describe the main workflow of our model and validation in terms of electroosmosis.

Figure 2 :
Figure 2: (a) Micro-CT image of the wet chalcopyrite-silica system.The system is a mixture of silica and chalcopyrite powders saturated with sodium chloride solution.(b) Segmentation process identifying 4 phases, including void, solution, silica, and chalcopyrite.Solution used herein is the 0.1 M NaCl.Multi-phase segmentation is based on CNN.(c) Segmented 3D synthetic ore system.(d) 3D distributions of silica and chalcopyrite in the system.

Fig. 5 (
Fig. 5 (a) shows the parabolic-like velocity profile that is the average velocity along the z-direction.The velocity profile of the LBPM and COMSOL show good agreement in all cases where σ = −1 × 10 −4 C/m 2 , −2 × 10 −4 C/m 2 and E = 1 × 10 5 V /m, 2 × 10 5 V /m, 3 × 10 5 V /m.The EOF is more significant under a large electric potential and surface charge, and thus results in a larger fluid velocity.The parabolic-like velocity profile matches the typical velocity profile

Figure 5 :
Figure 5: (a) Fully resolvable EOF validation results for the LBPM against the COMSOL under varied conditions.Good agreements are achieved in all cases.(b-c) The EOF velocity profile under the condition that EDL is fully resolvable for σ = −2e − 4C/m 2 , E = 1e + 5V /m.(b) The velocity profile generated from LBPM.(c) The velocity profile generated from COMSOL.Good agreement is found in the velocity at all locations of the domain.

Figure 6 :
Figure 6: Left: EOF validation results when the EDL is unresolvable for LBPM against COMSOL, under varied electric potential gradient.Right: EOF validation results when the EDL is unresolvable for LBPM against COMSOL under varied ζ.Good agreements are achieved in all cases.

Figure 7 :
Figure 7: Single-phase flow of lixiviant in a heterogeneous micro-channel system under only (a) hydraulic pressure gradient or (b) electric field.Velocity profile shows that heterogeneous flow is obtained with hydraulic pressure, whereas a homogeneous flow is obtained with electric potential.

Figure 8
Figure 8 shows the electric potential, velocity profiles for ζ = 0.002V and ζ = −0.03V, and the difference between the resulting velocity profiles.Electroosmosis is related to the ζ value of each mineral.For a negative ζ, the fluid moves towards the cathode, whereas for a positive ζ, the fluid moves towards the anode.The major mineral in the subdomain was silica with ζ = −0.005V .The fluid surrounding the silica flowed towards the cathode (positive velocity value).For ζ of chalcopyrite of 0.002V , the fluid flowed towards the anode (negative velocity value), indicating that

Figure 12 :
Figure 12: F e 3+ distribution is divided into seven concentration ranges.The voxel frequency refers to the number of voxels that are connected to the voxels of chalcopyrite and also fall into the satisfied concentration range.The results in early and late timesteps are compared for ψ = 0.4V and ψ = 1.2V .For the range > 501mol/m 3 , except for the late timestep for ψ = 1.2V , the voxel frequency for others are 1.

1 CNN Architecture and Training Schedule 1
Training data for U-ResNet [?] (architecture was shown in Fig.S1contained the grayscale micro-CT slice and its 2 corresponding ground truth from the WEKA segmentation tool with a size of 2013 × 2013 pixels.It was cropped into 3 patches with a size of 100 × 100 pixels.Overall, 400 patches were obtained, and split 80/20 for training and testing.

5
solver[?].The learning rate was reduced by a factor of 0.5 each time the loss reaches a plateau after ten epochs.The 6 training was implemented in PyTorch, using an Nvidia RTX 3090 GPU and 16 cores' CPU (AMD 5950X) with a 7 memory of 128 GB.The approximate training time per epoch was 5 seconds, the model was trained in 9 minutes for 8 100 epochs.After training, other 2969 full-size 2D slices were passed to the CNN model one-by-one for segmentation.

Figure S1 : 13 Fig. S2 .
Figure S1: Architecture of the U-ResNet used for segmentation, containing downsampling process as image feature extraction and upsampling layer for useful feature decoding.
on t u and t Di .For example, if t u is the largest among the sets { t u , t D1 , t D2 , ..., t Dn }, then N tu is set to 1 and N t D i can be determined using N t D i = t u / t Di , rounded up to the nearest integer.The open-source code for our EK model is accessible on GitHub (https://github.com/OPM/LBPM).Our EK model was first benchmarked with COMSOL Multiphysics for EOF under both conditions, where the EDL is resolvable and unresolvable in a two-dimensional microchannel.The validation results can be found in Supplementary Materials Section 3.After validation, we use a simple EK model to evaluate the benefits of EK transport compared to

Table 1 :
XRF and XRD results of the chalcopyrite powder.For XRF, powder was preoxidized before being fused into glass bead.Therefore, the XRF results show the element oxide weight percent.The XRD pattern shows the minerals occurring in the powder.Both results confirm that the most abundant mineral in the powder is chalcopyrite.SiO 2 SO 3 K 2 O CaO TiO 2 Fe 2 O 3 CuO

Table 2 :
Physical parameters used in the EK numerical simulation.
3.1.1Validation of EOF for the EK model

Table 4 :
Mineral compositions of all three sub-domains used in this study.Domain1 is the real scanned sub-domain.Domain2 and Domain3 are the manually generated sub-domains.