Jocelyne Erhel

Bastien Hamlat

Tangi Migot

Keywords: chemical kinetics, discontinuous ODEs, differential inclusions, projected dynamical systems, Filippov's theory, numerical simulation, analysis of positivity and existence

Une page se ferme en écrivant les dernières lignes de ce manuscrit... Je souhaite à tous les futurs étudiants doctorants d'arriver à s'épanouir pendant et après la thèse. Je garde en souvenir une tendresse pour chacune de mes rencontres pendant ces années de recherche et de rédaction.

Jocelyne ERHEL, Anthony MICHEL et Thibault FANEY, je tiens à vous exprimer mes plus sincères remerciements pour votre implication durant la thèse. Je vous remercie de m'avoir accordé votre confiance et beaucoup de liberté scientifique. Votre expérience et vos conseils m'ont été précieux. Vous avez tous fait part d'une parfaite bienveillance à mon égard et je vous en suis reconnaissant. Jocelyne, je te remercie d'avoir accepté de diriger ma thèse. Ce fut un réel plaisir d'apprendre à te connaître et de travailler avec toi. Mes visites à Rennes ont toujours été des moments agréables et très fructueux scientifiquement. Tes retours et tes apports m'ont permis de progresser au cours de ce long travail. Je garderai un très bon souvenir de nos échanges. Anthony, je te remercie de m'avoir recruté et pour avoir eu confiance en moi dès notre première rencontre. Notre travail collaboratif parfois quotidien a été moteur dans cette thèse. Tu as toujours été disponible lorsque j'avais besoin d'aide. Thibault, ton soutien et tes conseils scientifiques m'ont été indispensables durant la thèse. Je te remercie pour ta bonne humeur communicative, surtout durant les périodes difficiles.

Pour avoir accepté de faire partie du jury de la thèse je souhaiterais aussi remercier sincèrement la présidente du jury Danielle HILHORST, les rapporteurs Étienne AHUS-BORDE et Juliette LEBLOND pour leurs retours qui ont permis l'amélioration du manuscrit et les examinateurs Mounir HADDOU et Michel KERN pour l'intérêt qu'ils ont accordé à mes travaux de thèse.

Je remercie Jocelyne de m'avoir présenté Tangi MIGOT et de nous avoir permis de travailler ensemble. J'ai eu un réel plaisir à travailler avec toi, tu es toujours enthousiaste et la difficulté ne te fait pas peur bien au contraire. Merci d'avoir partagé avec moi tes compétences scientifiques. Bien évidement, je te remercie d'avoir travaillé à mes côtés et je te souhaite le meilleur pour ta carrière scientifique. Et puisque tout est permis : que Rennes soit champion de France. Durant ma thèse, j'ai eu le plaisir et la chance de faire mes entretiens annuels d'avancement du projet de thèse, organisés par l'école doctorale MathSTIC que je remercie très sincèrement, avec Michel KERN et Jean-Raynald DE DREUZY. Cela était bien plus qu'un simple entretien où l'on présente ses travaux. Nos discussions m'ont permis d'explorer des axes différents dans mes travaux et parfois de me rassurer lorsque les inquiétudes me tenaillaient. Ma thèse est une collaboration entre IFP Energies nouvelles et INRIA. Je souhaite remercier IFP Energies nouvelles de m'avoir recruté en thèse, de m'avoir octroyé tout le confort matériel et financier dont a besoin un jeune chercheur pour effectuer son travail. Je remercie également INRIA de Rennes d'avoir financé plusieurs de mes déplacements et de m'avoir toujours parfaitement bien accueilli lors de mes venues. Je souhaite remercier l'école doctorale MathSTIC. Merci à l'université Paris 2 Panthéon-Assas d'avoir financé la dernière année de ma thèse en m'octroyant un contrat ATER.

Je souhaite remercier la région Bretagne pour le financement de mon voyage à Milan dans le cadre de mes travaux de recherche. Je remercie le MOX de Politecnico di Milano pour m'avoir accueilli dans leur laboratoire. Merci à Luca FORMAGGIA et Luca BONAVENTURA de m'avoir permis de travailler dans leur équipe. Je remercie plus particulièrement Anna SCOTTI qui m'a supervisé pendant mon séjour (mathématiques, activités, problème de dents,...). J'ai été chanceux de la rencontrer, un grand merci. Je souhaite remercier Abramo AGOSTI pour les discussions scientifiques et l'aide qu'il a pu m'apporter. Une pensée pour tous les stagiaires, doctorants et post-doctorants : Federica, Ludovica, Francesco, Fannie, Laura, les deux Daniele, Alberto, Mattia, Matteo, Abele. Nous avons passé tant de bons moments ensemble (aperitivo + football) et vous m'avez même appris à faire le Tiramisu ! Enfin, je remercie Florent CHAVE, mon compère Montpelliérain, sans toi l'aventure Milanaise n'aurait pas été la même. Je te souhaite de réussir ce que tu entreprends.

Je tiens à remercier les membres de OSUR pour la summer school de Cargèse et plus particulièrement les doctorants que j'ai eu le plaisir de revoir plusieurs fois lors de mes séjours rennais.

Je souhaite remercier Terence BAYEN pour les discussions scientifiques que l'on a eu durant le congrès SMAI. J'ai été parfaitement accueilli dans l'équipe Fluminance de INRIA Rennes. Merci à Etienne MEMIN le chef de projet. Je tiens particulièrement à remercier Huguette BE-CHU, avec qui j'ai beaucoup échangé, pour sa disponibilité, son efficacité et d'avoir fait de son mieux pour me permettre de me déplacer dans les meilleures conditions. Merci aux doctorants, post-doctorants et permanents de l'équipe pour les moments passés ensembles.

Je remercie tous les membres de R115 et R114. Une pensée particulière pour Zakia BENJELLOUN TOUIMI qui a accordé une réelle attention à ma thèse et qui a été pour moi un soutien indispensable, simplement merci. Je tiens à remercier Carole GUENAUD et Sylvie HOGUET qui m'ont toujours aidé et qui font un travail extraordinaire. Elles ont toujours été bienveillantes et agréables avec moi. Sylvie, je garde un souvenir particulier de nos discussions, j'ai appris à connaître une personne formidable. Je remercie Huy TRAN QUANG pour ses conseils scientifiques et son aide, les doctorants de R115/114 ont de la chance de pouvoir compter sur toi. J'ai une pensée pour chaque membre permanent de IFP Energies nouvelles que j'ai eu le plaisir de rencontrer : Abdallah, Ani, Anthony, Aziz, Benjamin, Bruno, Catherine, Chakib, Christophe, Delphine, Elsie, Françoise, Frédéric, Guillaume, Henry, Ibtihel, Isabelle, Jef, Jean-Marc, Jean-Yves, Jean-Louis, Julien, Huy, Léo, Long, Laurent, Mathieu, Miguel, Nina, Olivier, Paul, Pierre, Raphaël, Soleiman, Stéphane, Sylvain, les trois Sylvie, Thibault, les deux Thomas, Xavier, Zakia. C'était un plaisir de partager des instants avec chacun d'entre vous pendant les repas ou durant les pauses de convivialité. J'ai eu le plaisir d'animer le club Belote de IFP Energies nouvelles et pour cela je remercie mes collègues et le CE pour nous avoir permis de nous réunir tous les vendredi midi pour jouer à la belote. J'y ai rencontré des gens formidables et très sincèrement vous m'avez permis de tenir dans les moments difficiles. Vous êtes et resterez dans mon coeur, sachez que même si nos routes se séparent vous pourrez toujours compter sur moi. En tant que premier président du club, je vous demande de ne jamais arrêter ! Jean-Yves comment oublier nos discussions sur la musique, notre goût pour Motörhead nous a tout de suite rapprochés, tu m'as quand même fait redécouvrir Scorpions. Cela m'a touché quand tu m'as invité à joueur aux cartes avec ta famille dans ma région natale. Je suis fier de te compter parmi mes proches. Christophe, j'ai adoré parler de tout avec toi (musique, démineur, maths,...). Nous avons partagé ensemble notre passion commune des problèmes mathématiques et cela était très plaisant, tu es une personne formidable et j'ai aimé chaque instant passé à tes cotés. Sylvie P., cette page ne suffirait pas pour te témoigner ma gratitude. Tu as été un modèle qui m'a permis d'avancer. Tu es une de mes plus belles rencontres durant la thèse. Delphine, la ! maman " des thésards. Tu es un bonheur quotidien. Ta bonne humeur et ta joie de vivre m'ont fait passé tellement de moments agréables durant cette thèse. Bruno, l'indispensable à la bonne ambiance. J'ai adoré jouer à MTG avec toi et discuter de tout. Tu es une super rencontre. Thomas, j'adore ton humour. Franchement tu es au top, ne change pas ! Xavier et Raphaël, le tandem à abattre à la coinche. Xavier, j'ai adoré discuter avec toi et rire de tout et de rien, tu es un gars super sympa. Raph, c'était vraiment un plaisir de partager pleins de moments avec toi, tu m'as beaucoup fait rire. Fred le Nîmois, nos racines communes gardoises nous ont rapprochés, j'ai passé de bons moments dans votre bureau (bureau de Delphine et Fred). Benoît, tu étais le seul à ne pas faire partie de notre département. Pourtant, tu nous as régalés avec tes règles bizarroïdes qui rendent un jeu que nous aimons tous différent mais encore plus fun ; tu n'étais pas non plus le dernier pour mettre l'ambiance. Long, ta bonne humeur me manquera. Sylvie W., je te remercie pour ton aide, j'ai eu beaucoup de chance de te compter parmi mes soutiens. J'ai une pensé pour tous les stagiaires qui ont croisé ma route à IFPEN : Irene, Nawel, Anis, Lynda , Marianne, Imane, Hassan, Olivier, Jeremy,... Une pensée pour Gouranga qui était en Post-doc à IFPEN. J'ai partagé beaucoup de moments avec toi et j'ai adoré nos discussions sur nos différences culturelles. Comment oublier mes collègues doctorants... Je tiens d'abord à tous vous dire que vous êtes de magnifiques rencontres. Mon année de recrutement a vu 3 doctorants être accueillis dans notre département : Julien, Karine et moi-même. Mes deux compères avec qui nous avons traversé tant d'épreuves. Karine, ce fut un plaisir de discuter avec toi. Nous avons passé des bons moments ensemble (sans parler de la pétanque !). Je te souhaite le meilleur pour la suite, tu le mérites. Julien, une sincère amitié est née entre nous. Tu fais partie des personnes importantes dans ma vie aujourd'hui. On a partagé un bureau pendant de longues années... mais j'espère qu'on continuera à partager une amitié tout le reste de notre existence. J'espère que tu t'épanouiras dans tes futures réalisations. Je tiens à remercier les doctorants qui m'ont accueilli. En particulier, Riad avec qui j'ai partagé tellement de bons moments ! à la pétanque " entre autres. Zakariae, nous avons partagé le même bureau durant toute la durée de ma thèse. J'ai trouvé en toi un ami et tu fais par-tie, aujourd'hui, des personnes importantes de ma vie. Tu es une personne exceptionnelle et je m'estime chanceux de t'avoir rencontré. Thói, j'ai eu le plaisir de partager le bureau avec toi. Nos échanges français/anglais et culturels étaient très enrichissants. Tu étais le ! smile " de notre bureau. Tu es une très belle rencontre et je te souhaite le meilleur. Sabrina, tu as été pour moi la meilleure arrivée dans notre groupe de doctorants. Tu t'es parfaitement intégrée et tu as permis d'accroître la bonne ambiance collective. Toujours souriante, prête à aider, active pour organiser des activités, . . . Tu es une personne sur qui on peut compter. Bien que je ne sois plus à IFP, je serai toujours là si tu as besoin. Guissel, tu as été une super rencontre. Ta joie de vivre et ton humour nous ont permis de passer des moments inoubliables. Alexis, pourtant un peu plus discret que les autres boute-en-trains, nous avons cependant beaucoup discuté et échangé sur des sujets divers et variés. Je te souhaite de réaliser tes futurs projets. Joëlle, j'ai aimé discuter avec toi. Tu es une super rencontre et je te souhaite sincèrement de réussir ce que tu entreprends. Son, nous n'avons que peu échangé mais c'était toujours un plaisir de parler avec toi. Ruben, chacune de tes venues était un plaisir. Narges, nous avons partagé une superbe expérience en Corse ensemble. Ce fut un plaisir de te rencontrer. Pour tous les doctorants qui n'ont pas encore soutenu je vous souhaite que ce jour arrive le plus rapidement possible et je suis sûr que vous ferez de bons docteurs.

Bien évidemment cette aventure n'aurait pas été la même sans mes proches. Une première pensée va à mes parents et mes soeurs. Mathilde, je suis extrêmement fier de toi. Tu es la soeur dont rêvent tous les frères. Marie, tu es encore jeune mais je te souhaite de réussir et surtout de t'épanouir. Madie, j'ai une pensée tendre pour toi. J'ai passé un agréable moment en ta compagnie lors de ma venue sur Lyon. Papa et Christelle, merci pour vos conseils et votre soutien. Maman et Yannick, merci pour votre entière dévotion et vos conseils. Maman merci pour l'amour que tu m'as apporté tout au long de ma vie. Je tiens à remercier mes grand-mères et les autres membres de ma famille. Une pensée particulière pour ma tante Annie qui est venue voir ma soutenance (cela m'a énormément touché). Je souhaite remercier Nicolas de m'avoir accueilli chez lui et de me considérer comme un membre de sa famille. Je suis triste de ne plus venir sur Rennes... Merci à mes amis proches qui m'ont soutenu durant cette épreuve : Alicia, Maïlys, Thibault, Julien, Olivier. Vous savez tous à quel point vous êtes importants dans ma vie.

Mes dernières pensées vont à la femme qui m'a épaulé, supporté et aimé. Amina nous y sommes arrivés... Je ne parle pas que de la thèse. Je suis heureux d'être celui que tu as choisi, tu combles de bonheur chacun de mes jours.

Résumé

Le sujet de thèse s'inscrit dans la continuité des travaux de thèse antérieurs réalisés à IFP Energies nouvelles ( [START_REF] Haeberlein | Time Space Domain Decomposition Methods for Reactive Transport, Application to CO2 Geological Storage[END_REF] et [START_REF] Tillier | Couplage réactions-transport pour la modélisation et la simulation du stockage géologique de CO2[END_REF]) et à INRIA ( [START_REF] De Dieuleveult | Un modèle numérique global et performant pour le couplage géochimie-transport[END_REF]). Tous ces travaux tentent d'apporter des réponses en terme de modélisation mathématique et de méthode numérique pour la simulation du transport réactif en milieu poreux.

En transport réactif les apparitions et les disparitions de phases et d'espèces engendrent un problème de modélisation des cinétiques chimiques. Des travaux sur les cinétiques chimiques ont déjà été réalisés pour les systèmes chimiques comprenant une unique phase minérale pure par réaction. Notre objectif est d'élargir la gamme des systèmes chimiques que l'on peut modéliser et de proposer des méthodes numériques de résolution.

Nous considérons un système chimique avec une phase aqueuse et des phases pures, où les réactions chimiques sont cinétiques. Un modèle avec des vitesses de réaction continues convient en phase aqueuse mais ne convient pas dès qu'il y a une phase pure, car les quantités peuvent devenir négatives. Une difficulté est de maintenir les quantités de matière positives lorsqu'une espèce disparaît.

Dans le premier chapitre, un modèle de type EDOs discontinues pour la cinétique chimique avec apparitions et disparitions d'espèces pour un nombre quelconque de minéraux est proposé. Ce travail a été réalisé en collaboration avec Jocelyne Erhel, Thibault Faney et Anthony Michel. Une version régularisée du modèle permet de prouver la positivité et l'existence de solutions. De plus, une analyse explicite est effectuée dans le cas contenant une espèce réactive intermédiaire. Enfin, une méthode numérique basée sur la formulation régularisée a été validée sur des exemples de systèmes réactifs.

Dans le deuxième chapitre, une reformulation du modèle de cinétique chimique proposé dans le premier chapitre à l'aide de la théorie de Filippov est fournie. Une preuve de l'existence et de la positivité des solutions de la formulation inclusion différentielle ainsi obtenue est réalisé. De plus, dans le cas des surfaces de discontinuité de codimension 1, une étude des configurations des champs fournit un résultat d'unicité et de caractérisation des trajectoires.

Dans le troisième chapitre, un modèle de cinétique chimique de type systèmes dynamiques projetées est proposé. Ce travail est effectué en collaboration avec Jocelyne Erhel et Tangi Migot. Une analyse de l'existence des solutions de ce modèle, des liens avec d'autres types de formulations et une méthode de résolution numérique adaptée sont présentés. Enfin, une illustration des résultats numériques obtenus est réalisée pour des systèmes de cinétique chimique.

Mots-clés : cinétiques chimiques, EDOs discontinues, inclusions différentielles, systèmes dynamiques projetées, théorie de Filippov, simulations numériques, analyse de positivité et d'existence.

Abstract

The thesis is a continuation of the previous thesis work carried out at IFP Energies nouvelles ( [START_REF] Haeberlein | Time Space Domain Decomposition Methods for Reactive Transport, Application to CO2 Geological Storage[END_REF] et [START_REF] Tillier | Couplage réactions-transport pour la modélisation et la simulation du stockage géologique de CO2[END_REF]) and at INRIA ( [START_REF] De Dieuleveult | Un modèle numérique global et performant pour le couplage géochimie-transport[END_REF]). All of these works attempt to provide answers to efficient modeling and numerical methods for resolving reactive transport in porous media.

In reactive transport, the appearance and disappearance of phases and species give rise to a problem in modeling chemical kinetics. Works on chemical kinetics has already been carried out for chemical systems comprising a single pure mineral phase per reaction. Our aim is to extend models and numerical methods to general chemical systems.

We consider a chemical system with an aqueous phase and pure phases, where the chemical reactions are kinetic. A model with continuous reaction rates is suitable in the aqueous phase but no longer appropriate as soon as there is a pure phase, because the quantities can become negative. One difficulty is to maintain the quantities of matter positive when a species disappears.

In the first chapter, a discontinuous ODEs model for chemical kinetics with appearance and disappearance of species for any number of minerals is proposed. This work is made in collaboration with Jocelyne Erhel, Thibault Faney and Anthony Michel. A regularized version of the model allows to prove the positivity and the existence of solutions. In addition, an explicit analysis is performed in a case with an intermediate reactive species. Finally, a numerical method based on the regularized formulation has been validated on examples of reactive systems.

In the second chapter, a reformulation of the chemical kinetics model proposed in the first chapter using Filippov's theory is provided. A proof of the existence and the positivity of the solutions of the differential inclusion formulation thus obtained is produced. In addition, in the case of the discontinuity surfaces of codimension 1, a study of the configurations of the fields provides a result of uniqueness and characterization of the trajectories.

In the third chapter, a model of chemical kinetics of the projected dynamical system is proposed. This work is carried out in collaboration with Jocelyne Erhel and Tangi Migot. An analysis of the existence of solutions of this model, a link with other types of formulations and an adapted numerical resolution method are provided. An illustration of the numerical results obtained is made for chemical kinetic systems.

aqueuses A i , i 1, ¤ ¤ ¤ , N a minérales M i , i 1 N a , ¤ ¤ ¤ , N a N p -Réactions chimiques : R j , j 1, ¤ ¤ ¤ , N r -Stoechiométrie : S ps i,j q 1¤i¤Na Np
¡ j ti t1, ¤ ¤ ¤ , N a u, s i,j 0u et M ¡ j ti t1 N a , ¤ ¤ ¤ , N a N p u, s i,j 0u -Produits de R j (indices) : A j ti t1, ¤ ¤ ¤ , N a u, s i,j ¡ 0u et M j ti t1 N a , ¤ ¤ ¤ , N a N p u, s i,j ¡ 0u -Réaction R j (voir (1.1)) : i A ¡ j |s ij |A i ķM ¡ j |s kj |M k Õ i A j |s ij |A i ķM j |s kj |M k .
-Quantités (concentrations) : c pc i q 1¤i¤Na Np pc a , c p q avec comme contrainte à respecter c ¥ 0 espèces aqueuses c a pc i q 1¤i¤Na espèces minérales c p pc i q 1 Na¤i¤Na Np -Activités :

minéraux a i 1, 1 N a ¤ i ¤ N a N p espèces aqueuses a i pc a q, 1 ¤ i ¤ N a avec (Assumption 1.2.1)
c a ¥ 0 ñ apc a q ¥ 0 1 c a,i 0 ñ a i pc a q 0.
-Loi d'action de masse ou taux "maximal" de réaction : τ j pour la réaction R j ,

j 1, ¤ ¤ ¤ , N r (voir (1.3)) τ j pc a q α j ¹ iA ¡ j a i pc a q |s ij | ¡ β j ¹ iA j a i pc a q |s ij | , j 1, . . . , N r .
-Vitesse de réaction : r j pour la réaction R j , j 1, ¤ ¤ ¤ , N r -Système cinétique (voir (1.2)) conservation de la masse dc dt Srpcq contrainte de positivité des quantités c ¥ 0, condition initiale cp0q c 0

-Invariants : Q base de kerpS T q (S T Q Q T S 0) c solution (satisfait la conservation de la masse) implique Q T c Q T c

Introduction

Cette thèse s'inscrit dans la continuité des travaux de thèse antérieurs réalisés à IFP Energies nouvelles ( [START_REF] Haeberlein | Time Space Domain Decomposition Methods for Reactive Transport, Application to CO2 Geological Storage[END_REF] et [START_REF] Tillier | Couplage réactions-transport pour la modélisation et la simulation du stockage géologique de CO2[END_REF]) et à l'INRIA ( [START_REF] De Dieuleveult | Un modèle numérique global et performant pour le couplage géochimie-transport[END_REF]). Ces travaux tentent d'apporter une réponse aux difficultés liées à la modélisation et à la simulation numérique du transport réactif en milieu poreux. Ce domaine a un grand nombre d'applications industrielles : EOR-CO 2 (production assistée d'hydrocarbures par injection de CO 2 ), stockage du gaz (injection de gaz dans le sous-sol afin de le conserver et de l'utiliser pour répondre à une demande ultérieure) ou encore géothermie (circulation d'un fluide dans le sous-sol afin de récupérer de l'énergie sous forme de chaleur). Ces applications sont basées sur l'écoulement d'une ou plusieurs phases fluides (eau, gaz, huile) dans une roche poreuse composée de phases minérales pouvant engendrer des interactions fluide-roche. Il est ainsi essentiel de disposer de formulations mathématiques bien posées et de solveurs numériques efficaces permettant de modéliser la complexité physique des systèmes étudiés.

Modèle de transport réactif multiphasique

Dans cette thèse, nous proposons des modèles mathématiques permettant de modéliser la cinétique chimique des interactions fluide-roche. Nous présentons ici la structure d'un problème de transport réactif multiphasique. L'objectif est de décrire les inconnues et les paramètres du système, ainsi que les équations à résoudre.

Description du système physique

Dans l'étude du transport réactif multiphasique en milieu poreux, il faut distinguer les lois et paramètres physiques issus de la physique de l'écoulement et ceux issus des interactions chimiques. Le système compositionnel d'un milieu poreux, noté Σ, permet de décrire la composition chimique du milieu poreux. Il est composé de deux sous-systèmes caractérisés par les états de matière : le sous-système fluide F Σ, comprenant les liquides et les gaz, et le sous-système solide S Σ, comprenant les solides. Chaque soussystème est lui-même constitué de phases notées α qui sont des ensembles de constituants chimiques pour lesquels les propriétés physiques sont homogènes. Les phases du soussystème solide (i.e. α S) considérées dans ce document sont pures, c'est à dire qu'elles ne contiennent qu'un unique constituant chimique. En d'autre termes, chaque minéral constitue une phase du sous-système solide. Chaque phase est composée de plusieurs constituants chimiques appelés espèces chimiques. On note E l'ensemble des espèces chimiques du système compositionnel Σ et E α l'ensemble des espèces chimiques contenues dans la phase α. Lorsqu'on utilise des espèces abstraites E E, au sens où l'on ne

Lois de fermeture

Avant de préciser les équations qui modélisent les phénomènes physiques et chimiques dans un milieu poreux, on introduit les équations de fermeture. Ces équations sont structurelles, elles permettent de garantir la cohérence du système d'inconnues. Equation de fermeture du système :

φ S φ F 1, (1)
l'équation de fermeture du système (1) permet de maintenir constante la somme de la proportion de fluide et de solide dans le système.

Equation de fermeture du fluide :

αF S α 1, (2)
Equation de fermeture du solide :

αS F α 1, (3)
les équations de fermeture du fluide (2) et du solide (3) permettent de maintenir la somme des proportions des phases qui les composent à 1. Equations de fermeture des phases :

iα x α i 1, dα S F, (4)
les équations de fermeture des phases (4) permettent de maintenir la somme des fractions molaires de chaque phase à 1. Le nombre d'équations de fermeture correspond au cardinal de l'ensemble des phases auquel on additionne trois. Il reste donc le cardinal de l'ensemble des espèces plus une équation à ajouter pour égaliser le nombre d'inconnues.

Physique de l'écoulement en milieu poreux

Les problèmes de transport réactif en milieu poreux font intervenir plusieurs lois physiques qui modélisent l'écoulement du fluide. Le système d'équations résultant est un système d'équations différentielles aux dérivées partielles (EDP). Il s'agit plus précisément d'une équation d'évolution avec un terme d'advection et de diffusion. Les paramètres suivants sont nécessaires pour définir la vitesse de Darcy [START_REF] Darcy | Les fontaines publiques de la ville de Dijon[END_REF] dans la phase α notée v α (m ¤ s ¡1 ) : le tenseur de perméabilité absolu K (m 2 ), la perméabilité relative associée à une phase kr α (sans unité), la viscosité d'une phase µ α (Pa ¤s), l'accélération de la gravité g (m ¤ s ¡2 ), la masse volumique d'une phase ρ α (kg ¤ m ¡3 ), la masse molaire d'une phase M α (kg ¤ mol ¡1 ) et la pression capillaire d'une phase P c α (Pa).

On peut définir la loi de Darcy comme suit, dα F,

v α ¡K kr α µ α ∇pP P c,α q ¡ ρ α g ¨. (5)
Un exemple pour les lois que suivent la perméabilité relative et la pression capillaire sont les lois d'hydrostatique de Brooks-Corey et sont définies dans le cas biphasique (phase aqueuse w et phase gazeuse g) comme suit,

S w S w ¡ S w,r 1 ¡ S w,r ¡ S g,r , (6)
où S α,r est la saturation résiduelle de la phase α tw, gu.

P c,g P e,g pS w q ¡ 1 λ , (7)

kr w S 2 3λ λ w , (8)
kr g p1 ¡ S w q 2 S 2 λ λ w , (9)
où P e,g est une constante de pression et λ est égal à 2 dans la plupart des cas. La loi de Darcy permet de définir la vitesse d'advection d'un problème de transport réactif en milieu poreux. L'opérateur de transport L, dépendant de la concentration des espèces chimiques, est défini pour chaque espèce dans une phase. Il est noté L α i pour la ième espèce de la phase α. On note c α i (mol ¤m ¡3 ) la concentration de la i-ème espèce dans la phase α, cette grandeur est proportionnelle à la fraction molaire de l'espèce considérée,

c α i x α i ρ α M α , (10)
et on note D α i (m 2 ¤ s ¡1 ) le coefficient de diffusion de la i-ème espèce dans la phase α.

Opérateur de transport advectif et diffusif :

L α i pc α i q div c α i v α ¡ D α i ∇c α i ¨, α F, (11)
L α i pc α i q 0, α S. (12)
Un autre phénomène physique à prendre en compte est l'apport de matière dû aux puits (producteurs et injecteurs) notés p P, ce terme source qui vient s'additionner au transport est noté Q. Le terme source Q est défini par espèce, on le note Q α i , il est dépendant de la pression et de la concentration.

Termes sources de type puits :

Q α i pP c α i,p T p,α pP P c α ¡ P p,α q d p,i pc α i ¡ c α p,i q ¨, (13)
tel que c α i,p (mol ¤ m ¡3 ) est la concentration de la i-ème espèce dans le p-ème puits, T p,α est la transmissivité, P p,α (Pa) est la pression de sa phase α et d p,i est un coefficient de diffusivité.

Interactions chimiques

Les interactions chimiques considérées sont des réactions chimiques de deux natures différentes : les réactions d'équilibre et les réactions de cinétique. Les réactions d'équilibre modélisent des interactions chimiques instantanées par rapport à l'échelle de temps considérée du phénomène physique. Mathématiquement, les considérer revient à résoudre un problème algébrique issu d'une élimination des variables dites secondaires du système d'EDP initial. Pour plus de détails se reporter aux différents travaux concernant les tables de Morel [START_REF] Morel | Principles and Applications of Aquatic Chemistry[END_REF]. Les cinétiques modélisent les transferts chimiques lents. Mathématiquement, la prise en compte des cinétiques se fait par l'ajout d'un terme source de vitesse de réaction au système de transport. Ces deux types d'interactions peuvent coexister et mathématiquement il n'y a aucune contradiction à considérer les deux simultanément. En effet, il suffit d'inclure les seconds membres liés à la cinétique puis d'effectuer l'élimination des variables et l'ajout des équations algébriques pour les équilibres.

Le système réactif est constitué de l'ensemble des réactions chimiques du problème. Une réaction chimique modélise les transformations de la matière. Dans la suite, l'indice i est rattaché à l'espèce et l'indice j à la réaction. Une réaction chimique se définit par une formule mathématique associant les espèces chimiques aux coefficients de stoechiométrie notés s i,j . Le coefficient stoechiométrique est une proportion d'une espèce consommée ou produite dans une réaction.

Les réactions d'équilibre sont symbolisées par le signe tandis que les réactions cinétiques, considérées comme réversibles, par le signe é.

Exemple 0.1.1 (Système réactif).

2W 1 W 2 M 1 W 3 G 1 3W 4 M 2 é M 3 (14) (15) (16)
Ce système réactif est composé de 8 espèces dont 4 espèces aqueuses, 3 espèces minérales et 1 espèce gazeuse. Il est constitué de 3 réactions dont les deux premières sont des réactions d'équilibre et la dernière est une réaction de cinétique.

Pour définir le processus réactionnel les constantes d'équilibre, notées K i , pour la ième espèce, doivent être fournies. La vitesse de chaque réaction cinétique j est notée τ j . Les vitesses de réaction sont modélisées par la loi d'action de masse cinétique :

τ j k r,j ¹ it1,¤¤¤ ,nu s i,j 0 rE i s |s i,j | ¡ k p,j ¹ it1,¤¤¤ ,nu s i,j ¡0 rE i s |s i,j | . ( 17
)
où k p,j , k r,j sont des constantes positives dépendantes des constantes d'équilibre des espèces intervenant dans la réaction et de la constante de vitesse de la réaction. rE i s sont les activités des espèces chimiques. L'activité d'une espèce chimique est une mesure de l'écart des propriétés de l'espèce chimique entre un état standard et un mélange réel à la même température. Cette quantité est modélisée par une loi fonction de la fraction molaire,

rE i s γ i x i , (18)
où γ i est le coefficient d'activité de l'espèce chimique E i . L'activité des espèces en phase pure est égale à 1.

Les réactions chimiques peuvent être représentées par la matrice de stoechiométrie S. Cette matrice contient les réactions chimiques. Pour cela, elle stocke en colonne les coefficients stoechiométriques de chaque espèce dans la réaction considérée. Ces coefficients sont signés par leur présence à droite ou à gauche de la réaction. Pour la suite on considérera la convention suivante, les espèces réactives (présentes à gauche de la réaction) auront des coefficients négatifs et les espèces produites (présentes à droite de la réaction) auront des coefficients positifs. La matrice de stoechiométrie est la concaténation de la matrice S kin , matrice de stoechiométrie pour les réactions de cinétique, et de la matrice S eq , matrice de stoechiométrie pour les réactions d'équilibre, S pS eq |S kin q.

Exemple 0.1.2 (Matrice de stoechiometrie). La matrice du système réactif de l'exemple 0.1.1 est définie comme suit, on choisit de ranger les espèces de la manière suivante (W 1 ,

W 2 , W 3 , W 4 , M 1 , M 2 , M 3 , G 1 ), S eq ¤ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¥ ¡2 0 ¡1 0 0 ¡1 0 0 1 0 0 0 0 0 0 1 , S kin ¤ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¥ 0 0 0 ¡3 0 ¡1 1 0 and S ¤ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¥ ¡2 0 0 ¡1 0 0 0 ¡1 0 0 0 ¡3 1 0 0 0 0 ¡1 0 0 1 0 1 0 .
La matrice de stoechiométrie nous permet de déduire les quantités de matière conservées. La matrice des quantités de matière conservées est notée Q et se déduit du noyau de la transposé de S : vecpQq kerpS T q. Respectivement, on définit Q kin et Q eq de la même manière.

