Multiscale simulation of colloids ingressing porous layers with evolving internal structure

We report on a reaction-diffusion model posed on multiple spatial scales that accounts for diffusion, aggregation, fragmentation, and deposition of populations of colloidal particles. The model is able to account for the heterogeneity of the internal porous structure of the layer. For simplicity, we represent the microstructures as discs with prescribed initial random distribution of radii. As microstructures grow due to the deposition of populations of colloidal particles, local clogging becomes possible, that is neighbouring disks may touch each other. We investigate how distributions of evolving microstructures influence the transport and storage properties of porous layers. As working tool, we propose a FD-FEM discretization of the multiscale model. We illustrate numerically local clogging effects on the dispersion tensor and quantify herewith the layer’s performance with respect to both the efficiency of the transport and the storage capacity. The presented model and numerical approach can be extended in a rather straightforward way to handle slightly more complex geometrical settings like thin porous structures with multi-layers in 2D, or single layers in 3D.


Introduction
We are interested in developing mathematical models and simulation tools that can mimic numerically the following Gedankenexperiment: we imagine a porous media having a rather regular internal structure.In 2D, we think of a locally periodic distribution of solid discs immersed into an easy-to-penetrate, compressible background; see e.g.Fig. 1.In 3D, the type of geometries that fit to our story are either collections of balls or of unidirectional parallel fibers [see e.g.Fig. 2 in Printsypar et al. (2019)], both embedded in a compressible environment.
We are now exposing this porous medium geometry to a large population of colloidal particles (small monomers) that like to ingress through the internal structure at hand (e.g.trough the white region in Fig. 1).Besides their spatial dynamics (driven by flow or diffusion), the colloidal particles tend to aggregate to reach bigger sizes (dimers, trimers, etc.) which then finally (i.e. when they are sufficiently big) deposit on the boundary of the microstructures they meet during their walking.This process of colloidal particle aggregation is called flocculation.As a result of this deposition process, the microstructures (black discs in Fig. 1) grow inside the porous medium.Conversely, these microstructures can, in principle, release selected population sizes thereby eventually also shrinking.Such a scenario extends naturally when allowing for chemically-reactive microstructures or when involving first-order phase transitions.The situation becomes numerically much more challenging and mathematically more interesting if the colloids are hot, or if the geometry of the regular medium has stochastic elements.Such extensions are outside of the focus of this work, but it is worth investigating them at some point.
As this thought experiment is rather generic, many porous media applications can be connected to it, at least remotely.Hence, there is a lot of literature available, however results seem to be bound to quite specific scenarios, and, as a consequence, aiming for generality is daunting.We refer here to some of our own previous work (Eden 2019;Krehel et al. 2015;Nikolopoulos 2018;van Noorden and Muntean 2011), as well as to that of others to see how colloidal populations are expected to affect the macroscopic transport in various classes of porous media like biological tissues [for drug delivery tasks (Ray et al. 2013)], natural oil reservoirs [to be unclogged via high pressure water containing inertial particles (Kokubun et al. 2019)], soils [encapsulating organic nutrients (Johnson and Elimelech 1995;Ray et al. 2012Ray et al. , 2018) ) or contaminants (Boccardo et al. 2018)], active membranes [expected to perform filtration (Iliev et al. 2014;Davis 2021;Printsypar et al. 2019;Iliev et al. 2020;Kiradjiev et al. 2021)], old textiles [exposed to the ingress of tiny particles inside museums or outdoors (Bonetti et al. 2019)].The excellent review (Ryan and Elimelech 1996) contains further relevant literature hints for a similar context.It turns out that the involved models are mostly phenomenological, where constitutive laws are assumed to hold and then efforts are paid in estimating to which extent such models can come close to reality.Other modeling ideas attempt to use some level of details known about the media's microstructure (shape, orientation, chemical composition, etc.) and scale up this information to an observable, macroscopic level.The upscaling procedures are typically done at the level of balance laws and deliver, usually under strong technical assumptions on both the averaging techniques and choice of shapes and arrangements of microstructures, Fig. 1 Poisson-disc sampling within a rectangle.Example of regular 2D porous medium that can be approximated via a locally periodic representation the structure of the expected effective balance equations as well as explicit formulas to compute the effective transport coefficients.This second investigation route is in line with our approach in Muntean and Nikolopoulos (2020), where we allowed for the possibility of local clogging, which, particularly in 2D or for thin 3D layers, has dramatic effects.1One of the resulting upscaled models derived in Muntean and Nikolopoulos (2020), was recently studied from the mathematical analysis perspective in Eden et al. (2021).
As main contribution in the context of this paper, we explore further the model from Eden et al. (2021) by means of numerical simulations; the details on the model equations and assumed geometry and parameters are described in Sect. 2. The main feature of our model is multifold: the coupling in the model equations goes over two spatial scales (macro and micro); the characteristic macro scale is fixed, while the characteristic micro scale length changes in time and depends implicitly on the solution itself.Mathematically, the object we are handling is an parabolic-elliptic-ODE system posed on multiple space scales, out of which one involves an a priori unknown moving boundary.This structure fits to what Showalter refers to as partial differential equations with distributed evolving-in-time microstructures; see Showalter (1992) for more context on this type of models.Relying on multiscale simulations, we illustrate in a couple of numerical experiments to which extent local clogging is able to affect macroscopic transport and storage for a particular class of porous media.We are able to trace changes in the effective (tensorial) transport coefficients induced by the evolution of the growing or shrinking microstructures, either receiving or releasing populations of colloidal particles.
The paper is organized as follows: In Sect.2, we describe the geometry of the heterogeneous medium we have in view as well as the set of model equations.Section 3 gives a glimpse on what is now regarding the mathematical structure of the evolution equations involved (from where they come for, well-posedness, approximation aspects), while in Sect. 4 we suggest suitable macroscocpic and microscopic space discretizations.We present in Sect. 5 a couple of numerical illustrations that show the ability of our approach to cope with variations in the of the porous medium's micro-scopic heterogeneities.Finally, we list in Sect.6 our conclusions and a few ideas for further related studies.

