# Derivation and mathematical study of a sorption-coagulation equation

###### Abstract.

This work is devoted to the derivation and the matematical study of a new model for water-soluble polymers and metal ions interactions, which are used in chemistry for their wide range of applications. First, we motivate and derive a model that describes the evolution of the configurational distribution of polymers. One of the novelty resides in the configuration variables which consider both, the size of the polymers and the quantity of metal ions they captured through sorption. The model consists in a non-linear transport equation with a quadratic source term, the coagulation. Then, we prove the existence of solutions for all time to the problem thanks to classical fixed point theory. Next, we reformulate the coagulation operator under a conservative form which allows to write a finite volume scheme. The sequence of approximated solutions is proved to be convergent (toward a solution to the problem) thanks to a stability principle. Finally, we illustrate the behaviour of the solutions using this numerical scheme and we intend to discuss on the long-time behaviour.

###### Key words and phrases:

Polymers; Metal ions; Coagulation equation; Existence of solutions; Finite volume scheme; Weak stability###### 2010 Mathematics Subject Classification:

Primary: 65R20, 82C05; Secondary: 35Q82, 35Q92###### Contents

- 1 Introduction
- 2 Rescaling and existence of global solutions
- 3 Numerical approximation
- 4 Numerical illustration and long-time behaviour
- 5 Conclusion

## 1. Introduction

### 1.1. Motivations

A class of macromolecules, the polymers, has emerged for their potential applications in various fields, as superconducting materials, liquid crystal, and biocompatible polymers. Also, in environmental science, they can be used for instance to remove pollutant from aqueous solutions, or bacteria, fungi and algae [Rivas2011, Rivas2003]. We are particularly interested in polymers which perform interactions with metal ions such as copper ions (Cu), lead ions (Pb), and many others.

One of the applications of this concept, being very promising, lies in membrane separation process. This technique is used to obtain highly purified water in industrial process, or to the contrary, to “wash” water after an industrial process before release it in the environment [Zagorodni2007, Rivas2003]. The idea is to take advantage the interactions of particular polymers with metal ions in order to retain free metal ions from an aqueous solution (water). The solution is then filtrate through a membrane, and due to the high molecular weight of the polymers, they cannot cross the membrane. It results that, on one side, ions are retained by the polymers, while on the other side, water has crossed the membrane and is free of ions. We refer to [Rivas2006] and the reviews [Rivas2003, Rivas2011] for a more precise description of this process.

Nevertheless, the development of this technique leads to some technical difficulties and gives rise to important questions in order to produce efficient methods. Among these, they are the role of fouling effect (adhesion of polymers to the membrane), aggregation phenomenon, concentration effect, interaction with the wall of the cell (recipient) and interaction with the fluid. For recent findings, techniques and models on the subject we refer to the works in [Palencia2010, Palencia2011, Palencia2009, Palencia2009b, Moreno2002, Moreno2006, Rivas2004].

To go further in this direction of a better understanding of such processes, we decided to propose a new model that accounts for the polymers - metal ions and polymers - polymers interactions. The model we develop here differs from the one derived for instance in [Moreno2006, Moreno2002, Rivas2006] where they used macroscopic quantities and equilibrium theory. Indeed, here we propose a new approach, ever used in other physical problems such as Ostwald ripening or coagulation-fragmentation equation, but, up to our knowledge, has never been formalised in this way for water-polymers in interactions with metal ions. The particularity of this problem reveals mathematical difficulties as for the existence and uniqueness of solution, the numerical approximation, and the long-time behaviour of the solution. In this work, we start to answer to some of these questions and give some hints on the behaviour of the solution.

The remainder of this introduction set the theory and derive the model studied in the present work.

### 1.2. Theory and model

A polymer is a macromolecule made up of many repeat units called monomers. To illustrate this definition, an example of polymer is the poly(acrylic acid), a synthetic polymers. Indeed, the acrylic acid which has the molecular formula plays the role of the monomer, whereas the poly(acrylic acid) is a succession of acrylic acid, forming a chain, and has the molecular formula , where stands for the number of occurrences of the monomer.