Nous allons préciser la contribution des réactions d'équilibre, comme énoncé en début de section ces réactions contribuent par l'ajout d'un système algébrique à résoudre. Pour cela il suffit de définir les contraintes à respecter. Elles dépendent du potentiel chimique,

µ (J ¤ mol ¡1 ). µ i µ 0 i RT lnprE i sq, (19)
où µ 0 i est le potentiel chimique de l'espèce E i dans l'état standard à la même température.

Le système d'équations induit par les équilibres chimiques s'obtient en multipliant les équations de conservations par la matrice Q eq et par l'ajout des contraintes linéaires S eq µ 0.

Système d'équations

Dans l'objectif de définir le système d'équations modélisant le transport réactif multiphasique, nous définissons, ci-dessous, le vecteur des quantités de matière n (mol),

n α i φ F S α ρ α M α x α i , α F, (20)
n α i φ S F α ρ α M α x α i , α S. (21)
En assemblant tous les concepts définis précédemment, nous définissons le système d'équations aux dérivées partielles et algébriques modélisant le transport réactif multi-phasique comme suit :

φ S φ F 1, (22) αF S α 1, (23) αS F α 1, ( 24
) iα x α i 1, dα S F, (25)
Q eq ¡ fn ft Lpcq Q S kin τ © 0, (26)
S eq µ 0.

(

) 27

Etat de l'art

Les formulations globalement implicites pour le transport réactif en milieu poreux sont apparues pour des problèmes de transport couplés avec des équilibres de phases. Ces formulations, voir [START_REF] Coats | An equation of state compositional model[END_REF] et [START_REF] Lauser | A new approach for phase transitions in miscible multi-phase flow in porous media[END_REF], permettent de résoudre de manière implicite des problèmes de transport multiphasique avec des phases fluides mélanges et des équilibres de phases. Des formulations incluant les cinétiques ont été proposées dans [START_REF] De Dieuleveult | A global strategy for solving reactive transport equations[END_REF], [START_REF] De Dieuleveult | A global approach to reactive transport: application to the momas benchmark[END_REF], [START_REF] Haeberlein | Time Space Domain Decomposition Methods for Reactive Transport, Application to CO2 Geological Storage[END_REF], [START_REF] Hoffmann | A general reduction scheme for reactive transport in porous media[END_REF], [START_REF] Kräutle | Existence of global solutions of multicomponent reactive transport problems with mass action kinetics in porous media[END_REF] et [START_REF] Tillier | Couplage réactions-transport pour la modélisation et la simulation du stockage géologique de CO2[END_REF]. Ces formulations globalement implicites permettent l'ajout d'un minéral, pour les équilibres et pour les cinétiques, en interaction chimique avec la phase mobile, fluide.

Les formulations de transport réactif multiphasique en milieu poreux sont des EDPs couplées avec des équations algébriques et des contraintes d'inégalités. Une liste non exhaustive des articles proposant un modèle de transport réactif, voir [START_REF] Bildstein | Modelisation geochimique des interactions eau-gaz-roche application a la diagenese minerale dans les reservoirs geologiques[END_REF], [START_REF] Coats | An equation of state compositional model[END_REF], [START_REF] Friedly | Solute transport with multiple equilibrium-controlled or kinetically-controlled chemical reactions[END_REF], [START_REF] Hoffmann | Existence and uniqueness of a global solution for reactive transport with mineral precipitation-dissolution and aquatic reactions in porous media[END_REF], [START_REF] Knabner | An analysis of crystal dissolution fronts in flows through porous media. Part 1: Compatible boundary conditions[END_REF], [START_REF] Kräutle | Existence of global solutions of multicomponent reactive transport problems with mass action kinetics in porous media[END_REF], [START_REF] Rubin | Transport of reacting solutes in porous media: Relation between mathematical nature of problem formulation and chemical nature of reactions[END_REF], [START_REF] Van Duijn | Travelling wave behaviour of crystal dissolution in porous media flow[END_REF], [START_REF] Van Duijn | An analysis of crystal dissolution fronts in flows through porous media part 2: incompatible boundary conditions[END_REF] et [START_REF] Van Duijn | Crystal dissolution and precipitation in porous media: pore scale analysis[END_REF]. L'ajout du multiphasique aux transferts de matière entraîne des possibles disparitions et apparitions de phases. Le système couplé doit maintenir les quantités de matière positives à l'aide des contraintes d'inégalité. La gestion de ces inégalités pose de réelles difficultés dans les simulateurs. Des travaux sur la prise en compte de ces contraintes ainsi que sur l'efficacité et la robustesse des schémas numériques ont été réalisés dans le but de simuler numériquement le transport réactif [START_REF] Amir | A global method for coupling transport with chemistry in heterogeneous porous media[END_REF], [START_REF] Brunner | Analysis of a modified second-order mixed hybrid bdm finite element method for transport problems in divergence form[END_REF], [START_REF] Buchholzer | The semismooth Newton method for the solution of reactive transport problems including mineral precipitation-dissolution reactions[END_REF], [START_REF] Carrayrou | Comparison of numerical methods for simulating strongly nonlinear and heterogeneous reactive transport problems-the MoMas benchmark case[END_REF], [START_REF] Corvisier | Implémentation des phénomènes de germination / mûrissement / croissance des phases solides secondaires dans un modèle de transport-réactif en milieu poreux géologique[END_REF], [START_REF] De Dieuleveult | A global strategy for solving reactive transport equations[END_REF], [START_REF] De Dieuleveult | A global approach to reactive transport: application to the momas benchmark[END_REF], [START_REF] Erhel | Analysis of a global reactive transport model and results for the MoMas benchmark[END_REF], [START_REF] Erhel | Solving partial differential algebraic equations and reactive transport models[END_REF], [START_REF] Haeberlein | Time Space Domain Decomposition Methods for Reactive Transport, Application to CO2 Geological Storage[END_REF], [START_REF] Hoffmann | Reactive transport and mineral dissolution/precipitation in porous media : efficient solution algorithms, benchmark computations and existence of global solutions[END_REF], [START_REF] Hoffmann | A general reduction scheme for reactive transport in porous media[END_REF] [50], [START_REF] Kräutle | Robust simulation of mineral precipitation-dissolution problems with variable mineral surface area[END_REF], [START_REF] Kumar | Convergence analysis for a conformal discretization of a model for precipitation and dissolution in porous media[END_REF], [START_REF] Lauser | A new approach for phase transitions in miscible multi-phase flow in porous media[END_REF], [START_REF] Leal | An overview of computational methods for chemical equilibrium and kinetic calculations for geochemical and reactive transport modeling[END_REF], [START_REF] Lichtner | Reactive transport in porous media[END_REF], [START_REF] Nourtier-Mazauric | Modélisation géochimique et numérique des interactions entre des solutions solides et une solution aqueuse. Extension du logiciel de réaction-transport Archimède et application à la diagenèse des réservoirs[END_REF], [START_REF] Sabit | Numerical methods for reactive transport[END_REF] et [START_REF] Tillier | Couplage réactions-transport pour la modélisation et la simulation du stockage géologique de CO2[END_REF].

La chimie des milieu poreux est de nature complexe et variée. Pour une introduction à la géochimie on se réfère aux livres suivants: [START_REF] Bethke | Geochemical reaction modeling[END_REF] et [START_REF] Lasaga | Kinetic Theory in the Earth Sciences[END_REF]. Les systèmes réactifs peuvent contenir des réactions de type cinétique et de type équilibre. Elles peuvent aussi contenir des espèces présentes dans différentes phases, par exemple des minéraux et des espèces aqueuses. Dans [START_REF] Angeli | A tutorial on chemical reaction network dynamics[END_REF], l'auteur explique la diversité des systèmes réactifs et propose une représentation de ceux-ci. La vitesse des réactions cinétiques est modélisée par la loi d'action de masse [START_REF] Chellaboina | Modeling and analysis of mass-action kinetics[END_REF] et [START_REF] Lund | Guldberg and waage and the law of mass action[END_REF]. Cette loi dans le cas de réaction purement aqueuse assure à elle seule la positivité des concentrations des espèces.

Dans cette thèse, on se focalise sur les difficultés engendrées par la prise en compte des réactions de cinétiques chimiques contenant des minéraux. Une des principales difficultés est de gérer l'apparition et la disparition des phases en maintenant la positivité des concentrations des espèces. Dans le cas des équilibres chimiques, un travail récent [START_REF] Erhel | Characterizations of solutions in geochemistry: existence, uniqueness, and precipitation diagram[END_REF] a notamment permis de caractériser les contextes de présence des minéraux et d'assurer les propriétés mathématiques du modèle proposé. En cinétique, des modèles de nature différente ont été développés pour garantir la positivité des concentrations lors d'un processus de dissolution ou de précipitation d'un minéral. A notre connaissance, ces formulations ne peuvent prendre en considération plus d'un minéral par réaction.

Dans [START_REF] Bouillard | Diffusion with dissolution and precipitation in a porous medium: Mathematical analysis and numerical approximation of a simplified model[END_REF], [START_REF] Hoffmann | Reactive transport and mineral dissolution/precipitation in porous media : efficient solution algorithms, benchmark computations and existence of global solutions[END_REF], [START_REF] Hoffmann | Existence and uniqueness of a global solution for reactive transport with mineral precipitation-dissolution and aquatic reactions in porous media[END_REF], [START_REF] Kräutle | Existence of global solutions of multicomponent reactive transport problems with mass action kinetics in porous media[END_REF], [START_REF] Van Duijn | An analysis of crystal dissolution fronts in flows through porous media part 2: incompatible boundary conditions[END_REF] et [START_REF] Van Duijn | Crystal dissolution and precipitation in porous media: pore scale analysis[END_REF], les auteurs proposent un modèle d'EDOs discontinues pour modéliser les réactions cinétiques. Ces modèles utilisent une vitesse de réaction modifiée de telle sorte que la réaction ne peut plus consommer le minéral lorsqu'il disparaît. La nouvelle vitesse de réaction est obtenue en multipliant le terme qui consomme le minéral par un indicateur de présence du minéral. Les auteurs démontrent l'existence, la positivité et l'unicité des solution par l'utilisation de résultats classiques d'analyse (Lemme de Grönwall et Théorème du point fixe).

Une autre alternative est fournie dans [START_REF] Hoffmann | Reactive transport and mineral dissolution/precipitation in porous media : efficient solution algorithms, benchmark computations and existence of global solutions[END_REF] et [START_REF] Hoffmann | Existence and uniqueness of a global solution for reactive transport with mineral precipitation-dissolution and aquatic reactions in porous media[END_REF] : les auteurs proposent un modèle de type complémentarité pour gérer les apparitions et disparitions des espèces. Une nouveauté qu'apporte cette formulation est la gestion des contextes de présence des espèces directement dans le système d'équations. Cela permet l'utilisation de méthodes numériques plus efficaces [START_REF] Kanzow | New ncp-functions and their properties[END_REF] et [START_REF] Qi | A survey of some nonsmooth equations and smoothing Newton methods[END_REF] que celles basées sur la gestion de chaque contexte de présence.

Dans [START_REF] Agosti | Analysis of a model for precipitation and dissolution coupled with a darcy flux[END_REF] et [START_REF] Hoffmann | Reactive transport and mineral dissolution/precipitation in porous media : efficient solution algorithms, benchmark computations and existence of global solutions[END_REF], les auteurs proposent une formulation de type inclusion différentielle pour les cinétiques chimiques. Dans [START_REF] Hoffmann | Reactive transport and mineral dissolution/precipitation in porous media : efficient solution algorithms, benchmark computations and existence of global solutions[END_REF], les auteurs démontrent l'équivalence des trois types de formulation (EDOs discontinue, complémentarité et inclusion différentielle) au sens faible. Dans [START_REF] Agosti | Analysis of a model for precipitation and dissolution coupled with a darcy flux[END_REF], les auteurs utilisent des résultats classiques de la théorie des inclusions différentielles pour prouver l'existence et l'unicité de la solution de son modèle, ils utilisent une généralisation du théorème de Cauchy-Lipschitz qui nécessite un second membre one-sided Lipschitz.

Une importante contribution concernant l'analyse des EDOs discontinues a été fournie par Filippov [START_REF]Differential Equations with Discontinuous Righthand Sides[END_REF], Utkin [START_REF] Utkin | Sliding modes in control and optimization[END_REF], Aubin et Cellina [START_REF] Aubin | Differential Inclusions[END_REF]. Dans [START_REF]Differential Equations with Discontinuous Righthand Sides[END_REF] et [START_REF] Filippov | Differential equations whose right-hand is discontinuous on intersecting surfaces[END_REF], l'auteur reformule les problèmes d'EDOs discontinues comme des problèmes de contrôle ou comme des problèmes d'inclusion différentielle. En particulier, il prouve des résultats d'existence et d'unicité généraux. Dans le cas particulier où la surface de discontinuité est de codimension 1, l'auteur prouve un résultat d'unicité ne nécessitant pas de montrer la propriété de one-sided Lipschitz pour le second membre. Dans le cas des codimensions supérieures, le résultat proposé nécessite la non-oscillation des solutions et la non-tangence à la surface de discontinuité du champ de vecteur discontinu. Ces propriétés peuvent être complexes à obtenir en pratique. Un complément à cette théorie a été réalisé en étudiant les trajectoires possibles des solutions suivant la codimension et suivant la configuration du champ de vecteur discontinu [START_REF] Alexander | Sliding modes in intersecting switching surfaces[END_REF], [START_REF]Sliding modes in intersecting switching surfaces, i: Blending[END_REF], [START_REF] Dieci | A Filippov sliding vector field on an attracting co-dimension 2 discontinuity surface, and a limited loss-of-attractivity analysis[END_REF], [START_REF] Dieci | Regularizing piecewise smooth differential systems: Co-dimension 2 discontinuity surface[END_REF], [START_REF] Dieci | Sliding motion in Filippov differential systems: Theoretical results and a computational approach[END_REF], [START_REF]Sliding motion on discontinuity surfaces of high co-dimension. a construction for selecting a Filippov vector field[END_REF], [START_REF]A survey of numerical methods for ivps of odes with discontinuous right-hand side[END_REF] et [START_REF] Guglielmi | Classification of hidden dynamics in discontinuous dynamical systems[END_REF]. En particulier, cela permet de définir une valeur du second membre discontinu sur les discontinuités.

Il existe des liens d'équivalence entre les systèmes de complémentarité, les inclusions différentielles et les systèmes de dynamique projetée, voir [START_REF] Brogliato | On the equivalence between complementarity systems, projected systems and differential inclusions[END_REF]. Les systèmes de dynamique projetée sont une autre manière de modéliser les EDOs discontinues qui, à notre connaissance, n'a pas été explorée dans le cadre de la cinétique. Cette théorie fournit des résultats d'existence et d'unicité de la solution [START_REF] Cojocaru | Projected dynamical systems and evolutionary variational inequalities via Hilbert spaces with applications[END_REF], [START_REF] Cojocaru | Infinite-dimensional projected dynamics and the 1-dimensional obstacle problem[END_REF], [START_REF] Cojocaru | Existence of solutions to projected differential equations in hilbert spaces[END_REF] et [START_REF] Nagurney | Projected dynamical systems and variational inequalities with applications[END_REF] et une large gamme de méthodes numériques [START_REF] Acary | Numerical Methods for Nonsmooth Dynamical Systems: Applications in Mechanics and Electronics[END_REF], [START_REF] Bauschke | On the convergence of von neumann's alternating projection algorithm for two sets[END_REF], [START_REF] Bauschke | On the method of cyclic projections for convex sets in hilbert space[END_REF] et [START_REF] Mangasarian | Nonlinear programming[END_REF]. Dans [START_REF] Bauschke | On the convergence of von neumann's alternating projection algorithm for two sets[END_REF] et [START_REF] Bauschke | On the method of cyclic projections for convex sets in hilbert space[END_REF], les auteurs proposent un algorithme qui permet de résoudre numériquement un système de dynamique projetée sur une intersection d'ensembles convexes.

Descriptions des travaux de thèse

En géochimie, on s'intéresse à l'étude des interactions chimiques fluide-roche dans le sous-sol, qui ont des échelles de temps très diverses. Les cinétiques chimiques permettent de modéliser les réactions chimiques lentes. Contrairement aux réactions d'équilibre, les réactions de cinétique nécessitent la connaissance de la loi de vitesse de réaction qui lui est associée, là où les équilibres ne demandent que la connaissance de l'état final de la réaction appelé équilibre chimique. Cette thèse traite de la modélisation mathématique et de la simulation numérique des cinétiques chimiques.

Les modèles introduits dans cette thèse pour les systèmes de cinétique chimique prennent en compte les apparitions et disparitions de phase lors de processus chimiques. Ils étendent les résultats de [START_REF] Bouillard | Diffusion with dissolution and precipitation in a porous medium: Mathematical analysis and numerical approximation of a simplified model[END_REF] et [START_REF] Hoffmann | Existence and uniqueness of a global solution for reactive transport with mineral precipitation-dissolution and aquatic reactions in porous media[END_REF] aux systèmes réactifs contenant plusieurs espèces en phase pure (minérale ou gazeuse) par réaction. Ces modèles permettent également de décrire des systèmes réactifs faisant intervenir des intermédiaires réactifs. La modélisation des systèmes de cinétique complexe, [START_REF] Angeli | A tutorial on chemical reaction network dynamics[END_REF] et [START_REF] Friedly | Solute transport with multiple equilibrium-controlled or kinetically-controlled chemical reactions[END_REF], est prise en charge par les modèles. De plus, les modèles proposés n'annulent pas nécessairement les vitesses des réactions chimiques lorsqu'une espèce disparaît.

On considère c pc a , c p q la concentration des espèces en phase aqueuse et en phase pure, le système d'équations relatif à la cinétique chimique est défini par dc dt Sr, cp0q c 0 ¥ 0, [START_REF] Deutsch | Best approximation in inner product spaces[END_REF] où S est la matrice de stoechiométrie associée aux réactions de cinétiques, c 0 la condition initiale et r représente le nouveau taux de réaction défini plus loin, voir [START_REF] Dieci | Sliding motion in Filippov differential systems: Theoretical results and a computational approach[END_REF].

Soit Q une base du noyau de S T , si c est une solution de [START_REF] Deutsch | Best approximation in inner product spaces[END_REF] alors

Q T c Q T c 0 . (29)
Ce sont les lois de conservation intrinsèques au système chimique. La loi de cinétique considérée est la loi d'action de masse cinétique définie pour une réaction j, quelconque, comme suit,

τ j : k f,j ¹ i : S i,j 0 rE i s K i ¨|S i,j | ¡ k b,j ¹ i : S i,j ¡0 rE i s K i ¨|S i,j | , (30)
i : k r,j k f,j ± i:S i,j 0 K ¡|S i,j | i et k p,j k b,j ± i:S i,j ¡0 K ¡|S i,j | i .
L'activité d'une espèce aqueuse est une fonction de la fraction molaire des espèces (on peut la réécrire en fonction de la concentration ou de la molalité). L'activité d'une espèce minérale en phase pure est égale à 1. Par conséquent, la loi de cinétique ne dépend pas des concentrations des phases pures. Comme montré dans [START_REF] Bouillard | Diffusion with dissolution and precipitation in a porous medium: Mathematical analysis and numerical approximation of a simplified model[END_REF], sans modification de la loi de vitesse de réaction, une espèce chimique en phase pure peut être amenée à prendre des valeurs de concentrations négatives.

Dans le premier chapitre, nous proposons un modèle discontinu de cinétique chimique permettant de traiter les apparitions et disparitions d'un nombre quelconque d'espèces minérales. Nous proposons une analyse détaillée d'un cas contenant une espèce réactive intermédiaire. Nous introduisons un modèle régularisé permettant de prouver un résultat d'existence et de positivité des solutions. Ce même modèle permet aussi de proposer une méthode numérique de résolution pour les cinétiques chimiques. Enfin, nous validons cette méthode numérique sur des exemples de systèmes réactifs.

Dans le deuxième chapitre, nous proposons une reformulation du modèle de cinétique chimique proposé dans le premier chapitre à l'aide de la théorie de Filippov. Nous réalisons une étude de l'existence et de la positivité des solutions de la formulation inclusion différentielle ainsi obtenue. Nous menons une étude de l'unicité et des trajectoires de la solution dans le cas des surfaces de discontinuité de codimension 1.

Dans le troisième chapitre, nous proposons un nouveau modèle de cinétique chimique de type système dynamique projetée. Nous réalisons une analyse de l'existence des solutions de ce modèle, des liens avec d'autres types de formulations et une méthode de résolution numérique adaptée. Nous illustrons les résultats numériques obtenus pour des systèmes de cinétique chimique.

Chapitre 1 : Modèles de cinétique chimique avec des phases pures

Ce travail fait l'objet d'un article en cours de rédaction en vue d'une soumission dans un journal. Les auteurs sont : Bastien Hamlat, Jocelyne Erhel, Anthony Michel et Thibault Faney.

Afin d'assurer la positivité de la concentration des espèces, on introduit une surface réactive dans la loi de cinétique. En particulier, une espèce disparue peut continuer à être consommée par le système réactif au plus à son taux de production. Par conséquent, la vitesse de réaction est définie par

r j : h j τ j , (31)
avec h j r0, 1s un facteur permettant de limiter la vitesse de réaction lorsque la surface réactive est constante ou d'ordre 0. Nous avons étendu le concept de surface réactive d'ordre 0 développé dans [START_REF] Bouillard | Diffusion with dissolution and precipitation in a porous medium: Mathematical analysis and numerical approximation of a simplified model[END_REF] et [START_REF] Hoffmann | Existence and uniqueness of a global solution for reactive transport with mineral precipitation-dissolution and aquatic reactions in porous media[END_REF] à des systèmes chimiques généraux comprenant des phases pures minérales et gazeuses ainsi qu'une phase aqueuse. Précisons notre choix de modélisation pour h : Définition 0.3.1. Les fonctions multivaluées U j , j 1, . . . , N r sont définies dans r0, Vq Np par :

U j pc p q t1u si dk M ¡ j , c k ¡ 0, ou M ¡ j r, r0, 1s sinon.
Les fonctions multivaluées V j , j 1, . . . , N r sont définies dans r0, Vq Np par :

V j pc p q t1u si dk M j , c k ¡ 0, ou M j r, r0, 1s sinon.

Le taux de réaction modifié est alors défini comme suit :

Définition 0.3.2. Soit T ¡ 0 et c pH 1 p0, T qq Na Np tel que cptq ¥ 0, 0 ¤ t ¤ T.
Le taux de réaction modifié r est défini comme suit :

r j u j τ j pcq ¡ v j τ ¡ j pcq, j 1, . . . , N r , (32)
où u j , v j L V p0, T q et u j ptq U j pc p ptqq, v j ptq V j pc p ptqq, 0 t T.

Nous obtenons ainsi le modèle d'inclusion différentielle suivant :

trouver c pH 1 p0, T qq Na Np , u pL V p0, T qq Nr , v pL V p0, T qq Nr tel que dcptq dt Srpt, cptqq, t p0, T q, cptq ¥ 0, cp0q c 0 , u j ptq U j pc p ptqq, v j ptq V j pc p ptqq, t p0, T q, r j u j τ j pc a q ¡ v j τ ¡ j pc a q, j 1, . . . , N r ,

où c 0 ¥ 0 est la condition initiale, U j et V j sont introduits en Définition 0.3.1.

On peut noter que les lois de conservation (29) sont respectées par ce modèle. Les objectifs sont alors de caractériser la surface réactive u j et v j , de prouver l'existence d'une solution et de construire un algorithme de résolution.

Une première contribution du chapitre est la construction d'un modèle régularisé. Nous appliquons une fonction de Heaviside régularisée H à chaque phase pure, puis une fonction mélange aux phases pures impliquées dans une réaction. Nous définissons ainsi des surfaces réactives régularisées u et v , ainsi qu'une vitesse de réaction régularisée r : Définition 0.3.3. Le taux de réaction régularisé est défini comme suit :

r j, pcq u j, pcqτ j pcq ¡ v j, pcqτ ¡ j pcq, (34)
où u j, and v j, sont définis par si M ¡ j r, alors u j, pcq 1, sinon u j, pcq min

kM ¡ j pH pc k qq, si M j r, alors v j, pcq 1, sinon v j, pcq min kM j pH pc k qq. ( 35
)
dt Sr pcptqq, t ¥ 0, cp0q c 0 .

(36)

Ce système [START_REF] Erhel | Analysis of a global reactive transport model and results for the MoMas benchmark[END_REF] admet une unique solution c . Théorème 0.3.1. On suppose qu'il existe un vecteur q tel que q ¡ 0 et S T q 0. le système [START_REF] Erhel | Analysis of a global reactive transport model and results for the MoMas benchmark[END_REF] admet une unique solution c dans r0, Vr qui satisfait

0 ¤ c ¤ q T c 0 min k pq k q .
Une deuxième contribution du chapitre est un résultat d'existence des solutions faibles de [START_REF]A survey of numerical methods for ivps of odes with discontinuous right-hand side[END_REF] grâce à une analyse de convergence de [START_REF] Erhel | Analysis of a global reactive transport model and results for the MoMas benchmark[END_REF].

Théorème 0.3.2. On suppose qu'il existe un vecteur q tel que q ¡ 0 et S T q 0. Soit T ¡ 0. La suite c , solution de [START_REF] Erhel | Analysis of a global reactive transport model and results for the MoMas benchmark[END_REF], converge (à une sous-suite près) fortement dans C 0 p0, T, R pNa Npq q. De plus, toute limite c vérifie 0 ¤ c ¤ q T c 0 min k pq k q . Pour tout j 1, . . . , N r , les contrôles régularisés u j, pc q et v j, pc q admettent des limites au sens faible-* dans L V p0, T, R Nr q notés u j et v j respectivement. De plus, pour tout j 1, . . . , N r et 0 t T , les limites c, u, v satisfont u j ptq U j pc p ptqq et v j ptq V j pc p ptqq. Enfin c, u, v sont des solutions faibles de [START_REF]A survey of numerical methods for ivps of odes with discontinuous right-hand side[END_REF].

Une troisième contribution du chapitre est l'analyse d'un cas de cinétique chimique :

M 1 é M 2 et M 2 é M 3 .
Dans ce contexte, le système discontinu est défini comme suit :

dc 2 dt r 1 ¡ r 2 , dc 3 dt r 2 , cp0q c 0 , c ¥ 0, (37) avec r 1 u 1 τ 1 ¡ v 1 τ ¡ 1 et r 2 u 2 τ 2 ¡ v 2 τ ¡
2 . Après une analyse des différents contextes possibles, on obtient la formulation du modèle discontinu limite. Le taux de réaction discontinu dans ce cas est défini comme suit :

4 r 1 Hpc 1 qτ 1 ¡ Hpc 2 qτ ¡ 1 ¡ δpc 2 qHpc 3 qτ ¡ 2 , r 2 Hpc 2 qτ 2 δpc 2 qHpc 1 qτ 1 ¡ Hpc 3 qτ ¡ 2 . ( 38
) avec δ : x Ñ 5 0 si x R ¦ 1 si x 0 et H la fonction de Heaviside.
La dernière contribution de ce chapitre est la proposition d'un schéma numérique basé sur la formulation régularisée [START_REF] Erhel | Analysis of a global reactive transport model and results for the MoMas benchmark[END_REF], qu'on illustre sur plusieurs cas chimiques :

A é M , M 1 é M 2 é M 3 , A 1 é M 2 é A 3 et d'autres cas plus complexes.

Chapitre 2 : Analyse d'un modèle de type EDOs discontinues pour la cinétique chimique multiphasique

Dans ce chapitre, nous utilisons le cadre des équations différentielles discontinues de la théorie de Filippov [START_REF]Differential Equations with Discontinuous Righthand Sides[END_REF] pour analyser le système [START_REF]A survey of numerical methods for ivps of odes with discontinuous right-hand side[END_REF].