Setting of the model equations
Here, we give a short summary of and motivation for our full model.This model was originally derived in Muntean and Nikolopoulos (2020) and further analyzed in Eden et al. (2021), we refer to those works for additional details regarding the problem and its properties.Please note that the focus of this present work, the clogging scenario, was not considered in Eden et al. (2021) and only tangentially touched upon in Muntean and Nikolopoulos (2020).
For a given T > 0, we consider S = (0, T ) to be the time interval during which we observe the overall physical process.The macroscopic region Ω ⊂ R2 is a bounded domain whose Lipschitz boundary admits the outer unit normal vector n = n(x) for almost all x ∈ ∂Ω.It represents the homogeneous description of the target porous material.In addition, let N be a given natural number indicating the maximal allowed size of an aggregate of colloid particles.The word size refers here to the number of primary particles making up the aggregate.For each i = {1, ..., N }, let u i : S × Ω → [0, ∞) (we set u = (u 1 , ..., u N )) denote the molar concentration of aggregates of size i that can be found at point x ∈ Ω and time t ∈ S.
We take the function v : S × Ω → [0, ∞) to represent the mass density of the adsorbed material i.e. the material deposited in the solid matrix (mass that is present in the system but currently does not take part in the diffusion and agglomeration processes); this mass can dissolve again allowing colloidal populations to re-enter the pore space.The processes of adsorption and desorption 2 is modeled in this context via an Robin-type exchange term in the form of Henry's law [see e.g.Helmig (1997); Krehel et al. (2015)] with positive exchange coefficients a i and b i : (1) Here the radius function r : S ×Ω → (0, r max ) (for some r max > 0) acts as a measure of the clogginess of the porous media [see ( 2)].
The exchange coefficient 2πr 1−πr 2 is given as the ratio of the perimeter of the ball (2πr ) and the pore volume; this relationship is justified via the upscaling in Muntean and Nikolopoulos (2020).Note that equilibrium between u i and v is accomplished when u i = a i b i v. To describe the aggregation processes taking place inside the pore space of the medium, we use a truncated variant of the Smoluchowski formulation [we point to Aldous (1999) for a review] given here by with positive coefficients γ jl .This is in fact a truncated version of the discrete Smoluchowski coagulation term [see, e.g., Aldous (1999)] where we only consider colloids up to a maximal size of N > 0 (i.e., γ jl = 0 for j + l > N ).The index i refers to the ith population (i ∈ {1, . . ., N }).The first sum in (2) accounts for the formation via coagulation under the assumption that colloids of size i can be formed when two smaller colloids of sizes j and l with the property j +l = i meet.In particular, the first sum is zero in R 1 (particles of size 1 can not form via coagulation) and 1 2 γ 11 u 2 1 in R 2 (particles of size 2 can only form via coagulation of two size-1 colloids).The second sum accounts for the loss of i-sized colloids by coagulating with a different colloid of size j to form a new one of size i + j ≤ N .The second sum is therefore empty in R N .Note that this is conceptually very similar to viewing the overall production process as a chain of second order chemical reactions.The truncation to a fixed population size holds particularly in the case of colloids in soils.If no truncation is taken into account, then fragmentation terms are needed to prevent the size of the population to grow indefinitely.In that case, the functional setting changes as for instance in Canizo et al. (2010) and references cited therein.Additionally, note that there are other ways to model the coagulation and fragmentation processes.We stick here to (2), an alternative way would be to consider instead, for instance, a discrete Becker-Döring dynamics.
It is important to observe that in the context of porous media the colloidal populations involve a finite size chain of the cluster, i.e. there will be a population of N -mers where N takes the maximum cluster size.For that reason, we deal with a truncated finite sum here.Interestingly, for many applications a good choice of such N is rather low; see e.g.Krehel et al. (2015).
The diffusion-reaction system for the different N populations of aggregates is then given via We refer to this model component as (P i ).Note that the effective porosity is wrapped into the diffusion tensor D i (r ) as this is apparent from Eq. ( 10) as well as in the second term in the right-hand side of Eq. ( 4), the latter being derived by the homogenization process [see the statements in Eden et al. (2021)].
1 Page 6 of 19 GEM -International Journal on Geomathematics (2023) 14 :1 The effective diffusion matrix D i (r ) ∈ R 2×2 can be calculated using any solution Here Y = (0, 1) 2 denotes the unit cell, B(r ) is the closed ball with radius r and center point y = ( 1 2 , 1 2 ), and e k the k-th unit normal vector.Finally, ν = ν(y) denotes the outer unit normal vector of Y\B(r ).Note that w k is only unique up to a constant: this is fine as we are only interested in ∇w k later.
We then have where φ(r ) := 1−πr 2 |Y | denotes the porosity of the medium and d i > 0 are known constants.In our case, since |Y | = 1, we get φ(r ) = 1 − πr 2 .The structure of both the cell problem and of the effective diffusivity are derived in Muntean and Nikolopoulos (2020) using asymptotic homogenization techniques adapted to the case of locallyperiodic media.
Finally, the evolution of v is governed by an ODE parametrized in x ∈ Ω, say (P iii ), and is given via and the radius function is governed by the following ODE parametrized in x ∈ Ω referred as (P iv ) where α > 0 is an additional proportionality coefficient.We have labeled these subproblems (P i ) − (P iv ) for easier referencing.
It is important to keep in mind that changes in time of v and r are closely connected via 2πα∂ t v = ∂ t r or, more explicitly, via For that reason, Problem (P iv ), that is ( 13)-( 14), can be substituted with (15).A possible initial choice for the initial distribution of radii r 0 is depicted in Fig. 2. We point out there also what will happen at the final time T ; more details on the parameter setup are given in the simulation sections.What concerns the modeling of the deposition of the colloidal populations, our choice is similar to one reported in Johnson and Elimelech (1995).
Such choice accounts for the simple observation that the adsorbed material leads to the clogging of the pore under the fundamental assumption of the growth of the radius is proportional to the amount of material that is adsorbed.For a more concrete argumentation for this particular structure, we again point to Muntean and Nikolopoulos (2020).
We can now exploit the structure of the effective diffusion tensor D = [(D i ) jk ] given cf.(10) to point out the corresponding structure of the tortuosity tensor, denoted here by T = [(T i ) jk ].Classical descriptions of the structure of porous media [cf.e.g.Bear (1988)] allow for the relationship where is in fact here φ(r ) the volumetric porosity, while T is the tortuosity tensor.
Together with ( 16), (10) ensure the following structure to the tortuosity tensor: 1 Page 8 of 19 GEM -International Journal on Geomathematics (2023) 14 :1 The relationship ( 16) is supposed to hold for general 3D porous domains, just that in general T is hard to estimate.For thin porous domains, alternatives have been found cf.e.g.Patasius et al. (2005), Grisan et al. (2003) and Idris et al. (2022) where simple expressions can deliver an approximate tortuosity scalar.As we treat in this framework a 2D case, such arguments can in principle deliver here as well a corresponding tortuosity scalar. 3However, the main advantage of our setting is that arrays of regularly-shaped microstructures are computationally tractable, hence we have a direct access to the entries of T simply by computing numerically (weak) solutions to families of cell problems (P ii ), i.e. ( 7)-( 9).This is simply a direct benefit of the fact that we could apply homogenization-type methods for the case of locally periodic porous media to derive the precise structure of the cell problems and corresponding effective diffusion tensor.

