General Local Reactive Boundary Condition for Dissolution and Precipitation Using the Lattice Boltzmann Method

A general and local reactive boundary condition (RBC) for studying first‐order equilibrium reactions using the lattice Boltzmann method is presented. Its main characteristics are accurate reproduction of wall diffusion, invariance to the wall and grid orientation, and absence of nonphysical artifacts. The scheme is successfully tested for different benchmark cases considering diffusion, advection, and reactions of fluids at solid‐liquid interfaces. Unlike other comparable RBCs from the literature, the novel scheme is valid for a large range of Péclet and Damköhler numbers, and shows realistic pattern formation during precipitation. In addition, quantitative results are in good accordance with analytical solutions and values from literature. Combining the new RBC with the rest fraction method, Péclet‐Reynolds ratios of up to 1,000 can be achieved. Overall, the novel RBC accurately models first‐order reactions, is applicable for complex geometries, and allows efficiently simulating dissolution and precipitation phenomena in fluids at the pore scale.


Introduction
Dissolution and precipitation of solid solutes occur in a broad range of natural phenomena and technological applications.Typical examples are dissolution of minerals in subsurface hydrology (Přikryl et al., 2017;Baqer & Chen, 2022;Andrews et al., 2023), biofilm growth in nutrient rich environments (Jung & Meile, 2021), etching into substrates (Cui et al., 2019), and conversion of active material in energy storage systems (Fang et al., 2021;Danner & Latz, 2019).
Despite their omnipresence, a detailed experimental analysis of such processes is often lacking due to the dominant small length and short time scales on which they occur.As such, computational methods and more specifically continuum approaches are used to describe these phenomena at large.But even they are stretched to their limits, when dissolution and precipitation occur in structurally complex porous media where the relevant physics happen at the pore scale.
A computational method that has proven to give insight into such mesoscopic phenomena is the lattice Boltzmann method (LBM).It can be applied to predict flow, transport, and reactions in porous media (Lautenschlaeger et al., 2023;Guiltinan et al., 2021;Liu et al., 2021).Especially over the last decade, LBM has gained importance both technically and application-wise.On the one hand, the ease of meshing and parallelization makes LBM extremely favorable for high-performance computing (Latt et al., 2021;Krause et al., 2021;Kellers et al., 2023).On the other hand, it can be applied to a broad variety of application fields.In the context of porous media these are for example: oil and gas flow in underground formations (Ren et al., 2015;H. Li et al., 2015), oil recovery with in-situ combustion (Lei & Luo, 2022), pore structure evolution in cement manufacturing (Patel et al., 2021(Patel et al., , 2014)), combustion of porous solid rocket fuel (Wang & Zhu, 2018), water transport and reactions in fuel cells (Sarkezi-Selsky et al., 2022), multi-phase flow in batteries (Lautenschlaeger, Prifling, et al., 2022;Lautenschlaeger, Weinmiller, et al., 2022;Danner et al., 2016), and dissolution and precipitation reactions in stone formations (Chen et al., 2014;Kang & Lichtner, 2013;Tian & Wang, 2017;L. Zhang et al., 2019).
In most of these cases, dissolution and precipitation play an important role.Using LBM, such reactions are typically modelled as heterogeneous reactive boundary conditions (RBC).Different approaches to translate macroscopic reaction behavior into mesoscopic schemes are known from literature.They range from simple modified bounce-back schemes (Kang et al., 2002(Kang et al., , 2007)), over pseudo-homogeneous reactions (Patel et al., 2014), and schemes based on interpolation using the wall normal (Walsh & Saar, 2010;Chen et al., 2013;Xie et al., 2021;L. Li et al., 2017), up to general local boundary schemes for first-order (Patel, 2016;Verhaeghe et al., 2006;Ju et al., 2020) and higherorder reactions (Kashani et al., 2022;Hiorth et al., 2013).Two important scheme characteristics are generality and locality.Here, general means that the underlying scheme is the same regardless of the surrounding cells, i.e. corners, edges and flat boundaries use the same scheme.Whereas local means that only information of the current boundary cell is used.In combination, they are a necessary basis for versatile and computationally efficient RBCs applicable to complex geometries.
Several approaches have been considered to derive such general local boundary schemes.Three of which are discussed in more detail in this paper: 1) Verhaeghe et al. (2006) was amongst the first to develop such a model for multi-component LBM using the momentum transfer analysis of Bouzidi, Firdaouss, and Lallemand (2001).2) Patel (2016) developed an approach that bases on bouncing back the non-equilibrium part and including the flux as a concentration gradient while correctly capturing the macroscopic wall diffusion.3) Ju et al. (2020) extended an approach, which was originally developed for finite differences (T.Zhang et al., 2012) and then converted to a local scheme (Meng & Guo, 2016), to also consider the correct macroscopic wall diffusion.All aforementioned RBC schemes show robust behavior in cases where the RBC normal is aligned with the lattice.However, cases different from these benchmarks are either missing in the corresponding papers or they show a nonphysical behavior of the corresponding scheme there.
Therefore, in this paper, a novel RBC scheme is presented which is based on the work of Verhaeghe et al. (2006), Ju et al. (2020), andPatel (2016).It is shown to overcome most deficiencies of these methods.It captures the correct wall diffusion and wall normal behavior, and is applicable to a broad range of reaction regimes.It is shown that the new RBC scheme can be combined with the rest fraction method (Sullivan, Johns, et al., 2005) to enable numerically stable simulations even at large ratios between the solute diffusivity and solvent viscosity.
The new RBC scheme as well as the RBC schemes of Verhaeghe et al. (2006), Ju et al. (2020), and Patel (2016) are tested using numerous benchmark cases that vary in complexity.For the sake of comparability, all schemes are reformulated to a common notation.The test cases include: 1) A robust reaction-diffusion verification case for which analytical solutions exists.2) Pattern formation in a precipitation process.3) A sophisticated benchmark case from the literature (Molins et al., 2020), considering reactions in channel flow.All simulations are conducted in 2D.They can, however, be easily extended to 3D.
The paper is organized as follows: In Section 2, LBM basics, the rest fraction method, as well as the coupling of LBM to dissolution and precipitation are introduced.This section also covers the reformulation of the aforementioned RBC schemes.The verification and validation of the different schemes are compared and discussed in Section 3. Finally, the findings are summarized and conclusions are given in Section 4.