The model we investigate below concerns highly soluble polymers, particularly water-soluble polymers. Rather, it aims at modelling the interactions between these polymers and metal ions in solution. These interactions are responsible for the retention of metal ions by the polymer and take place in specific sites of the chain, called functional groups. The functional groups are repeat subunits of the polymer able to hold one metal ion each. It is a group of atoms are molecules and it can be, or the monomer itself, or else composed of several monomers. Its composition fluctuates according to: the structure of the chain; the metal ions used; as well as the type of interaction committed. The reader can refer to [Rivas2003, Rivas2006] to get some examples of such polymers, but, among these we mention polyelectrolytes and polychelatogens. The former is a class of polymers whose functional groups are charged, so metal ions bind to them by an electrostatic interaction. For instance, the monomer of the poly( acrylic acid) polymer above-mentioned, loose a proton () in solution leading to a negatively charged polymer. While for the latter metal ions bind to functional groups forming coordinate bonds. However, some polymers combine both interactions, even sometimes others weakest interactions can take place, see [Rivas2011]. Furthermore, all depends on both the polymer and the metal ion considered. Nevertheless, here we follow [Rivas2006] and we assume the interactions modelled might be formulated as reversible chemical reactions. Thus, we acknowledge two processes here: the binding process which consists in the association between a metal ion and a functional group of a polymer also named adsorption; and the opposite process which consists in the dissociation of a metal ion from a functional group or desorption. Both processes formed the sorption phenomenon.

Modeling such interactions, between water-soluble polymers and metal ions, recently received a considerable attention. Several models has been developed, analysed and compared to experimental data, we would mention again the works made in [Moreno2006, Moreno2002, Rivas2006]. The authors mainly deal with macroscopic quantities, which are the concentrations of metal ions bound to polymers and free, the concentration of polymers, etc, together with equilibrium theory to build suitable models.

In the present work, we derive a model that accounts for the evolution of a system of interacting particles consisting of water-soluble polymers and metal ions in solution. Its establishment resides in a chemical formulation of both microscopic processes (association and dissociation). For this, let us label by a polymer consisting of functional groups, such that of them are associated to a metal ion. The variable can reach any numbers (at least theoretically), indeed the polymers can be as big as the technology permits it. Then, the number belongs to by definition of a functional groups, since each functional group reacts with one and only one metal ion, thus a polymer made up of functional groups cannot bind more than metal ions. Remark, it makes clear that if is the number of functional groups occupied by a metal ion, it is also the number of metal ion bound to the polymer. The couple is named the configuration of the polymer. Finally, we label by a free metal ions, free meaning that it is dissolved and not bound to a polymer.

The reversible chemical reaction between one free metal ion and one polymer is

(1.1) |

where is the binding rate (or association rate) at which a polymer of configuration interacts with a metal ions to get a new configuration , while is the dissociation rate at which a metal ions is removed from a polymer , and this latter gets the new configuration . Both rates depend on the configuration of the polymer and inherently on the surrounding conditions (pH of the solution, temperature, etc) supposed to be fixed.

In addition to polymer - metal ion interactions, we take into account polymer -polymer interactions. Actually, under experimental conditions, the presence of other polymers leads to inter-polymer complexes and can produce gels or precipitates, as pointed out in [Rivas2003]. We incorporate in the model the formation of inter-polymer complexes through a binary coagulation, which is the formation of a biggest polymer from two smaller ones. This reaction can be written as follows:

(1.2) |

where the coagulation rate is the rate at which a polymer with the configuration and an other one with configuration will produce a new biggest macromolecule. The reaction preserves the number of functional groups and metal ions bound, so the new polymers gets the configuration .

We devote the next section to establish the set of equations studied in the remainder of this paper. To this end we would finish here by some comments on the mathematical objects used.

First, we use evolution equation which means solutions to the system will be time-dependent functions and are not at the equilibrium yet.

Second, the experiments use generally an high number of particles (both polymers and metal ions). That is why quantities are expressed in term of their molar or mass concentration. So that, we introduce two new convenient variables, and , respectively for the mole of functional groups and the mole of associated (to metal ions) functional groups, where the inverse of the Avogadro’s number. Both variables and can reach a continuum of values, so we consider them as continuous variables on the contrary to variables and which are natural number. Thus, the states of the system at a time will be given by a configurational distribution function where the quantity

provides the molar concentration of polymers in the system having a configuration in the subset of the configurational space (defined latter).

