A new analytical formulation of retention effects on particle diffusion processes

The ultimate purpose of this paper is to present a new analytical formulation to simulate diffusion with retention in a reactive medium under stable thermodynamic conditions. The analysis of diffusion with retention in a continuum medium is developed after the solution of an equivalent problem using a discrete approach. The new law may be interpreted as the reduction of all diffusion processes with retention to a unifying phenomenon that can adequately simulate the retention effect namely a circulatory motion. It is remarkable that the governing equation requires a fourth order differential term as suggested by the discrete approach. The relative fraction of diffusion particles β is introduced as a control parameter in the diffusion-retention law as suggested by the discrete approach. This control parameter is essential to avoid retention isolated from the diffusion process. Two matrices referring to material properties are introduced and related to the real phenomenon through the circulation hypothesis. The governing equation may be highly non-linear even if the material properties are constant, but the retention effect is a function of the concentration level, that is, β is a function of the concentration.


INTRODUCTION
Spreading of particles or microorganisms immersed in a given medium or deployed on a given substrata is frequently modeled as a diffusion process. The mathematical formulation of the classical diffusion problem in a three-dimensional space x = [x 1 , x 2 , x 3 ] T can be written as follows: Consider the process depicted in Figure 1. The rule governing the redistribution of contents of each cell indicates that a fraction of the contents αp n is retained in the n th cell and the exceeding volume is evenly transferred to the neighboring cells, that is, 0.5βp n to the left, (n − 1) th cell, and 0.5βp n to the right, (n + 1) th cell, at each time step, where β = (1 − α). This means that the dispersion runs slower than for the classical diffusion problem. Note that if β = 1, the problem is reduced to the classical Gaussian distribution.
Translating this rule into algebraic expressions we get: /2 /2 /2 /2 /2 /2 /2 /2 n n+1 n+2 n+3 n 1 n 2 n 3 n n+1 n+2 n+3 n 1 n 2 n 3 n n+1 n+2 n+3 n 1 n 2 n 3 It is now possible to manipulate Eq. (1a) and Eq. (1b) in order to obtain finite difference terms. It is important to notice that the rearrangement of terms must be done carefully otherwise we could obtain an equation that doesn't reproduce the required solutions for the limits β = 0 and β = 1. That is, the differential terms involved in the final expression may not reproduce the behavior required by the fundamental set of equations for the entire range of the parameters. Therefore, it is necessary to test the intermediate expressions at critical steps to verify if the underlying assumptions related to Eq. (1a) and Eq. (1b) are preserved. This means that, for β = 1, corresponding to no retention, the classical Gaussian distribution is obtained, that is, the contents of each cell is equally redistributed to the neighboring ones at each time step. If β = 0 on the other hand the solution is stationary that is, the total content of each cell remains the same all time.
Subtracting Eq. (1a) from Eq. (1b), after reducing all terms to the time t − 1, leads to: p t+1 n − p t n = β − 1 2 p t−1 n−1 + p t−1 n+1 + (1 − β) p t−1 n−1 + p t−1 n+1 − p t−1 n + 1 4 β p t−1 n−2 + 2 p t−1 n + p t−1 n+2 . (2) After a relatively lengthy and careful algebraic manipulation -see appendix -the following expression is finally obtained: Note that taking β = 1 Eq. (3) reduces to the classical diffusion problem, that is no retention, and with β = 0 Eq. (3) represents a stationary solution since the right hand side term vanishes. Consequently the time rate variation of the contents p n equals zero for all t for all cells.
Calling K 2 = L 2 0 /2T 0 and K 4 = L 4 1 /4T 0 , and taking the limit as x → 0 and t → 0, we obtain: The fourth order term with negative sign introduces the retention effect. Clearly this term comes in naturally, without any artificial assumption, as an immediate consequence of the temporary retention imposed by the redistribution law. We would like to call the attention that similar equations can be obtained by introducing non-linear correlations in the Fick's law (Cohen and Murray 1981). These equations, however, are derived to take into account physical phenomena that hardly could be associated to the retention effect. The result presented above justifies the following statement: PROPOSITION. Suppose a one-dimensional diffusion process consisting of N particles moving in a homogeneous, isotropic supporting medium according to the Fick's law. If a fraction f (N ) = (1 − β) of the diffusing particles is temporarily trapped along their trajectories due to some mechanical, physical, chemical or physicochemical interaction with the medium, the governing equation must include a fourth order differential term of the form: is the fraction of the particles temporarily trapped, with 0 < β ≤ 1 constant, and p(x, t) stands for the particle concentration in the medium.
It will be shown in the sequel that this proposition may be generalized in order to allow β and to be time and space dependent parameters. In this case, however, a more complex differential equation is obtained including non-linear terms, and it would be necessary a more detailed discussion that we are going to postpone for a future paper.