Soit i l'indice d'un minéral, on note Σ i tc i 0u une surface de discontinuité.

Ces surfaces séparent l'espace des phases pures en des domaines G k qui correspondent à des portions de l'espace avec un signe constant pour les concentrations des minéraux.

Chaque domaine G k est caractérisé par un vecteur e k : pour tout i t1, ¤ ¤ ¤ , N p u,

e k i 5 1, si c i ¡ 0, 0, si c i 0. On note G l'union disjointe des G k et Σ l'union des Σ i .
Pour construire le modèle EDOs discontinues, l'idée est de prolonger le modèle [START_REF]A survey of numerical methods for ivps of odes with discontinuous right-hand side[END_REF] dans les domaines G k où au moins une concentration des minéraux est négative. Dans notre modèle, les limiteurs de réaction valent 0 dans le domaine G k si au moins un des minéraux à concentration négative est consommé par la réaction j. Le premier résultat de ce chapitre assure que toute fonction absolument continue solution de [START_REF] Haeberlein | Time Space Domain Decomposition Methods for Reactive Transport, Application to CO2 Geological Storage[END_REF] est à valeur positive ou nulle.

On définit cette extension comme suit, r

k : G k Þ Ñ R Nr tel que r k j pcq ¹ iM ¡ j e k i ¨τ j pcq ¡ ¹ iM j e k i ¨τ ¡ j pcq, dk t1, ¤ ¤ ¤ , 2 Np u, dj t1, ¤ ¤ ¤ , N r u.
Théorème 0.3.3. Soit c : r0, T s Þ ÝÑ R Na Np une solution de (42). Alors dt r0, T s, cptq C 0 tc ¦ pR q Na Np , Q T c ¦ Q T c 0 u.
De plus, grâce à la théorie de Filippov [START_REF]Differential Equations with Discontinuous Righthand Sides[END_REF], nous montrons l'existence d'une solution de [START_REF] Haeberlein | Time Space Domain Decomposition Methods for Reactive Transport, Application to CO2 Geological Storage[END_REF]. Théorème 0.3.4. Sous l'hypothèse qu'il existe q ¡ 0 tel que q KerpS T q, pour tout c 0 pR q Na Np , il existe une solution c : r0, VrÑ C 0 de [START_REF] Haeberlein | Time Space Domain Decomposition Methods for Reactive Transport, Application to CO2 Geological Storage[END_REF].

Le deuxième résultat concerne l'unicité de la solution de [START_REF] Haeberlein | Time Space Domain Decomposition Methods for Reactive Transport, Application to CO2 Geological Storage[END_REF] ainsi que l'étude des trajectoires dans le cas où un seul minéral disparaît. Pour ce faire, l'idée est d'effectuer une classification et une étude de toutes les configurations possibles de champs signés afin d'en déduire les trajectoires des solutions de [START_REF] Haeberlein | Time Space Domain Decomposition Methods for Reactive Transport, Application to CO2 Geological Storage[END_REF]. En particulier, l'étude montre que le cas répulsif n'est pas possible en cinétique chimique.

Dans le contexte de l'étude des discontinuités de codimension 1 on définit

G 1 tpc a , c p q, c a ¥ 0, c p ¡ 0u et G 2 tpc a , c p q, c a ¥ 0, c p,1 0, c p,i ¡ 0, i ¥ 2u les deux ensembles autour de la surface de discontinuité Σ 1 tpc a , c p q, c a ¥ 0, c p,1 0, c p,i ¡ 0, i ¥ 2u. De la même manière, f 1 et f 2 sont les second-membres discontinus de (42) dans les domaines G 1 et G 2 .
Le second-membre multivalué F de (42) est défini comme suit

F pcq 6 9 8 9 7 tf 1 pcqu, c G 1 , tf 2 pcqu, c G 2 , tλf 1 pcq p1 ¡ λqf 1 pcq, λ r0, 1su, c Σ 1 . (43)
Notre étude des configurations des champs de vecteur de cinétique chimique nous conduit aux trois cas possibles décrits ci-dessous :

Soit c Σ 1 , -Attractif : rf 1 1 pcq 0 pet f 2 1 pcq ¥ 0qs ou rf 1 1 pcq 0 et f 2 1 pcq ¡ 0s.
-Transverse :

f 1 1 pcq ¡ 0 (et f 2 1 pcq ¡ 0).
-Tangent :

f 1 1 pcq 0 et f 2 1 pcq 0.
Dans les cas attractif et transverse on prouve l'unicité de la solution de (42) en s'appuyant sur des résultats de la théorie de Filippov dans le contexte des surfaces de discontinuité de codimension 1.

Chapitre 3 : Modèle avec dynamique projetée pour les réactions de dissolution et de précipitation des minéraux en géochimie

Ce travail fait l'objet d'un article en cours de rédaction en vue d'une soumission dans un journal. Les auteurs sont : Jocelyne Erhel, Bastien Hamlat et Tangi Migot.

Dans ce travail, réalisé en collaboration avec Jocelyne Erhel et Tangi Migot, nous proposons un autre modèle de cinétique chimique, une analyse d'existence des solutions et une méthode numérique de résolution. La première contribution de ce chapitre est de proposer un modèle pour lequel la dynamique de chaque espèce est limitée par un facteur λ i r0, 1s de sorte que dc i dt λ i pSτq i . On note Dpcq l'ensemble qui vérifie cette contrainte, défini ci-dessous,

Dpcq : tdiagpSτpcqqλ with 0 ¤ λ ¤ 1u.
De plus, les concentrations des espèces chimiques doivent respecter les lois de conservation et la positivité. On note C 0 l'ensemble associé à ces contraintes, défini comme suit,

C 0 : tc ¥ 0 : Q T pc ¡ c 0 q 0u.
Cette approche aboutit à un système projectif de façon à garantir la positivité de c, à respecter les lois de conservation et à maintenir la vitesse des concentrations proportionnelle à la vitesse non limitée. Le modèle proposé et étudié dans ce chapitre est le suivant

6 9 8 9 7 9 cptq proj pDpcptqq T pC 0 , cptqq; Sτ pcptqqq , cptq C 0 , cp0q c 0 , (44)
M 1 é M 2 et M 2 é M 3 et dynamique projeté sans la contrainte 9 c Dpcq dans le cas M 1 é M 2 et M 2 é M 3 .
La troisième contribution de ce chapitre est la démonstration de l'existence des solutions de [START_REF] Hoffmann | Reactive transport and mineral dissolution/precipitation in porous media : efficient solution algorithms, benchmark computations and existence of global solutions[END_REF]. Pour cela, nous proposons une étude d'un schéma discrétisé de [START_REF] Hoffmann | Reactive transport and mineral dissolution/precipitation in porous media : efficient solution algorithms, benchmark computations and existence of global solutions[END_REF] ainsi qu'une analyse de convergence. Nous utilisons un schéma explicite en temps transformant [START_REF] Hoffmann | Reactive transport and mineral dissolution/precipitation in porous media : efficient solution algorithms, benchmark computations and existence of global solutions[END_REF] en un problème non linéaire à chaque pas de temps. Soit h t f N h où N h est le nombre de pas de temps et

t f le temps final. Pour tout k 0, ¤ ¤ ¤ , N h ¡ 1, nous définissons c k 1 par 5 c k 1 K k , c k 1 projpK k ; c k hSτ pc k qq, (45) où K k C 0 pc k hDpc k qq. Il est à noter que grâce aux propriétés de K k , c k 1 existe et est unique.
Nous prouvons l'existence d'une solution de [START_REF] Hoffmann | Reactive transport and mineral dissolution/precipitation in porous media : efficient solution algorithms, benchmark computations and existence of global solutions[END_REF] en construisant une fonction c h , définie pour t r0,

t f s et k t h par c h ptq c k h pc k 1 h ¡ c k h q t ¡ kh h . ( 46
)
On a le résultat suivant Théorème 0.3.5. c h converge uniformément vers c solution de [START_REF] Hoffmann | Reactive transport and mineral dissolution/precipitation in porous media : efficient solution algorithms, benchmark computations and existence of global solutions[END_REF].

La preuve repose sur des résultats d'analyse classique comme le théorème d'Arzela-Ascoli et le lemme de Mazur. La dernière contribution de ce chapitre concerne la simulation numérique. Nous utilisons, l'algorithme itératif avec double projection pour calculer c k 1 pour une tolérance ¡ 0 :

c k 1 0 c k hSτ pc k q, e 0 e ¡1 0, j 0,
Répéter en incrémentant j :

c k 1 2j 1 projpC 0 ; c k 1 2j e 2j¡1 q, e 2j 1 c k 1 2j e 2j¡1 ¡ c k 1 2j 1 , c k 1 2j 2 c k hproj £ Dpc k q; c k 1 2j 1 e 2j ¡ c k h , e 2j 2 c k 1 2j 1 e 2j ¡ c k 1 2j 2 , jusqu'à }c k 1 2j ¡ c k 1 2j¡1 } ¤ (47) alors c k 1 c k 1 2j .
Cet algorithme est programmé en Matlab afin de résoudre ( 44) numériquement. Une validation de la méthode numérique est effectuée sur des cas tests.

Chapter 1 Kinetic models with pure phases

This work corresponds to a preprint in view of a submission in a journal. The authors are: Bastien Hamlat, Jocelyne Erhel, Anthony Michel and Thibault Faney.

Introduction

In this paper, we are interested in reactive transport modeling. This kind of models are a coupled ODEs-PDEs, where the PDEs describe the fluid transport and the ODEs describe the kinetic reaction, for reactive transport models including minerals see [START_REF] Knabner | An analysis of crystal dissolution fronts in flows through porous media. Part 1: Compatible boundary conditions[END_REF], [START_REF] Rubin | Transport of reacting solutes in porous media: Relation between mathematical nature of problem formulation and chemical nature of reactions[END_REF], [START_REF] Van Duijn | Travelling wave behaviour of crystal dissolution in porous media flow[END_REF], [START_REF] Van Duijn | An analysis of crystal dissolution fronts in flows through porous media part 2: incompatible boundary conditions[END_REF] and [START_REF] Van Duijn | Crystal dissolution and precipitation in porous media: pore scale analysis[END_REF]. Algebraic constraints are also included in order to take into account equilibrium reactions in the model. In the present work, we focus on ODEs models for chemical kinetic systems.

In chemical kinetics, the speed of a reaction is described by the reaction rate. The reaction rate is modeled by the law of mass action [START_REF] Chellaboina | Modeling and analysis of mass-action kinetics[END_REF], a nonlinear kinetic law, see Section 1.2.4. As proved in several papers, [START_REF] Bouillard | Diffusion with dissolution and precipitation in a porous medium: Mathematical analysis and numerical approximation of a simplified model[END_REF] and [START_REF] Hoffmann | Existence and uniqueness of a global solution for reactive transport with mineral precipitation-dissolution and aquatic reactions in porous media[END_REF], this law is not mathematically adapted to take into account the disappearance of a mineral. In order to handle this, we consider kinetic-controlled models [START_REF] Agosti | Analysis of a model for precipitation and dissolution coupled with a darcy flux[END_REF], [START_REF] Bouillard | Diffusion with dissolution and precipitation in a porous medium: Mathematical analysis and numerical approximation of a simplified model[END_REF], [START_REF] Friedly | Solute transport with multiple equilibrium-controlled or kinetically-controlled chemical reactions[END_REF], [START_REF] Hoffmann | Existence and uniqueness of a global solution for reactive transport with mineral precipitation-dissolution and aquatic reactions in porous media[END_REF], [START_REF] Knabner | An analysis of crystal dissolution fronts in flows through porous media. Part 1: Compatible boundary conditions[END_REF], [START_REF] Kräutle | Existence of global solutions of multicomponent reactive transport problems with mass action kinetics in porous media[END_REF], [START_REF] Kräutle | Robust simulation of mineral precipitation-dissolution problems with variable mineral surface area[END_REF] and [START_REF] Van Duijn | Crystal dissolution and precipitation in porous media: pore scale analysis[END_REF]. A first possibility is to consider a discontinuous reaction rate, see Section 1.3.1 for an introduction of this model and for more details [START_REF] Bouillard | Diffusion with dissolution and precipitation in a porous medium: Mathematical analysis and numerical approximation of a simplified model[END_REF] and [START_REF] Friedly | Solute transport with multiple equilibrium-controlled or kinetically-controlled chemical reactions[END_REF]. A second solution is to consider a set-va lued model, see Section 1.3.2 for an introduction and for more details [START_REF] Agosti | Analysis of a model for precipitation and dissolution coupled with a darcy flux[END_REF], [START_REF] Knabner | An analysis of crystal dissolution fronts in flows through porous media. Part 1: Compatible boundary conditions[END_REF] and [START_REF] Van Duijn | Crystal dissolution and precipitation in porous media: pore scale analysis[END_REF]. A third possibility is to consider complementarity conditions in order to change the value of the reaction rate depending on the presence of the mineral, see Section 1.3.3 for an introduction and [START_REF] Hoffmann | Existence and uniqueness of a global solution for reactive transport with mineral precipitation-dissolution and aquatic reactions in porous media[END_REF] for more details. All of these models allow to consider one mineral in dissolution or precipitation in a kinetic reaction.

The aim of this paper is to investigate the case of several minerals in a kinetic reaction or a mineral involved in several reactions. The physical motivation is to consider complex chemical systems as described in [START_REF] Angeli | A tutorial on chemical reaction network dynamics[END_REF]. The main difficulty is to propose a kinetic model which allows to maintain the reaction process when one of the minerals reactant disappear. In this paper, we generalize the set-valued formulation for kinetic problems containing several minerals in each reaction, see [START_REF] Friedly | Solute transport with multiple equilibrium-controlled or kinetically-controlled chemical reactions[END_REF] and [START_REF] Rubin | Transport of reacting solutes in porous media: Relation between mathematical nature of problem formulation and chemical nature of reactions[END_REF] for example of reaction containing several minerals. This generalization is obtained by introducing two limiter set-valued functions for each reaction (one for the reactant side and another for the product side). Each of this limiters can depend on several minerals, i.e. each reactant limiter depends on all the reactants in its reaction. We follow the idea of [START_REF] Van Duijn | Crystal dissolution and precipitation in porous media: pore scale analysis[END_REF] by proposing a regularized model to prove the existence of solutions of the set-valued model. Moreover, one novelty of this work is to use the regularized formulation for numerical experiments. Keeping in mind that the aim is to couple the kinetic model in a fully implicit reactive transport solver [START_REF] Kräutle | Robust simulation of mineral precipitation-dissolution problems with variable mineral surface area[END_REF] and [START_REF] Leal | An overview of computational methods for chemical equilibrium and kinetic calculations for geochemical and reactive transport modeling[END_REF].

This paper is divided into five main parts. The first one is a general introduction to chemical kinetics. The second one gathers the state of the art in the case of one mineral in one reaction. In the third part, we present our set-valued model, its regularization and associated convergence and existence results. In the fourth part, we investigate the case M 1 é M 2 and M 2 é M 3 . The last part concerns numerical experiments: simulation, convergence and an example of non-uniqueness.

Mathematical model of kinetic reactions

We consider a geochemical system with an aqueous phase and several mineral phases. The aqueous phase is composed of water and several diluted species whereas a mineral phase contains only one species. We do not consider source or sink terms. We assume that the pressure and temperature are constant during the reactions.

Kinetic reactions in this system are governed by reaction rates and conservation laws. Reactions are considered reversible, with forward and backward directions. The number of reactions is N r and each reaction is noted R j , j 1, . . . , N r . The number of aqueous species is N a ¡ 0 and each aqueous species is noted

A i , i 1, . . . , N a . The number of minerals is N p ¥ 0 and if N p ¡ 0, each mineral is noted M k , k N a 1, . . . , N a N p .
The stoichiometric matrix S Z pNa Npq¢Nr contains the coefficients of reactions.

We set an arbitrary forward direction of a reversible reaction R j . Species on the left are called reactants and are associated to negative stoichiometric coefficients, whereas species on the right, called products, correspond to positive coefficients [START_REF] Leal | An overview of computational methods for chemical equilibrium and kinetic calculations for geochemical and reactive transport modeling[END_REF]. We denote by A ¡ j and M ¡ j the sets of aqueous and mineral reactants in reaction R j , with negative stoichiometric coefficients s ij 0; by A j and M j the sets of aqueous and mineral products, with positive coefficients s ij ¡ 0.

A reaction R j is written

i A ¡ j |s ij |A i ķM ¡ j |s kj |M k Õ i A j |s ij |A i ķM j |s kj |M k .
(1.1)

Examples

We provide here two examples of such reactions in geology.

In the first example, calcite CaCO 3 reacts with calcium oxide CaO which reacts with calcium hydroxide CapOHq 2 :

R 1 : CaCO 3 Õ CaO CO 2 , R 2 : CaO H 2 0 Õ CapOHq 2 .
The second example still concerns calcite, which can react with dolomite CaM gpCO 3 q 2 and siderite F eCO 3 when the solution contains magnesium and irons:

R 1 : 2CaCO 3 M g Õ CaM gpCO 3 q 2 Ca , R 2 : 3CaCO 3 M g F e Õ CaM gpCO 3 q 2 F eCO 3 2Ca .
Finally, it can be noted that calcite is a form of calcium carbonate which can interact with another form, aragonite, under some conditions.

In the domain of electrochemistry, reactions between solute and minerals occur in batteries.

Quantities of species and reaction rates

The unknowns in such a geochemical system are the number of moles of each species, grouped into a vector of functions c pc i q, i 1, . . . , N a N p of the time variable t with t ¥ 0. We denote by c a pc i q, i 1, . . . , N a the vector of quantities for aqueous species, and by c p pc k q, k N a 1, . . . , N a N p the vector of mineral quantities. The quantities c should satisfy a positivity constraint cptq ¥ 0, t ¥ 0. We also require that c pH 1 p0, T qq Na Np with T ¡ 0.

In this paper, we propose a model of differential equations, which satisfy mass conservation laws and various invariants. The right-hand side uses reaction rates r pr j q, functions in L 2 p0, T q, based on mass action laws. Our objective is to define a model of reaction rates such that the differential system has a solution satisfying the requirements given above.

Conservation laws

Mass conservation laws are grouped into a system of differential equations: dcptq dt Srpt, cptqq, t p0, T q cptq ¥ 0, cp0q c 0 ,

(1.2)
where c 0 ¥ 0 is a given initial condition, and reaction rates r are discontinuous functions, defined thanks to mass action laws and set-valued functions. We define a weak solution of system (1.2).

Definition 1.2.1. Let T ¡ 0, let c H 1 p0, T q Na Np . It is a weak solution of system (1.2) if c ¥ 0 and dt r0, T s, c i ptq c 0,i » t 0 pSrq i ps, cpsqqds, i 1, . . . , N a N p .
It can be noted that this system implies invariants, which correspond to various balance equations. Indeed, let Q be a basis of the kernel of S T , such that Q T S S T Q 0. Then, any solution of (1.2) satisfies dt ¥ 0, Q T cptq Q T c 0 . Functions Q T c are thus constant in time and are conservative variables. For example, if some species are ions, then one conservative variable is the charge balance.

If a conservative variable q T c is defined by a positive vector q ¡ 0, then quantities of species are bounded: Lemma 1.2.1. Let q be a vector such that q ¡ 0. If c satisfies c ¥ 0 and q T c q T c 0 then c ¤ C with C 1 min k pq k q q T c 0 . Proof. Let k t1, . . . , N a N p u. First, q T c q k c k °i$k q i c i ¥ q k c k since c i ¥ 0 and

q i ¡ 0. Then q k c k ¤ q T c 0 and c k ¤ 1 q k q T c 0 ¤ C.
Physically, this result prohibits the chemical system from creating matter. In other words, the total mass of each atom is an invariant of the chemical system.

Maximal reaction rates

Reactions R j involve a forward reaction, with reactants and products, and a backward reaction, where products become reactants and vice versa. Each direction has its own reaction rate, and the global rate is the difference between both. We call it the maximal reaction rate.

When a reaction involves only aqueous species, its rate is given by a mass action law, which is a function of species' activities [START_REF] Chellaboina | Modeling and analysis of mass-action kinetics[END_REF]. When a pure phase exists, its activity is equal to 1; we extend it to 1 when the quantity of pure phase is zero. Various models of aqueous activities exist in the literature, see [START_REF] Carrayrou | Modélisation des phénomènes de transport réactif en milieu poreux saturé[END_REF] for example. Throughout the paper, we make the following assumptions. Assumption 1.2.1. The activity a i of a solute A i is a C 1 function of c a which satisfies:

c a ¥ 0 ñ apc a q ¥ 0 c a,i 0 ñ a i pc a q 0.
The forward reaction rate of reaction R j is given by

τ f,j k f,j ¹ iA ¡ j p a i K i q |s ij | ¹ kM ¡ j p 1 K k q |s kj |
and its backward reaction rate is

τ b,j k b,j ¹ iA j p a i K i q |s ij | ¹ kM j p 1 K k q |s kj | If A ¡ j is empty, the product ± iA ¡ j p a i K i q |s ij | is equal to 1.
The same property applies to other subsets. The maximal reaction rate is the difference between forward and backward reaction rates:

τ j τ f,j ¡ τ b,j .
Finally, we will use the following formula, where constants are grouped together:

τ j pc a q α j ¹ iA ¡ j a i pc a q |s ij | ¡ β j ¹ iA j a i pc a q |s ij | , j 1, . . . , N r . (1.3)
where α j and β j are strictly positive constants of reactions.

Thermodynamic equilibrium of reaction R j is obtained when forward and backward reactions are balanced, in other words when the reaction rate is zero: τ j 0. It corresponds to the mass action law for chemistry at equilibrium.

Continuous model

We consider the continuous differential system with the vector τ of maximal reaction rates:

5 dcptq dt Sτ pc a ptqq, t ¥ 0, cp0q c 0 .

(1.4) Proposition 1.2.1. The geochemical system (1.4) has a unique solution satisfying c pH 1 p0, Vqq Na Np .

Proof. Thanks to Assumption 1.2.1, the maximal rate τ is C 1 . Since the right-hand side is regular, we can apply theorem of Cauchy-Lipschitz so that the system has a unique solution in p0, Vq.

In order to study the non-negativity of quantities c i , we introduce the notion of essential non negativity [START_REF] Chellaboina | Modeling and analysis of mass-action kinetics[END_REF].

Definition 1.2.2. Let f be a function from R Na to R Na . It is essentially non-negative if, for any k 1, . . . , N a , c i ¥ 0, i 1, . . . , N a and c k 0 ñ f k pcq ¥ 0.
Since the reaction rates τ j depend only on the aqueous species, we can define a function f from R Na to R Na by f k pc a q pSτpc a qq k , k 1, . . . , N a .

(1.5) Proposition 1.2.2. The function f defined by (1.5) is locally Lipschitz continuous and essentially non negative.

Proof. Since τ is regular, f is at least locally Lipschitz continuous.

Let c a ¥ 0 with c k 0 for some k. Then apc a q ¥ 0 and a k pc k q 0. If s kj ¡ 0 then A k A j and τ j pc a q ¥ 0 thus s kj τ j ¥ 0. If s kj 0 then A k A ¡ j and τ j pc a q ¤ 0 thus s kj τ j ¥ 0. Finally, f k °j:

s kj $0 s kj τ j pc a q ¥ 0.
Thanks to this property, we can prove that the quantities of all aqueous species remain positive.

Proposition 1.2.3. If c 0 ¥ 0 then the unique solution c of system (1.4) satisfies c a ¥ 0.

Proof. Let us consider the system of ODE with only the first N a equations, in other words with the function f , which is essentially non negative and locally Lipschitz continuous by Proposition 1.2.2. By using Theorem 1 in [START_REF] Chellaboina | Modeling and analysis of mass-action kinetics[END_REF], we get the result wanted.

On the other hand, denoting by c the solution of (1.4) it may exist a mineral M k , N a 1 ¤ k ¤ N a N p and a time t ¦ such that c k pt ¦ q 0 and dt ¡ t ¦ , c k ptq 0. Models of reaction rates r aim at limiting the reaction so that the solution of (1.2) is non-negative. We first consider a specific case, well studied in the literature.

State of the art: one mineral in one reaction

Here, we consider the particular case of two aqueous species A 1 and A 2 which precipitate into a mineral M :

s 1 A 1 s 2 A 2 Õ M
with s 1 ¡ 0, s 2 ¡ 0. An example is the precipitation-dissolution of silver chloride: Ag pwq Cl ¡ pwq é AgCl psq . Such kinetic models have been studied by several authors [START_REF] Agosti | Analysis of a model for precipitation and dissolution coupled with a darcy flux[END_REF][START_REF] Bouillard | Diffusion with dissolution and precipitation in a porous medium: Mathematical analysis and numerical approximation of a simplified model[END_REF][START_REF] Hoffmann | Existence and uniqueness of a global solution for reactive transport with mineral precipitation-dissolution and aquatic reactions in porous media[END_REF].

Since the system of three equations has two invariants, it is sufficient to consider the equation for c M : dc M dt r.

The maximal reaction rate given by (1.3) is

τ αa s 1 1 a s 2
2 ¡ β, where a 1 , a 2 are activities of aqueous species A 1 , A 2 . It seems natural to choose r τ as long as c M is not fully dissolved (c M ¡ 0). Different choices of r appear in the literature when c M 0. With a discontinuous model, r is such that the reaction stops as soon as the mineral is fully dissolved [START_REF] Bouillard | Diffusion with dissolution and precipitation in a porous medium: Mathematical analysis and numerical approximation of a simplified model[END_REF]. With a set-valued model, r belongs to a given set [START_REF] Hoffmann | Existence and uniqueness of a global solution for reactive transport with mineral precipitation-dissolution and aquatic reactions in porous media[END_REF]; moreover, the value at c M 0 can be specified [START_REF] Agosti | Analysis of a model for precipitation and dissolution coupled with a darcy flux[END_REF]. A third model is based on a complementarity problem [START_REF] Hoffmann | Existence and uniqueness of a global solution for reactive transport with mineral precipitation-dissolution and aquatic reactions in porous media[END_REF]. Under some assumptions, the three models are equivalent in the sense that they have the same unique weak solution [START_REF] Hoffmann | Existence and uniqueness of a global solution for reactive transport with mineral precipitation-dissolution and aquatic reactions in porous media[END_REF]. Another model defines a projected dynamical system, which limits the derivative of each mineral species [START_REF] Erhel | A projected dynamical system approach to mineral precipitation-dissolution reactions in geochemistry[END_REF].

Discontinuous model

The discontinuous model implies r 0 after dissolution [START_REF] Bouillard | Diffusion with dissolution and precipitation in a porous medium: Mathematical analysis and numerical approximation of a simplified model[END_REF]:

4 if c M ¡ 0 or τ ¥ 0 then r τ, if c M 0 and τ ¤ 0 then r 0. (1.6)
Let us introduce the positive part τ pc a q and the negative part τ ¡ pc a q of the continuous reaction rate so that τ pc a q τ pc a q ¡ τ ¡ pc a q.

In order to ensure the positivity of c M , the reaction rate r is defined by multiplying the negative part of τ by a control function.

rpcq τ pc a q ¡ Hpc M qτ ¡ pc a q, (1.7)
where H is the Heaviside function defined in r0, Vr by

4 if x ¡ 0 then Hpxq 1, if x 0 then Hpxq 0.
Clearly, the choice (1.7) satisfies (1.6). It can be noted that the Heaviside function H thus the reaction rate r can be defined for negative values: Hpxq 0 if x 0. The kinetic model (1.2) is then a system of differential equations, where the reaction rate is a piecewise discontinuous function. It has a unique weak solution [START_REF] Bouillard | Diffusion with dissolution and precipitation in a porous medium: Mathematical analysis and numerical approximation of a simplified model[END_REF].

Set-valued model

The set-valued model does not specify the value of r with total dissolution. The idea is to limit the backward reaction rate β by a set-valued function:

r αa s 1 1 a s 2 2 ¡ w M β, where w M W pc M q with W puq t1u if u ¡ 0, r0, 1s if u 0, t0u if u 0.
(1.8) Thus, r τ as long as c M ¡ 0 and τ ¤ r ¤ αa s 1 1 a s 2 2 . When c M 0, the value of w M can be chosen such that r 0 [START_REF] Agosti | Analysis of a model for precipitation and dissolution coupled with a darcy flux[END_REF]. This model can be extended to a set of aqueous and mineral reactions, provided that a mineral reaction contains exactly one mineral which participates only in this reaction of precipitation and dissolution [START_REF] Hoffmann | Existence and uniqueness of a global solution for reactive transport with mineral precipitation-dissolution and aquatic reactions in porous media[END_REF]. An aqueous reaction R j is written

i A ¡ j |s ij |A i Õ i A j |s ij |A i , whereas a mineral reaction R j involving the mineral M k is written i A ¡ j |s ij |A i Õ i A j |s ij |A i M k .
Let R j be a reaction with mineral M k , with τ j defined by (1.3). Then the set-valued model for reaction rate r j is written

r j α j ¹ iA ¡ j a i pc a q |s ij | ¡ w k β j ¹ iA j a i pc a q |s ij | ,
where w k W pc k q with the set-valued function W defined by (1.8).

Under some assumptions, the reactive transport model coupling advection-diffusion with the geochemical system has a unique weak solution [START_REF] Hoffmann | Existence and uniqueness of a global solution for reactive transport with mineral precipitation-dissolution and aquatic reactions in porous media[END_REF]. The proof uses a regularized chemical system where the set-valued function W is replaced by a regularized Heaviside function.

Complementarity model

Here, the value of r when c M 0 is not given, since the function r is defined by the following complementarity problem:

c M pr ¡ τ q 0, c M ¥ 0, r ¥ τ,
which has a unique weak solution [START_REF] Hoffmann | Existence and uniqueness of a global solution for reactive transport with mineral precipitation-dissolution and aquatic reactions in porous media[END_REF].

Reference example

In order to illustrate the appearance and disappearance of a pure phase, we study a very simple model with one reaction involving one aqueous species and one mineral:

A é M . System (1.
2) is then written, with c pc A , c M q: 6 9 9 9 9 9 8 9 9 9 9 9 7

dc A ptq dt ¡rpcptqq, dc M ptq dt rpcptqq, c A ptq ¥ 0, c M ptq ¥ 0, c A p0q c A,0 , c M p0q c M,0 .
(1.9)