###### Remark 1.1.

It could be possible to keep a discrete version of the problem considering the concentration of all the with and the number of functional groups (as natural number). Nevertheless, it is sometimes mathematically less tractable due to the huge number of equations involved. But, the continuous model can be seen as a limit of the discrete one, through an appropriate rescaling and here we have in mind the parameters . The reader can make himself an idea, on what would be the discrete model here and how it is linked to the continous one, in [Laurencot2002c, Vasseur2002] where rigorous derivations of the continuous Lifshitz-Slyozov equation from the Becker-Döring model is made. Both model being closed to the one presented here.

### 1.3. Equations

Below, we use the two continuous variables for the polymer configuration. First, the quantity in mole of functional groups. Second, the quantity in mole of functional groups associate each to a metal ion. We name, the variable : the quantity of occupied functional groups (in mole). Thus, the set of admissible configurations is:

Indeed, the polymer can reach any size, so its number of functional groups can be as large as we want. While its number of occupied functional groups can belong to , i.e up to the total number of functional groups of the polymer. Then, we define the configurational distribution function of polymers, denoted by , as a function of time over the configuration space . The system governing the evolution of is

(1.3) |

where is the coagulation operator and denotes the rate of association-dissociation or by other means the sorption rate. The latter determines the mechanism of ions transfer with a polymer (association and dissociation of ions), it is the continuous operator describing reaction (1.1), while the former rules the polymer - polymer interaction given by (1.2), where both are defined below. Equation (1.3) on the configurational distribution function is not enough to characterised all the system. To complete the model, it remains to introduce a second equation on the molar concentration of free metal-ions (being in the solution but unbound to polymers). We denote this concentration by , as a function of time , and it is given by the constraint of metal ions conservation (bound and free), namely

(1.4) |

where is a constant which stands for the total quantity of metal ions in the system (bound and unbound to polymers). Indeed, as approaches the molar concentration of polymers with configuration in the limit of and small, if we multiply this quantity by the number of occupied functional groups (= the number of metal ions bounded) in mole, thus we get the molar concentration of metal ions bounded onto the polymers with such a configuration. Thus, the integral term in (1.4) account for the molar concentration of metal ions associated to polymers and the balance equation (1.4) expresses the conservation of matter, between associate and free metal ions in the system.

Let us now define more precisely the rate of metal ions association-dissociation. As a general form for we consider the following chemical reaction rate

where is the association rate, or adsorption, at which a free monomers bind to a polymer (depending on the type of interaction and the diffusion rate of the particles) and is the dissociation rate, or desorption (depending on the strength of the interaction). The association rate is multiply by , as a classical law of mass action and for sake of simplicity we restrict ourselves to which is the order of the reaction. A relevant example of association-dissociation rate would be an analogy to the Langmuir’s law (for the adsorption of metal ions onto a surface), namely

(1.5) |

with parameters and geometrical factor, see [Somorjai2010]. Indeed, the association rate depends on the quantity of available functional groups while the dissociation rate depends only on the quantity of metal ions bound to the polymers (the more metal ions, the more the probability of a dissociation is great).

Next we explicit the coagulation operator with the coagulation rate defined as a nonnegative function over satisfying the symmetry assumption

(1.6) |

It gives the rate at which two polymers with the configuration and will coagulate. The symmetry assumption follows form the impossibility in the system to distinguish the coagulation of a with a or the coagulation of a with a because it is the same reaction. Then, can be decomposed by a gain term and a depletion term that is

where

(1.7) | ||||

(1.8) |

The gain term accounts for the production of a polymer with a configuration thanks to the coagulation of a with and a with . Likewise, the depletion term account for the disappearance in the system of a polymer when coagulate to any other for the benefit of a new polymer with configuration .

### 1.4. Contents of the paper and related works

The remainder of this paper is devoted to the existence of global solutions to problem (1.3)-(1.4), its numerical approximation, simulations and a discussion on the long-time behaviour.

