Agnès Labande

Jean-Claude Daran

Nicholas J Long

Andrew J P White

Rinaldo Poli

Professor Didier Astruc

Rhodium(iii

published or not. The documents may come

Introduction

N-Heterocyclic carbenes (NHCs) have become widely used ligands in many areas from transition metal catalysis, 1 organocatalysis, 2 to biochemistry and medicine. 3 They are generally described as organophosphine analogues, although their -donor character is more pronounced. Unlike phosphines, they generally lead to air-stable, robust complexes which do not easily undergo ligand dissociation. However, the M-C bond stability could be detrimental to the catalytic activity. Functionalized Nheterocyclic carbene ligands combine the robust nature of the resulting metal-NHC bond with the flexibility of the secondary donor function (hemilability, stereoelectronic control on the active site) and have been successfully used in many catalytic reactions. [4][START_REF] Bierenstiel | Selected papers on S-NHC ligands[END_REF][6] However, any alteration of steric or electronic properties of a ligand implies structural modifications and thus more synthetic work.

Fine tuning of the electronic properties can also be achieved by introduction of redox-active moieties in the ligands. [START_REF] Bildstein | For ferrocene-containing NHC ligands[END_REF] This has been demonstrated in the field of small molecule sensing 8,9 and in catalysis. 10 In the latter domain, transition metal complexes containing either substitutionally inert 11 or hemilabile 12 redoxactive ligands have shown charge-dependent behaviour.

Our group has a strong interest in the chemistry and catalytic activity of functionalised N-heterocyclic carbene complexes. We have shown that NHCs bearing thioether or phosphine donors give very efficient and robust catalysts. 13,14 The work reported here describes the preparation of redox-active NHC-thioether ligands incorporating the ubiquitous ferrocene group and the study of their coordination chemistry and electrochemical behaviour.

Results and discussion

50

Synthesis

Scheme 1 Synthesis of imidazolium salts 3a and 3b.

The synthesis of the imidazolium salt 3b, i.e. the NHC-thioether precursor, was achieved in three steps from bromoferrocene 1. 55 Although the preparation of intermediate 2 and its 1,1'disubstituted analogue have been described by Mirkin et al., 15,16 we decided to follow a different synthetic route. Mirkin's reported procedures involved the monolithiation of ferrocene and afforded the products in very low yields, by reaction with S8 and 60 tosylate TsOCH2CH2Cl (compound 2, 15% yield), 15 or with disulfide (ClCH2CH2S)2 (1,1'-disubstituted analogue, 13% yield). 16 We decided to start from bromoferrocene, which could be selectively and efficiently transformed into lithioferrocene. The reaction of the latter with bis(2-chloroethyl) disulfide was 65 then carried out at the THF reflux temperature for 16 h and allowed us to dramatically improve the average yield of 2, ranging from 70 to 79%. The intermediate imidazolium salt 3a was prepared by reaction of 2 with N-mesitylimidazole in the presence of NaI to facilitate the substitution (54% yield). This compound was characterised by 1 H NMR to confirm the presence of the acidic proton, characteristic of an imidazolium salt, at 9.87 ppm. However, due to the potential hygroscopic nature of salts bearing halogen counteranions, and to the non-innocent behaviour of these anions in coordination chemistry, a metathesis with NaBF4 was carried out to afford 3b with an overall yield of 52% from 2. 1 H NMR of 3b showed the characteristic imidazolium proton signal at 8.78 ppm, shifted upfield relative to 3a.

The NHC-thioether rhodium(III) and ruthenium(II) complexes 4 and 5 were then prepared by a one-pot synthesis in the presence of Ag2O as transmetallating agent. 6,17 The reaction was carried out with 3b in order to avoid a potential mixture of Cl 20 and I on the Rh and Ru centres. It was conducted in the presence of excess Ag2O in dichloromethane at reflux. Forty hours were generally necessary to drive the reaction to completion and consume all the imidazolium salt, as confirmed by the disappearance of the proton signal at 8.78 ppm. The synthesis of silver carbene complexes from imidazolium salts containing noncoordinating counteranions is known but usually requires harsher reaction conditions, such as high temperatures, longer reaction times and the use of molecular sieves, than those generated from halide-containing imidazolium salts. 18 At least one example describes the use of such silver complexes in transmetallation reactions to prepare Rh(I) complexes. 19 Cationic complexes 4 and 5 were obtained in very good yields (4, 86%; 5, 85%) as deep red and orange solids respectively. Both complexes are air-stable and can be purified by column chromatography on neutral alumina.

35

Characterisation of complexes 4 and 5

The 1 H and 13 C NMR spectra of complex 4 confirmed the coordination of both NHC and thioether moieties. The coordination of the carbene on rhodium was evidenced by 13 C 40 NMR by a doublet at δ 166.3, with a coupling constant JRh-C of 54 Hz, typical of rhodium(III)-NHC complexes. [13][14]20 All three methyl groups from the mesityl moiety are now distinct, whereas the two ortho methyl groups had the same chemical shift in 3b. This means that the rotation of the mesityl group along the C-N 45 bond is blocked, presumably because of close contact with the Cp* ligand. The CH2S signal shifts upfield on going from 3b (δ 38.1) to 4 (δ 35.2). However, the effects of rhodium coordination were best evidenced by 1 H NMR. All the protons of the SCH2CH2Im bridge become non equivalent and appear as a 50 complex ABCD system composed of two doublets of doublets and two doublets of triplets, with the coupling pattern confirmed by a 2D-COSY 1 H-1 H experiment. This specific pattern is set by the conformational rigidity imposed by the presence of several bulky groups. Indeed, the measured coupling constants are in 55 good agreement with those predicted from the measured dihedral angles from the X-ray structure of complex 4, allowing us to assign these resonances more precisely as shown in figure 1. As in the 13 C NMR spectrum, differentiation of the two aromatic protons (slight) and ortho, ortho' methyl substituents (large) is 60 observed in the 1 H NMR spectrum.

The 1 H NMR spectrum of the ruthenium complex 5 showed a very different pattern, as all signals were broad even at low temperatures. The formation of the Ru-NHC bond, however, was confirmed by 13 C NMR with a signal at δ 164.5, similar to other 65 NHC ruthenium(II) complexes bearing a p-cymene substituent. 6 Most signals were assigned by 13 C NMR: the o,o'-CH3 substituents of the mesityl group were barely distinct, whereas both CH(CH3)2 signals from the p-cymene moiety were well separated at δ 20.4 and 24.8. The signals for the CH2Im carbon 70 appeared in the expected region at δ 49.2, as well as the CH2S signal at δ 33.7. The 13 C NMR signals for ferrocene were situated in the expected region between δ 72 and 68 but appeared very broad. The 1 H NMR spectrum analysis was more complicated. The proton signals of the mesityl and p-cymene groups were 75 assigned without ambiguity and all six methyl groups were differentiated. However, only one proton of CH2Im could be assigned at δ 4.86. The other one was not found even with the help of 2D 1 H-1 H or 1 H- 13 C correlation experiments, as were the CH2S protons. Those might be overlapping with other signals, or 80 broadened out. However, the mass spectrometry and X-ray diffraction analyses confirmed the expected structure for 5.

Complexes 4 and 5 were analysed by single crystal X-ray diffraction. Both structures consist of isolated cations and anions with no short interionic contacts. In compound 4, there is one

S1-C2-C3-N4 -50.5(2) -73.3(4)
Cg1 defines the centroid of the cycle above the metal.

coordination site, the overall coordination geometry about each metal centre is pseudo-tetrahedral or typical piano-stool geometry. However, in the case of ruthenium p-cymene complexes, it is usual to describe them as pseudo-octahedral, the  6 -ligand occupying one face of the octahedron. The other three sites in complex 5 are occupied by the S1, C5 and Cl1 atoms. Whatever the description of the coordination about the metal, the M1-S1, M1-C5 and M1-Cl1 bond lengths are comparable for the two compounds (Table 1). The only marked differences concern the angles. Indeed, the Cg1-M1-Cl1 angle is much larger for the ruthenium complex 5, 125.09(6)°, than for the rhodium one 4, 119.92(3)°. Such a difference may be explained by the bulkier p-cymene ligand with the substituted methyl roughly above the Cl position. Another difference is the larger value observed for the S1-Ru1-C5 angle -91.54( 12)° compared to 85.68( 5)° for S1-Rh1-C5. This is certainly related to the difference in conformation for the chelate ring, which is boat for the rhodium complex and half-chair for the ruthenium one. The 30 difference is reflected in the S1-C2-C3-N4 torsion angles, -50.5(2)° for 4 and -73.3(4)° for 5. This difference in conformation results in a different dihedral angle between the coordinated Cp of the ferrocene moiety and the imidazol-2ylidene ring: 61.03 (12)° for 4 and 39.9(2°) for 5.

35

The Rh1-C5 distance of 2.0550(17) Å compares well with related NHC-Cp* rhodium complexes. [20][21] The Ru1-C5 distance of 2.091(4) Å is slightly long but in the range of other NHCruthenium complexes, which are found between 2.028(4) and 2.0828( 14) Å. 6,21 The average M1-Cl1 and M1-S1 distances are 40 in agreement with other structurally characterised cationic  5 -Cp* rhodium or  6 -p-cymene ruthenium complexes reported in the Cambridge Structural Database. 22 Both ferrocene moieties have eclipsed conformation and usual geometry. 22

Electrochemistry

45

The imidazolium salts 3a, 3b and complexes 4 and 5 were analysed by cyclic voltammetry in dichloromethane. The voltammogram of 3a exhibited several waves, as the iodide anion itself is redox-active. 23 The redox potential for the ferrocene group is observed at E1/2 = +0.64 V/SCE with a chemically 50 reversible behaviour (Ep is 120 mV, however, and may indicate slow electron transfer at the surface of the electrode). The voltammogram of 3b only shows the chemically reversible wave for the ferrocene unit at +0.64 V/SCE (Ep = 111 mV), see Figure 3. The Fe II /Fe III couple in the rhodium complex 4 (Figure 4) is anodically shifted compared to the imidazolium salts with a potential of +0.77 V/SCE (Ep = 111 mV) and exhibits chemically reversible behaviour. The anodic shift signals a weak electron-withdrawing effect of coordination to rhodium and 4 indicates electronic communication between iron and rhodium through the sulfur atom. The cyclic voltammogram of the ruthenium complex 5 also presents a chemically reversible redox wave at +0.76 V/SCE (Ep = 84 mV) attributed to the Fe II /Fe III couple, see Figure 4. A very weak oxidation wave was detected at approximately +1.7 V/SCE and was tentatively assigned to the Ru II /Ru III couple. Analysis by square-wave voltammetry more clearly revealed this wave at +1.68 V/SCE (see Supporting Information), however it is far less intense than the Fe II /Fe III wave. A similar Ru(II) p-cymene complex, bearing a chelating 10 NHC-thioether ligand with methyl in place of ferrocenyl, exhibits a Ru II /Ru III redox wave at +1.55 V/SCE. 6 The anodic shift of the oxidation potential for the Ru II /Ru III couple is attributed to the stronger electron-withdrawing effect of the oxidised form of ferrocene compared to CH3. The aspect of the voltammogram also suggests that the iron oxidation has an effect on the redox behaviour of ruthenium, whose oxidation becomes irreversible. 25 The electrochemical behaviour of both complexes suggests that it should be possible to alter selectively and reversibly the oxidation state of iron by a careful choice of oxidizing and reducing

Conclusions

Two half-sandwich rhodium(III) and ruthenium(II) complexes 30 bearing a new redox-active NHC-thioether ligand have been prepared and their electrochemical properties investigated. Cyclic voltammetry experiments have shown that the ferrocenyl unit can be oxidised and reduced reversibly. The catalytic properties of these complexes and the influence of the redox state of iron on 35 the catalytic performances are under investigation and the results will be reported in due course.

Experimental General Procedures

All reactions were performed under an atmosphere of dry 2). [15][16] A solution of bromoferrocene (300 mg, 1.13 mmol, containing 92% FcBr, 7% ferrocene and 1% FcBr2) in THF (15 mL) was cooled to -78°C. n-70 BuLi (800 L, 1.24 mmol) was added dropwise and the solution was stirred for 30 min at -78°C. A solution of 2chloroethyldisulfide (260 mg, 1.36 mmol) in THF (5 mL) was added dropwise and the mixture was allowed to warm to room temperature for 2 h. The mixture was then heated at 65°C for 16 75 h. It was allowed to cool to room temperature and H2O (10 mL) was added. The phases were separated and the aqueous phase was extracted with Et2O. The organic phases were combined, washed with H2O, dried (MgSO4), filtered and concentrated in vacuo. The yellow liquid was purified by column chromatography on neutral 80 alumina (grade I). The first band contained mainly ferrocene and the desired product was obtained as a yellow liquid in the second band (250 mg, 79%). 1 H NMR (400 MHz, CDCl3)  4.34 (t, JHH = 1.8 Hz, 2H, Cp), 4.26 (t, JHH = 1.8 Hz, 2H, Cp), 4.24 (s, 5H, Cp'), 3.60-3.56 (m, 2H, CH2Cl), 2.89-2.85 (m, 2H, CH2S). A mixture of 2 (250 mg, 0.89 mmol), N-mesitylimidazole (249 mg, 1.34 mmol) and NaI (267 mg, 1.78 mmol) in degassed acetonitrile (20 mL) was heated at 80°C for 16 h. After cooling to room temperature, the solvent was evaporated and the residue taken up into CH2Cl2. The solution was filtered through Celite and the solvent evaporated to a minimum. Excess Et2O was added and the yellow precipitate was filtered on Celite, rinsed with Et2O then taken up into CH2Cl2. Evaporation of the solvent afforded 3a as a yellow solid (268 mg, 54%). 1 H NMR (400 MHz, CDCl3)  9.87 (t, JHH = 1.5 Hz, 1H, NCHN), 7.85 (t, JHH = 1.7 Hz, 1H, HC=CH Im + ), 7.20 (t, JHH = 1.8 Hz, 1H, HC=CH Im + ), 7.05 (s, 2H, CH Mes), 4.79 (pseudo t, JHH = 5.9 Hz, 2H, CH2Im + ), 4.34 (d, JHH = 1.8 Hz, 2H, Cp), 4.27 (d, JHH = 1.8 Hz, 2H, Cp), 4.20 (s, 5H, Cp'), 3.21 (pseudo t, JHH = 5.9 Hz, 2H, CH2S), 2.38 (s, 3H, p-CH3 Mes), 2.15 (s, 6H, o-CH3 Mes).

1-(-(ferrocenylthio)-ethyl)-3-(2,4,6-trimethylphenyl) imidazolium tetrafluoroborate (3b).

Compound 2 (224 mg) was treated as described previously and crude 3a was dissolved in degassed CH2Cl2 (15 mL). A solution of NaBF4 (438 mg, 3.99 mmol) in degassed H2O (10 mL) was added and the mixture was stirred vigorously for 16 h. The two phases were separated and the aqueous phase was extracted with CH2Cl2. The combined organic phases were washed with H2O, dried (Na2SO4), filtered and concentrated in vacuo to ca. 2 mL. Excess Et2O was added and the yellow precipitate was filtered through Celite, rinsed with Et2O then taken up into CH2Cl2. Evaporation of the solvent afforded 3b as a yellow solid (216 mg, 52%). 1 H NMR (400 MHz, CDCl3)  8.78 (s, 1H, NCHN), 7.73 (s, 1H, HC=CH Im + ), 7.23 (s, 1H, HC=CH Im + ), 7.03 (s, 2H, CH Mes), 4.48 (t, JHH = 5.9 Hz, 2H, CH2Im + ), 4.33 (d, JHH = 1.7 Hz, 2H, Cp), 4.25 (d, JHH = 1.7 Hz, 2H, Cp), 4.18 (s, 5H, Cp'), 3.07 (t, JHH = 5.9 Hz, 2H, CH2S), 2.36 (s, 3H, p-CH3 Mes), 2.08 (s, 6H, o-CH3 Mes). 13 55.60; H, 5.25; N, 5.41; found C, 55.48; H, 5.28; N, 5.37.

Rhodium complex (4). Imidazolium salt 3b (56.4 mg, 0.109 mmol), [RhCp*Cl2]2 (59.3 mg, 0.096 mmol) and Ag2O (14.8 mg, 0.064 mmol) were dissolved in dry CH2Cl2 (50 mL) and the Schlenk tube was protected from light. The mixture was heated at reflux for 40 h then filtered through Celite and the solvent was evaporated in vacuo. The residue was purified by column chromatography on neutral alumina (grade I, deactivated with 6-8% H2O; eluent: CH2Cl2/MeOH 95:5) to give a deep red solid (74 mg, 86%). X-Ray quality crystals were obtained by slow diffusion of cyclohexane in a saturated CH2Cl2 solution. 1 H NMR (400 MHz, CDCl3)  7.94 (d, 3 JHH = 1.8 Hz, 1H, HC=CH Im), 6.83 (d,3 JHH = 1.8 Hz,1H,HC=CH Im),2H,C6H2(CH3)3), 5.33 (dd, 2 JHH = 14.6 Hz, 3 JHH = 3.0 Hz, 1H, CHCHDIm), 4.35 (app. s, 6H, 5xC5H5 + C5H4S), 4.26 (m, 1H, C5H4S), 4.20 (m, 1H, C5H4S), 3.91 (dd, 2 JHH = 14.2 Hz, 3 JHH = 2.4 Hz, 1H, CHAHBS), 3.78 (td, 2 JHH = 3 JHH = 14.2 Hz, 3 JHH = 3.7 Hz, 1H, CHCHDIm), 3.55 (td, 2 JHH = 3 JHH = 13.5 Hz, 3 JHH = 4.2 Hz, 1H, CHAHBS), 3.34 (app. s, 1H, C5H4S), 2.29 (s, 3H, C6H2(CH3)3), 1.98 (s, 3H, C6H2(CH3)3), 1.76 (s, 3H, C6H2(CH3)3), 1.63 (s, 15H, C5(CH3)5). 13 C NMR (100 Hz, CDCl3)  166.29 (d, JCRh = 54.3 Hz, NCN), 139.10, 138.61, 135.44, 133.57 (4x Cquat C6H2(CH3)3), 129.16, 127.32 (2xCH C6H2(CH3)3), 126.08, 124.75 (2x HC=CH Im), 100.03 (d, JCRh = 6.2 Hz, (C5(CH3)5), 80.35 (Cquat C5H4S), 72.92 (CH C5H4S), 70.21 (5xCH C5H5), 69.87, 69.02, 66.50 (3xCH C5H4S), 50.03 (CH2Im), 35.17 (CH2S), 21.11, 20.08, 18.74 (3xC6H2(CH3)3), 9.47 (C5(CH3)5). MS (ESI + ) m/z 703.5 (M + -BF4, 100%). Elem. Anal. C34H41BClF4FeN2RhS (790.50): calc. C, 51.61; H, 5.23; N, 70 3.54; found C, 51.51; H, 5.16; N, 3.46. Ruthenium complex (5). Imidazolium salt 3b (52.0 mg, 0.1 mmol), [Ru(p-cymene)Cl2]2 (54.0 mg, 0.088 mmol) and Ag2O (13.6 mg, 0.059 mmol) were dissolved in dry CH2Cl2 (30 mL) and the Schlenk tube was protected from light. The mixture was 75 heated at reflux for 40 h then filtered through Celite and the solvent was evaporated to ca. 2 mL. Excess hexane was added and the precipitate was filtered through Celite, rinsed with hexane and recovered by dissolving in CH2Cl2. The solid was further purified by column chromatography on neutral alumina (grade I, ), 141.02, 137.90, 137.45, 135.45 (4x Cquat C6H2(CH3)3), 129.74, 129.67 (2xCH C6H2(CH3)3), 126.63, 122.67 (2x HC=CH Im), 107.77, 103.52 (2xCquat p-cym), 97.68, 93.16, 88.51, 79.78 (4xCH p-cym)

X-ray crystallography.

A single crystal of each compound was 100 mounted under inert perfluoropolyether at the tip of a cryoloop and cooled in the cryostream of an Agilent Technologies OD XCALIBUR 3 CCD diffractometer for 4 or a GEMINI EOS CCD diffractometer for 5. Data were collected using the monochromatic MoK radiation (= 0.71073). The structures 105 were solved by direct methods (SIR97) 27 and refined by leastsquares procedures on F 2 using SHELXL-97. 28 All H atoms attached to carbon were introduced in idealised positions and treated as riding on their parent atoms in the calculations. In compound 4, the BF4 was found to be disordered, and three 110 orientations of ca. 81, 10 and 9% occupancy were identified. The geometries of all three orientations were optimised, the thermal parameters of adjacent equivalent atoms were restrained to be similar, and only the non-hydrogen atoms of the major occupancy orientation were refined anisotropically (the others were refined 115 isotropically). The drawing of the molecules was realised with the help of ORTEP3. 29 Crystal data and refinement parameters: see footntotes. Selected bond distances and angles are given in Table 1. Crystallographic data (excluding structure factors) have been deposited with the Cambridge Crystallographic Data Centre as 120 supplementary publication no. CCDC 816175-816176. Copies of the data can be obtained free of charge on application to the Director, CCDC, 12 Union Road, Cambridge CB2 1EZ, UK (fax: (+44) 1223-336-033; e-mail: deposit@ccdc.cam.ac.uk).

Notes and references

Scheme 2

 2 Scheme 2 Synthesis of rhodium(III) and ruthenium(II) complexes.

Figure 1 4 :

 14 Figure 1 Assignment of the methylene protons in the 1 H NMR spectrum of 4: ABCD system, simulated (top) and experimental (bottom).

Figure 2 . 5  6 -Figure 2

 2562 Figure 2. The structures are closely related, they both contain the same thioether-NHC ligand chelating the metal to form a sixmembered ring and a Cl atom but differ by the -ligand: Cp* for

  24

55

 55

Figure 3

 3 Figure 3 Cyclic voltammograms of imidazolium salts 3a (dashed) and 3b (plain), 1 mM in CH2Cl2 at a scan rate of 0.2 V/s.

Figure 4

 4 Figure 4 Cyclic voltammogram of complexes 4 and 5, 1 mM in CH2Cl2 at a scan rate of 0.2 V/s.

  45nitrogen or argon using standard Schlenk techniques, unless otherwise stated. Solvents were dried and distilled by standard techniques prior to use. NMR spectra were recorded on Bruker AV300, AV400 or AV500 spectrometers. NMR chemical shifts were determined by reference to the residual 1 H or13 C solvent 50 peaks in the deuterated solvents. The NMR simulation experiment was done with SpinWorks 3. Electrospray (ES) mass spectra were recorded either at Imperial College on a Micromass LCT Premier spectrometer or at the Université Paul Sabatier by the Service Commun de Spectrométrie de Masse on a MS/MS 55 API-365 (Perkin Elmer Sciex). Elemental analyses were performed either by Mr. Stephen Boyer at London Metropolitan University or by the Service d'Analyses du Laboratoire de Chimie de Coordination in Toulouse. Cyclic voltammetry (CV) and square-wave voltammetry (SQW) experiments were carried 60 out with a PGSTAT-302N potentiostat (Metrohm) at 20°C using a Pt disk working electrode, a Pt wire counter electrode and a saturated calomel reference electrode (SCE); a 0.1 M n-Bu4BF4 solution was used as the supporting electrolyte. All electrochemical data are referenced versus SCE. Bromoferrocene 65 1 was prepared by a previously described procedure.

  C NMR (100 Hz, CDCl3)  141.36 (Cquat Mes), 136.94 (NCHN), 134.49 (2xCquat Mes), 130.64 (Cquat Mes), 129.83 (2xCH Mes), 123.70 (C=C Im + ), 123.52 (C=C Im + ), 73.94 (2xCH Cp), 69.99 (2xCH Cp), 69.61 (5xCH Cp'), 49.29 (CH2-Im + ), 38.07 (CH2S), 21.14 (p-CH3 Mes), 17.32 (2x o-CH3 Mes). Cquat Cp not found. 19 F NMR (282 MHz, CDCl3)  -151.6 (BF4 -). MS (ESI + ) m/z 431 (M + , 100%). Elem. Anal. C24H27BF4FeN2S (518.03): calc. C,

  , 71.08 (C5H5), 71.5-68.0 (C5H4), 49.20 (CH2Im), 95 33.69 (CH2S), 29.51 ((CH3)2CH p-cym), 24.77 ((CH3)2CH pcym), 21.53 (C6H2(CH3)3), 20.41 ((CH3)2CH p-cym), 18.63, 18.62 (C6H2(CH3)3), 17.95 (CH3 p-cym). MS (ESI + ) m/z 701 (M + -BF4, 100%).

6 a

 6 CNRS ; LCC (Laboratoire de Chimie de Coordination) ; 205, route de Narbonne, F-31077 Toulouse, France, and Université de Toulouse ; UPS,INPT ; LCC ; F-31077 Toulouse, France. Fax: +33 561553003; Tel: +33 561333158; E-mail: agnes.labande@☺lcc-toulouse.fr, corresponding author. b Department of Chemistry, Imperial College London, London, SW7 2AZ, United Kingdom. Fax: +44 2075945804; Tel: +44 2075945781; E-mail: n.long@imperial.ac.uk. c Institut Universitaire de France, 103, bd Saint-Michel, 75005 Paris, France. † Electronic Supplementary Information (ESI) available: Square-wave voltammogram of complex 5, crystal data and structure refinement and selected bonds and angles for complexes 4 and 5. See DOI: 10.1039/b000000x/ ‡ Crystal data for 4: [C34H41ClFeN2RhS](BF4)•CH2Cl2, M = 875.69, triclinic, (no. 2), a = 10.6786(2), b = 11.8064(3), c = 15.8586(3) Å, α = 96.8086(18), β = 102.4583(19), γ = 106.009(2)°, V = 1842.52(7) Å 3 , Z = 2, Dc = 1.578 g cm -3 , μ(Mo-Kα) = 1.164 mm -1 , T = 173 K, orange tablets, Oxford Diffraction Xcalibur 3 diffractometer; 12107 independent measured reflections (Rint = 0.0239), F 2 refinement, R1(obs) = 0.0331, 20 wR2(all) = 0.0856, 10327 independent observed absorption-corrected reflections [|Fo| > 4σ(|Fo|), 2θmax = 66°], 484 parameters. CCDC 816175. Crystal data for 5: [C34H40ClFeN2RuS](BF4)•CH2Cl2, M = 787.92, monoclinic, P21/n (no. 14), a = 11.8820(3), b = 10.2423(2), c = 27.5100(7) Å, β = 94.493(2), V = 3337.65(14) Å 3 , Z = 4, Dc = 1.568 g cm -3 , μ(Mo-Kα) = 1.080 mm -1 , T = 180 K, orange tablets, Oxford Diffraction Gemini Eos diffractometer; 7574 independent measured reflections (Rint = 0.063), F 2 refinement, R1(obs) = 0.0512, wR2(all) = 0.1073, 5590 independent observed absorption-corrected reflections [|Fo| > 4σ(|Fo|), 2θmax = 55°], 412 parameters. CCDC 816176.

Table 1

 1

	M1-Cg1	1.8334(8)	1.7460(15)
	M1-Cl1	2.3785(4)	2.3944(9)
	M1-S1	2.3622(4)	2.3501(10)
	M1-C5	2.0550(17)	2.091(4)
	S1-C2	1.822(2)	1.798(4)
	S1-C18	1.7631(9)	1.774(4)
	C2-C3	1.524(3)	1.510(6)
	C3-N4	1.459(2)	1.454(6)
	C7-C8	1.341(3)	1.331(7)
	Cg1-M1-Cl1	119.92(3)	125.09(6)
	Cg1-M1-S1	123.28(3)	119.72(6)
	Cg1-M1-C5	133.87(6)	134.56(12)
	Cl1-M1-S1	90.407(16)	91.66(4)
	Cl1-M1-C5	91.64(5)	81.79(11)
	S1-M1-C5	85.68(5)	91.54(12)

Selected bond lengths (Å), angles (°) and torsion angle (°) for complexes 4 and 5 4 (M = Rh) 5 (M = Ru)

Institut National Polytechnique for a grant to R.P. and A.L. (Soutien à la Mobilité Internationale).

Acknowledgements

The authors thank the CNRS, the Agence Nationale de la Recherche for a grant to A.L. (ANR-07-JCJC-0041) and the 1 For selected books and recent reviews, see: S. Díez-González, N-Heterocyclic Carbenes : From Laboratory Curiosities to Efficient Synthetic Tools, RSC Publishing, Cambridge, UK, 2011 ; F.Glorius, N-Heterocyclic Carbenes in Transition Metal Catalysis, Springer-Verlag, Heidelberg, Germany, 2007 ; V. César, S. Bellemin-

For a comprehensive review on chemical redox agents, see: N. G. Connelly and W. E. Geiger, Chem. Rev., 1996, 96


Bastien Hamlat

Jocelyne Erhel

Anthony Michel

Thibault Faney

Analysis and numerical computation of geochemical systems with kinetic precipitation and dissolution reactions involving several minerals

In this paper, we study underground media where chemical interactions involve aqueous and mineral species. In most kinetic models, a mineral species is involved in a single precipitation or dissolution reaction. In this work, we consider more general precipitation and dissolution kinetic reactions, which can include several minerals, which can each participate in several reactions. We propose a differential inclusion model, with reaction rates defined by set-valued functions, in order to ensure positive quantities of species. We also develop a regularized model, with reaction rates defined by regularized Heaviside functions and blendings if several minerals dissolve in the same reaction. We prove that the unique solution of this system of ODEs converges towards a solution of the differential inclusion model. We run numerical experiments by using a classical ODEs solver for computing the regularized solution. We get good numerical approximations of the quantities of species and the discontinuous reaction rates.

Introduction

Underground media provide both natural resources and potential sites for waste storage. The conflicting resource and storage objectives reveal a wide panel of chemical and biological interactions. Evolution is characterized by many complex processes, such as precipitation and dissolution with or without kinetics and ion exchanges involving evolution of pH and redox conditions. Risk assessment and resource management rely on numerical simulations at large scale. Relevant and efficient reactive transport models are essential to assess the complexity of physical and chemical processes at large time and space scales [START_REF] De Dieuleveult | A global strategy for solving reactive transport equations[END_REF][START_REF] Erhel | Solving partial differential algebraic equations and reactive transport models[END_REF].

In this paper, we consider geochemical reactions involving minerals [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], [START_REF] Leal | An overview of computational methods for chemical equilibrium and kinetic calculations for geochemical and reactive transport modeling[END_REF] and [START_REF] Van Duijn | Crystal dissolution and precipitation in porous media: pore scale analysis[END_REF]. If reactions are sufficiently fast, they can be considered in a local equilibrium, mathematically modeled by algebraic nonlinear equations [START_REF] Erhel | Characterizations of solutions in geochemistry: existence, uniqueness, and precipitation diagram[END_REF]. On the other hand, kinetic reactions are governed by mass balance equations, written as ordinary differential equations (ODEs), where the right-hand side involves stoichiometry of the geochemical system and reaction rates. These ODEs must take into account precipitation or dissolution of minerals. When a mineral participates in only one reaction, with only aqueous species, the kinetic rate can be described by a discontinuous function, a setvalued function, or a complementarity problem, [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]. These three models are equivalent under some conditions [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 some complex chemical systems, several minerals participate in the same kinetic reaction, or a mineral is involved in several reactions, [START_REF] Angeli | A tutorial on chemical reaction network dynamics[END_REF], [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]. Reaction rates must be limited as previously when dissolving minerals disappear. However, the extension of the single mineral limiting process is far from straightforward when considering a chemical network of several mineral reactions. It becomes even more intricate if several minerals dissolve in the same reaction, but are still produced in other reactions.

In this paper, we propose a generalized model dealing with these cases. A key idea of this paper is to control reaction rates by set-valued functions, which limit the rate when dissolving minerals disappear. This results in a differential inclusion problem. In the analysis of differential inclusion systems, a first difficulty is to prove the existence of a solution, while a second difficulty is to provide a numerical method to solve the system. For geochemical systems with one mineral per reaction and one reaction per mineral, Filippov's theory [START_REF] Filippov | Differential Equations with Discontinuous Righthand Sides[END_REF] provides a strong theoretical basis for the study of such systems and allows the use of efficient custom even driven methods to solve the resulting system of ODEs [START_REF] Agosti | Analysis of a model for precipitation and dissolution coupled with a darcy flux[END_REF]. Unfortunately, characterizing the set-valued functions for the systems considered in this paper is still an open question [START_REF] Dieci | Regularizing piecewise smooth differential systems: Co-dimension 2 discontinuity surface[END_REF] [10] [START_REF] Dieci | Piecewise smooth systems near a co-dimension 2 discontinuity manifold: Can one say what should happen?[END_REF]. It is also possible to limit the derivatives of each mineral species in a projected dynamical system, which is another topic under consideration [START_REF] Erhel | A projected dynamical system approach to mineral precipitation-dissolution reactions in geochemistry[END_REF].

In this work, we overcome these difficulties using regularized reaction rates to obtain a regular system of ODEs. A second key idea is the choice of an appropriate blending function to uniquely define the regularized limiters. This leads to a unique solution of the regularized system for a given blending function.Then, we prove convergence of the regularized quantities of species and limiters towards a solution of the differential inclusion problem. Since the regularized system is easy to solve, it provides a numerical approximation of the desired result.

This paper is organized as follows. Section 2 is a general introduction to geochemical kinetics with aqueous and mineral species. Section 3 summarizes the state of the art in the case of one mineral per reaction and one reaction per mineral. In section 4, we discuss the discontinuous model applied to a single reaction between a solute and a mineral. Then section 5 is devoted to our general differential inclusion problem, with the corresponding regularized model. We prove convergence and existence results. In order to illustrate our model, we consider a system of two reactions where a mineral is a product of the first reaction and a reactant of the second one. A detailed analysis of convergence is given in section 6. Numerical experiments, developed in Matlab, show the efficiency of regularization to get an approximate solution. In section 7, we propose several test cases, with either simple or complex geochemical systems. We conclude and discuss future work in last section.

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 > 0 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 (Na+Np)×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]. In reaction R j , we denote by A - i and M - k the sets of aqueous and mineral reactants, by A + i and M + k the sets of aqueous and mineral products.

A reaction R j is written

i∈A - j |s ij |A i + k∈M - j |s kj |M k i∈A + j |s ij |A i + k∈M + j |s kj |M k . (1)
In this paper, we only consider balanced chemical equations, with as many atoms on the reactants side than on the products side.

Examples of reactions with more than one mineral

We provide here three examples of geochemical reactions, where a mineral precipitates into another one.

In the first example, calcite CaCO 3 is decarbonated to produce calcium oxide CaO which undergoes hydrolysis to form calcium hydroxide Ca(OH) 2 :

R 1 : CaCO 3 CaO + CO 2 , R 2 : CaO + H 2 O Ca(OH) 2 .
In the second example, calcite either reacts with Magnesium ions to form dolomite CaM g(CO 3 ) 2 or with Magnesium and Iron ions to form siderite F eCO 3 :

R 1 : 2CaCO 3 + M g ++ CaM g(CO 3 ) 2 + Ca ++ , R 2 : 3CaCO 3 + M g ++ + F e ++ CaM g(CO 3 ) 2 + F eCO 3 + 2Ca ++ .
The last example is the phase transition of calcite into aragonite:

R 1 : (CaCO 3 ) α (CaCO 3 ) β
where α and β represent respectively the calcite and aragonite crystal forms of calcium carbonate.

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 = (c i ), i = 1, . . . , N a + N p of the time variable t with t ≥ 0. We denote by c a = (c i ), i = 1, . . . , N a the vector of quantities for aqueous species, and by c p = (c k ), k = N a + 1, . . . , N a + N p the vector of mineral quantities. The quantities c should satisfy a positivity constraint c(t) ≥ 0, t ≥ 0. We also require that c ∈ (H 1 (0, T )) 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 = (r j ), functions in L 2 (0, T ), 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:

     dc(t) dt = Sr(t), t ∈ (0, T ) c(t) ≥ 0, c(0) = c 0 , (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 (2).

Definition 2.1. Let T > 0, let c ∈ H 1 (0, T ) Na+Np . It is a weak solution of system (2) if c ≥ 0 and ∀t ∈ [0, T ], c i (t) = c 0,i + t 0 (Sr) i (s)ds, i = 1, . . . , N a + N p .
It can be noted that this system implies invariants, which correspond to balance equations and conservative variables. Indeed, let Q be a basis of it, or equivalently a basis of the subspace orthogonal to the image of S. We get Q T S = S T Q = 0. Then, any solution of (2) satisfies ∀t ≥ 0, Q T c(t) = Q T c 0 . Functions Q T c are thus constant in time and are conservative variables.

In particular, let q be the vector of molar mass of species, which satisfies q > 0. Since we assume that each equation is mass balanced, q satisfies S T q = 0, so that total mass q T c is an invariant of the system. This vector provides an upper bound on the solutions of (2) using the following lemma: Lemma 2.1. Let q be the vector of molar mass of species involved in the chemical system.

If c satisfies c ≥ 0 and q T c = q T c 0 then c ≤ C with C = 1 min k (q k ) q T c 0 . Proof. Let k ∈ {1, . . . , N a + N p }. 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.
Other notable invariants in usual chemical systems can be the charge balance or the total mass of each chemical element.

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 net rate is the difference between both. 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 2.1. The activity a i of a solute A i is a C 1 function of c a which satisfies:

c a ≥ 0 ⇒ a(c a ) ≥ 0 c a,i ≤ 0 ⇒ a i (c a ) = 0.
The forward reaction rate of reaction R j is given by

τ f,j = k f,j i∈A - j ( a i K i ) |s ij | k∈M - j ( 1 K k ) |s kj |
and its backward reaction rate is

τ b,j = k b,j i∈A + j ( a i K i ) |s ij | k∈M + j ( 1 K k ) |s kj | If A - j is empty, the product i∈A - j ( a i K i ) |s ij | is equal to 1.
The same property applies to other subsets. The net 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 (c a ) = α j i∈A - j a i (c a ) |s ij | -β j i∈A + j a i (c a ) |s ij | , j = 1, . . . , N r . (3)
where α j and β j are strictly positive constants. Thermodynamic equilibrium of reaction R j is obtained when forward and backward reactions are balanced, in other words when the net reaction rate is zero: τ j = 0. It corresponds to the mass action law for chemistry at equilibrium.

In our differential inclusion model, we control reaction rates r j with τ j being an upper bound (in absolute value). From now on, we call the vector τ the maximal reaction rates.

Continuous model

We consider the continuous differential system with the vector τ of maximal reaction rates: dc(t) dt = Sτ (c a (t)), t ≥ 0,

c(0) = c 0 . (4)
Proposition 2.1. The geochemical system (4) has a unique solution satisfying c ∈ (H 1 (0, +∞)) Na+Np .

Proof. Thanks to Assumption 2.1, the maximal rate τ is C 1 . Since the right-hand side is regular, we can apply Cauchy-Lipschitz theorem so that the system has a unique solution in (0, +∞).

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 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 (c) ≥ 0.
Since the reaction rates τ j depend only on the quantities c a of aqueous species, we can define a function f from R Na to R Na by

f k (c a ) = (Sτ (c a )) k , k = 1, . . . , N a . (5)
Proposition 2.2. The function f defined by ( 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 a(c a ) ≥ 0 and a k (c k ) = 0. If s kj > 0 then A k ∈ A + j and τ j (c a ) ≥ 0 thus s kj τ j ≥ 0. If s kj < 0 then A k ∈ A - j and τ j (c a ) ≤ 0 thus s kj τ j ≥ 0. Finally, f k (c a ) = j: s kj =0 s kj τ j (c a ) ≥ 0.

Thanks to this property, we prove that the quantities of all aqueous species remain positive.

Proposition 2.3. If c 0 ≥ 0 then the unique solution c of system (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 2.2. By using Theorem 1 in [START_REF] Chellaboina | Modeling and analysis of mass-action kinetics[END_REF], we get the result wanted.

Denoting by c the solution of (4), it may exist a mineral M k , N a + 1 ≤ k ≤ N a + N p and a time t * such that c k (t * ) < 0. Mathematical models aim at limiting reaction rates r so that the solution of (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 and dissolution of calcite: Ca ++ + CO -- 3 CaCO 3 . 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 ( 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]. If τ (c a ) > 0 then τ (c a ) = τ + (c a ) and the product mineral M precipitates. There is no need to limit the reaction. On the other hand, if τ (c a ) < 0 then τ (c a ) = τ -(c a ) and the product mineral M dissolves. The reaction rate r must be limited in order to ensure the positivity of c p : it is defined by multiplying the negative part of τ by a control function.

Discontinuous model

r(c) = τ + (c a ) -H(c M )τ -(c a ), (6)
where H is the Heaviside function defined in [0, +∞[ by

if x > 0 then H(x) = 1, if x = 0 then H(x) = 0.
It can be noted that the Heaviside function H thus the reaction rate r can be defined for negative values: H(x) = 0 if x < 0. Clearly, the reaction stops (r = 0) when M is fully dissolved and continues at the maximal rate (r = τ ) otherwise. In other words, the choice (6) satisfies the following requirements:

if c M > 0 or τ ≥ 0 then r = τ, if c M = 0 and τ ≤ 0 then r = 0. ( 7
)
The kinetic model ( 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 (c M ) with

W (u) = {1} if u > 0, = [0, 1] if u = 0, = {0} if u < 0. (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]. Therefore 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 .
It can be noted that it is more restrictive than general reactions (1) considered in this paper.

Let R j be a reaction with mineral M k , with τ j defined by [START_REF] Angeli | A tutorial on chemical reaction network dynamics[END_REF]. Then the set-valued model for reaction rate r j is written

r j = α j i∈A - j a i (c a ) |s ij | -w k β j i∈A + j a i (c a ) |s ij | ,
where w k ∈ W (c k ) with the set-valued function W defined by [START_REF] Dieci | Piecewise smooth systems near a co-dimension 2 discontinuity manifold: Can one say what should happen?[END_REF].

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, but the function r is defined by the following complementarity problem:

c M (r -τ ) = 0, c M ≥ 0, r ≥ τ.
Then system (2) with r defined above 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].

Discontinuous solution of system A M

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 (2) is then written, with c = (c A , c M ):

             dc A (t) dt = -r(t), dc M (t) dt = +r(t), c A (t) ≥ 0, c M (t) ≥ 0, c A (0) = c A,0 , c M (0) = c M,0 . (9)
The conservative variable is (c

A + c M ), since Q = (1, 1) T .
Here, we choose a simple activity model with a(c A ) = c A , so that the mass action law (3) becomes τ (c A ) = αc A -β, which we rewrite (with α = 1 and

β = K) τ (c A ) = c A -K,
where K is a strictly positive constant.

The continuous problem ( 4) with the reaction rate τ (c A ) has a unique solution:

c * A (t) = K + τ 0 exp(-t), t ≥ 0, c * M (t) = c M,0 + τ 0 (1 -exp(-t)), t ≥ 0. ( 10
)
where τ 0 = c A,0 -K. Thus τ (t) = τ 0 exp(-t) and c * A (t) > 0, t > 0. If τ 0 = 0, then τ (t) = 0 and r(t) = 0, thus the system is at equilibrium at any time. From now on, we assume that τ 0 = 0.

The discontinuous reaction rate (6) becomes here:

r(c A , c M ) = τ + (c A ) -H(c M )τ -(c A ). ( 11
)
The discontinuous model [START_REF] Dieci | Regularizing piecewise smooth differential systems: Co-dimension 2 discontinuity surface[END_REF] with reaction rate (11) has a unique weak solution, satisfying definition 2.1. In what follows, we compute this solution, in the three situations depicted in Figure 1.

Precipitation

If τ 0 > 0, then ∀t > 0, c * A (t) > 0, c * M (t) > 0, τ (t) > 0.
The solute precipitates but does not disappear and the quantity of mineral increases. Here, c * is the weak solution of system (9) and formula [START_REF] Erhel | A projected dynamical system approach to mineral precipitation-dissolution reactions in geochemistry[END_REF] gives r = τ.

Partial dissolution

If -c M,0 ≤ τ 0 < 0, then ∀t > 0, c * A (t) > 0, c * M (t) > 0.
The mineral dissolves but does not disappear and the quantity of solute increases. Here too, c * is the weak solution of system (9) and r = τ.

Full dissolution

If τ 0 < -c M,0 , then c * M (t M ) = 0 where t M ≥ 0 is given by t M = log( τ 0 τ 0 + c M,0 ) ( 12
)
The mineral is fully dissolved. Let (c A , c M ) be defined by:

       c A (t) = c * A (t), t ≤ t M , c M (t) = c * M (t), t ≤ t M , c A (t) = c A,0 + c M,0 , t ≥ t M , c M (t) = 0, t ≥ t M ,
It is the weak solution of ( 9) with r = H(c M )τ .

A differential inclusion model for general geochemical systems

In this section, we consider the general geochemical system described in section 2. We develop an original new model based on a differential inclusion [START_REF] Filippov | Differential Equations with Discontinuous Righthand Sides[END_REF]. In this model, each reaction rate r j is equal to τ j if the quantities of dissolving minerals are strictly positive, and must be controlled if one mineral is fully dissolved. It is important to note that 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 . A dissolving mineral M k in reaction R j is such that s kj τ j < 0. For example, if τ j > 0, then τ + j = τ j and τ - j = 0; reactant minerals dissolve, so that the limiter takes into account reactant minerals. On the other hand, if τ j < 0, the limiter is based on dissolving product minerals. As in the set-valued model, limiters belong to set-valued functions.

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:

∀j = 1, . . . , N r , |r j | ≤ |τ j | and r j τ j ≥ 0. ( 13
)
For any j ∈ {1, . . . , N r }, we introduce set-valued functions U j and V j of variables c p ∈ (R + ) Np . Then we introduce two control functions u j , v j ∈ L ∞ (0, T ) in order to define the reaction rates.

Definition 5.1. The set-valued functions U j , j = 1, . . . , N r are defined in [0, +∞) Np by:

U j (c p ) = {1} if ∀k ∈ M - j , c k > 0, or M - j = ∅, = [0, 1] else.
The set-valued functions V j , j = 1, . . . , N r are defined in [0, +∞) Np by:

V j (c p ) = {1} if ∀k ∈ M + j , c k > 0, or M + j = ∅, = [0, 1] else.
For general geochemical systems, we can now define our model of reaction rates. Definition 5.2. Let T > 0 and c ∈ (H 1 (0, T )) Na+Np such that c(t) ≥ 0, 0 ≤ t ≤ T. Reaction rates r are defined by

r j = u j τ + j (c) -v j τ - j (c), j = 1, . . . , N r , (14)
where u j , v j ∈ L ∞ (0, T ) and u j (t) ∈ U j (c p (t)), v j (t) ∈ V j (c p (t)), 0 < t < T.

Clearly, the vector of reaction rates r satisfies [START_REF] Erhel | Solving partial differential algebraic equations and reactive transport models[END_REF]. 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] Filippov | Differential Equations with Discontinuous Righthand Sides[END_REF]:

find c ∈ (H 1 (0, T )) Na+Np , u ∈ (L ∞ (0, T )) Nr , v ∈ (L ∞ (0, T )) Nr such that              dc(t) dt = Sr(t), t ∈ (0, T ), c(t) ≥ 0, t ∈ (0, T ), c(0) = c 0 , u j (t) ∈ U j (c p (t)), v j (t) ∈ V j (c p (t)), t ∈ (0, T ), r j = u j τ + j (c a ) -v j τ - j (c a ), j = 1, . . . , N r , (15)
where c 0 ≥ 0 is a given initial condition, τ is defined by (3), U j and V j are introduced in Definition 5.1.

We are now faced to the two issues stated in the introduction: prove existence of a solution and compute it. Here, 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].

Regularized model

Our 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 (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.

Let us first define a regularized Heaviside function H :

   if x ≥ then H (x) = 1, if 0 ≤ x ≤ then H (x) = γ(x/ ), if x = 0 then H (x) = 0, ( 16
)
where γ is an increasing Lipschitz continuous function defined in the interval [0, 1] such that γ(0) = 0 and γ(1) = 1, so that the function H is Lipschitz continuous and bounded.

For example, we can choose γ(x) = x.

We can now define regularized reaction rates. Let us first consider a reaction, where only one mineral M k participates. If it is a reactant, then u j, (c) = H (c k ) and v j, (c) = 1. On the other hand, if it is a product, then u j, (c) = 1 and v j, (c) = H (c k ). Then either

r j = τ j or r j = H (c k )τ j .
It becomes more difficult in a general case with a set of reacting minerals, since we have to deal with several potential dissolutions. n this case, we define the reaction limiter by using a blending of the regularized Heaviside functions. Let us consider a reaction R j with τ j > 0 such that reactant minerals dissolve. We limit τ + j by a blending function, as introduced in [START_REF] Alexander | Sliding modes in intersecting switching surfaces, i: Blending[END_REF], which is equal to 0 as soon as one of the reactant minerals is fully dissolved and is equal to 1 as long as all of them are present in sufficient quantity. Similarly, we limit τ - j by a blending function applied to product minerals. Examples of such blendings are the minimum and the product functions. For the sake of simplicity, we choose the min function in what follows.

Definition 5.3. For general geochemical systems, the regularized reaction rates are defined by r j, (c) = u j, (c)τ + j (c) -v j, (c)τ - j (c), [START_REF] Knabner | An analysis of crystal dissolution fronts in flows through porous media. Part 1: Compatible boundary conditions[END_REF] where u j, and v j, are defined by

     if M - j = ∅, then u j, (c) = 1, else u j, (c) = min k∈M - j (H (c k )), if M + j = ∅, then v j, (c) = 1, else v j, (c) = min k∈M + j (H (c k )). ( 18
)
It can be noted that these regularized reaction rates satisfy the requirements [START_REF] Erhel | Solving partial differential algebraic equations and reactive transport models[END_REF]. Now, we define a regularized system of kinetic reactions:

dc(t) dt = Sr (c(t)), t ≥ 0, c(0) = c 0 . (19)
It is easy to show that the regularized system (19) has a unique solution. Proof. By using Lemma 5.1 and linearity, f is locally Lipschitz continuous. Let c a vector of species quantities such that c ≥ 0.

First, let A k , 1 ≤ k ≤ N a be an aqueous species such that c k = 0 for some k ∈ {1, . . . N a }. We proved in Proposition 2.2 that s kj τ j ≥ 0, so using 13, s kj r j, ≥ 0 and (Sr ) 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 (Sr j, ) k ≥ 0. Proposition 5.2. If c 0 ≥ 0 then the solution of system [START_REF] Kräutle | Robust simulation of mineral precipitation-dissolution problems with variable mineral surface area[END_REF] satisfies c ≥ 0 in the maximal interval.

Proof. As in Proposition 2.3, it is a direct consequence of Lemmas 5.1 and 5.2.

Convergence of the regularized solutions

In this section, we prove convergence of the regularized solutions, towards a non-negative and bounded weak solution of the differential inclusion problem. Recall that the vector q of molar mass satisfies q > 0 and S T q = 0. Using Lemma 2.1, we prove that the regularized solution is bounded and can be defined in [0, +∞[: Proposition 5.3. Let q be the vector of molar mass. System [START_REF] Kräutle | Robust simulation of mineral precipitation-dissolution problems with variable mineral surface area[END_REF]

has a unique solution c in [0, +∞[ which satisfies 0 ≤ c ≤ C,
where C is defined in Lemma 2.1.

Proof. Let c be the maximal solution of [START_REF] Kräutle | Robust simulation of mineral precipitation-dissolution problems with variable mineral surface area[END_REF]. Since S T q = 0, it satisfies q T c = q T c 0 . Thanks to Proposition 5.2, it also satisfies c ≥ 0. Thus, using Lemma 2.1, it also satisfies c ≤ C. So, the maximal interval of the non-negative solution is [0, +∞[.

It also implies that the sequence c converges.

Theorem 5.1. Let T > 0. The sequence c converges strongly in C 0 (0, T, R (Na+Np) ).

Any limit c of a subsequence satisfies 0 ≤ c ≤ C. Also for any j = 1, . . . , N r , subsequences of the regularized controls u j, (c ) and v j, (c ) converge weakly-* in L ∞ (0, T, R Nr ) 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 (t) ∈ U j (c p (t)) and v j (t) ∈ V j (c p (t)).

Proof. Let c be the solution of [START_REF] Kräutle | Robust simulation of mineral precipitation-dissolution problems with variable mineral surface area[END_REF] in [0, +∞[. Since 0 ≤ c ≤ C, the norms c L 2 (0,T ) and Sr (c ) L 2 (0,T ) are bounded, thus c H 1 (0,T ) is also bounded.

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

Let 1 ≤ j ≤ N r . Since 0 ≤ u j, (c ) ≤ 1, a subsequence of u j, (c ) converges weakly-* in L ∞ (0, T ) and the limit u j is bounded below by 0 and above by 1. The result is similar for v j . Thus ∀j = 1, . . . , N r , the control satisfies u j ∈ L ∞ (0, T ) and 0 ≤ u j (t) ≤ 1, 0 < t < T. Now, let 0 < t < T . We want to show that u j (t) ∈ U j (c p (t)). In other words, we want to show that if M - j = ∅ and if ∀k ∈ M - j , c k (t) > 0, then u j (t) = 1. Since c is continuous, there exists an interval I containing t such that ∀s ∈ I, ∀k ∈ M - j , c k (s) > 0. Since a subsequence of c converges towards c in C 0 (0, T ), the same property is true for sufficiently small: ∀s ∈ I, ∀k ∈ M - j , c k, (s) > 0. Consequently, ∀s ∈ I, u j, (c )(s) = 1. Since u j, (c ) is a constant function of time in I, it converges towards 1 and u j (t) = 1. This concludes the proof. Now, we can prove the existence of a weak solution of system (15). Theorem 5.2. Let T > 0 and c a limit of c in C 0 (0, T, R (Na+Np) ). Let u and v weakly-* limits of a subsequence of u (c ) and v (c ) in L ∞ (0, T, R Nr ). Then (c, u, v) is a weak solution of system [START_REF] Friedly | Solute transport with multiple equilibrium-controlled or kinetically-controlled chemical reactions[END_REF]. Proof. First, the limits u j and v j satisfy the requirements of definition 5.2. Also, c(0) = c 0 and c ≥ 0. Since c weakly converges towards c in H

1 (0, T ), ∀φ ∈ (C ∞ c (0, T )), lim →0 T 0 dc (s) dt - dc(s) dt φ(s) ds = 0. ( 20
)
Thanks to Lemma 5.1, u j, , u j ∈ L ∞ (0, T ) and u j, ∞ ≤ 1. Also, since τ j and c , c are continuous, τ + j (c), τ + j (c ) ∈ L 1 (0, T ). Let us decompose r j, -r j in the following way:

r j, -r j = u j, (c )(τ + j (c ) -τ + j (c)), + (u j, (c ) -u j )τ + j (c), -v j, (c )(τ - j (c ) -τ - j (c)), -(v j, (c ) -v j )τ - j (c).
Then ∀φ ∈ C ∞ c (0, T ),

| T 0 (r j, (c (s)) -r j (s, c(s)))φ(s) ds| ≤ τ + j (c ) -τ + j (c) L 1 (0,T ) φ L ∞ (0,T ) + | T 0 (u j, (c (s)) -u j (s))τ + j (c(s))φ(s) ds| + τ - j (c ) -τ - j (c)) L 1 (0,T ) φ L ∞ (0,T ) + | T 0 (v j, (c (s)) -v j (s))τ - j (c(s))φ(s) ds|.
Since c converges towards c in C 0 (0, T ) and τ j is continuous, we get lim →0 τ + j (c )τ + j (c) L 1 (0,T ) = 0. Since a subsequence of u j, (c ) weakly-* converges to u j in L ∞ (0, T ) and τ + j (c)φ ∈ L 1 (0, T ), we also obtain lim

Using both weak convergence results ( 20) and ( 21), we conclude that c, u, v are a weak solution of system [START_REF] Friedly | Solute transport with multiple equilibrium-controlled or kinetically-controlled chemical reactions[END_REF].

Test case with a mineral involved in two reactions

The system we consider here is a set of two reactions with a mineral which participates in both reactions. This system is similar to the two reactions between calcite, calcium oxide and calcium hydroxide described in Section 2.1. Other species are also minerals in order to get constant reaction rates, so that reactions are written M 1 M 2 and M 2 M 3 . The objective is to analyze reaction rates after dissolution of one or several minerals. We can expect that the first reaction stops if M 1 dissolves completely, as for the reference example. Similarly, the second reaction stops if M 3 is fully dissolved. But the case where M 2 precipitates in the first reaction and dissolves in the second one is more complex. We analyze below the different situations and compute the limit of the regularized solution.

The two maximal reaction rates defined by the mass action laws (3) are τ = (τ 1 , τ 2 ). System [START_REF] Friedly | Solute transport with multiple equilibrium-controlled or kinetically-controlled chemical reactions[END_REF] becomes

                   dc 1 dt = -r 1 , dc 2 dt = r 1 -r 2 , dc 3 dt = r 2 , c(0) = c 0 , c ≥ 0, (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 = (1, 1, 1) T is a basis vector satisfying q > 0. The conservative variable (up to a constant) is c 1 +c 2 +c 3 so that it is possible to consider only two equations out of the three in [START_REF] Shampine | The matlab ode suite[END_REF].

The unique solution c * of the continuous model ( 4) is a linear function of time:

   c * 1 (t) = c 1,0 -τ 1 t, c * 2 (t) = c 2,0 + (τ 1 -τ 2 )t, c * 3 (t) = 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 . ( 23
)
Let c * ,i be the solution of the regularized system [START_REF] Kräutle | Robust simulation of mineral precipitation-dissolution problems with variable mineral surface area[END_REF], where the function γ is identity. Since there is only one dissolving mineral in each reaction, there is no blending in the regularized model. In other words, the regularized functions are H (c k ), where c k is the dissolving mineral. We analyze convergence in different cases and show stronger results than Theorem 5.2. Indeed, for almost any t > 0, the regularized controls u (c (t)) and v (c (t)) converge towards discontinuous functions u(c(t)) and v(c(t)).

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

r 1, = H (c 1 )τ 1 , r 2, = H (c 3 )τ 2 . ( 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

c 1, (t) = c 1,0 -τ 1 t = c * 1 (t), t ≤ t 1, , c 1, (t) = exp(- τ 1 (t -t 1, )), t ≥ t 1, ,
and the regularized control is

H (c 1, (t)) = 1, t ≤ t 1, , H (c 1, (t)) = exp(- τ 1 (t -t 1, )), t ≥ t 1, .
It is easy to verify that ∀t ≥ 0, lim

→0 c 1, (t) = c 1 (t)
where c 1 is defined by

c 1 (t) = c 1,0 -τ 1 t = c * 1 (t), t ≤ t * 1 , c 1 (t) = 0, t ≥ t * 1 .
Also, the control converges: ∀t ≥ 0,

t = t * 1 lim →0 H (c 1, (t)) = H(c 1 (t))
. Clearly, results are the same for the variable c 3 and c 2 is deduced from the conservative quantity. Thus c is solution of ( 22) where

r 1 = H(c 1 )τ 1 , r 2 = H(c 3 )τ 2 .

Dissolution of M 2

Here, we assume that τ 1 < 0 and τ 2 > 0 and c 2,0 > 0. The pure phase M 2 dissolves and produces M 1 and M 3 . The regularized reactions rates are given by

r 1, = H (c 2 )τ 1 , r 2, = H (c 2 )τ 2 . ( 25
)
In that case, the second equation is independent and similar to the previous case. In the same way as t 1, , we define, with 0 < < c 2,0 , 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 (τ 2 -τ 1 ) instead of τ 1 ), and the control H (c 2, (t)) converges towards H(c 2 (t)) for t = t * 2 . It is slightly more complicated to compute the solution c 1, :

c 1, (t) = c 1,0 -τ 1 t = c * 1 (t), t ≤ t 2, , c 1, (t) = - τ 1 τ 2 -τ 1 (1 -exp(- τ 2 -τ 1 (t -t 2, ))) + c * 1 (t 2, ), t ≥ t 2, .
Therefore, for all t ≥ 0, the regularized solution converges towards c 1 defined by

c 1 (t) = c 1,0 -τ 1 t = c * 1 (t), t ≤ t * 2 , c 1 (t) = c * 1 (t 2 ), t ≥ t * 2 .
The third quantity c 3 is obtained thanks to the conservation law. The limit c is solution of [START_REF] Shampine | The matlab ode suite[END_REF] where reaction rates are the discontinuous functions

r 1 = H(c 2 )τ 1 , r 2 = H(c 2 )τ 2 .
As illustrated in Figure 3, minerals M 1 and M 3 precipitate until M 2 is fully dissolved and both reactions R 1 and R 2 stop.

Dissolution of M 1 then M 2

Here, we assume that τ 2 > 0 and τ 1 > 0, also c 1,0 > 0 and c 2,0 > 0; moreover, if τ 2 > τ 1 , we assume that t * 1 < t * 2 . In these conditions, 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.

With these assumptions and with sufficiently small,c * 2 (t) ≥ for t ≤ t 1, . As long as c 2, (t) ≥ , it is given by

c 2, (t) = c 2,0 -(τ 2 -τ 1 )t = c * 2 (t), t ≤ t 1, , c 2, (t) = (1 -exp(- τ 1 (t -t 1, ))) -τ 2 (t -t 1, ) + c * 2 (t 1, ), t ≥ t 1, .
Since c 2, (t 1, ) ≥ and lim t→+∞ c 2, (t) = -∞, by continuity, there exists t2, ≥ t 1, such that c 2, ( t2, ) = . Then, for t ≥ t2, the equation becomes

dc 2 dt = c 1 τ 1 - c 2 τ 2 , t ≥ t2, with c 2 ( t2, ) = . The solution is c 2, (t) = (a exp(- τ 1 (t -t 1, )) + (1 -ab ) exp(- τ 2 (t -t2, ))) where a = τ 1 τ 2 -τ 1 and b = exp(- τ 1 ( t2, -t 1, )).
Finally, the regularized quantity c 2, satisfies the equation

       c 2, (t) = c 2,0 -(τ 2 -τ 1 )t = c * 2 (t), t ≤ t 1, , c 2, (t) = (1 -exp(- τ 1 (t -t 1, ))) -τ 2 (t -t 1, ) + c * 2 (t 1, ), t 1, ≤ t ≤ t2, , c 2, (t) = (a exp(- τ 1 (t -t 1, )) + (1 -ab ) exp(- τ 2 (t -t2, ))) t ≥ t2, .
Let us now analyze the convergence of t2, and c 2, (t) when tends towards 0. The equality c 2, ( t2, ) = can be rewritten

-exp(- τ 1 ( t2, -t 1, ))) -τ 2 ( t2, -t 1, ) + c * 2 (t 1, ) = 0.
Since t2, ≥ t 1, , it implies that t2, has a limit which we denote by t2 , and which satisfies

τ 2 ( t2 -t * 1 ) = c * 2 (t * 1 ) = 0. Then ∀t ≥ 0, lim →0 c 2, (t) = c 2 (t) given by    c 2 (t) = c * 2 (t), t ≥ t * 1 , c 2 (t) = -τ 2 (t -t * 1 ) + c * 2 (t * 1 ), t * 1 ≤ t ≤ t2 , c 2 (t) = 0, t ≥ t2 . Moreover, ∀t ≥ 0, t = t * 1 , t = t2 , lim →0 H (c 2, (t)) = H(c 2 (t)).
Finally, c is solution of the discontinuous system [START_REF] Shampine | The matlab ode suite[END_REF] with

r 1 = H(c 1 )τ 1 , r 2 = H(c 2 )τ 2 .
In Figure 4 where τ 1 > τ 2 , mineral M 1 dissolves at rate τ 1 and mineral M 2 precipitates at rate τ 1 -τ 2 , then reaction R 1 stops and M 2 dissolves at rate τ 2 until reaction R 2 stops. Up to now, the regularized Heaviside functions converged towards the Heaviside function and reactions stopped when a mineral was fully dissolved. This case is more remarkable as it illustrates the capacity of the regularized model to continue the reactions 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 (c 1 )τ 1 -H (c 2 )τ 2 ,
so that the regularized solution c 2, satisfies

               c 2, (t) = c 2,0 -(τ 2 -τ 1 )t = c * 2 (t), t ≤ t 2, , c 2, (t) = τ 2 (τ 1 + (τ 2 -τ 1 ) exp(- τ 2 (t -t 2, ))), t 2, ≤ t ≤ t 1, , c 2, (t) = ( τ 2 -τ 1 τ 2 exp(- τ 2 (t -t 2, )) - τ 2 1 τ 2 (τ 2 -τ 1 ) exp(- τ 2 (t -t 1, ))+ τ 1 τ 2 -τ 1 exp(- τ 1 (t -t 1, ))), t ≥ t 1, .
Thus the regularized solution c 2, converges, for all t ≥ 0, towards c 2 defined by

c 2 (t) = c 2,0 -(τ 2 -τ 1 )t = c * 2 (t), t ≤ t * 2 , c 2 (t) = 0, t ≥ t * 2 .
Also, the regularized control H (c 2, ) has a special limit, for any t ≥ 0,

t = t * 1 , t = t * 2 : lim →0 H (c 2, (t)) = H(c 2 (t)) + H(c 1 (t))δ(c 2 (t)) τ 1 τ 2 ,
where δ(x) is the Dirac function at point x. The limit c of the regularized solution satisfies the discontinuous differential equations [START_REF] Shampine | The matlab ode suite[END_REF] with

r 1 = H(c 1 )τ 1 , r 2 = H(c 2 )τ 2 + H(c 1 )δ(c 2 )τ 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 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 6.1. The limit c of the regularized solution c is solution of system [START_REF] Shampine | The matlab ode suite[END_REF] where the reaction rates are discontinuous functions defined by:

r 1 = H(c 1 )τ + 1 -H(c 2 )τ - 1 -δ(c 2 )H(c 3 )τ - 2 , r 2 = H(c 2 )τ + 2 + δ(c 2 )H(c 1 )τ + 1 -H(c 3 )τ - 2 . (26)
Moreover, the regularized limiters u (c ) and v (c ) strongly converge, up to a subsequence, in

L p (0, T ), p ∈ [1, ∞[.
Proof. Existence is proved above by construction. We proved in all cases that u (c ) converges almost everywhere to u in [0, T ]. Thanks to the dominated convergence theorem applied to |u (c (t)) -u(t)| 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 (t * ) = 0 then c i (t) = 0, t ≥ t * so that dc i dt = 0, t > t * . Finally, the reaction rates satisfy (Sr) i = 0, t > t * . For example, when c 2 becomes fully dissolved, then r 2 = r 1 . This result is a so-called sliding mode of differential inclusion problems [START_REF] Filippov | Differential Equations with Discontinuous Righthand Sides[END_REF].

Numerical experiments

In this section we describe and discuss the results of five numerical experiments done with synthetic geochemical systems. The first one is the reference example A M defined in Section 4. The second one, with two reactions M 1 M 2 and M 2 M 3 , is analyzed in Section 6. Then we replace M 1 and M 3 by aqueous species A 1 and A 3 in the third test case. In order to introduce competition between two reactions dissolving M 2 , we add a third reaction M 2 A 4 in the fourth test case. These four examples do not require blending in the regularized model, in contrary to the last one, which involves two reactant minerals in the same reaction.

Numerical context

We have developed a Matlab code for solving the regularized ODE [START_REF] Kräutle | Robust simulation of mineral precipitation-dissolution problems with variable mineral surface area[END_REF], where we vary the value of the regularizing parameter in order to analyze numerical convergence. We have implemented several regularizing functions defined 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], namely a linear function, an exponential smoothing, and a polynomial one. We show results with the linear function, other results were similar. For the last test case, we implemented two types of blending: the min function defined in [START_REF] Kräutle | Existence of global solutions of multicomponent reactive transport problems with mass action kinetics in porous media[END_REF], and the prod function defined below:

     if M - j = ∅, then u j, (c) = 1, else u j, (c) = k∈M - j (H (c k )), if M + j = ∅, then v j, (c) = 1, else v j, (c) = k∈M + j (H (c k )). (27)
The maximal reaction rates are implemented in a Matlab function, using data of the physical model, mainly stoichiometry and constants in [START_REF] Angeli | A tutorial on chemical reaction network dynamics[END_REF]. 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 (c a ) = c i c H 2 O + Na l=1 c a,l , (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 constant. 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 [0, T ] and the initial condition c 0 0 are parameters.

We define a set of discrete times in [0, T ]:

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 (t k ), where c (t k ) is a vector of size N a + N p .

In the first two tests, we know the limit c of the regularized model, exact solution of the differential inclusion model. Therefore we can compute a discrete error for each species i:

e c i (t k ) = |c i (t k ) -c ,i (t k )|, k = 0, • • • , N f ,
and a global error:

e c (t k ) = Na+Np i=1 e c i (t k ). ( 29
)
We also analyze the errors in the limiting discontinuous functions u ,j and v ,j for each reaction R j . In the first two tests, the number N l of discontinuity points is finite and the corresponding times t * l , l = 1, • • • , N l , are known. As above, we compute discrete errors for each reaction R j : (30)

e u j (t k ) = |u j (t k ) -u ,j (t k )|, k = 0, • • • , N f , e v j (t k ) = |v j (t k ) -v ,j (t k )|, k = 0, • • • , N f ,
Then we compute discrete L p -norms with p = 1, 2, ∞. We also define a discrete error at times t k outside [t * l -δ, t * l + δ], ∀l = 1, • • • , N l , where δ > 0. Then we compute the discrete L ∞ norm out of the jumps defined by these intervals.

Reference example A -M , total dissolution

This test case corresponds to the reference example of Section 4:

A M,
with A a solute in the water phase and M a mineral. We consider the total dissolution case, as in Section 4.3, with data given in table 1 . The activity of A is defined by (28) and the reaction rate by [START_REF] Angeli | A tutorial on chemical reaction network dynamics[END_REF], so that

a(c A ) := c A c A + c H 2 O and τ (c A ) := αa(c A ) -β.
Since there is only one mineral product in the chemical reaction, the limiter is v = H (c M, ).

Reaction A M S -1 1 (α, β) (2, 2) c H 2 O 2 c 0 (0.2 0.6)
Table 1: Test case A -M : physical parameters Figure 6 shows the regularized solution and the reaction rate computed with = 0.001. We observe that the limit of the regularized quantities is the discontinuous solution of Section 3. Physical data are such that the mineral is fully dissolved at time t * and for t ≥ t * , the quantity of mineral c M remains equal to 0. The reaction stops: r = τ until time t * then r = 0. Now we run the regularized model for a decreasing sequence of parameters and compute errors. We choose T = 0.5, dt = T 10 4 = 5 × 10 -5 and δ = 0.1. We observe in Figure 7 that are localized around the discontinuity. As expected, when tends towards 0, the error e c converges towards 0 at any time, and the error e u converges towards 0 except at the time t * of discontinuity. In Figure 8, we plot L 1 and L ∞ norms of errors e c and e u in function of the parameter . As proved in Lemma 5.1, the error e c converges in both norms whereas the error e u converges only in L 1 norm. Nevertheless, for this test case, the error e u converges numerically in the L ∞ norm, out of jumps.

Test case

M 1 -M 2 -M 3
Now, we investigate the case where one mineral participates in two reactions, as a reactant in one of them and a product in the other. The geochemical system is described in Section 6:

M 1 M 2 and M 2 M 3 ,
where M 1 , M 2 and M 3 are three pure phase minerals, so that activities are equal to 1 and maximal reaction rates τ 1 and τ 2 are constant. We consider the case of Section 6.4, where the mineral M 2 precipitates in the first reaction and is dissolved in the second one.

Table 2 provides a summary of physical data.

Reactions

R 1 : M 1 M 2 R 2 : M 2 M 3 S   -1 0 1 -1 0 1   (τ 1 , τ 2 ) (1, 4) c 0 (1, 0.6, 0.4) Table 2: Test case M 1 -M 2 -M 3 : physical parameters
When we apply solver ode15s directly to the discontinuous system (26), with various values of , the solver gives either a false result or no result at all. This failure highlights the need for a numerical regularization. In order to justify this observation, we analyze numerical convergence for a decreasing sequence of . We choose T = 2, dt = 2 10 4 = 2 × 10 -4 and δ = 0.1, as defined in Section 7. In Figure 10, we plot the regularized quantities (c 1 , c 2 , c 3 ) and controls (u 1 , u 2 ), as well as the errors e c and e u , in function of . We observe, as proved in Theorem 5.2, that the quantities c converge strongly towards c in C 0 , with an error e c maximal around the times of discontinuity. Also, the controls H (c 1 ) and H (c 2 ) converge at any time except those of discontinuity towards respectively H(c 1 ) and H(c 2 ) + δ(c 2 )H(c 1 ) τ 1 τ 2 , as expected from Section 6.4. Figure 11 represents the norms of the errors in the quantities e c and in the controls e u . Numerically, the rate of convergence of e c is almost 1 for both norms L 1 and L ∞ . In Theorem 5.1, we proved strong convergence of e u in L 1 norm. Numerically, the L 1 -norm of e u converges with a rate almost equal to 1. Also, we observe a fast convergence of e u in the L ∞ norm, out of jumps.

Test case

A 1 -M 2 -A 3
This test case is very similar to the previous one, except that species M 1 and M 3 are replaced by aquesous species A 1 and A 3 :

A 1 M 2 and M 2 A 3 ,
where A 1 and A 3 are aqueous species and M 2 a pure phase mineral. Physical data, given in Table 3, are such that mineral M 2 precipitates in the first reaction and dissolves in the second one. Activities of water solute A 1 and A 3 are given by ( 28) and maximal reaction rates by [START_REF] Angeli | A tutorial on chemical reaction network dynamics[END_REF]. They are non constant functions of c a :

a i (c a ) = c A i c A 1 + c A 3 + c H 2 O , i = 1 or 3, τ 1 (c a ) = α 1 a 1 (c a ) -β 1 and τ 2 (c a ) = α 2 -β 2 a 3 (c a )
Here the first reaction is not limited so that r 1 = τ 1 , whereas r 2, = H (c 2, ).

Reactions Figure 12 shows the regularized solution and reaction rates computed with = 0.001. We get a good approximation of the limit, solution of the discontinuous equations [START_REF] Friedly | Solute transport with multiple equilibrium-controlled or kinetically-controlled chemical reactions[END_REF]. We observe, like in the previous examples, that once the mineral M 2 is fully dissolved, it remains so. Thus the derivative dc 2 /dt is equal to 0, in other words r 2 = τ 1 , meaning that H (c 2, ) converges towards H(c 2 ) + δ(c 2 ) τ 1 τ 2 , as in the previous example.

R 1 : A 1 M 2 R 2 : M 2 A 3 S   -1 0 1 -1 0 1   (α 1 , β 1 ) (α 2 , β 2 ) (30, 3) (4.5, 18.75) c 0 (2, 5, 0) c H 2 O 20 Table 3: Test case A 1 -M 2 -A 3 : physical parameters

Test case with competition

A 1 -M 2 -A 3 and M 2 -A 4
In this test case, we still consider the two reactions of the previous example and add a reaction of dissolution:

A 1 M 2 and M 2 A 3 and M 2 A 4 ,
where M 2 is a pure phase mineral, A 1 , A 3 and A 4 are species in aqueous phase. Physical data are such that M 2 gets fully dissolved at a time of discontinuity t * , see Table 4. The first reaction is not limited thus r 1 = τ 1 whereas r 2 = u 2 (c 2 )τ 2 and r 3 = u 3 (c 2 )τ 3 , with u 2 = u 3 = 1 if c 2 > 0. The activity of water solute c a is given by ( 28) and maximal reaction rates by (3): The regularized limiters are u 2, = u 3, = H (c 2, ) thus, at convergence, u 2 = u 3 . Figure 13 shows the regularized solution and reaction rates computed with = 0.001. We get an approximation of the discontinuous solution. Again, we observe that once M 2 is dissolved at time t * , the quantity c 2 remains equal to 0. Thus dc 2 /dt = 0 for t > t * , in other words r 2 + r 3 = r 1 , so that u 2 = u 3 = τ 1 τ 2 +τ 3 . In some sense, the rate r 1 is evenly distributed between r 2 and r 3 , as if A 1 produces A 3 and A 4 .

a i (c a ) = c A i c A 1 + c A 3 + c A 4 + c H 2 O , i = 1, 3 or 4, τ 1 (c a ) = α 1 a 1 (c a ) -β 1 , τ 2 (c a ) = α 2 -β 2 a 3 (c a ) and τ 3 (c a ) = α 3 -β 3 a 4 (c a ). 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     (α 1 , β 1 ) (α 2 , β 2 ) (α 3 , β 3 ) ( 30

Test case with blending

Up to now, there was only one dissolving mineral in a kinetic reaction. We investigate cases where at least two minerals dissolve in the same reaction, so that our regularized model requires using a blending function. We consider the following example: where M i , i = 1, • • • , 7, are pure phase minerals, so that maximal reaction rates are constant. Table 5 contains a summary of physical parameters.

Reactions Table 5: Test case with blending: physical parameters

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
We use two blendings in the regularized model, the min function as in [START_REF] Kräutle | Existence of global solutions of multicomponent reactive transport problems with mass action kinetics in porous media[END_REF] and the prod function as in (27). As long as M 1 and M 4 are not fully dissolved, u 1, = u 2, = 1. Since u 3, = min(u 4, (c 1, ), u 5, (c 2, )), at convergence u 3 = min(u 4 , u 5 ) (similarly with prod).

In Figure 14, we plot the controls u j, obtained with the two blendings, when eps=0.001. As expected, we observe convergence except at times of discontinuity. Under these conditions, M 1 and M 4 are not yet fully dissolved so that u 1 = u 2 = 1, whereas M 2 is first fully dissolved, followed by M 1 . Again, after dissolution of a mineral, its time derivative is equal to 0. This observation gives rise to equations satisfied by the controls u 4 and u 5 . Between both dissolutions, u 4 = 1 and r 2 = r 3 + r 5 , yielding u Consequently, the discontinuous reaction rates are different and the quantities c are also different with these two blendings. For example, dc 6 dt = u 4 , dc 7 dt = u 5 . So, we observe that the regularized limiters and quantities depend on the blending, as well as their limits. Thus the choice of the blending has a direct impact on the differential inclusion problem. In a future work, we plan to compare in more details the different blendings.

Discussion and conclusion

In this paper, we consider general precipitation dissolution reactions, where minerals can interact with other minerals and participate in several reactions. We propose a differential inclusion model with reaction rates defined by set-valued functions. Indeed, when a mineral is dissolved, reactions can stop or continue more slowly, thanks to limiting functions taking value in the interval [0, 1]. We prove existence of a continuous and positive solution to this problem by regularizing the reaction rates. This regularized model allows us also to compute numerical approximations of a discontinuous solution.

In our tests, we observe that once a mineral is dissolved, it remains so, at least during a time interval. This provides equations to get unique values of the limiting functions and the discontinuous reaction rates. We plan to further analyze this sliding mode, in order to characterize the set-valued functions in the differential inclusion problem. We also plan to study whether the limit of the regularized solution is the unique solution of the differential inclusion problem. When blending is required, regularized reaction rates as well as discontinuous ones depend on the choice of blending. We plan to compare these different models from a chemical point of view. In our differential inclusion model, the geochemical system is closed, with no external source term. Also, geochemical interactions are often coupled with transport processes such as dispersion or conduction. Future work will aim at adding source terms and at studying reactive transport models with these kinetic reactions.

  Let us introduce the positive part τ + (c a ) and the negative part τ -(c a ) of the continuous reaction rate so that τ (c a ) = τ + (c a ) -τ -(c a ).

Figure 1 :

 1 Figure 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 5 . 1 .Proposition 5 . 1 .

 5151 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. For any T > 0, system[START_REF] Kräutle | Robust simulation of mineral precipitation-dissolution problems with variable mineral surface area[END_REF] with non-negative initial condition c 0 , has a unique maximal solution c in [0, T [. Proof. We apply Cauchy-Lipschitz theorem. Now, in order to prove the non negativity of the solution, we start by showing that the right-hand side is essentially non negative.Lemma 5.2. The function f defined by f (c) = Sr (c) is locally Lipschitz continuous and essentially non negative.

→0 T 0 0 (

 00 (u j, (c ) -u j )τ + j (c)φ ds = 0. Similar results apply to v j and τ - j . So, ∀φ ∈ C ∞ c (0, T ), lim →0 T Sr (c (s)) -Sr(s, c(s)))φ(s) ds = 0.

Figure 2

 2 Figure 2 illustrates this case. Mineral M 2 precipitates at rate (τ 1 -τ 2 ) 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 2 :

 2 Figure 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 3 :

 3 Figure 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 4 :

 4 Figure 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 5 :

 5 Figure 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.

e

  u j (t k ), e v (t k ) = Nr j=1 e v j (t k ).

Figure 6 :

 6 Figure 6: Case A -M with dissolution of M . Results with = 0.001.

Figure 7 :

 7 Figure 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 c (left) and error e u (right). The interval [t * -δ, t * + δ] is represented in pink.

Figure 8 :

 8 Figure 8: Numerical convergence for the case A -M with dissolution of M .

Figure 9 :

 9 Figure 9: Case M 1 -M 2 -M 3 with dissolution of M 2 then M 1 . Results with = 0.001.

Figure 9

 9 Figure9shows the regularized solution and both limiters obtained with = 0.001. Qualitatively, the quantities computed are close to the discontinuous ones given in Section 6.4.

Figure 10 :

 10 Figure 10: Case M 1 -M 2 -M 3 with dissolution of M 2 then M 1 . Results and errors for different values of . Top: quantities (left) and controls u (right). Bottom: errors e c (left) and e u (right). Both intervals [t * -δ, t * + δ] are represented in pink.

Figure 11 :

 11 Figure 11: Numerical convergence for the case M 1 -M 2 -M 3 in dissolution of M 2 then M 1 .

Figure 12 :

 12 Figure 12: Case A 1 -M 2 -A 3 with total dissolution of M 2 . Results with = 0.001.

20 Table 4 :Figure 13 :

 20413 Figure 13: Case Competition with aqueous species. Results with = 0.001.

M 3 M 1 and M 4 M 2

 342 and M 1 + M 2 M 5 and M 1 M 6 and M 2 M 7 ,

5 = τ 2 τ 3 1 τ 3 + τ 4 and u 5 = τ 2 -τ 3 u 4 τ 5 .

 5313525 + τ 5 = 0.35. After dissolution of the two minerals, r 1 = r 3 + r 4 and r 2 = r 3 + r 5 . With the min blending, we get min(u 4 , u 5 ) = u 4 , so that u 4 = τ Numerical values are u 4 = 0.25, u 5 = 0.45, u 3 = min(u 4 , u 5 ) = 0.25. , u 3 = u 4 u 5 0.1726. (32)

Figure 14 :

 14 Figure 14: Case with two different blendings. Results with = 0.001.