The conservative variable is pc A c M q, since Q p1, 1q T . Here, we choose a simple activity model with apc A q c A , so that the mass action law (1.3) becomes τ pc A q αc A ¡ β, which we rewrite (with α 1 and β K)

τ pc A q c A ¡ K,
where K is a strictly positive constant.

The continuous problem (1.4) with the reaction rate τ pc A q has a unique solution:

4 c ¦ A ptq K τ 0 expp¡tq, t ¥ 0, c ¦ M ptq c M,0 τ 0 p1 ¡ expp¡tqq, t ¥ 0.
(1.10) where τ 0 c A,0 ¡ K. Thus τ ptq τ 0 expp¡tq and c ¦ A ptq ¡ 0, t ¡ 0. If τ 0 0, then τ ptq 0 and rptq 0, thus the system is at equilibrium at any time. From now on, we assume that τ 0 $ 0.

Solution of the discontinuous model

The discontinuous reaction rate (1.7) becomes here:

rpc A , c M q τ pc A q ¡ Hpc M qτ ¡ pc A q.

Precipitation

If τ 0 ¡ 0, then dt ¡ 0, c ¦ A ptq ¡ 0, c ¦ M ptq ¡ 0, τ ptq ¡ 0. The solute precipitates but does not disappear and the quantity of mineral increases. Here, c ¦ is the weak solution of system (1.9) and formula (1.11) gives r τ.

1.3.4.3 Partial dissolution If ¡c M,0 ¤ τ 0 0, then dt ¡ 0, c ¦ A ptq ¡ 0, c ¦ M ptq ¡ 0.
The mineral dissolves but does not disappear and the quantity of solute increases. Here too, c ¦ is the weak solution of system (1.9) and r τ.

Full dissolution

If τ 0 ¡c M,0 , then c ¦ M pt M q 0 where t M ¥ 0 is given by

t M logp τ 0 τ 0 c M,0 q (1.12)
The mineral is fully dissolved. Let pc A , c M q be defined by:

6 9 9 8 9 9 7 c A ptq c ¦ A ptq, t ¤ t M , c M ptq c ¦ M ptq, t ¤ t M , c A ptq c A,0 c M,0 , t ¥ t M , c M ptq 0, t ¥ t M ,
It is the weak solution of (1.9) with r Hpc M qτ.

A differential inclusion model for general geochemical systems

In this section, we consider the general geochemical system described in section 1.2. We develop an original new model based on a differential inclusion [START_REF]Differential Equations with Discontinuous Righthand Sides[END_REF]. Another model defined in [START_REF] Erhel | A projected dynamical system approach to mineral precipitation-dissolution reactions in geochemistry[END_REF] is based on a projected dynamical system. In our model, the reaction rate r j is equal to τ j if the quantities of all reactant minerals in the reaction are strictly positive and it must be controlled if the quantity of one reactant mineral is zero. However, these controls may depend on the vector c p of all mineral quantities. As in the discontinuous model, we limit the positive and negative parts of τ j ; nevertheless, these controls belong to set-valued functions as in the set-valued model.

Our objective is to define reaction rates r in order to ensure a non-negative solution. We propose a model where reaction rates r j are governed by maximal reaction rates τ j in the following way: dj 1, . . . , N r , |r j | ¤ |τ j | and r j τ j ¥ 0.

(1.13) For any j t1, . . . , N r u, we introduce set-valued functions U j and V j of variables c p pR q Np . Then we introduce two control functions u j , v j L V p0, T q in order to define the reaction rates. Definition 1.4.1. The set-valued functions U j , j 1, . . . , N r are defined in r0, Vq Np by:

U j pc p q t1u if dk M ¡ j , c k ¡ 0, or M ¡ j r, r0, 1s else.
The set-valued functions V j , j 1, . . . , N r are defined in r0, Vq Np by: V j pc p q t1u if dk M j , c k ¡ 0, or M j r, r0, 1s else.

For general geochemical systems, we can now define our model of reaction rates.

Definition 1.4.2. Let T ¡ 0 and c pH 1 p0, T qq Na Np such that cptq ¥ 0, 0 ¤ t ¤ T.

Reaction rates r are defined by

r j u j τ j pcq ¡ v j τ ¡ j pcq, j 1, . . . , N r , (1.14)
where u j , v j L V p0, T q and u j ptq U j pc p ptqq, v j ptq V j pc p ptqq, 0 t T.

Clearly, the vector of reaction rates r satisfies (1.13). The reaction rate r j , j 1, . . . , N r is equal to the maximal value τ j if all the dissolving minerals are present. If, for example, τ j ¡ 0 and a reactant mineral is fully dissolved, then 0 ¤ r j ¤ τ j . The reaction is not stopped but it may continue at a limited rate.

The problem to solve is then a differential inclusion system [START_REF]Differential Equations with Discontinuous Righthand Sides[END_REF]: find c pH 1 p0, T qq Na Np , u pL V p0, T qq Nr , v pL V p0, T qq Nr such that dcptq dt Srpt, cptqq, t p0, T q, cptq ¥ 0, cp0q c 0 , u j ptq U j pc p ptqq, v j ptq V j pc p ptqq, t p0, T q, r j u j τ j pc a q ¡ v j τ ¡ j pc a q, j 1, . . . , N r ,

where c 0 ¥ 0 is a given initial condition, τ is defined by (1.3), U j and V j are introduced in Definition 1.4.1. In order to analyze and to solve this system, we introduce a regularization method, as in [START_REF] Hoffmann | Existence and uniqueness of a global solution for reactive transport with mineral precipitation-dissolution and aquatic reactions in porous media[END_REF]. A regularized model is controlled by a parameter , with a continuous reaction rate. We get a system of stiff ODEs and prove that it has a unique non negative solution. Then we prove that the differential inclusion (1.15) has a weak solution, limit of a convergent sequence of regularized functions. Moreover, the regularized problem can be solved by a classical ODE solver and provides a numerical method to approximate a solution of the discontinuous system.

Regularized model

In order to analyze the inclusion differential system, we introduce a regularized problem, as in [START_REF] Hoffmann | Existence and uniqueness of a global solution for reactive transport with mineral precipitation-dissolution and aquatic reactions in porous media[END_REF]. Let us first define a regularized Heaviside function H :

6 8 7 if x ¥ then H pxq 1, if 0 ¤ x ¤ then H pxq γpx{ q, if x 0 then H pxq 0, (1.16)
where γ is an increasing Lipschitz continuous function defined in the interval r0, 1s such that γp0q 0 and γp1q 1, so that the function H is Lipschitz continuous and bounded. For example, we can choose γpxq x.

We can now define regularized reaction rates. If only one mineral M k participates as a reactant in reaction R j then u j, pcq H pc k q and v j, pcq 1. If it is a product, then u j, pcq 1 and v j, pcq H pc k q. But if two minerals participate as reactants in the same reaction R j , we use a blending function to define the limiter, which takes values in the interval r0, 1s. This blending of the reactant minerals must be equal to 0 as soon as one of them disappears and be equal to 1 as long as all of them are present in sufficient quantity. Similar properties apply to product minerals. Examples of such blending are the minimum and the product functions. In what follows, we choose the min function.

Definition 1.4.3. For general geochemical systems, the regularized reaction rates are defined by r j, pcq u j, pcqτ j pcq ¡ v j, pcqτ ¡ j pcq,

(1.17

)
where u j, and v j, are defined by if M ¡ j r, then u j, pcq 1, else u j, pcq min kM ¡ j pH pc k qq, if M j r, then v j, pcq 1, else v j, pcq min kM j pH pc k qq.

(1.18)

It can be noted that these regularized reaction rates satisfy the requirements (1.13).

We study the regularized system of kinetic reactions:

5 dcptq dt Sr pcptqq, t ¥ 0, cp0q c 0 . (1.19)
It is easy to show that the regularized system (1.19) has a unique solution. hT ¡ 0 such that system (1.19) with non-negative initial condition c 0 , has a unique maximal non-negative solution c in r0, T r.

Proof. We apply the theorem of Cauchy-Lipschitz. Now, in order to prove the non negativity of the solution, we start by showing that the right-hand side is essentially non negative. First, let A k , 1 ¤ k ¤ N a be an aqueous species such that c k 0 for some k t1, . . . N a u. We proved in Proposition 1.2.2 that s kj τ j ¥ 0, so using 1.13, s kj r j, ¥ 0 and pSr q k ¥ 0.

Second, let M k , N a 1 ¤ k ¤ N a N p be a mineral such that c k 0. If M k is a reactant of reaction R j , then s kj 0; since c k 0 we get u j, 0 and r j, ¤ 0. On the contrary, if M k is a product of reaction R j , then s kj ¡ 0; since c k 0 we get v j, 0 and r j, ¥ 0. Finally pSr j, q k ¥ 0.

Convergence of the regularized solutions

In this section, we prove convergence of the regularized solutions, towards a nonnegative and bounded weak solution of the differential inclusion problem. From now on, we make the following assumption: Assumption 1.4.1. We assume the existence of a vector q such that q ¡ 0 and S T q 0. This assumption implies that the regularized solution is bounded and can be defined in r0, Vr: Proof. Let c be the maximal solution of (1.19). Since S T q 0, it satisfies q T c q T c 0 . Thanks to Theorem 1.4.2, it also satisfies c ¥ 0. Thus, using Lemma 1.2.1, it also satisfies c ¤ C. So, the maximal interval of the non-negative solution is r0, Vr.

It also implies that the sequence c converges:

Proposition 1.4.3. Let T ¡ 0. Under Assumption 1.4.1, the sequence c converges strongly in C 0 p0, T, R pNa Npq q. Any limit c of a subsequence satisfies 0 ¤ c ¤ C. Also for any j 1, . . . , N r , subsequences of the regularized controls u j, pc q and v j, pc q converge weakly-* in L V p0, T, R Nr q towards respectively u j and v j which satisfy 0 ¤ u j ¤ 1 and 0 ¤ v j ¤ 1. Moreover, for any j 1, . . . , N r and any 0 t T , limits c, u, v defined above satisfy u j ptq U j pc p ptqq and v j ptq V j pc p ptqq. Proof. Let c be the solution of (1.19) in r0, Vr. Since 0 ¤ c ¤ C, the norms }c } L 2 p0,Tq and }Sr pc q} L 2 p0,Tq are bounded, thus }c } H 1 p0,Tq is also bounded.

By compact injection of H 1 p0, T q in C 0 p0, T q, a subsequence of c converges in C 0 p0, T q and the limit c is bounded below by 0 and above by C.

Since 0 ¤ u j, pc q ¤ 1, a subsequence of u j, pc q converges weakly-* in L V p0, T q and the limit u j is bounded below by 0 and above by 1. The result is similar for v j . Thus dj 1, . . . , N r , the control satisfies u j L V p0, T q and 0 ¤ u j ptq ¤ 1, 0 t T. Now, we want to show that u j ptq U j pc p ptqq. Let 1 ¤ j ¤ N r and 0 t ¦ T such that M ¡ j $ r and dk M ¡ j , c k pt ¦ q ¡ 0. Since c is continuous, there exists an interval I containing t ¦ such that dt I, dk M ¡ j , c k ptq ¡ 0. Since a subsequence of c converges towards c in C 0 p0, T q, for sufficiently small, the same property is satisfied:

dt I, dk M ¡ j , c k, ptq ¡ δ k c k ptq 2 ¡ 0.
Consequently, dt I, u j, pc qptq 1, so that lim Ñ0 u j, pc q 1 in C 0 pIq, and dt I, u j ptq 1. This concludes the proof. Now, we can prove the existence of a weak solution of system (1.15):

Theorem 1.4.2. We make Assumption 1.4.1. Let T ¡ 0 and c a limit of c in C 0 p0, T, R pNa Npq q. Let u and v weakly-* limits of a subsequence of u pc q and v pc q in L V p0, T, R Nr q. Then c, u, v are a weak solution of system (1.15).

Proof. First, the limits u j and v j satisfy the requirements of definition 1.4.2. Also, cp0q c 0 and c ¥ 0. Since c weakly converges towards c in H 1 p0, T q, dφ pC V c p0, T qq, (1.20)

Thanks to Lemma 1.4.3, u j, , u j L V p0, T q and }u j, } V ¤ 1. Also, since τ j and c , c are continuous, τ j pcq, τ j pc q L 1 p0, T q.

Let us decompose r j, ¡ r j in the following way:

r j, ¡ r j u j, pc qpτ j pc q ¡ τ j pcqq, pu j, pc q ¡ u j qτ j pcq, ¡ v j, pc qpτ ¡ j pc q ¡ τ ¡ j pcqq, ¡ pv j, pc q ¡ v j qτ ¡ j pcq.

Then dφ C V c p0, T q, | ³ T 0 pr j, pc psqq ¡ r j ps, cpsqqqφpsq ds| ¤ }τ j pc q ¡ τ j pcq} L 1 p0,Tq }φ} L V p0,Tq | ³ T 0 pu j, pc psqq ¡ u j psqqτ j pcpsqqφpsq ds| }τ ¡ j pc q ¡ τ ¡ j pcqq} L 1 p0,Tq }φ} L V p0,Tq | ³ T 0 pv j, pc psqq ¡ v j psqqτ ¡
j pcpsqqφpsq ds|. Since c converges towards c in C 0 p0, T q and τ j is continuous, we get lim Ñ0 }τ j pc q ¡ τ j pcq} L 1 p0,Tq 0. Since a subsequence of u j, pc q weakly-* converges to u j in L V p0, T q and τ j pcqφ L 1 p0, T q, we also obtain lim Ñ0 ³ T 0 pu j, pc q ¡ u j qτ j pcqφ ds 0. Similar results apply to v j and τ ¡ j .

So, dφ C V c p0, T q, lim Ñ0 » T 0 pSr pc psqq ¡ Srps, cpsqqqφpsq ds 0.

(1.21)

Using both weak convergence results (1.20) and (1.21), we conclude that c, u, v are a weak solution of system (1.15).

Test case with three minerals

The problem we consider here is a set of two reactions with three minerals: M 1 é M 2 and M 2 é M 3 . This system is similar to the two reactions between calcite, calcium oxide and calcium hydroxide. The two continuous reaction rates defined by the mass action laws (1.3) are constant:

τ pτ 1 , τ 2 q.
System (1.15) becomes 6 9 9 9 9 9 9 9 9 8 9 9 9 9 9 9 9 9 7

dc 1 dt ¡r 1 , dc 2 dt r 1 ¡ r 2 , dc 3 dt r 2 , cp0q c 0 , c ¥ 0, (1.22) with r 1 u 1 τ 1 ¡ v 1 τ ¡ 1 and r 2 u 2 τ 2 ¡ v 2 τ ¡ 2 .
The kernel of S T is of dimension 1 and q p1, 1, 1q T is a basis vector satisfying q ¡ 0.

The conservative variable is c 1 c 2 c 3 so that it is possible to consider only two equations out of the three in (1.22).

The unique solution c ¦ of system (1.4) is a linear function of time:

6 8 7 c ¦ 1 ptq c 1,0 ¡ τ 1 t, c ¦ 2 ptq c 2,0 pτ 1 ¡ τ 2 qt, c ¦ 3 ptq c 3,0 τ 2 t.
It is easy to compute the first times t ¦ i ¥ 0 when the quantities c ¦ i are zero (and then negative): If