In Section 2, we establish in Theorem 2.2 the existence of global solutions. The proof is based on fixed point theorems two treat the non-linear terms, one for the coagulation operator (1.7)-(1.8), the second for the constraint (1.4). The technique used, implying hypotheses on the regularity of the coefficients, is based on the works on Lifshitz-Slyozov (LS) equation in [Collet2000] and LS with encounters (coagulation) in [Collet1999]. Techniques also adapted in [Helal2013] for biological polymers. The LS equation is a size structured model for clusters (polymers or more general) formation by addition-depletion of monomers [Lifshitz1961], while LS with encounters also take into account merging clusters . The coagulation (only size-dependent) is part of the class of coagulation-fragmentation (CF) equation, where fragmentation is the reverse operator (break-up, splitting of clusters), and it has been studied from a mathematical point of view for instance in [Dubovskii1996, Laurencot2000], also in [Laurencot2002, Amann2005] for the CF equation with space diffusion, and in [Broizat2010] which generalized CF equation with a kinetic approach. The particularity of the model here, is the boundary conditions and the conservation involves. Both necessitates a careful attention in the proof of existence. The former requests to treat the characteristics that come from the boundary, the latter need a particular attention due to the nature of the configurational space.

In Section 3, we propose a finite volume scheme to construct numerical solutions approaching the problem. The numerical scheme is in the spirit of the works made in [Bourgade2008] and [Filbet2008], where the authors propose a reformulation of the CF equation in a manner well-adapted to a finite volume scheme, namely a conservative form. Here, we write a similar conservative form, but as a cross derivative with respect to both variables. For suitable numerical scheme, we also refer to [Filbet2004] for LS, to [Goudon2013a] for LS with encounters and to [Goudon2012] for a model with space diffusion. Then the rest of the section is devoted to the proof of the convergence result given in Theorem 3.4. It is based on weak compactness of sequences of approximations. We emphasis that the numerical scheme introduced here takes its originality in the introduction of the cross-derivative for the coagulation operator. The main difficulty reside in its approximation which involves new conservation. Moreover, we get a better regularity in the weak stability principle result than in [Bourgade2008], which may also apply to their case.

In section 4, we present a numerical simulation which is produced by the numerical scheme. Based on this simulation, we discuss and interpret the long-time behaviour of the model and we show how it might be related to the behaviour of a non-autonomous coagulation equation. Here, it might be very interesting to connect this problem to [Herrmann2012, Collet2002, Gabriel2012] on the LS and related model or [Escobedo2005, Desvilettes2009] for the coagulation equation.

## 2. Rescaling and existence of global solutions

### 2.1. Rescaled problem

The nature of the configurational space is not really convenient for computations both for the theoretical point of view and the numerical implementation. Thus, we decide to rescale the problem and operate a change of variable in the distribution with respect to a physically relevant new variable . We call it the ion ratio (without dimension) and its definition is, for ,

We introduce then the new unknown defined, over the new configuration space , by

Then, we operate a a change of variable in the association-dissociation and coagulation rates, by introducing define over and over , such that

Thus, they satisfy

(2.1) |

with and . Also, the symmetry assumption (1.6) becomes

(2.2) |

A formal computation leads to the assessment

Finally, letting such that

(2.3) | ||||

(2.4) |

with , we get

In the following and for the rest we drop tildes in the rescaled problem, for sake of clarity. Now, we are able to reformulate the problem which is to find the distribution satisfying

(2.5) |

with the constraint

(2.6) |

and boundary condition (1.9) remains given by

(2.7) |

while the initial conditions is only a change of variables in (1.10):

(2.8) |

### 2.2. Hypotheses and result

The study of the problem (2.5)-(2.6) with (2.7) and (2.8) requires hypothesis whether they naturally arise in the problem or technicals. Namely, we assume that:

H1. The initial distribution is nonnegative and such that

(2.9) |

H2. The coagulation rate is nonnegative, satisfies (2.2) and

(2.10) |

H3. The rates functions are both nonnegatives and

(2.11) |

and for all ,

(2.12) |

H4. For all and ,

(2.13) |

and

(2.14) |

Here, denotes a constant. Note that (2.13) ensures the characteristics remain in the set and allows us to prescribe the boundary (2.7). In fact, it is equivalent with respect to (2.1) and (H4) to assume

(2.15) |

for all .

###### Remark 2.1.

With such variables, we note that example (1.5) becomes

And, hypothesis (H4) is consistent with this example.

###### Definition 2.1 (weak solution 1).

