Characterization of gas diffusion electrodes for metalair batteries
Abstract
Gas diffusion electrodes are commonly used in high energy density metalair batteries for the supply of oxygen. Hydrophobic binder materials ensure the coexistence of gas and liquid phase in the pore network. The phase distribution has a strong influence on transport processes and electrochemical reactions. In this article we present 2D and 3D RothmanKeller type multiphase LatticeBoltzmann models which take into account the heterogeneous wetting behavior of gas diffusion electrodes. The simulations are performed on FIBSEM 3D reconstructions of an Ag model electrode for predefined saturation of the pore space with the liquid phase. The resulting pressuresaturation characteristics and transport correlations are important input parameters for modeling approaches on the continuum scale and allow for an efficient development of improved gas diffusion electrodes.
keywords:
LatticeBoltzmann method, gas diffusion electrodes, FIBSEM tomography, metalair batteries, multiphase flowsquare, compress
1 Introduction
Metalair batteries possess a very high theoretical energy density which makes them interesting for both mobile and stationary applications Rahman:MetalAirReview (); Cheng2012:MetalAirReview (). At the negative electrode, metals like Al Egan:AluminumAirReview (), Li Girishkumar2010:LiAirPromiseChallenges (), Mg Zhang2015:MgAir (), Na Palomares:NaReview (); Slater:NaReview (), and Zn Lee2011:ZnAir (); Li2014:ZnAir () were suggested in the literature Kim:MetalElectrodesReview (). In recent years Liair batteries received the most attention in the battery community Christensen2012:LiO2Review (); ShaoHorn:BridgingMechanisticUnderstanding (). However, the only system which successfully reached the stage of mass production is the primary Znair battery. At the positive electrode oxygen is reduced and evolved during discharge and charge, respectively. Sufficient supply of O during discharge is accomplished by the concept of porous gas diffusion electrodes (GDEs). The electrodes are commonly made of carbon materials. However, in aqueous systems carbon is known to dissolve (’carbon corrosion’) and carbonfree GDEs were proposed Wittmaier2014:Ag (); Wittmaier2014:Ni (); Wittmaier2015:AgCo (). Hydrophobic binder materials ensure the coexistence of gas and liquid phase in the porous structure of the GDE. The saturation behavior is characteristic for the porous material and can be described by capillary pressure saturation () curves. The amount and distribution of the liquid phase has a strong influence on transport processes. The transport in the gas phase ensures a good supply of O and, thus, allows to draw high current densities. Moreover, the binder improves the mechanical stability of the electrode.
Due to their application in alkaline fuel cells the concept of gas diffusion electrodes was studied already in the 1960s. Design optimizations of the electrodes were mainly done by intensive experimental studies. In recent years the improvements in computational efficiency made computer simulations a common design tool in the engineering disciplines. However, traditional CFD (computational fluid dynamics) tools like the volume of fluid (VOF) method Hirt1981:VOF (); Ferreira2015:VOF () have their limitations in the simulation of multiphase flow in complex geometries. In recent years the lattice Boltzmann method (LBM) McNamara1988:LBM (); Sukop2006:LBM (); Succi2001:LBM () became increasingly popular for this class of problems because it is easy to implement and scales favorably. In LBM a probability distribution of discrete particle velocities is propagated on a computational lattice. Interactions between particles, boundaries, etc. are modeled by suitable collision operators. Several multiphase models were suggested in the literature for the simulation of immiscible fluids Sukop2015:MultiphaseLBM (). The most prominent ones are the ShanChen model Shan1993:SC1 (); Shan1993:SC2 (), the freeenergy model Swift1995:FE1 (); Swift1996:FE2 (), and the color gradient or RothmanKeller (RK) model Rothman1988:RK (); Gunstensen1991:RK1 (); Gunstensen1992:RK2 (). A common problem of the methods is the numerical stability and accuracy in the simulation of systems with high density and viscosity ratios which also includes the airwater system. Recent publications present modifications of the models which are able to overcome this limitation Inamuro2004:FE1 (); Zheng2006 (). This extends the applicability of the method to technical systems like gas diffusion media of polymer electrolyte fuel cells and the number of publications in the field increased rapidly Niu2007:LBM (); Park2008:LBM (); Mukherjee2011:Review (); Hao2012:CapillaryPressure (); Nabovati2014:PorosityEffect (); Dong2015:LBM_compression (). However, to our knowledge this is the first publication of porescale LBM simulations of gas diffusion electrodes for metalair batteries. In our work we use an RK type multiphase model to simulate multiphase flow in two and three dimensions. The model is based on the work of Leclaire Leclaire2011:ColorGradient (); Leclaire2012:Recoloring (); Leclaire2013:MultiPhase () and Liu et al. Liu2012:RK3D (). In their publications the authors successfully demonstrated the simulation of high density ratios. In our study we focus on aqueous electrolytes, however, the presented methodology is not limited to this case. An important feature of our model is that we explicitly take into account the nonhomogeneous wetting properties of the GDE which consists of hydrophilic electrode particles and hydrophobic binder fibers. This is an important step and has been barely pursued in previous porescale studies of electrochemical devices Hao2012:CapillaryPressure (); Molaeimanesh:PTFE (); Hao2010:Permeabilities (); Hao2010:WaterTransport (); Zhou2010:WaterTransport (); Mukherjee2009:WaterManagement ().
The focus of our article is set on the development of a methodology for the characterization of gas diffusion electrodes for metalair batteries. First, the structure of porous Ag model electrodes is reconstructed based on FIBSEM tomography as explained in Section 2. The reconstruction serves as simulation domain for multiphase LBM simulations. Model equations and computational details are summarized in Sections 3 and 4. 2D and 3D simulations are performed to simulate the evolution of the phase distribution towards an energetic minimum (Section 5.1). The results are evaluated to obtain characteristic pressuresaturation curves (Section 5.2) and saturation dependent transport parameters (Section 5.3).
2 Electrode reconstruction
We use samples of Ag model electrodes which were characterized regarding their electrochemical performance in our previous article Danner2014:Halfcell (). The focus of this work is set on the structural characterization and investigation of transport processes on the porescale. In this section we explain in detail the methodology which was developed for the reconstruction of gas diffusion electrodes using FIBSEM tomography. The suggested procedure is schematically shown in Figure 1.
2.1 FibSem
To obtain the microstructure for our simulations a dualbeam ZEISS NVision 40 SEMFIB instrument is used. The SEM has a fieldemission electron gun with an acceleration voltage between 5 and 30 kV and a vertical electronoptic axis. The FIB is based on a Ga primary ionbeam with a 30 kV acceleration voltage. The ionoptic axis is at an angle of 54° to the electronoptic axis. The standard vacuum levels for the electrongun and the sample chamber are 10 and 10 mbar, respectively. The serial sectioning of the sample is achieved using FIB and the images are acquired using the SEM at specified milling intervals.
Before the SEM imaging of the sample, it must be ensured that the SEM and FIB images correspond to the same regionofinterest of the sample. This condition is obtained in the following way: Firstly, we set the sample stage to the eucentric tilt position to avoid an offset in sample position with sample tilt. Secondly, the coincidence of the electron and ion beam has to be ensured. This is established by tilting the sample by 54° with a constant workingdistance of about 5 mm and adjusting the Zaxis until the SEM image of the sample comes into focus. A trench is cut into the sample in the vicinity of the regionofinterest with a relatively intense ionbeam current of 6.5 nA such that the crosssectional (CS) plane becomes visible. The CS plane is then gently polished with an ionbeam current of 300 pA. Subsequently, the total region for 3D tomography is selected for ionmilling. The two samples were filled with a low viscosity epoxy resin in order to improve the contrast of the images. This is an important step to facilitate the reconstruction process. The left panel of Figure 1 shows a representative SEM image. We note that the resin has rather well impregnated the pore space of the GDEs. Horizontal and vertical dimension are in the following named and , respectively. The direction perpendicular to the  plane is denoted by and represents the direction of the FIB cut. The total thickness of the cut is 10 µm. In order to optimize the total duration of milling and imageacquisition, the SEM images are obtained at every fifth FIBslice which results in 84 images representing slices of 0.12 µm thickness. Note, that (, ) values of the voxelsize are related by the relation = /sin(54°), where the value of is determined by the magnification of the image and the 54° origin from the tilt of the sample.
2.2 Structure generation
The second panel of Figure 1 shows a representative binarized and cropped image of the electrode microstructure (black). The images of different slices were aligned to account for the sampledrift during imaging. The perspective correction and pixelsize adjustments were done in the software package IMOD. The epoxy resin improves the contrast in the images and helps to identify solid particles. In spots where the impregnation of the electrode is incomplete the phases are assigned manually. Further details of the methodology of reconstruction is discussed in detail elsewhere Eswara2014:FIBSEM (). Finally, the images were stacked to a virtual structure in the commercial software GeoDict GeoDict:UserGuide (). The resulting geometry can be seen in the third panel of Figure 1.
The reconstructions are mirrored at the plane in order to increase the simulation domain in the direction of FIB sampling (direction). This step avoids systematic errors due to periodic boundary conditions which are used in the LBM simulations. The resulting geometry is subsequently coarsened in order to decrease the computational load. The dimensions of the samples are summarized in Table 1.
As outlined above the hydrophobic binder material is important for the distribution of the liquid phase and, thus, performance of the device. Unfortunately, the binder distribution is lost in the imaging and reconstruction process. In our reconstructions we distribute binder fibers by a random walk algorithm on the electrode surface. Binder fibers crossing the void space between electrode particles are neglected. The resulting binder distribution is shown in red color in the right panel of Figure 1. It is in reasonable optical agreement with distributions observed in SEM images Danner2014:Halfcell ().
3 Lattice Boltzmann Method
3.1 Model description
In this work we perform two (D2Q9) and three (D3Q19) dimensional simulations of twophase flow in porous gas diffusion electrodes. The LBM multiphase models employed in this study are based on the colorgradient approach and follow the recent publications of Liu and Leclaire et al. Leclaire2011:ColorGradient (); Leclaire2012:Recoloring (); Leclaire2013:MultiPhase (); Liu2012:RK3D (). A detailed derivation and discussion of the model equations can be found in their publications.
In RK type models the state of each phase (=gas, liquid) is described with the help of a probability distribution function , where is the index representing discrete directions in velocity space (see Eq. (27) and (28)). The macroscopic properties at a lattice node are given by the moments of the probability distribution. For instance, the density is calculated by the 0th moment which is simply the sum of the probabilities in all lattice directions
(1) 
The fluid velocity follows as the 1st moment which can be regarded as average of the discrete velocities weighted by their probability
(2) 
Finally, the pressure can be calculated according to
(3) 
where is a simulation parameter given in A.
In general the temporal evolution of the probability distribution is described by
(4) 
where the collision operator is in RK type models a result of three sub operators Tolke2002:RK (); Reis2007:RK ()
(5) 
The operators are applied successively to the probability distribution and are defined by