τ 1 ¡ 0 then t ¦ 1 c 1,0 τ 1 , If τ 2 ¡ τ 1 ¡ 0 then t ¦ 2 c 2,0 τ 2 ¡ τ 1 , If τ 2 0 then t ¦ 3 ¡ c 3,0 τ 2 .
(

Let c ¦

,i be the solution of the regularized system (1.19), where the function γ is identity. We analyze its convergence in different cases. We will show stronger results than Theorem 1.4.2. Indeed, c converges towards c and for any t ¡ 0, the regularized controls u pc ptqq and v pc ptqq converge towards discontinuous functions upcptqq and vpcptqq. Moreover, c is the unique solution of a discontinuous system.

Dissolution of M 1 and M 3

Here, we assume that τ 1 ¡ 0 and τ 2 0, that c 1,0 ¡ 0 and c 3,0 ¡ 0. The pure phase M 1 dissolves and produces M 2 while M 3 dissolves and produces M 2 . The regularized reaction rates are given by 4 r 1, H pc 1 qτ 1 , r 2, H pc 3 qτ 2 .

(1.24)

The regularized solutions are easy to compute in that case, since the first and third equations are independent.

Taking sufficiently small, let t 1, c 1,0 ¡ τ 1 so that 0 t 1, t ¦

1 and lim t 1, t ¦

1 . The solution c 1, is then

5 c 1, ptq c 1,0 ¡ τ 1 t c ¦ 1 ptq, t ¤ t 1, , c 1, ptq expp¡ τ 1 pt ¡ t 1, qq, t ¥ t 1, . c 1 ptq c 1,0 ¡ τ 1 t c ¦ 1 ptq, t ¤ t ¦ 1 , c 1 ptq 0, t ¥ t ¦ 1 .
Also, the control converges: dt ¥ 0, t $ t ¦

1 lim

Ñ0

H pc 1, ptqq Hpc 1 ptqq.

Clearly, results are the same for the variable c 3 and c 2 is deduced from the conservative quantity. Thus c is solution of (1.22) where

r 1 Hpc 1 qτ 1 , r 2 Hpc 3 qτ 2 .

Dissolution of M 2

Here, we assume that τ 1 0 and τ 2 ¡ 0, that c 2,0 ¡ 0. The pure phase M 2 dissolves and produces M 1 and M 3 . The regularized reactions rates are given by 4 r 1, H pc 2 qτ 1 , r 2, H pc 2 qτ 2 .

(1. [START_REF] De Dieuleveult | Un modèle numérique global et performant pour le couplage géochimie-transport[END_REF] In that case, the second equation is independent and similar to the previous case. In the same way as t 1, , we define t 2, c 2,0 ¡ τ 2 ¡ τ 1 so that 0 t 2, t ¦ 2 and t 2, converges to t ¦

2 . Thus the quantity c 2, converges for all t ¥ 0 towards c 2 defined as c 1 above (with pτ 2 ¡ τ 1 q instead of τ 1 ), and the control H pc 2, ptqq converges towards Hpc 2 ptqq for t $ t ¦ 2 .

36

It is slightly more complicated to compute the solution c 1, :

5 c 1, ptq c 1,0 ¡ τ 1 t c ¦ 1 ptq, t ¤ t 2, , c 1, ptq ¡ τ 1 τ 2 ¡ τ 1 p1 ¡ expp¡ τ 2 ¡ τ 1 pt ¡ t 2, qqq c ¦ 1 pt 2, q, t ¥ t 2, .
Therefore, for all t ¥ 0, the regularized solution converges towards c 1 defined by

4 c 1 ptq c 1,0 ¡ τ 1 t c ¦ 1 ptq, t ¤ t ¦ 2 , c 1 ptq c ¦ 1 pt 2 q, t ¥ t ¦
2 . The third quantity c 3 is obtained thanks to the conservation law. The limit c is solution of (1.22) where reaction rates are the discontinuous functions

r 1 Hpc 2 qτ 1 , r 2 Hpc 2 qτ 2 .
As illustrated in Figure 1.3, minerals M 1 and M 3 precipitate until M 2 is fully dissolved and both reactions R 1 and R 2 stop. Here, we assume that τ 2 ¡ 0 and τ 1 ¡ 0, the mineral M 1 dissolves and produces M 2 while M 2 dissolves and produces M 3 . As previously, the regularized solution c 1, converges towards c 1 defined above. We also assume that t 2, 0 or t 1, ¤ t 2, . With these assumptions, we get, as long as c 2, ptq ¥ H pc 2, q 1 and

5 c 2, ptq c 2,0 ¡ pτ 2 ¡ τ 1 qt c ¦ 2 ptq, t ¤ t 1, , c 2, ptq p1 ¡ expp¡ τ 1 pt ¡ t 1, qqq ¡ τ 2 pt ¡ t 1, q c ¦ 2 pt 1, q, t ¥ t 1, .
On the other hand, using that c 2, pt 1, q ¥ , lim

tÑ V p1 ¡ expp¡ τ 1 pt ¡ t 1, qqq ¡ τ 2 pt ¡ t 1, q c ¦ 2 pt 1, q ¡V and t Þ Ñ p1 ¡ expp¡ τ 1 pt ¡ t 1, qqq ¡ τ 2 pt ¡ t 1, q c ¦ 2 pt 1,
q is decreasing for all t ¥ t 1, , there exists a unique t2, ¥ t 1, such that c 2, p t2, q . Then, for t ¥ t2, the equation becomes dc 2 dt c 1 τ 1 ¡ c 2 τ 2 , t ¥ t2, with c 2 p t2, q . The solution is c 2, ptq pa expp¡ τ 1 pt ¡ t 1, qq p1 ¡ ab q expp¡ τ 2 pt ¡ t2, qqq where a τ 1 τ 2 ¡ τ 1 and b expp¡ τ 1 p t2, ¡ t 1, qq.

We get

6 9 9 8 9 9 7 c 2, ptq c 2,0 ¡ pτ 2 ¡ τ 1 qt c ¦ 2 ptq, t ¤ t 1, , c 2, ptq p1 ¡ expp¡ τ 1 pt ¡ t 1, qqq ¡ τ 2 pt ¡ t 1, q c ¦ 2 pt 1, q, t 1, ¤ t ¤ t2, , c 2, ptq pa expp¡ τ 1 pt ¡ t 1, qq p1 ¡ ab q expp¡ τ 2 pt ¡ t2, qqq t ¥ t2,
Let us now analyze the convergence of t2, and c 2, ptq.

We recall below the definition of t2, c 2, p t2, q , direct computations yield the following relation for all ¡ expp¡ τ 1 p t2, ¡ t 1, qqq ¡ τ 2 p t2, ¡ t 1, q c ¦ 2 pt 1, q 0, since t2, ¡ t 1, for all we have that expp¡ τ 1 p t2, ¡ t 1, qq is uniformly bounded for all , hence the first term vanishes when Ñ 0. This ensures the existence of lim Ñ0 t2, since t 1, and c ¦ 2 pt 1, q admit a limit. We set then lim Ñ0 t2, t2 . Moreover, in order to compute the limit c 2 we show that t2, ¤ t2 . Indeed, the function gpt, q : p1 ¡ expp¡ τ 1 pt ¡ t 1, qqq ¡ τ 2 pt ¡ t 1, q c ¦ 2 pt 1, q is decreasing in for all t ¥ t ¦

1 . This ensures that for all ¡ I and t ¥ t ¦ 1 , gpt, q gpt, I q. Moreover, note that gp¤, q c 2, on rt 1, , t2, s (which decreases toward ¡V when t goes to V), then necessarily t2, t2, I for all ¡ I . Passing through the limit I Ñ 0 we get t2, ¤ t2 for all ¡ 0. We get lim Ñ0 t2, t2 which satisfies t2, ¤ t2 and

τ 2 p t2 ¡ t ¦ 1 q c ¦ 2 pt ¦ 1 q. Then dt ¥ 0, lim Ñ0 c 2, ptq c 2 ptq given by 6 8 7 c 2 ptq c ¦ 2 ptq, t ¥ t ¦ 1 , c 2 ptq ¡τ 2 pt ¡ t ¦ 1 q c ¦ 2 pt ¦ 1 q, t ¦ 1 ¤ t ¤ t2 , c 2 ptq 0, t ¥ t2 .
Moreover, dt ¥ 0, t $ t ¦ 1 , t $ t2 , lim Ñ0 H pc 2, ptqq Hpc 2 ptqq.

Finally, c is solution of the discontinuous system (1.22) with r 1 Hpc 1 qτ 1 , r 2 Hpc 2 qτ 2 .

In Up to now, the regularized Heaviside functions converged towards the Heaviside function and reactions stopped when the reactant mineral was fully dissolved. Here, we will see that reactions can continue although M 2 is fully dissolved.

We assume that τ 2 ¡ τ 1 ¡ 0 and t ¦ 2 t ¦ 1 so that, for sufficiently small, 0 t 2, ¤ t 1, .

The second regularized equation is

dc 2 dt H pc 1 qτ 1 ¡ H pc 2 qτ 2 .
It is still possible to compute the regularized solution c 2, , although we do not provide all the details here:

c 2, ptq c 2,0 ¡ pτ 2 ¡ τ 1 qt c ¦ 2 ptq, t ¤ t 2, , c 2, ptq τ 2 pτ 1 pτ 2 ¡ τ 1 q expp¡ τ 2 pt ¡ t 2, qqq, t 2, ¤ t ¤ t 1, , c 2, ptq pb 1, expp¡ τ 2 pt ¡ t 2, qq b 2, expp¡ τ 2 pt ¡ t 1, qq b 3, expp¡ τ 1 pt ¡ t 2, qq b 4, expp¡ τ 1 pt ¡ t 1, qqq, t ¥ t 1, ,
where b i, are bounded epsilon-dependant constants. Thus the regularized solution c 2, converges, for all t ¥ 0, towards c 2 defined by

4 c 2 ptq c 2,0 ¡ pτ 2 ¡ τ 1 qt c ¦ 2 ptq, t t ¦ 2 , c 2, ptq 0, t ¥ t ¦ 2 .
Also, the regularized control H pc 2, q has a special limit, for any t ¥ 0,

t $ t ¦ 1 , t $ t ¦ 2 : lim Ñ0 H pc 2, ptqq Hpc 2 ptqq Hpc 1 ptqqδpc 2 ptqq τ 1 τ 2 ,
where δpxq is the Dirac function at point x. The limit c of the regularized solution satisfies the discontinuous differential equations (1.22) with r 1 Hpc 1 qτ 1 , r 2 Hpc 2 qτ 2 Hpc 1 qδpc 2 qτ 1 .

Mineral M 3 precipitates at rate τ 2 until M 2 is fully dissolved. Then, reaction R 2 does not stop but continues at rate τ 1 , so that M 2 remains fully dissolved, until M 1 is fully dissolved. In some sense, the two reactions are replaced by the unique reaction M 1 Ñ M 3 at rate τ 1 . This behaviour is shown in Figure 1.5.

Summary of test case

M 1 ¡ M 2 ¡ M 3
The other cases can be easily deduced from the four cases considered above. Finally, the regularized solution converges towards c which is solution of a discontinuous differential system. Moreover, the regularized controls converge, except at a finite number of discontinuity points, to the discontinuous ones. Theorem 1.5.1. The limit c of the regularized solution c is solution of system (1.22) where the reaction rates are discontinuous functions defined by:

4 r 1 Hpc 1 qτ 1 ¡ Hpc 2 qτ ¡ 1 ¡ δpc 2 qHpc 3 qτ ¡ 2 , r 2 Hpc 2 qτ 2 δpc 2 qHpc 1 qτ 1 ¡ Hpc 3 qτ ¡ 2 .
(1.26)

Moreover, the regularized limiters u pc q and v pc q strongly converge, up to a subsequence, in L p p0, T q, p r1, Vr.

Proof. Existence is proved above by construction. We proved in all cases that u pc q converges almost everywhere to u in r0, T s. Thanks to the dominated convergence theorem applied to |u pc ptqq ¡ uptq| p , we get strong convergence in L p .

We observe, in all cases of this example, that when a mineral disappears at a given time, then it remains fully dissolved afterwards. In other words, if c i pt ¦ q 0 then c i ptq 0, t ¥ t ¦ so that dc i dt 0, t ¡ t ¦ . Finally, the reaction rates satisfy pSrq i 0, t ¡ t ¦ . For example, when c 2 becomes fully dissolved, then r 2 r 1 .

Numerical experiments

In this section we describe and discuss the results of numerical experiments done with synthetic geochemical systems. We developed a Matlab code for solving the regularized ODE (1.19) and we study numerical convergence by varying the value of the parameter.

We define six test cases. The first one is the reference example A ¡ M defined in Section 1.3.4. The second one is M 1 ¡ M 2 ¡ M 3 analyzed in Section 1.5. Then we replace M 1 and M 3 by aqueous species A 1 and A 3 in the third test case. The fourth one (respectively the fifth one) corresponds to the previous second test case (respectively previous third test case) for which we add one mineral species M 4 and one reaction M 2 é M 4 (respectively one aqueous species A 4 and one reaction M 2 é A 4 ). With this construction, both reactions consuming M 2 in both, fourth and fifth, cases are in competition. Finally, the last test case contains seven minerals and five reactions such that one of them contains two minerals in the reactive side. In this configuration, we show numerically the non-uniqueness of the limiter u according the the choice of the regularized limiter u .

Numerical context 1.6.1.1 ODE solver

In our experiments, we used the linear function defined in (1.16). We also implemented an exponential smoothing and a polynomial one with similar results. The value of is a parameter that varies to analyze the numerical convergence.

For cases with more than one mineral per reaction we introduced a blending function (1.18) in the regularized model. In the first experiments, there is no need of blending. In the last one, we compare the min and the prod functions, in order to observe the effect of blending. The blending product is defined as follows We used the Matlab solver ode15s which implements an implicit scheme with variable order in order to solve stiff differential equations (see [START_REF] Shampine | The matlab ode suite[END_REF] and Matlab documentation for ode15s). An adaptive time step uses tolerance options RelTolSolver and AbsTolSolver, both taken as the default value 10 ¡6 . The Jacobian matrix of the nonlinear right-hand side is computed by a finite difference method. The time interval r0, T s and the initial condition c 0 ¥ 0 are parameters.

6 9 8 9 7 if M ¡ j r, then u j, pcq 1, else u j, pcq ± kM ¡ j pH pc k qq, if M j r, then v j, pcq 1, else v j, pcq ±
We define a set of discrete times in r0, T s:

t k k dt, k 0, ¤ ¤ ¤ , N f with t N f T and dt T N f
where dt is an external time step. We get as output of ode15s a set of discrete values c pt k q, where c pt k q is a vector of size N a N p .

Physical model

The maximal reaction rates are implemented in a Matlab function, using data of the physical model, mainly stoichiometry and constants of equilibrium. In our tests, the activity of a mineral species is equal to 1 and the activity of an aqueous species is equal to the molar fraction:

a i pc a q c i c H 2 O Na °i1 c a,i , (1.28)
where c H 2 O is the quantity of water. This definition assumes that the coefficient of activity is equal to 1. Moreover, we consider a dilute solution such that c H 2 O is considered constant.

Discrete error norms

In some tests, we know the limit c of the regularized model, so that c is an exact solution of the differential inclusion model. Therefore we can compute a discrete error for each species i:

e i pt k q |c i pt k q ¡ c ,i pt k q|, k 0, ¤ ¤ ¤ , N f ,
and a global error:

ept k q Na Np i1 e i pt k q.
Then we compute the L p -norm with p 1, 2, V.

We also analyze the errors in the limiting discontinuous functions u ,j and v ,j for each reaction R j . In some tests, we also know the limits u j and v j of u ,j and v ,j . The number N l of discontinuity points is finite and the corresponding times t ¦ k , k 1, ¤ ¤ ¤ , N l , are known. As above, we compute discrete error for each reaction R j and L p -norm:

e u j pt k q |u j pt k q ¡ u ,j pt k q|, k 0, ¤ ¤ ¤ , N f , e v j pt k q |v j pt k q ¡ v ,j pt k q|, k 0, ¤ ¤ ¤ , N f .
We also compute the discrete error at times t k outside of rt ¦ k ¡δ, t ¦ k δs, dk 1, ¤ ¤ ¤ , N l , where δ ¡ 0. We call out-of-jumps norms the L p -norms computed in this subset.

Reference example A ¡ M , total dissolution

In this section we consider the reference example of Section 1.3.4:

A é M,
with A a solute in the water phase w : tH 2 O, Au and M a mineral. The stoichiometric matrix is S ¡1 1 ¨T .

We consider the total dissolution case, as in Section 1.3.4.4, with a different reaction rate. The initial state is c 0 p0.2 0.6q and c H 2 O 2. In particular, the water is not involved in the reaction hence the quantity of water is constant (its value is 2). The activity of M is equal to 1 (as every mineral). The activity of A is defined by (1.28) and the reaction rate by (1.3), so that

apc A q : c A c A c H 2 O
and τ pc A q : αapc A q ¡ β Table 1.1 contains a summary of physical parameters of this test case. We can observe that errors are localized around the discontinuity. For the concentrations, the maximum value of the error and the length of the interval around t ¦ for which the error is not negligible decreases when tends towards 0. For the limiter u, the area decreases when tends towards 0 but the error remains large at the discontinuity point. Indeed, lim

Reaction A é M S ¢ ¡1 1 pα, βq p2, 2q c H 2 O 2 c 0 p0.2 0.6q Table 1.1 -Test case A ¡ M : physical parameters

Ñ0

H pc M, ptqq Hpc M ptqq outside a neighbourhood of t ¦ .

In Figure 1.8, we show the convergence of the concentration in L 1 and L V , as proved in Lemma 1.4.3. The order of convergence is about 1.85 for the L 1 -norm and 1 for the L V -norm. For the limiter, the numerical result shows in Figure 1.8 are more general in this case. In fact, there is convergence in L 1 and in L V outside the discontinuity. We do not represent the values of L 1 -norm without the discontinuity but the convergence rate is faster than L 1 -norm with the discontinuity. The convergence rates of the error for the seven first values of are almost 0.93 for the L 1 -norm, this value increases towards 1 if we consider the five last values of .

Test case

M 1 ¡ M 2 ¡ M 3
The following example is described in Section 1.5, M 1 é M 2 and M 2 é M 3 , with M 1 , M 2 and M 3 are three pure phase minerals. The stoichiometric matrix is S ¢ ¡1 1 0 0 ¡1 1 T . We consider the case described in Section 1.5.4, where a mineral, here M 2 , precipitates in one reaction and is dissolved in a second one.

To compute the reaction rate we use α p4 12q T and β p3 8q T , both quantities are defined in (1.3). Then τ 1 1 and τ 2 4. In other words, this test case corresponds to the following choice for physical quantities: K p 1 4 1 3 1 2 q T and k f k b p1 4q T . In this case, both reaction rates are positive. The initial state is c 0 p1 0.6 0.4q T , with all of this parameters, we are in the configuration of Section 1.5.4. Table 1.2 contains a summary of physical parameters of this test case.

Reactions

R 1 : M 1 é M 2 R 2 : M 2 é M 3 S ¤ ¥ ¡1 0 1 ¡1 0 1 pα 1 , β 1 q pα 2 , β 2 q p4, 3q p12, 8q c 0 p1, 0.6, 0.4q Table 1.2 -Test case M 1 ¡ M 2 ¡ M 3 : physical parameters
First, we apply solver ode15s directly to the discontinuous system (1.26), with various input parameters. In any case, the solver fails, with a false result or no result at all. This example shows the need for a numerical regularization. We can see qualitatively that the solution obtained is close to the solution computed in Section 1.5.4. In order to justify this observation, we investigate the error of the regularized solution computed with our method, for a decreasing sequence of , according to the real one. Let us consider T 2, dt 2 10 4 2 ¢ 10 ¡4 and δ 0.1, defined in Section 1.6. In the graph representing the concentration, we observe a convergence when tends to 0. Moreover, the graph representing the error of the concentration shows that the error is maximal at the discontinuity times and the error converges to 0 when tends to 0.

Both graphs representing the limiter indicate that the error is localized around the discontinuities. The difference with the concentration is that the maximal values of the error do not seem to decrease. For the limiter, Theorem 1.4.2 ensures only a weakly-* convergence in L V but in this case we can prove the strong convergence in L 1 . Numerically, the L 1 -norm converges with a rate almost equal to 1 but the uniform convergence is not satisfied.

Without considering the discontinuities, the limiter converges uniformly. Numerically, we obtain a really fast convergence outside the jumps.

Test case

A 1 ¡ M 2 ¡ A 3
The following example has the same structure as the previous example. The difference here is the consideration of the aqueous species which results in the non-constant value of the reaction rates,

A 1 é M 2 and M 2 é A 3 ,
where A 1 and A 3 are aqueous species (in this case the water phase is w :

tH 2 O, A 1 , A 3 u)
and M 2 a pure phase mineral. The stoichiometric matrix is

S ¢ 0 ¡1 1 0 0 0 ¡1 1 T ,
with the first line corresponds to H 2 O. In the sequel, we consider only the case where the reaction rates are positive and the M 2 species disappears before the first reaction reaches equilibrium.

The reaction rates are computed using α p30 4.5q T and β p3 18.75q T , both quantities are defined in (1.3). In other word, this test case corresponds to the following choice for physical quantities: K p 1 1 1 30 1 3 2 25 q T and k f k b p1 1.5q T . The initial state is c 0 p2 5 0q T , the quantity of the solvent H 2 O is more than 90 percent of the total initial quantity of the water phase (c H 2 O 20).

The activity of water solute A 1 and A 3 is given by (1.28) and reaction rates by (1.3):

a i pc a q c A i c A 1 c A 3 c H 2 O
, i 1 or 3, τ 1 pc a q α 1 a 1 pc a q ¡ β 1 and τ 2 pc a q α 2 ¡ β 2 a 3 pc a q Table 1.3 contains a summary of physical parameters of this test case.

Reactions We observe that, as in Section 1.6.3, when the central species M 2 disappears, the value of the limiter u 2 is such that r 2 r 1 τ 1 , i.e. u 1 1 and u 2 τ 1 τ 2 .

R 1 : A 1 é M 2 R 2 : M 2 é A 3 S ¤ ¥ ¡1 0 1 ¡1 0 1 pα 1 , β 1 q pα 2 , β 2 q p30, 3q p4.5, 18.75q c 0 p2, 5, 0q c H 2 O 20 Table 1.3 -Test case A 1 ¡ M 2 ¡ A 3 : physical parameters
This example illustrates the behavior of limiters in the case where the reaction rates are not constant functions.

Test case with competition

The aim of the following examples are to observe reactions in competition. We want to specify the solution given by our model in the case where a reaction produces a mineral species and the two other reactions consume this one. The interesting case is when the species dissolved by the competition disappears. In fact, in the case M 1 ¡ M 2 ¡ M 3 we prove that reaction does not stop; the same result holds true in the case of competition. So, the main issue is to determine the new reaction rate value after the disappearance. We numerically prove that our model gives as value to the limiter the evenly distributed theoretical value (1.29). We consider the following example, M 1 é M 2 and M 2 é M 3 and M 2 é M 4 with M 1 , M 2 , M 3 and M 4 four pure phase minerals. The stoichiometric matrix is

S ¤ ¥ ¡1 1 0 0 0 ¡1 1 0 0 ¡1 0 1 T .
In the sequel we consider only the case where the reaction rates are positive and the M 2 species disappears first.

The reaction rates are computed using α p4 1.5 6q T and β p3 1 4q T , both quantities are defined in (1.3), so that τ 1 1, τ 2 0.5 and τ 3 2. In other word, this test case corresponds to the following choice for physical quantities:

K p 1 4 1 3 1 2 1 2 q T and k f k b p1 0.5 2q T . The initial state is c 0 p1 0.6 0.3 0.2q T .
Table 1.4 contains a summary of physical parameters of this test case.

Reactions We recall that the limiters for the second and the third reaction are the same, as defined in (1.17). Now, let us consider the linear regularization with 0.001. The results represented in Figure 1.13 allow us to extend those observed for the case

R 1 : M 1 é M 2 R 2 : M 2 é M 3 R 3 : M 2 é M 4 S ¤ ¦ ¦ ¥ ¡1 0 0 1 ¡1 ¡1 0 1 0 0 0 1 pα 1 , β 1 q pα 2 , β 2 q pα 3 , β 3 q p4, 3q p1
M 1 ¡ M 2 ¡ M 3 .
Here, mineral M 2 is fully dissolved at time t ¦ 2 0.4, but reaction R 1 still continues until time t ¦ 1 1, when mineral M 1 is fully dissolved. We observe a numerical convergence of the limiters. Between times t ¦ 2 and t ¦ 1 , we get u 1 1 and

u 2 u 3 u 1 τ 1 τ 2 τ 3 u 1 1 0.5 2 0.4u 1 , if c 2 0.
(1.29) After time t ¦ 1 , we get u 1 u 2 u 3 0. In this example, reactions R 2 and R 3 are limited proportionally to the ratio of consumption and production. Let us consider a similar test case including aqueous species as described below,

A 1 é M 2 and M 2 é A 3 and M 2 é A 4 ,
with M 2 a pure phase mineral and A 1 , A 3 and A 4 are species in aqueous phase w : tH 2 O, A 1 , A 3 , A 4 u. The stoichiometric matrix is the same as before. The reaction rates are computed using α p30 4.5 4.5q T and β p3 18.75 45q T , both quantities are defined in (1.3). The activity of water solute A 1 , A 3 and A 4 is given by (1.28) and reaction rates by (1.3):

a i pc a q c A i c A 1 c A 3 c A 4 c H 2 O
, i 1, 3 or 4, τ 1 pc a q α 1 a 1 pc a q ¡ β 1 , τ 2 pc a q α 2 ¡ β 2 a 3 pc a q and τ 3 pc a q α 3 ¡ β 3 a 4 pc a q.

In this case, the three reaction rates are positive. The initial state is c 0 p5 2 0 0q T and the intial quatity of

H 2 O is c H 2 O 20.
Table 1.5 contains a summary of physical parameters of this test case.

Reactions

R 1 : A 1 é M 2 R 2 : M 2 é A 3 R 3 : M 2 é A 4 S ¤ ¦ ¦ ¥ ¡1 0 0 1 ¡1 ¡1 0 1 0 0 0 1 pα 1 , β 1 q pα 2 , β 2 q pα 3 , β 3 q p30, 3q p4.5, 18.75q p4.5, 45q c 0 p5200q c H 2 O 20
Table 1.5 -Test case with competition and aqueous species: physical parameters Numerically, we obtain a similar convergence of the limiters. Indeed, they converge towards u 2 and u 3 given by: u 2 u 3 u 1 τ 1 τ 2 τ 3 , if c 2 0, with r i u i τ i , for i 1, 2, 3 and u 1 1.

Test case with blending (non-uniqueness)

The aim of this section is to show the non-uniqueness of the limiter. We consider the following example,

M 3 é M 1 and M 4 é M 2 and M 1 M 2 é M 5 and M 1 é M 6 and M 2 é M 7 ,
with M i , i 1, ¤ ¤ ¤ , 7, are pure phase minerals. In the sequel we consider only the case where reaction rates are positive and M 1 and M 2 disappears first. The reaction rates are computed using α p1.5 2.1 4 2 2q T and β p1 1.4 3 1 1q T , both quantities are defined in (1.3), τ 1 0.5, τ 2 0.7, τ 3 1, τ 4 1 and τ 5 1. Table 1.6 contains a summary of physical parameters of this test case. Here, reaction R 3 involves two minerals as reactants, so that the regularized limiter u 3, is defined with a blending function. We use here the min function as in (1.18) and the prod function as in (1.27).

Reactions Chapter 2

R 1 : M 3 é M 1 R 2 : M 4 é M 2 R 3 : M 1 M 2 é M 5 R 4 : M 1 é M 6 R 5 : M 2 é M 7 S ¤ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¥ 1 0 ¡1 ¡1 0 0 1 ¡1 0 ¡1 ¡1 0 0 0 0 0 ¡1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 pα 1 , β 1 q pα 2 , β 2 q pα 3 , β 3 q pα 4 , β 4 q pα 5 , β 5 q p1.5, 1q p2

Analysis of a discontinuous ODEs model for multiphase chemical kinetics 2.1 Introduction

In this chapter we present a differential inclusion model describing kinetic reactions. The idea is to consider the right-hand side as a multivalued map, see [START_REF] Filippov | Differential equations whose right-hand is discontinuous on intersecting surfaces[END_REF], [START_REF]Differential Equations with Discontinuous Righthand Sides[END_REF] and [START_REF] Utkin | Sliding modes in control and optimization[END_REF] for an introduction to this theory. In [START_REF] Agosti | Analysis of a model for precipitation and dissolution coupled with a darcy flux[END_REF], [START_REF] Bouillard | Diffusion with dissolution and precipitation in a porous medium: Mathematical analysis and numerical approximation of a simplified model[END_REF] and [START_REF] Hoffmann | Existence and uniqueness of a global solution for reactive transport with mineral precipitation-dissolution and aquatic reactions in porous media[END_REF], authors consider discontinuous right-hand side for kinetic reactions: a discontinuous ODEs formulation in [START_REF] Bouillard | Diffusion with dissolution and precipitation in a porous medium: Mathematical analysis and numerical approximation of a simplified model[END_REF] and [START_REF] Hoffmann | Existence and uniqueness of a global solution for reactive transport with mineral precipitation-dissolution and aquatic reactions in porous media[END_REF], a complementarity formulation in [START_REF] Hoffmann | Existence and uniqueness of a global solution for reactive transport with mineral precipitation-dissolution and aquatic reactions in porous media[END_REF] and a differential inclusion formulation in [START_REF] Agosti | Analysis of a model for precipitation and dissolution coupled with a darcy flux[END_REF] and [START_REF] Hoffmann | Existence and uniqueness of a global solution for reactive transport with mineral precipitation-dissolution and aquatic reactions in porous media[END_REF]. In all of these papers, authors prove existence, positivity and uniqueness of the solution. All of these models assume that reactions cannot contain more than one mineral and every mineral can be present in only one chemical reaction.

In differential inclusion theory, the difficulty is that the right-hand side of the ODE system is no more uniquely defined on the discontinuity surfaces: it is a set-valued function. Hence, an important step of the analysis is to describe the trajectory of a solution. This is done by investigating the properties of the right hand side of the ODE system, which can be considered as a vector field. In [START_REF]Sliding modes in intersecting switching surfaces, i: Blending[END_REF], [START_REF] Dieci | A Filippov sliding vector field on an attracting co-dimension 2 discontinuity surface, and a limited loss-of-attractivity analysis[END_REF], [START_REF] Dieci | Regularizing piecewise smooth differential systems: Co-dimension 2 discontinuity surface[END_REF], [START_REF] Dieci | Sliding motion in Filippov differential systems: Theoretical results and a computational approach[END_REF], [START_REF]Sliding motion on discontinuity surfaces of high co-dimension. a construction for selecting a Filippov vector field[END_REF] and [START_REF] Guglielmi | Classification of hidden dynamics in discontinuous dynamical systems[END_REF], authors characterize different vector field configurations on discontinuity surfaces. In this chapter, we use the result presented in these papers to define the right-hand side of the differential inclusion model for kinetic systems for which only one mineral can disappear.

In chapter 1 (Section 1.4), we defined a differential inclusion system for cptq ¥ 0 : (2.1)

The reaction rate rpt, cq is a continuous function on the domain r0, T s ¢ G where G tc R Na Np | di t1, ¤ ¤ ¤ , N a N p u, c i ¡ 0u is the positive domain. However, r is only defined as a set-valued function on the domain r0, T s ¢ fG , which is the domain where at least one mineral is completely dissolved. The objective of this chapter is to better characterize our problem using Filippov's theory [START_REF]Differential Equations with Discontinuous Righthand Sides[END_REF], in order to use the relevant results and theorems to characterize the solution. First, we extend the model by defining the reaction rate for strict negative mineral concentrations. This allows us to use Filippov's theory [START_REF]Differential Equations with Discontinuous Righthand Sides[END_REF] to more precisely define the possible values of r on r0, T s ¢ fG . With this definition, we can then prove the existence of non-negative solutions to the differential inclusion system. Second, we will study the different cases resulting from various values of the neighbouring reaction rates around a discontinuity. We precise the cases for which Filippov's original definition leads to uniquely determined reaction rates on r0, T s¢fG and show that this leads to the existence of a unique non-negative solution to the resulting discontinuous ODEs system that we characterize.

Differential inclusion model for chemical kinetics 2.2.1 Domains and surfaces of discontinuity

Let c pc a , c p q, where c a R Na is the vector of aqueous species concentrations and c p R Np is the vector of mineral species concentrations.

We define Σ it1,¤¤¤ ,Npu

Σ i , with Σ i tc pc a , c p q R Na Np | c p,i 0u
e k i 5 1, if dc G k , c p,i ¡ 0, 0, if dc G k , c p,i 0. Let G : kt1,¤¤¤ ,2 Np u G k .
We have G Σ R Na Np , although G and Σ do not form a partition of R Na Np since the discontinuity surfaces Σ i intersect. For c R Na Np , let d ¥ 0 be the codimension of c, that is the number of mineral components equal to 0. The intersections of several discontinuity surfaces Σ i are characterized by the codimension of each point at the intersection.

A point c Σ of codimension d ¡ 0 has 2 d neighbouring domains G k . We denote by Gpcq the set of indices of these domains. A point c in a domain G k is of codimension 0 and Gpcq tku.

Example 2.2.1 (Space decomposition for two minerals). For a system containing only two minerals, M 1 and M 2 , the space is partitioned in four domains G k , k t1, 2, 3, 4u and an interface Σ, as shown in Figure 2.1. The point c p0, 0q lies at the intersection of Σ 1 and Σ 2 . It is of codimension 2 and has 4 neighbouring spaces. The point c I p0, ¡1.3q is in Σ 1 , it is of codimension 1 and has 2 neighbouring spaces. c P p¡2.5, ¡1q is in G 3 and it is of codimension 0.

Σ C M 2 C M 1 G 1 , e 1 p1, 1q T G 2 , e 2 p0, 1q T G 3 , e 3 p0, 0q T G 4 , e 4 p1, 0q T c P Gpc P q t4u c I Gpc I q t3, 4u c Gpcq t1, 2, 3, 4u

Extension of reaction rates

The law of mass action in Chapter 1 (Section 1.2.4) is only defined for non-negative aqueous species concentrations. The activity a of an aqueous species i satisfies the following property:

5 c a,i ¥ 0 ñ a i pc a q ¥ 0. c a,i ¤ 0 ñ a i pc a q 0. (2.2)
We extend this law to negative concentrations by setting their activity to zero in that case: Definition 2.2.1 (LMA for negative concentrations of aqueous species). dj t1, ¤ ¤ ¤ , N r u, dc R Na ¢ pR ¦ q Np , τ j pcq :

α j ¹ iA ¡ j a i pc a q ¨|S i,j | ¡ β j ¹ iA j a i pc a q ¨|S i,j | , (2.3)
with α j and β j are positive constants of reactions.

We now extend the definition of the reaction rate (1.14) Hpc i q ¨τ ¡ j pcq, dj t1, ¤ ¤ ¤ , N r u,

(2.4)

where

Hpxq 5 1, if x ¡ 0, 0, if x 0.

57

This definition means that we cancel the reaction rate of any reaction that consumes a mineral whose concentration is negative. In the positive domain c p ¡ 0, all the kinetic reactions proceed as usual : rpcq τ pcq.

The extended reaction rates r, defined in (2.4), are uniquely defined in G and are undefined in Σ.

Differential inclusion kinetic model

The idea behind the discontinuous ODE theory developed in [START_REF]Differential Equations with Discontinuous Righthand Sides[END_REF] and [START_REF] Utkin | Sliding modes in control and optimization[END_REF] originates from physical considerations: the motion of an ODE solution on a discontinuous surface Σ where the vector field is not provided must depend on the well-defined vector fields in the neighbouring spaces G k . A natural condition is to define the vector field on Σ as a convex combination of the vector fields in the neighbouring spaces G k . Filippov's theory then states that the set of admissible values of the vector field in Σ is the convex hull of the limit values of neighbouring vector fields.

The limit values of neighbouring vector fields in Σ are obtained by extending by

continuity each function r |G k : G k Þ Ñ R Nr to G k . Let r k : G k Þ Ñ R Nr be this extension such that dk t1, ¤ ¤ ¤ , 2 Np u, dj t1, ¤ ¤ ¤ , N r u and dc G k , r k j pcq min iM ¡ j e k i ¨τ j pcq ¡ min iM j e k i ¨τ ¡ j pcq. (2.5)
We can now define the reaction rates on Σ as the convex hull of the neighbouring r k : Definition 2.2.3 (Extended reaction rates in R Na Np ). The reaction rates r : r0,

1s 2 Np ¢ R Na Np Þ Ñ R Nr are defined by, dc R Na Np , rpλ, cq ķGpcq λ k r k pcq, λ Λpcq, (2.6)
with

Λpcq tλ r0, 1s 2 Np , ķGpcq λ k 1, λ i 0 if i Gpcqu. (2.7)
The Filippov set-valued vector field is then F pcq 2 Srpλ, cq, λ Λpcq @ . Definition 2.2.4 (Kinetic Differential Inclusion system). We define a kinetic model as a Differential Inclusion system:

6 8 7 dc dt 2 Srpλ, cq, λ Λpcq @ , cp0q c 0 ¥ 0, (2.8)
with r defined by (2.6).

Remark 2.2.1. If c G k then Λpcq is reduced to a singleton, Λpcq t1u, and rpλ, cq r k pcq so that F pcq tSr k pcqu. In particular, rpλ, cq τ pcq if c p ¡ 0.

We introduce below the definition of a solution to the problem (2.8).

Definition 2.2.5 (Kinetic Differential Inclusion solution). For a given initial condition c 0 ¥ 0 and a final time T ¡ 0, c : r0, T s Þ Ñ R Na Np is a solution of (2.8) if and only if c is an absolutely continuous function and there exists a measurable function λ : r0, T s Þ Ñ r0, 1s 2 Np such that for all t r0, T s, λptq Λpcptqq and

cptq c 0 S » t 0 rpλpsq, cpsqq ds. ( 2.9)
We note that Filippov's definition satisfies the conditions of the limited model (1.15) introduced in Chapter 1.

Proposition 2.2.1 (Link between (1.15) and (2.8)). Let r be defined by (2.6). dc ¥ 0 and dλ Λpcq,

r j pλ, cq u j τ j pcq ¡ v j τ ¡ j pcq, dj t1, ¤ ¤ ¤ , N r u with u j ķGpcq λ k min iM ¡ j e k i ¨ U j pcq and v j ķGpcq λ k min iM j e k
i ¨ V j pcq, where U j and V j are defined in Chapter 1 (Definition 1.4.1).

Proof. Let c ¥ 0, let λ Λpcq and let j t1, ¤ ¤ ¤ , N r u,

r j pλ, cq °kGpcq λ k r k j pcq °kGpcq λ k min iM ¡ j e k i ¨τ j pcq ¡ °kGpcq λ k min iM j e k i ¨τ ¡ j pcq u j τ j pcq ¡ v j τ ¡ j pcq, with u j ķGpcq λ k min iM ¡ j e k i ¨and v j ķGpcq λ k min iM j e k i ¨.
It remains to prove that u j U j pcq and v j V j pcq. Let us prove it for u j (a similar proof can be done for v j ): by construction u j r0, 1s because ķGpcq λ k 1 and min

iM ¡ j e k i r0, 1s. If c R Na Np such that di M ¡ j , c i ¡ 0 then dk Gpcq and di M ¡ j , e k i 1. So, u j °kGpcq λ k 1.
Remark 2.2.2. The Differential Inclusion model (2.8) is more restrictive than model (1.15). Indeed, the limiters u and v are now linearly dependent of the parameter λ.

Example 2.2.2 (Differential Inclusion formulation for the reference example A é M ). This example is described in Section 1.3.4. Here, tSτpcqu,

G 1 : tc R 2 , c M ¡ 0u, G 2 : tc R 2 , c M 0u and Σ : tc R 2 , c M 0u.
if c G 1 , tSτ pcqu, if c G 2 , tS λ 1 τ pcq λ 2 τ pcq ¨, λ 1 λ 2 1, λ r0, 1s 2 u, if c Σ, cp0q c 0 ¥ 0. (2.10)
We show that (2.8) leads to desirable properties of physical solutions such as existence and non-negativity.

Essential positivity property

Let us show that the extended reaction rates in the negative domain maintain the solution in the positive domain.

Proposition 2.2.2 (Essential Positivity). dc R Np Na , di t1, ¤ ¤ ¤ , N p N a u such that c i 0, pSrpλ, cqq i ¥ 0, dλ Λpcq. Proof. Let i t1, ¤ ¤ ¤ , N p N a u, let c R Na Np such that c i 0.
We recall that for every λ Λpcq, rpλ, cq

ķGpcq λ k r k pcq with ķGpcq λ k 1, λ k r0, 1s, dk Gpcq. Hence, pSrpλ, cqq i jt1,¤¤¤ ,Nru ķGpcq λ k S i,j ¡ ¹ lM ¡ j e k l © τ j pcq ¡ jt1,¤¤¤ ,Nru ķGpcq λ k S i,j ¡ ¹ lM j e k l © τ ¡ j pcq,
with τ j pcq defined in (2.3). Let k Gpcq and j t1, ¤ ¤ ¤ , N r u. There are two possibilities for the indices i, the species is a mineral or is an aqueous species.

1. Let's consider the aqueous case i t1, ¤ ¤ ¤ , N a u, 1.1. If i A j , we have S i,j ¡ 0, using the assumption (2.2),

c i 0 ñ a i pc a q 0 ñ τ j pcq α j ¹ i I A ¡ j a i Ipc a q ¨|S i I ,j | ¥ 0.
Then τ j pcq τ j pcq and τ ¡ j pcq 0.

So,

6 9 9 8 9 9 7 λ k S i,j ¡ ± lM ¡ j e k l © τ j pcq ¥ 0, λ k S i,j ¡ ± lM j e k l © τ ¡ j pcq 0. 1.2. If i A ¡
j then S i,j 0, using the assumption (2.2),

c i 0 ñ a i pc a q 0 ñ τ j pcq ¡β j ¹ i I A j a i Ipc a q ¨|S i I ,j | ¤ 0.
Then ¡τ ¡ j pcq τ j pcq and τ j pcq 0.

So,

λ k S i,j ¡ ± lM ¡ j e k l © τ j pcq 0, ¡λ k S i,j ¡ ± lM j e k l © τ ¡ j pcq ¥ 0. 60 1.3. Else, S i,j 0 obviously λ k S i,j ¡ ± lM ¡ j e k l © τ j pcq λ k S i,j ¡ ± lM j e k l © τ ¡ j pcq 0. 2. Let's consider the mineral case i tN a 1, ¤ ¤ ¤ , N p N a u, 2.1. If i M j then S i,j ¡ 0. c i 0 ñ ¹ lM j e k l 0.
So,

6 9 9 8 9 9 7 λ k S i,j ¡ ± lM ¡ j e k l © τ j pcq ¥ 0, λ k S i,j ¡ ± lM j e k l © τ ¡ j pcq 0. 2.2. If i M ¡ j then S i,j 0. c i 0 ñ ¹ lM ¡ j e k l 0. So, 6 9 9 8 9 9 7
λ k S i,j ¡ ± lM ¡ j e k l © τ j pcq 0, ¡λ k S i,j ¡ ± lM j e k l © τ ¡ j pcq ¥ 0. 2.3. Else, S i,j 0 obviously λ k S i,j ¡ ± lM ¡ j e k l © τ j pcq λ k S i,j ¡ ± lM j e k l © τ ¡ j pcq 0.
Summing over all terms, we have pSrpλ, cqq i ¥ 0, for every λ Λpcq.

Existence of solutions

As in Chapter 1, let Q be a basis of the kernel of S T such that Q T S 0 and let C 0 tc pR q Na Np , Q T c Q T c 0 u. We have the following result: Theorem 2.2.1 (Mass conservation and positivity). Let c : r0, T s Þ Ñ R Na Np be a solution of (2.8). Then dt r0, T s, cptq C 0 . Proof. Let c : r0, T s Þ Ñ R Na Np be a solution of (2.8).

By multiplying (2.9) by Q T , we get:

Q T cptq Q T c 0 » t 0 Q T Srpλpsq, cpsqq ds Q T c 0 .
There remains to prove that c is positive. By contradiction, let us assume that there exists t I s0, T s and i t1, ¤ ¤ ¤ , N a N p u such that c i pt I q 0. Let T i tt r0, t I s, c i ptq 0u. This set is non-empty because c i is continuous and c i p0q ¥ 0. It is closed and bounded because c i is continuous and T i r0, t I s c ¡1 i pt0uq.

Therefore, T i is an non-empty compact set and admits a maximum element t ¦ . Using the continuity of c i , we deduce that dt st ¦ , t I s, c i ptq 0.

On the other hand, using (2.9), we have:

c i pt I q c i pt ¦ q » t I t ¦ dc i psq dt ds » t I t ¦ f i |tc i 0u pλpsq, cpsqq ds, where f i pSrq i .
Hence, by using the essential positivity Proposition 2.2.2, we get c i pt I q ¥ 0 which contradicts our initial assumption.

As in Chapter 1, we assume that there exists a conservative variable q T c defined by a positive vector q. Hypothesis 2.2.1 (Positive conservative variable).

hq KerpS T q, q ¡ 0.

(2.11)

We enunciate below an existence and boundedness result for the solution of (2.8)

Lemma 2.2.1 (Bounded constraints). Under Hypothesis 2.2.1, C 0 is bounded. Proof. See Lemma 1.2.1.
Theorem 2.2.2 (Existence). Under Hypothesis 2.2.1, for all c 0 pR q Na Np , there exists at least one solution c of (2.8). Moreover, c : r0, VrÑ C 0 .

Proof. Let c 0 pR q Na Np and let F the right-hand side of (2.8). Let B 0 be the unit ball of R Na Np . By construction, for all c C 0 B 0 , F pcq is closed, convex, non-empty and bounded. Moreover, using Lemma 3 in [START_REF]Differential Equations with Discontinuous Righthand Sides[END_REF] page 67, F is upper semicontinuous. The assumptions of Theorem 1 in [START_REF]Differential Equations with Discontinuous Righthand Sides[END_REF] page 77 are then satisfied.

Therefore, there exists T r0, Vr and c : r0, T s Ñ C 0 B 0 an absolutely continuous function solution of dcptq

dt F |C 0 B 0 pcptqq, cp0q c 0 , with T ¥ 0.5 sup Bpc 0 ,0.5q |F pcq| .
Combining the mass conservation property and the positivity proved in Theorem 2.2.1, every solution of (2.8) starting in C 0 remains in C 0 . Hence, c : r0, T s Ñ C 0 . Moreover, we have T ¥ 0.5 sup C 0 B 0 |F pcq| ¡ 0. Using that the existence time T is bigger than a strict positive constant and by repeating the same process considering cpT q as the initial time (by construction the set C 0 is the same), we prove that there exists c : r0, 2T s Ñ C 0 solution of (2.8). This concludes the proof.

Uniqueness of the solution when only one mineral disappears

We study the case where only one mineral can be dissolved. Without loss of generality, we denote M 1 this mineral and c 1 c p,1 the concentration of

M 1 . Let G 1 tpc a , c p q, c a ¥ 0, c p ¡ 0u, G 2 tpc a , c p q, c a ¥ 0, c p,1 0, c p,i ¡ 0, i ¥ 2u and Σ 1 tpc a , c p q, c a ¥ 0, c p,1 0, c p,i ¡ 0, i ¥ 2u. Thus any c Σ 1 is of codimension 1. η 1 ° 1 G 1 e 1 1 1 G 2 e 2 1 0

Figure 2.2 -Space decomposition in codimension 1

Let f 1 (respectively f 2 ) be the extended right-hand side of (2.8) in G 1 (respectively G 2 ).

The set-valued function F is defined as follows

F pcq 6 9 8 9 7 tf 1 pcqu, c G 1 , tf 2 pcqu, c G 2 , tλf 1 pcq p1 ¡ λqf 2 pcq, λ r0, 1su, c Σ 1 , (2.12)
where f 1 pcq Sτ pcq and f 2 pcq Sr 2 pcq.

The normal components f 1 1 and f 2 1 are given by:

6 9 9 8 9 9 7 f 1 1 pcq jt1,¤¤¤ ,Nru S 1,j τ j pcq, c G 1 f 2 1 pcq jt1,¤¤¤ ,Nru ¡ S ¡ 1,j τ ¡ j pcq S 1,j τ j pcq © , c G 2 . (2.13)
Before describing the possible configurations of the kinetic fields we present the following result Proposition 2.3.1 (Partial attractiveness of kinetic fields). The vector field satisfies the following assertions:

-dc G 2 Σ 1 , f 2 1 pcq ¥ 0, -dc Σ 1 , f 2 1 pcq 0 ñ f 1 1 pcq ¤ 0. Proof.
By sum of positive terms and using Proposition 2.2.2, the assertion dc G 2 Σ 1 , f 2 1 pcq ¥ 0 holds. Let us prove the second assertion. Let c Σ 1 such that f 2 1 pcq 0,

f 1 1 pcq f 1 1 pcq ¡ f 2 1 pcq ¡ jt1,¤¤¤ ,Nru ¡ S ¡ 1,j τ j pcq S 1,j τ ¡ j pcq © ¤ 0.
That concludes the proof.

Let c Σ 1 , there are three possible configurations for the kinetic fields:

-Attractive: rf 1 1 pcq 0 pand f 2 1 pcq ¥ 0qs or rf 1 1 pcq 0 and f 2 1 pcq ¡ 0s.

-Transverse: f 1 1 pcq ¡ 0 (and f 2 1 pcq ¡ 0, thanks to Proposition 2.3.1).

-Tangent: f 1 1 pcq 0 and f 2 1 pcq 0.

In order to prove the uniqueness of the solution of (2.8) in the case of pointwise attractive or transverse field in codimension 1, we recall a classical result proved by Filippov in [START_REF]Differential Equations with Discontinuous Righthand Sides[END_REF].

Theorem 2.3.1 (Theorem 2 page 110 in [START_REF]Differential Equations with Discontinuous Righthand Sides[END_REF]). Let G R n be separated by a smooth surface Σ C 2 , with η a normal vector, into domains G 1 and G 2 . Let f and ff fx i , i 1, ¤ ¤ ¤ , n be continuous in domains r0, t f s ¢ G 1 and r0, t f s ¢ G 2 up to the boundary. Let the function

f 1 ¡ f 2 C 1 in r0, t f s ¢ Σ.
If for all pt, xq r0, t f s¢Σ at least one of the inequalities f 1 pt, xq ¤η 0 or f 2 pt, xq ¤η ¡ 0 is fulfilled, the right uniqueness of the solution in G of the following problem occurs,

5 dxptq dt F pt, xptqq, xp0q x 0 G, (2.14) with F pt, xptqq 6 9 8 9 7 tf 1 pt, xptqqu, xptq G 1 , tf 2 pt, xptqqu, xptq G 2 , tλf 1 pt, xptqq p1 ¡ λqf 2 pt, xptqq, λ r0, 1su, xptq Σ.
The following uniqueness result is a direct consequence of Theorem 2.3.1.

Corollary 2.3.1 (Uniqueness). Let G tpc a , c p q ¥ 0, c p,i ¡ 0, di ¥ 2u and let c 0 G.

If for all c in Σ 1 , f 1 1 pcq $ 0 or f 2 1 pcq $ 0 then the system (2.8) has a unique solution in

G.

Proof. Let c Σ 1 such that f 1 1 pcq $ 0 or f 2 1 pcq $ 0 then the field in c is in one of the two configurations:

rf 1 1 pcq 0 and f 2 1 pc ¥ 0qs or rf 1 1 pcq 0 and

f 2 1 pcq ¡ 0s, -f 1 1 pcq ¡ 0 and f 2 1 pcq ¡ 0. Hence dc Σ 1 , f 1 1 pcq 0 or f 2 1 pcq ¡ 0. The result holds using Theorem 2.3.1 in G tpc a , c p q R Na Np , c p,i ¡ 0, di ¥ 2u combined with Theorem 2.2.1.
Remark 2.3.1. In particular, Corollary 2.3.1 shows the uniqueness in the case where the field is attractive (or tangent) for all the points in Σ 1 .

The following result shows the positivity of mineral, other that M 1 , components. Proposition 2.3.2 (Positive mineral concentrations). Let t ¦ r0, T r and c be a solution of (2.8) such that cpt ¦ q Σ 1 . There exists t I st ¦ , T s such that for all t rt ¦ , t I s, c p,i ptq ¡ 0, di ¥ 2.

Proof. The proof is done using the continuity of c.

In the next sections we describe in more details the three possible configurations for the kinetic fields. In particular, we describe the trajectories of the solutions. such that ds rt ¦ , t I s, cpsq V G 1 . Hence, ds rt ¦ , t I s, f 1 1 pcpsqq ¤ 0. Let's prove by contradiction that ds rt ¦ , t I s, c 1 psq 0. We assume that there exists z rt ¦ , t I s such that c 1 pzq ¡ 0. Let A tt rt ¦ , zs, c 1 ptq 0u be a non-empty and compact set and let z 1 be its maximal element. By definition z 1 z and we have:

c 1 pzq c 1 pz 1 q » z z 1 f 1 1 psq ds ¤ 0,
which is a contradiction and that concludes the proof.

We then determine a unique sliding field on Σ 1 in the attractive case.

Corollary 2.3.2 (Sliding field). Under the same assumptions of Theorem 2.3.2, there exists t I ¡ t ¦ such that for almost every t rt ¦ , t I s, F pcptqq 2 f pλptq, cptqq, λptq

f 2 1 pcptqq f 2 1 pcptqq ¡ f 1 1 pcptqq @ is a singleton.
Proof. Using Theorem 2.3.2, there exists t I ¡ t ¦ such that for almost every t rt ¦ , t I s, F 1 pcptqq t0u. So, for almost every t rt ¦ , t I s, the right-hand side satisfies the following system:

6 9 8 9 7 f 1 pλptq, cptqq 0, f pλptq, cptqq λptqf 1 pcptqq p1 ¡ λptqqf 2 pcptqq, λptq r0, 1s, (2.16)
which has a unique solution

λptq f 2 1 pcptqq f 2 1 pcptqq ¡ f 1 1 pcptqq .
By hypothesis f 1 1 pcptqq ¤ 0, f 2 1 pcptqq ¥ 0 and f 1 1 pcptqq ¡ f 2 1 pcptqq $ 0 so λptq r0, 1s.

Figure 2.4 illustrates the attractive field in codimension 1. The two vector fields f 1 and f 2 satisfy the attractive condition. We draw in green the normal component of both fields and in red the convex hull F and in orange the sliding value f Σ 1 of the right-hand side of (2.16). We discuss below the relevance of Hypothesis 2.3.1. Indeed it is important to emphasize that this hypothesis can be proved in the case c

Σ 1 G 1 G 2 Σ 1 G 1 G 2 c f 2 pcq f 1 pcq F pcq f Σ 1 pcq f 1 1 pcq f 2 1 pcq c f Σ 1 pcq f 1 pcq f 2 pcq F pcq f 2 1 pcq
Σ 1 if f 1 1 pcq 0. Proposition 2.3.3 (Attractive field in a neighbourhood in the case f 1 1 0). If f is an attractive field at c Σ 1 such that f 1 1 pcq 0 then there exists a neighbourhood V of c such that dc ¦ G 1 V, f 1 1 pc ¦ q 0.
So, Hypothesis 2.3.1 holds true.

Proof. Since f 1 1 is continuous in G 1 we get the desired result. The two examples below (Example 2.3.1 and 2.3.1) correspond to attractive situations of the case M 1 M 2 M 3 described in Chapter 1, when M 2 disappears first. So here the kinetic field is discontinuous on Σ 2 defined by

f 1 pcq ¤ ¥ ¡τ 1 τ 1 ¡ τ 2 τ 2 , f 2 pcq ¤ ¥ ¡τ 1 τ 1 0 , with τ 1 ¡ 0 and τ 2 ¡ 0.
Using the analysis developed in this chapter, we get the same value for the right-hand side on the discontinuity surface.

Example 2.3.1 (case M 1 M 2 M 3 when M 2 disappears first in case f 1 1 0). Let consider the case τ 1 ¡τ 2 0. In this case the field is attractive because f 1 2 pcq 0. Then Corollary 2.3.2 yields λ τ 1 τ 2
and the unique possible value of the right-hand side F on Σ 2 tc 2 0, c 1 ¡ 0u is f Σ 2 defined as follow,

f Σ 2 ¤ ¥ ¡τ 1 0 τ 1 .
We recover the limited right-hand side defined in Chapter 1.

In the case where f 1 1 pcq 0 and f 2 1 pcq ¡ 0 at c Σ 1 , Hypothesis 2.3.1 is not necessarily satisfied: a difficult case arises when f 1 1 pcq ¡ 0 for c G 1 but f 1 1 pcq 0 for c Σ 1 , i.e. only the limit of the vector field is tangent to Σ

0, c 1 ¡ 0u is f Σ 2 defined as follow, f Σ 2 ¤ ¥ ¡τ 1 0 τ 2 .
So, using Lemma 2.3.1, dt st ¦ , t I s, c 1 ptq ¡ 0.

The following example shows a chemical system in this situation.

Example 2.3.3 (case M 1 M 2 M 3 when M 2 start at a 0 concentration in the transverse case). In the case M 1 M 2 M 3 when M 2 is the only one to disappear, both reactions have a positive orientation, τ 1 ¡ 0, τ 2 ¡ 0 and τ 1 ¡ τ 2 ¡ 0 and the concentration value of M 2 is equal to zero, c 0,2 0 (this corresponds chemically to consider M 2 absent at the start of chemical process), we get,

f 1 pcq ¤ ¥ ¡τ 1 τ 1 ¡ τ 2 τ 2 , f 2 pcq ¤ ¥ ¡τ 1 τ 1 0 .
In this case the field is transverse because f 1 2 pcq ¡ 0 and f 2 2 pcq ¡ 0. Thus, for all possible choices of F at the initial point, the solution are governed by f 1 after the initial time (until the disappearance of M 1 ):

cptq c 0 » t 0 f 1 pcpsqq ds, dt ¡ 0 and c 1 ptq ¡ 0.

Tangent field

The tangent field configuration corresponds to the case where the first component of both fields are equal to zero on a point of the discontinuity surface Σ 1 . See Figure 2.6 for an illustration. (2.18)

A direct consequence of this definition is F 1 pcq t0u for all c satisfying the tangent condition. In the tangent case we have the following Proposition.

f 1 pcq f 2 pcq F pcq c Σ 1 G 1 G 2
Proposition 2.3.4 (Tangent conditions). Let c Σ 1 such that f is tangent in c. Then the reaction rates of reactions containing M 1 are equal to zero, i.e. dj t1, ¤ ¤ ¤ , N r u, if S 1,j $ 0 then τ j pcq 0.

Proof. Let c Σ 1 such that f is tangent in c, i.e. S 1,j τ j pcq 0,

f 2 1 pcq jt1,¤¤¤ ,Nru ¡ S ¡ 1,j τ ¡ j pcq S 1,j τ j pcq © 0. Hence, f 2
1 pcq 0 ñ dj t1, ¤ ¤ ¤ , N r u, S 1,j τ j pcq 0 and S ¡ 1,j τ ¡ j pcq 0, using that

f 1 1 pcq jt1,¤¤¤ ,Nru ¡ S ¡ 1,j τ ¡ j pcq S 1,j τ j pcq © ¡ jt1,¤¤¤ ,Nru ¡ S 1,j τ ¡ j pcq S ¡ 1,j τ j pcq © , f 1
1 pcq 0 and f 2 1 pcq 0 ñ dj t1, ¤ ¤ ¤ , N r u, S 1,j τ ¡ j pcq 0 and S ¡ 1,j τ j pcq 0. Gathering previous properties we get dj t1, ¤ ¤ ¤ , N r u, S 1,j τ j pcq 0. Corollary 2.3.3 (Continuity of the field in the tangent case). Let c Σ 1 such that f is tangent in c. Then, F pcq tf 1 pcqu tf 2 pcqu. This concludes the proof.

Proof. Let c Σ 1 such that f is tangent in c. Let k t2, ¤ ¤ ¤ , N a N p u.
In the tangent case, the assumptions of Theorem 2.3.1 are not satisfied. However, the previous corollary means that the field f is actually continuous in the tangent case. Therefore, if we assume that these tangent conditions hold true in an open neighbourhood of a point c Σ 1 , then we have a unique solution of (2.8) due to the Lipschitz continuity of f in that neighbourhood. Unfortunately, this assumption does not always hold as shown in the following example: Example 2.3.4 (Reference example A é M in tangent configuration). Let denote τ αc a ¡ β the reaction rate with α ¡ 0 and β ¡ 0. We remind the discontinuous system in this case: 2 pc ¦ q are equal to zero, so that the tangent conditions are satisfied. In this case the other components are also equal zero:, f 1 1 pc ¦ q f 2 1 pc ¦ q 0. Hence, F pc ¦ q 3 ¢ 0 0 A . Uniqueness results for this particular example in the tangent case are obtained in [START_REF] Agosti | Analysis of a model for precipitation and dissolution coupled with a darcy flux[END_REF] and [START_REF] Hoffmann | Existence and uniqueness of a global solution for reactive transport with mineral precipitation-dissolution and aquatic reactions in porous media[END_REF].

In [START_REF] Agosti | Analysis of a model for precipitation and dissolution coupled with a darcy flux[END_REF], the authors use the mass conservation to prove the one-sided Lipschitz continuity of the right-hand side, which yields the uniqueness. In fact, Theorem 2.3.1 is proved using two main arguments: one-sided Lipschitz property and a suitable change of variable allowing to consider the right orthogonal space. In view of [START_REF]Differential Equations with Discontinuous Righthand Sides[END_REF], a future research direction would be to identify the right scalar product to prove the one-sided Lipschitz property of the right-hand side of (2.8) in a neighbourhood of a point in the tangent case for a general kinetic system.

Introduction

Many environmental studies rely on modeling geochemical reactions such as CO 2 sequestration, water injection, and geothermal energy. In some of these applications, when the reaction rate is relatively slow compared to the environment, the reaction does not necessarily reach an equilibrium. In this case, we consider kinetically controlled chemical reactions. In several of the aforementioned applications, the chemical reactions are coupled with hydrodynamic processes leading for instance to a reactive transport problem [START_REF] Lichtner | Reactive transport in porous media[END_REF][START_REF] Carrayrou | Comparison of numerical methods for simulating strongly nonlinear and heterogeneous reactive transport problems-the MoMas benchmark case[END_REF][START_REF] Hoffmann | A general reduction scheme for reactive transport in porous media[END_REF][START_REF] Erhel | Analysis of a global reactive transport model and results for the MoMas benchmark[END_REF][START_REF] De Dieuleveult | A global strategy for solving reactive transport equations[END_REF]. It is of great interest to get an in-depth knowledge of the geochemical properties involved as it is very often the source of non-linearity or discontinuity in the system. In particular, this is true when considering heterogeneous reactions such as precipitation-dissolution reactions [START_REF] Rubin | Transport of reacting solutes in porous media: Relation between mathematical nature of problem formulation and chemical nature of reactions[END_REF]. The difficulty here comes from the change of phase, which implies discontinuities in the reaction rate [START_REF] Friedly | Solute transport with multiple equilibrium-controlled or kinetically-controlled chemical reactions[END_REF]. We also refer

A Projected Dynamical System Model

In this section, we first recall some characteristics of kinetic chemical systems. Then we give some constraints and introduce our new model.

Kinetic system

We denote by n the number of species taking part into m reactions, where we assume that n ¡ m. The vector x R n contains the species concentrations. We consider chemical reactions composed of aqueous components and minerals. There is no restriction on the number of minerals per reaction on either side of the reaction. We denote by W t1, ¤ ¤ ¤ , nu the set of indices of the aqueous species.

S Z n¢m denotes the stoichiometry matrix, where, by convention, the stoichiometric coefficients of the reactants are taken as negative integers. The sign convention is required here to define the rate of each reaction.

Let KerpS T q be the kernel of S T . A basis Q of KerpS T q is a collection of vectors of so-called conservative quantities. Assuming that S is a full rank matrix, Q R n¢pn¡mq is also a full rank matrix satisfying ImpQq KerpS T q.

We assume that the activity of an aqueous species is equal to its concentration, whereas the activity of a mineral species is equal to 1. In these conditions, we use the kinetic mass action law to define the reaction rates τ pxq R m . We get, for each j 1, . . . , m: τ j pxq : k r,j

¹ iW s i,j 0 x ¡s i,j i ¡ k p,j ¹ iW s i,j ¡0 x s i,j i , (3.1)
where we denote by k r (resp. k p ) the constant vector of the left side (resp. of the right side) of reactions corresponding to the reactants (resp. to the products).

The table below summarizes notations used throughout the paper.

Name

Model

Reaction rates are driving the evolution of the system such that 9

x Sτ pxq and

Q T 9
x 0 as long as x ¡ 0. However, this evolution is perturbed (but not necessarily stopped) when a mineral species is vanishing, as illustrated for instance in [40, Exemple 3].

In order to design a model taking into account vanishing minerals, we first assume that the derivative of the concentration vector remains proportional to the reaction rate.

In other words, for each species i t1, . . . , nu, there exists λ i r0, 1s such that 9

x i pSτpxqq i λ i . This can be rewritten as 9

x Dpxq : tdiagpSτpxqqλ with λ r0, 1s n u. 1Moreover, concentrations must satisfy contraints. Starting from a non-negative initial condition x 0 R n , the vector xptq should remain non-negative for all t R . Besides, conservative variables Q T xptq should be invariant. Said otherwise, the concentration vector satisfies at any time

xptq C 0 : tx R n : Q T px ¡ x 0 q 0u.
The following lemma links the conserved quantities to the derivative of x. Lemma 3.1.1. Let x : r0, VrÑ R n be an absolutely continuous function such that xp0q x 0 . For almost every t r0, Vr, it holds true that Q T pxptq ¡ x 0 q 0 ðñ Q T 9 xptq 0. Proof. The function yptq : Q T pxptq¡x 0 q is also absolutely continuous with 9 yptq Q T 9 xptq almost everywhere and yp0q 0. Then, using derivation properties, we get, almost everywhere, yptq y 0 ðñ 9 yptq 0 and the desired result.

Before introducing our model, we recall two classical definitions from convex analysis.

The tangent cone to a convex set C at x C is defined by T pC; xq td : hd ν Ñ d and hτ ν Ó 0 with x τ ν d ν C dνu.

(

The projection of a vector v onto C is defined by ProjpC; vq : tw C : }v ¡ w} 2 inf zC }v ¡ z} 2 u.

