## 3 Numerical Simulation
### Geometry and Sensitivity Analysis  
The computational domain is a rectangle with a dimensionless height of \(x_{2}=1\) and dimensionless length of \(x_{1}=17.5\), in which 10 half cylinders are embedded. The distance between the centers of the cylinders in the \(x_{1}\) direction is 1.5 dimensionless units, and the radius of the cylinders is 0.4 dimensionless units.  
The computational domain \(\Omega_{f}\) is triangulated using the grid generator Gmsh [42]. The script for preparing the geometry is written in Python. A demonstration of the computational grid is shown in Figure 2.  
For the FEM approximation of the velocity, the pressure, and the respective test functions, the following finite dimensional subspaces are selected: \(\mathbf{V}_{h}\subset\mathbf{V}\), \(\mathbf{\hat{V}}_{h}\subset\mathbf{\hat{V}}\) and \(Q_{h}\subset Q\). Taylor-Hood \(P_{2}-P_{1}\) elements [43] are used here. These are continuous \(P_{2}\) Lagrange elements for the velocity components and continuous \(P_{1}\) Lagrange elements for the pressure field. The computations are carried out using the computing platform for partial differential equations FEniCS (website fenicsproject.org) [44; 45].  
Figure 2: Computational grid.  
As mentioned above, slow flows are mainly of interest for the current study. The convergence of the solution with respect to the refinement of the grid is illustrated in Figure 2(a). The convergence of the solution with respect to the time step is illustrated in Figure 2(b). We have used three computational grids: a basic grid with 18,743 nodes and 35,958 triangles, a coarse grid with 4760 nodes and 8754 triangles and a fine grid with 72,745 nodes and 142,460 triangles. From the results, it can be concluded that the coarse grid provides good enough accuracy for the numerical solution. Taking into account that we will solve a multitude of direct problems, it was decided to conduct research on a coarse grid.  
A sensitivity study is carried out in order to identify how the change of different parameters leads to changes in the breakthrough curves:  
* \(\mathrm{Pe}=[1,10,100]\);
* \(\mathrm{Da}_{a}=[50,100,200]\);
* \(\mathrm{Da}_{d}=[0.5,1,2]\);
* \(\mathrm{M}=[100,1000,100,000]\).  
Such sensitivity studies are often tge first stage for optimization or parameter identification procedures. The dependence of the average outlet concentration from the Peclet number is shown in Figure 2(c).  
The rate of adsorption is characterized by \(\mathrm{Da}_{a}\). The influence of this parameter on the outlet concentration is illustrated in Figure 2(e). Increasing \(\mathrm{Da}_{a}\), as expected, leads to more intensive adsorption and a larger amount of the deposited substance. The influence of the parameter \(\mathrm{Da}_{d}\) on the outflow concentration is illustrated in Figure 2(f).  
The Henry isotherm usually describes the initial stages of the adsorption well. In cases when only a limited mass can be adsorbed at a surface, the Langmuir isotherm should be used to reflect the decay of the adsorption rate close to the saturation. In this case, an additional parameter appears, namely \(\mathrm{M}\). The influence of \(\mathrm{M}\) on the average output concentration is shown in Figure 2(d).  
The results are illustrated in Figures 2(c)-f. The curve in these figures is a breakthrough curve, which describes how much solute has crossed the outflow boundary, and it is given by  
\[c_{out}(t)=\frac{\int_{\Gamma_{out}}c(\mathbf{x},t)d\mathbf{x}}{\int_{\Gamma_{out}}d \mathbf{x}}. \tag{20}\]  
Figure 3: _Cont._