PHENOMENOLOGICAL CONSIDERATIONS OF DIFFUSION PROCESSES
Before stating the laws that govern retention, it is worthwhile describing some basic ideas sustaining the formulation that will be introduced later on. Scientists and particularly engineers have frequently to deal with a class of problems that, behind an apparent simplicity, hides complex events. A problem with hidden complexity results from the simultaneous convergence of several phenomena involving a relatively large number of transformation processes. Most of the problems belonging to this class that appear in engineering practice require the determination of certain variables, that we will call here macrostate variables, representing complex physical, chemical or physicochemical interactions at micro scales. For the engineering point of view, within certain limits, it is not necessary to analyze the phenomenon at the smallest possible scale, because frequently the basic phenomena are not yet fully understood. It suffices knowing the behavior of the macro-state variables. That is, for a problem with hidden complexity, in several circumstances, the main task is to determine the relationship between the input variables and the output variables flowing in and out of a "black box".
The determination of a method or a theory adequate to establish the behavior of the macro-variables is by no means an easy task. The main difficulty arises from the fact that the theory must play a unification role. That is, the theory must as far as possible translate, with a relatively simple set of rules and laws, the response of the macro-state variables to the perturbations introduced in the system under consideration. More than that, a successful theory should give adequate response to several similar phenomena belonging to the same class. For instance, dispersion of gas, liquid or solid particles in different supporting media should be accurately modeled by the same theoretical framework because they belong to the same class of problems independently of the particle nature and of the medium. The macro-variable for this class of phenomenon is the particle concentration denoted by p(x, t). The theory proposed in this paper is derived from the following ideal process.
Consider a collection of particles immersed in some medium and moving in a preferred direction driven by repulsion forces such that they displace from high concentration regions towards low concentration regions free from constraints. The motion is not subjected to any critical barrier that could induce a significant perturbation in the process. Under this circumstance it is plausible to admit that the average speed of the particles is proportional to the concentration gradient, the particles being driven from regions with high concentration levels, strong repulsive forces, to regions with low concentration levels, weak repulsive forces. This is the fundamental idea supporting Fick's law. For the one-dimensional problem in a homogeneous medium, the absolute value of the diffusion speed may be defined as Dgrad( p(x, t)), where D is the diffusion coefficient.
The diffusion coefficient cannot be considered as a parameter characterizing all physicochemical phenomena associated to the diffusion processes. It is instead a coefficient associate to an ideal event that simulates satisfactorily a wide range of phenomena. The diffusion coefficient may have different physical interpretations depending on the particular phenomenon under consideration. Another way to see this methodology is to consider a "black box" that transforms input variables into output variables according to well-established laws, irrespectively from the physics that may be going on inside the box. In this sense, the Fick's law is a decomplexification process. Now suppose that other events interfere in the diffusion mechanism. In this case the law must be modified accordingly. For instance, the effects of sources or sinks in a diffusion process cannot be modeled by some artificial modification of the flow velocity or distorting the interpretation of the diffusion coefficient. A new term is necessary to be incorporated in the governing equation. Diffusion with sinks or sources belongs to a different class of phenomena.
To avoid possible misunderstandings referring to more complex diffusion problems, we would like to remark that it is possible within a same class of phenomena to improve the respective laws incorporating extra terms. For the case of diffusion, for instance, Fick's law may be extended to incorporate higher order terms in the definition of the flux b(x, t). For the one-dimensional case Cohen and Murray, starting from a generalized law for the internal energy, arrived to the following expression for the flux: This means that the driving forces inducing the particle motion are not related only to the gradient of the concentration grad(Ap(x, t)) as in Fick's law, but require the consideration of higher order terms. Since the corresponding governing equation includes for this case two new parameters B and k, it would be presumably possible to obtain a better representation of experimental results by adjusting properly the values of these two parameters. This is a refinement of the theoretical framework, but the basic event remains the same, belonging to the class of diffusion problems without taking into account any other phenomenon as retention for instance. If however there is some essential modification of the phenomenology governing a given process, it is necessary to review the theoretical framework to find some universal model adequate to decomplexify the corresponding class of phenomenon. Unfortunately attempts to extend a theory fitted to represent a given class of events to a new phenomenological class are not seldom found in the literature. The results are certainly not universal and may fail to express a reliable response of the macro-state variables. This is the case of diffusion with temporary retention. The phenomenon of temporary retention cannot be adequately modeled by Fick's law or any refinement of the classical diffusion law. Diffusion with temporary retention belongs to another class of phenomena. A new suitable law could be stated axiomatically. It would be better to find a kind of universal event that may encompass a large number of phenomena and may also provide a consistent evaluation of the macro-state variables.
The main purpose of this paper is to propose a model that can be used to simulate with good accuracy the retention effect in a diffusion process. Based on that model, a new universal theoretical framework to decomplexify the temporary process is derived.