Singlephase collision (SRTBGK)
(6) In our model we apply the standard SRTBGK approximation Bhatnagar1954:BGKOperator () to model fluid interactions in the same phase
(7) where is the relaxation time of the collision process and the local equilibrium distribution of fluid velocities. The relaxation time can be related to the kinematic viscosity by
(8) In our simulations we use a single relaxation time for both phases which is calculated from a density average of the kinematic viscosity
(9) The local equilibrium distribution is given by the Maxwell distributions
(10) where (Eq. (31) and (32)) determines the compressibility of the fluid and is a lattice specific weighting parameter (see Eq. (29) and (30)).

Twophase collision (Perturbation)
(11) The effect of surface tension between the two phases is modeled by the perturbation operator. In order to distinguish the phases it is convenient to introduce a colorfield
(12) Regions of gas and liquid phase are marked by values of close to 1 and 1, respectively. The perturbation operator takes the form Reis2007:RK (); Liu2012:RK3D ()
(13) is a parameter controlling the surface tension (see Eq. (34)), is the color gradient in the twophase region, and is a parameter which has to be chosen in order to recover the NavierStokes Equations Reis2007:RK (); Liu2012:RK3D () (see Eq. (36) and (37)). The color gradient determining the surface normal in the twophase region is approximated by higherorder isotropic discretization schemes which significantly reduce spurious velocities Sbragaglia2007 (); Liu2012:RK3D (); Leclaire2011:ColorGradient ().

Twophase collision (Recoloring)
(14) The perturbation operator models the effect of surface tension, however, it does not enforce phase separation. This is ensured by the recoloring operator
(15) where is the total density, the total probability distribution, and the angle between the color gradient and the lattice direction vector
(16) The operator allows a moderate mixing of the two phases at the interface which additionally reduces spurious currents. The interface thickness is controlled by the parameter . Moreover, it was shown that the proposed operator circumvents the problem of ’latticepinning’ LatvaKokko2005:LatticePinning ().