We remark here that regularity (2.16), where means continuous from to a Banach space with respect to the weak topology of . Hypothesis (H1) to (H3) suffice to define (2.17). Particularly, (2.3)-(2.4) entail, as we will see later, that belongs to .

We can now state the main result:

###### Theorem 2.2 (Global existence).

Proof of Theorem 2.2 relies on 2 main steps, which are similar to the ones used for instance in [Collet1999] and [Collet2000] for LS equation. The first step consists in the construction of a mild solution of equation (2.5) for a given . This is achieved through a fixed point theorem by virtue of the contraction property of the coagulation operator. The second step follows a second fixed point which associates (2.5) to the constraint (2.6) on .

Since the method is rather classical, we only provide in the next section the key arguments of the proof, by highlighting the differences between our problem and LS equation with encounter. Particularly, the treatment of the characteristics.

### 2.3. Existence of solutions

#### 2.3.1. The autonomous problem.

We start the analysis of the problem for a given nonnegative with , i.e. we avoid the difficulty induced by the constraint (2.6). A well-know approach is to construct the characteristics of the transport operator. These are the curves parametrized by and associated to given, for any , by the solution of

According to (2.11)-(2.12) and (2.15), there exists a unique solution . We only consider the characteristic while they are defined, i.e . We would first remark that (H4) ensures the characteristics remain into for any . Moreover, we define the origin time . So, from the characteristics curves we construct the so-called mild-formulation which is solution of

(2.18) |

for all and a.e. and

Note that, for an enough regular solution, the boundary condition (2.7) is satisfy by (2.18) since .

Then, to prove Theorem 2.2 we need to recover the notion of weak solution from the one of mild solution. This is given by the following the result:

###### Lemma 2.3.

This result is well known and comes from a change of variable and an identification process, we refer to [Collet1999, Collet2000] for LS equation or [Diperna1989] for Boltzmann equation. The only delicate point remain in the treatment of the origin time and to this end we should proceed as in [Helal2013]. The monotonicity hypothesis (2.14) is crucial to separate continuously the characteristics coming from and and hence to be able to construct the weak solution from (2.18). Now, with the help of this lemma, it is sufficient to prove the existence of a mild solution.

Before claiming the existence of mild solution for a given , let us introduce some a priori properties of the coagulation operator. Namely, for any and both belongs to , we have

(2.19) | |||

(2.20) |

These two estimates ensure that maps into itself and is Lipschitz on any bounded subset of . Finally, we would remark that for any and , it holds

(2.21) |

with , which is called weak formulation of the coagulation operator, see [Desvilettes2009, Escobedo2003, Laurencot2002]. We obtain this identity by inversion of integrals applying Fubini’s theorem, then changes of variable. In particular, when ,

(2.22) |

And, if moreover , then

(2.23) |

Now, for a given , we define the set:

and we are ready to claim the next proposition.

###### Proposition 2.4.

Let and with belongs to the associated set . If , then there exists a unique nonnegative mild solution, i.e. satisfying (2.18), with

Moreover, for all we have

(2.24) |

and

(2.25) |

###### Proof.

Here we only give a sketch of the proof.

Step 1. Existence and uniqueness. The local existence of a unique nonnegative solution , for some small enough, readily follows from the Banach fixed point theorem applied to the operator that maps to the right-hand side of (2.18) on a bounded subset of . To that, we follow line-to-line [Collet1999] using properties (2.19) and (2.20).

Then, the global existence, for any time , is obtained using estimation (2.24), indeed by a classical argument we construct a unique solution on intervals , , etc. So, it remains to prove (2.24), which directly follows from the integration of (2.18), using that and (2.22).

Step 2. Mass conservation. It remains to prove (2.25) which needs and too. But, identity (2.21) holds only for and a priori is not integrable. It is enough to follow [Collet1999]*Lemma 4 which involves a regularization procedure using in (2.18) and construct a sequence of approximation . Then, computing the norm of with the mild formulation yields to the strong convergence of toward the solution in when , since a.e. . Finally, because has a compact support, is integrable and from (H2) we get the uniform bound (in ):

for some constant obtained by a Gronwall’s lemma on the mild formulation. We get and is integrable. We conclude coming back to (2.18) and integrating it against , which yields (2.25) thanks to (2.23). ∎