Clearly, C 0 is a non-empty convex set defined by linear inequality and equality constraints. Thus, the tangent cone to C 0 at x is non-empty, closed and convex and has the following explicit form: T pC 0 ; xq td : Q T d 0, d i ¥ 0 for all i such that x i 0u.

(3.3) Hence T pC 0 ; xq Dpxq is a closed and convex set. Now, we can define our model. Combining the two conditions described above, we get the following projected dynamical system: This model is similar to a classical projected dynamical system, with projection onto the tangent cone, which has a unique solution [START_REF] Nagurney | Projected dynamical systems and variational inequalities with applications[END_REF][START_REF] Cojocaru | Projected dynamical systems and evolutionary variational inequalities via Hilbert spaces with applications[END_REF]. However, here, projection is onto the intersection of T pC 0 , xptqq with the set Dpxptqq.

Therefore, we must prove the existence of a solution satisfying our projected dynamical system (3.4). Definition 3.1.1. Let x 0 R n . A function x : r0, t f s Ñ C 0 is a solution of (3.4) if x is absolutely continuous and satisfies (3.4) for almost every t r0, t f s.

Link with Other Models

Link with a complementarity system

In the context of multiphase models, it is recurrent to use complementarity problems in the presence of heterogeneous reactions. In [START_REF] Hoffmann | Existence and uniqueness of a global solution for reactive transport with mineral precipitation-dissolution and aquatic reactions in porous media[END_REF], the authors propose a complementarity formulation of kinetic systems, assuming one mineral per reaction and one reaction per mineral.

We focus on a reference example with an aqueous species W , a mineral species M , and one reaction: W é M, so that S p¡1 1q T and Q p1 1q T . The reaction rate is given by τ pxptqq k r x 1 ptq¡k p , with two constants k r and k p .

The complementarity model introduced in [START_REF] Hoffmann | Existence and uniqueness of a global solution for reactive transport with mineral precipitation-dissolution and aquatic reactions in porous media[END_REF] is written:

6 9 8 9 7
x 1 ptq 9

x 2 ptq 0, 0 ¤ x 2 ptq u 9

x 2 ptq ¡ τ pxptqq ¥ 0, xp0q x 0 .

(3.5)

For this example, we prove that our model (3.4) is equivalent to the complementarity model (3.5). Theorem 3.2.1. The projected dynamical system (3.4) is pointwise equivalent to the complementarity system (3.5).

Proof. Let x : r0, t f s Ñ R 2 be a Lipschitz continuous function. We can split r0, t f s into a finite number of subintervals such that, in each open interval I, x is differentiable and satisfies one of the three following conditions: either dt I, x 1 ptq 0 or dt I, x 2 ptq 0 or dt I, x 1 ptq ¡ 0 and x 2 ptq ¡ 0.

In the first case, since x 1 ptq 0 over I, it follows that dt I, 9

x 1 ptq 0. Hence, noticing that the constraint 9

x 1 ptq 9

x 2 ptq 0 is common to both problems, we get dt I, 9

x 1 ptq 9

x 2 ptq 0. Therefore, both problems are equivalent in this case. The second case is similar using x 2 instead of x 1 .

In the third case, the complementarity problem (3.5) can be rewritten dt I, 9 xptq p¡τpxptqq, τ pxptqqq T .

On the other hand, ProjpDpxptqq T pC 0 ; xptqq; p¡τpxptqq, τ pxptqqq T q p¡τpxptqq, τ pxptqqq T , so that the projected dynamical system is equivalent. Finally, both models are equivalent in each open interval I, thus almost everywhere in r0, t f s.

Link with a discontinuous model

In some kinetic models, positivity is ensured by limiting each reaction rate [START_REF] Bouillard | Diffusion with dissolution and precipitation in a porous medium: Mathematical analysis and numerical approximation of a simplified model[END_REF][START_REF] Hamlat | Kinetic models with pure phases[END_REF]. In the model proposed in this paper, we choose to limit each time derivative.

For the reference example W é M above, discontinuous models proposed in [START_REF] Bouillard | Diffusion with dissolution and precipitation in a porous medium: Mathematical analysis and numerical approximation of a simplified model[END_REF][START_REF] Hamlat | Kinetic models with pure phases[END_REF] are equivalent to the complementarity problem [START_REF] Hoffmann | A general reduction scheme for reactive transport in porous media[END_REF]. So our model is also equivalent to discontinuous ones.

As said above, we also consider kinetic systems where reactions can involve several minerals and where a mineral can participate in several reactions. We consider here a kinetic system with two reactions involving three minerals [START_REF] Hamlat | Kinetic models with pure phases[END_REF]:

M 1 é M 2 and M 2 é M 3 .
The stoichiometry matrix is S ¢ ¡1 1 0 0 ¡1 1 T and a basis of KerpS T q is Q p1 1 1q T . The reaction rate is constant and given by τ pxq k r ¡ k p . For this example, we choose k r and k p so that τ p2 6q T and the initial condition is x 0 p1 1 1q T . Then Sτ p¡2 ¡ 4 6q T .

The discontinuous model proposed in [START_REF] Hamlat | Kinetic models with pure phases[END_REF] has a unique solution, which is piecewise affine. It is plotted in Figure 3.1 and defined below: p1 1 1q T , 0 ¤ t 1{4, p 1 2 0 5 2 q T , 1{4 ¤ t 1{2, p0 0 3q T , 1{2 ¤ t ¤ t f . The first conditions are easy to check: dt, xptq C 0 and 9 xptq Dpxptqq T pC 0 ; xptqq.

xptq ¤ ¥ x ¦,1 ptq ¡ 2λ 1 ptqpt ¡ zptqq x ¦,2 ptq ¡ 4λ 2 ptqpt ¡ zptqq x ¦,3 ptq 6λ 3 ptqpt ¡ zptqq with λptq 6 9 8 9 7 p1 1 1q T , 0 ¤ t 1{4, p1 0 1 3 q T , 1{4 ¤ t 1{2, p0 0 0q T , 1{2 ¤ t ¤ t f ,

It remains to verify that 9

xptq is indeed the projection yptq of Sτ pxptqq onto the set Dpxptqq T pC 0 ; xptqq for all t. This is clearly true for t r0, 1 4 r: 9 xptq yptq. Let us consider t s 1 4 , 1 2 r. Since yptq Dpxptqq there exists µ r0, 1s 3 such that yptq diagpSτ qµ p¡2µ 1 , ¡4µ 2 , 6µ 3 q T . Besides, since 9 xptq T pC 0 ; xptqq and using (3.3), it follows that µ 2 ptq 0 (recall that x 2 p 1 4 q 0 and pSτpxptqqq 2 0). Moreover, Qyptq 0 yields µ 1 3µ 3 . Since µ r0, 1s 3 , we get µ p1 0 1 3 q T . Thus 9 xptq yptq. Finally, for t s 1 2 , 1s, Dpxptqq T pC 0 ; xptqq contains only p0, 0, 0q T , since x 1 p 1 4 q x 2 p 1 4 q 0 and pSτpxptqqq 1 0, pSτpxptqqq 2 0.

In the numerical examples, we will also consider a similar system, with two aqueous species and one mineral:

W 1 é M and M é W 2 .

Link with a classical projected dynamical system

In this section, we illustrate by an example why the projection should contain explicitly the constraint 9 xptq Dpxptqq.

We introduce a classical projected dynamical system, which has a unique solution [START_REF] Nagurney | Projected dynamical systems and variational inequalities with applications[END_REF]: x Dpxptqq.

At time t 1 4 , the second mineral M 2 is fully dissolved. Then 9

x 1 ptq Sτ 1 λ 1 , with λ 1 ¡ 1 thus 9 xptq Dpxptqq. In other words, model (3.6) does not satisfy the constraint for the time derivative in our model.

Existence of a Solution

In this section, we prove the existence of a solution to (3.4). We give a constructive proof by introducing a numerical scheme. For a given time step h ¡ 0, we consider forward difference approximation of the derivative. Then, an approximate solution is defined by a projection and an affine interpolation. We prove that this explicit scheme converges as h Ñ 0 to a solution of (3.4).

Throughout this section, we denote by tau the largest integer smaller than a.

Explicit scheme and interpolation

Consider a fixed step size h ¡ 0. For k N such that kh ¤ t f , let tx k u k be the sequence defined such that x 0 is the initial point of (3.4) and

x k 1 Proj pK k ; x k hSτ px k qq , (3.7)
where K k C 0 px k hDpx k qq. Our first observation is that this sequence is well-defined. Proposition 3.3.1. Let h ¡ 0. For all k 0, . . . , t t f h u, the projection in (3.7) exists and is unique.

Proof. Note that for any k, the set K k is non-empty as x k K k . Moreover, the set K k is closed and convex. Therefore, the projection in (3.7) exists and is unique.

We now sum up in the following lemma some properties of the sequence tx k u k computed in (3.7).

Lemma 3.3.1. Let h ¡ 0, k N such that hk ¤ t f , and x k C 0 . Then z K k ðñ z ¡ x k h Dpx k q 3 d : Q T d 0, d ¥ ¡ x k h A .
Moreover, (3.7) is equivalent to

5 x k 1 x k hd k , x k 1 ¡x k h Proj Dpx k q 2 d : Q T d 0, d ¥ ¡ x k h @ ; Sτ px k q ¨. (3.8)
Proof. Let z and d be related by z x k hd. Since

x k C 0 , it holds that z C 0 ðñ Q T d 0, d ¥ ¡ x k h and z x k hDpx k q ðñ d Dpx k q. Thus z K k if and only if d Dpx k q td : Q T d 0, d ¥ ¡ x k h u.
Also, the projection in (3.7) can be rewritten as

}x k 1 ¡ x k ¡ hSτ px k q} ¤ }z ¡ x k ¡ hSτ px k q}, dz K k . Hence, dividing by h ¡ 0 and defining d k by x k 1 x k hd k , previous inequality is equivalent to }d k ¡ Sτ px k q} ¤ }d ¡ Sτ px k q}, dd Dpx k q td : Q T d 0, d ¥ ¡ x k h u.
The result follows again by definition of a projection.