Streaming
(17) Finally, the modified probability distributions are streamed to their new positions on the computational grid and the corresponding boundary conditions are applied (see Section 3.2).
3.2 Boundary conditions
There are two types of boundaries in microstructure resolved simulations of porous media: boundaries of the computational domain and boundaries at the fluidsolid interface. For the first type we employ simple periodic boundary conditions. This choice has implications on the simulation procedure and will be discussed in Section 4.1.
At the fluidsolid interface we use a simple bounce back scheme. The contact angle is adjusted by assigning an effective density to the probability distribution at solid nodes LatvaKokko2005:ContactAngle (). This has an effect on the local colorgradient and determines if the surface is wetting or nonwetting. In our model the densities on solid nodes were adjusted to match the contact angles which are observed experimentally (see Figure 2). Parameters can be found in Table 2.
3.3 Model parameterization
In this study we parameterize our models for the airwater system in order to represent metalair batteries with aqueous electrolytes. It has to be noted that other electrolyte systems using organic solvents or ionic liquids can be treated in the same framework. An analysis of dimensionless numbers is helpful to evaluate the relevant forces of the physical problem. Fluid flow in porous media is mainly determined by gravitational, viscous, inertial, and surface forces. The dimensionless numbers describing the ratio between these forces are

Bond number
(18) 
Capillary number
(19) 
Reynolds number
(20)
The characteristic length is in our studies given by the mean pore diameter µm of the electrode. The flow velocity can be approximated by the flow rate of volume replacement experiments and is typically less than the m s used in above calculations Dwenger2011:CapillaryPressure (). Moreover, in the limit of fluid mechanical equilibrium it will tend to zero. The small values of the dimensionless numbers demonstrate that gravitational and viscous forces are negligible compared to the strong influence of surface forces. Therefore, we use unity density and viscosity ratios in our simulations to improve the numerical stability. In a post processing step suitable scaling rules are applied to the simulation results Schaap2007:Rescaling (). The Laplace equation
(21) 
allows to deduce a relationship between the capillary pressure of the ’physical’ () and simulated system () according to
(22) 
where is one lattice unit, the size of one voxel and and the surface tension of the physical system and LBM simulation, respectively. The parameters of this study are summarized in Table 2.
3.4 Model validation
In order to validate our model we present two simple test cases.
Bubble test
The pressure difference across the interface of a steady bubble can be calculated by the Laplace equation (Eq. (21)). In order to validate our model we perform 2D simulations on a lattice with 100x100 nodes and varying bubble diameter. Figure 2 a) shows simulation results for diameters ranging from 10 to 90 lattice units. The size of one lattice unit corresponds to the voxel size of the reconstruction (see Table 1). The simulated pressures reproduce favorably the Laplace equation. At small diameters a minor deviation of simulation results is observed. This can be assigned to the increasing ratio of surface to bulk voxels which is important to notice for the calculation of curves.
Static contact angle
The second case is the simulation of static contact angles as presented in Figure 2 b). The fluid density of solid nodes was adjusted to reproduce the contact angle of water on Ag Zhao:ContactAngleAgPTFE () and PTFE surfaces Mukherjee2009:WaterManagement (). The Figure demonstrates that our model is able to simulate the wetting characteristics of the Ag GDE sample.
4 Simulation methodology
4.1 Initial conditions
The curves of porous media are commonly determined with the method of standard porosimetry (MSP) Gostick2006:CapillaryPressure (); Kumbur2007:LeverettIII () or dynamic volume replacement experiments Fairweather2007:CapillaryPressure (); Gostick2008:CapillaryPressure (); Nguyen2008:CapillaryPressure (); Dwenger2011:CapillaryPressure (). After filling the porous medium to a defined saturation level the corresponding average capillary pressure follows as difference between the pressure in the gas and the liquid phase. It was found that the direction of the process (injection/imbibition or removal/drainage of the liquid phase) causes a hysteresis in the resulting curves. It is therefore expected that a different protocol for loading the porous medium with liquid or gas phase will lead to slightly different curves. Since the conditions for the standard dynamic volume experiments do not match the conditions in the GDE under dynamic operations, we chose a different simulation set up to determine curves. This setup corresponds to establishing a force balance in the volume of the porous medium after initializing the simulation with a maximally wetting or nonwetting condition respectively for a given volume fraction of liquid in the GDE. We argue, that the thus obtained pressure saturation curves are closer to the real operating condition of a GDE in a metalair cell, where electrolyte is pushed out of the GDE due to the occurrence of solid reaction products in the bulk of the porous medium and not due to the application of external pressure forces Danner2014:Halfcell (). The different initial configurations within our LBM model with periodic boundary conditions are shown in the first row of Figure 3. During drainage the electrolyte will remain preferentially on the hydrophilic electrode particles. We maximize the interface between liquid phase and hydrophilic Ag particles by randomly placing liquid droplets on the electrode surface or, if all surface sites are occupied, in contact with another liquid droplet (configuration I). In the second configuration the liquid phase is introduced in the pore space as one solid block (Figure 3 right column  configuration II) in order to minimize the interfacial area between gas and liquid phase. This situation is comparable to the imbibition process in the experimental setup where the electrolyte is forced into the porous structure and also occupies areas which are energetically not favorable. Our simulations are able to reproduce the hysteretic behavior which is also observed in the forced drainageimbibition experiments. These results will be discussed in more detail in Section 5.2.
4.2 Pressuresaturation curves
For the determination of curves we perform independent 2D and 3D simulations at various saturation levels of the GDE sample. The number of iterations was chosen sufficiently large to achieve fluid mechanical equilibrium. At the end of the simulation the capillary pressure follows as , where is averaged over the simulation domain. In the case of 2D simulations we conduct at each saturation 10 independent simulations on randomly chosen twodimensional slices of the reconstruction in order to get statistically significant results.
Pressure saturation curves are commonly described with the so called Leverett function Leverett1940:J (); Udell1985:J ()
(23) 
where is the permeability and the porosity of the electrode. In this study we use a modification of the standard formulation according to Hao et al. Hao2012:CapillaryPressure () to account for the heterogeneous wetting properties of the GDE
(24) 
4.3 Transport parameters
The final phase distribution at the end of the 3D simulations is used as input for the calculation of saturationdependent effective transport parameters and surface areas in the software GeoDict GeoDict:UserGuide (). The final fluid distribution is basically regarded as ’frozen’ in the porous structure. For the determination of effective transport parameters in the gas phase the voxels of the liquid phase are then treated as solids and vice versa. In their study Garcia et al. Garcia2015:Diffusivity () took a similar approach to determine the diffusivity of partially saturated GDLs. However, they obtained the distribution of the liquid phase from tomographic images instead of twophase simulations.
In the modeling of electrochemical devices the Bruggeman correlation is often used for the approximation of effective transport parameters. The diffusivity which is the ratio between effective and bulk diffusion coefficients is in this approach described by
(25) 
where is the volume fraction of the transporting phase and the socalled Bruggeman coefficient. The volume fraction of the liquid and gas phase are calculated according to and , respectively.
5 Results and discussion
5.1 Electrolyte distribution
Figure 3 shows the phase distribution during 3D simulations of configuration I (left) and configuration II (right) at a saturation of the pore space with liquid electrolyte of 50%. The first row illustrates the initial conditions which were applied to mimic the drainage and imbibition of the liquid phase. (cf. Section 4.1). In configuration I the electrolyte is initially distributed on the electrode surface. During the simulation the distribution evolves towards a local minimum of the free energy. Our simulations reach a quasi stationary state after about 110 iterations. Further changes are only marginal. This is also reflected in the corresponding pressure signal which is shown as function of iterations in Figure 4 b). A visually similar steadystate is also observed for configuration II. This indicates that the initial condition which we apply in our simulations has only a minor influence on the saturation dependence of effective transport parameters. This effect will be discussed in Section 5.3.
5.2 Pressuresaturation curves
Figure 4 shows the capillary pressure signal during 2D and 3D simulations representing the drainage process. In both cases the pressure signal approaches a constant value at the end of the simulations which indicates that the system is in a local minimum of the free energy. The 3D simulations give a clear trend of increasing capillary pressure with increasing liquid phase saturation. It has to be noted that the 2D simulations shown in Figure 4 a) were not performed on the same slice of the reconstruction. Therefore, the fluctuations in pressure which occur are due to different simulation domains and a clear trend can not be deduced from the graph. In our approach we take the average of 10 simulations to capture this effect statistically. Figure 5 shows curves for configuration I (left) and configuration II (right). The results of the 2D simulations are shown in red color and the error bars represent the corresponding standard deviation. The error bars are larger at low and high saturation which can be explained by a reduced configurational freedom. Meaning, reconfigurations of the phases are suppressed and they remain in their initial (random) configuration. The calculated pressures of 3D simulations are included as blue triangles. There is a systematic deviation at high saturation, however, the general agreement between 2D and 3D simulations is favorable. This is an important result which shows that a series of computationally efficient 2D simulations can be used for a screening of different electrode structures. In a previous publication Danner2014:Halfcell () we predicted a liquid phase saturation of 50 % at by a fit of continuum simulations to electrochemical data. This is in excellent agreement with the results presented in this work and indicates the validity of our approach.
Pressure saturation curves are commonly presented as dimensionless Leverett functions (see Eq. (24)). Hao et al. showed that an expression of the form
(26) 
is able to give a good representation of data. We fit the parameters of Eq. (26) to the curves of our 2D simulations. Figure 6 shows Leverett functions of configuration I (left) and configuration II (right). Similar to the experimental work on GDLs we observe a hysteresis at intermediate saturation Fairweather2007:CapillaryPressure (); Dwenger2011:CapillaryPressure (). This can be explained by looking at the experimental procedure and taking into account microstructural effects in the GDE. Some areas are connected to the pore network only through narrow pores. According to the Laplace equation they are filled with the liquid phase only at high capillary pressures. On the reverse process a higher negative pressure difference is needed to withdraw the electrolyte from this part of the pore network. A similar reasoning can be made for the simulation approach presented in this study. In the simulations starting with configuration I all parts of the pore network are accessible for the liquid phase. This results in a lower capillary pressure compared to simulations with initial configuration II where the access is restricted to surrounding pores. This result justifies our choice of initial conditions to mimic the two processes. However, it has to be noted that the hysteresis is within the standard deviation of the simulations. Experimental studies on the electrodes will be needed to validate our simulation methodology. The parameters of Eq. (26) resulting from the fit to the simulation data are summarized in Table 3.
5.3 Effective transport parameters
Figure 7 shows the diffusivity in the gas and liquid phase for configuration I (left) and configuration II (right). Results are presented individually for the three spatial coordinates , , and . The graphs show that the transport behavior is anisotropic. Especially, in direction which is the direction of FIB sampling the transport parameters are smaller. This can be attributed to a lower spatial resolution and indicates that the anisotropy might be an artifact of the reconstruction. The overall trend of the diffusivity with saturation is the same for both initial configurations. The diffusivity of the gas phase decreases with an increasing saturation of the pore space with the liquid phase. This corresponds to the behavior predicted by the Bruggeman correlation (cf. Eq. (25)) presented in Figure 8 a).
The main transport direction in the GDE during operation is in our notation denoted by and we limit our discussion in the remainder of this section to this most important case. Figure 8 a) shows the calculated diffusivity in the gas and liquid phase for simulations with initial configuration I and II. As mentioned in Section 5.1 the final liquid phase distribution is similar for both cases and, thus, also saturation dependent transport coefficients are comparable. The dotted line represents the Bruggeman correlation with a coefficient of =1.7. The graph demonstrates that the correlation gives a reasonable representation of the simulated values although the coefficient is slightly higher than the standard value of 1.5 which is usually assumed in the literature.
An important parameter for the performance of the GDE is the length of the triple phase boundary (TPB). Here, gas, liquid, and solid phase are in close contact and reaction kinetics are facile. The dependence of the TPB on the saturation is shown in Figure 8 b) for both, configuration I and configuration II. The length of the TPB first increases with saturation. Then, it reaches a plateau because large parts of the pore space are completely filled with electrolyte. The same effect leads to a decrease in TPB length at high saturation. This shows that at a high to moderate saturation the performance of the GDE is improved. The deviation between configuration I and configuration II origins from the more homogeneous distribution of the liquid phase in configuration I.
The correlations determined above are important input for simulations on the continuum scale Horstmann2013:Precipitation (); Danner2014:Halfcell () and allow an efficient improvement of gas diffusion electrodes. An overview of this multiscale approach can also be found in Ref. Gruebl2014:Methodology ().
6 Conclusions
In this article we present a methodology for the characterization of gas diffusion electrodes for metalair batteries. First, we take FIBSEM images of an Ag model electrode. In a subsequent step the images are binarized and stacked to a virtual reconstruction of the electrode. On this geometry we performed 2D and 3D LBM simulations with a RK type multiphase model. To the best of our knowledge this study is the first publication investigating the complex multiphase transport mechanism in GDEs using LBM on the porescale. An important feature of our model is that we take into account the nonhomogeneous wetting behavior of hydrophilic Ag particles and hydrophobic binder. The model is parametrized to represent an aqueous electrolyte system, however, other electrolytes (e.g. organic solutions or ionic liquids) can be treated in the same framework. Results of the simulations are evaluated for the determination of curves and saturation dependent transport parameters. Our simulations demonstrate that a series of 2D simulations is an efficient tool for the screening of characteristics of GDEs. Transport parameters and specific surface areas are determined from the final phase distributions of 3D simulations. The diffusivity generally follows the trends predicted by the Bruggeman correlation with a coefficient of 1.7. The results of this study provide important input parameters for simulations on the continuum scale. Moreover, the methodology presented here can be used as a tool for the optimization of GDEs for metal air batteries.
Acknowledgments Funding from the European Union’s innovation programme Horizon 2020 (20142020) under grant agreement n°646186 (ZAS) is gratefully acknowledged. The authors would also like to thank Norbert Wagner, Institute for Engineering Thermodynamics, DLR Stuttgart for providing the samples of Ag model electrodes and Prasanth Balasubramanian, Ulm University for help in tomogram reconstruction. Santhana Eswara would like to thank Jörg Bernhard and Prof. Ute Kaiser of Ulm University for fruitful discussions. Timo Danner would like to thank Prof. Wolfgang Bessler, University of Applied Sciences Offenburg, for discussions and comments.
References
 (1) M. A. Rahman, X. Wang, C. Wen, High Energy Density MetalAir Batteries: A Review, J. Electrochem. Soc. 160 (10) (2013) A1759–A1771. doi:10.1149/2.062310jes.
 (2) F. Cheng, J. Chen, Metalair batteries: from oxygen reduction electrochemistry to cathode catalysts, Chem. Soc. Rev. 41 (2012) 2172–2192. doi:10.1039/C1CS15228A.
 (3) D. Egan, C. Ponce de León, R. Wood, R. Jones, K. Stokes, F. Walsh, Developments in electrode materials and electrolytes for aluminiumair batteries, J. Power Sources 236 (2013) 293–310. doi:10.1016/j.jpowsour.2013.01.141.
 (4) G. Girishkumar, B. McCloskey, A. C. Luntz, S. Swanson, W. Wilcke, Lithium  Air Battery: Promise and Challenges, J. Phys. Chem. Lett. 1 (14) (2010) 2193–2203. doi:10.1021/jz1005384.
 (5) T. Zhang, Z. Tao, J. Chen, Magnesiumair batteries: from principle to application, Mater. Horiz. 1 (2014) 196–206. doi:10.1039/C3MH00059A.
 (6) V. Palomares, P. Serras, I. Villaluenga, K. B. Hueso, J. CarreteroGonzález, T. Rojo, Naion batteries, recent advances and present challenges to become low cost energy storage systems, Energy Environ. Sci. 5 (3) (2012) 5884. doi:10.1039/c2ee02781j.
 (7) M. D. Slater, D. Kim, E. Lee, C. S. Johnson, SodiumIon Batteries, Adv. Func. Mater. 23 (8) (2013) 947–958. doi:10.1002/adfm.201200691.
 (8) J.S. Lee, S. Tai Kim, R. Cao, N.S. Choi, M. Liu, K. T. Lee, J. Cho, Metalair batteries with high energy density: Liair versus znair, Adv. Energy Mater. 1 (1) (2011) 34–50. doi:10.1002/aenm.201000010.
 (9) Y. Li, H. Dai, Recent advances in zincair batteries, Chem. Soc. Rev. 43 (15) (2014) 5257–5275. doi:10.1039/C4CS00015C.
 (10) H. Kim, G. Jeong, Y.U. Kim, J.H. Kim, C.M. Park, H.J. Sohn, Metallic anodes for next generation secondary batteries., Chem. Soc. Rev. 42 (23) (2013) 9011–34. doi:10.1039/c3cs60177c.
 (11) J. Christensen, P. Albertus, R. S. SanchezCarrera, T. Lohmann, B. Kozinsky, R. Liedtke, J. Ahmed, A. Kojic, A Critical Review of Li/Air Batteries, J. Electrochem. Soc. 159 (2) (2012) R1–R30. doi:10.1149/2.086202jes.
 (12) Y.C. Lu, B. M. Gallant, D. G. Kwabi, J. R. Harding, R. R. Mitchell, M. S. Whittingham, Y. ShaoHorn, Lithiumoxygen batteries: bridging mechanistic understanding and battery performance, Energy Environ. Sci. 6 (2013) 750–768. doi:10.1039/C3EE23966G.
 (13) D. Wittmaier, N. Wagner, K. A. Friedrich, H. M. Amin, H. Baltruschat, Modified carbonfree silver electrodes for the use as cathodes in lithiumair batteries with an aqueous alkaline electrolyte, J. Power Sources 265 (2014) 299 – 308. doi:10.1016/j.jpowsour.2014.04.142.
 (14) D. Wittmaier, S. Aisenbrey, N. Wagner, K. A. Friedrich, Bifunctional, carbonfree nickel/cobaltoxide cathodes for lithiumair batteries with an aqueous alkaline electrolyte, Electrochim. Acta 149 (2014) 355 – 363. doi:10.1016/j.electacta.2014.10.088.
 (15) D. Wittmaier, N. A. Cañas, I. Biswas, K. A. Friedrich, Highly Stable CarbonFree Ag/CoOCathodes for LithiumAir Batteries: Electrochemical and Structural Investigations, Adv. Energy Mater.doi:10.1002/aenm.201500763.
 (16) C. Hirt, B. Nichols, Volume of fluid (vof) method for the dynamics of free boundaries, J. Comput. Phys. 39 (1) (1981) 201 – 225. doi:10.1016/00219991(81)901455.
 (17) R. B. Ferreira, D. S. Falcao, V. B. Oliveira, A. M. F. R. Pinto, Numerical simulations of twophase flow in proton exchange membrane fuel cells using the volume of fluid method  A review, J. Power Sources 277 (2015) 329–342. doi:10.1016/j.jpowsour.2014.11.124.
 (18) G. R. McNamara, G. Zanetti, Use of the boltzmann equation to simulate latticegas automata, Phys. Rev. Lett. 61 (1988) 2332–2335. doi:10.1103/PhysRevLett.61.2332.
 (19) M. C. Sukop, D. T. Thorne, Lattice Boltzmann Modeling: An Introduction for Geoscientists and Engineers, 1st Edition, Springer Berlin Heidelberg, 2006.
 (20) S. Succi, The Lattice Boltzmann Equation: For Fluid Dynamics and Beyond, Numerical Mathematics and Scientific Computation, Clarendon Press, 2001.
 (21) H. Huang, M. Sukop, X. Lu, Multiphase Lattice Boltzmann Methods: Theory and Application, Wiley, 2015.
 (22) X. Shan, H. Chen, Lattice boltzmann model for simulating flows with multiple phases and components, Phys. Rev. E 47 (1993) 1815–1819. doi:10.1103/PhysRevE.47.1815.
 (23) X. Shan, H. Chen, Simulation of nonideal gases and liquidgas phase transitions by the lattice boltzmann equation, Phys. Rev. E 49 (1994) 2941–2948. doi:10.1103/PhysRevE.49.2941.
 (24) M. R. Swift, W. R. Osborn, J. M. Yeomans, Lattice boltzmann simulation of nonideal fluids, Phys. Rev. Lett. 75 (1995) 830–833. doi:10.1103/PhysRevLett.75.830.
 (25) M. R. Swift, E. Orlandini, W. R. Osborn, J. M. Yeomans, Lattice boltzmann simulations of liquidgas and binary fluid systems, Phys. Rev. E 54 (1996) 5041–5052. doi:10.1103/PhysRevE.54.5041.
 (26) D. H. Rothman, J. M. Keller, Immiscible cellularautomaton fluids, J. Stat. Phys. 52 (34) (1988) 1119–1127. doi:10.1007/BF01019743.
 (27) A. K. Gunstensen, D. H. Rothman, S. Zaleski, G. Zanetti, Lattice boltzmann model of immiscible fluids, Phys. Rev. A 43 (1991) 4320–4327. doi:10.1103/PhysRevA.43.4320.
 (28) A. K. Gunstensen, D. H. Rothman, Microscopic modeling of immiscible fluids in three dimensions by a lattice boltzmann method, Europhys. Lett. 18 (2) (1992) 157.
 (29) T. Inamuro, T. Ogata, S. Tajima, N. Konishi, A lattice Boltzmann method for incompressible twophase flows with large density differences, J. Comput. Phys. 198 (2) (2004) 628–644. doi:10.1016/j.jcp.2004.01.019.
 (30) H. Zheng, C. Shu, Y. Chew, A lattice Boltzmann model for multiphase flows with large density ratio, J. Comput. Phys. 218 (1) (2006) 353–371. doi:10.1016/j.jcp.2006.02.015.
 (31) X. Niu, T. Munekata, S. Hyodo, K. Suga, An investigation of watergas transport processes in the gasdiffusionlayer of a PEM fuel cell by a multiphase multiplerelaxationtime lattice boltzmann model, J. Power Sources 172 (2) (2007) 542 – 552. doi:10.1016/j.jpowsour.2007.05.081.
 (32) J. Park, X. Li, Multiphase microscale flow simulation in the electrodes of a PEM fuel cell by lattice Boltzmann method, J. Power Sources 178 (1) (2008) 248–257. doi:10.1016/j.jpowsour.2007.12.008.
 (33) P. P. Mukherjee, Q. Kang, C.Y. Wang, Porescale modeling of twophase transport in polymer electrolyte fuel cellsprogress and perspective, Energy Environ. Sci. 4 (2) (2011) 346–369. doi:10.1039/b926077c.
 (34) L. Hao, P. Cheng, Capillary pressures in carbon paper gas diffusion layers having hydrophilic and hydrophobic pores, Int. J. Heat Mass Trans. 55 (13) (2012) 133–139. doi:10.1016/j.ijheatmasstransfer.2011.08.049.
 (35) A. Nabovati, J. Hinebaugh, A. Bazylak, C. H. Amon, Effect of porosity heterogeneity on the permeability and tortuosity of gas diffusion layers in polymer electrolyte membrane fuel cells, J. Power Sources 248 (2014) 83–90. doi:10.1016/j.jpowsour.2013.09.061.
 (36) D. H. Jeon, H. Kim, Effect of compression on water transport in gas diffusion layer of polymer electrolyte membrane fuel cell using lattice Boltzmann method, J. Power Sources 294 (2015) 393–405. doi:10.1016/j.jpowsour.2015.06.080.
 (37) S. Leclaire, M. Reggio, J.Y. Trépanier, Isotropic color gradient for simulating very highdensity ratios with a twophase flow lattice Boltzmann model, Comput. Fluids 48 (1) (2011) 98–112. doi:10.1016/j.compfluid.2011.04.001.
 (38) S. Leclaire, M. Reggio, J.Y. Trépanier, Numerical evaluation of two recoloring operators for an immiscible twophase flow lattice Boltzmann model, Appl. Math. Model. 36 (5) (2012) 2237 – 2252. doi:10.1016/j.apm.2011.08.027.
 (39) S. Leclaire, M. Reggio, J.Y. Trépanier, Progress and investigation on lattice Boltzmann modeling of multiple immiscible fluids or components with variable density and viscosity ratios, J. Comput. Phys. 246 (2013) 318–342. doi:10.1016/j.jcp.2013.03.039.
 (40) H. Liu, A. J. Valocchi, Q. Kang, Threedimensional lattice Boltzmann model for immiscible twophase flow simulations, Phys. Rev. E 85 (4) (2012) 046309. doi:10.1103/PhysRevE.85.046309.
 (41) G. R. Molaeimanesh, M. H. Akbari, Impact of PTFE distribution on the removal of liquid water from a PEMFC electrode by lattice Boltzmann method, Int. J. Hydrogen Energy 39 (16) (2014) 8401–8409. doi:10.1016/j.ijhydene.2014.03.089.
 (42) L. Hao, P. Cheng, Porescale simulations on relative permeabilities of porous media by lattice Boltzmann method, Int. J. Heat Mass Tran. 53 (910) (2010) 1908–1913. doi:10.1016/j.ijheatmasstransfer.2009.12.066.
 (43) L. Hao, P. Cheng, Lattice Boltzmann simulations of water transport in gas diffusion layer of a polymer electrolyte membrane fuel cell, J. Power Sources 195 (12, SI) (2010) 3870–3881. doi:10.1016/j.jpowsour.2009.11.125.
 (44) P. Zhou, C. W. Wu, Liquid water transport mechanism in the gas diffusion layer, J. Power Sources 195 (5) (2010) 1408–1415. doi:10.1016/j.jpowsour.2009.09.019.
 (45) P. P. Mukherjee, C.Y. Wang, Q. Kang, Mesoscopic modeling of twophase behavior and flooding phenomena in polymer electrolyte fuel cells, Electrochim. Acta 54 (27) (2009) 6861–6875. doi:10.1016/j.electacta.2009.06.066.
 (46) T. Danner, B. Horstmann, D. Wittmaier, N. Wagner, W. G. Bessler, Reaction and transport in Ag/AgO gas diffusion electrodes of aqueous LiO batteries: Experiments and modeling, J. Power Sources 264 (0) (2014) 320 – 332. doi:10.1016/j.jpowsour.2014.03.149.
 (47) S. K. EswaraMoorthy, P. Balasubramanian, W. van Mierlo, J. Bernhard, M. Marinaro, M. WohlfahrtMehrens, L. Jörissen, U. Kaiser, An in situ semfibbased method for contrast enhancement and tomographic reconstruction for structural quantification of porous carbon electrodes, Microsc. Microanal. 20 (2014) 1576–1580. doi:10.1017/S1431927614012884.

