A Deep Bed Filtration Model of Two-Component Suspension in Dual-Zone Porous Medium

: In the paper, a mathematical model for the ﬁltration of two-component suspensions in a dual-zone porous medium is considered. The model consists of the mass balance equations, the kinetic equations for active and passive zones of porous medium for each component of the suspension and Darcy’s law. To solve the problem, a numerical algorithm for computer experimentation is developed on the basis of ﬁnite difference method. Based on numerical results, the main characteristics of suspension ﬁltration in a porous medium are established. Inﬂuences of model parameters on transport and deposition of suspended particles of two-component suspension in porous media are analysed. It is shown that the polydispersity of suspension and multistage nature of the deposition kinetics can lead to various effects that are not characteristic for the transport of one-component suspensions with one-stage particle deposition kinetics. In particular, in distribution of the concentration of suspended particles in a moving ﬂuid non-monotonic dynamics are obtained at individual points in the medium. is shown that at the points of the medium near to the input section, the concentration of deposited particles can reach partial capacities in the passive zone.


Introduction
Many scientific and practical studies in such fields as fluid and gas mechanics, hydrogeology, and ecology lead to problems in the theory of filtration of inhomogeneous liquids. Filtration models of heterogeneous liquids in porous media and methods for their numerical research serve as the theoretical basis for the processes of treating drinking and wasting water [1], increasing the efficiency of oil and gas production [2], studying colloid transport in chemical engineering [3], etc. Therefore, the development of adequate mathematical models for filtering heterogeneous liquids and solute transport in porous media, as well as their study are important tasks of mathematical modeling.
As mentioned in [4], during the filtration of inhomogeneous liquids particles can deposit either on the filter surface or in the filter on the grains of the medium. Hence, by the method of particle deposition, filtration of suspension is divided into two groups: cake filtration and deep bed filtration [4]. Cake filtration usually is used in the cases when the suspension has large particles with high concentration whereas deep bed filtration is used for low concentration suspension [5]. In some cases, these two types are used together, when the suspension has particles with different sizes [6].
Various approaches have been developed to simulate the particle transport and deposition process, which can be arbitrarily divided into the so-called network models, particle trajectory analysis models, and empirical models [7].
A simulation of particle deposition in a porous medium based on the network representation was proposed in [8]. Network models were used to describe the change in permeability during particle deposition, using the method of random movement of a particle through a lattice with a given probability of particle capture. Approaches using the theory of an effective medium were developed in [9][10][11][12]. Most studies considered two mechanisms of particle deposition: blockage of small capillaries by large particles and the deposition of small particles on the walls of pore channels with the formation of a sediment layer. Taking these mechanisms into account, attempts were made in [9] to predict the behaviour of the concentration of deposited particles, the filtration coefficient, and the change in permeability.
In the models of particle trajectory analysis [13,14], the skeleton of a porous medium is represented as a set of single sections, which can be spherical, cylindrical, or other simple shapes. The trajectory of each particle during the flow of fluid through separate sections of the medium is characterized by a stream function determined from the problem of flowing around, for example, a sphere or cylinder under the condition that various particles act on the particle (gravitational, hydrodynamic, van der Waals forces, etc.) [15]. Trajectory analysis leads to significant progress in understanding the deposition process at the individual pore level, although it is rather complicated. This approach allows to quite successfully predicting the value of the filtration coefficient and the change in the concentration of deposited particles. However, it has some limitations associated with determining changes in the permeability of a porous medium.
In empirical (phenomenological) models [7,[16][17][18][19], the process of filtering a suspension in a porous medium is described using equations of continuum mechanics, which are supplemented by empirical relations. In classical deep bed filtration models system of partial differential equations describing the process consists of mass balance equation end kinetic equation of deposition [16,17]. In one-dimensional case, mass balance equation takes the form [20][21][22][23][24] ∂c ∂t + v ∂c ∂x where c is the concentration of suspended particles in the suspension, ρ is the concentration of deposition, D is diffusion (hydrodynamic dispersion) coefficient, v is the velocity, t is time, x is coordinate. In some cases the diffusive (dispersive) fluxes are ignored for large particles and in large scale approximation, since they are insignificantly smaller than the advective particle fluxes [25,26], so mass balance negligible diffusion takes the form [9,10,25] ∂c ∂t + v ∂c ∂x Mainly, kinetic equations can be divided into two groups: irreversible and reversible deposition according to their characteristics. In irreversible deposition [16,17], deposited particles never get back to flow again and kinetic equation is in the form where β is the coefficient characterizing intensity of deposition, which may be changed during the filtration. In reversible deposition [18][19][20][21][22][23][24][25][26], both attachment and detachment of suspended particles may occur in parallel and kinetic equation has the form where β att and β det are empirical attachment and detachment coefficients.
In the literature, we can find various types of kinetic equations, because there are many factors affecting the process of suspended particles transportation and deposition in porous medium such as particle size distribution [24,[27][28][29], grain size distribution [28][29][30], hydrodynamic forces [30,31], temperature [32,33] and others. In [7,19], a multistage deposition kinetics is used, with ripening, transient and breakthrough stages. Studies of the flow in deep bed filters revealed the existence of two components: the jet stream and stagnant zones in the areas of grain contact [34]. The deposition formed in the pores depends on the place of formation. It was noted [34] that on the frontal parts of grains, that is, in the areas washed by the jet component, the deposition has a dynamic character: the attached particles quite often they break down in a stream and again reattach. Deposition in stagnant zones is more passive. Its formation occurs slowly in comparison with deposition on the frontal parts of grains, but reattachment from stagnant zones is not observed either. According to this, in [34][35][36][37], the porous media is assumed to have two types of zones: active and passive. It is assumed that active and passive zones evenly distributed in porous medium and at each point of the porous medium portion of active and passive zones proportional to their capacity. In active zones, the deposition is formed in reversible form, whereas in passive zones it is formed in irreversible form. The kinetic equations are given by the equations, respectively, for each zone, where ρ a is the concentration of deposition in the active zone (m 3 /m 3 ), ρ p is the concentration of deposition in the passive zone (m 3 /m 3 ), c 0 is the inlet suspension concentration, β a is the coefficient characterizing the kinetics in the active zone (1/s), β p0 is the kinetic coefficient (1/s) of passive zone and associated with the effect of compaction (aging) of the deposition, ρ p1 is the value of ρ p at which "aging" begins, ρ a0 and ρ p0 are the capacities of the active and passive zones, respectively, the total filter capacity is ρ 0 = ρ a0 + ρ p0 . Here, at the beginning of deposition formation, kinetic coefficient of passive zone β p = β p0 . Starting with a certain concentration of ρ p1 , the value of β p becomes less than β p0 , and a further decrease in β p goes inversely with the amount of deposition ρ p . Finally, when the deposition concentration is close to saturation, β p decreases more rapidly. This section is approximated by a step (similar to the approximation of strongly convex isotherms by a rectangular isotherm in the dynamics of sorption) [34]. So the function ρ p (x, t) is continuous, but its first derivation has first type discontinuity at point ρ p = ρ p0 . Most of real suspensions are multicomponent. The solid particles suspended in the suspension can differ in their properties, that is to say, they can have different geometric, physicochemical and electrokinetic characteristics [34,38,39]. Many recent works takes into account the polydispersity of suspension [6,[40][41][42][43][44][45]. Two-component suspensions consisting of two disperse components with different concentrations and kinetic properties represent the simplest form of a multicomponent suspension. In most cases two component suspension have been studied [41][42][43][44][45][46].
As we mentioned above, several varieties of phenomenological relationships have been obtained to describe the flow of various types of disperse systems, but suspension polydispersity and multistage kinetics have not been taken into account together. In this paper, we formulate a problem of two component suspension filtration in a dual-zone porous medium with multistage deposition kinetics. The problem is numerically solved. Influences of suspension polydispersity, dynamical factors and other model parameter on filtration characteristics are studied.

Mathematical Model of Two-Component Suspension Filtration
Let us consider a semi-infinite homogeneous porous media with initial porosity m 0 , filled with a homogeneous fluid. At the point x = 0 starting from t = 0 a two-component suspension with 0 are concentrations of first and second type of particles in inlet suspension, respectively) is injected into the porous media with the velocity v (t) = v = const.
The mathematical model of filtration of two-component suspensions dual-zone porous media is a generalization of the corresponding one-component model [34,[42][43][44]. The system of equations consists of balance equations for each component, Darcy's law and kinetic equations for each zone of media.
The balance equation with a given flow velocity regime, including the diffusion term is where m is current porosity of media, c (i) is the concentration of the ith component of the suspension p are respectively the concentration of deposition in the active and passive zones of the ith component of the suspension (m 3 /m 3 ), D (i) is diffusion coefficient for each component (m 2 /s), i = 1, 2 corresponds to the numbers of the components. These mass balance equations derived on the base of assumption "accumulation rate = (flow in) − (flow out)" [16] for each component. The first term on left hand side describes changes in concentration of suspended particles, the second is convective term, and the third and fourth terms describes changes in concentration of deposition in active and passive zone, respectively. The expression on right hand side describes diffusion term.
For the changes of porosity and permeability of porous media, we use Darcy's law and the Carman-Cozeny equation where K (m) is filtration coefficient, |∇p| is the modulus of pressure gradient, k 0 = const. We generalized the kinetic equation of the process of attachment and detachment of particles in the active zone for a two-component suspension reflecting the dynamic factors [18,37] as follows a0 is the capacity of the active zones for the ith component of the suspension, β (i) a is the coefficient characterizing the kinetics in the active zone of the ith component of the suspension, γ (i) , ω (i) are constant coefficients characterizing intensity of influence of pressure gradient on particle attachment and detachment process in active zone.
The equations of kinetics of the formation of irreversible deposition in passive zone are as follows where β (i) p is the coefficient characterizing the kinetics in the passive zone of the ith component of the suspension.
We denote the partial capacity of the passive zones by ρ (i) p0 . As in [34], we can assume that the total filter capacity is mostly determined by the first component so that ρ (1) p0 > ρ (2) p0 . Consequently, deposition process in passive zone ceases when the total deposition reaches full capacity of passive zone, that is ρ (1) p + ρ (2) p = ρ p0 . In addition, the accumulation of deposition of the second component ceases when the partial capacity is reached, ρ (2) p = ρ (2) p0 , even if the full capacity is not exhausted, i.e., ρ (1) p + ρ (2) p < ρ p0 . In addition, the "kinetic quality" of the first component is higher than second one, i.e., β characterizing the effect of compaction (aging) of the deposition are described in accordance with [34], by the equations where ρ (i) p1 are the limiting concentrations that is from ρ p1 begins compaction (aging) of the reversible form of the deposition for the corresponding component of the suspension.
The initial and boundary conditions are in the form
In the area D = {0 ≤ x < ∞, 0 ≤ t ≤ T}, we consider the net ω hτ = x k , t j , x k = kh, k = 0, 1, ..., t j = jτ, j = 0, 1, ..., J, τ = T/J , where, J is the number of nodes of time. Instead of functions c (i) (t, x), ρ  The balance Equation (7) are approximated by the equations Combining the Equations (8) with (9), (10), |∇p| can be written in grid form as follows The difference scheme for the Equation (11) is For Equations (12)- (14), we have Conditions (15) are also reduced to a grid form as follows where K is a sufficiently large number for which c (i) j K = 0 approximately. After simple transformations, the schemes (18) and (19) take the form We transform the scheme (16) to obtain the following system of linear equations where Using the values ρ where

From boundary conditions we have
Numerical algorithm is follows. First by using the initial conditions we find values c (i)

Results and Discussion
Numerical experiment carried out with program developed by the authors in Python.
The following values of parameters: c p0 . Due to this fact, the capacity of the passive zone for the second component is saturated not only with concentration ρ (2) p , but also with ρ (1) p . Approximately at t ≈ 900 s, the deposition concentration of the second component increases to some x and then decreases. Over time, the interval of increasing of ρ (2) p expands.
With increasing values of the parameters ρ (1) p1 and ρ (2) p1 (Figure 2) near the point x = 0, the capacity of the passive zones is completely saturated with the deposition more quickly, at t ≈ 450 s. Here, unlike the above capacity of the second component is fully saturated with itself an exception is only a few points near x = 0. Furthermore, quick filling the capacity of passive zone leads to increase of c (i) , ρ (i) a . In Figure 3 some results are presented to estimate the effect of diffusion term. A comparison of the graphs with corresponding graphs above shows that diffusion leads to smearing of profiles. Concentration of suspended particles will flow further in this case. This also leads to formation of more deposition at the points away from inlet point in both active and passive zones. Only in passive zone deposition reaches the maximum value and stationary situation is established in this zone at points near to x = 0. But in profiles of c (i) and ρ (i) a , the stationary situation has not established yet in studied time.
From Figure 4 it is seen that dynamic factors affect mainly the formation of deposition in active zone. As in the case of a single-component suspension [37], the potential capacity of the active zone is not completely saturated with deposition. As the parameters characterizing intensity of gradient of pressure for the first component of suspension γ (1) , ω (1) are greater than the parameters of second one γ (2) , ω (2) capacity of active zone for first component is less filled with deposition than that for the second one. This also leads to a little increase in c (i) , ρ (i) p . Figure 5 presents the dynamics of ρ  a and ρ (2) a at the fixed points x = 0 m and x = 0.02 m with and without dynamical factors to estimate the influence of parameters γ (i) , ω (i) . We can see some non-monotonic changes in ρ (1) a . From Figure 6a we can conclude that when dynamic factors are not taken into account over time, the capacity of active zone for each component ρ (i) a0 is completely saturated with deposition at the point x = 0. In the case when dynamical factors are taken into account at the point x = 0 over time first ρ (i) a increases, then decreases, and finally the dynamical equilibrium occurs before capacity of active zone is completely saturated with deposition. Due to increasing of depositions this leads to decrease of porosity, and by increasing the pressure gradient intensity the detachment of deposited particles also increases. From Figure 6b at the point x = 0.02 dynamical equilibrium is not occurred in studied time. Including dynamic factors in kinetic equation leads to increase of concentration of deposition in active zone. This is more noticeable in dynamics of ρ      (1) = ω (2) = 0, D (1) = 2.0 · 10 −6 m 2 /s, D (2) = 1.5 · 10 −6 m 2 /s.   (2) = ω (2) = 0.05, D (1) = 2.0 · 10 −6 m 2 /s, D (2) = 1.5 · 10 −6 m 2 /s.

Conclusions
A mathematical model of the process of two-component suspension filtration in homogeneous porous media with the formation of deposition in passive and active zones has been developed. The system of filtration equations consists of the mass balance equation, the kinetics taking into account the dynamic factors in the active zones, multistage deposition formation in the passive zone for each component of the suspension and Darcy's law. To solve the problem, a numerical algorithm has been developed. Numerical experiments carried out using a computer code developed by authors in Python. The concentration fields and the concentration of deposition are determined for both passive and active zones for each component of suspension. The multistage dynamics of deposition, which characterizes the various features of each stage and their parameters, is considered. The two components of the suspension have different kinetics of attachment and detachment from pores (for active zones). It has been shown that polydisversity of suspension and multistage deposition kinetics can lead to various effects that are not characteristic of the transfer of one-component suspensions with one-stage particle deposition kinetics. The dynamics of particle deposition from both components of the suspension at various points in the medium has also been analysed. It has been shown that at the points of the medium close to the inlet, the concentration of deposited particles in the passive zone can reach partial capacities. It has been determined that near the free surface of the medium, over time, the concentration of the first component of the suspension becomes greater than the capacity of the passive zone for the first component. Due to this the capacity of the passive zone for the second component is saturated not only with this component, but also with the first component. The increase of parameters estimating "aging", leads to a more intense saturation of the passive zone of deposition in the initial stages of the process. It has been established that for some set of initial data used in the calculations, the concentration of deposition of the second component can relatively quickly reach the partial capacities. For the concentration of deposition in active zone, non-monotonic dynamics have been obtained at individual points in the medium. This is explained by the multi-stage kinetics of particle deposition and effect of dynamical factors.