Numerical Methods
All LBM simulations presented in this work were conducted using an extended version of the LBM simulation package Palabos (Latt et al., 2021).There, the rest fraction method, reactive boundary conditions, and dissolution and precipitation were implemented.

LBM Fundamentals
For classical fluid flow, LBM solves the Boltzmann equation discretized in the phase-space mostly on a uniform grid.The velocity space is represented by the so-called lattice.It consists of a set of velocities {e i }, weights {w i } and a specific lattice speed of sound c s .It forms the framework to distribute a set of particles that are represented by a particle distribution function, typically called population.In the following, it is indicated by either f i or g i , where i is the lattice direction.Fluid flow emerges through local collision and subsequent streaming of the populations.For reactive flows, the hydrodynamic carrier fluid and the advected scalar field are distinguished using separate lattices and equations.They are coupled by the fluid velocity, determined from the carrier fluid, which enters the advected scalar field.

Hydrodynamic Equations
The dynamics of the carrier fluid described by the Navier-Stokes (NS) equations are solved using the lattice Boltzmann equation combined with the well-known BGK collision operator (Krüger et al., 2017) where τ NS is the relaxation time which is determined by the fluid viscosity ν = c 2 s (τ NS − 0.5) and f eq i is the equilibrium population The fluid density ρ and the fluid velocity u are determined from f i as In the following, the D2Q9 lattice is applied to solve the hydrodynamic equations in 2D.Further details are given in Appendix A.
At high relaxation times, the BGK collision operator introduces numerical slip at bounce-back walls leading to inaccuracies and instabilities.To overcome these artefacts, the two relaxation time (TRT) collision operator (Ginzburg et al., 2008;d'Humières & Ginzburg, 2009) is used in such cases.It is given as (5) Here, ī is the direction opposite to i, i.e. eī = −e i .The values of τ + NS and τ − NS are defined as from the viscosity of the fluid ν and the so-called "magic" parameter Λ.In the following, Λ = 3 16 , if not specified otherwise.It ensures that walls implemented via the bounce-back scheme are located exactly halfway between a solid and a fluid node.

Advection-Diffusion Equations
Advection and diffusion (AD) of a scalar field are described in a similar manner as the hydrodynamic equations given by Eq. (1).In the following, AD populations are indicated by g i .
Here, τ AD is determined by the diffusivity of the scalar D = c 2 s (τ AD − 0.5) and g eq i is Here, u is the local advection velocity of the corresponding fluid field (cf.Eq. ( 3)).From Eq. ( 7), the concentration C of the scalar field is determined as Note, that using this local advection velocity, the first moment of g i , i.e. momentum, is not conserved.Therefore, Eq. ( 7) solves only the AD equation but not the NS equations (cf.Eq. ( 1)).With fewer degrees of freedom, it is sufficient to use a reduced D2Q5 lattice (cf.Appendix A) and a linear equilibrium function (cf.Eq.( 8)) to solve Eq. ( 7) which significantly decreases computational efforts.
The scalar field can also be solved using the TRT collision operator, for both accuracy and stability reasons.The corresponding set of equations reads Analogous to the hydrodynamic equations, τ + AD and τ − AD depend on both the diffusivity and the magic parameter, Λ = 3 16 by default:

Rest Fraction Method
Coupling NS and AD via the advection velocity requires the same spatial and temporal discretization of both fields.Together with the two facts that 1) τ j should be within the range of τ j ≈ [0.6, 1.4] to ensure numerical stability and accuracy (He et al., 1997;Krüger et al., 2017, p.188) and 2) τ j is proportional to ν or D (and therefore also to the Reynolds or Péclet number) this strongly limits the range of applicable Reynolds-to-Péclet ratios.
Using the TRT method where τ j ≈ [0.6, 4.0] can slightly extend this range.However, a technique to significantly increase the range is the rest fraction method.It decouples the discretization and relaxation time by redefining the equilibrium function of the scalar field (Sullivan, Sani, et al., 2005;Sullivan, Johns, et al., 2005).Following (Looije et al., 2018), here, the original formulation is interpreted as a method to change c 2 s as well as the weights w i by the rest fraction J 0 where "dim" is 2 for D2Q5 and 3 for D3Q7.This approach does not affect the other equations, i.e.Eqs.(1)−( 12), unlike the original formulation.Note that setting the free parameter J 0 = 1/3 for D2Q5 or J 0 = 1/4 for D3Q7 results in the commonly used values for c 2 s and w i .The rest fraction method allows for tuning the relation between diffusivity and relaxation time, with D = (1−J0) dim (τ AD − 0.5).Therefore, Reynolds-to-Péclet ratios of up to 1000 can be achieved by adjusting solely J 0 , while keeping τ ≈ 1.
General Local RBC for Dissolution and Precipitation using LBM