(48)
Math2Market GmbH  Becker, J.; Wiegmann, A.,
GeoDict (2013).
URL http://www.geodict.com  (49) J. Tölke, M. Krafczyk, M. Schulz, E. Rank, Lattice Boltzmann simulations of binary fluid flow through porous media., Philos. T. Roy. Soc. A 360 (1792) (2002) 535–45. doi:10.1098/rsta.2001.0944.
 (50) T. Reis, T. N. Phillips, Lattice Boltzmann model for simulating immiscible twophase flows, J. Phys. AMath. Theor. 40 (14) (2007) 4033–4053. doi:10.1088/17518113/40/14/018.
 (51) P. L. Bhatnagar, E. P. Gross, M. Krook, A model for collision processes in gases. i. small amplitude processes in charged and neutral onecomponent systems, Phys. Rev. 94 (1954) 511–525. doi:10.1103/PhysRev.94.511.
 (52) M. Sbragaglia, R. Benzi, L. Biferale, S. Succi, K. Sugiyama, F. Toschi, Generalized lattice Boltzmann method with multirange pseudopotential, Phys. Rev. E 75 (2) (2007) 026702. doi:10.1103/PhysRevE.75.026702.
 (53) M. LatvaKokko, D. Rothman, Diffusion properties of gradientbased lattice Boltzmann models of immiscible fluids, Phys. Rev. E 71 (5) (2005) 056702. doi:10.1103/PhysRevE.71.056702.
 (54) M. LatvaKokko, D. Rothman, Static contact angle in lattice Boltzmann models of immiscible fluids, Phys. Rev. E 72 (4) (2005) 046701. doi:10.1103/PhysRevE.72.046701.
 (55) S. Dwenger, G. Eigenberger, U. Nieken, Measurement of Capillary Pressure–Saturation Relationships Under Defined Compression Levels for Gas Diffusion Media of PEM Fuel Cells, Transport Porous Med. 91 (1) (2011) 281–294. doi:10.1007/s1124201198444.
 (56) M. G. Schaap, M. L. Porter, B. S. B. Christensen, D. Wildenschild, Comparison of pressuresaturation characteristics derived from computed tomography and lattice Boltzmann simulations, Water Resour. Res. 43 (12) (2007) W12S06. doi:10.1029/2006WR005730.
 (57) Q. Zhao, Y. Liu, C. Wang, Development and evaluation of electroless AgPTFE composite coatings with antimicrobial and anticorrosion properties, Appl. Surf. Sci. 252 (5) (2005) 1620–1627. doi:10.1016/j.apsusc.2005.02.098.
 (58) J. T. Gostick, M. W. Fowler, M. A. Ioannidis, M. D. Pritzker, Y. Volfkovich, A. Sakars, Capillary pressure and hydrophilic porosity in gas diffusion layers for polymer electrolyte fuel cells, J. Power Sources 156 (2) (2006) 375–387. doi:10.1016/j.jpowsour.2005.05.086.
 (59) E. C. Kumbur, K. V. Sharp, M. M. Mench, Validated Leverett Approach for Multiphase Flow in PEFC Diffusion Media  III. Temperature Effect and Unified Approach, J. Electrochem. Soc. 154 (12) (2007) B1315 – B1324. doi:10.1149/1.2784286.
 (60) J. D. Fairweather, P. Cheung, J. StPierre, D. T. Schwartz, A microfluidic approach for measuring capillary pressure in PEMFC gas diffusion layers, Electrochem. Commun. 9 (9) (2007) 2340–2345. doi:10.1016/j.elecom.2007.06.042.
 (61) J. T. Gostick, M. A. Ioannidis, M. W. Fowler, M. D. Pritzker, Direct measurement of the capillary pressure characteristics of waterairgas diffusion layer systems for pem fuel cells, Electrochem. Commun. 10 (10) (2008) 1520 – 1523. doi:10.1016/j.elecom.2008.08.008.
 (62) T. V. Nguyen, G. Lin, H. Ohn, X. Wang, Measurement of Capillary Pressure Property of Gas Diffusion Media Used in Proton Exchange Membrane Fuel Cells, Electrochem. Solid St. 11 (8) (2008) B127–B131. doi:10.1149/1.2929063.
 (63) M. Leverett, Capillary behavior in porous solids, T. Am. I. Min. Met. Eng. 142 (1941) 152–169. doi:10.2118/941152G.
 (64) K. S. Udell, Heat transfer in porous media considering phase change and capillaritythe heat pipe effect, Int. J. Heat Mass Tran. 28 (2) (1985) 485 – 495. doi:10.1016/00179310(85)900821.
 (65) P. A. GarciaSalaberri, J. T. Gostick, G. Hwang, A. Z. Weber, M. Vera, Effective diffusivity in partiallysaturated carbonfiber gas diffusion layers: Effect of local saturation and application to macroscopic continuum models, J. Power Sources 296 (2015) 440 – 453. doi:10.1016/j.jpowsour.2015.07.034.
 (66) B. Horstmann, T. Danner, W. G. Bessler, Precipitation in aqueous lithium–oxygen batteries: a modelbased analysis, Energy Environ. Sci. 6 (2013) 1299–1314. doi:10.1039/c3ee24299d.
 (67) D. Grübl, T. Danner, V. P. Schulz, A. Latz, W. G. Bessler, Multimethodology modeling and design of lithiumair cells with aqueous electrolyte, ECS Trans. 62 (1) (2014) 137–149. doi:10.1149/06201.0137ecst.
 (68) Y. Qian, D. D’Humières, P. Lallemand, Lattice BGK Models for NavierStokes Equation , Europhys. Lett. 17 (6BIS) (1992) 479–484. doi:10.1209/02955075/17/6/001.