Using the sequence tx k u k , for h ¡ 0 and x 0 C 0 , we define an interpolation function

x h : r0, t f s Ñ R n by x h ptq x k px k 1 ¡ x k q t ¡ kh h , with k t t h u.
(3.9) Here, k is chosen such that t rkh, pk 1qhr.

The function x h is piecewise affine and differentiable almost everywhere:

9 x h ptq x k 1 ¡ x k h ,
d : Q T d 0, d ¥ ¡ x k h @ ; Sτ px k q ¨, x h ptq C 0 , x h p0q x 0 , (3.11)
where k t t h u.

Proof. For any t r0, t f s,

x h ptq is a convex combination of x k C 0 and x k 1 C 0 . Since C 0 is a convex set, x h ptq C 0 . By definition, x h p0q x 0 .
Finally, by (3.10) and (3.8), the derivative satisfies the first equation of (3.11).

Bounded constraints

In the special case where C 0 is bounded, a solution of the model is Lipschitz continuous. We introduce an assumption limiting the production of species. More precisely, we assume that the columns of Q T are positively linearly independent. Assumption 3.3.1.

Q T d 0, d ¥ 0 ñ d 0. (3.12)
Using Motzkin alternative theorem [START_REF] Mangasarian | Nonlinear programming[END_REF], (3.12) is equivalent to the following condition:

td | Qd ¡ 0u $ r.
Since Q is a basis of KerpS T q, it can also be written tq KerpS T q, | q ¡ 0u $ r.

This last condition is used for example in [START_REF] Hoffmann | Existence and uniqueness of a global solution for reactive transport with mineral precipitation-dissolution and aquatic reactions in porous media[END_REF]Assumption 4.1] The following lemma shows that, thanks to Assumption 3.3.1, the production of concentration is necessarily finite. Proof. Remind that C 0 tx : Q T px ¡ x 0 q 0, x ¥ 0u. Assume by contradiction that there exists a sequence tx n u with x n C 0 for all n N such that lim nÑV }x n } V V.

Let d n : xn¡x 0 }xn}V . This sequence is bounded, therefore it converges, up to a subse- quence, to some limit:

lim nÑV d n lim nÑV x n }x n } V d ¥ 0 with d $ 0.
We obtain d $ 0 using that if lim nÑV }x n } V V then there exists a component of x n such that, up to a subsequence, the limit of this component is infinite. So, there exists a non-zero component of d.

Besides, since x n C 0 it holds true that Q T d n 0 for all n and then Q T d 0. However, Assumption 3.3.1 gives that d 0, which is a contradiction.

The following example illustrates a system where production may arbitrarily grow and Assumption 3.3.1 does not hold.

Example 3.3.1. Consider the following system of two reactions with three aqueous species: W 1 é W 2 and W 2 é 2W 1 W 3 . This system seems to create an infinite quantity of species, since whenever the first reaction consumes a quantity of W 1 , then the second reaction doubles this quantity.

In this case, the stoichiometry matrix and a basis of the kernel of its transpose are given by A direct consequence of this lemma is that any element in the set Dpxq is also bounded under Assumption 3.3.1. Lemma 3.3.4. Assume that (3.12) holds true. There exists ρ ¡ 0 such that, for any x C 0 and any d Dpxq, }d} ¤ ρ. Proof. Let x C 0 and d Dpxq. By definition of D, there exists λ r0, 1s n such that d diagpSτ pxqqλ. Therefore, it holds }d} }diagpSτpxqqλ} ¤ }λ}}Sτpxq} ¤ }Sτpxq}.

S ¤ ¥ ¡1 2 1 ¡1 0 1 and Q ¤ ¥ 1
Since C 0 is bounded by Lemma 3.3.3, and τ is continuous, then τ pxq is uniformly bounded. Hence, there exists ρ ¡ 0 such that }d} ¤ }Sτpxq} ¤ ρ.

Convergence of the discrete function

Throughout this section, we make Assumption (3.3.1), such that C 0 is bounded. Thanks to this assumption, we will prove convergence of x h towards a function x. Proposition 3.3.2. There exists a subsequence H and a Lipschitz continuous function

x : r0, t f s Ñ C 0 R n such that x h Ñ H x uniformly.
Proof. By Lemma 3.3.3, the sequence tx h u h is equi-bounded. Also, by Lemma 3.3.4, the sequence tx h u h is equi-Lipschitz with constant ρ.

Then, all the assumptions of the Arzela-Ascoli theorem are satisfied. Therefore, there exists a subsequence H, such that h Ñ H 0, and a Lipschitz mapping x : r0, t f s Ñ R n such that xp0q x 0 and tx h u hH converges uniformly to x in r0, t f s. Besides, since C 0 is closed, then xptq C 0 for all t r0, t f s.

Since k depends on h and t, the sequence x k converges. Lemma 3.3.5. Let t r0, t f s. The sequence x k with k t t h u converges towards xptq:

x k Ñ H xptq.
Proof. By definition of x h in (3.9) and using Lemma 3.3.4, it follows

}x h ptq ¡ x k } ¤ }px k 1 ¡ x k q t ¡ kh h } ¤ ρpt ¡ khq ¤ ρh.
Passing to the limit as h Ñ 0 , and using convergence of x h ptq in Proposition 3.3.2, we get convergence of x k towards xptq.

Another consequence of Proposition 3.3.2 is the weak convergence of the derivative. Lemma 3.3.6. 9

x h Ñ H 9

x weakly in L 1 pr0, t f s, R n q.

Proof. The result is a consequence of a straightforward application of [47, Lemma 2.2].

Even though the convergence is weak, there exists a sequence of functions which converge strongly towards the derivative of x. Lemma 3.3.7. There exists a sequence of functions ty h u h Ñ H 9

x and satisfying

y h ļ J h α l 9
x h l , with ļ J h α l 1, and α l ¥ 0, dl J h , where J h is a finite set of indices such that ph l q lJ h H. Proof. The sequence t 9

x h u hH converges weakly to 9

x, by Lemma 3.3.6. Then, applying Mazur's Lemma, there exists a sequence ty h u h , with values in the convex hull of the set t 9

x h u hH , which converges to 9

x.

The following lemma will also be useful when passing to the limit in the constraint set. The set convergence has to be understood in the Painlevé-Kuratowski sense [64, Chapter 5.B]. Lemma 3.3.8. Let z C 0 and tz h u h C 0 be a sequence such that z h Ñ hÑ0 z. Then

lim hÑ0 Dpz h q td : Q T d 0, d ¥ ¡ z h h u Dpzq td : Q T d 0, d ¥ ¡ lim hÑ0 z h h u. Proof. By Lemma 3.3.
T d 0, d ¥ ¡ z h
h u are both non-empty and compact for all h ¡ 0. Therefore, by [START_REF] Rockafellar | Variational analysis[END_REF]Theorem 5.7], it follows that lim sup hÑ0 Dpz h q td : Q

T d 0, d ¥ ¡ z h h u Dpzq td : Q T d 0, d ¥ ¡ lim hÑ0 z h h u. On the other hand, let ξ Dpzq td : Q T d 0, d ¥ ¡ lim hÑ0 z h h u. Since x Þ Ñ Sτ pxq is continuous, there exists a sequence ξ h Dpz h q td : Q T d 0, d ¥ ¡ z h h u which converges to ξ. Thus Dpzq td : Q T d 0, d ¥ ¡ lim hÑ0 z h h u lim inf hÑ0 Dpz h q td : Q T d 0, d ¥ ¡ z h h u.
All in all, the limit exists and the result follows.

Existence of a solution

Again, we make Assumption 3.3.1 throughout this section. Denote by I : tt r0, t f s : 9 xptq existsu. By Sard's theorem, the set r0, t f szI has a zero measure.

We will prove that the limit x of x h is a solution of the model. First, we show that it satisfies the constraints. Lemma 3.3.9. For all t I, 9 xptq Dpxptqq td :

Q T d 0, d ¥ ¡ lim H x k l h l u, with l J h and k l t t h l u.
Proof. Let t I and h ¡ 0. Considering y h defined in Lemma 3.3.7, and using (3.11), we get

y h ptq conv £ ¤ lJ h Dpx k l q td : Q T d 0, d ¥ ¡ x k l h l u ,
where conv denotes the convex hull and k l t t h l u.

We note that h l Ñ H 0 and, using Lemma x k l Ñ H xptq. Thus the sequence x k l ptq satisfies assumptions of Lemma 3.3.8. Passing to the limit in y h , we get the result. Lemma 3.3.10. For all t I, 9 xptq T pC 0 ; xptqq.

Proof. Recall that the tangent cone to C 0 is given by (3.3). By Lemma 3.3.9, Q T 9 xptq 0.

By Lemma 3.3.9, 9

xptq Dpxptqq so that 9

x i ptq pSτpxptqqq i λ i for i 1, . . . , n, with 0 ¤ λ i ¤ 1. Let i such that x i ptq 0. We must show that 9

x i ptq ¥ 0. If pSτpxptqqq i ¥ 0, then 9

x i ptq ¥ 0. On the other hand, if pSτpxptqqq i 0, by continuity of τ and x, there is an open interval U centered in t such that pSτpxpsqqq i 0 so that 9

x i psq ¤ 0 for all s U I. Also, x i ptq 0 and x i psq ¥ 0. Since x is absolutely continuous, we get 9

x i ptq 0. This concludes the proof. Lemma 3.3.11. The limit x satisfies the inclusion below, for t I and l J h :

Dpxptqq T pC 0 ; xptqq Dpxptqq td : Q T d 0, d ¥ ¡ lim hH x k l h l u. Proof. Since x k l ¥ 0, the limit satisfies ¡ lim hH x k l h l ¤ 0. Let i such that x i ptq ¡ 0, then ¡ lim hH px k l q i h l ¡V. So, if d T pC 0 ; xptqq, it satisfies d ¥ ¡ lim hH x k l h l .
We are now able to prove the existence of a solution to our model. Theorem 3.3.1. The limit function x is a solution of (3.4).

Proof. The function x : r0, t f s Ñ C 0 is a Lipschitz continuous function such that xp0q x 0 . It satisfies the constraints of the model, so that it remains to verify the projection property.

Considering y h defined in Lemma 3.3.7, and using the triangle inequality, we get }Sτpx k q ¡ y h ptq} ¤ max lJ h }Sτpx k q ¡ 9 x h l ptq}.

(3.13) Let j J h which maximizes this inequality and d Dpx k j q td : Q T d 0, d ¥ ¡ x k j h j u.

Since x h j is solution of (3.11), it satisfies }Sτpx k j q ¡ 9 x h j ptq} ¤ }Sτpx k j q ¡ d}, therefore }Sτpx k q ¡ y h ptq} ¤ }Sτpx k q ¡ Sτ px k j q} }Sτpx k j q ¡ d}.

T d 0, d ¥ ¡ lim hH x k j h j u.
Finally, Lemma 3.3.11 provides the result.

An Alternating Projection Algorithm

We conclude this analysis by presenting an alternating projection method to compute the sequence x k . This iterative algorithm consists in computing a sequence ty j u j , starting from x k hSτ px k q, which converges towards x k 1 . The sequence is defined as follows:

y 0 x k hSτ px k q, e 0 e ¡1 0, y 2j 1 ProjpC 0 ; y 2j e 2j¡1 q, e 2j 1 y 2j e 2j¡1 ¡ y 2j 1 , y 2j 2 x k hProj Proof. Using [28, Proposition 2.7 (ii) and (iv)], the fourth line of (3.14) can be rewritten

¡ Dpx k q; y 2j 1 e 2j ¡ x k h © , e 2j
y 2j 2 Projpx k hDpx k q; y 2j 1 e 2j q.
The algorithm is a direct application of the Boyle-Dykstra Theorem [START_REF] Deutsch | Best approximation in inner product spaces[END_REF]Theorem 9.24].

We already noted that C 0 and x k hDpx k q are closed convex set and that their intersection K k is non-empty. Therefore, the sequence y j converges and the limit is nothing else than x k 1 .

Numerical Experiments

We developed a software in Matlab in order to compute the discrete solution given by (3.9). For each k 0, . . . , t t f h u, and given x k , the new iterate x k 1 is computed with the alternating projection algorithm (3.14), implemented with an error tolerance of ¡ 0 and a maximal number of iterations j max as in Algorithm 1.

Data: starting vector x k ; iteration threshold ; 1 Begin ; 2 Set j : 0, y 0 : x k hSτ px k q, e 0 : e ¡1 : 0; Algorithm 1: Alternating projection algorithm to (3.7). Projections are computed using Matlab lsqlin function. In all experiments, the threshold is set to 10 ¡7 , while the step size h is 10 ¡2 in the time interval r0, 1s.

The sequel is divided into five examples illustrating the new approach discussed in this article. In all these examples, Algorithm 1 converges to a solution.

Example W 1 W 2

The first example is a single reaction with two aqueous species:

W 1 é W 2 .
The stoichiometric matrix is S p¡1 1q T , and a basis of kerpS T q is Q p1 1q T . The reaction rate is τ pxq k r x 1 ¡ k p x 2 . Figure 3.3 shows the results on two different instances. In the first case (left), the second species W 2 is initially absent, then increases as the first species W 1 decreases. The second case (right) illustrates an increase of W 1 with a decrease of W 2 . In both cases, the quantities tend to an equilibrium at infinite time. Asymptotic values are given by τ pxq 0 and Q T x Q T x 0 .

Example W M

The following classical example, considered earlier in Section 3.2.1, is a single reaction with one aqueous species W and one mineral M :

W é M.
As in the previous test case, S p¡1 1q T , and, Q p1 1q T . The reaction rate is τ pxq k r x 1 ¡ k p . The following example is now two reactions involving three mineral:

M 1 é M 2 and M 2 é M 3 .
The stoichiometric matrix is S ¢ ¡1 1 0 0 ¡1 1 T , and, Q p1 1 1q T is a non-zero element in the kernel of S T . The reaction rate is τ pxq k r ¡ k p . This exemple is used in Section 3.2.2 to illustrate a solution of (3.4) with discontinuous reaction rates. Figure 3.5 shows the evolution of the three solutions. In this more sophisticated example, we see a three stages reaction. Each stage is affected by the disappearance of one mineral and then follow a different reaction rate. The mineral M 2 , which appears in both reactions vanishes first. Then, the mineral M 1 reduces thanks to the constraint on the conservation until only the mineral M 3 remains and the reaction stops. The following example considers again two reactions, but involving two minerals and one aqueous species:

W é M 1 and M 1 é M 2 .
The stoichiometric matrix is S Figure 3.6 shows the dynamics of the chemical system. The aqueous species W decreases and produces the first mineral M 1 which is dissolved into the second mineral M 2 . Under the conditions chosen, M 2 is fully dissolved in a finite time. However, because the first reaction continues, the second reaction can also continue, at the reaction rate τ 1 instead of τ 2 . In some sense, the system behaves as if W precipitates into mineral M 2 , at rate τ 1 . Asymptotic equilibrium, at infinite time, is given by x 2 0, Q T x Q T x 0 and τ 1 px 1 q 0.

Concluding Discussion

The study presented in this article introduces a new model for precipitation-dissolution reactions in the form of a projected dynamical system. The main benefit of this new approach is to offer the possibility to consider reactions with several minerals. This new model extends the complementarity system model for one mineral introduced in [START_REF] Hoffmann | Existence and uniqueness of a global solution for reactive transport with mineral precipitation-dissolution and aquatic reactions in porous media[END_REF]. We presented a numerical scheme to find a solution to the differential system and proved its convergence. The existence of solutions of the equation follows as a consequence of the convergence. The numerical approach and the model have been validated on some examples.

Further researches may focus on the numerical perspectives offered by this work. The actual proposed method is based on an explicit discretization. The perspective to use an implicit scheme is an interesting path. Although in this case, the projection on D becomes more difficult.

The motivation to consider kinetic reactions is to assume that the reaction is too slow to attain its equilibrium. However, assuming that after a long time, the system reaches a stable state, we should obtain a connection with models for equilibrium reactions [START_REF] Erhel | Characterizations of solutions in geochemistry: existence, uniqueness, and precipitation diagram[END_REF].

Finally, an essential motivation in this work is to apply this new technique in reactive transport modeling.

Conclusion et perspectives

Résumé des contributions de la thèse

Dans cette thèse nous apportons des réponses pour la modélisation des réactions de cinétiques chimiques multiphasiques. Nous proposons plusieurs types de modèle de limitation des vitesses de réaction cinétique et plusieurs manières de limiter la vitesse des concentrations des espèces chimiques. Dans le chapitre 1, nous proposons une approche par régularisation d'un modèle de surface réactive discontinue amont et un modèle limite de type inclusion différentielle. Dans le chapitre 2, nous proposons une approche discontinue d'un modèle de surface réactive amont, en utilisant la théorie de Filippov. Ces formulations sont basées sur une limitation amont des vitesses de réaction. Dans le chapitre 3, nous proposons une nouvelle approche de limitation en modulant directement la vitesse de chaque espèce à l'aide d'une dynamique projetée. Ce troisième modèle est différent des deux premiers car la limitation est appliquée espèce par espèce et non réaction par réaction, la stoechiométrie globale étant maintenue par l'ajout de contraintes de conservation.

Pour chacun des trois modèles, nous avons prouvé des résultats d'existence et de positivité des solutions. De plus, dans le chapitre 2, nous avons montré l'unicité des solutions du modèle développé pour le cas où un seul minéral disparaît et sous l'hypothèse que le champ ne soit pas tangent. Nous montrons que le modèle développé dans le chapitre 2 est une version plus précise du modèle développé dans le chapitre 1. Dans le chapitre 3, nous montrons sur des exemples choisis les différences et liens entre le modèle développé et d'autres modèles de cinétiques (modèle de type complémentarité, modèle discontinue du chapitre 1 et modèle projeté sans la contrainte sur la vitesse des réactions).

Enfin, nous proposons deux méthodes numériques de résolution des cinétiques : une basée sur le modèle régularisé du chapitre 1 et l'autre sur le modèle avec dynamique projeté du chapitre 3. Dans le chapitre 1, nous avons prouvé la convergence des solutions du modèle régularisé vers une solution d'un modèle de type inclusion différentielle. Dans le chapitre 3, nous avons montré la convergence du schéma numérique basé sur un algorithme à double projections. Nous illustrons, pour les deux méthodes, les résultats obtenus sur des cas tests unitaires variés issus de la cinétique chimique.

Travaux en cours et perspectives

Les trois premières perspectives font l'objet de travaux en cours.

Une première perspective est l'étude de la formulation de type inclusion différentielle 91 du chapitre 2 dans le cas où plusieurs minéraux disparaissent. Une étude préliminaire de classification des cas a été réalisée (configurations des champs possibles). Une étude du cas attractif en codimension 2 montre que le convexifié de Filippov ne détermine pas une valeur unique du second membre. Une sélection du second membre s'impose alors pour espérer l'unicité de la solution. Il existe plusieurs méthodes de sélection de λ dont l'une des plus connues est le "sigmoid blending" [START_REF]Sliding modes in intersecting switching surfaces, i: Blending[END_REF] and [START_REF]Sliding motion on discontinuity surfaces of high co-dimension. a construction for selecting a Filippov vector field[END_REF]. En chimie, la nature même du problème nous restreint au domaine des concentrations positives; en effet nous avons montré dans les chapitre 1 et 2 que le champs sont essentiellement positifs. Nous pouvons donc proposer une nouvelle sélection récursive du second membre adaptée aux problèmes de cinétique chimique comme combinaison convexe des champs de même codimension dans le domaine positif et du champ dans le domaine négatif (le champ du domaine pour lequel tous les minéraux ont une concentration négative). Cependant, une difficulté réside dans la démonstration de l'unicité de la solution après sélection du second membre. Une piste intéressante serait alors l'utilisation du théorème général d'unicité de Filippov [39, Theorem 4 Chapter 2 Section 10]. Plus généralement, on peut se poser la question pour un nombre quelconque de minéraux. On peut définir une version plus générale des conditions d'attractivité, mieux adaptée aux problèmes de cinétiques (voir chapitre 1 pour le cas avec un minéral), et démontrer une version générale du Théorème 2.3.2 (sliding mode).

Une deuxième perspective est d'établir le lien entre la formulation régularisée du chapitre 1 (selon le choix de la fonction de mélange pour u ) et la formulation au sens de Filippov du chapitre 2 (selon le choix du paramettre du convexifié λ). Dans le chapitre 1, nous avons montré la convergence faible-* de u vers u. Une première difficulté est d'expliciter la formule de u en fonction du choix de la fonction de mélange pour u . Des premiers travaux numériques nous permettent de conjecturer que si u ,j min iM ¡ j pH pc ,i qq alors u j min iM ¡ j pHpc i qq (de même pour le produit), avec H la fonction de Heaviside indéterminée en 0. Une idée serait d'effectuer une analyse de convergence forte afin de prouver cette conjecture. En supposant la conjecture vraie, on peut montrer que si le second membre limite du chapitre 1, f pu, v, cq, et le second membre du chapitre 2, f pλ, cq, sont égaux alors il existe une application linéaire inversible (dépendant seulement du nombre de minéraux N p ) entre H et λ, i.e. H A 2 Np λ où H est défini par H 1, Hpc 1 q, ¤ ¤ ¤ , Hpc Np q, minpHpc 1 q, Hpc 1 qq, ¤ ¤ ¤ , minpHpc Np¡1 q, Hpc Np qq, ¤ ¤ ¤ , minpHpc 1 q, ¤ ¤ ¤ , Hpc Np qq La propriété est également vraie dans le cas où H est défini avec la fonction produit. Plus généralement, cette relation nous permet de calculer des valeurs de H connaissant λ. En particulier, on retrouve le λ sigmoid blending si on fait le choix produit pour le mélange. De plus, la détermination des valeurs de λ nous permet d'obtenir les bonnes valeurs en zéro des fonctions de Heaviside Hpc i q, i 1, ¤ ¤ ¤ , N p .

Une troisième perspective serait de proposer une formulation de type complémentarité générale. Dans [START_REF] Hoffmann | Existence and uniqueness of a global solution for reactive transport with mineral precipitation-dissolution and aquatic reactions in porous media[END_REF], les auteurs proposent une formulation de type complémentarité pour des réactions contenant un unique minéral. Nous avons mené des premiers travaux et obtenu une formulation pour des réactions contenant un nombre quelconque de minéraux 92 : pour une réaction de cinétique générale, Eb 1 ¡ Eb i 0, di A ¡ M ¡ zt1u

Eb 1 Eb i 0, di A M medp min iM ¡ pc i q , ¡ min iM pc i q , Eb 1 q 0, (3.15) avec med la fonction qui renvoie la valeur médiane des trois quantitées (la deuxième valeur une fois ordonné). Cependant, considérer un nombre quelconque de réactions n'est pas direct. Une piste serait d'étudier les liens entre les formulations de types EDOs discontinues, complémentarité et dynamique projetée, afin d'obtenir une formulation de type complémentarité équivalente aux modèles du chapitre 1 ou 3. Une autre piste serait de travailler directement sur les contraintes du système.

Une quatrième perspective serait de proposer un modèle général de transport réactif couplant équilibre et cinétique chimique dans l'objectif d'effectuer de la simulation numérique. Nous avons implémenté la loi de cinétique régularisée du chapitre 1 dans le logiciel de transport réactif IFPEN Geoxim combinée avec une formulation multiphasique de type Coats existante. Pour cela, nous avons dû adapter notre modèle à la structure du logiciel et réaliser le couplage avec la partie déjà existante, transport réactif avec équilibres chimiques. Cela constitue une première implémentation globalement implicite d'un modèle de transport réactif cinétique avec un nombre de minéraux quelconque par réaction. Des premières simulations numériques sur un maillage 2D ont été effectués pour valider l'approche sur des cas unitaires présentés lors du congrès CMWR 2018 (Computational Methods in Water Resources). Un premier travail envisagé porterait sur l'amélioration et la validation de la méthode numérique développée. Un objectif supplémentaire serait de considérer des problèmes issus de la physique ou de l'industrie. Pour cela, une première étape serait de valider et d'adapter les méthodes de résolution numérique pour les cinétiques chimiques (voir chapitre 1 et 3) à des problématiques appliquées en géosciences ou dans d'autres domaines avec des processus et des données réalistes. Un second travail envisagé porterait sur une analyse théorique des propriétés du système de transport-chimie couplé, en discret en en continu. Résumé : Cette thèse concerne la modélisation mathématique des réactions cinétiques comprenant des phases pures. Dans le premier chapitre, un modèle de type EDOs discontinues pour la cinétique avec apparitions et disparitions d'espèces pour un nombre quelconque de minéraux est proposé. Une version régularisée du modèle permet de prouver la positivité et l'existence. Une analyse explicite plus approfondie dans le cas contenant une espèce réactive intermédiaire est menée. Dans le deuxième chapitre, une reformulation du modèle de cinétique chimique utilisant la théorie de Filippov est proposée. Une preuve de l'existence et de la positivité des solutions est réalisée. De plus, dans le cas des surfaces de discontinuité de codimension 1, une étude des configurations des champs fournit un résultat d'unicité et de caractérisation des trajectoires. Dans le troisième chapitre, un modèle de cinétique chimique de type systèmes dynamiques projetés est proposé. Une analyse de l'existence des solutions de ce modèle, des liens avec d'autres types de formulations et une méthode de résolution numérique adaptée sont proposés. Enfin, une illustration des résultats numériques obtenus est réalisée pour des systèmes de cinétique chimique.

Title: Mathematical modeling of multiphase kinetic reactions in geochemistry

Keywords: chemical kinetics, discontinuous ODEs, differential inclusions, projected dynamical system, Filippov's theory, numerical simulations, existence analysis.