Reactive Boundary Conditions
All RBC schemes discussed in this paper macroscopically describe a first-order equilibrium reaction given as Here, J R is the reaction rate, interpretable as a population flux.k r is the reaction rate coefficient and C eq and C wall are the equilibrium concentration and concentration at the wall, respectively.The AD scalar field can also represent temperature, in which case the population flux J R describes a heat flux, instead of a reaction.
The macroscopic flux J R cannot be directly used in LBM due to its mesoscopic character, where the flux is determined using the populations.In addition, C wall is unknown and needs to be estimated using an appropriate mesoscopic scheme.Details of how to transform the macroscopic first-order equilibrium reaction (cf.Eq. ( 15)) into a mesoscopic boundary condition (cf.Eq. ( 16)) are given in the supporting information (SI) Text S1.1.
Fig. 1 schematically shows the functioning principle of the set of RBC schemes discussed in this paper.It starts with the population g * i = g i (x, t 0 ) pointing towards the wall and before streaming.After streaming this population becomes gi = g i (x + e i ∆t, t 0 + ∆t).The actual main step of the RBC scheme determines the unknown population gī from the known gi .For the next streaming step gī is required, and it is the final outcome of an RBC scheme.The RBC schemes of Verhaeghe et al. (2006), Patel (2016) and Ju et al. (2020), and the new RBC scheme are presented in the following.For brevity, in the following only the first authors are mentioned and the corresponding schemes are abbreviated to S.V., S.P.J. and S.New, respectively.

General Local RBC for Dissolution and Precipitation using LBM
For comparability they are reformulated to a similar notation.
It is shown in detail in the SI (cf.Text S1) that they can be expressed in the form The difference of the schemes only enters by the definition of the term and more specifically by the incorporation of the wall normal term (eī • n) and the correction factor γ = τ AD /(τ AD − 0.5).This reformulation approach also allows cross-implementation of already developed features from one RBC to the other, e.g.moving walls.Here, it is used to highlight the differences of the schemes.
In the scheme of Verhaeghe the wall normal term is included in the numerator and γ is not considered.In the schemes of Patel and Ju the wall normal term is included in the denominator and γ is considered.Compared to the scheme of Verhaeghe, in the schemes of Patel and Ju the γ can be interpreted as a shift in the effective reaction rate.The novel scheme proposed here is only slightly different and combines elements from both the scheme of Verhaeghe and the schemes of Patel and Ju.This, however, significantly improves both accuracy and applicability for complex geometries, and overcomes the τ -dependent wall diffusion as is shown in Section 3.
For all three schemes, the limits of the reaction rate are similar (Verhaeghe et al., 2006;Ju et al., 2020;Patel, 2016): At high reactivity (k r → ∞), the schemes simplify to the anti bounce-back scheme (anti-BB) gī = 2w i C eq − gi , which sets the wall concentration to the equilibrium concentration.If no reactions occur (k r = 0), the schemes simplify to the common bounce-back method, i.e. gī = gi .