List of Tables
 1 Dimensions and resolution of the electrode sample in direction. The reconstructed electrodes are based on cubic voxels. For the LBM simulations the structure is mirrored in direction to increase the computational domain.
 2 Physical parameters of the system (SI units) and corresponding input parameters of the LBM simulations (lattice units).
 3 Parameters of the Leverett function resulting from a fit to Eq. (26).
Dimensions / µm  Voxels /   Voxel size / nm  
x/y/z  x/y/z  x/y/z  
FIBSEM  27.8/23.9/10.0  504/346/84  55.25/69.1l/119 
Reconstruction  27.8/23.9/10.0  504/432/180  55.25 
LBM  27.8/23.9/20.0  126/108/90  221 
Gas  Liquid  Gas  Liquid  GasLiquid  Ag  Binder  
Physical  1.18  997  1.58  8.93  7.28  67 °Zhao:ContactAngleAgPTFE ()  140 °Mukherjee2009:WaterManagement () 
LBM  1  1  1/6  1/6  0.1  1.195/0.805†  0.357/1.643† 
† Density on wall nodes (gas/liquid) mu lu 
configuration I  0.813  0.233  5.079  2.32310  13.47 

configuration II  4.55710  3.79710  8.138  0.01  10.98 
List of Figures
 1 Methodology to reconstruct the LBM simulation domain from FIBSEM images. Full details regarding the dimensions and voxel sizes are given in Table 1. Scale bar of SEM image is 1 µm.
 2 Validation simulations of the LBM model. Left: Bubble test. Right: Influence of nonwetting and wetting surfaces.
 3 Electrolyte distribution in the porous GDE for configuration I (left) and configuration II (right).
 4 Pressure as function of iterations during simulations employing configuration I. Left: 2D simulations. Right: 3D simulations.
 5 Pressuresaturation curves during for configuration I (left) and configuration II (right). Results of the 2D simulations are shown in red color and error bars represent the standard deviation of the simulations. Results of the 3D simulations are represented by blue triangles.
 6 Plot of the dimensionless Leverett Function over saturation for configuration I (blue) and configuration II (red). Dashed lines are a result of the fit to Eq. (26). Symbols represent corresponding simulation results.
 7 Diffusivity in the gas (open circles) and liquid (open triangles) phase. Values are calculated based on the final electrolyte distribution obtained in the 3D LBM simulations. Plots for the three spatial coordinates (,,) are shown separately for configuration I (left) and configuration II (right).
 8 Diffusivity in the main transport direction (left) and specific length of the triple phase boundary (right) as function of the saturation for configuration I (blue) and configuration II (red). The Bruggeman correlation (Eq. (25), =1.70) is included as gray line in Figure 8 a).
 1 Methodology to reconstruct the LBM simulation domain from FIBSEM images. Full details regarding the dimensions and voxel sizes are given in Table 1. Scale bar of SEM image is 1 µm.
 2 Validation simulations of the LBM model. Left: Bubble test. Right: Influence of nonwetting and wetting surfaces.
 3 Electrolyte distribution in the porous GDE for configuration I (left) and configuration II (right).
 4 Pressure as function of iterations during simulations employing configuration I. Left: 2D simulations. Right: 3D simulations.
 5 Pressuresaturation curves during for configuration I (left) and configuration II (right). Results of the 2D simulations are shown in red color and error bars represent the standard deviation of the simulations. Results of the 3D simulations are represented by blue triangles.
 6 Plot of the dimensionless Leverett Function over saturation for configuration I (blue) and configuration II (red). Dashed lines are a result of the fit to Eq. (26). Symbols represent corresponding simulation results.
 7 Diffusivity in the gas (open circles) and liquid (open triangles) phase. Values are calculated based on the final electrolyte distribution obtained in the 3D LBM simulations. Plots for the three spatial coordinates (,,) are shown separately for configuration I (left) and configuration II (right).
 8 Diffusivity in the main transport direction (left) and specific length of the triple phase boundary (right) as function of the saturation for configuration I (blue) and configuration II (red). The Bruggeman correlation (Eq. (25), =1.70) is included as gray line in Figure 8 a).