Abstract: This thesis focuses on the modeling chemical kinetics for reactions involving pure phases. In the first chapter, a discontinuous ODEs model for reactions with appearance and disappearance of species for any number of minerals is proposed. A regularized version of the model can prove positivity and existence. An explicit analysis in the case containing an intermediate reactive species is investigated. In the second chapter, a reformulation of the chemical kinetics model using Filippov's theory is proposed. A proof of the existence and the positivity of the solutions is given. In addition, in the case of discontinuity surfaces of codimension 1, a study of the configurations of the vector fields provides a result of uniqueness and a characterization of the trajectories. In the third chapter, a model of chemical kinetics of the projected dynamical system is proposed. An analysis of the existence of solutions of this model, links with other types of formulations and an adapted numerical resolution method are provided. An illustration of the numerical results obtained is made for chemical kinetic systems.

  1¤j¤Nr et s i,j Z -Phase pure est une phase contenant une unique espèce (dans nos travaux il s'agit exclusivement des minéraux) -Multiphasique : ce dit d'un système comprenant plusieurs phases (dans nos exemples il y a toujours une phase aqueuse et une ou plusieurs phases pures minérales) -Réactifs de R j (indices) : A

0 -k e k i 5 1,Notions et notations du chapitre 3 - 9 x-

 0539 Modèle continue : avec r τ (voir (1.4)) -Limiteurs u et v (voir (1.14)) -Modèle limité (inclusion différentielle) (voir (1.15)) -Modèle régularisé (voir (1.19)) -Blending : la fonction de mélange des limiteurs pour le modèle régularisé (voir (1.18) et (1.27)) -Système réactif avec compétition : Si la consommation d'une espèce du système est due à plusieurs réactions du système -Surfaces de discontinuité : Σ it1,¤¤¤ ,Npu Σ i avec Σ i tc pc a , c p q R Na Np | c p,i 0u les hyperplans correspondant à la valeur de concentration nulle d'un minéral -Domaines : G k le plus grand ensemble ouvert caractérisé par une valeur unique du signe des concentrations des minéraux -Vecteur représentant de G k : e if dc G k , c p,i ¡ 0, 0, if dc G k , c p,i 0. -Ensemble des indices des domaines voisins à un point c : Gpcq -Modèle de type Filippov (inclusion différentielle) (voir (2.8)) -Codimension d'un point c : le nombre de composente minérale égale à zéro -Champs attractif : voir (2.15) -Champs transverse : voir (2.17) -Champs tangent : voir (2.18) Espèces chimiques : nombre d'espèces totales n aqueuses W i et minérales M i -Nombre de réactions : m -Stoechiométrie : S Z n¢m -Quantités (concentrations) : x -Espèce aqueuses (ensemble d'indice) : W -Taux de réaction : τ (voir (3.1)) -Invariants : Q base de kerpS T q (S T Q Q T S 0) Projection : ProjpC; vq : tw C : }v ¡ w} 2 inf zC }v ¡ z} 2 u -Cône tangent :(voir (3.2)) T pC; xq td : hd ν Ñ d and hτ ν Ó 0 with x τ ν d ν C dνu. -Contraintes sur les concentrations : C 0 tx R n : Q T px ¡ x 0 q 0u -Contraintes sur la vitesse des concentrations : T pC 0 ; xq Dpxq avec Dpxq : tdiagpSτpxqqλ with λ r0, 1s n u -Système avec dynamique projetée : voir (3.4)

  connaît pas leur constitution atomique, on utilise la convention suivante : W (ou A) complété d'un indice dénote une espèce en phase aqueuse (i.e. phase dans laquelle l'eau est le solvant), M complété d'un indice dénote un minéral pur et G complété d'un indice dénote une espèce gazeuse. Par conséquent, une fois ordonnée, une espèce dans une phase est entièrement caractérisée par son indice dans la phase, cela permet d'utiliser la notation i α pour référencer la i-ème espèce dans la phase α.

Figure 1 -

 1 Figure 1 -Description du système compositionnel

( 39 )F pcq 2 f

 392 Afin de définir le modèle d'inclusions différentielles permettant l'utilisation de la théorie de Filippov, on introduit les notations suivantes : Gpcq est l'ensemble des indices des domaines G k adhérent à c. Le taux de réaction multivalué utilisant la convexification de Filippov est défini comme suit rpλ, cq ķGpcq λ k r k pcq, λ Λpcq, (40) avec Λpcq tλ r0, 1s card Gpcq , ķGpcq λ k 1u. (41) Il en résulte le système d'inclusion différentielle suivant : pλ, cq, λ Λpcq @ avec f pλ, cq Srpλ, cq, le second membre multivalué au sens du convexifié de Filippov.

( 1 . 11 )

 111 The discontinuous model (1.9) with reaction rate (1.11) has a unique weak solution, satisfying definition 1.2.1. In what follows, we compute this solution, in the three situations depicted in Figure1.1.

Figure 1 . 1 -

 11 Figure 1.1 -Quantities of solute c A and mineral c M for reaction A é M with K 10. Top left: precipitation with c A,0 11, c M,0 8. Top right: partial dissolution with c A,0 7, c M,0 4. Bottom: full dissolution with c A,0 5, c M,0 4.

Lemma 1 . 4 . 1 .

 141 The function r is locally Lipschitz continuous. Proof. The functions u j, and v j, are bounded and Lipschitz continuous. The function τ is C 1 and locally bounded. Therefore the function r is locally Lipschitz continuous. Proposition 1.4.1.

Lemma 1 . 4 . 2 .

 142 The function f defined by f pcq Sr pcq is locally Lipschitz continuous and essentially non negative. Proof. By using Lemma 1.4.1 and linearity, f is locally Lipschitz continuous. Let c a vector of species quantities such that c ¥ 0.

Proposition 1 . 4 . 2 .

 142 If c 0 ¥ 0 then the solution of system (1.19) satisfies c ¥ 0 in the maximal interval. Proof. As in Proposition 1.2.3, it is a direct consequence of Lemmas 1.4.1 and 1.4.2.

Theorem 1 . 4 . 1 .

 141 Under Assumption 1.4.1, system (1.19) has a unique solution c in r0, Vr which satisfies 0 ¤ c ¤ C, where C is defined in Lemma 1.2.1.

Figure 1 .

 1 Figure 1.2 illustrates this case. Mineral M 2 precipitates at rate pτ 1 ¡ τ 2 q until M 3 is fully dissolved, then reaction R 2 stops and M 2 precipitates at rate τ 1 until M 1 is fully dissolved and reaction R 1 stops.

Figure 1 . 2 -

 12 Figure 1.2 -Reactions M 1 é M 2 and M 2 é M 3 with maximal rates τ 1 1, τ 2 ¡3 and initial conditions c 1,0 5, c 2,0 0, c 3,0 5.

Figure 1 . 3 -

 13 Figure 1.3 -Reactions M 1 é M 2 and M 2 é M 3 with maximal rates τ 1 ¡1, τ 2 3 and initial conditions c 1,0 c 3,0 0, c 2,0 10.

Figure 1 .4 where t ¦ 2 0

 12 , mineral M 1 dissolves at rate τ 1 and mineral M 2 dissolves at rate τ 2 ¡ τ 1 then reaction R 1 stops and M 2 dissolves at rate τ 2 until reaction R 2 stops.

Figure 1 . 4 -

 14 Figure 1.4 -Reactions M 1 é M 2 and M 2 é M 3 with maximal rates τ 1 2, τ 2 1 and initial conditions c 1,0 2, c 2,0 3, c 3,0 0.

Figure 1 . 5 -

 15 Figure 1.5 -Reactions M 1 é M 2 and M 2 é M 3 with maximal rates τ 1 1, τ 2 3 and initial conditions c 1,0 5, c 2,0 5, c 3,0 0.

Figure 1 .

 1 Figure 1.6 shows the regularized solution and the reaction rate computed with 0.001.

Figure 1 . 6 -

 16 Figure 1.6 -Case A ¡ M with dissolution of M . Results with 0.001.

Figure 1 . 7 -

 17 Figure 1.7 -Case A ¡ M with dissolution of M . Results and errors for different values of . Top: Quantity of mineral (left) and function u (right). Bottom: error e (left) and error e u (right). The interval rt ¦ ¡ δ, t ¦ δs is represented in pink.

Figure 1 . 8 -

 18 Figure 1.8 -Numerical convergence for the case A ¡ M with dissolution of M .

Figure 1 .

 1 Figure 1.9 shows the regularized solution and both limiters computed with 0.001.

Figure 1 . 9 -

 19 Figure 1.9 -Case M 1 ¡M 2 ¡M 3 with dissolution of M 2 then M 1 . Results with 0.001.

Figure 1 . 10 -

 110 Figure 1.10 -Case M 1 ¡ M 2 ¡ M 3 with dissolution of M 2 then M 1 . Results and errors for different values of . Top: concentrations (left) and function u (right). Bottom: error e (left) and error e u (right). Both intervals rt ¦ ¡ δ, t ¦ δs are represented in pink.

Figure 1 .

 1 Figure 1.11 represents errors in function of for c and u. In fact, Theorem 1.4.2 ensures that c converges strongly towards c in C 0 . Numerically, the rate of convergence is almost 1 for both norms.

Figure 1 . 11 -

 111 Figure 1.11 -Numerical convergence for the case M 1 ¡ M 2 ¡ M 3 in dissolution of M 2 then M 1 .

Figure 1 .Figure 1 . 12 -

 1112 Figure 1.12 shows the regularized solution and reaction rates computed with 0.001,

Figure 1 . 13 -

 113 Figure 1.13 -Case with competition with only minerals. Results with 0.001.

Figure 1 .Figure 1 . 14 -

 1114 Figure 1.14 shows the regularized solution and reaction rates computed with 0.001,

  , the surface where at least one mineral has zero concentration. The set Σ divides the space R Na Np into 2 Np open subspaces G k , k t1, ¤ ¤ ¤ , 2 Np u, uniquely characterized by the sign of the concentration of the mineral species. It is convenient to use a vector representation e k R Np of the domains G k :

Figure 2 . 1 -

 21 Figure 2.1 -Space decomposition for two minerals

Figure 2 . 4 -

 24 Figure 2.4 -Right-hand side in the attractive case at c of codimension 1

Definition 2 . 3 . 3 ( 1 1

 2331 Tangent field). A kinetic discontinuous field f is tangent at c Σ 1 if and only if f pcq 0 and f 2 1 pcq 0.

Figure 2 . 6 -

 26 Figure 2.6 -Field on a discontinuity surface in the tangent case

Sf 1 k

 1 k,j τ j pcq jt1,¤¤¤ ,Nru S 1,j ¡0 S k,j τ j pcq ¡ jt1,¤¤¤ ,Nru S 1,j 0 S k,j τ ¡ j pcq.Using Proposition 2.3.4, dj t1, ¤ ¤ ¤ , N r u such that S 1,j $ 0, τ j pcq 0. So, pcq f 2 k pcq jt1,¤¤¤ ,Nru S 1,j 0 S k,j τ j pcq.

  a ¡ βq ¡ Hpc p qpαc a ¡ βq ¡ ¨, dc p dt pαc a ¡ βq ¡ Hpc p qpαc a ¡ βq ¡ .So, the values of the discontinuous right-hand side in all domains are:-c p ¡ 0 then f 1 pcq ¢ ¡pαc a ¡ βq pαc a ¡ βq ¡ pαc a ¡ βq ¡ pαc a ¡ βq ¡ ¢ ¡pαc a ¡ βq pαc a ¡ βq , -c p 0 then f 2 pcq ¢ ¡pαc a ¡ βq pαc a ¡ βq .

Figure 2 .

 2 Figure 2.7 illustrates the field configuration described above.

Figure 2 . 7 - 1 2 pc ¦ q and f 2

 2712 Figure 2.7 -Field configuration of the reference example in the tangent case

Q

  R n¢pn¡mq Conserved quantities matrix m N Number of reactionsk p R m Products species constant S Z n¢m Stoichiometry matrix k r R m Reactants species constant x R n Concentration τ R mReaction rate

9 xptq

 9 Proj pDpxptqq T pC 0 , xptqq; Sτ pxptqqq , xp0q x 0 . (3.4)

  zptq

Figure 3 . 1 -

 31 Figure 3.1 -Solution of the case M 1 M 2 M 3 considering τ p2 6q T and x 0 p1 1 1q T .

5 9xptq

 5 Proj pT pC 0 , xptqq; Sτ pxptqqq , xp0q x 0 .(3.6) We consider the same example as above:M 1 é M 2 and M 2 é M 3 .

Figure 3 .

 3 Figure 3.2 shows the solution of (3.6).

3 Figure 3 . 2 -

 332 Figure 3.2 -Solution of the case M 1 M 2 M 3 considering τ p2 6q T and x 0 p1 1 1q T without the constraint 9

Lemma 3 . 3 . 3 .

 333 If Q satisfies (3.12) then C 0 is bounded.

1 ¡1.

 1 Assumption 3.3.1 fails to hold here as Q T d 0 for d p1 1 2q T ¡ 0.

  3 and Lemma 3.3.4, the sets C 0 and Dpz h qtd : Q

2 y 2j 1 e 2j ¡ y 2j 2 .( 3 . 14 )

 2314 The following result proves the convergence of the sequence (3.14). Theorem 3.4.1. The sequence ty j u j computed by (3.14) converges to x k 1 satisfying (3.7).

3 repeat 4 y 2j 1 5 e 6 d 7 y 8 e 9 j : j 1 ;

 341567891 ProjpC 0 ; y 2j e 2j¡1 q ; 2j 1 y 2j e 2j¡1 ¡ y 2j 1 ; j : Proj¡ Dpx k q; y 2j 1 e 2j ¡ x k h © ; 2j 2 : x k h d j ; 2j 2 y 2j 1 e 2j ¡ y 2j 2 ;10 until }y 2j ¡ y 2j¡1 } ¤ ; 11 return: x k 1 : y 2j .

Figure 3 . 3 -

 33 Figure 3.3 -Example W 1 W 2 . On the left: k r 2, k p 1 and x 0 p1 0q T ; on the right: k r 1, k p 5 and x 0 p1 1q T .

Figure 3 . 4 -Figure 3 . 4

 3434 Figure 3.4 -Example W M , with k r 1, k p 2 and x 0 p1 0.5q T .

3 Figure 3 . 5 -

 335 Figure 3.5 -Example M 1 M 2 M 3 . This example considers k r p3 7q T , k p p1 1q T and x 0 p1.5 1 0q T .

3 Figure 3 . 6 -

 336 Figure 3.6 -Example W M 1 M 2 , with k r p3 7q T , k p p1 1q T and x 0 p1.5 1 0q T .

  i ķM |s k |M k . Soit 1 M ¡ , Eb i dc i dt |s i |τ, di A ¡ M ¡ et Eb i dc i dt ¡ |s i |τ, di A M .

Titre:

  Modélisation mathématique de réactions cinétiques multiphasiques en géochimie Mot clés : cinétique chimique, EDO discontinue, inclusion différentielle, système dynamique projeté, théorie de Filippov, simulation numérique, analyse d'existence.

Table des matières

 des Introduction 0.1 Modèle de transport réactif multiphasique . . . . . . . . . . . . . . . . . 0.1.1 Description du système physique . . . . . . . . . . . . . . . . . . Etat de l'art . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 0.3 Descriptions des travaux de thèse . . . . . . . . . . . . . . . . . . . . . . Maximal reaction rates . . . . . . . . . . . . . . . . . . . . . . . . 1.2.5 Continuous model . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 State of the art: one mineral in one reaction . . . . . . . . . . . . . . . . 1.3.1 Discontinuous model . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.2 Set-valued model . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.3 Complementarity model . . . . . . . . . . . . . . . . . . . . . . . 1.3.4 Reference example . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.4.1 Solution of the discontinuous model . . . . . . . . . . . . 1.3.4.2 Precipitation . . . . . . . . . . . . . . . . . . . . . . . . 1.3.4.3 Partial dissolution . . . . . . . . . . . . . . . . . . . . . .4 A differential inclusion model for general geochemical systems . . . . . . 1.4.1 Regularized model . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4.2 Convergence of the regularized solutions . . . . . . . . . . . . . . 1.5 Test case with three minerals . . . . . . . . . . . . . . . . . . . . . . . . 1.5.1 Dissolution of M 1 and M 3 . . . . . . . . . . . . . . . . . . . . . . 1.5.2 Dissolution of M 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.5.3 Dissolution of M 1 then M 2 . . . . . . . . . . . . . . . . . . . . . . 1.5.4 Dissolution of M 2 then M 1 . . . . . . . . . . . . . . . . . . . . . . 1.5.5 Summary of test case M 1 ¡ M 2 ¡ M 3 . . . . . . . . . . . . . . . . 1.6 Numerical experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.6.1 Numerical context . . . . . . . . . . . . . . . . . . . . . . . . . . 1.6.1.1 ODE solver . . . . . . . . . . . . . . . . . . . . . . . . . 1.6.1.2 Physical model . . . . . . . . . . . . . . . . . . . . . . . 1.6.1.3 Discrete error norms . . . . . . . . . . . . . . . . . . . . 1.6.2 Reference example A ¡ M , total dissolution . . . . . . . . . . . . 1.6.3 Test case M 1 ¡ M 2 ¡ M Link with Other Models . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.1 Link with a complementarity system . . . . . . . . . . . . . . . . 3.2.2 Link with a discontinuous model . . . . . . . . . . . . . . . . . . 3.2.3 Link with a classical projected dynamical system . . . . . . . . . 3.3 Existence of a Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . Existence of a solution . . . . . . . . . . . . . . . . . . . . . . . . 3.4 An Alternating Projection Algorithm . . . . . . . . . . . . . . . . . . . . 3.5 Numerical Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5.1 Example W 1 W 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5.2 Example W M . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5.3 Example M 1 M 2 M 3 . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5.4 Example W M 1 M 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.6 Concluding Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . .

	Glossaire
	Conclusion et perspectives Bibliography Notions et notations des chapitres 1 et 2
	-Espèces chimiques :
	iii

0.1.2 Inconnues du système . . . . . . . . . . . . . . . . . . . . . . . . . 0.1.3 Lois de fermeture . . . . . . . . . . . . . . . . . . . . . . . . . . . 0.1.4 Physique de l'écoulement en milieu poreux . . . . . . . . . . . . . 0.1.5 Interactions chimiques . . . . . . . . . . . . . . . . . . . . . . . . 0.1.6 Système d'équations . . . . . . . . . . . . . . . . . . . . . . . . . 0.2 0.3.1 Chapitre 1 : Modèles de cinétique chimique avec des phases pures 0.3.2 Chapitre 2 : Analyse d'un modèle de type EDOs discontinues pour la cinétique chimique multiphasique . . . . . . . . . . . . . . . . . 0.3.3 Chapitre 3 : Modèle avec dynamique projetée pour les réactions de dissolution et de précipitation des minéraux en géochimie . . . . .

1 Kinetic models with pure phases 1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Mathematical model of kinetic reactions . . . . . . . . . . . . . . . . . . 1.2.1 Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.2 Quantities of species and reaction rates . . . . . . . . . . . . . . . 1.2.3 Conservation laws . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.4 1.3.4.4 Full dissolution . . . . . . . . . . . . . . . . . . . . . . . 13 . . . . . . . . . . . . . . . . . . . . . . . 1.6.4 Test case A 1 ¡ M 2 ¡ A 3 . . . . . . . . . . . . . . . . . . . . . . . 1.6.5 Test case with competition . . . . . . . . . . . . . . . . . . . . . . 1.6.6 Test case with blending (non-uniqueness) . . . . . . . . . . . . . . 2 Analysis of a discontinuous ODEs model for multiphase chemical kinetics 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Differential inclusion model for chemical kinetics . . . . . . . . . . . . . . 2.2.1 Domains and surfaces of discontinuity . . . . . . . . . . . . . . . . 2.2.2 Extension of reaction rates . . . . . . . . . . . . . . . . . . . . . . 2.2.3 Differential inclusion kinetic model . . . . . . . . . . . . . . . . . 2.2.4 Essential positivity property . . . . . . . . . . . . . . . . . . . . . 2.2.5 Existence of solutions . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Uniqueness of the solution when only one mineral disappears . . . . . . . 2.3.1 Attractive field . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.2 Transverse field . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.3 Tangent field . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 A projected dynamical system approach to mineral precipitation-dissolution reactions in geochemistry 3.1 A Projected Dynamical System Model . . . . . . . . . . . . . . . . . . . 3.1.1 Kinetic system . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1.2 Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 3.3.1 Explicit scheme and interpolation . . . . . . . . . . . . . . . . . . 3.3.2 Bounded constraints . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.3 Convergence of the discrete function . . . . . . . . . . . . . . . . ii 3.3.4

  où rE i s est l'activité d'une espèce, K i est la constante d'équilibre d'une espèce, k f,j est la constante de vitesse du côté réactif et k b,j est la constante de vitesse du côté produit.

	Remark 0.3.1. La formule (30) fait apparaître les dépendances de k r et k p , définies
	en (17), par rapport aux constantes d'équilibre K

  où T pC 0 ; cq est le cône tangent de C 0 en c. La projection d'un vecteur v sur un ensemble E est définie par projpE; vq : tw E : }v ¡ w} 2 inf zE }v ¡ z} 2 u.

La deuxième contribution de ce chapitre est un lien entre le modèle

[START_REF] Hoffmann | Reactive transport and mineral dissolution/precipitation in porous media : efficient solution algorithms, benchmark computations and existence of global solutions[END_REF]

et d'autres modèles : complémentarité dans le cas A é M , discontinu (33) dans le cas

Table 1 . 4

 14

	c 0	.5, 1q p6, 4q p1, 0.6, 0.3, 0.2q

-Test case with competition: physical parameters

Table 1 .

 1 6 -Test case with blending: physical parameters Figure1.15 shows numerical results obtained with eps=0.001 for these two blendings. We observe that limiters converge, but to different values.Indeed, with min function, we get numerical convergence to u 4 0.25, u 5 0.45, u 3 minpu 4 , u 5 q 0.25.

			c 0	.1, 1.4q p4, 3q p2, 1q p2, 1q p0.2, 0.1, 2, 2, 2, 2, 2q
					(1.30)
	With prod function, numerical limits of limiters are
	u 4	43 2 ¡ 3 5	0.32736, u 5	43 2 ¡ 2 5	0.52736, u 3 u 4 u 5 0.1726. (1.31)

Consequently, solutions are different since their dynamics are governed by different limiters. For example, dc 6 dt u 4 , dc 7 dt u 5 .

  introduced in Chapter 1 to negative concentrations of mineral species. Definition 2.2.2 (Extended reaction rates in G). The reaction rates r : G Þ Ñ R Nr are defined by, dc G, r j pcq min

	iM ¡ j	Hpc i q ¨τ	j pcq ¡ min iM j

  1 . The proof of Theorem 2.3.2 doesn't apply and we can't determine if the solution remains on Σ 1 or enters G 1 . Nevertheless, there exist chemical systems with tangent vector fields satisfying Hypothesis 2.3.1: Example 2.3.2 (case M 1 M 2 M 3 when M 2 disappears first in the case f 1 1 0 and f 2 1 ¡ 0). Let consider the case τ 1 ¡ τ 2 0. In this case the field is attractive because f 1 2 pcq 0 and f 2 2 pcq ¡ 0. Hypothesis 2.3.1 holds true. Then Corollary 2.3.2, in this context, yields λ 1 and the unique possible value of the right-hand side F on Σ 2 tc 2

  with k t t h u, and t $ lh, l 0, . . . , t t f Lemma 3.3.2. Let h ¡ 0 and x 0 C 0 . The discrete function x h is solution of the following discrete projected dynamical system: for all t $ lh, l 0, . . . , t t f

			h	u.	(3.10)
	The discrete function x h satisfies an approximate model:
			h	u,
	6 9 8	9 x h ptq Proj Dpx k q	2
	9 7

For every vector a R n and b R n , diagpaqb diagpbqa pa i b i q it1,¤¤¤ ,nu .

Remerciements

Using the limit values of the limiters obtained, system (1.15) has not a unique solution in this case. Indeed for t ¥ maxpt ¦ 1 , t ¦ 2 q, both choices, (1.30) and (1.31), give different solutions of (1.15) satisfying dc 1 ptq dt 0 and dc 2 ptq dt 0.

Attractive field

The attractive field configuration occurs if the field in G 1 points towards the discontinuity surface Σ 1 with at most one of both vector fields equal to 0. See Figure 2.3 for an illustration. Definition 2.3.1 (Attractive field). A kinetic discontinuous field f is attractive at c Σ 1 if and only if rf 1 1 pcq 0 pand f 2 1 pcq ¥ 0qs or rf 1 1 pcq 0 and f 2 1 pcq ¡ 0s.

(2.15)

3 -Field on a discontinuity surface in the attractive case. Left: strict attractive case. Right: attractive case with a tangent vector field f 1 .

We study the trajectory of the unique solution in the attractive case. Let introduce the following assumption. Hypothesis 2.3.1 (Attractive neighbourhood). Let c Σ 1 such that the field is attractive. There exists V an open neighbourhood of c such that dc ¦ G 1 V, f 1 1 pc ¦ q ¤ 0.

The following result ensures that in the attractive case, the solution remains on Σ 1 .

Theorem 2.3.2 (Attractive sliding motion). Let c : r0, T s Þ Ñ C 0 be a solution of (2.8).

Suppose that ht ¦ r0, T r such that cpt ¦ q Σ 1 , the attractive conditions are satisfied at cpt ¦ q and we assume that Hypothesis 2.3.1 holds true. Then there exists t I ¡ t ¦ such that for all t rt ¦ , t I s, the trajectory remains on Σ 1 . Proof. Let c : r0, T s Þ Ñ C 0 be a solution of (2.8). Let t ¦ r0, T r such that cpt ¦ q Σ 1 and the attractive conditions are satisfied at cpt ¦ q. Using Proposition 2.3.2, we know that all minerals components c p,i , di ¥ 2, are strictly positive. It remains to prove that there exists t I ¡ t ¦ such that for all t rt ¦ , t I s, c 1 ptq 0. Using Hypothesis 2.3.1, there exists V an open neighbourhood of cpt ¦ q such that dc ¦ G 1 V, f 1 1 pc ¦ q ¤ 0. Using the continuity and the non-negativity of c, there exists t I ¡ t ¦

Transverse field

The transverse field configuration is defined by the two neighbouring vector fields pointing in the same direction. See Figure 2.5 for an illustration. Definition 2.3.2 (Transverse field). A kinetic discontinuous field f is transverse at c Σ 1 if and only if

Proof. Using Proposition 2.3.1, we get f 2 1 pcq ¡ 0. By continuity, as in the proof of Proposition 2.3.3, we conclude the proof.

All possible values of the vector field F point toward G 1 . Therefore, it seems clear that any solution will leave Σ 1 to enter G 1 . The result below shows that remarkably, the choice of the vector field on the discontinuity has no impact on the solution.

Theorem 2.3.3 (Solutions leaving the surface). Let c a solution of (2.8). Let t ¦ r0, T r such that cpt ¦ q Σ 1 and f 1 1 pcpt ¦ qq ¡ 0. There exists t I st ¦ , T r such that ds st ¦ , t I s, c 1 psq ¡ 0. Proof. Let c be a solution of (2.8) and let t ¦ r0, T r such that cpt ¦ q Σ 1 and f 1 1 pcpt ¦ qq ¡ 0.

There exists t I st ¦ , T r such that ds rt ¦ , t I s, cpsq V, neighbourhood of cpt ¦ q defined in Lemma 2.3.1. Moreover, Using Theorem 2.2.1, ds rt ¦ , t I s, cpsq pG 1 Σ 1 q V. Let t st ¦ , t I s. Using Definition 2.2.5, there exists a function λ : r0, T s Þ Ñ r0, 1s, satisfying ds rt ¦ , t I s such that cpsq G 1 , λpsq 1, -

the reader to the recent survey on computational methods for geochemistry and reactive transport modeling in [START_REF] Leal | An overview of computational methods for chemical equilibrium and kinetic calculations for geochemical and reactive transport modeling[END_REF] for more insights.

We focus here on kinetically-controlled reactions involving aqueous and mineral species. The concentration evolution of each species is governed by an ideal rate, depending on the law of mass action, while satisfying non-negativity and conservation constraints. As mentioned in [START_REF] Friedly | Solute transport with multiple equilibrium-controlled or kinetically-controlled chemical reactions[END_REF], the resulting differential system might be discontinuous. This constitutes a real challenge both for the design of numerical methods and the study of the properties of the system. More sophisticated models can deal simultaneously with both kinetically controlled and equilibrium controlled reactions. We decided to focus on the former here for more simplicity. We refer to [START_REF] Hoffmann | A general reduction scheme for reactive transport in porous media[END_REF][START_REF] Erhel | Analysis of a global reactive transport model and results for the MoMas benchmark[END_REF][START_REF] De Dieuleveult | A global strategy for solving reactive transport equations[END_REF][START_REF] Amir | A global method for coupling transport with chemistry in heterogeneous porous media[END_REF] or recently to [START_REF] Erhel | Characterizations of solutions in geochemistry: existence, uniqueness, and precipitation diagram[END_REF] for more details on this important subject.

Despite the difficulty of the problem, the past decade has seen some new developments in the study of such geochemical model [START_REF] Hoffmann | A general reduction scheme for reactive transport in porous media[END_REF][START_REF] Bouillard | Diffusion with dissolution and precipitation in a porous medium: Mathematical analysis and numerical approximation of a simplified model[END_REF][START_REF] Agosti | Analysis of a model for precipitation and dissolution coupled with a darcy flux[END_REF][START_REF] Hoffmann | Existence and uniqueness of a global solution for reactive transport with mineral precipitation-dissolution and aquatic reactions in porous media[END_REF][START_REF] Kumar | Convergence analysis for a conformal discretization of a model for precipitation and dissolution in porous media[END_REF]. These new developments are following up some early discussions in [START_REF] Friedly | Solute transport with multiple equilibrium-controlled or kinetically-controlled chemical reactions[END_REF][START_REF] Knabner | An analysis of crystal dissolution fronts in flows through porous media. Part 1: Compatible boundary conditions[END_REF][START_REF] Van Duijn | Travelling wave behaviour of crystal dissolution in porous media flow[END_REF][START_REF] Van Duijn | An analysis of crystal dissolution fronts in flows through porous media part 2: incompatible boundary conditions[END_REF][START_REF] Van Duijn | Crystal dissolution and precipitation in porous media: pore scale analysis[END_REF]. In the literature, the discontinuity induced by the phase change has been handled in several ways: an Heaviside function in [START_REF] Agosti | Analysis of a model for precipitation and dissolution coupled with a darcy flux[END_REF][START_REF] Kumar | Convergence analysis for a conformal discretization of a model for precipitation and dissolution in porous media[END_REF][START_REF] Knabner | An analysis of crystal dissolution fronts in flows through porous media. Part 1: Compatible boundary conditions[END_REF][START_REF] Van Duijn | Travelling wave behaviour of crystal dissolution in porous media flow[END_REF][START_REF] Van Duijn | An analysis of crystal dissolution fronts in flows through porous media part 2: incompatible boundary conditions[END_REF][START_REF] Van Duijn | Crystal dissolution and precipitation in porous media: pore scale analysis[END_REF], a complementarity condition in [START_REF] Hoffmann | Existence and uniqueness of a global solution for reactive transport with mineral precipitation-dissolution and aquatic reactions in porous media[END_REF], and an ordinary differential equation with discontinuous reaction rate in [START_REF] Friedly | Solute transport with multiple equilibrium-controlled or kinetically-controlled chemical reactions[END_REF][START_REF] Bouillard | Diffusion with dissolution and precipitation in a porous medium: Mathematical analysis and numerical approximation of a simplified model[END_REF]. In a recent article [START_REF] Hoffmann | Existence and uniqueness of a global solution for reactive transport with mineral precipitation-dissolution and aquatic reactions in porous media[END_REF], the authors discuss the equivalence between these three approaches. In [START_REF] Hoffmann | A general reduction scheme for reactive transport in porous media[END_REF], based on a dynamical system with complementarity condition (in other words, a complementarity system), the authors study a semi-smooth Newton method to solve reactive transport problems with kinetic reactions.

A common point in the aforementioned articles is that the authors assume that there is at most one mineral per reaction. To the best of our knowledge, the only exception is [START_REF] Friedly | Solute transport with multiple equilibrium-controlled or kinetically-controlled chemical reactions[END_REF], in which the authors focus on some examples. Our main motivation in this paper is to propose a new model that allows several minerals in the reactions.

Our new model is a projected dynamical system, which is a differential equation that simultaneously deals with the dynamics and the set of constraints. The difficulty here is that the right-hand side of the differential equation is discontinuous at the boundary of the set. To tackle this difficulty, we propose a numerical scheme based on an explicit discretization of the equation and a projection-type method. We prove the convergence of this numerical method, which proves the existence of a solution constructively. We show on several examples the results obtained by an implementation of this new approach. Moreover, we show through a classical example that our new model extends the complementarity system from [START_REF] Hoffmann | Existence and uniqueness of a global solution for reactive transport with mineral precipitation-dissolution and aquatic reactions in porous media[END_REF].

The paper is organized as follows. In Section 1, we introduce the projected dynamical model for geochemistry and describe our assumptions. In Section 2, we study the connection between our new formulation and existing models used in the literature on two examples. In Section 3, we introduce a discretization of this projected dynamical system. We prove the convergence of the numerical scheme and, as a consequence, the existence of a solution of our model. In Section 4, we present an alternating projection algorithm to compute the discretized solution. In Section 5, we numerically solve several kinetic test cases. In particular, these examples illustrate that our model allows solving numerically kinetic systems with more than one mineral per reaction. Finally, we discuss some avenues for future work in Section 6.