Dissolution and Precipitation
Finally, simulating phase changes such as dissolution and precipitation at solid surfaces is built up on the RBC scheme.After the determination of the reaction flux using the RBC scheme, the flux enters a method that determines the amount of phase being changed as well as the direction of this process.In this paper, the so-called adjacent growth method (Kang et al., 2010;Kang & Lichtner, 2013) is used.It keeps track of the solid fraction φ for each cell during the simulation, where the change of φ, i.e. ∆φ, is determined as Here, V m is the dimensionless molar volume and a m is the specific surface area.The summation term represents the reaction flux determined from the RBC scheme (cf.Fig. 1).In this study and as proposed by Kang, Lichtner, and Zhang (2006), ∆t = 1 and a m = 1 are used.
During the run time of the simulation φ is dynamically updated.If φ in a cell exceeds a certain threshold, the phase type of this cell is changed from solid to fluid for dissolution, or vice versa for precipitation.In this study, threshold values of Kang et al. (2006) and Pedersen, Jettestuen, Vinningland, and Hiorth (2014) are used.They used φ = 0 for dissolution and φ = 1 for precipitation.This mimics a phase change hysteresis.For precipitation, the local concentration is converted to φ and any excess is distributed to the surrounding fluid cells.The amount being distributed is weighted by the reaction rate of the cell.For dissolution, an initial concentration equal to the average concentration of the surrounding fluid cells is set and an the corresponding amount of φ is removed from the surrounding solid cells.

Results
The RBC schemes of Verhaeghe, Patel and Ju, and the new scheme are compared and verified using a variety of simulation tests.These concern the accuracy of the schemes using a robust 2D reactiondiffusion verification case with analytical solution, pattern formation during precipitation, and an encompassing case study with a reactive circular object inside a channel flow proposed by Molins et al. (2020).
The different simulation tests are conducted for a broad range of advection, diffusion, and reaction regimes.They are characterized by the non-dimensional Reynolds (Re), Péclet (Pe), Damköhler (Da), and Péclet-Damköhler (PeDa) number, which are given as Here, U and L are the characteristic velocity and length, respectively.Note that the Péclet-Damköhler number, which opposes reaction and diffusion rates, defines the behavior of pure reaction-diffusion problems.For PeDa → ∞ diffusion is limiting and for PeDa → 0 it is reaction.

Analytical Reaction-Diffusion Verification
A simple 2D reaction-diffusion problem is studied, for which also an analytical solution exists (cf.Eq. ( 22)).The corresponding simulation setup is shown in Fig. 2. The bulk is solved using the BGK collision operator.The simulation considers pure diffusion of a concentration field within a rectangular domain, thus u = 0.At the left boundary, a constant concentration C 0 is defined.The bottom and right boundaries are set to ∇C = 0.At the top boundary a first-order equilibrium reaction is modelled using either the scheme of Verhaeghe, the schemes of Patel and Ju, or the new scheme.
The analytical solution to this problem (Carslaw & Jaeger, 1986, p. 167) using the formulation of Kang et al. (2006) is where C (analyt) (x, y) is the local analytical concentration, C eq the equilibrium concentration at the top boundary, and a and b are the width and height, respectively.In this paper, the sum in Eq. ( 22) is evaluated up to n=100.The parameters N 2 n is given as and β n are the solutions to the following transcendental Simulations are conducted for PeDa numbers in the range of PeDa = [10 −1 , 10 5 ] and for the two

Configuration A: Aligned System
In configuration A, the term e i • n in Eqs. ( 17)−( 19) simplifies to 1 for e 2 pointing towards the simulation domain, and 0 otherwise.Thus, the implementation of the wall normal does not affect k i (cf.Eqs. ( 17)−( 19)) such that the schemes of Patel and Ju, and the new scheme are identical; only the scheme of Verhaeghe differs by the factor γ.
The simulation results of the concentration field are shown in Fig. 3 for PeDa = 1 and PeDa = 100.The colors show the concentration.The analytical solution is depicted by the black contour lines for the distinct values given in the plot.In addition, the corresponding absolute error is shown in Fig. 4. For the predominantly diffusion-limited case (PeDa = 100), shown in Figs. 3 and 4 d) to f), all RBC schemes are at least in good agreement with the analytical solution.While for the scheme of Verhaeghe the absolute error is < 10 −2 over large areas of the simulation domain, the schemes of Patel and Ju, and the new scheme are even more accurate (absolute error < 10 −3 ).This is, however, different for the predominantly reaction-limited case (PeDa = 1), shown in Figs. 3 and 4   There, the scheme of Verhaeghe shows a significant deviation from the analytical solution, while the schemes of Patel and Ju, and the new scheme are again in almost perfect agreement with the analytical solution.
For configuration A, the difference from the scheme of Verhaeghe compared to the schemes of Patel and Ju, and the new scheme is solely attributed to the absence of γ.As described in the literature (Ju et al., 2020;Patel, 2016), including γ ensures that the macroscopic diffusivity at the wall is captured correctly.Multiplying the PeDa number in Eq. ( 24) by a factor of 1/γ results in no error when using the scheme of Verhaeghe, confirming the shift in effective PeDa number simulated.
Additionally, this setup is used to confirm that integrating RBC with the rest fraction method and TRT collision operator is valid.The resultant MAE increase is less than 1% (cf.SI Text 2), and thus does not result in significant additional errors.