Appendix A LBM parameters
a.1 Computational lattice
The D2Q9 and D3Q19 lattices Qian1992:Nomencalture () in cartesian coordinates are given by
(27) 
and
(28) 
respectively.
a.2 Singlephase collision
For the D2Q9 and D3Q19 lattices is given by
(29) 
and
(30) 
respectively.
The parameter is related to the compressibility of the fluid and, thus, to the speed of sound and hydrostatic pressure in phase . For the D2Q9 and D3Q19 lattice is given by Leclaire2013:MultiPhase (); Liu2012:RK3D ()
(31) 
and
(32) 
respectively, where one of the is a free parameter setting the pressure level in the system. The values of are related by
(33) 
in order to guarantee a stable interface (). In our simulations we set to 4/9.
a.3 Twophase collision: Perturbation
The parameter controlling the surface tension is defined as Liu2012:RK3D ()
(34) 
where is the surface tension and the density averaged relaxation time. In our simulations we use a fourthorder isotropic approximation for the calculation of the color gradient Liu2012:RK3D ()
(35) 
In the literature approximations of higher order are reported Sbragaglia2007 (), however, the computational load increases significantly.
The parameters depend on the lattice and have to be chosen in order to ensure the conservation of mass in the perturbation step. For the D2Q9 and D3Q19 lattice the parameters are given by Reis2007:RK (); Liu2012:RK3D ()
(36) 
and
(37) 
respectively.
a.4 Twophase collision: Recoloring
The only additional parameter appearing in the recoloring operator is which determines the thickness of the interfacial layer.