THE DIFFUSION CONSTITUTIVE LAW
The diffusion law adopted here is the classical Fick's law. Define D as: where n is the unit vector normal to ∂ V and d is the diffusivity vector. The diffusivity vector consistently with the Fick's law is defined as: Let us call D the diffusion matrix that is related to the physical properties of the medium, and p(x, t) is the mass concentration function. The components of D characterize the diffusion properties of the mass concentration in the embedding medium. For the general case, the matrix D is full (MacElroy 1996) and the diffusivity vector reads: The minus sign comes from the well-known fact that under normal conditions the motion direction is oriented towards the regions of decreasing mass concentration.

THE RETENTION CONSTITUTIVE LAW
In this section we will introduce the key contribution of this paper, which is the derivation of the function representing the retention rate of the particle concentration trapped on the boundary layer in the vicinity of ∂ V . The reasoning that will follow was guided by the result obtained with the discrete approach in the sense that the final equation governing diffusion with retention should include fourth order differential terms.
Instead of using an axiomatic approach through the definition of analytical expressions that would lead to the expected results, we will explore a plausible physical interpretation along with the development of the mathematical formulation. The physics of the problem postulated here is not unique, but is adequate to simulate a large range of cases. A consistent interpretation can also be given to the material parameters that will be introduced later on. The reader must have in mind that our purpose is to achieve a general law, that is, a law that would fit possibly all cases of diffusion including temporary retention. This means that the basic phenomenon introduced here to explain retention, has to be seen as a universal phenomenon or, more precisely, a universal virtual phenomenon.
The study of particular cases trying to simulate as close as possible the physical and chemical interactions governing the retention effect belongs to other approaches focusing very specific problems. Each problem requires a particular solution that cannot be extrapolated to other cases.
Two main ideas are enclosed in the definition of R : i) Recall that the terms representing the retention effect according to the previous derivation must be multiplied by β(1 − β) as shown in Eq. (4). To be consistent with the discrete approach, the effective mass temporarily trapped per unit of time R must be proportional to the relative mass of the diffusing particles β as will become explicit in the next section from Eq. (5a). Clearly the factor β(1 − β) is a control parameter imposing variation limits on the blocking effect to avoid unacceptable mass accumulation on the boundary layer next to ∂ V , where otherwise the density might grow beyond physically acceptable limits. Moreover it is also physically consistent to admit the retention effect directly proportional to the diffusing fraction. Indeed, the probability of particles being trapped increases with increasing intensity of particles diffusing in the medium.
ii) Now let us introduce a vector b for the retention state similar to the vector d for the diffusion state. The effective mass trapped in the boundary layer per unit of time, that is R , is proportional to the component of the vector b orthogonal to the boundary surface ∂ V . Considering the proportionality law proposed in i) above we may write: An Acad Bras Cienc LUIZ BEVILACQUA, AUGUSTO C.N.R. GALEÃO and FLAVIO P. COSTA The vector b will be called blocking vector. We introduce now the fundamental hypothesis to derive the formation law of the blocking vector.
First of all let us recall that we are dealing with temporary retention. Permanent retention corresponds to a sink that belongs to another class of phenomenon. Temporary retention may occur under different forms that are very difficult to represent analytically. As a matter of fact, if specific formulations for each kind of phenomenon are to be derived -temporary interstitial retention, temporary blocking by mechanical obstacles, chemical reactions, tortuosity of the substrata -it would be necessary to develop very complex theories, each one corresponding to a specific class of phenomenon.
What we intend to do here is to detect the key phenomenon that is common to all or at least to the majority of the diffusion processes under the effect of retention. That is, we are seeking for a unified model that represents the fundamental effects of the retention process for a large range of cases. The retention coefficients appearing in the evolution equation will in some cases represent the parameters corresponding to a virtual phenomenon simulating the real one.
The unifying phenomenon that will be assumed to simulate the retention effect for all cases is the circulation in a boundary layer in the neighborhood of the surface ∂ V . That is, a fraction of the particles follows a helicoidal trajectory consisting of n spirals distributed on very thin boundary layer s in the neighborhood of the surface ∂ V (Fig. 2a). While travelling along this helicoidal trajectory, the particle will be considered trapped in the boundary layer in the vicinity of the surface ∂ V . The retention effect is therefore simulated by a temporary circulation in the neighborhood of the surface enclosing the volume V . The time delay for a particle to cross the boundary layer is proportional to the number of loops n and the angular velocity ω. In order to simplify the analysis, consider the loops compressed into a single closed trajectory contained in a plane tangent to the boundary. That is, the particle moves with velocity v on a circumferential path with radius a on the plane . To simulate the displacement of the particle along the direction parallel to n after the completion of n loops covering a distance s in a time t, the plane is assumed to move inwards the volume V with a velocity u. The retention effect is therefore simulated by a circulation with angular velocity ω = |v| a n on the plane coupled with a translation inwards the volume V with velocity u. The modulus of u may be easily calculated |u| = s t = ωλ 2π , and therefore u = ωλ 2π n. The modulus of u is very small since it is assumed that the pitch λ of the helicoidally motion is very small consistent with the hypothesis of retention within a very thin boundary layer. That is, the fraction of the particles moving on a temporarily closed trajectory (Fig. 2b) coupled with a very slow motion u of the circulation plane inwards the volume V will be considered trapped on the surface ∂ V.Thus, the retention phenomenology defined here requires the particles to be temporarily confined in a boundary layer in the vicinity of the surface ∂ V .
In order to describe the equivalent circulation motion, two different matrices of material parameters will be introduced. The first one denoted by C characterizes the velocity, of the equivalent circulatory motion. The second material matrix called R expresses the effectiveness of the blocking effect through the projection of the enclosed area by one loop on the tangent plane (t 1 xt 2 ) of ∂ V multiplied by the number of loops, such that the time delay equivalent to the temporary retention is achieved. Cases where a fraction of the dispersing particles is physically trapped for a certain amount of time on some material obstacle, as is the case of porous media, can be interpreted as a singularity where the radius of the circulation tends to zero and the motion is reduced to a spin. Description of some specific physical phenomena referring to the blocking time can be found in (Deleersnijder et al. 2006, Green 1996, Mota et al. 2004). Now given the above hypothesis let us introduce the retention law.
1. The retention phenomenon occurring simultaneously with diffusion in a reactive medium under a state of stable thermodynamic equilibrium may be simulated by a circulatory motion defined by the rotational of an appropriate vector t given by: where C is a matrix of equivalent physical constants. That is the components c i j of C refer to a circulation effect that adequately simulates the real phenomenon with respect to the dynamical and geometric characteristics of the circulatory motion.
2. The blocking matrix B is given by: and the blocking vector b representing the retention rate is given by: (Cgrad( p(x, t))) .