Brief review of known results
The precise structure of the problem studied here has been derived via homogenization methods for locally periodic media in Muntean and Nikolopoulos (2020), allowing for the possibility of microstructures to produce local clogging.Our original motivation to look into this class of problems is described in detail in the PhD thesis (Krehel 2014), where mathematical models have been devised to describe how hot colloids are supposed to be transported through heterogeneous, possibly chemically reactive, porous media; see, for instance, the related works (Ray 2013;Ray et al. 2012).
The overall problem we are considering in this work is then given by Problems (P i )−(P iv ) [or Eqs. ( 4)-( 14)].Given 0 < s < T , consider the time interval (0, s) ⊂ S. A weak solution to our problem is a set of functions (u, v, w, r ) with the regularity that satisfies (4)-( 14) in the standard weak Sobolev setting.
The model is a two-scale (macro-micro), quasilinear reaction-diffusion system posed in a 2D porous medium which undergoes microstrucural problems.The mathematical structure of this model is rather well-understood; see Eden et al. (2021) where the authors of this paper studied the weak solvability question.The analysis relies on a suitable application of Schauder's fixed point theorem which also provides a convergent algorithm for an iteration method to compute finite element/finite difference approximations of smooth solutions to our multiscale model (as suggested in Sect.4).The existence of solutions is guaranteed by Theorem 13 in Eden et al. (2021), while their uniqueness is expected to hold (cf.Remark 14 in loc.cit.).Interestingly, one can 3 Such tortuosity scalar is not supposed to comply with (17).It is rather seen as a ratio of two (gas) permeabilities (or of equivalent formulations); see e.g.Idris et al. (2022).
prove directly a continuous dependence result (in suitable norms) of the changes in radius with respect to small variations in the colloidal populations [cf.Lemma 7 in Eden et al. (2021)].The case of potential clogging, the primary focus in this work, was not considered in Eden et al. (2021).