###### Remark 2.2.

The boundedness of is a key ingredients in the estimation used in the proof above. Indeed, without it the control of the mass could break down, i.e. in finite time , known as gelation phenomena. Nevertheless, the condition could be relaxed, generally up to a sub-linear coagulation kernel, see for instance [Escobedo2003].

We close this section by stating an additional regularity on the pseudo-moment of the mild solution, key argument to couple the constraint (2.6) to (2.5).

###### Corollary 2.5.

###### Proof.

We are now ready to apply the second fix point to connect and in the system.

#### 2.3.2. Fix point on .

Again we follow [Collet1999, Collet2000], i.e we let and we define the map

where is the positive part and the unique mild solution on associated to thanks to Proposition 2.4. It follows that maps into itself. Moreover, using (2.26)

then, since is Lipschitz, it holds that with

a.e. , see [Ziemer1989]*Theorem 2.1.11. Thus, for any by using hypothesis on the rates (H3), it yields

Next, we invoke Ascoli theorem to claim that is relatively compact in . Now, a Schauder fix point theorem would achieve the proof of Theorem 2.2. It remains to prove the continuity of the map . Let be a sequence of converging to for the uniform norm. We need to prove that

which is done by estimating

Indeed, the right hand side of this inequality goes to zero following line-to-line the proof of [Collet1999]*Lemma 5 and 6 to conclude on the one hand the continuity and on the other that, in fact,

to drop . Thus, there exist such that

This achieves the proof of Theorem 2.2.

## 3. Numerical approximation

### 3.1. A conservative truncated formulation

The discretization of the problem (1.3)-(1.4) gives rise to three main difficulties. First, the unboundedness of the space. Indeed, one of the two variables has been reduced to the interval , with a physical meaning, but the -variable can reach any size in . Thus, we decided to proceed as in [Bourgade2008] and carry out a truncation of the problem considering a “maximal reachable size”, or cut-off, . The link between both problems, truncated and full, when is not taken in consideration here. The reader can refer to step 2 in the proof of Proposition 2.4 to get some hints related to this topic. The purpose is to provide a converging numerical approximation of a truncated problem for a fixed . The second issue arises when we look toward conservations of the system. In [Bourgade2008], the authors propose a reformulation of the coagulation operator into a divergence form. We are inspired by this method and adapt it to our problem. Indeed, this formulation appears natural for finite volume scheme and has the advantage to provide exact conservations at the discrete level.

The starting point is the weak formulation of , see (2.21). One can take for some and we formally get an expression of the form

where the coagulation reads now

(3.1) |

where . Now the coagulation operator has been reformulated in a manner well adapted to a finite volume scheme (the volumes averages of are brought out directly). It remains to truncate the problem. It can be achieved in two different ways as mentioned in [Bourgade2008, Filbet2004], the authors discuss about conservative and non-conservative form. These two options can be derived respectively by taking or . The first option avoids the formation of clusters larger than thus it will preserve the mass, while the second induces a loss of polymers due to the creation of larger clusters than . This latter is convenient to study gelation phenomenon, see [Escobedo2003] for a review on coagulation. Here, we restrict ourself to the conservative form and obtain the truncated operator by taking

Then, replacing by in (3.1) it yields, for any , to

(3.2) |

The last main issue lies in the discretization of (1.4), i.e. the algebraic constraint driving . We will not be able to properly derive an approximation of in the coagulation operator which would allow us to control the sign of . Once again, we reformulate this constraint obtaining an evolution equation on by a time derivation of it. Thus the problem (1.3)-(1.4) reads now

(3.3) |

and

(3.4) |

The boundary condition reads now

(3.5) |

and the initial data are

(3.6) |

It is now appropriate to introduce the technical assumptions used through this section:

H1’. The initial distribution is nonnegative and , with

(3.7) |

H2’. The coagulation rate is nonnegative and

(3.8) |

H3’. The rate functions are nonnegatives and

(3.9) |

H4’. The function is a non-increasing function:

(3.10) |

###### Remark 3.1.

We emphasize that hypotheses H2’ and H3’ are not so restrictive in front of the truncation, it could allow unbounded rate on the full configuration space locally bounded which seems reasonable.