A multiscale quasilinear system for colloids deposition in porous media: Weak solvability and numerical simulation of a near-clogging scenario
AA multiscale quasilinear system for colloids deposition inporous media: Weak solvability and numerical simulation of anear-clogging scenario
Michael Eden † , Christos Nikolopoulos (cid:63) , and Adrian Muntean ‡† Zentrum f ¨ur Technomathematik, Department of Mathematics and Computer Science,University of Bremen, Germany (cid:63)
Department of Mathematics, School of Sciences, University of The Aegean, Greece ‡ Department of Mathematics and Computer Science, Karlstad University, SwedenFebruary 9, 2021
Abstract
We study the weak solvability of a quasilinear reaction-diffusion system nonlinearly coupledwith an linear elliptic system posed in a domain with distributed microscopic balls in D . The sizeof these balls are governed by an ODE with direct feedback on the overall problem. The systemdescribes the diffusion, aggregation, fragmentation, and deposition of populations of colloidalparticles of various sizes inside a porous media made of prescribed arrangement of balls. Themathematical analysis of the problem relies on a suitable application of Schauder’s fixed pointtheorem which also provides a convergent algorithm for an iteration method to compute finitedifference approximations of smooth solutions to our multiscale model. Numerical simulationsillustrate the behavior of the local concentration of the colloidal populations close to cloggingsituations. Keywords:
Colloidal transport and deposition, reactive porous media, weak solutions to stronglynonlinear parabolic systems, two-scale finite difference approximation, clogging
MSC2020:
We study a two-scale system modeling the effective diffusive transport as well as the aggrega-tion, fragmentation, and deposition of populations of colloidal particles inside porous media. Suchsituations arise, for instance, in membrane filtration scenarios [12, 26], papermaking [22], immobi-lization of colloids in soils [7], or transport of colloidal contaminants in groundwater [32].We are particularly interested in situations where micro-structural changes due to the depositionor dissolution of colloids are allowed to take place. This can locally change both the transportpatterns and storage capacity of the medium; see [4, 9, 14, 19, 28, 34] for related cases. This varietyof technological and natural processes is based on the transfer of colloidal particles from liquidsuspension onto stationary surfaces [16]. From this perspective, one can perceive that the porousmedia we are considering here behave like materials with reactive internal microstructures (see [8]for a periodic setting) and, based on [31], they are sometimes classified as media with distributedmicrostructures. Additional motivation for this work comes from our own research on reactive flowin porous media and is linked very much with the work of P. Ortoleva and J. Chadam (see e.g. [6]and follow up papers), but it is worth mentioning that quite related aspects arise in pharmacy andmedicine like drug delivery, thrombosis formation on arterial walls, evolution of Alzheimer’s disease.We refer the reader, for instance, to [5, 33, 13] for works in this direction.1 a r X i v : . [ m a t h . A P ] F e b enoting with u = ( u , ..., u N ) ( i = 1 , ..., N ) the molar concentrations of colloids of size i (with N ∈ N the maximal size), its time evolution can be modelled by a quasi-linear parabolic system inthe form of ∂ t u i − div( D i ( u ) ∇ u i ) = F i ( u ) , (1)where F i ( u ) accounts for the aggregation, segregation, and adsorption processes and D i ( u i ) thechanging permeability as consequence of the micro-structural changes (like clogging) inside theporous medium itself. While equation (1) is purely macroscopic, the computation of the effectivepermeability D i ( u i ) is done on the micro-scale therefore leading to the two-scale nature of our prob-lem. This system is a compact and abstract reformulation of a two-scale model for colloidal transportderived in [21] via asymptotic homogenization (more details are given in Section 2). Structurally sim-ilar (two-scale model with geometrical changes) models were investigated in, e.g., [11, 25].In this work, we take a D -cross-section of a porous medium and assume the solid matrix of thecross section to be made up of circles of not-necessarily uniform radius. The growth and shrinkageof these circles, which represent the underlying micro-structural changes of the porous medium, aremodelled via a scalar quantity governed by an additional ODE. For a similar geometrical setup see,e.g., [25]. The model and the resulting mathematical problem gets more complicated if we were toallow for more general geometries (e.g., evolving C -interfaces) that can not be represented by ascalar quantity like the radius in our setting. We treat our geometries in D mainly for the sake ofsimplicity of inequalities and transformations and also because the simulation work is easier to behandled in D compared to D , there is no fundamental element in the analysis that is sensitive todimensions (like Sobolev embeddings would be for example). As a consequence, the mathematicalanalysis part can be extended to D with suitable modifications on the upper and lower a priori bounds on the radii of the balls-like microstructure.The quasilinear structure of the problem together with the multiscale coupling is non-standard.Here, we point out that D i and F i are non-linear operators that are not defined via point wise eval-uation (in the sense of D i ( u )( t, x ) = D i ( t, x, u ( t, x )) ). In particular, it does not fit directly to theframework elaborated in, e.g., [2] and it requires an approach that utilizes the underlying couplingpresent in the model equations behind the abstract system. A similar two-scale problem allowing formicro-structural changes was investigated in [20].In Section 2 we explain our working strategy to prove the existence of weak solutions to theoverall problem. To keep things simple, we consider that the local porosity φ ( r ) does not degen-erate. Note however, that it is technically possible to include in the analysis simple degeneracies(like neighboring microstructures touching in single points [30]), a complete (local) clogging beinghowever out of reach. Besides the non-degeneracy of the effective parameter, another simplificationis included – the absence of the flow. Note that if the colloidal populations would be immersed in afluid flow, then, most likely, besides the balance equations of the linear momentum one would alsohave to take into account the charge transport taking place between oppositely charged populationsof particles; see e.g. [15, 27] for more information in this direction.The paper is organized as follows: In Section 2 we present the model and outline our strategyfor the analysis of our problem. We list the needed mathematical details of the problem so that wecan prove in Theorem 11 the existence of a weak solution. In Section 4, we solve numerically ourmultiscale quasilinear problem and discuss the obtained numerical results for realistic parameterregimes. We add in Section 5 a detailed discussion of the potential of our problem, expected results,and related aspects. In the following, let S = (0 , T ) be the time interval of interest and Ω ⊂ R a bounded Lipschitzdomain. In addition, let N ∈ N be a given number indicating the maximal possible size of anaggregate of colloid particle, where size refers to the number of primary particles making up the2ggregate. For each i = { , ..., N } , let u i : S × Ω → [0 , ∞ ) (we set u = ( u , ..., u N ) ) denote themolar concentration density of aggregates of size i at point x ∈ Ω at time t ∈ S . We take the function v : S × Ω → [0 , ∞ ) to represent the mass density of absorbed material (mass that is in the systembut currently not part of the diffusion and agglomeration process); this mass can be dissolved againby a Robin-type exchange allowing colloidal populations to re-enter the pore space. This processof absorption and dissolution is modelled in this context via an Robin-type exchange term (see e.g.[18]) in the form of πr − πr ( a i u i − β i v ) . Here, the radius function r : S × Ω → (0 , r max ) (for some r max > ) acts as a measure of the clogginess of the porous media and πr − πr is the ratio of the size of the boundary between fluidspace and pore to the pore volume.To describe the aggregation and fragmentation processes taking place inside the pore space ofthe medium, we use the Smoluchowski formulation (we point to [1] for a review) given here by R i ( u ) = 12 (cid:88) j + l = i γ jl u j u l − u i N − i (cid:88) j =1 γ jl u j . It is important to note that in the context of porous media the colloidal populations involve a finite sizechain of the cluster, i.e. there will be a population of N -mers where N takes the maximum clustersize. As a result, we deal with a finite sum here. Interestingly, for many applications a good choiceof such N is rather low; see e.g. [18].The diffusion-reaction system for the different aggregates is then given via ∂ t u i − div( D i ( r ) ∇ u i ) = R i ( u ) − πr − πr ( a i u i − β i v ) in S × Ω , (2a) − D i ( r ) ∇ u i · n = 0 on S × ∂ Ω , (2b) u i (0) = u i in Ω . (2c)The effective diffusion matrix (including diffusion, dispersion, and tortuosity effects) D i ( r ) ∈ R × can be calculated using any solution w k , k = 1 , , of the cell problem − ∆ w k = 0 in S × ( Y \ B ( r )) , (2d) −∇ w k · n = e k · n on S × ∂B ( r ) , (2e) y (cid:55)→ w ( · , · , y ) is Y -periodic . (2f)Here, Y = (0 , denotes the unit cell, B ( r ) is the closed ball with radius r and center point a = ( / , / ) , and e k the k -th unit normal vector. We have ( d i > are known constants) ( D i ) jk = d i φ ( r ) (cid:90) Y \ B ( r ) ( ∇ w k + e k ) · e j d z where φ ( r ) = − πr | Ω | denotes the porosity density of the medium. For more details regarding the cellproblem and the effective diffusivity, we refer to [21] where they are established via homogenization.Finally, the evolution of v is governed by an ODE parametrized in x ∈ Ω ∂ t v = N (cid:88) i =1 ( α i u i − β i v ) in S × Ω , (2g) v (0) = v in Ω . (2h)3nd the radius function is governed by the following ODE parametrized in x ∈ Ω ∂ t r = 2 πα N (cid:88) i =1 ( a i u i − β i v )) in S × Ω , (2i) r (0) = r in Ω . (2j)A possible initial choice for the radii r is depicted in Figure 1. We point out there also what willhappen 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 toone reported in [16]. -0.2 0 0.2 0.4 0.6 0.8 1 1.2-0.200.20.40.60.81 Schematic representation of r(x ,x ,0) -0.2 0 0.2 0.4 0.6 0.8 1 1.2-0.200.20.40.60.81 Schematic representation of r(x ,x ,T) Figure 1:
Example of r ( x , x , t = 0) with corresponding r ( x , x , t = T ) of the same simulation. Theparameter setting is as discussed in Figure 12. Regions with larger circles correspond to low porosity andpermeability. This accounts for the simple observation that the absorbed material leads to the clogging of thepore under the fundamental assumption of the growth of the radius is proportional to the amount ofmaterial that is absorbed. For a more concrete argumentation for this particular structure, we againpoint to [21].The overall problem we are considering in this work is then given by equations (2a–2j). Regard-ing our concept of a weak solution of this system:
Definiton 1 (Weak solution) . For a time interval (0 , s ) ⊂ S , a weak solution to the problem is givenby a set of functions ( u, v, w, r ) with the regularity u i ∈ L ((0 , s ); H (Ω)) ∩ L ∞ ((0 , s ) × Ω) such that ∂ t u i ∈ L ((0 , s ) × Ω) ,w ∈ L ((0 , s ) × Ω; H ( Y )) , v ∈ W , ((0 , s ); L (Ω)) , r ∈ W , ((0 , s ); L (Ω)) that satisfies equations (2a–2j) in the standard weak Sobolev setting. Solution strategy.
Without yet caring about regularity issues (like smoothness, integrability, mea-surability) and possible singularities, we want to suggest our solution strategy for the problem givenby equations (2a–2j) and show how it relates to the abstract quasi-linear PDE System 1.We start with a few comments regarding the particular structure of our problem where we referto the subproblems ( i ) - ( iv ) for u, w, v, r , viz.(A) The problem is strongly coupled: ( i ) depends on u, w, v, r , ( ii ) on w, r , ( iii ) on u, v , and ( iv ) on r, u, v . 4B) Problem ( i ) is parabolic in u , ( ii ) elliptic in w , ( iii ) and ( iv ) are first order ODEs in v and r .(C) Problem ( i ) is nonlinear in u and r , ( ii ) is nonlinear in r , and ( iii ) and ( iv ) are linear.(D) Problem ( ii ) is not a real free boundary problem, as the underlying domain Y \ B ( r ) dependsonly on ( t, x ) while the derivatives are w.r.t. y .As a consequence of points (A)–(D), a natural strategy is to first tackle the ODEs and to use them toinform the cell problem and the parabolic system. In the following, we outline the intermediate stepsinvolved in getting to the abstract fixed-point problem that will be the starting point for our analysis inTheorem 11:Step ( a ) : Looking at the linear ODE vor v (given by equations (2g) and (2h)), we find the char-acterization of v in terms of u via (setting b = (cid:80) Ni =1 β i ) v ( t, x ) = e − bt (cid:32) v ( x ) + N (cid:88) i =1 α i (cid:90) t e bτ u i ( τ, x ) d τ (cid:33) . With this in mind, we can eliminate v for u in our problem by setting v = L v ( u ) , where L v is theabstract solution operator for the v -problem.Step ( b ) : Similarly, looking at the second ODE (problem ( iv ) ), we have r ( t, x ) = r ( x ) + 2 πα N (cid:88) i =1 (cid:90) t ( a i u i ( τ, x ) − β i v ( τ, x )) d τ With this characterization, we can introduce the corresponding solution operator (cid:102) L r via r = L r ( u, v ) = L r ( u, L v ( u )) = (cid:102) L r ( u ) . Step ( c ) : Looking at the cell problem ( k = 1 , − ∆ w k = 0 in S × ( Y \ B ( r )) , −∇ w k · n = e k · n on S × ∂B ( r ) ,y (cid:55)→ τ ( · , · , y ) is Y -periodic , we expect to get solutions for every given r > such that B ( r ) ∩ ∂Y = ∅ . We introduce thecorresponding solution operator via w = L w ( r ) = (cid:16) L w ◦ (cid:102) L v (cid:17) ( u ) = (cid:102) L w ( u ) . Step ( d ) : Putting everything together, we can rewrite the parabolic problem ∂ t u i − div( D i ( r, w ) ∇ u i ) = R i ( u ) − πr − πr ( a i u i − β i v ) into ∂ t u i − div (cid:16) D i (cid:16)(cid:102) L r ( u ) , (cid:102) L w ( u ) (cid:17) ∇ u i (cid:17) = R i ( u ) − π (cid:102) L r ( u )1 − π ( (cid:102) L r ( u )) ( a i u i − β i L v ( u )) This highly nonlinear system of PDEs is now given only in terms of the unknown function u . On anabstract level, we therefore want to investigate parabolic system like ∂ t u i − div (cid:16)(cid:99) D i ( u ) ∇ u i (cid:17) = F i ( u ) in S × Ω , (3a) − (cid:99) D i ( u ) ∇ u i · n = 0 on S × ∂ Ω , (3b) u i (0) = u i in Ω (3c)where F i ( u ) = R i ( u ) − π (cid:102) L r ( u )1 − π ( (cid:102) L r ( u )) ( a i u i − β i L v ( u )) . The exact setting regarding function spaces will be settled in the following section.5
Analysis
In this section, we present the detailed fixed-point argument (as outlined in Section 2) for thenon-linear problem given via equations (3a–3c):The strategy of our proof is a three-step process:1) For a given function ˜ u (of sufficient regularity), we establish well-posedness and estimates forthe linear problem given by ∂ t u i − div (cid:16)(cid:99) D i (˜ u ) ∇ u i (cid:17) = F i (˜ u ) in S × Ω , (4a) − (cid:99) D i (˜ u ) ∇ u i · n = 0 on S × ∂ Ω , (4b) u i (0) = u i in Ω . (4c)This is established in Lemma 7.2) We show that there is a set such that the solution operator for equations (4a–4c) maps thatset onto itself, see Lemma 9. This result is local in time, since we need to keep t small in orderto control the norm of the solution.3) Finally, we employ Schauder’s fixed point theorem to establish the existence of at least onesolution, see Theorem 11.For some arbitrary (later to be fixed) M > and s ∈ (0 , T ) , let T s,M = { u ∈ L ((0 , s ) × Ω) N : (cid:107) u i (cid:107) ∞ ≤ M ( i = 1 , ..., N ) } . For ease of notation, for any given u of sufficient regularity we will write v u = L v ( u ) , r u = (cid:102) L r ( u ) , w u = (cid:102) L w ( u ) for the corresponding solution given for the particular subproblem and Y u = Y \ B ( r u ) . We start by collecting some important auxiliary results and estimates that will be needed in theconstruction of the actual fixed-point argument.Function Assumption Reason r / ≤ r ( x ) ≤ / Room for growth and shrinkage u i ≤ u i ( x ) ≤ M / Keeping the solution in T s,M v ≤ v ( x ) ≤ C v Bounding v u Table 1:
Assumptions regarding the initial data.
In a first step, we establish some sufficient conditions for the diffusivity matrix to not degenerate.Note that at this point it is not clear that this condition can be satisfied; this is shown in Lemma 3.
Lemma 2 (Diffusivity) . If u ∈ L ((0 , s ) × Ω) is chosen such that ≤ r u ≤ − ε for some small ε > , we find that (cid:99) D i ( u ) is symmetric and positive definite, i.e., (cid:99) D i ( u ) ξ · ξ ≥ c i | ξ | where theconstants c i > do not depend on u and ξ ∈ R . In addition, (cid:99) D i ( u ) ∈ L ∞ ((0 , s ) × Ω) .Proof. Its entries are given by ( (cid:99) D i ( u )) jk = d i φ ( r u ) (cid:90) Y u ( ∇ w u,k + e k ) · e j d z where r u = (cid:102) L iv ( u ) , Y u = Y \ B ( r u ) , and w u = ( w u, , w u, ) = L ii ( r u ) . The D i are symmetric since (cid:90) Y u ( ∇ w u,k + e k ) · e j d z = (cid:90) Y u ( ∇ w u,k + e k ) · ( ∇ w u,j + e j ) d z
6y way of w u,k solving the cell problem.Via that representation, non negativity is also straightforward to show (we refer to [23, Section12.5] for a similar argument) as long as φ ( r u ) is non negative. For the positivity, we have to ensurethat there is some c i > such that φ ( r u ) , | Y u | ≥ c i for all ( t, x ) ∈ S × Ω . Both hold true if r u isbounded away from / , i.e, if there is some ε > such that r u ≤ − ε for all ( t, x ) ∈ S × Ω .Now, regarding the boundedness of D i , we first see that | φ ( r u ) | ≤ | Ω | − when ≤ r u ≤ − ε is satisfied. Due to | Y u | ≤ | Y | = 1 , boundedness of D i is clear.In the following, we will try to establish sufficient conditions for a function u ∈ L ((0 , s ) × Ω) toguarantee that the condition r u ≤ − ε is met. Setting a u ( t, x ) = 2 πα N (cid:88) i =1 ( a i u i − β i v u ) , we get r u ( t, x ) = r ( x ) + (cid:90) t a u ( τ, x ) d τ. (5) Lemma 3 (Bounds for r ) . If M, ε , ε > satisfy M (cid:16) e bt − (cid:17) ≤ b παa min { − r − ε , inf r − ε − παbt sup v } (6) for all t ∈ (0 , s ) , it holds ε ≤ r u ≤ − ε for all u ∈ T s,M .Proof. For every u ∈ T s,M , we find that v u ( t, x ) = e − bt (cid:32) v ( x ) + N (cid:88) i =1 a i (cid:90) t e bτ u i ( τ, x ) d τ (cid:33) . As a consequence, − ab M ( e bt − ≤ v u ( t, x ) ≤ v ( x ) + ab M ( e bt − . This implies a u = 2 πα N (cid:88) i =1 ( a i u i − β i v u ) ≤ πα (cid:16) aM + aM ( e bt − (cid:17) = 2 παaM e bt as well as a u ≥ − πα (cid:16) aM + bv ( x ) + aM ( e bt − (cid:17) = − πα (cid:16) bv ( x ) + aM e bt (cid:17) . Therefore, inf r − πα (cid:18) tb sup v + aMb ( e bt − (cid:19) ≤ r u ( t, x ) ≤ sup r + 2 π αaMb ( e bt − . As a consequence, ε < r u < − ε can be ensured by the following two conditions: M (cid:16) e bt − (cid:17) ≤ b παa (1 − r − ε ) ,M (cid:16) e bt − (cid:17) ≤ b παa (inf r − ε − παbt sup v ) . emark 4. The condition 6 required in Lemma 3 can always be met (over some possibly smalltime interval (0 , s ) ) for M, ε , ε small enough as long as the initial radius distribution satisfies ε < r ( x ) < − ε . Connecting Lemma 3 with Lemma 2 leads to well behaved diffusivities for u ∈ T s,M .The additional bound from below in the form of ε is needed for the transformation for the cell problemfor w k . Now, looking at the r.h.s. of our reaction diffusion equation, we have for u ∈ T s,M (setting γ = max i,j γ ij ): − M γ (cid:18) N − k + 12 (cid:19) ≤ R k ( u ) ≤ M γ (cid:18) N − k + 12 (cid:19) (1 ≤ k ≤ N ) . (7)Due to r u ≤ / and πr u − πr u ≤ π − π / ≤ we arrive at πr u − πr u ( a i u i − β i v u ) ≤ (cid:16) a i M + ab β i M ( e bt − (cid:17) , (8)and πr u − πr u ( a i u i − β i v u ) ≥ − (cid:16) a i M + β i (cid:16) v ( x ) + ab M ( e bt − (cid:17)(cid:17) . (9)As a consequence, for every u ∈ T s,M , we find that F i ( u ) ∈ L ∞ ( S × Ω) for all i = 1 , ..., N . Inparticular, we find that sup {(cid:107) F i ( u ) (cid:107) ∞ : u ∈ T s,M } = C (10)where the constant C depends only s, M . Lemma 5 (Estimates for the radius) . For u (1) , u (2) ∈ T s,M let r (1) , r (2) be the corresponding solu-tions of the radius ODE problem. Then, (cid:12)(cid:12)(cid:12) r (1) − r (2) (cid:12)(cid:12)(cid:12) ≤ C (cid:90) t (cid:18)(cid:12)(cid:12)(cid:12) u (1) − u (2) (cid:12)(cid:12)(cid:12) + (cid:90) τ e bs (cid:12)(cid:12)(cid:12) u (1) − u (2) (cid:12)(cid:12)(cid:12) d s (cid:19) d τ. where the constant C > is independent of the particular choice of u ( k ) ( k = 1 , )Proof. The radius ODE can be solved by integration ( k = 1 , ): r ( k ) ( t, x ) = r ( x ) + 2 πα N (cid:88) i =1 (cid:90) t a i u ( k ) i ( τ, x ) − β i v ( k ) ( τ, x ) d τ where v ( k ) are given via v ( k ) ( t, x ) = e − bt (cid:32) v ( x ) + N (cid:88) i =1 a i (cid:90) t e bτ u ( k ) i ( τ, x ) d τ (cid:33) . Consequently, we can estimate (cid:12)(cid:12)(cid:12) r (1) − r (2) (cid:12)(cid:12)(cid:12) ≤ πα N (cid:88) i =1 (cid:90) t a i (cid:12)(cid:12)(cid:12) u (1) i − u (2) i (cid:12)(cid:12)(cid:12) + β i N (cid:88) j =1 a j (cid:90) τ e bs (cid:12)(cid:12)(cid:12) u (1) j − u (2) j (cid:12)(cid:12)(cid:12) d s d τ ≤ C (cid:90) t (cid:18)(cid:12)(cid:12)(cid:12) u (1) − u (2) (cid:12)(cid:12)(cid:12) + (cid:90) τ e bs (cid:12)(cid:12)(cid:12) u (1) − u (2) (cid:12)(cid:12)(cid:12) d s (cid:19) d τ. C > is independent of the particular choice of u ( k ) ( k = 1 , ). Lemma 6 (Estimates for the cell problem) . Let ε ≤ r ≤ r ≤ / (1 − ε ) and let w ( i ) k , k, i = 1 , ,solve − ∆ w ( i ) k = 0 in S × Y ( i ) , −∇ w ( i ) k · n = e k · n on S × Σ ( i ) , (cid:90) Y ( i ) w ( i ) k ( y ) d y = 0 ,y (cid:55)→ w ( i ) k ( y ) is Y -periodic . Then, the following estimate holds: (cid:12)(cid:12)(cid:12)(cid:12)(cid:90) Y (1) ∇ w (1) k · e j d y − (cid:90) Y (2) ∇ w (2) k · e j d y (cid:12)(cid:12)(cid:12)(cid:12) ≤ C | r (1) − r (2) | , where the constant C > might dependent on e and e but not on the particular choice of r (1) and r (2) . Here, we have set Y ( j ) = Y \ B ( r ( j ) ) and Σ ( j ) = ∂B ( r ( j ) ) .Proof. We prove this statement in three steps. First, we introduce a coordinate transform that allowsus to compare the different solutions and, second, go on proving some important energy estimates.Finally, we use these energy estimates to proof the desired result.
Step1: Transformation:
We set a = ( / , / ) and introduce the transformation ξ : Y → Y givenby ξ ( y ) = y, | y − a | ≥ / , (1 − χ ( | y − a | )) y + χ ( | y − a | ) (cid:0) r (1) / r (2) ( y − a ) + a (cid:1) , r (2) ≤ | y − a | ≤ / , r (1) / r (2) ( y − a ) + a, | y − a | ≤ r (2) Here, χ : [ r (2) , / ] → [0 , is a smooth cut-off function with compact support (i.e., χ ∈ C ∞ ( r (2) , / ) )satisfying χ ( r (2) ) = 1 , χ ( / ) = 0 , as well as − / ε ≤ χ (cid:48) ( z ) ≤ . As a result, ξ is a smooth functionas well and satisfies ξ ( Y (2) ) = Y (1) and n Σ (1) ( ξ ( y )) = n Σ (2) ( y ) for all y ∈ Σ (2) . Y (2) B ( r (2) ) Y (1) B ( r (1) ) ξ : Y (2) → Y (1) Figure 2:
Sketch of the transformation connecting reference cells for different radii r (1) and r (2) . Calculating the Jacobi matrix for ξ , we see that Dξ = I for | y − a | ≥ / and Dξ = (cid:0) r (1) / r (2) (cid:1) I for | y − a | ≤ r (2) . For the transition part, i.e., r (2) ≤ | y − a | ≤ / , we calculate9 y i ξ j ( y ) = ∂ y i (cid:2) y (cid:55)→ (1 − χ ( | y − a | )) y j + χ ( | y − a | ) (cid:0) r (1) / r (2) ( y j − / ) + / (cid:1)(cid:3) = δ ij (cid:18) r (1) / r (2) − χ ( | y − a | ) (cid:19) + (cid:0) r (1) / r (2) ( y j − / ) + / − y j (cid:1) y i − / | y − a | χ (cid:48) ( | y − a | ) As a consequence, we find that the Jacobian is given by the symmetric matrix Dξ ( y ) = a ( | y − a | ) (cid:18) (cid:19) + b ( | y − a | ) (cid:18) ( y − / ) ( y − / )( y − / )( y − / )( y − / ) ( y − / ) (cid:19) (11)where (setting r = r (2) − r (1) ≥ ) a ( z ) = (cid:18) − rr (2) χ ( z ) (cid:19) , b ( z ) = − χ (cid:48) ( z ) z rr (2) . We can calculate the determinant as det Dξ ( y ) = a ( | y − a | ) (cid:0) a ( | y − a | ) + b ( | y − a | ) (cid:0) y − y + y − y + 1 (cid:1) (cid:1) . Since a ( | y − a | ) > , b ( | y − a | ) ≥ and y − y + y − y + 1 > for all y = ( y , y ) ∈ Y , we findthat det Dξ ( y ) ≥ inf r (2) ≤| y − a |≤ / a ( | y − a | ) = (cid:32) r (1) r (2) (cid:33) . This shows that ε ≤ (cid:18) ε / (1 − ε ) (cid:19) ≤ det Dξ ( y ) ≤ which implies invertibility of Dξ . Step 2: Energy estimates . In the following, we set F ( y ) = Dξ ( y ) and J ( y ) = | det F ( y ) | . Westart with the the weak forms (cid:90) Y ( i ) ∇ w ( i ) k · ∇ η ( i ) d z = (cid:90) Σ ( i ) e k · n Σ ( i ) η ( i ) d σ (cid:16) η ( i ) ∈ H ( Y ( i ) ) , i = 1 , (cid:17) . We take the difference of these two weak forms: (cid:90) Y (1) ∇ w (1) k · ∇ η (1) d y − (cid:90) Y (2) ∇ w (2) k · ∇ η (2) d y = e k · (cid:20)(cid:90) Σ (1) n Σ (1) η (1) d σ − (cid:90) Σ (2) n Σ (2) η (2) d σ (cid:21) . and transform the surface integral on the right-hand side in order to arrive at (cid:90) Σ (1) n Σ (1) η ( i ) d σ − (cid:90) Σ (2) n Σ (2) η (2) d σ = (cid:90) Σ (2) n Σ (1) ( ξ ( y )) η (1) ( ξ ( y )) | det Dξ ( y ) | d σ − (cid:90) Σ (2) n Σ (2) η (2) d σ. By construction, we have n Σ (1) ( ξ ( y )) = n Σ (2) ( y ) for all y ∈ Σ (2) leading to (cid:90) Σ (1) n Σ (1) η (1) d σ − (cid:90) Σ (2) n Σ (2) η (2) d σ = (cid:90) Σ (2) (cid:18) η (1) ( ξ ( y )) det Dξ ( y ) − η (2) ( y ) (cid:19) n Σ (2) d σ = (cid:90) Σ (2) (cid:18) η (1) ( ξ ( y )) − η (2) ( y ) (cid:19) det Dξ ( y ) n Σ (2) d σ + (cid:90) Σ (2) (cid:18) det Dξ ( y ) − (cid:19) η (2) ( y ) n Σ (2) d σ (cid:90) Y (1) ∇ w (1) k · ∇ η (1) d y − (cid:90) Y (2) ∇ w (2) k · ∇ η (2) d y = (cid:90) Y (2) det Dξ ( Dξ ) − ∇ w (1) k ( ξ ) · ∇ η (1) ( ξ ) − ∇ w (2) k · ∇ η (2) d y and, as a consequence, (cid:90) Y (2) det Dξ ( Dξ ) − ∇ w (1) k ( ξ ) · ∇ η (1) ( ξ ) − ∇ w (2) k · ∇ η (2) d y = (cid:90) Σ (2) (cid:18) η (1) ( ξ ( y )) − η (2) ( y ) (cid:19) det Dξ ( y ) n Σ (2) d σ + (cid:90) Σ (2) (cid:18) | det Dξ ( y ) | − (cid:19) η (2) ( y ) n Σ (2) d σ. Now, choosing (cid:101) η = η (2) = (cid:101) w (1) k − w (2) k =: w k , this leads to (cid:107)∇ w k (cid:107) L ( Y (2) ) ≤ (cid:90) Y (2) (cid:12)(cid:12) det Dξ ( Dξ ) − − I (cid:12)(cid:12) (cid:12)(cid:12)(cid:12) ∇ (cid:101) w (1) k (cid:12)(cid:12)(cid:12) · |∇ w k | d y + (cid:90) Σ (2) (cid:12)(cid:12)(cid:12)(cid:12) det Dξ ( y ) − (cid:12)(cid:12)(cid:12)(cid:12) | w k | d σ. For y ∈ Σ (2) , i.e., | y − a | = r (2) , we have − det Dξ ( y ) = 1 − (cid:32) r (1) r (2) (cid:33) = ( r (2) ) − ( r (1) ) ( r (2) ) ≤ rr (2) . Now, for y ∈ Y with | y − a | ≥ / , we have det Dξ = 1 and Dξ = I and, in the case that r (2) ≤ | y − a | ≤ / , (cid:12)(cid:12) det Dξ ( Dξ ) − − I (cid:12)(cid:12) ≤ | det Dξ − || Dξ | + (cid:12)(cid:12) ( Dξ ) − − I (cid:12)(cid:12) | Dξ | + (cid:12)(cid:12) ( Dξ ) − − I (cid:12)(cid:12) Since | Dξ | ≥ det Dξ ≥ ε and − det Dξ ( y ) ≤ r / r (2) : (cid:12)(cid:12) det Dξ ( Dξ ) − − I (cid:12)(cid:12) ≤ r r (2) ε + (cid:18) ε (cid:19) (cid:12)(cid:12) ( Dξ ) − − I (cid:12)(cid:12) Finally, via (cid:12)(cid:12) ( Dξ ) − − I (cid:12)(cid:12) ≤ (cid:12)(cid:12) ( Dξ ) − (cid:12)(cid:12) | I − Dξ | ≤ ε | I − Dξ | we arrive at (looking at equation (11)) (cid:12)(cid:12) det Dξ ( Dξ ) − − I (cid:12)(cid:12) ≤ rr (2) (cid:18) ε + 2 ε + 1 + 1 ε r (2) (cid:19) Therefore we find that (cid:107)∇ w k (cid:107) L ( Y (2) ) ≤ C ( ε , ε ) r (cid:18)(cid:90) Y (2) (cid:12)(cid:12)(cid:12) ∇ (cid:101) w (1) k (cid:12)(cid:12)(cid:12) · |∇ w k | d y + (cid:90) Σ (2) | w k | d σ (cid:19) . Applying Poincar ´e’s inequality (possible due to the zero average condition) and the trace theoremleads to the energy estimate (cid:107) w k (cid:107) H ( Y (2) ) ≤ ˜ C ( ε , ε ) r, (12)11here the constant ˜ C ( ε , ε ) > is independent of r (1) and r (2) . Step 3: Proving the result . Using equation (12), we go on by estimating the following key expres-sion: (cid:12)(cid:12)(cid:12)(cid:12)(cid:90) Y (1) ∇ w (1) k · e j d y − (cid:90) Y (2) ∇ w (2) k · e j d y (cid:12)(cid:12)(cid:12)(cid:12) ≤ (cid:12)(cid:12)(cid:12)(cid:12)(cid:90) Y (1) ∇ w (1) k d y − (cid:90) Y (2) ∇ w (2) k d y (cid:12)(cid:12)(cid:12)(cid:12) = (cid:12)(cid:12)(cid:12)(cid:12)(cid:90) Y (2) det Dξ ( Dξ ) − ∇ (cid:101) w (1) k − ∇ w (2) k d y (cid:12)(cid:12)(cid:12)(cid:12) ≤ (cid:98) C ( ε , ε ) r, Now, let ε , ε , M ∗ , s ∗ > and initial conditions r , v be chosen such that ε ≤ r u ( t, x ) ≤ − ε for all ( t, x ) ∈ (0 , s ∗ ) × Ω and all u ∈ T s ∗ ,M ∗ (this is possible due to Lemmas 2 and 3). Also,let ≤ u i ( x ) ≤ M ∗ / . These choices imply F ( u ) = ( F ( u ) , ..., F N ( u )) ∈ L ∞ ((0 , s ∗ ) × Ω) N for all u ∈ T s ∗ ,M ∗ (see equation (10)). In the following, let s ∈ (0 , s ∗ ) and M ∈ (0 , M ∗ ) .We will now look at the linearized problem: For some ˜ u ∈ T s,M , we try to find a function u ∈ W ((0 , s ); H (Ω)) solving ∂ t u i − div (cid:16)(cid:99) D i (˜ u ) ∇ u i (cid:17) = F i (˜ u ) in S × Ω , (13a) − (cid:99) D i (˜ u ) ∇ u i · n = 0 on S × ∂ Ω , (13b) u i (0) = u i in Ω . (13c) Lemma 7 (Existence result for linearized problem) . For each ˜ u ∈ T s,M , there is a unique u ∈ W ((0 , s ); H (Ω)) solving the problem given by equations (13a–13c). Moreover, the following apriori estimates are satisfied (cid:107) ∂ t u (cid:107) L ( S ; H (Ω) ∗ ) + (cid:107) u (cid:107) L ∞ ((0 ,s ); L (Ω)) + (cid:107)∇ u (cid:107) L ((0 ,s ) × Ω) ≤ C (cid:32) (cid:107) u (cid:107) L (Ω) + N (cid:88) i =1 (cid:107) F i (˜ u ) (cid:107) L ∞ ((0 ,s ) × Ω) (cid:33) where the constant C > does not depend on ˜ u , s , and M . Please note that the above estimateimplies boundedness in W ((0 , s ); H (Ω)) as well.Proof. Since ˜ u ∈ T s,M , we have F i (˜ u ) ∈ L ∞ ((0 , s ) × Ω) ( i = 1 , ..., N ). Also, the diffusivity matrix (cid:99) D i (˜ u ) is uniformly positive definite (i.e., there is c i > such that (cid:99) D i (˜ u )( t, x ) ξ · ξ ≥ c i | ξ | for all ( t, x ) ∈ (0 , s ) × Ω and all ξ ∈ R ). Finally, as the D i are also bounded, the existence of a uniquesolution follows by standard theory of parabolic PDE.To search for the needed a priori estimates, we test the weak form with u i . Hence, we are led to (cid:107) u i ( t ) (cid:107) L (Ω) + 2 c i (cid:90) t (cid:107)∇ u i (cid:107) L (Ω) d τ ≤ (cid:107) u i (cid:107) L (Ω) + 2 (cid:90) t (cid:90) Ω | F i (˜ u ) u i | d x d τ ( t ∈ (0 , s )) . From here, summing over i = 1 , ..., N and applying Gr ¨onwall’s inequality leads to the desired esti-mate for u and ∇ u . Similarly, taking a test function ϕ ∈ L ((0 , s ); H (Ω)) such that (cid:107) ϕ (cid:107) ≤ , wefind that 12 ∂ t u i , ϕ (cid:105) L ((0 ,s ); H (Ω) ∗ ) ≤ (cid:90) Ω | F i (˜ u ) ϕ | d x + (cid:90) Ω | (cid:99) D i (˜ u ) ∇ u i ∇ ϕ | d x thus completing the estimate.With the solvability of the linarized problem established, we want to investigate under what cir-cumstances we can ensure that u ∈ T s,M as well; as this would then naturally lead to a potentialfixed-point scheme. As a first point, any (cid:101) u ∈ T s,M leads to a solution u ∈ W ((0 , s ); H (Ω)) whichagain leads to the corresponding solution operator L : T s,M → W ((0 , s ); H (Ω)) N . We now need to show, that s ∈ (0 , s ∗ ) and M ∈ (0 , M ∗ ) can be chosen such that L [ T s,M ] ⊂ T s,M .With the following lemma, we first establish L [ T s,M ] ⊂ L ∞ ((0 , s ) × Ω) N . Lemma 8 (Boundedness) . For every ˜ u ∈ T s,M , the solution of the linearized equation is boundedby − t ess sup( F i (˜ u )) − ≤ u i ≤ ess sup u i + t ess sup F i (˜ u ) . In particular, we have u ∈ L ∞ ((0 , s ) × Ω) N .Proof. By the linearity of the problem, we can decompose the solution u i = π i + ω i , where ∂ t π i − div (cid:16)(cid:99) D i (˜ u ) ∇ π i (cid:17) = 0 in S × Ω , − (cid:99) D i (˜ u ) ∇ π i · n = 0 on S × ∂ Ω ,π i (0) = u i in Ω , ∂ t ω i − div (cid:16)(cid:99) D i (˜ u ) ∇ ω i (cid:17) = F i (˜ u ) in S × Ω , − (cid:99) D i (˜ u ) ∇ ω i · n = 0 on S × ∂ Ω ,ω i (0) = 0 in Ω . Estimating the π i -problems via ( π i − L i ) + for L i = ess sup u i , we find that π i ≤ L i . Using Duhamel’sprinciple, we get ω i ( t, x ) = (cid:82) t h i ( τ, t, x ) d τ where the τ -parametrized function h i solves ∂ t h i − div (cid:16)(cid:99) D i (˜ u ) ∇ h i (cid:17) = 0 in S × Ω , − (cid:99) D i (˜ u ) ∇ h i · n = 0 on S × ∂ Ω ,h i (0) = F i (˜ u ( τ, · )) in Ω . This implies h i ≤ ess sup( F i (˜ u )) + and, as a consequence ω i ≤ t ess sup( F i (˜ u )) + . Finally, we have u i ≤ ess sup u i + t ess sup( F i (˜ u )) + . Now, since u i ≥ , we find that π i ≥ as well. Testing with ( h i + ess sup( F i (˜ u )) − ) − , we arrive at h i ≥ − ess sup( F i (˜ u )) − and, as a consequence ω i ≥ − t ess sup( F i (˜ u )) − . This shows u i ≥ − t ess sup( F i (˜ u )) − . In particular, we find that u i ∈ L ∞ ((0 , s ) × Ω) with (cid:107) u i ( t ) (cid:107) L ∞ (Ω) ≤ (cid:107) u i (cid:107) L ∞ (Ω) + t (cid:107) F i (˜ u )( t ) (cid:107) L ∞ (Ω) . Now, in order to get concrete bounds for the solution u = ( u , ..., u N ) , we have to take a closerlook at the right-hand sides: For the F i (˜ u ) , we have the estimates (given our assumptions on r , s ∗ ,and M ∗ and using equations (7–9)): F i (˜ u ) ≤ M (cid:18) M γ (cid:18) N + k + 12 (cid:19) + 15 (cid:16) a i + ab β i ( e bt − (cid:17)(cid:19) + 15 β i v ( x ) ,F i (˜ u ) ≥ − M (cid:18) M γ (cid:18) N + k + 12 (cid:19) + 15 (cid:16) a i + ab β i ( e bt − (cid:17)(cid:19) (cid:107) F i (˜ u )( t ) (cid:107) L ∞ (Ω) ≤ β i (cid:107) v (cid:107) L ∞ (Ω) + M (cid:18) M γ (cid:18) N + k + 12 (cid:19) + 15 (cid:16) a i + ab β i ( e bt − (cid:17)(cid:19) . (14)With this estimate at hand, we are now able to establish that L is a self-mapping for a suitable choiceof ( s, M ) . Lemma 9 (Fixed-point operator) . For any M ∈ (0 , M ∗ ) there is s ∈ (0 , s ∗ ) such that for every ˜ u ∈ T s,M the solution u of the linearized problem also satisfies u = L (˜ u ) ∈ T s,M .Proof. For any given M ∈ (0 , M ∗ ) , we find that lim t → t (cid:107) F i (˜ u ) (cid:107) ∞ → i = 1 , ..., N ) . uniformly for ˜ u ∈ T s,M (see inequality 14). As a consequence, it is possible to find s ∈ (0 , s ∗ ) such that s (cid:107) F i (˜ u ) (cid:107) ∞ ≤ M / for all i = 1 , ..., N and for all ˜ u ∈ T s,M . This implies u ∈ T s,M viaLemma 8.Please note that T s,M is a closed subset of L ((0 , s ) × Ω) . In the following lemma we investigatecontinuity of the fixed point operator Lemma 10 (Continuity) . The operator L : T s,M → L ((0 , s ) × Ω) is continuous with respect to the L -norm.Proof. Now let ˜ u, ˜ u ( k ) ∈ T s,M such that ˜ u ( k ) → ˜ u in L ((0 , s ) × Ω) for k → ∞ . In addition, let u = L (˜ u ) and u ( k ) = L (˜ u ( k ) ) ( k ∈ N ) be the corresponding unique solutions to the linearizedproblem (see Lemma 7).Now, the sequence u ( k ) is bounded in W ((0 , s ); H (Ω)) since ≤ ˜ u ( k ) ≤ M and the a prioriestimates given by Lemma 7. Since W ((0 , s ); H (Ω)) is a reflexive Banach space and since it iscompactly embedded in L ((0 , s ) × Ω) ( Lions-Aubin lemma ), there is a subsequence (for ease ofnotation, still denoted by u ( k ) ) and a limit function u ∗ such that u ( k ) converges to u ∗ strongly andweakly in L ((0 , s ) × Ω) and W ((0 , s ); H (Ω)) , respectively. Without loss of generality, we also have u ( k ) → u pointwise almost everywhere over (0 , s ) × Ω (possibly by choosing a further subsequence).In the following, we show continuity by establishing that u ∗ = u . The components of u ( k ) satisfy (for all ϕ ∈ H (Ω) and t ∈ (0 , s ) ) (cid:104) ∂ t u ( k ) i , ϕ (cid:105) H (Ω) ∗ + (cid:90) Ω (cid:99) D i (˜ u ( k ) ) ∇ u ( k ) · ∇ ϕ d x = (cid:90) Ω F i (˜ u ( k ) ) ϕ d x. Now, since ˜ u ( k ) → ˜ u in L ((0 , s ) × Ω)) , it holds (cid:90) Ω F i (˜ u ( k ) ) ϕ d x → (cid:90) Ω F i (˜ u ) ϕ d x ( ϕ ∈ H (Ω) , i = 1 , .., N ) . For the diffusion term, we take a look at (cid:90) Ω (cid:16)(cid:99) D i (˜ u ) ∇ u ∗ − (cid:99) D i (˜ u ( k ) ) ∇ u ( k ) (cid:17) · ∇ ϕ d x = (cid:90) Ω (cid:99) D i (˜ u ) ∇ (cid:16) u ∗ − u ( k ) (cid:17) · ∇ ϕ d x + (cid:90) Ω (cid:16)(cid:99) D i (˜ u ) − (cid:99) D i (˜ u ( k ) ) (cid:17) ∇ u ( k ) · ∇ ϕ d x. Due to this resulting statement:
Every subsequence has a further subsequence converging to u . u ( k ) to u ∗ in W ((0 , s ); H (Ω)) . Looking at the second term, we recall (cid:16)(cid:99) D i (˜ u ) − (cid:99) D i (˜ u ( k ) ) (cid:17) lm = d i (cid:18) φ ( r (0) ) (cid:90) Y (0) ( ∇ w (0) l + e l ) · e m d z − φ ( r ( k ) ) (cid:90) Y ( k ) ( ∇ w ( k ) l + e l ) · e m d z (cid:19) , which can be estimated using Lemmas 5 and 6 (cid:12)(cid:12)(cid:12)(cid:99) D i (˜ u ) − (cid:99) D i (˜ u ( k ) ) (cid:12)(cid:12)(cid:12) ≤ C (cid:90) t (cid:18)(cid:12)(cid:12)(cid:12) ˜ u − ˜ u ( k ) (cid:12)(cid:12)(cid:12) + (cid:90) τ e bs (cid:12)(cid:12)(cid:12) ˜ u − ˜ u ( k ) (cid:12)(cid:12)(cid:12) d s (cid:19) d τ. Here, we have used for the porosity that (cid:12)(cid:12)(cid:12) φ ( r (0) ) − φ ( r ( k ) ) (cid:12)(cid:12)(cid:12) ≤ π | Ω | (cid:12)(cid:12)(cid:12) r (0) − r ( k ) (cid:12)(cid:12)(cid:12) . Now, since ˜ u ( k ) → ˜ u almost everywhere over (0 , s ) × Ω , dominated convergence leads to (cid:90) Ω (cid:16)(cid:99) D i (˜ u ) − (cid:99) D i (˜ u ( k ) ) (cid:17) ∇ u ( k ) · ∇ ϕ d x → As a consequence, u ∗ = u . Theorem 11 (Existence) . The operator L : T s,M → L ((0 , s ) × Ω) has at least one fixed-point u ∗ ∈ W ((0 , s ); H (Ω)) .Proof. T s,M is a non-empty, closed, and convex subset of L ((0 , s ) × Ω) and L is continuous withrespect to the L ((0 , s ) × Ω) norm (Lemma 10). Moreover, we have L [ T s,M ] ⊂ T s,M via Lemma 9.Finally, since L [ T s,M ] ⊂ W ((0 , s ); H (Ω)) which is compactly embedded in L ((0 , s ) × Ω) by virtueof Lions-Aubin’s lemma, we can employ Schauder’s fixed point thorem to conclude the existence ofat least one fixed-point u ∗ ∈ W ((0 , s ); H (Ω)) ∩ T s,M . Remark 12.
Relying for instance on techniques from [10], we expect the weak solution given byTheorem 11 to be of higher regularity provided that data (boundary of Ω , initial conditions) aresufficiently smooth. This could change, however, if we were to allow actual clogging of the porousmedium. The aim is to solve numerically the two-dimensional macroscopic model problem for the speciesconcentration u i ( i ∈ { , . . . , N } ) and v . To focus the attention on physically relevant choices ofparameters, we use the setup described in [16]; see also [18, 21] for more details. Essentially, welook at a theoretical model describing the dynamics of colloid deposition on collector surfaces, whenboth inter-particle, and particle-surface electrostatic interactions are assumed to be negligible. Thenumerical range of the used parameters fit to the situations that can relate to the immobilization ofbio-colloids in soils.The simulation output we are looking after includes approximated space and time concentrationprofiles of colloidal populations, spatial distribution of microstructures for given time slices, and esti-mated amount of deposited colloidal mass. This information helps us detect in a posteriori way thelocations in Ω where deposition-induced clogging is likely to happen.15e have ∂ t u i ( x, t ) = D ijk ( x, t )∆ x u i ( x, t ) + R i ( u ) − L ( x, t ) A ( x, t ) ( a i u i ( x, t ) − β i v ( x, t )) , describing the diffusion of u i in the macroscopic domain Ω .The effective diffusion tensor has the form D ijk ( x, t ) = d i φ ( x, t ) τ jk ( x, t ) , where the entries τ jk ( x, t ) = (cid:90) Y ( x,t ) (cid:0) δ j,k + ∇ y j w k ( z, t ) (cid:1) dz, for all i = 1 , . . . , N , j, k = 1 , .In addition, the length L and area A functions related to the motion of the boundary (for r < / )are: L ( x, t ) = (cid:90) Γ( x,t ) ds = 2 πr ( x, t ) , A ( x, t ) = (cid:90) Y ( x,t ) dy = 1 − πr ( x, t ) , (in 2D) (15a) R i ( u ) = 12 (cid:88) i + j = k α i,j β i,j u i u j − u k ∞ (cid:88) i =1 α k,i β k,i u i . (15b)Moreover, the cell functions w := ( w ( x, y, t ) , w ( x, y, t )) , assumed to have constant mean,satisfy − ∆ y w i = 0 , i = 1 , in Y ( x, t ) , (15c) − n ( x, t ) · ∇ y w i = 0 , on ∂Y, − n ( x, t ) · ∇ y w i = n i ( x, t ) , on ∂B ( r ) . (15d)with Γ e := ∂Y being the boundary of the cell n ( x, t ) = ( n ( x, t ) , n ( x, t )) is the correspondingnormal vector.Equation (15a) needs to be complemented with corresponding initial and boundary conditions.In the sequel of this section, we focus the discussion on the case of a two dimensional macroscopicdomain, i.e. x = ( x , x ) ∈ [0 , × [0 , .We set Robin conditions at the one side of the square ∂u i ∂n ( x , , t ) + b r u i ( x , , t ) = (cid:26) u bi ( x ) > t ∈ [0 , t ] , t > t , , x ∈ [0 , , (15e)while we impose Neumann boundary conditions for the rest of the boundary ∂u i ∂n ( x , x , t ) = 0 , (15f)for ( x , x ) such that ≤ x ≤ with x = 0 , or ≤ x ≤ with x = 0 and with initial conditions u i ( x,
0) = u ai ( x ) ≥ . (15g)Moreover, we have ∂ t v ( x, t ) = N (cid:88) i =1 α i u i ( x, t ) − βv ( x, t ) , (15h)with some initial condition v ( x,
0) = v a ( x ) ≥ , (15i)and r ( x, t ) ∂ t r ( x, t ) = α (cid:32) N (cid:88) i =1 a i u i ( x, t ) − βv ( x, t ) (cid:33) L ( x, t ) , (15j)together with some initial distribution r ( x,
0) = r a ( x ) > , (15k)for x ∈ [0 , × [0 , . We discuss in Section 4.2 additional choices of suitable initial and boundaryconditions. 16 .2 Discretization schemes To treat problem (15) numerically, we need to obtain firstly a numerical approximation for thecell problems (15c) and determine the shape of the corresponding cell functions w , w posed in Y ( x, t ) .More specifically, we proceed for the various values of r , for r a ≤ r ( x, t ) ≤ / . We take apartition of width δr , r a = r , r = r + δr, . . . , r M = 1 / .Then since Y is determined as the area contained inside the square cell and outside the circleof radius r , we obtain a sequence of solutions for the cell problem (15c) for each Y i correspondingto the radius r i of the partition.We use a finite element scheme to solve these cell problems. To be precise, we use the MATLABfinite element package ” Distmesh ” (see details in [24]) to triangulate the domain Y i = Y ( r i ) .Furthermore, a solver has been implemented to handle this specific problem (equations (15c)); itworks in a similar fashion as applied in [21].In Figure 3, we illustrate the numerical solution for this problem for a particular choice of r i .Specifically, we choose to look at r i = . . -0.051 0.8 0.6 y y w Figure 3:
Numerical solution of the cell problem (15c) and specifically for w with r i = . . Having available the numerical evaluation of the cell functions w as approximate solutions to thecell problems (15c) and (15d), the entries of the diffusion tensor D ijk = (cid:82) Y ( x,t ) d i (cid:0) δ j,k + ∇ y j w k (cid:1) , i = 1 , . . . , N , j, k = 1 , can be calculated directly and for each ( x, t ) and consequently for thecorresponding value for r ( x, t ) and thus for Y ( x, t ) . Then the corresponding value of D ijk ( x, t ) isapproximated via linear interpolation.Next, we solve the system of equations (15a)-(15k). We use a finite difference scheme to solvethe two-dimensional version of the field equation (15a), together with its boundary and initial condi-tions. More specifically we consider a square domain Ω = [0 , × [0 , .For this purpose we implement a forward finite difference scheme and for this purpose initiallywe consider a uniform partition of the domain Ω , with x = ( x , x ) ∈ Ω , ≤ x ≤ , ≤ x ≤ ,of ( M + 1) × ( M + 1) points with spacial step δx = δx = δx , with x (cid:96) = (cid:96) δx , (cid:96) = 0 , , . . . M , x (cid:96) = (cid:96) δx , (cid:96) = 0 , , . . . M .Additionally, we take a partition of N T points in the time interval [0 , T ] , where T is the maximumtime of the simulation, with step δt and t n = nδt , i = 0 , . . . N T − .Let U in(cid:96) ,(cid:96) the numerical approximation of the species i of the solution of equation (15a) at thepoint ( x (cid:96) , x (cid:96) , t n ) of Ω T = Ω × [0 , T ] , that is u i ( x (cid:96) , x (cid:96) , t n ) (cid:39) U in(cid:96) ,(cid:96) . Moreover we denote17y D in(cid:96) ,(cid:96) the corresponding approximation of the diffusion coefficients D ijk ( x (cid:96) , x (cid:96) , t n ) (cid:39) D in(cid:96) ,(cid:96) and similarly by V in(cid:96) ,(cid:96) the approximation for the species v , v ( x (cid:96) , x (cid:96) , t n ) (cid:39) V n(cid:96) ,(cid:96) . Finite difference scheme for the model equations.
Initially we focus on the appropriate dis-cretization of the terms in (15a). For the spatial derivatives ∂∂x s (cid:16) D i ( x, t ) ∂u i ∂x s (cid:17) , where s = 1 , weapply a discretization of the form ∂∂x (cid:16) D i ( x, t ) ∂u i ∂x (cid:17) (cid:39) ∆ ( u i ( D i u ix )) x := δx (cid:104) D in(cid:96) + ,(cid:96) (cid:16) U in(cid:96) ,(cid:96) − U in(cid:96) ,(cid:96) δx (cid:17) − D in(cid:96) − ,(cid:96) (cid:16) U in(cid:96) ,(cid:96) − U in(cid:96) − ,(cid:96) δx (cid:17)(cid:105) ∂∂x (cid:16) D i ( x, t ) ∂∂x (cid:17) (cid:39) ∆ ( u i ( D i u ix )) x := δx (cid:104) D in(cid:96) ,(cid:96) + (cid:16) U in(cid:96) ,(cid:96) − U in(cid:96) ,(cid:96) δx (cid:17) − D in(cid:96) ,(cid:96) − (cid:16) U in(cid:96) ,(cid:96) − U in(cid:96) ,(cid:96) − δx (cid:17)(cid:105) D i(cid:96) + ,(cid:96) = D i(cid:96) ,(cid:96) +D i(cid:96) ,(cid:96) , D i(cid:96) − ,(cid:96) = D i(cid:96) ,(cid:96) +D i(cid:96) − ,(cid:96) , D i(cid:96) ,(cid:96) + = D i(cid:96) ,(cid:96) +D i(cid:96) ,(cid:96) , D i(cid:96) ,(cid:96) − = D i(cid:96) ,(cid:96) +D i(cid:96) ,(cid:96) − . Moreover we use a standard forward in time discretization for the time derivative and we concludewith a finite difference scheme of the form for the species u i ’s, U in +1 (cid:96) ,(cid:96) = U in(cid:96) ,(cid:96) + δt ∆ ( U i ( Du ix )) x + δt ∆ ( U i ( Du ix )) x + δtR in(cid:96) ,(cid:96) − δtF n(cid:96) ,(cid:96) and for the species v V n +1 (cid:96) ,(cid:96) = V n(cid:96) ,(cid:96) + δt N (cid:88) i =1 α i U in(cid:96) ,(cid:96) − βV n(cid:96) ,(cid:96) , where R in(cid:96) ,(cid:96) = 12 (cid:88) p + q = s α p,q β p,q U pn(cid:96) ,(cid:96) U pn(cid:96) ,(cid:96) − U sn(cid:96) ,(cid:96) ∞ (cid:88) p =1 α s,p β s,p U pn(cid:96) ,(cid:96) , and F n(cid:96) ,(cid:96) = L n(cid:96) ,(cid:96) A n(cid:96) ,(cid:96) (cid:0) a i U in(cid:96) ,(cid:96) − β i V n(cid:96) ,(cid:96) (cid:1) , are the approximations of the source terms at the point ( x (cid:96) , x (cid:96) , t n ) .In addition, the functions for the length L ( r ) and for the area A ( r ) , are approximated, for r ≤ / by the relations: L n(cid:96) ,(cid:96) = 2 πr n(cid:96) ,(cid:96) , A n(cid:96) ,(cid:96) = 1 − π ( r n(cid:96) ,(cid:96) ) , (in 2D) . Furthermore, we have the approximate value r n(cid:96) ,(cid:96) of the radius r given by r n +1 (cid:96) ,(cid:96) = r n(cid:96) ,(cid:96) + δt r n(cid:96) ,(cid:96) α (cid:32) N (cid:88) i =1 a i U in(cid:96) ,(cid:96) − βV n(cid:96) ,(cid:96) (cid:33) L n(cid:96) ,(cid:96) . In the first set of simulations we consider homogeneous Neumann boundary conditions at thethree edges of the square Ω , namely at x = 0 , x = 1 for ≤ x ≤ and at x = 1 , ≤ x ≤ .At the edge x = 0 , ≤ x ≤ we impose Robin boundary conditions given by equation (15e).That is we consider a scenario of having inflow at this side of Ω for a particular time period, [0 , t ] which stops after some time t , and we want mainly to observe the deposition process of the colloidspecies around the solid cores of the cells. The later can be apparent by the variation in time of theradius r . 18e take zero distributions as initial conditions ( t = 0 ) for the colloidal populations, while weconsider various specific initial distributions for the radius r .We consider N = 3 mobile species u i and one immobile species v . Our model needs aquite large number of parameters. We take them as follows: κ = 1 , ( d , d , d ) = ( . , . , . , ( a , a , a ) = ( . , . , . , ( β , β , β ) = (1 , , , α i,j = . , β i,j = 100 , i, j = 1 , . . . , u ia ( x ) =0 , v a ( x ) = 0 , r a ( x ) = . , ≤ x ≤ .Regarding the choice of boundary condition at ( x , , we take the function u bi to be defined as ( u b , u b , u b ) = ( u b x ( i − x ) , , with u b = 25 for t ∈ [0 , t ] and zero for t > t , with t = 2 . Moreover, we let b r = 0 . , v ( x , x ,
0) =0 , and r ( x , x ,
0) = 0 . .In addition, we take as final simulation time T = 3 and set the remaining parameters to be M = 41 , R := δt/δx = 0 . . Approximated concentration profiles.
In the first of the following graphs, i.e. in Figure 4, concen-tration profiles of the colloidal population u are plotted against space. Similar profiles are exhibitedby the other colloidal populations as well. As general rule, we keep the discussion about whathappens with u only as here the effects are more visible. This corresponds also to the physicalsituation when most of the mass is contained in the monomer population, while the amount of ob-servable dimer, trimer, 4-mer populations is considerably lower; see e.g. [18] and references citedtherein.In the first two frames we have t < t ; hence we can see that there is an inflow in Ω throughone edge and so we can observe the diffusion of u taking place in the x direction. In the last twoframes taken at times after t (hence here the inflow has stopped) we see that the concentration of u near the edge drops possibly due to an activation of the reaction mechanisms. Especially, thedeposition activates and consumes monomers initially involved in diffusion. x x t=0.75 u x x t=1.5 u u t=2.25 x x u t=3 x x Figure 4:
Concentration profiles at different time steps for the species u . In Figure 5, we present similar graph for the concentration of u . As expected, the behaviour is19imilar as for the species u . Moreover, for the third species u during the simulation we notice nodifference in its qualitative behaviour. x x t=0.75 u x x t=1.5 u u t=0.25 x x u t=3 x x Figure 5:
Concentration profiles at different time steps for the species u . Regarding the behaviour of the immobile species v pointed out in Figure 6, we observe an initialdistribution in the first two frames t = 0 . , t = 1 . , following the form of the mobile species u i and an increase inside the domain Ω . After the inflow stops, for instance, see the last two frames t = 1 , , t = 3 , the distribution of the mass of the deposited species appears to be stationary.Focusing now in the behaviour of r , we present in Figure 7 time frames of contour plots ofthe radius at times t i = 0 . , , , , , , . We observe the expected increase of the radius withrespect to time. Even for t > t = 2 , after the inflow has stopped to happen, we still have a slightincrease of the radius due to the accumulation of the immobile species around the spherical coresof the cells.As final remarks regarding this numerical experiment, the main observables u , u , u , and v are plotted in Figure 8 against time for fixed locations inside the domain Ω ; see specifically the points (0 , . , the center (0 . , . , (0 . , and at the corner (0 , . Approximations with non uniform initial radius.
In the following experiment we consider for thesame scenario of initial and boundary conditions, (15e), (15f), (15g), a non uniform distribution forthe initial values of the radius r = r ( x , x , . Specifically, we consider larger values of the radiusin the form of two peaks centered at the points (0 . , . and (0 . , . and with r a having the form r a = r c + r exp (cid:2) − c ( x − . − c ( x − . (cid:3) + r exp (cid:2) − c ( x − . − c ( x − . (cid:3) . In this context, we take r c = 0 . , r = 0 . , c = 60 so that the maximum radius at these twopoints is quite large but smaller than one ( max r ( x , x , (cid:39) . ) as it can be seen in the yellowarea shown in Figure 10. Here we also set M = 41 for the spatial partition and R = 0 . The restof the parameters values are the same as in the previous numerical experiment.The effect of the non-uniform initial radius distribution is apparent in the evolution of the speciesof the model; particularly, this non-uniformity effect can be traced back in the evolution of the popu-lation u as exhibited in Figure 9. 20 x t=0.75 x v x x t=1.5 v x x t=2.25 v x t=3 x v Figure 6:
Mass at different time steps for the deposited species v . Due to the inflow from the edge x = 0 , we have now high values in the u concentration aroundthis edge (yellow area) of the domain, while inside the domain we have lower value (blue areas);this behavior can be seen in the first two frames of the simulation ( t = 0 . , t = 1 . ). We notice agradual increasing perturbation of the symmetric form of u around the point (0 . , . due to thefact that, precisely at this point, we have large values of r . In the next frames, at ( t = 2 . , t = 3) and particularly at t = 2 . , we observe the concentration of u after the time that the inflow in thedomain has stopped ( t > t and ∂u i ∂n ( x , , t ) + b r u i ( x , , t ) = 0 ). The dominant mechanisms noware the diffusion and the surface reaction, i.e. the deposition of material around the cores of thecells. Thus we observe lower values of u (blue and green areas) around the points with larger r (close to the two initial peaks of r ) where there the material has been deposited and higher values(yellow areas) in between the aforementioned peak points where the values of r are smaller anddeposition is slower. Essentially due to the same mechanism, at the final frame t = 3 at the end ofthe simulation, the values of u decrease and tend to zero with slower speed within the area closeto the corner (0 , .In Figure 10, we present the contour plot of the initial value of r for this experiment. In Figure 11,we point out the spatial distribution of the radius r = r ( x , x , T ) , where T is the final time of thesimulation. In this case, we observe a behaviour consistent with what happens with the profile of thecolloidal population u towards the end of the simulation, i.e. around t = 3 . This effect is shown inFigure 9.Higher values of r equal to . , where clogging occurs, are taken in the lower part of the domainnear the edge x = 0 as well as in the neighbour of the points (0 . , . and (0 . , . ; observe theyellow areas in Figure 11. In the rest of the domain Ω the radius r attains lower values. This is inline with the observed behaviour of the concentration profiles of u around the end of the simulation.The evolution of the diffusivity during the experiment is also apparent in Figure 12. We noticeinitially low values of it in the areas (blue regions) around the two peaks and higher values in theintermediate area (yellow region), in the first frame for t = 0 . . As r gradually increases thecorresponding areas with low diffusivity expand as we can see in the second and third frame for t = 1 . , . , and finally, also for t = 3 at the end of the simulation where we obtain the final map21 = . . . . . . x x t=1.5 . . . . . . x x t=2.25 . . . . . . . x x t=3 . . . . . . . x x Figure 7:
Contour plots of the radius r = r ( x , x , t i ) for the time steps t i = 0 . , , , , , , . of the diffusivity. This contains also information on the tortuosity of the material. The latter frame isin fact a ”reverse” image of Figure 11 as very low values of D are linked to clogging around the blueareas where r is large.It is worthwhile to note that the spatial distribution of the balls-like microstructure that corre-sponds to the vizualization shown in Figure 11 of the effective transport coefficient is pointed out inFigure 1. The unavoidable occurrence of clogging is pointed out in all these representations. We have proven the existence of a weak solution to a specific coupled multiscale quasilinearsystem describing the diffusion, aggregation, fragmentation, and deposition of populations of col-loidal particles in porous media. The structure of the system was originally derived in [21] and wekept it here.Tracking numerically the x -dependence in the shape of the microstructures rises serious compu-tational problems especially in 3D or even in 2D when working with low-regular shapes. Because ofthe strong separation between the macroscopic length scale and the microscopic length scale, suchsetting is parallelizable; see [29] for a prestudy in this direction done for a micro-macro reaction-diffusion problem with x -dependent microstructure arising in the context of transport of nutrients inplants. The approach used in [29] is potentially applicable here as well. Moreover, what concernsthe discretization techniques used in this framework, a more advanced finite difference scheme,such as an appropriate version of Du Fort Frankel scheme, can give in principle more flexibility andaccuracy in the numerical computations, e.g. by allowing larger time steps.Our multiscale model can allow for further relevant extensions in at least twofold direction:(1) For instance, a particularly interesting development would be to allow for some amount ofstochasticity in the balance laws. In this spirit, the ODE for the growth of the balls induced by thedeposition of the species v could have not only a random distribution of initial positions but alsosome suitably scaled ”Brownian noise” in the production term mimicking an additional contribution This is tractable with the current form of the model. igure 8: Concentration profiles of the species u i , v versus time at different spatial points in the squaredomain. eventually due to a non-uniform deposition of colloids on the boundary of the microstructures (com-pare with the setting from [3]). The difficulty in this case is that, due to the strong coupling in thesystem, the overall problem becomes a quasilinear SPDE, which is much more difficult to handlemathematically and from the simulation point of view compared with our current purely deterministicsetting.(2) Another development that would be interesting to follow in the deterministic setup is to at-tempt a computational efficient hybrid-type modeling. In this context, one idea would be to couplecontinuum population models for colloidal dynamics with discrete network models describing themechanics of the underlying material (see e.g. the approach proposed in [17] having paper as tar-get material). Relevant questions would be: What is the counterpart of our equation for the radiusgrowth of a ball B ( r ) , when the ball is replaced by a point? How does ”continuum” deposition takeplace on ”discrete” fixed locations? Are points able to absorb matter in D and D ?We expect that the non-standard type of couplings suggested in (1) and (2) (i.e. deterministic-stochastic and continuum-discrete) can potentially be posed in terms of measured-valued balanceequations. We will investigate some of these ideas in follow-up works.23 =0.75 . . . . . . x x t=1.5 . . . . . . x x t=2.25 . . . . . . . . x x t=3 . . . . . . x x Figure 9:
Contour plots at different time steps for the concentration of the species u for the case of nonuni-form initial radius distribution. . . . . . . . . . . . . . . . . . . . . r=0.42 r=0.42 Figure 10:
Contour plot steps for initial radius distribution r a = r ( x , x , . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure 11:
Contour plot for the radius distribution r a = r ( x , x , T ) at the end of the simulation. t=0.75 . . . . . . . . x x t=1.5 . . . . . . . . . . . . . x x t=2.25 . . . . . . . . . . . . . x x t=3 . . . . . . . . . . . . x x Figure 12:
Contour plots at different time steps for the effective diffusivity D ( x, t ) for the case of nonuniforminitial radius distribution. cknowledgments AM is partially supported by the grant VR 2018-03648 ”
Homogenization and dimension reductionof thin heterogeneous layers ”. We thank R. E. Showalter (Oregon) and O. Richardson (Karlstad) foruseful discussions on closely related topics.
References [1] D. J. Aldous. Deterministic and stochastic models for coalescence (aggregation and coagula-tion): A review of the mean-field theory for probabilists.
Bernoulli , pages 3–48, 1999.[2] H. W. Alt and S. Luckhaus. Quasilinear elliptic-parabolic equations.
Math. Z. , 183:311–341,1983.[3] H. Bessaih, Y. Efendiev, and R. F. Maris. Stochastic homogenization of a diffusion-reactionmodel.
DCDS - Series A , 39(9), 2019.[4] G. Boccardo, E. Crevacore, R. Sethi, and M. Icardi. A robust upscaling of the effective particledeposition rate in porous media.
Journal of Contaminant Hydrology , 212:3–13, 2018.[5] G. Bonacucina, M. Cespi, M. Misici-Falzi, and G. F. Palmieri. Colloidal soft matter as drugdelivery system.
Journal of Pharmaceutical Sciences , 89(1):1–42, 2009.[6] J. Chadam and P. Ortoleva. A mathematical problem in geochemistry: The reaction-infiltrationinstability.
Rocky Mountain J. Math. , 21(2), 1991.[7] Y. Chen, J. Ma, X. Wu, L. Weng, and Y. Li. Sedimentation and transport of different soil colloids:Effects of Goethite and humic acid.
Water , 12:980, 2020.[8] C. Conca, J. I. Diaz, and C. Timofte. On the homnogenization of a transmission problem arisingin chemistry.
Romanian Reports in Physics , 56(4):613–622, 2004.[9] M. P. Dalwadi, Y. Wang, J. R. King, and N. P. Minton. Upscaling diffusion through first-ordervolumetric sinks: A homogenization of bacterial nutrient uptake.
SIAM J. Appl. Math. , 78:1300–1329, 2018.[10] E. DiBenedetto.
Degenerate Parabolic Equations . Springer Verlag, Berlin, 1993.[11] M. Eden. Homogenization of a moving boundary problem with prescribed normal velocity.
Adv.Math. Sci. Appl , 28(2):313–341, 2019.[12] A. Fasano and A. Mikeli´c. On the filtration through porous media with partially soluble perme-able grains.
Nonlinear differ. equ. appl. , 7:91–105, 2000.[13] B. Franchi, M. Heida, and S. Lorenzani. A mathematical model for Alzheimer’s disease: Anapproach via stochastic homogenization of the Smoluchowski equation.
Communications inMathematical Sciences , 18(4):1105–1134, 2020.[14] B. Hallak, E. Specht, F. Herz, R. Gr ¨opler, and G. Warnecke. Influence of particle size distributionon the limestone decomposition in single shaft kilns.
Energy Procedia , 2017.[15] R. J ¨ager. Erosion and deposition in porous media. Master’s thesis, ETH Z ¨urich, Switzerland,2020.[16] P. R. Johnson and M. Elimelech. Dynamics of colloid deposition in porous media: Blockingbased on random sequential adsorption.
Langmuir , 11(3):801–812, 1995.[17] G. Kettil, A. M ˚alqvist, A. Mark, M. Fredlund, K. Wester, and F. Edelvik. Numerical upscaling ofdiscrete network models.
BIT Numerical Mathematics , 60:67–92, 2020.2618] O. Krehel, A. Muntean, and P. Knabner. Multiscale modeling of colloidal dynamics in porousmedia including aggregation and deposition.
Advances in Water Resources , 86:209–216, 2015.[19] J. Maes and C. Soulaine. A unified single-field volume-of-fluid-based formulation for multi-component interfacial transfer with local volume changes.
Journal of Computational Physics ,402:109024, 2020.[20] S. A. Meier. Global existence and uniqueness of solutions for a two-scale reaction-diffusionsystem with evolving pore geometry.
Nonlinear Anal. , 71(1-2):258–274, 2009.[21] A. Muntean and C. Nikolopoulos. Colloidal transport in locally periodic evolving porous media—An upscaling exercise.
SIAM J. Appl. Math. , 80(1):448–475, 2020.[22] A. Nyfl ¨ott, E. Moons, C. Bonnerup, G. Carlsson, L. J ¨arnstr ¨om, and M. Lestelius. The influenceof clay orientation in dispersion barrier coatings on oxygen permeation.
Applied Clay Science ,126:17–24, 2016.[23] G. Pavliotis and A. Stuart.
Multiscale Methods : Averaging and Homogenization . Springer, NewYork, 2008.[24] P. O. Persson and G. Strang. A simple mesh generator in MATLAB.
SIAM Review , 46(2):329–345, 1998.[25] M. A. Peter. Coupled reaction-diffusion processes inducing an evolution of the microstructure:analysis and homogenization.
Nonlinear Anal. , 70(2):806–821, 2009.[26] G. Printsypar, M. Bruna, and I. Griffiths. The influence of porous-medium microstructure onfiltration.
Journal of Fluid Mechanics , 86:484–516, 2019.[27] N. Ray.
Colloidal Transport in Porous Media-Modeling and Analysis . PhD thesis, University ofErlangen, Germany, 2013.[28] N. Ray, A. Rupp, R. Schultz, and P. Knabner. Old and new approaches predicting the diffusionin porous media.
Transport in Porous Media , 124:803–824, 2018.[29] O. M. Richardson, O. Lakkis, A. Muntean, and C. Venkataraman. Parallel two-scale finite el-ement implementation of a system with varying microstructures. Technical report, KarlstadUniversity, Sweden, 2021.[30] R. Schulz, N. Ray, F. Frank, H. Mahato, and P. Knabner. Strong solvability up to clogging ofan effective diffusion-precipitation model in an evolving porous medium.
European Journal ofApplied Mathematics , pages 1–29, 2016.[31] R. E. Showalter. Distributed microstructure models of porous media. In U. Hornung, editor,
Flow in Porous Media , pages 153–163. Oberwolfach, 1992.[32] N. Suciu, F. A. Radu, S. Attinger, L. Sch ¨uller, and P. Knabner. A Fokker-Planck approach forprobability distributions of species concentrations transported in heterogeneous media.
Journalof Computational and Applied Mathematics , 289:241–252, 2015.[33] C. Valladolid, M. Martinez-Vargas, N. Sekhar, F. Lam, C. Brown, T. Palzkill, A. Tischer, M. Auton,K. V. Vijayan, R. E. Rumbaut, T. C. Nguyen, and M. A. Cruz. Modulating the rate of fibrinformation and clot structure attenuates microvascular thrombosis in systemic inflammation.
Blood Advances , 4(7):1340–1349, 2020.[34] T. L. van Noorden and A. Muntean. Homogenisation of a locally periodic medium with areas oflow and high diffusivity.