Multiscale spatial discretization
We present a discretization of the problem consisting of Eqs.(4-14) for a two dimensional macroscopic domain Ω = (0, 1) × (0, 1), which is allowing us to solve numerically the problem.The aim is to focus on the behaviour of the weak solution when we have an initial random distribution of the porous radii as well as in the case when we have variations of the relative sizes of the parameters a i s and b i s as they arise in Eq. ( 11).The former set of experiments will give us an idea of how the porous radius evolves from an initial distribution towards clogging, while the latter can enlighten the significance of the mass interchange between colloidal and deposited species in the process.Naturally, Eq. ( 4) needs to be complemented with corresponding initial and boundary conditions.As we assume a two-dimensional macroscopic domain, we take without loss of generality x = (x 1 , x 2 ) ∈ [0, 1]×[0, 1].Furthermore, we set Robin boundary conditions at the one side of the square for the u i s, that is for the particular scenario when we have inflow of colloids in the domain from only one side of the domain.We impose Neumann boundary conditions (hence, a perfect reflection) for the rest of the boundary for (x 1 , x 2 ) such that 0 ≤ x 2 ≤ 1 with x 1 = 0, 1 or 0 ≤ x 1 ≤ 1 with x 2 = 0.As initial conditions, we set We also consider an arbitrarily prescribed initial distribution v 0 for the deposited species and r 0 for the radius given by Eqs. ( 12) and ( 14), respectively.Solution strategy We follow the same solution strategy as proposed in our recent work (Eden et al. 2021).Initially, we need to obtain a numerical approximation for the (weak) solution to the cell problems, i.e. we approximate the solutions to the Eqs.( 7)-( 9), and then determine the shape of the corresponding cell functions w 1 , w 2 posed in Y \ B(r ).Given that r ∈ [r 0 , 1/2], we take a partition of width δr , r a = r 0 , r 1 = r 0 + δr , . . ., r M 1 = 1/2.Next, we obtain a sequence of solutions for the cell problem (7) for each Y \ B(r i ) corresponding to the radius r i of the partition.These sequence of problems are solved with a finite element scheme [see Muntean 123and Nikolopoulos (2020) for the technical details].Having available the numerical evaluation of the cell functions w k (k = 1, 2) as approximate solutions to the cell problems ( 7)-( 9), the entries of the diffusion tensor (D i ) jk in Eq. ( 10), that is the the integral of the derivatives of the cell functions, together with the porosity φ times the constant diffusion coefficient d i , can now be calculated directly for each r i of the aforementioned partition.Consequently, for each r (x, t), and thus, for each Y \ B(r (x, t)), the corresponding value of (D i ) jk can be approximated directly via a standard linear interpolation.
As a next step, we proceed with the solving of the system of Eqs. ( 5)-( 7).To this end, we apply a finite difference scheme 4 to approximate the field Eq. ( 5) together with its boundary and initial conditions.
We consider a uniform partition of the domain Ω, Additionally, we take a partition of N T points of the time interval S = (0, T ), where T is the maximum time of the simulation, with step δt and t n = nδt, i = 0, . . .N T −1.
Let U i n 1 , 2 be the numerical approximation of the ith species of the solution of Eq. ( 4) at the point The discretization of the terms in (5) has the following form: Moreover, by applying a standard forward in time discretization scheme for the time derivative, we obtain the following finite difference scheme for tracking the evolution 4 The reader may wonder why we have chosen for a combined FEM-FD scheme to handle our problem.This is simply done for the convenience of the implementation.Other options are though available.For instance, to facilitate a rigorous numerical analysis of the discretization schemes, one would need to opt for a fully two-scale FEM approach, of for a fully two-scale FD approach.This will be done in a follow-up approach.
of both species u i and v, namely where β := N i=1 b i and and Finally, the approximate value r n 1 , 2 of the radius r is given by