An Acad Bras Cienc
LUIZ BEVILACQUA, AUGUSTO C.N.R. GALEÃO and FLAVIO P. COSTA The components of the matrix R are also related to material properties that define the retention effectiveness that is the time delay imposed by the blocking effect. The matrix R stands for the control mechanism regulating the circulation intensity that is the number of loops in the equivalent spiral motion. The experimental determination of C and R defines the equivalent circulation associated with the effective actual phenomenon.
3. The net retention rate R corresponding to the particles temporarily trapped on the surface ∂ V is given by the projection of b on the normal to the surface ∂ V multiplied by the particles diffusing fraction β.
That is: Now taking the physical interpretation explained before for the retention effect as a circulation plus a translation, it is obvious that b = u. Note that we are dealing with a process evolving under a thermodynamic stable equilibrium. For more complex cases involving thermodynamic coupling the energy conservation principle must be taken into account. The theory for those cases introducing the retention effect has still to be developed.
The matrices C and R containing the physical parameters c i j and r i j depend on the interaction properties between supporting medium and diffusing substance. The determination of the dimensions of the variables intervening in the retention law may help to find the physical meaning of the components of C and R. Let the physical meaning of t be defined as the angular momentum of the trapped particles per unit of volume. The corresponding dimension is written as: The notations M, L , T mean the dimensions of mass, length and time, respectively. Combining the expression above for the dimension of t with the definition of rot(t) introduced in the retention law we have: From which follows the dimension of the components of the matrix C; Note that this is the same dimension as for the diffusion coefficients composing the diffusion matrix D. Now the dimension of B = grad(rot(t)) comes immediately: Now combining the dimension of b required by the principle of mass conservation with the definition of the blocking vector b introduced in the retention law, b = R div (B), we obtain: T . This is a plausible definition of the dimensions of the variables related to the retention law. From the dimensions of rot(t) and b, the dimensions of the components of C and R can be readily obtained, namely [L 2 /T ] and [L] 2 , respectively, which is a an expansion speed and an area as suggested by the phenomenological assumption introduced before. The dimension of the components of the blocking matrix B shows that B may be interpreted as the time rate of the specific mass variation. It would be possible to generate other slightly different combinations of the fundamental dimensions but, if the concentration and concentration gradient are preserved as key elements any other alternative would not introduce any essential deviation of the interpretation presented here.