Configuration B: Rotated System
In configuration B, primarily the influence of the wall normal on k i (cf.Eqs. ( 17)−( 19)) is tested.The simulation setup is similar to configuration A, with the only difference that it is rotated clockwise by 45 • .This leads to the term e i • n in Eqs. ( 17)−( 19) being equal to − √ 2/2 for e 1, 2 , i.e. pointing towards the simulation domain.Thus, none of the three schemes compared here are identical anymore.Deviations from the configuration A predominately originate from the implementation of the wall normal.
Again, for PeDa = 1 and PeDa = 100, the simulation results and corresponding absolute error are shown in Figs. 5 and 6, respectively.The meaning of the colors and lines is identical to those from Figs. 3 and 4.
The simulation results of the scheme of Verhaeghe and the new scheme are similar to those from configuration A. Again, multiplying the PeDa number in Eq. ( 24) by a factor of 1/γ, to match the analytical PeDa to the simulated PeDa, results in no error of the the scheme of Verhaeghe.However, a large deviation from the analytical solution is observed for the schemes of Patel and Ju especially at low PeDa values.

Configuration A & B for Extended PeDa Range
The MAE (cf.Eq. ( 25)) of the RBC schemes are summarized in Fig. 7 for both configuration A and B over a vast range of PeDa values.When PeDa → ∞, the RBC schemes approach the anti-BB method.This was shown in the literature (Verhaeghe et al., 2006;Ju et al., 2020) and is further elaborated in the SI Text S1.3.Note, that the anti-BB sets the equilibrium concentration at the boundary.Therefore, this approximation is expected to have significant errors at low PeDa values.Still, it is included as reference in the extended MAE plots.For configuration A (cf. Fig. 7 a)), it clearly shows how the scheme of Verhaeghe has significant errors for 10 −1 PeDa 10 2 , with an decreasing error for when PeDa→ ∞.On the other hand, the schemes of Patel and Ju, and the new scheme are barely affected by the choice of PeDa.This highlights the impact of γ.The MAE decreases for all schemes in the reaction-limited case, i.e PeDa < 1, which has the trivial solution of C(x, y) = C 0 , with the RBC approaching bounce-back behavior.
For configuration B (cf. Fig. 7 b)), the scheme of Verhaeghe and the new scheme behave the same as in configuration A (cf. Fig. 7 a)).However, the scheme of Patel and Ju shows a significant error, due to their wall-normal implementation.In contrast, the new scheme developed in the current study, is almost not affected by the choice of PeDa for both configuration A and B.
The MAE values for RBC and anti-BB schemes converge for PeDa > 10 3 in configuration A and for PeDa > 10 4 in configuration B. Thus, simulations with PeDa > 10 4 show similar results for all aforementioned RBC schemes, i.e. anti-BB, the schemes of Verhaeghe, Patel and Ju, and the new scheme.

Pattern Formation During Precipitation
This test case is motivated by the work of Pedersen et al. (2014).They investigated pattern formation of grain growth during precipitation and its dependence on the grid orientation.For their study they used the scheme of Verhaeghe.
The simulation domain is circular with a diameter of 300 grid cells.This is realized by adding boundary cells along this circle in a square domain of size 300 lm × 300 lm cells.On the boundary the concentration is constant, i.e.C circ = 1.An initial seed is placed in the center of the domain.Its General Local RBC for Dissolution and Precipitation using LBM shape is defined by the parametric equation x(s) = [0.1 + 0.02 cos(8πs)] cos(2πs), y(s) = [0.1 + 0.02 cos(8πs)] sin(2πs), ( 26) with s as the parameter.Two variants of this setup are simulated: One where the initial seed is aligned with the grid orientation, and another where the initial seed is rotated clockwise by 19 • .The boundaries of the seed are defined by a RBC with C eq = 0 and k r = 0.0014 (cf.Eq. ( 16)).The bulk is solved using the BGK collision operator.Pure diffusion is considered such that u = 0 within the domain.To compute the wall normal, the simple and efficient isotropic finite difference method (Kumar, 2004) is used.
The results determined with the scheme of Verhaeghe, the schemes of Patel and Ju, and the new scheme are shown in Fig. 8.The resulting patterns when using the scheme of Verhaeghe and the new scheme are similar, as shown in Fig. 8 a) & d) and Fig. 8 c) & f), respectively.However, for the scheme of Verhaeghe the growth is less pronounced which is due to the effective PeDa shift.A correction of this shift further improves the accordance of both schemes.In contrast, the pattern resulting from the schemes of Patel and Ju shows dendritic behavior and grid dependence, shown in Fig. 8 b) & e).The reason is due to the implementation of the wall normal in k ( S.P.J.) i (cf.Eq. ( 18)), leading to faster reactions in diagonal walls compared to straight walls.