Simulation case studies
Our main interest is to capture numerically the effect of the deposition process of the colloidal species around the solid cores of the cells.The effect is expected to be visible as soon variations in time of the radius r occur.
To illustrate numerically our model, we consider N = 3 mobile species u i and one immobile species v.Note though that N can be in practice much higher, e.g. for soils it is not unusual to reach N = 40.We take zero distributions as initial conditions (t = 0) for the colloidal populations, while we consider various specific initial distributions for the radius r .
Our model needs a quite large number of parameters.These dimensionless parameters were chosen in a range indicated by the model evaluation in Johnson and Elimelech (1995), Krehel et al. (2015) and Muntean and Nikolopoulos (2020).We take them as follows: We also take the function u b i to be defined as with u b 10 = 25 for t > 0. In this way we simulate an inflow of u 1 stronger in the center of the line segment x 2 = 0.Moreover, we let b r = 0.5, v(x 1 , x 2 , 0) = 0.In addition, we take as final simulation time T = 3 and set the remaining parameters to be M = 41, R := δt/δx 2 = 0.3.
Approximated concentration profiles with random initial radius distribution In this set of simulations a random initial distribution is assumed.We take r 0 to have random values of uniform distribution in the interval [0.01, 0.25].The form of the initial random distribution and corresponding porosity is presented in Fig. 3.
Next, in the first row of the panels shown in Fig. 4, concentration profiles of the colloidal population u 1 are plotted against space.A similar behavior of the concentration profiles is exhibited by the colloidal populations u 2 and u 3 as we can see in the second and third row of the same figure.For presentation purposes, we focus on the subdomain of Ω where variations are significant, that is on [0, 1]×[0, 0.5].We notice that there is no significant effect on the shape of the concentration profiles due to the randomness of the radius initial distribution.Small random perturbations smooth out with time.The inflow at the one segment of the square, x 2 = 0, causes an increase of the deposition, and consequently, also of the radius.This leads to a slowing down of the diffusion near the boundary.The latter also results in not having significant (and visible in Fig. 4) variations of the u i 's with respect to time, especially at the later stages of the simulation.This effect is also apparent in the next set of contour plots of the radius distribution as we can see in Fig. 5 (the perturbations due to randomness are visible).
In addition, we refer the reader to the next set of panels in Fig. 6.Here we show a sequence of line plots of the porosity and of the tortuosity given by Eq. ( 10) against time at specific points in the domain Ω.These specific points are: (x 1 , x 2 ) = (0, 5, 0), (0, 5, 0.25), (0, 5, 0.75), (1, 0.5).It us worth noting that in all these points, the porosity decreases with time and tends to stabilize at the end of the simulation, while the tortuosity (accounting for the specific internal geometry of our model) fluctuates around the value T ∼ 1.Moreover, it is instructive to see in Fig. 7 the distribution of the colloidal species u i , i = 1, 2, 3 against time for these specific spatial points.
Finally, in Fig. 8, we calculate the flow of one of the colloidal species, u 1 through a line segment of our domain vertical to the direction of the inflow from the boundary segment x 2 = 0. Specifically, we calculate the quantity (x 2 , t) = − initial bump at the first stage till t ∼ 1 and a stationary behaviour at larger times, decreasing as we move deeper in the domain.The later is due to clogging occurring near the boundary for larger times.Approximated concentration profiles with uniform initial radius distribution and variation on the relative size of a i and b i In the next series of simulations, we investigate the effect of the relative sizes of the parameters a i and b i on the flow of the colloidal species through a fixed frame with respect to time.Having as a baseline the values of the parameters chosen in the previous experiment we test the behaviour of Fig. 6 Line plots of the porosity (blue) and tortuosity (red) against time at the points (x 1 , x 2 ) = (0, 5, 0), (0, 5, 0.25), (0, 5, 0.75), (1, 0, 5).Due to the isotropic nature of the diffusivity in our case (direct result of our particular choice of microstructures), the tortuosity tensor is here a scaled unit matrix (colour figure online) the model for the cases that the a's or b's have different order of magnitude testing in this way the transfer mechanism of Henry's law.
More specifically, as we can see in Fig. 9, for a fixed point x 0 2 = 0.5 we calculate the total amount of u Then U 1 is plotted with respect to time for the case when (a) 3), a i = 30 b i , (blue line).The rest of the parameters are set to be the same as in the previous numerical experiment.These choises for the parameters a i , b i correspond to Damköhler numbers Da much larger as order of magnitude in cases (a) and (c) compared to case (b).
Additionally, in the next set of four panels presented in Fig. 10, we focus the attention on the flow of one of the colloidal species (u 1 ) through a line segment as in Fig. 8.We plot the flow (x 2 , t) at the points x 2 = 0, 0.25, 0.5, 0.75 for the three cases (a) a i b i , (blue line), (b) a i = b i , (red line), (c) a i b i , (orange line).In all these cases, it seems that the flow becomes stationary for larger times.Near the boundary (x 2 = 0, 0.25) when a i b i (case (c)) the flow is larger compared with cases (a) and (b) while it takes smaller values deeper in the square domain Fig. 7 Line plots of the colloidal species u i against time at the points (x 1 , x 2 ) = (0, 5, 0), (0, 5, 0.25), (0, 5, 0.75), (1, 0, 5) x 2 = 0.5, 0.75.In conclusion, the flow is choked for case (c) where a >> b, while in cases (a) and (b) we can observe some flow.