An Acad Bras Cienc
The physical components c i j and r i j must be determined experimentally. These coefficients depend exclusively on the physicochemical or mechanical properties of the reactive medium interacting with the dispersing particles. One has to keep in mind that the experimental determination of C and R based on the present theory fits to a circulation effect simulating the actual retention phenomenon. The coefficients c i j and r i j may represent the real physical reaction depending on how the circulation hypothesis fits to the real case under consideration. This simplification, or, as stated before, this decomplexification process, is by no means a disadvantage of the present theory as compared with other approaches available in the specialized literature. The other theories are also based on a set of hypotheses that rarely explains the connection with the actual phenomenon. So the experimental coefficients obtained using other theories equally fails to represent those associated to the real case. The present formulation has the advantage of making it explicit the kind of simulation adopted namely the circulatory motion. Also the same equivalent phenomenon is taken for all cases that allows for the comparison of the coefficients c i j and r i j associated to different cases since all of them are modeled through the same universal representation.
It was mentioned before that the dimensions of the matrices D and C are the same. Also both appear in the formulations of the diffusion law and the retention law as well, as factors of the same quantity composing similar expressions, Dgrad( p(x, t)) and Cgrad( p(x, t)). This suggests that it would be plausible, in the absence of other more significant reason, to assume that they are equal, that is C = D. If this is true, the retention coefficients determining the temporary blocking process are reduced to those composing the matrix R. This would simplify considerably the experimental investigation. For sake of generality we will keep both independent in the sequel since there is no significant simplification from the theoretical point of view to make them equal.
Maybe most important of all is the introduction of higher order differential terms resulting from the present formulation of the retention law. This law is clearly inspired by the results obtained with the discrete approach without which it would be hard to establish a consistent representation of temporary retention. The presence of fourth order terms, consistent with the discrete approach, naturally follows from the law proposed above. Now for the sake of completeness let us write the expanded expressions for the rot(t), B and b according to the above functional representation. In expanded notation we have successively:

An Acad Bras Cienc
LUIZ BEVILACQUA, AUGUSTO C.N.R. GALEÃO and FLAVIO P. COSTA Up to now we have left the theory as general as possible within the framework established for the representation of diffusion with retention. In this sense, the matrices D, C and R have been assumed to be complete. It is quite reasonable, however, for several cases, to assume that D, C and R are diagonals matrices, that is d i j = 0, c i j = 0 and r i j = 0 for i = j. This hypothesis holds for the case of orthotropic processes. For the more general case of skewed anisotropy the matrices D, C and R cannot be reduced to the diagonal terms.
For the sake of simplicity we will restrict ourselves in the sequel to orthotropic diffusion processes with retention. This will not impose any critical limitations to the theory developed here.
The diffusivity vector can now be written: Similarly for the variables associated to the retention process we get: For this case, the matrix B coincides with the Jacobian of rot(t), that is: The vector b reads:

THE CONTINUUM APPROACH FORMULATION OF DIFFUSION WITH RETENTION
Now once established the basic mechanisms for the diffusion with retention, let us derive the governing equation assuming, as said before, that the process takes place under stable thermodynamic conditions. Consider a sufficient large set of elementary particles scattered in some reactive medium, that is, a medium with the property of delaying the motion of a certain number of particles for a given period of time as reported by some authors in physicochemical experiments (Atsumi 2002, D'Angelo et al. 2003, Deleersnijder et al. 2006, Green 1996, Liu and Thompson 2002, Muhammad 2004. Consider the open set V bounded by the surface ∂ V at a given time t (Fig. 3). The total number of particles in V ∪ ∂ V is divided into two groups. The first group comprises the particles belonging to V , and the set of particles distributed in the

An Acad Bras Cienc
A NEW FORMULATION FOR DIFFUSION-RETENTION PROBLEMS neighborhood of ∂ V constitutes the second group. The specific mass of the dispersing particles, considered as a continuum, at that time t referred to an arbitrary point x = x 1 x 2 x 3 T is defined by p(x, t). For the sake of simplicity we will use the Cartesian system for the analytical development. Call M 1 the net mass variation inside the volume V in the time interval t. That is: Now, suppose that the medium interacts with the elementary particles in such a way that a fraction of the particles on ∂ V scatters through the medium, and the complementary fraction is trapped in a state that we will call temporarily retention state. Particles belonging to this state are captured back into the volume V crossing the boundary at very low speed and remaining temporarily within a thin boundary layer next to ∂ V as explained in the previous section. In this sense we say that, while in the retention state, the particles are trapped. The time rate of mass variation per area unit given by D characterizes the diffusion process, that is the outward particle flux through ∂ V , and the time rate of mass variation per area unit given by R characterizes the retention effect, that is the inward particle flux through ∂ V . Both processes are active simultaneously on ∂ V . Let us assume that diffusion and retention are independent processes homogeneously distributed on the surface ∂ V covering two disconnected partitions of ∂ V , that is, β∂ V for diffusion and (1 − β)∂ V for retention. This means that the particles on ∂ V are either spreading -diffusion -or trapped -retention. Therefore, calling M 2 the mass variation in the volume V due to the net flux through the boundary ∂ V in the time interval t, we have: An Acad Bras Cienc LUIZ BEVILACQUA, AUGUSTO C.N.R. GALEÃO and FLAVIO P. COSTA The principle of mass conservation in the absence of sources or sinks reads ( M 1 + M 2 ) t = 0 or recaling Eq. (5a) and Eq. (5b): As stated above, the variables D and R are the effective diffusion rate and the effective retention rate, respectively. Using the constitutive relations for diffusion and retention introduced in the previous sections, Eq. (5c) may be rewritten: Now assuming the function p(x, t) sufficiently smooth and using the divergence theorem we get: Introducing the expressions of b and d previously derived and taking the limit as t → 0 we obtain: Since V is arbitrary and expanding the representation of B: (C grad( p(x, t)))) + div βdiv (D p(x, t)) .
Or using the component notation: Equation (7b) is the governing equation for the case of diffusion with retention. It differs significantly from some others approaches found in the technical literature to simulate this type of phenomenon. The theory developed above is consistent with the straightforward discrete derivation of the population dispersion with temporarily retention to colonize the occupied territory. Equation (7a) is in fact the generalized version of Eq. (4) corresponding to the discrete derivation. The determination of the physical components R, C and D depends on experimental efforts. Anyway, starting from the Eq. (7a) it should be possible to devise meaningful experimental strategies since it is clear what should be detected and measured. For one-dimensional diffusion with retention along the x axis, the Eq. (7a) reduces to a simpler expression, namely: Despite that the above equation is considerably simpler, it keeps some very fundamental information that applies in general for more complex mass diffusion problems. Particularly interesting is the role of the An Acad Bras Cienc parameter β. Note first that the diffusion term grows steadily with β, while the retention term varies according to the logistic equation (Fig. 4). As a consequence, low intensity diffusion processes do not mean high intensity retention. The control of the retention process introduced in the governing equation avoids the existence of pure retention as expected, but allows for diffusion only, as it is also consistent with the physics of the process. Now it is also reasonable to expect that the retention mechanism saturates after some time. This idea may be introduced by defining the parameter β as a function of p(x, t). Let us assume β, a power function of the saturation quotient q(x, t) = p(x, t) p * , that is β = [q(x, t)] m , where p * is the saturation limit of the mass concentration. With this expression for β and with r , c and k constants, the Eq. (7c) gives: where the coefficients φ i , i = 1, 2, 3, 4 are given by: Figure 5 illustrates the variation of β with q(x, t) for several values of m. For m < 1 the retention effect reaches the maximum, that is β = 0.5 for very low saturation quotient. That is, for relatively low levels of mass concentration, the capacity of the supporting medium to block diffusion is close to the upper limit. For m > 1 the situation is reversed. The maximum blocking effect is reached for relatively high saturation quotient (Fig. 5). For the case of linear correlation between β and q, that is m = 1, the Eq. (7d) reduces to:

An Acad Bras Cienc
Close to the maximum retention effect that is when q(x, t) is in the neighborhood of 0.5, we may write q(x, t) ≈ 0.5 + ε(x, t). Then, close to this point, the Eq. (7e) can be written:

CONCLUDING REMARKS
The phenomenon of diffusion with retention has been already observed and reported by several authors. Analytical formulations have also been proposed by some authors based on certain assumptions focused on the particular problem under consideration. All of them assume consistently the Fick's law as the analytical representation of the diffusion process even in the presence of retention effects. Consequently, the governing equations for diffusion with retention appearing in the literature so far, to the best of our knowledge, belong to the class of the second order partial differential equation as in the classical diffusion approach. They are as a matter of fact modified versions of the classical diffusion equation. As a consequence it would be virtually impossible to derive a general theory since each problem requires a particular approach. The classical solution could apply with the introduction of sources and sinks, such that the overall effect on the number of particles present in the medium remains unaltered after sufficiently long time. But as explained before, this is a different phenomenon since it requires energy exchange with the surrounding medium.
We believe that the diffusion through the brain intracellular space (Nicholson andSykova 1998, Mota et al. 2004) may be analyzed with the theory proposed here. The intracellular spaces, that is, the substratum where neuroactive substances and molecules circulate, introduce a three-dimensional diffusion process due to the tortuosity of the free trajectories. If the tortuosity acts as an obstacle introducing only a kind of circulation orthogonal to the main diffusion direction, the process may be simulated by a temporary delay and the approach introduced here may represent this type of phenomenon. This is a typical example of a complex problem being simulated by a general decomplexfication model. This paper shows that the second order differential equation is not adequate to represent the physics of diffusion with retention. The retention process requires a fourth order term multiplied by a control factor β(1 − β) that is essential to avoid the formation of unattainable mass concentration. It is remarkable that the solution of the discrete problem of population invasion with temporary territory occupation provided the fundamental clue to derive the solution of the respective continuous problem. It has been shown very clearly that a fourth order term appears in a straightforward manner to represent the retention phenomenon. We believe that this argument suffices to prove the need to introduce this term to achieve an accurate model for diffusion with retention. Problems with a source term are even more sensitive to the presence of the fourth order term. Indeed the stability conditions cannot be correctly formulated without the consideration of the retention coefficients r i j and c i j . Moreover, the redistribution parameter β plays a very important role on the stability conditions. Indeed, for β ≈ (1 − ε) where ε is small, solutions of the type ζ = ζ 0 exp 2πi x λ − ωt may cause instability for a very large range of frequencies independently of the values of the diffusion and retention coefficients (Bevilacqua et al. in press). This means that a state of relative intense diffusion coupled with an external source is very sensitive to small retention perturbations for a very large frequency range.
Fourth order terms may also appear when the Fick's law is expanded to introduce nonlinear terms (Cohen and Murray 1981). This approach, however, has no connection with the retention phenomenon. Moreover, the nonlinear expansion of the Fick's law introduces always nonlinear terms in the governing equation and ignores the control parameter multiplying the fourth order term. The approach presented here doesn't require necessarily the introduction of nonlinear terms in the governing equation.
It is worthwhile mentioning the occurrence of fourth order terms in other engineering and physical problems as those related to the determination of thin films profiles. Examples are molecular beam epitaxy (Barabási and Stanley 1995) and surface coating with viscous fluids (Schwartz and Roy 2004, Myers 1998, Mullins 1957. Fourth order term also appears in the partial differential equation for the determination of the entropy evolution connected to the Shannon information density (Broadbridge 2008).
The theory presented here is satisfactory for moderate temperatures and small temperature variations. The extended theory including thermodynamic coupling has still to be developed. Classical approaches can be found in the current literature (Kirkland et al. 1990, Eon et al. 1997, Fitch et al. 1996, Schimpf and Giddings 1989, but there is no satisfactory theory including retention effects. The determination of the functional dependence of β on mass concentration, on temperature and on other intervening molecular forces is, for instance, a challenge for further investigation. We believe, however, that, before extending the theory to include energy equations, the experimental validation of the blocking effect as introduced here, for thermodynamically stable processes, should be done. To concentrate efforts to find experimental evidence of the functional dependence of β on the saturation quotient is of particular importance. The terms of the matrices C and R have also to be determined experimentally. As a first step it is advisable to focus the attention on homogeneous media and one-dimensional flow. For these cases, the terms of matrices R and C come together (rc) which makes the experimental analysis easier despite being limited regarding the phenomenology inherent in the process.
Analytical solution of the equations derived above could be found maybe only for very particular cases. We believe that efforts should be concentrated on qualitative behavior, asymptotic and numerical solutions. Similar equations have been studied by several authors (You and Kavet 2000, Lu et al. 2007, Rubinstein et al. 1989, Bernoff and Bertozzi 1995, Myers 1998.
We believe that the theory presented here is consistent, compatible with a very robust discrete approach and represents satisfactorily the retention phenomenon in reactive media. We have introduced a general retention matrix for the sake of completeness. Possibly the retention matrix is also a diagonal matrix, which simplifies considerable the formulation.
Finally other type of phenomena in reactive media as counter-diffusion may occur (Brandani et al. 2000). For this case a reactive medium is one that has the property of triggering a mechanism that partially reverses the flow direction. This means that particle dispersion may occur simultaneously as diffusion processes in opposite directions. This mechanism may be thought of as being induced by the presence of ionic channels distributed in the support medium that are able to partially reverse the flow direction. The governing equations for this case may be derived in a similar way. Fourth order terms, however, may not appear in the governing equations.

ACKNOWLEDGMENTS
This research was partially supported by the Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq) Research Fellowship granted to the first two authors, by the CNPq Research Grant "Universal 14/2009" and by Fundação de Amparo à Pesquisa do Estado da Bahia (FAPESB).