Reaction in a Channel Flow
This case study considers a reactive object in a channel flow, i.e. under the influence of advection.It is motivated by the extensive comparison paper by Molins et al. (2020) in which they compared five General Local RBC for Dissolution and Precipitation using LBM different simulation methods (Chombo-Crunch, dissolFoam, LBM, OpenFoam, and Vortex).Both Molins et al. (2020) and Ju et al. (2020) used the same dissolution experiment of Soulaine, Roman, Kovscek, and Tchelepi (2017) to validate their methods.
The simulation setup is shown in Fig. 9.It is 2D and consists of a channel with height H = 0.05 cm and width W = 0.1 cm, simulated with resolution 600 lm × 1200 lm and a lattice velocity of U LB = 0.001 lm ls −1 .The channel is initially filled with a fluid at rest and with constant concentration equal to that of the inlet.A circular reactive object with diameter D = 0.02 cm is placed in the middle of the channel.The boundary conditions for the carrier fluid are 1) bounce-back walls at the top, bottom, and at the object's surface, 2) plug velocity profile u in = 0.12 cm/s at the inlet on the left, and 3) fixed density ρ out = 1 at the outlet on the right.The boundary conditions for the concentration field are 1) bounce-back walls at the top and bottom, 2) constant concentration at the inlet C in = 10 mol/m 3 , 3) zero concentration gradient at the outlet, and 4) a RBC at the object's surface.The RBC is varied between the three schemes and the anti-BB.Large ratios of Pe to Re are realized using the TRT collision operator and rest fraction method for the concentration lattice.The following values are chosen: J 0 = 0.99 for Pe=600 and J 0 = 1/3 for Pe=6.First, the flow of the carrier fluid, i.e. the NS lattice, was simulated until steady state was achieved.Then, the concentration field, i.e. the AD lattice, was solved and the average reaction rate of the object was determined as Here the integral is performed on the outlet.Simulations were conducted for Re = 0.6 and different combinations of Pe = [6, 600] and Da = [0.178,17800].These are: case 1) Pe = 600, Da = 178 (diffusion limited); case 2) Pe = 600, Da = 17800 (diffusion limited); case 3) Pe = 6, Da = 178 (advection-diffusion limited); and case 4) Pe = 6, Da = 0.178 (advection-reaction-diffusion limited).As is shown in the following, the cases 1) to 3) are not suited to differentiate the schemes, since they are not reaction limited.
A qualitatively comparison of the results of case 4) for the schemes is given in Fig. 10.A compar- ison between the cases 1) to 4) for the new scheme is shown in Fig. 1 in the SI.In Case 4) when using the scheme of Verhaeghe (cf.Fig. 10 a)), the simulation is strongly reaction limited, as indicated by the high concentration around the object's surface.In contrast, the anti-BB scheme (cf.Fig. 10 d)) shows the case where reaction is not limiting at all.Both the new scheme and the schemes of Patel and Ju (cf.Table 2 shows that for the cases 1) to 3), the simulations are not reaction limited, resulting in the RBCs behaving similar to the anti-BB.Thus, cases in which the PeDa ≥ 1000 are not suited for validation of RBC schemes.The experiment of Soulaine et al. (2017) has a PeDa = 3000.In case 4), the reaction is not limiting anymore, as indicated by the much higher anti-BB result.Here, the result of the various RBC schemes diverge.Compared to the new scheme, the scheme of Verhaeghe underestimates and the schemes of Patel and Ju overestimates R avg by the significant factors of 0.43 and 1.6 respectively.

Conclusion
A new general local reactive boundary condition (RBC) scheme was presented that accurately captures first-order equilibrium reactions for a wide range of Péclet (Pe) and Damköhler (Da) numbers for complex geometries.It combines aspects of previous RBCs from Verhaeghe et al. (2006), Ju et al. (2020), andPatel (2016), but overcomes their deficiencies.This means, the new RBC scheme does not suffer from a τ -dependent wall diffusion and accurately considers the wall normal for complex geometries and surfaces that are not aligned with the simulation grid.
In this study, all aforementioned RBC schemes were first reformulated to a similar notation.Then, all RBC schemes were tested and compared using three different verification cases.These were chosen such that they were relevant, representative, and covered important features of RBCs and their applications: 1) A robust 2D reaction setup purely driven by diffusion.The problem has an analytical solution and the RBC schemes were studied over a wide range of PeDa values as well as under grid rotation.2) Pattern formation in a precipitation process.Here, especially the impact of rotational variance and the behavior of the wall normal was studied.3) An advanced reaction setup including advection around a reactive object in a channel flow.The RBC schemes were coupled to the rest fraction method, a wide range of PeDa values was studied, and the results were compared to results from the literature.
The test cases demonstrate a broad applicability and high accuracy of the new RBC scheme.For the 2D reaction setup, the new RBC scheme shows the best accordance with the analytical solutions.Additionally, it was shown that it correctly simulates wall diffusion and is invariant with respect to relaxation time and grid orientation.For the simple precipitation case, the new RBC scheme showed physically sound pattern formations.Moreover, in contrast to the other RBC schemes, it was not affected by grid rotation.For the advection case, large ratios of Péclet and Reynolds number (≤ 1000) were successfully simulated.Reaction rates determined using the new RBC scheme are in general accordance with results from the literature (Molins et al., 2020).
All in all, the new RBC accurately describes first-order reactions and can be applied to simulate precipitation and dissolution phenomena even for complex geometries.This is a strong advantage over other RBCs previously described in the literature.In addition, as a general and local scheme, it is easy to implement for both 2D and 3D simulations, it is computationally efficient and facilitates parallel computation.The new RBC scheme is missing some advanced features, e.g.moving walls, or interpolated sub-grid wall locations.However, the presented general reformulation allows those already developed features to be incorporated easily.
The new RBC can be used to study reaction processes in a broad range of research fields.Potential applications might be reactive flows through porous media, with and without dissolution and precipitation (Pereira, 2022;L. Zhang et al., 2019;D. Zhang et al., 2021;C. Zhang et al., 2021;Xu et al., 2018;Jiang et al., 2021), pore structure evolution in cement manufacturing (Patel et al., 2021), or morphological changes due to the conversion of active material in energy storage systems (Fang et al., 2021).Especially the last topic will be in the focus of our future work.