Concluding remarks and outlook
After reviewing our model equations, we did illustrate numerically clogging effects.The obtained results confirm our expectations from the modeling and analysis point of view, but are not yet verified experimentally.It would be very beneficial to get additional trust by getting the chance to compare our numerical results with pointwise measurements around detected clogging points.
Furthermore, additional work needs to be done regarding the numerical approximation of the target problem.A unified (two-scale) finite element approach can allow for a more efficient tackling of the problem for a large variety of choices of shapes of microstructures, possibly also allowing for a high performance computing approach like in Richardson et al. (2021).
It would be very interesting as well to consider random spacing of the porous radii.Playing numerically with variations in the positioning of the balls centering inside the unit cell adds extra difficulty in solving the cell problems, and consequently, in calculating efficiently the entries of the diffusion tensor.A successful handling of such aspect would enlarge considerably the types of porous structures where our methodology is in principle applicable.What concerns the current two-dimensional microscopic setting, it is worth noting that allowing for r to be greater than one [see Muntean and Nikolopoulos (2020)] would be an interesting addition both analytically and numerically to our current approach.Last but not least, from the viewpoint of practical applications, being able to handle a three-dimensional setting for the macroscopic problem coupled with a threedimensional setting for the microscopic problem is needed to obtain more realistic results near the fully clogged or partially clogged regions.

Fig. 3
Fig. 3 Left: random initial distribution of radii.Here r 0 takes random values following a uniform distribution in the interval [0.01, 0.25] at the first (left) graph.Right: corresponding porosity

Fig. 4 Fig. 5
Fig. 4 Concentration profiles at different time steps for the species u i , i = 1, 2, 3 in the case that we have random initial distribution of the radii.Note the accumulation localized around 0.5, more visible for the populations u 2 and u 3

Fig. 10
Fig. 10 Plot of (x 2 , t) for a a i b i , (orange lines), b a i = b i , (red lines), c a i b i , (blue lines), for x 2 = 0, 0.25, 0.5, 0.75 (colour figure online)