Modeling Zinc Electrowinning for Current Efﬁciency Prediction Based on Nernst-Plank Equation and Electrode Gas Evolution Reaction Kinetics

Zinc electrowinning is by nature complex, involving the electrodeposition of Zn ions on cathode, oxygen evolution on the anode, as well as the hydrogen evolution on the cathode as a side reaction which decreases the current efﬁciency and increases the energy consumption while improving local mass transport. Because of the capital-intensive features and the low efﬁciencies of experimental research approaches, modeling and simulation have become a promising way to study the electrowinning processes. The coupling of gas-liquid two-phase ﬂow from gas evolution with the electrochemistry process, however, is challenging to simulate. In this paper, an approach of building a comprehensive Zn electrowinning model based on COMSOL Multiphysics was introduced. Further, the basic model results in velocity ﬁeld, concentration distribution, gas fractions, and current efﬁciency were demonstrated and discussed. In addition, validation of the current efﬁciency predictions from the model was performed through a comparison with experimental test data.©The

Zinc, which is the third most common non-ferrous metal in the world, has a wide range of significant applications in alloys, corrosion protection, battery manufacturing, etc. 1 Industrial zinc production involves approaches of hydrometallurgy, pyrometallurgy or a combination of these two processes. Approximately 90% of the world's zinc production utilizes hydrometallurgical processing methods, which provide important technical and economic advantages. 2,3 Several important electrode reactions are involved in the Zn electrowinning process, including the reduction of Zn ions and the side reaction of hydrogen evolution on the cathode as well as the decomposition of water at the anode, as shown in Equations 1 through 3.
Zn 2+ + 2e − → Zn [1] 2H + + 2e − → H 2 [2] 2H 2 O → 4H + + O 2 + 4e − [3] Due to the nature of the electrochemical reactions, Zn electrowinning is an energy-intensive process that consumes a large amount of energy, hence increasing current efficiency and reducing energy consumption become one of the most important objectives of the operation. In order to optimize the operating parameters and achieve high energy efficiency, a great deal of experimental or numerical work has been performed or developed. However, due to a large amount of electrolyte involved and the hysteresis involved (there is a time delay in the performance change with the parameter adjustment) of the experimental approach, investigating the Zn electrowinning process experimentally is not always practical and economically feasible. Therefore, modeling and simulation have become an important way to achieve optimization with relatively low cost, high efficiency, good flexibility, and increased accuracy.
In general, established models can be divided into four categories: current efficiency prediction equations, pertinent parameter empirical * Electrochemical Society Student Member. * * Electrochemical Society Member. z E-mail: zongliang.zhang@utah.edu; Joshua.werner@uky.edu; michael.free@utah. edu equations, fundamental models, and CFD models, as exemplified below. For current efficiency prediction equations, Wark 4 formulated a semi-empirical relationship between current efficiency and the ratio between the zinc sulfate concentration and the sulfuric acid concentration. Cruz et al. 5 developed a more sophisticated empirical model to predict the current efficiency and cathodic overpotential of Zn electrodeposition based on multilinear regression. Nine variables of electrolyte composition in 3 levels of their concentration were used to study their effects on current efficiency. Some works 6-8 based on experimental studies have been done to decide the empirical equations for the pertinent physical or chemical constants. For example, Guerra gave several reliable empirical equations of the electrolyte density and viscosity. 7 The equilibrium exchange current density i 0 is a function of concentration, temperature, and rate constant. The associated rate constant depends on the substrate. i 0 is generally determined by experimental methods or with empirical equations. [9][10][11] Barton 16 summarized the work that has been done on Computational Fluid Dynamics (CFD) models for zinc processing. These models contributed a lot to the understanding of the Zn electrowinning process and the development of models that describe the Zn electrowinning process more accurately. Simplifications of some significant features of the real electrowinning process, nevertheless, leave the models incomplete and less accurate than what is needed. To achieve more reliable simulations, a recently developed, more comprehensive Zn electrowinning model with the combination of two-phase flow, electrochemistry, mass transfer, gas evolution reactions, and current efficiency prediction is described in this paper.