Open Research
An extended version of the Parallel Lattice Boltzmann Solver (version 2.3.0)-short Palaboswas used for all LBM simulations in this study.The original version of Palabos is preserved at General Local RBC for Dissolution and Precipitation using LBM https://palabos.unige.ch/,available via GNU Affero General Public License version 3 without login required and developed openly at https://gitlab.com/unigespc/palabos(Latt et al., 2021).The reactive boundary conditions and the rest fraction method were implemented separately.

Introduction
This supporting information provides details on the reactive boundary conditions schemes, their reformulation and behavior in the limit.Additionally, the supporting information contains a section on the aggregation of Molins et al. (2020) simulation data.

Text S1 Reactive Boundary Condition Expressions
In the following, the derivation of the reactive boundary condition (RBC) of Verhaeghe, Arnout, Blanpain, and Wollants (2006) is given in detail.In addition, the approach to reformulate the RBC schemes of Ju, Zhang, and Guo (2020) and Patel (2016) to the notation of Verhaeghe et al. is shown.As in the main text, only the first authors are mentioned in future discussions.The nomenclature and symbols are chosen in a consistent manner with respect to the current paper.
Text S1.1 Derivation of the RBC scheme of Verhaeghe The scheme of Verhaeghe is based on a Dirichlet boundary condition implemented using the anti bounce-back method (cf.Eq. (S.1)), a wall flux analogous to moving walls (cf.Eq. (S.2)) (Bouzidi et al., 2001), and a definition of the reactive flux (cf.Eq. (S.3)).The corresponding equations read: S.1 with the outgoing population g i , the incoming population gī and the wall concentration C w , and where momenta flux due to the wall movement is interpreted as concentration flux along the wall normal where J R is the reactive flux defined by a first-order equilibrium reaction.From the aforementioned equations, the scheme of Verhaeghe is derived as given in the following.Solving Eq. (S.1) for C w and inserting it into Eq.(S.3) gives Inserting Eq. (S.4) into Eq.(S.2) results in Finally, applying c i = −cī in the definition of k i results in kī = −k i , and finally leads to This final equation (cf.Eq. (S.8)) is similar to the result from Pedersen, Jettestuen, Vinningland, and Hiorth (2014).

Text S1.2 Reformulation of the RBC schemes of Ju and Patel
Starting from the original formulations of the schemes of Ju and Patel it is shown under which conditions they can be reformulated similar to the scheme of Verhaeghe (cf.Eq. (S.8)).
The original formulation of Patel (2016, p.69) is Here, V and A s are the volume and surface area, respectively, of the solid at the boundary.They are unity when non-dimensionalized, resulting in .

S.10
The original formulation of

S.12
The reformulation of Eqs.(S.10) & (S.11) are performed with the relation 2w i = c 2 s .It can be applied only for reduced lattices along all directions, except of the rest velocity i = 0.The latter, however, relates to g 0 which is not modified by any of the mentioned boundary schemes.
The reformulation steps required for the scheme of Patel are as follows: Using the second definition of Eq. (S.12), 1 γ in Eq. (S.10) gives gī = k r C eq − gi S.17 Here, the w i C w − g i term is the discrepancy between the incoming and outgoing population, which is indicative of the gradient at the boundary itself.
The RBC schemes to approach the anti-BB scheme, i.e. ∆g (ABB vs RBC) ī → 0, under the following conditions.Either for 2 / (1 + kī) → 0 which holds for increasing values of kī and Da, or for (w i C w − gi ) → 0 which holds for a decreasing gradient of the wall concentration and increasing value of Pe.Finally, both conditions together show that the RBC schemes approach the anti-BB scheme as PeDa → ∞.

