## 2 Pore-Scale Reactive Transport
### Species Transport  
When a problem needs to be solved for a range of parameters (as is our goal here), working with the dimensionless form of the equations give definitive advantages. The height of the computational domain \(\Omega_{f}\), namely \(l\), is used for scaling spatial sizes. The scaling of the velocity is done by the inlet velocity \(\bar{u}\), and the scaling of the concentration is done by the inlet concentration \(\bar{c}\). The concentration of the solute in the fluid is denoted by \(c(\mathbf{x},t)\). The unsteady solute transport under the conditions of heterogeneous reactions is governed by the convection-diffusion equation in its dimensionless form:  
\[\frac{\partial c}{\partial t}+\
abla(\mathbf{u}c)-\frac{1}{\text{Pe}}\
abla^{2}c= 0,\quad\mathbf{x}\in\Omega_{f},\quad t>0, \tag{7}\]  
Figure 1: Sketch of the pore scale domain. where  
\[\text{Pe}=\frac{l\bar{a}}{D}\]  
is the Peclet number. The boundary conditions are written as follows:  
\[c=\bar{c},\quad\mathbf{x}\in\Gamma_{in}, \tag{8}\] \[\
abla c\cdot\mathbf{n}=0,\quad\mathbf{x}\in\Gamma_{sim}\cup\Gamma_{out},\] (9) \[\frac{\partial m}{\partial t}=-\
abla c\cdot\mathbf{n},\quad\mathbf{x} \in\Gamma_{s}. \tag{10}\]  
The concentration of the solute at the inlet is assumed to be known. Zero diffusive flux of the solute at the outlet and on the external boundaries of the domain is prescribed. Note that convective flux via the outlet is implicitly allowed by the above equations. The surface reactions that occur at the obstacles' surface, \(\Gamma_{s}\), satisfy the mass conservation law, in this particular case meaning that the change in the adsorbed surface concentration is equal to the flux from the fluid to the surface (10), where \(m\) is the surface concentration of the adsorbed solute [36]. A mixed kinetic-diffusion adsorption description is used:  
\[\frac{\partial m}{\partial t}=f(c,m). \tag{11}\]  
For reactive boundaries, the choice of \(f\) and its dependence on \(c\) and \(m\) are critical for a correct description of the reaction dynamics at the solid-fluid interface. A number of different isotherms (i.e., different functions \(f(c,m)\)) exist for describing these dynamics, which are dependent on the solute attributes, the order of the reaction, and the interface type. The simplest of these is the Henry isotherm, which assumes a linear relationship between the near-surface concentration and the surface concentration of the adsorbed particles. In its dimensionless form, this is written as follows:  
\[\frac{\partial m}{\partial t}=\text{Da}_{a}c-\text{Da}_{d}m,\quad\mathbf{x}\in \Gamma_{s}, \tag{12}\]  
where the adsorption and desorption Damkoler numbers are given by  
\[\text{Da}_{a}=\frac{k_{a}}{\bar{u}},\quad\text{Da}_{d}=\frac{k_{d}l}{\bar{u}}.\]  
Here, \(k_{a}\geq 0\) is the rate of adsorption, measured in unit length per unit time, and \(k_{d}\geq 0\) is the rate of desorption, measured per unit time. We used this isotherm in our previous work in our diffusion-dominated case investigation [37] and we will come back to it later when comparing the operation of the MBC algorithm in the penultimate section. In this work, we use the Langmuir adsorption isotherm, which is a more complicated three-parameter model:  
\[\frac{\partial m}{\partial t}=\text{Da}_{a}c\Big{(}1-\frac{m_{\infty}}{\text {M}}\Big{)}-\text{Da}_{d}m,\quad\mathbf{x}\in\Gamma_{s}, \tag{13}\]  
where the dimensionless parameter \(M\) is given by  
Figure 1:0  
Here, \(m_{\infty}>0\) is the maximal possible adsorbed surface concentration. In comparison to the Henry isotherm (12), the Langmuir isotherm (13) predicts a decrease in the rate of adsorption, as the adsorbed concentration increases due to the reduction in the available adsorption surface. The formulation of the initial boundary value problem in addition to the governing Equations (1), (2) and (7), boundary conditions (3)-(6), (8)-(10), and specified isotherm (12) or (13) requires specification of the initial conditions:  
Figure 1:1  
Figure 1:2