Model Description
Zn electrowinning involves a variety of complex physical and chemical processes. To simulate the Zn electrowinning process as accurately as possible, these processes, including mass transfer, electrochemistry, fluid flow, and electrode reactions, etc., must be considered as completely as possible into the model. Furthermore, the interaction among these processes makes it harder to simulate the process. For example, as the electrochemical reactions consume or produce substances, these source terms affect the mass transfer process. In addition, the mass transfer process will affect the fluid flow in the electrowinning cell because of the density gradient created. These complexly coupled processes can be described by the corresponding governing equations, and they can be characterized and described by a model based on the Computational Fluid Dynamics (CFD) module coupled with the tertiary current distribution electrochemistry module in the COMSOL Multiphysics package. The governing equations for the Zn electrowinning process are very similar to the copper electrowinning process, which have been discussed in the previous works. [17][18] Governing equations.-Mass transfer.-In the Zn electrowinning process, the mass transfer is comprised of three processes: diffusion, migration, and convection. In the electrolyte, the governing equation for mass transfer in solution is the Nernst-Plank equation, as shown in Eq. 4 [4] where N i , z i , u i , c i , and D i are the flux density, charge, mobility, concentration, and diffusivity of species i, F is Faraday's constant, ∇ l is an electric field, ∇c i is a concentration gradient, and v is the velocity vector. The material balance in the electrolyte is governed by Eq. 5.
Because of the electroneutrality of the electrolyte, the electrolyte current density is simplified and shown in Eq. 6.
where i e is the current density in the electrolyte, and other variables were defined above.
The continuity equation simplified via the low gas concentration assumption is shown in Eq 8.
where u l is the velocity vector (SI unit: (m/s), p is the pressure (SI unit: Pa), φ l is the phase volume fraction (SI unit: m 3 /m 3 ), ρ l is the density (SI unit: kg/m 3 ), g is the gravity vector (SI unit: m/s 2 ), F is any additional volume force (SI unit: N/m 3 ), I is the identity matrix.
Geometry and mesh.-To facilitate the convergence and calculation of the model, a 2D geometry configuration was used in the model. However, the 3D schematic diagram of the electrowinning cell and the electrode configurations are shown in Figure 1a in order to give an intuitive view of the simulated cell. Figure 1b shows the side view of the cell and the 2D model simulated area were marked with a red dash line. In this geometry, one anode and one cathode as well as the electrolyte between them were considered. Note that the size of the electrodes as well as the gap between them was exaggerated in the 3D plot to show a clearer set of anode and cathode and how the 2D geometry was obtained from the 3D plot. The final 2D geometry applied in the model with precise geometric dimensioning is shown in Figure  1c. The free triangular mesh was applied in the main geometry with user-defined finer mesh in the vicinity of the cathode and extremely fine mesh on the top part of the cathode to deal with the fast fluid flow in the top part. Boundary layers were applied at the cathode surface to facilitate the mass transport calculation in the boundary layer. The mesh size distribution in the entire cell and the details in the top part are shown in Figure 2.
Boundary conditions.-Electrodes.-Equilibrium potential.-In order to determine the cell voltage, the equilibrium potential of both the anodic and cathodic reactions need to be calculated. The equilibrium potential E is given by the Nernst equation as Eq. 9.
303RT n F log a reduced species a oxidi zed speices [9]  This equation uses the standard thermodynamic potential along with the reduced species activities and the oxidized species activities to determine the equilibrium potential for a half-cell reaction.
Electrode kinetics.-The local current density is calculated with an adapted Butler-Volmer Eq. 10.
where i loc is the local current density at the interface (also called the charge transfer current density), i 0 is equilibrium exchange current density, C R,S is the surface concentration of the reduced species, C R,B is the bulk concentration of the reduced species, C O,S is the surface concentration of the oxidized species, C O,B is the bulk concentration of the oxidized species, α a is the anodic symmetry factor, α c is the cathodic symmetry factor, z is the number of electrons transferred in the rate limiting step (nearly always one), F is the Faraday constant, R is the gas constant, T is the absolute temperature, and η is the overpotential.
The corresponding standard electrode potential, the exchange current density, and the charge transfer coefficient for each reaction are shown in Table I. 10,20 Inlet, outlet, and walls.-Liquid and gas phases have different inlets and outlets. For the liquid phase, the inlet and outlet are shown in the model geometry shown in Figure 1. The inlet fluid velocity was calculated from the experimental flow rate and was set to 3.2 × 10 −4 m/s with the initial concentrations of Zn and H 2 SO 4 shown in Table II. The outlet of the liquid phase was set as a liquid pressure boundary with flows going out from that boundary naturally. For the gas phase, the inlet is the anode (oxygen) and cathode (hydrogen), while the outlet is the entire top line and the liquid outlet. Other boundaries are set as no-slip boundary walls.   Table II. The main electrolyte parameters, i.e. electrolyte density and electrolyte viscosity, have an important effect on the electrolyte fluid flow behavior, and these parameters change with the operating conditions. Therefore, in order to get correct electrolyte parameters under different operation conditions, the relationship of electrolyte density and viscosity with a variety of operating parameters are taken into account. The density and viscosity of the electrolyte are generally affected by the ionic concentrations and the electrolyte temperature. A number of researchers [6][7][8][9] have studied the changes of density and viscosity within a certain range of ion concentration and temperature of ZnSO 4 -H 2 SO 4 -H 2 O electrolyte system. In the literature, some empirical equations have been developed based on the experimental data. Among these equations, the following two equations were selected based on Guerra's recommendations 7 to calculate the electrolyte density and viscosity under different conditions in the model, as shown in Equations 11 and 12. 6 where, C Zn 2+ is Zn ion concentration, mol/L; C H 2 SO 4 is sulfuric acid concentration, mol/L; T is electrolyte temperature, K; ρ is electrolyte density, kg/m 3 ; μ is electrolyte viscosity, Pa · s. Similarly, the diffusion coefficients of the major ions will also have an important effect on the mass transport process. The diffusion coefficient will change with operating parameters.   [13] where, D i is the diffusion coefficient of species "i"; k is the Boltzmann constant; T is the temperature; μ is the solution viscosity; r i is the radius of species "i".
where, T is electrolyte temperature, K; D Zn 2+ is Zn ion diffusion coefficient, m/s 2 .
Operational parameters.-Similar to a previous copper electrowinning model, 17 21 to facilitate further model validation.

Results and Discussion
As the model is established, it was first calculated based on the basic parameters shown in Table II and Table III. As the model was built as a transient model, therefore it gives results both on the transient stage and the steady stage. At the beginning of the test, the fluid flow and the concentration gradients of all the species in the electrolyte start to develop and change drastically with time. As all these processes build up, the electrowinning process reaches a pseudo-steady state when all the parameters are relatively stable but fluctuate with time slightly. Finally, the process is stabilized and reaches a steady state, which means no apparent changes in the parameters could be observed with time. One example of the parameter development process is shown in Figure 3, which demonstrates the Zn concentration   Figure 1c. In Figure 3, it can be seen that the Zn concentration at this point changes a lot with time in the first 180 s, and stabilized at around 200 s. In constant operation, the performance of the electrowinning cell at steady-state attracts more attention. As a consequence, the results at pseudo-steady-state or steady-state are of more interest when studying the influence of operating parameters on process performance under various conditions. By comparing the model results at different times, it is determined that the process reaches a pseudo-steady-state after around 180 s of operation. Therefore, the model was run to 180 s to reach a pseudo-steady state, and then run for 30 s more based on these results. The data at 210s were used to study the performance of various steady-state phenomena and model results.
Base model results.-The model using the operating parameters in Table II and Table III was referred as the base model and used to calculate the fluid flow, electrolyte potential, species concentration distributions, current density, and current efficiency in the Zn electrowinning cell. These results are discussed and shown below. Figure 4 shows the calculated velocity field in the Zn electrowinning cell for the base model. As it shows, there is a large vortex formed at the upper part of the cell because of the strong stirring effect of oxygen bubbles. The process for the generation of this fluid flow field can be summarized below. First, the bubbles are generated at the anode surface, and they start to rise because of the effect of buoyancy. As the bubbles accelerate and join more bubbles and coalesce as they rise to higher positions in the cell, the speed and momentum of the bubbles increase and approach the maximum value at the top of the cell. The interaction between the bubbles and the electrolyte solution generates a similar fluid field in the liquid. The difference is that when the bubbles reach the top surface of the electrolyte, most of them were released to the ambient environment. However, the liquid fluid flow continues, but it is forced to change direction at the surface and heads toward the cathode and downward, thereby generating a vortex because of the interaction between the remaining bubbles and the fluid electrolyte. This explains the reason why the fluid velocity at the lower part is much slower than at the upper part. This Zn electrowinning fluid flow field is very similar to the fluid flow field in the Cu electrowinning cell studied previously 17,18 except that the maximum velocity is higher than that in the copper model because of the high current density applied. Figure 5a and Figure 5b show the Zn concentration distribution in the cell and along the cathode, respectively. In bulk solution, the Zn concentration is uniform and close to the initial or inlet concentration because of the strong stirring effect of gas bubbles, making the cell a well-mixed vessel. A boundary layer formed near cathode because of the consumption of Zn ions on the cathode. Therefore, there exists a large concentration gradient in the boundary layer. As shown in Figure  5b, the Zn concentration on the cathode surface can go lower than 100 mol/m 3 . However, there exist two peaks on the Zn concentration distribution curve along the cathode. One is located at the top area of the cathode while the other one appears at the bottom part of the cathode. The peak on the top cathode can be easily explained by the influence of the vortex. Because of the existence of the vortex and the fast fluid flow in that area, the boundary layer thickness in this area is much thinner than in other areas, which means less concentration gradient is needed to achieve the set current density, yielding a higher Zn concentration on the top cathode. Similar explanations apply to the formation of the peak at the bottom cathode. As can be seen from Figure 4, besides the vortex at the top cell, there exists another fast flow stream along the cathode. Although it is not as strong as the flow in the top part, it is still faster than most of the electrolyte flow, flowing at around 0.02 m/s -0.03 m/s. As this stream of flow hits the cell bottom, it is forced toward the bottom part of the cathode, thereby reduces the thickness of the boundary layer and boosts the surface Zn concentration in this area.
Because of the resistance in the electrolyte, there is a voltage drop in the electrolyte. Figure 6 shows the electrolyte potential distribution in the Zn electrowinning cell. From calculation results, the maximum voltage drop in the electrolyte between electrodes is about 0.4 V.
There is gas evolution on the cathode and the anode. The buoyancy force of these gas bubbles is a strong agitation force for the fluid flow. Therefore, knowing how the gas distributed in the cell will be beneficial to the understanding of the fluid pattern in the electrowinning cell. As shown in the gas fraction distribution in Figure 7a, a great majority of the gas bubbles concentrate near the anode surface. However, a small fraction of hydrogen gas bubbles formed at the cathode surface in Zn electrowinning, as shown in Figure 7b, speeding up the mass transport of the Zn and other ions near cathode surface.
Based on the kinetics of the Zn deposition reaction and hydrogen evolution reaction on the cathode, the local current densities of Zn deposition reaction, hydrogen evolution reaction, as well as the overall current density were calculated by the model, as shown in Figure 8a. The overall current density was set to 500 A/m 2 , while the current density of the hydrogen evolution reaction is around 50 A/m 2 along the cathode except for both ends of the cathode. That means the local current efficiency in the middle cathode is a little lower than at other positions, which is shown in Figure 8b. Zn reaction rate is highest at the upper part of the cathode where the vortex exists.
Current efficiency is calculated from the ratio of Zn reaction current density to the overall current density. Figure 8b shows the time average current efficiency in the 30s-time interval and the overall current efficiency along the cathode in this 30 s period. As discussed before, the current efficiency is high at both ends of the cathode and low at the middle part of the cathode. The reasons for this phenomenon could be the fast mass transport of Zn ions in these regions due to fast fluid flow and the edge effect. The overall current efficiency is about 93%, which is a reasonable value in actual Zn electrowinning process.

Model validation and current efficiency prediction.-A number
of experimental tests have been done to study the Zn electrowinning process and the associated experimental data can be used to validate model predicted results. In the study of Scott et al., 21 the effects of Zn concentration, H 2 SO 4 concentration, temperature, and current density, et al. were investigated with the experimental method. Some of these experimental results were used for the validation of model   Figure 9. 21 From Figure 9, the model predicted general trend of the current efficiencies with Zn concentration shows good agreement with the experimentally measured current efficiencies, except the model predicted current efficiency is lower than experimental data when the Zn concentration is very low (around 20 g/L). But considering that generally low Zn concentration is not common and not practical in Zn electrowinning, this prediction could be correct, but not experimentally validated. Furthermore, even under the same Zn concentration, different experiments conducted by different researchers can vary over a reasonable range because of different operating conditions. Therefore, it is deemed that the model results are good enough to predict the current efficiency and other significant parameters in Zn electrowinning.
The current efficiencies distribution along the cathode under 3 different Zn concentrations mentioned above are shown in Figure 10.
In all 3 scenarios, the current efficiencies are higher at the upper part of the cathode and relatively stable at the middle part of the cathode.
To further understand these distributions of current efficiencies, the Zn concentrations distribution along the cathode at steady-state are plotted and shown in Figure 11. From this figure, the Zn concentration on the cathode varies under different bulk Zn concentrations. Generally, the current densities have a strong effect on the Zn concentration on the cathode. For 20 g/L Zn concentration condition, the Zn concentration on the cathode is almost 0 because the Zn electrowinning is running at the limiting current density. Higher Zn concentration means a lower resistivity. Therefore, the voltage drop on electrolyte will decrease as the Zn concentration increase. It is found from this plot that the cathode surface concentration increases with the increase of bulk Zn concentration when other parameters are identical.

Conclusions
A comprehensive Zn electrowinning model based on the Nernst-Plank equation and electrode gas evolution reaction kinetics was established with COMSOL Multiphysics software package. Various model results including the velocity field in the electrowinning cell, the concentration distribution of involved species, as well as the  current density and the current efficiency distributions were obtained, demonstrating valuable insights into the complex processes in the Zn electrowinning cell. The current efficiency values calculated using the model were very close to experimental values, thereby providing significant validation for this model. The establishment of the model and its application to the investigation of the influence of the Zn concentration as well as the subsequent model validation lead to several conclusions relating to this model and its use. (1) The involvement of the fluid flow, the Nernst-Plank equation described mass transport process, and the electrode reaction kinetics of metal deposition and gas evolution on the anode and the cathode allow the model to be accurate and effective in simulating and studying the Zn electrowinning process. (2) By studying the influence of bulk electrolyte Zn 2+ and H + concentrations on the current efficiency and its comparison to the experimental data, the adequacy of the model in investigating the current efficiency was verified, hence further studies on the current efficiency could be conducted with this model. Overall, the establishment of this comprehensive Zn electrowinning model will be beneficial for the  understanding of Zn electrowinning and associated process optimization in the future.