Text S2 Detailed error analysis of the rest fraction method
The additional error introduced by applying the rest fraction method is quantified using the MAE.Simulations are conducted using the new scheme for the simulation setup of configuration A of the 2D reaction-diffusion problem described in Section 3.1.It is tested against the simulation without rest fraction as baseline.Both BGK and TRT collision operators are used with the following parameter variations: J 0 = {0.1,0. 3, 0.5, 0.9} at PeDa=1, where k r is kept constant resulting in τ AD = {0.6 7, 0.74, 0.82, 2.1}.The results are shown in Table S1.As expected results for the J 0 = 0. 3 and the baseline are identical.For the BGK simulations, the J 0 = 0.9 simulation was unstable due to a too high τ AD , which can be mitigated with the TRT collision operator.For both BGK and TRT collision operators, the difference in the MAE between the baseline and the rest fraction simulations are negligibly higher.Overall, the rest fraction method and TRT collision operator can be used to simulate large Pe to Re ratios.

Figure 1 :
Figure 1: Schematic representation of the first-order RBC.Fluid nodes (circles) are depicted on the left side.Solid nodes (squares) are depicted on the right side.The top and the bottom rows only differ in time, but not in space.The streaming of population g * i = g i (x, t 0 ) to gi = g i (x + e i ∆t, t 0 + ∆t) are indicated by the orange arrows; the result of the RBC gī = gī(x + e i ∆t, t 0 + ∆t) with the blue arrow.In addition, the top right quadrant shows the D2Q5 velocity set (e i for i ∈ {0, 1, 2, 3, 4}).

Figure 2 :C
Figure 2: Simulation setup of the 2D reaction-diffusion problem.The boundary conditions are constant concentration (left), zero concentration gradient (right and bottom), and first-order equilibrium reaction (top).An exemplary concentration field and its analytical solution are indicated by the background and the contour lines, respectively.Two variants of this setup are studied; one where the boundaries are aligned with the D2Q5 grid (configuration A) and another which is rotated clockwise by a 45 • angle (configuration B).

Figure 3 :
Figure 3: Simulation results of configuration A. The concentration field is shown by the color code given in the legend.The analytical solution (cf.Eq. (22)) is indicated by the black contour lines for the distinct values given in the plot.
a) to c).General Local RBC for Dissolution and Precipitation using LBM

Figure 4 :
Figure 4: Absolute error of configuration A. The error is shown by the color code given in the legend.

Figure 5 :
Figure 5: Simulation results of configuration B. The concentration field is shown by the color code given in the legend.The analytical solution (cf.Eq. (22)) is indicated by the black contour lines for the distinct values given in the plot.

Figure 6 :
Figure 6: Absolute error of configuration A. The error is shown by the color code given in the legend.

Figure 7 :
Figure 7: MAE analysis of configuration A and B with extended PeDa range.The MAE of the schemes of Patel and Ju, and the new scheme are identical for configuration A.

Figure 8 :
Figure 8: Pattern formation during precipitation.The final grains are shown in black.The initial grain seeds and the simulation domain are indicated by the white line and the gray-shaded background, respectively.

Figure 9 :
Figure 9: Simulation setup of the 2D reactive object in channel flow from Molins et al. (2020).The boundary conditions for 1) the fluid lattice are constant plug flow (left), constant density (right) and bounce-back walls (top, bottom, object), and for 2) the scalar lattice are constant concentration (left), zero concentration gradient (right), bounce-back walls (top, bottom) and RBC (object).

Figure 10 :
Figure 10: Concentration field in the vicinity of a reactive object (grey circle) in a channel flow simulated for case 4).Using a) the scheme of Verhaeghe, b) the schemes of Patel and Ju, c) the new scheme, and d) the anti-BB scheme, which sets the concentration to the equilibrium concentration.The concentration is shown by the color code given in the legend.
Fig. 10 b) & c), respectively) fall between the two extremes.The numerical results of R avg are given in Table 2 together with the values from Molins et al. (2020).There, ranges are given as different methods were used to determine R avg (cf.SI Text S3).Overall, the RBC schemes agree well with the results of Molins et al., where available, indicating good integration of the RBC schemes with the rest fraction method.

Figure S1 :
Figure S1: Concentration field in the vicinity of a reactive object (grey circle) in a channel flow simulated using the new scheme (cf.3.3 Reaction in a Channel Flow in the main paper).a) & b) Cases 1) & 2) strongly diffusion limited and advection-reaction driven, causing the tail with low concentration behind the object.c) Case 3) advection-diffusion limited and reaction driven, indicated by a broader area of low concentration.d) Case 4) advection-reaction-diffusion limited, with high concentration at the object's surface.The concentration is shown by the color code given in the legend.

Table 1 :
Parameter set for the 2D reaction-diffusion case

Table S1 :
Impact of the Rest Fraction Method on Configuration A Simulations