Luis Loaeza

Ricardo Corona-Sánchez

Geraldyne Castro

Margarita Romero-Ávila

Rosa Santillan

Valérie Maraval

Remi Chauvin

Norberto Farfán
email: n.f.norberto.farfan@gmail.comr.c.chauvin@lcc-toulouse.fr

Synthesis and optical properties of 1-ethyl-indol-3-yl-substituted aza-BODIPY dyes at the 1,7positions

Keywords:

New aza-BODIPYs with significant bathochromic shifts in their UV-vis absorption spectrum were synthesized by the installation of 1-ethyl-indol-3-yl moieties at the 1,7-positions of the aza-BODIPY core bearing 4-substituted phenyl or thien-2-yl groups at the 3,5-positions. The compounds described herein show strong NIR absorptions up to 779 nm, and emissions up to 805 nm in chloroform solution, and present large molar extinction coefficients. The absorption and emission bands of aza-BODIPYs were significantly red-shifted, by ca 100 nm, relative to the 3,5-dianisyl-1,7diphenyl aza-BODIPY. TD-DFT calculation indicate an effective electronic conjugation between the aza-BODIPY core and the substituents. Absorption of the aza-BODIPYs corresponds to a charge transfer from donor substituents toward the boron-centered core. Moreover, these computational results show that all the compounds present HOMO-LUMO band gap values that are in the organic semiconductor range.

Introduction

In recent years, the interest in molecules absorbing and emitting in the near-infrared region (NIR) has increased, especially for the purpose of medicine and materials science. In the domain of medicine, the region between 700 and 1100 nm, denominated the "biological window", is attractive for numerous biological applications because in this range of wavelength the light has deep tissue penetration and does not damage living tissues. 1 In materials science, wavelengths extending in the NIR region are particularly important for optoelectronic applications. Indeed, about 51% of solar energy falls in the near-infrared region, 2 calling for the development of new dyes and materials absorbing in this region to increase the efficiency in solar energy conversion (e.g. by organic photovoltaic solar cells or water-splitting photocatalysts 3 ). Amongst the large variety of NIR chromophores, aza-BODIPY dyes have attracted significant attention due to their remarkable photophysical properties such as tunable absorption and emission wavelengths, high photostability, and large molar absorption coefficients. 4 These characteristics make them promising candidates for potential applications, e.g. in bioimaging, 5 photodynamic therapy, 6 sensing, 7 optoelectronics, 8 and general nonlinear optics. 9

Taking as reference the 1,3,5,7-tetraphenyl-aza-BODIPY which absorbs at 650 nm and emits at 672 nm in chloroform, 10 several structural modifications have been reported to shift the absorption and emission bands further to the NIR range. These approaches include the substitution by electron-donating groups 9,11 -or both electron-donating and electron-withdrawing groups- 12 on the phenyl rings, rigidification of rotatable moieties, 13 extension of the π-conjugation, 14 benzannulation of the pyrrole ring, 15 formation of intramolecular B-O six-membered rings, 16 and reduction of the torsion angles by replacing phenyl substituents by thienyls. 17 Nevertheless, one of the most popular general strategy to modify the optical properties of π-conjugated molecules consists in the introduction of resonating donor (D) and acceptor (A) motifs to form push-pull systems. 18 This approach applied to aza-BODIPYs has proven to contribute significantly to a bathochromic shift of their absorption to the NIR region, as, for example, in 1,7-difluorophoresubstituted aza-BODIPYs in which a naphthalene, pyrene, or N-substituted carbazole group operates as the donor unit and the aza-BODIPY core as the acceptor fragment. 19 A large variety of aza-BODIPYs with diverse substituents at the 1,7-positions was reported, but to our knowledge, strongly electron-donating 3-indolyl moieties 20 have not been envisaged at those positions yet, though it has been shown that such motifs at the 2,6-positions improve the optical properties of aza-BODIPYs. 21 Therefore, in this work, 1-ethyl-indol-3-yl substituents were introduced in the 1,7-positions with different substituents at the 3,5-positions of the aza-BODIPY core (Figure 1), with the view to studying their effect on the photophysical properties of the obtained chromophores. These 1-ethyl-indol-3-yl-substituted aza-BODIPYs are shown to exhibit red-shifted absorption and emission maxima. Additionally, TD-DFT computational studies were carried out to present a complementary spectroscopic insight with the purpose of estimating the wavelength of maximum absorption of this type of dyes.

Results and Discussion

Synthesis and Characterization

The synthesis of the target compounds 4a-c (Scheme 1) was performed using the O'Shea approach. 9,22 This strategy involves the formation of the α,β-unsaturated ketones 1a-c by an aldol condensation reaction using one equivalent of 1-ethylindole-3-carboxaldehyde and one equivalent of the selected methylketone in the presence of KOH. The enones 1a-c were found to precipitate as yellow solids in 84-89% yields. Then, the 1,3-diaryl-4-nitrobutan-1-ones 2a-c were synthesized in 80-85% yields via a 1,4 Michael addition reaction of nitromethane to the corresponding α,βunsaturated ketones 1a-c. The obtained adducts 2a-c were subsequently treated with an excess of ammonium acetate in butanol to give the aza-dipyrromethenes 3a-c in 25-34% yields. Finally, the the aza-BODIPYs 4a-c were obtained by reaction of the corresponding aza-dipyrromethene with BF 3 •Et 2 O in the presence of triethylamine. The targeted aza-BODIPYs 4a-c were isolated as dark blue-violet solids in 79-88% yields. Particularly, the 11 B chemical shift value of the thienyl substituted aza-BODIPY 4c was found slightly high frequency shifted, compared to 4a and 4b as a result of an increased delocalization of electronic density in the aza-BODIPY core of 4c. The 3,5-dianisyl-1,7-diphenyl aza-BODIPY 4aPh was also prepared according to a literature protocol 9 and was used to compare their optical properties versus the compounds 4a-c. The optical properties of the aza-BODIPYs 4a-c were evaluated in three solvents with different polarities: tetrahydrofuran, chloroform, and toluene. The three dyes 4a-c present absorption and emission in the NIR region with large molar extinction coefficients. The values of the wavelength of maximum absorption, emission, molar extinction coefficients, and relevant photophysical data are listed in Table 2, and compared to those of the selected 3,5-dianisyl-1,7-diphenyl-aza-BODIPY reference 4aPh. 9

Table 2. Optical properties of aza-BODIPYs 4a-c in different solvents a β = Kamlet-Taft solvent basicity parameter, ET N = solvent polarity parameter normalized vs tetramethylsilane (0.00) and water (1.00). 23 b Fwhmabs, Fwhmem: full-width at half-maximum. c Data from the literature. 9

The aza-BODIPYs 4-c show a dual absorption profile that extends from 500 to 850 nm (Figure 2, Table 2). The maximum absorption band located at 750 nm for 4a, 768 nm for 4b, and 779 nm for 4c in chloroform solutions is attributed to the S0→S1 transition. The corresponding absorption coefficients are in the range of 48 400-52 700 M -1 cm -1 . A second absorption band with medium intensity is observed in the 500-650 nm region (ε = 18 600-35 600 M -1 cm -1 ) and should correspond to the S0→S2 transition, being particularly intense for the dye 4c. Furthermore, the less-intense band around 300-400 nm can be assigned to the 1-ethyl-indol-3-yl moieties. The data in Table 2 show that the 1-ethyl-indol-3-yl substituents in 4a induced a 62 nm red-shift of the maximum absorption wavelength as compared to the phenyl substituents in the reference aza-BODIPY 4aPh (λabs = 688 nm in chloroform). 9 This is attributed to the electron-donating effect of the 1-ethyl-indol-3-yl moiety. The maximum absorption wavelength of the 3,5-dithienyl-

Compound Solvent (β, ET N ) a λabs, ε (nm, M -1 cm -1 )
Fwhmabs substituted aza-BODIPY 4c is red-shifted by 29 and 11 nm in comparison with 3,5-diarylsubstituted counterparts 4a and 4b, respectively. This significant bathochromic shift for 4c can be attributed to the smaller torsion angles between the pyrrole ring plane and the thienyl ring plane, versus 4-substituted phenyl ring planes. As observed in Table 2, the absorption wavelengths of the aza-BODIPYs 4a-c are only slightly solvent-dependent. For example, upon increasing the polarity of the solvent from toluene to THF, the most noticeable solvatochromic effect was observed for 4b, bearing the less electron-donating group (C6H4Br), with a 10 nm bathochromic shift of the maximum absorption band. A systematic trend for all the compounds is a slight increase of the λabs value through the sequence toluene ≈ chloroform < THF, i.e. through increasing solvent Lewis basicity as quantified by the Kamlet-Taft parameter β, suggesting a possible stabilization of the most polar zwitterionic valence bond form of the aza-BODIPY by interaction of 23 the THF oxygen lone pair with the indolium moiety (Et-N … CH + <-:O(C2H4)2). However, the full-width at halfmaximum of absorption (Fwhmabs) of the main absorption band is also broader for 4b (around 114 nm) than for 4a and 4c in all the solvents. Excitation of the aza-BODIPYs 4a-c at 700 nm induces emission in the range 771-813 nm depending on the solvent (Table 2). In chloroform solution, the emission bands of 4a-c are localized between 781 and 805 nm (Figure 3). In this solvent, the wavelength of maximum emission of the 3,5-dithienyl-substituted aza-BODIPY 4c is red-shifted by 90 nm as compared to the reference aza-BODIPY 4aPh (λem = 715 nm in chloroform), 9 but only by 24 nm vs 4a and 1 nm vs 4b. Solvatochromism studies revealed that the emission wavelengths of aza-BODIPYs 4a-c are only slightly solvent-dependent, with small shifts of 2-12 nm for 4a, 8-17 nm for 4b, and 2-7 nm for 4c, the most notable solvatochromic effect of 17 nm being observed for 4b when changing the solvent from tetrahydrofuran to toluene. The largest Stokes shift of 36 nm occurs for 4b in polar solvents (tetrahydrofuran and chloroform). For the three compounds a marked systematic solvent effect is observed for the Stoke shift: 5-8 nm greater in THF and chloroform than in toluene (thus grossly correlating with the normalized solvent polarity parameter ET N ). 23 It is finally noteworthy that the aza-BODIPY 4a exhibits the broadest full-width half-maximum of emission (Fwhmem) around 65 nm in all the solvents.

Computational Methodology

To elucidate the origin of the optical properties in aza-BODIPYs 4a-c, computational calculations were performed with a standard protocol reported in the literature for aza-BODIPY dyes. 24 The ground state geometries of 4a-c were optimized at the PBE0/6-311G(2d,p) level of theory, followed by calculation of the harmonic vibrational frequencies to ascertain that the optimized geometries correspond to minima of energy. The electronic transitions were calculated for the optimized geometries at the TD-DFT/BMK/6-311+G(2d,p) level of approximation. The solvation by chloroform has been estimated in the calculations under the PCM model. All the calculations were performed using Gaussian 09 and the GaussView 5.0 visualization suite. 25 The calculated dihedral angles (Table 3) between the thienyl ring planes and the mean plane of the aza-BODIPY core of 4c (ϕ = 26°) are smaller, thus improving the molecular planarity in comparison with the 4-substituted phenyl derivatives (ϕ = 37° for 4a and ϕ = 39° for 4b). The dihedral angles between the indolyl planes and the mean plane of the aza-BODIPY core are calculated to be much lower at 17.7 ± 0.5°, the steric hindrance being smaller on this "north" side of the aza-BODIPY core. The contour plots of the frontier molecular orbitals (FMOs) for the aza-BODIPYs 4a-c are presented in Figure 4. Considering the TD-DFT calculated excited states, the wavelength of maximum absorption of 4a exclusively accounts for a HOMO→LUMO electronic promotion, whereas those of 4b and 4c corresponds to a combination of HOMO→LUMO and HOMO-1→LUMO excitations with a dominance of the former. The delocalization of the FMOs is a consequence of an effective electronic conjugation between the aza-BODIPY core and substituents through the 1,3,5,7-positions. Likewise, the TD-DFT calculations indicate that for all the dyes 4a-c, the absorptions in the visible range correspond to charge transfers from the substituents to the central aza-BODIPY core. Due to similar topologies of the LUMO in molecules 4a-c, it is observed that the aza-BODIPY core acts as a strong accepting unit and the 1-ethyl-indol-3-yl substituent as a strong donating unit.

The λabs max values calculated for 4a-c are presented in Table 4 and compared to the experimental values. These absorption wavelengths reflect the HOMO-LUMO gap, as evidenced by the TD-DFT studies. The calculated λabs max value was underestimated by 9 nm for 4a, 54 nm for 4b, and 42 nm for 4c as compared to the experimental data, but remains in quite good agreement. The calculated HOMO and LUMO orbital energy levels, and the values of the HOMO-LUMO gaps are given in Table 4. The HOMO and LUMO eigenvalues of the aza-BODIPYs 4a-c are located between the HOMO and LUMO eigenvalues of PCBM (HOMO = -6.0 eV, LUMO = -4.3 eV) and P3HT (HOMO = -5.1 eV, LUMO = -2.9 eV), thus making diindolyl-aza-BODIPYs potential candidates for the elaboration of ternary inverted organic solar cells. 7b The HOMO-LUMO optical gaps of 4a-c were also calculated from the absorption onset in the UV-Vis studies, and were found to be quite low, between 1.48 and 1.53 eV.

Table 3. Computed values of the torsion angles (in degrees) defining the tilt of the substituents with respect to the aza-BODIPY core plane.

Conclusion

The synthesis, absorption and emission properties of unprecedented aza-BODIPYs bearing 1-ethylindol-3-yl substituents as donating motifs at the 1,7-positions and 4-substituted phenyl and thien- 2-yl rings at the 3,5-positions, have been disclosed. By comparison with the 3,5-dianisyl-1,7diphenyl-aza-BODIPY reference, the absorption and emission bands are shifted to the nearinfrared region thanks to the 1-ethyl-indol-3-yl moieties, in agreement with their strong electrondonating character. The introduction of such electron-rich fluorophores at the 1,7-positions thus appears to be a simple method to shift the absorption and emission bands to the NIR range, by avoiding "π-frustration" 26 through the electron-accepting and dipolar character of the aza-BODIPY core. Furthermore, the thienyl groups were shown to develop stronger electronic interactions with the aza-BODIPY core than the 4-substituted phenyl counterparts, and thus to induce a redshift of the corresponding electronic spectra. TD-DFT studies allowed a detailed analysis of the effective electronic interactions between the donor periphery and the acceptor core in all the synthesized 1-ethyl-indol-3-yl-aza-BODIPY derivatives. DFT-calculated HOMO-LUMO band gap values are found to lie in the organic semiconductor range, thus opening prospects of applications in organic photovoltaics.

Experimental Section

Materials, solvents and instrumentation

All starting materials were purchased from Aldrich and used without further purification, except for 1-ethylindole-3-carboxaldehyde that was synthesized according to a literature protocol. 27 Solvents were dried using standard methods or distilled prior to use. Reactions were monitored by TLC on precoated silica gel plates (Aldrich Silica gel on TLC plates with fluorescent indicator 254 nm) and revealed by exposure to a UV lamp. FT-IR spectra were measured on a PerkinElmer Frontier FT-IR spectrophotometer. High-resolution mass spectrometry was acquired using an Agilent Technologies ESI-TOF spectrometer.

Experimental UV-visible spectra were determined using a Thermo Scientific Evolution 220 spectrophotometer.

Experimental emission spectra were recorded using a FluoroMax-4 spectrofluorometer.

Synthesis

General procedure for α,β-unsaturated ketones (1a-c)

A mixture of 1-ethylindole-3-carboxaldehyde (1 equivalent, 8.66 mmol) and the corresponding methylketone (1 equivalent, 8.66 mmol) in ethanol (20 mL) was treated with KOH (2.5 equivalents, 21.65 mmol). During the course of the reaction, the product precipitated from the reaction mixture. The precipitate was filtered and washed with water. The products were purified through column chromatography on silica gel using a n-hexane/acetone eluting mixture. 163.0, 138.0, 137.3, 132.7, 132.0, 130.5, 126.4, 123.0, 121.5, 120.9, 116.7, 113.7, 113.1, 110.2, 55.5, 41.5, 15.3. FT-IR (ATR, cm -1 ): 3093,1648,1597,1583,1559,1372,1257,1213,1161,1026,813,727,594 139.5, 137.5, 138.1, 133.3, 131.9, 130.0, 127.1, 126.4, 123.3, 121.8, 121.0, 116.4, 113.2, 110.4, 41.6, 15.3. FT-IR (ATR, cm -1 ): 3104,1647,1570,1551,1524,1394,1279,1206,1184,1007,806,726,644 4, 146.5, 138.0, 137.4, 133.0, 132.7, 130.8, 128.1, 126.4, 123.1, 121.6, 120.9, 116.5, 113.0, 110.3, 41.5, 15.2. FT-IR (ATR, cm -1 ): 3096,1640,1562,1509,1412,1273,1207,1179,981,814,709,634. ESI-TOF HRMS for C 17 H 16 NOS calculated: 282.0947, found: [M + H] + 282.0947.

General procedure for 1,3-diaryl-4-nitrobutan-1-ones (2a-c)

The α,β-unsaturated ketone 1a-c (1 equivalent, 6.55 mmol) was dissolved in a mixture of methanol/triethylamine (3:2 v/v, 30 mL) and treated with nitromethane (5 equivalents, 32.75 mmol). The reaction mixture was refluxed for 24 hours. After cooling to room temperature, the solvent was removed in vacuo and the residue was dissolved in dichloromethane and washed with a 0.5 M solution of aqueous HCl. The combined organic layers were washed with water, dried over sodium sulfate and concentrated in vacuo. The products were purified through column chromatography on silica gel using a mixture of n-hexane/acetone as eluting system.

4.2.2.1. Compound 2a. The product was obtained as a yellowish oil in 85% yield. 1 H-NMR (400.16 MHz, CDCl 3 , δ/ppm): 7. 92 -7.86 (m, 2H), 7.65 -7.60 (m, 1H), 7.32 -7.28 (m, 1H), 7.24 -7.18 (m, 1H), 7.15 -7.09 (m, 1H), 7.02 (s, 1H), 6.91 -6.85 (m, 2H), 4.92 -4.78 (m, 2H), 4.50 (p, J = 6.8 Hz, 1H), 4.07 (q, J = 7.3 Hz, 2H), 3.82 (s,3H),2H),1.39 (t,J = 7.3 Hz,3H). 13 C-NMR (100.63 MHz, CDCl 3 , δ/ppm): 196. 2, 163.8, 136.2, 130.4, 129.7, 126.5, 125.1, 121.9, 119.4, 118.7, 113.9, 112.4, 109.8, 79.2, 55.5, 41.0, 40.8, 31.7, 15.4. FT-IR (ATR, cm -1 ): 3051,2974,2934,1672,1547,1376,1212,1166,1016,816,734 196.8, 136.3, 135.4, 132.1, 129.7, 128.7, 126.4, 125.1, 122.1, 119.6, 118.7, 112.0, 109.9, 79.1, 41.2, 41.1, 31.6, 15.4 143.9, 136.2, 134.2, 132.3, 128.3, 126.4, 125.2, 122.0, 119.5, 118.7, 111.9, 109.9, 79.0, 41.8, 41.0, 31.9, 15.4. FT-IR (ATR, cm -1 ): 3073,2965,2935,1671,1548,1374,1257,1167,1023,828,699

General procedure for the preparation of aza-dipyrromethenes (3a-c)

A mixture of the corresponding 1,3-diaryl-4-nitrobutan-1-ones 2a-c (1 equivalent, 2.33 mmol) and ammonium acetate (30 equivalents, 69.89 mmol) in butanol (15 mL) was refluxed overnight. After cooling to room temperature and concentration to 2 mL, the precipitate was filtered and washed with methanol. The products were purified through column chromatography on silica gel using dichloromethane as eluent. 1, 150.2, 137.8, 136.9, 130.6, 128.5, 127.8, 126.0, 122.4, 121.5, 120.7, 115.1, 111.4, 110.3, 110.2, 56.0, 41.9, 15.6 2, 143.9, 142.9, 132.8, 130.9, 130.4, 128.4, 128.2, 127.7, 122.6, 121.4, 121.0, 111.9, 110.4, 110.1, 42.0, 15.6 137.9, 136.9, 130.7, 129.2, 128.7, 127.7, 127.3, 122.5, 121.4, 120.9, 111.8, 110.4, 110.0, 41.9, 15.6

General procedure for the synthesis of aza-BODIPYs (4a-c)

A solution of the corresponding aza-dipyrromethene 3a-c (1 equivalent, 0.29 mmol) in dry dichloromethane (20 mL) was treated with triethylamine (2.2 equivalents, 0.64 mmol). The resulting mixture was stirred at room temperature for 10 minutes, followed by addition of BF 3 •OEt 2 (4.4 equivalents, 1.27 mmol). The reaction mixture was refluxed overnight. After cooling to room temperature, the reaction mixture was washed with water. The combined organic layers were washed with brine, dried over sodium sulfate and concentrated in vacuo. The products were purified through column chromatography on silica gel using a mixture of n-hexane/ dichloromethane as eluting system. 158.0, 145.8, 139.0, 137.2, 132.0, 131.7, 127.6, 125.4, 123.1, 121.6, 121.5, 115.4, 114.4, 110.7, 110.0, 56.0, 42.3, 15.7. 19 F-NMR (376.49 6, 140.1, 137.3, 132.5, 132.1, 131.9, 131.4, 128.4, 127.5, 125.2, 123.4, 121.9, 121.5, 115.4, 110.9, 110.2, 42.4, 15.6. 19 F-NMR (376.49 MHz, CD 2 Cl 2 , δ/ppm): -129.44 (q, J(B-F) = 31.9 Hz). 11 B-NMR (128.39 MHz, CD 2 Cl 2 , δ/ppm): 1.00 (t, J(B-F) = 31.9 Hz). FT-IR (ATR, cm -1 ): 3047,2973,1570,1509,1453,1425,1371,1196,1003,798,737. ESI-TOF HRMS for C 40 H 30 BBr 2 F 2 N 5 calculated: 787.0929, found: [M] + 787.0921.

Figure 1 .

 1 Figure 1. Aza-BODIPYs 4a-4c reported in this work (data in chloroform).

Scheme 1 .

 1 Scheme 1. Synthesis of 1-ethyl-indol-3-yl-substituted aza-BODIPYs 4a-c. (i) KOH, EtOH, rt, (ii) CH 3 NO 2 , MeOH/NEt 3 , reflux, (iii) NH 4 OAc, BuOH, reflux, (iv) NEt 3 , BF 3 •Et 2 O, CH 2 Cl 2 , reflux

Figure 2 .

 2 Figure 2. Normalized absorption spectra of the aza-BODIPYs 4a-c in chloroform

Figure 3 .

 3 Figure 3. Normalized emission spectra of the aza-BODIPYs 4a-c in chloroform

ac

  Calculated with PBE0/6-311G(2d,p) level of theory. b Optical band gap estimated from the onset wavelength of the UV-Vis spectra in CHCl 3 . Calculated with TD-DFT/BMK/6-311+G(2d,p) level of theory. f = oscillator strength. d Experimental maximum absorption wavelength in CHCl 3 .

Figure 4 .

 4 Figure 4. Contour plots of the frontier molecular orbitals for the aza-BODIPYs 4a-c (PBE0/6-311G(2d,p))

Table

). These compounds showed a typical quartet 19 F NMR signal and a triplet 11 B NMR signal due to 11 B-19 F coupling. The 19 F signals appeared between -129.44 ppm and -135.82 ppm, while 11 B signals appeared between 1.00 ppm and 1.27 ppm.

Table 1 .

 1 19 F and 11 B NMR chemical shift values recorded in CD 2 Cl 2

	Compound 19 F signal (δ/ppm) 11 B signal (δ/ppm)
	4aPh	-131.76	0.06
	4a	-130.11	1.16
	4b	-129.44	1.00
	4c	-135.82	1.27
	2.2. Spectroscopic characteristics

Table 4 .

 4 Energies of the HOMO and the LUMO, HOMO-LUMO gaps, and comparison between TD-DFTcomputed and experimental λabs max values.

	Compound ϕ 1	ϕ 2	ϕ 3	ϕ 4
	4a	18	17	37	37
	4b	18	17	39	39
	4c		18	17	26	26
	Compound HOMO a	LUMO a	E(H-L) a	E(H-L) b optical
	(eV)	(eV)	(eV)	(eV)

  H, 13 C,19 F, and 11 B NMR spectra were recorded using a Bruker 400 spectrometer. Chemical shifts (δ/ppm) are reported relative to Si(CH 3 ) 4 , CDCl 3 , or CD 2 Cl 2 .

  4.2.1.1. Compound 1a. The product was obtained as a yellow solid in 89% yield. Melting point: 112-113 °C.1 H-NMR (400.16 MHz, CDCl 3 , δ/ppm): 8.14 -8.05 (m, 3H), 8.04 -8.00 (m, 1H), 7.57(d,

	J = 15.5 Hz, 1H), 7.51 (s, 1H), 7.43 -7.37 (m, 1H), 7.35 -7.28 (m, 2H), 7.04 -6.94 (m, 2H), 4.19 (q, J
	= 7.3 Hz, 2H), 3.89 (s, 3H), 1.51 (t, J = 7.3 Hz, 3H). 13 C-NMR (100.63 MHz, CDCl 3 , δ/ppm): 189.0,

  . ESI-TOF HRMS for C 20 H 20 NO 2 calculated: 306.1488, found: [M + H] + 306.1488.

	4.2.1.2. Compound 1b. The product was obtained as a yellow solid in 84% yield. Melting point:
	136-137 °C. 1 H-NMR (400.16 MHz, CDCl 3 , δ/ppm): 8.10 (d, J = 15.5 Hz, 1H), 8.03 -7.98 (m, 1H),
	7.95 -7.90 (m, 2H), 7.67 -7.62 (m, 2H), 7.56 (s, 1H), 7.48 (d, J = 15.5 Hz, 1H), 7.43 -7.40 (m, 1H),
	7.37 -7.29 (m, 2H), 4.23 (q, J = 7.3 Hz, 2H), 1.54 (t, J = 7.3 Hz, 3H). 13 C-NMR (100.63 MHz, CDCl 3 ,
	δ/ppm): 189.6,

  . ESI-TOF HRMS for C 19 H 17 BrNO calculated: 354.0488, found: [M + H] + 354.0486.

4.2.1.3

. Compound 1c. The product was obtained as a yellow solid in 87% yield. Melting point: 132-133 °C.

1 H-NMR (400.16

MHz, CDCl 3 , δ/ppm): 8.12 (d, J = 15.4 Hz, 1H), 8.04 -7.98 (m, 1H), 7.88 (dd, J = 3.8, 1.2 Hz, 1H), 7.63 (dd, J = 5.0, 1.1 Hz, 1H), 7.53 (s, 1H), 7.48 -7.37 (m, 2H), 7.36 -7.29 (m, 2H), 7.18 (dd, J = 4.9, 3.7 Hz, 1H), 4.19 (q, J = 7.3 Hz, 2H), 1.51 (t, J = 7.3 Hz, 3H). 13 C-NMR (100.63 MHz, CDCl 3 , δ/ppm): 182.

  . ESI-TOF HRMS for C 21 H 23 N 2 O 4 calculated: 367.1652, found: [M + H] + 367.1653.

	4.2.2.2. Compound 2b. The product was obtained as a yellowish oil in 80% yield. 1 H-NMR (400.16
	MHz, CDCl 3 , δ/ppm): 7.79 -7.76 (m, 2H), 7.64 -7.60 (m, 1H), 7.59 -7.56 (m, 2H), 7.34 -7.31 (m,
	1H), 7.26 -7.21 (m, 1H), 7.14 (ddd, J = 8.0, 7.0, 1.1 Hz, 1H), 7.02 (s, 1H), 4.92 -4.81 (m, 2H), 4.50
	(p, J = 7.0 Hz, 1H), 4.11 (q, J = 7.3 Hz, 2H), 3.61 -3.48 (m, 2H), 1.42 (t, J = 7.3 Hz, 3H). 13 C-NMR
	(100.63 MHz, CDCl 3 , δ/ppm):

  . ESI-TOF HRMS for C 18 H 19 N 2 O 3 S calculated: 343.1110, found: [M + H] + 343.1108.

  . FT-IR (ATR, cm -1 ): 3058, 2928, 1601, 1456, 1244, 1151, 1024, 900, 736. ESI-TOF HRMS for C 42 H 38 N 5 O 2 calculated: 644.3020, found: [M + H] + 644.3023.

	4.2.3.2. Compound 3b. The product was obtained as a dark blue-black solid in 25% yield. Melting
	point: >300 °C. 1 H-NMR (400.16 MHz, CD 2 Cl 2 , δ/ppm): 8.11 -8.07 (m, 2H), 7.90 -7.86 (m, 6H), 7.72
	-7.68 (m, 4H), 7.45 -7.42 (m, 2H), 7.32 -7.27 (m, 2H), 7.24 -7.19 (m, 4H), 4.09 (q, J = 7.3 Hz, 4H),
	1.37 (t, J = 7.3 Hz, 6H). 13 C-NMR (100.63 MHz, CD 2 Cl 2 , δ/ppm): 153.5, 150.

  . FT-IR (ATR, cm -1 ):3059, 2921, 1588, 1531, 1453, 1233, 1006, 901, 788, 732. ESI-TOF HRMS for C 40 H 32 Br 2 N 5 calculated: 740.1018, found: [M + H] + 740.1016.

	4.2.3.3. Compound 3c. The product was obtained as a dark blue-black solid in 28% yield. Melting
	point: >300 °C. 1 H-NMR (400.16 MHz, CD 2 Cl 2 , δ/ppm: 8.08 (d, J = 7.9 Hz, 2H), 7.88 (s, 2H), 7.70 -
	7.66 (m, 2H), 7.56 -7.51 (m, 2H), 7.46 -7.40 (m, 2H), 7.32 -7.18 (m, 6H), 7.10 (s, 2H), 4.09 (q, J =
	7.3 Hz, 4H), 1.37 (t, J = 7.3 Hz, 6H). 13 C-NMR (100.63 MHz, CD 2 Cl 2 , δ/ppm): 150.2, 150.0, 138.1,

  The product was obtained as a dark blue-violet solid in yield. Melting point: >300 °C.1 H-NMR (400.16 MHz, CD 2 Cl 2 , δ/ppm): 8.15 -8.11 (m, 2H), 8.10 -8.06 (m, 4H), 8.04 (s, 2H), 7.48 -7.44 (m, 2H), 7.35 -7.30 (m, 2H), 7.28 -7.23 (m, 2H), 7.07 -7.02 (m, 6H), 4.19 (q, J = 7.3 Hz, 4H), 3.90 (s, 6H), 1.46 (t, J = 7.3 Hz, 6H). 13 C-NMR (100.63 MHz, CD 2 Cl 2 , δ/ppm): 162.0,

4.2.4.1. Compound 4a.

  MHz, CD 2 Cl 2 , δ/ppm): -130.11 (q, J(B-F) = 32.3 Hz). 11 B-NMR (128.39 MHz, CD 2 Cl 2 , δ/ppm): 1.16 (t, J(B-F) = 32.3 Hz). FT-IR (ATR, cm -1 ):3044, 2928, 1601, 1561, 1453, 1433, 1380, 1077, 966, 735. ESI-TOF HRMS for C 42 H 37 BF 2 N 5 O 2 calculated: 692. 3002, found: [M + H] + 692.3003.4.2.4.2. Compound 4b.The product was obtained as a dark blue-violet solid in 79% yield. Melting point: >300 °C.1 H-NMR (400.16 MHz, CD 2 Cl 2 , δ/ppm): 8.14 -8.10 (m, 2H), 8.06 (s, 2H), 7.96 -7.91 (m, 4H), 7.66 -7.62 (m, 4H), 7.49 -7.45 (m, 2H), 7.36 -7.24 (m, 4H), 7.02 (s, 2H), 4.20 (q, J = 7.3 Hz, 4H), 1.46 (t, J = 7.3 Hz, 6H). 13 C-NMR (100.63 MHz, CD 2 Cl 2 , δ/ppm): 157.

Acknowledgements

L. L. thanks CONACYT for the grant (576732) supporting his doctoral studies. The authors acknowledge the financial support from CONACYT, PAIP, PAPIIT (IN 222819), and the French-Mexican International Laboratory (LIA-LCMMC-CONACYT) programs.

Declaration of interests

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.


Manon Lincker
email: manon.lincker@mines-paristech.fr

Jennifer L Druhan

Sophie Guillon

Vincent Lagneau

Thermodynamic and kinetic constants for isotopic fractionation modeling with or without major isotope hypothesis

Keywords: geochemical modeling, isotopic fractionation, major isotope hypothesis

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Geochemical speciation and reactive transport (RT) codes were originally developed several decades ago to provide a means of considering the long term implications and management of nuclear waste disposal [START_REF] Lichtner | Multiphase multicomponent nonisothermal reactive transport in partially saturated porous media[END_REF], [START_REF] Steefel | Multicomponent reactive transport in discrete fractures: I. controls on reaction front geometry[END_REF], [START_REF] Steefel | Multicomponent reactive transport in discrete fractures: Ii: Infiltration of hyperalkaline groundwater at maqarin, jordan, a natural analogue site[END_REF]. Since then, these basic principles have been broadly expanded to other purposes like pollution remediation ( [START_REF] Wanner | A chromate-contaminated site in southern switzerland-part 2: Reactive transport modeling to optimize remediation options[END_REF], [START_REF] Wanner | Benchmarking the simulation of cr isotope fractionation[END_REF], [START_REF] De Windt | Modelling of long-term dynamic leaching tests applied to solidified/stabilised waste[END_REF], [START_REF] Avasarala | Reactive transport of u and v from abandoned uranium mine wastes[END_REF], [START_REF] Li | Mineral transformation and biomass accumulation associated with uranium bioremediation at rifle, colorado[END_REF]), mining [START_REF] Lagneau | Industrial deployment of reactive transport simulation: an application to uranium in situ recovery[END_REF], resources exploitation [START_REF] Sin | 2d simulation of natural gas reservoir by two-phase multicomponent reactive flow and transport -description of a benchmarking exercise[END_REF] and water-rock interactions [START_REF] Maher | Reactive transport processes that drive chemical weathering: From making space for water to dismantling continents[END_REF].

Over a comparable period of time, stable isotope data have been increasingly utilized [START_REF] Urey | The thermodynamic properties of isotopic substances[END_REF] for identification and quantification of water-rock interactions ( [START_REF] Banner | Water-rock interaction history of regionally extensive dolomites of the burlington-keokuk formation (mississippian): isotopic evidence[END_REF], [START_REF] Banner | Calculation of simultaneous isotopic and trace element variations during water-rock interaction with applications to carbonate diagenesis[END_REF], [START_REF] Fisher | Water-rock interaction in tertiary sandstones, san joaquin basin, california, usa: diagenetic controls on water composition[END_REF], ...), weathering ( [START_REF] Dasch | Strontium isotopes in weathering profiles, deep-sea sediments, and sedimentary rocks[END_REF], [START_REF] Pistiner | Lithium-isotope fractionation during continental weathering processes[END_REF], [START_REF] Lawrence | Hydrogen and oxygen isotope systematics in weathering profiles[END_REF], [START_REF] Georg | Silicon isotope variations accompanying basalt weathering in iceland[END_REF], [START_REF] Hatton | Investigation of subglacial weathering under the greenland ice sheet using silicon isotopes[END_REF], ...), vegetation and rhizosphere cycling of water, soil or minerals ( [START_REF] Cenki-Tok | The impact of water-rock interaction and vegetation on calcium isotope fractionation in soil-and stream waters of a small, forested catchment (the strengbach case)[END_REF], [START_REF] Sprenger | Illuminating hydrological processes at the soil-vegetationatmosphere interface with water stable isotopes[END_REF], ...), pollution ( [START_REF] Farmer | Stable lead isotope record of lead pollution in loch lomond sediments since 1630 ad[END_REF], [START_REF] Cloquet | Tracing source pollution in soils using cadmium and lead isotopes[END_REF], [START_REF] Sherman | The use of pb, sr, and hg isotopes in great lakes precipitation as a tool for pollution source attribution[END_REF], [START_REF] Nestler | Isotopes for improved management of nitrate pollution in aqueous resources: review of surface water field studies[END_REF], [START_REF] Rogers | Stable carbon and nitrogen isotope signatures indicate recovery of marine biota from sewage pollution at moa point, new zealand[END_REF], ...), contaminant degradation pathways such as chlorinated hydrocarbons ( [START_REF] Wanner | Quantification of degradation of chlorinated hydrocarbons in saturated low permeability sediments using compound-specific isotope analysis[END_REF], [START_REF] Hunkeler | Monitoring oxidation of chlorinated ethenes by permanganate in groundwater using stable isotopes: Laboratory and field studies[END_REF], [START_REF] Hunkeler | Assessing chlorinated ethene degradation in a large scale contaminant plume by dual carbon-chlorine isotope analysis and quantitative pcr[END_REF],...), oceanography ( [START_REF] Goldstein | Long-lived isotopic tracers in oceanography, paleoceanography, and ice-sheet dynamics[END_REF], [START_REF] Sigman | Ocean process tracers: nitrogen isotopes in the ocean[END_REF], [START_REF] Almogi-Labin | Quaternary paleo-oceanography, pteropod preservation and stableisotope record of the red sea[END_REF], [START_REF] Brand | Geochemistry of modern brachiopods: applications and implications for oceanography and paleoceanography[END_REF], ...), paleoreconstruction ( [START_REF] Lemarchand | Boron isotope systematics in large rivers: implications for the marine boron budget and paleo-ph reconstruction over the cenozoic[END_REF], [START_REF] Rozanski | Reconstruction of past climates from stable isotope records of palaeo-precipitation preserved in continental archives[END_REF], [START_REF] Thomas | Neodymium isotopic reconstruction of late paleocene-early eocene thermohaline circulation[END_REF], [START_REF] Gázquez | Quantification of paleo-aquifer changes using clumped isotopes in subaqueous carbonate speleothems[END_REF], [START_REF] Frei | A systematic look at chromium isotopes in modern shells-implications for paleo-environmental reconstructions[END_REF], ...) and biogeochemical cycling ( [START_REF] Group | Geotraces-an international study of the global marine biogeochemical cycles of trace elements and their isotopes[END_REF], [START_REF] Hindshaw | Experimental constraints on mg isotope fractionation during clay formation: Implications for the global biogeochemical cycle of mg[END_REF], [START_REF] Schmitt | Calcium biogeochemical cycle at the beech tree-soil solution interface from the strengbach czo (ne france): insights from stable ca and radiogenic sr isotopes[END_REF], [START_REF] Denk | The nitrogen cycle: A review of isotope effects and isotope modeling approaches[END_REF], [START_REF] Abadie | Iron isotopes reveal distinct dissolved iron sources and pathways in the intermediate versus deep southern ocean[END_REF], ...) .

RT software are now increasingly used to help interpret and fully leverage these isotopic data [START_REF] Druhan | Stable isotope fractionation by transport and transformation[END_REF].

As isotopes and fractionating pathways continue to be incorporated and refined into reactive transport simulators ( [START_REF] Druhan | A large column analog experiment of stable isotope variations during reactive transport: I. a comprehensive model of sulfur cycling and δ34s fractionation[END_REF], [START_REF] Druhan | A large column analog experiment of stable isotope variations during reactive transport: Ii. carbon mass balance, microbial community structure and predation[END_REF], [START_REF] Druhan | The influence of mixing on stable isotope ratios in porous media: A revised rayleigh model[END_REF], [START_REF] Wanner | Benchmarking the simulation of cr isotope fractionation[END_REF], [START_REF] Wanner | Assessing the control on the effective kinetic cr isotope fractionation factor: A reactive transport modeling approach[END_REF], [START_REF] Steefel | Reactive transport codes for subsurface environmental simulation[END_REF], [START_REF] Druhan | Stable and radioactive carbon isotope partitioning in soils and saturated systems: a reactive transport modeling benchmark study[END_REF]) their capacity to allow a deeper and more quantitative understanding of reactive transformations and system dynamics is becoming increasingly evident.

Reactive transport models which explicitly include stable isotope fractionation have been developed for a variety of geoscience and environmental engineering applications including contamination( [START_REF] Wanner | Benchmarking the simulation of cr isotope fractionation[END_REF], [START_REF] Van Breukelen | Reactive transport modelling of biogeochemical processes and carbon isotope geochemistry inside a landfill leachate plume[END_REF], [START_REF] Van Breukelen | Quantification of sequential chlorinated ethene degradation by use of a reactive transport model incorporating isotope fractionation[END_REF]), unsaturated zone water dynamics [START_REF] Singleton | Multiphase reactive transport modeling of seasonal infiltration events and stable isotope fractionation in unsaturated zone pore water and vapor at the hanford site[END_REF], flow and transport through heterogeneous porous media [START_REF] Druhan | A model linking stable isotope fractionation to water flux and transit times in heterogeneous porous media[END_REF], experimental analogs ( [START_REF] Druhan | A large column analog experiment of stable isotope variations during reactive transport: I. a comprehensive model of sulfur cycling and δ34s fractionation[END_REF], [START_REF] Druhan | A large column analog experiment of stable isotope variations during reactive transport: Ii. carbon mass balance, microbial community structure and predation[END_REF]), marine sediments pore fluids [START_REF] Maher | The mineral dissolution rate conundrum: Insights from reactive transport modeling of u isotopes and pore fluid chemistry in marine sediments[END_REF] and groundwater treatment systems [START_REF] Gibson | 34s/32s fractionation during sulfate reduction in groundwater treatment systems: Reactive transport modeling[END_REF], etc. These applications require a rigorous integration of isotopic and geochemical data for modeling purposes.

Geochemical speciation codes are commonly based on the primary component or basis species approach [START_REF] Lichtner | Continuum formulation of multicomponent-multiphase reactive transport[END_REF], in which the total number of chemical species tracked during transport is limited to the number that are linearly independent from one another. Practically speaking, there is no need to waste computational effort tracking the advection of both bicarbonate and carbonic acid if they are coupled through an equilibrium relationship. Implementation of isotopic fractionation then means that each isotope must be created as a primary component, with all necessary attendant equilibria described. This exercise is usually performed by duplicating entries within the thermodynamic database, so that each isotopologue is considered as a fully-fledged chemical species ( [START_REF] Druhan | Stable isotope fractionation by transport and transformation[END_REF], [START_REF] Druhan | Calcium isotope fractionation in groundwater: Molecular scale processes influencing field scale behavior[END_REF], [START_REF] Wanner | Benchmarking the simulation of cr isotope fractionation[END_REF], [START_REF] Druhan | A model linking stable isotope fractionation to water flux and transit times in heterogeneous porous media[END_REF], [START_REF] Druhan | The influence of mixing on stable isotope ratios in porous media: A revised rayleigh model[END_REF]).

Kinetically controlled reactions are modeled by reactive transport codes by defining the parameters of a selected kinetic law linking two or more components. When modeling kinetic isotopic fractionation, the common strategy is consistent with equilibrium reactions and involves duplication of isotopes as unique 'species' in the database so that they may be assigned different rate constants depending on the isotopologue involved in the reaction to reproduce kinetic isotopic fractionation. This approach brings the significant advantage that the basic framework of the pre-existing RT codes remains unchanged, requiring only appropriate identification of parameter values. Typically, this is accomplished by reasonably assuming that the thermodynamic and/or rate constant used for the overall or 'bulk' formation of a given compound are effectively equivalent to those describing the most abundant isotope. For example, the equilibrium constant describing the abundance of bicarbonate and water through the dissociation of carbonic acid would be used to describe the 12 C isotopologue of these compounds. The thermodynamic or rate constants of the rarer isotope(s) (e.g. 13 C) are then determined based on reported fractionation factors [START_REF] Druhan | Stable isotope fractionation by transport and transformation[END_REF] appropriate for each equilibrium or kinetic reaction at a given temperature.

As the capacity for highly precise mass ratio measurements continues to expand, an increasing number of isotope systems are being developed in application to geosciences and environmental engineering studies.

Unlike the 'traditional' isotope systems (e.g. H, C, N, O, S), many of these elements do not present a major isotope, for instance Cl [START_REF] Hendry | Chloride and chlorine isotopes (36cl and δ37cl) as tracers of solute migration in a thick, clay-rich aquitard system[END_REF], [START_REF] Hunkeler | Assessing chlorinated ethene degradation in a large scale contaminant plume by dual carbon-chlorine isotope analysis and quantitative pcr[END_REF], Mo [START_REF] Migeon | Molybdenum isotope fractionation during acid leaching of a granitic uranium ore[END_REF], [START_REF] Barling | Molybdenum isotope fractionation during adsorption by manganese oxides[END_REF], Br [START_REF] Shouakar-Stash | Geochemistry and stable isotopic signatures, including chlorine and bromine isotopes, of the deep groundwaters of the siberian platform, russia[END_REF], Cu [START_REF] Larson | Copper isotope ratios in magmatic and hydrothermal ore-forming environments[END_REF], [START_REF] Mathur | Exploration potential of cu isotope fractionation in porphyry copper deposits[END_REF], Zn [START_REF] Cloquet | Variation in the isotopic composition of zinc in the natural environment and the use of zinc isotopes in biogeosciences: a review[END_REF], [START_REF] Pichat | Zinc isotope variations in deep-sea carbonates from the eastern equatorial pacific over the last 175 ka[END_REF], etc. The issue of the major isotope was considered several decades ago by Thorstenson and Parkhurst [START_REF] Thorstenson | Calculation of individual isotope equilibrium constants for geochemical reactions[END_REF] , who offered direct derivations of the individual isotope equilibrium constants of some species, e.g. CO 2(g) , CO 2(aq) ,

HCO - 3 , CO 2- 3 , CaCO 3(aq) , H 2 O (g) , H 2 O (l) , H 3 O + , OH -, HCO - 3 and CO 2- 3
in ions pairs with Na + , K + , Ca 2+ and Mg 2+ following the earlier approach developed by [START_REF] Urey | The thermodynamic properties of isotopic substances[END_REF] [START_REF] Urey | The thermodynamic properties of isotopic substances[END_REF]. While the methodology does not require the definition of a major isotope, it cannot be readily generalized and requires an accounting of the degree of symmetry of exchanged atoms which the authors only provided for the aforementioned species [START_REF] Thorstenson | Calculation of individual isotope equilibrium constants for geochemical reactions[END_REF]. Therefore, continued expansion of isotope-enabled RT applications suggests that the common method by which isotopic partitioning is implemented in these models requires careful consideration, and, where necessary, a more generalized approach such that the benefits of RT codes can be applied to emerging isotopic compounds that do not have an obvious choice of major isotope.

The present study addresses this need by revisiting the classical mathematical formulation for equilibrium and kinetic isotopic fractionation as typically implemented in RT codes to produce a generalized method relaxing the constraint associated with the need for a major isotope. We will first walk through a series of modifications to the formal approach of equilibrium fractionation, then consider the extension to kinetic fractionation and finally offer a demonstration in the context of zinc isotopes in Mississippi Valley Type (MVT) ore deposits.

Classical method for equilibrium isotopic fractionation modeling with RT codes

Reactive transport codes developed following the basis species approach [START_REF] Lichtner | Continuum formulation of multicomponent-multiphase reactive transport[END_REF] consider isotopologues as individual components with their own thermodynamic constants. This means that two isotopes of a given element in a given compound exist as two separate 'species' that operate in tandem through parallel reactions and transport pathways. One of the most basic ways in which these 'species' are tied together as isotopes of the same element is through the slight differences in their thermodynamic constants, which stems from mass-dependent equilibrium fractionation. We begin with the mathematical development commonly used to determine the thermodynamic constants for reactions involving different isotopologues.

Duplication of isotopologues in the database

In the following equations, E i stands for a basis species (mol.kg -1 ), Ē for a derived species (mol.kg -1 ), γ for the activity coefficient, N c is the number of basis components and ν represents the stoichiometric coefficients. We hereby present the general formalism, which is appropriate to any isotope system of interest.

For the sake of clarity, we also offer an example: isotopes of sulfur in anhydrite and sulfates (where 32 S can be considered as a major isotope and 34 S as a minor isotope of sulfur). Critically, this choice is intended only for the sake of clarity and does not restrict the derivation to this context.

Let us consider the following chemical reaction :

Nc-1 i=0 ν i E i Ē, (1)
and its associated equilibrium constant:

K = γ • Ē i γ νi i E νi i (2)
For the specific case of CaSO 4 formation, it becomes :

Ca 2+ + SO 2- 4 CaSO 4(s) , (3)
K CaSO 4(s) = γ CaSO 4(s) CaSO 4(s) γ Ca 2+ Ca 2+ • γ SO 2- 4 SO 2- 4 . (4)
When multiple isotopes of one element are considered, they are denoted through an isotope index λ. For example, E 0 and E λ or for our particular case study 32 S and 34 S. We represent E 0 as the basis species which carries the isotopes, namely SO 2- 4 for the S isotopes example. In equation [START_REF] Abadie | Iron isotopes reveal distinct dissolved iron sources and pathways in the intermediate versus deep southern ocean[END_REF], this carrying species is isolated out of the sum leading to:

ν 0 E 0 + Nc-1 i=1 ν i E i Ē (5)
Let us assume only two isotopes for the moment: the major (most abundant one, e.g. 32 S) and a minor one (least abundant, e.g. 34 S). Furthermore, let us assume here that there is only one isotopically substituted element in the derived species. E 0 0 is substituted only once by E λ 0 in Ēλ :

(ν 0 -1)E 0 0 + E λ 0 + Nc-1 i=1 ν i E i Ēλ . (6)
The equilibrium constant can be written as:

K = γ Ēλ γ ν0 0 (E 0(ν0-1) 0 • E λ 0 ) • Nc-1 i=1 γ νi i E νi i , (7)
In the case of our example, only one atom of sulfur is considered, thus the stoichiometry of SO 2- 4 is 1

(ν 0 = 1). We introduce N iso the number of different isotopes considered and the isotopic ratios R λ = Ēλ / Ē0 and R λ 0 = Ēλ 0 / Ē0 0 . Substituting R λ and R λ 0 in (2) isolates the thermodynamic constant of the major isotope:

K = K 0 1 + Niso-1 λ=1 R λ 1 + Niso-1 λ=1 R λ 0 ν0 . (8)
Likewise for the CaSO 4(s) example:

K = K 32 1 + R 34 CaSO 4(s) 1 + R 34 SO 2- 4 (9)

Simplification based on major isotope hypothesis

Where the major isotope assumption holds true, the isotopic ratios are very small (compared to 1) and it is reasonable to assume that the thermodynamic constant K of the major isotope (K 32 in our example) equals the generally known constant for the bulk chemical species, as is often already available in geochemical databases. Physically, this leads to the assumption that there is one isotope whose abundance is predominant compared to others. This assumption, which is valid for many elements (like C, N, O, Ca, K, S, Fe, Li, ...), yields a condition under which R λ << 1 and R λ 0 << 1 so that λ R λ << 1 and λ R λ 0 << 1. As a consequence, :

K 0 K (10)
For our CaSO 4(s) example, we can safely assume that R 34 CaSO 4(s) << 1 and R 34 SO 2-

4

<< 1. This leads to

K 32 K (11)
The equilibrium fractionation factor is defined as follows:

α λ B-A = R λ B R λ A ( 12
)
with A standing for the reactant, B for the instantaneous product and λ the isotope. Using the fractionation factor between two species, the equilibrium constant for the reaction which involves a minor isotope can be written as a function of the equilibrium constant of the major isotope reaction :

α λ 0 = R λ R λ 0 = Ēλ Ē0 E λ 0 E 0 0 = K λ K 0 (13)
so that

K λ α λ 0 K. ( 14
)
Continuing with CaSO 4 example:

α 34 CaSO 4(s) -SO 2- 4 = K 34 K 32 (15)
and

K 34 α 34 CaSO 4(s) -SO 2- 4 • K (16)

Limitations

We now expand our consideration to the potential for error introduced when such an approach is applied to elements without significant difference in abundance(s) between the major and minor isotope(s) (Ni, Cu, Zn, Mg, Cl, Br, etc). These concerns motivate a generalization of the calculations described above for thermodynamic constants without the need to invoke a major isotope hypothesis.

Calculation: Circumvention of the major isotope hypothesis

Simple substitution

Here we propose to expand the calculations described above to avoid the use of the major isotope assumption. Instead of the major isotope hypothesis, we rely on characteristically small variations of enrichment factor (defined as 1 -α) and R. As a secure upper limit of enrichment in natural systems, we can assume that is smaller than 100 (i.e. α between 0.90 and 1.1). This is generally a safe assumption for the values of low-(e.g. C, N) and mid-mass (e.g. Ca) isotope systems [START_REF] Mook | Carbon isotope fractionation between dissolved bicarbonate and gaseous carbon-dioxide[END_REF], [START_REF] Bottinga | Calculated fractionation factors for carbon and hydrogen isotope exchange in the system calcite-carbon dioxide-graphite-methane-hydrogen-water vapor[END_REF], [START_REF] Oelkers | Using stable mg isotope signatures to assess the fate of magnesium during the in situ mineralisation of co2 and h2s at the carbfix site in sw-iceland[END_REF], [START_REF] Chang | Mg and ca isotope fractionation during caco3 biomineralisation[END_REF], [START_REF] Druhan | Calcium isotope fractionation in groundwater: Molecular scale processes influencing field scale behavior[END_REF] and even holds true for large differences in mass among the isotopes of a given element, as in the case of Li [START_REF] Bouchez | Modeling novel stable isotope ratios in the weathering zone[END_REF]. In the following mathematical developments, the second order and higher terms in R will therefore be neglected.

Once again we consider the following reaction where element E 0 has two isotopologues E 0 0 and E λ 0 :

E 0 + Nc-1 i=1 ν i E i Ē. ( 17
)
Therefore, we can define the following isotopic ratios R, fractionation factor α and enrichment factor :

R λ 0 = E λ 0 E 0 0 (18) Rλ = Ēλ Ē0 (19)
α λ = Rλ R λ 0 = 1 -λ . ( 20
)
The derivation of equilibrium constant K leads to:

K = γ Ē γ 0 E 0 • i (γ i E i ) νi = γ Ē0 • 1 + Ēλ Ē0 γ 0 E 0 • i (γ i E i ) νi • 1 + E λ 0 E 0 0 = K 0 • 1 + α λ R λ 0 1 + R λ 0 = K 0 • 1 + R λ 0 + λ R λ 0 1 + R λ 0 ( 21
)
where K 0 is the equilibrium constant of the reaction involving E 0 0 as E 0 isotopologue.

Thus

K 0 = K • 1 + R λ 0 1 + R λ 0 + λ R λ 0 (22)
and

K λ = α λ • K • 1 + R λ 0 1 + R λ 0 + λ R λ 0 . (23)
In the specific case where E 0 0 is a major isotope, this result is identical to section 2.

We call now T 0 the total of species E 0 such as T 0 = Ē + E 0 . Using the equilibrium constant K to express E 0 leads to:

E 0 = T 0 1 + K γ 0 γ i (γ i E i ) νi . (24)
We next investigate the difference in accuracy when using the simplification derived from the major isotope hypothesis but without making any assumption about the relative abundance of any isotopologues.

For this purpose, we will use the precise expression of the equilibrium constants of the reactions involving isotopologues E 0 0 and E λ 0 respectively (equations 22 and 23) and approximation of these constants: K 0 = K and K λ = α λ K, corresponding to equations 10 and 14. To evaluate these approximations we monitor the isotopic ratio Rλ 0 when using approximations relative to the isotopic ratio R λ 0 when using exact equilibrium constants (equations [START_REF] Druhan | Stable and radioactive carbon isotope partitioning in soils and saturated systems: a reactive transport modeling benchmark study[END_REF][START_REF] Druhan | A model linking stable isotope fractionation to water flux and transit times in heterogeneous porous media[END_REF]. This leads to:

R = Rλ 0 R λ 0 = Êλ 0 Ê0 0 E λ 0 E 0 0 = 1 + K γ 0 γ i (γ i E i ) νi T 0 0 • T λ 0 1 + α λ K γ 0 γ i (γ i E i ) νi 1 + K 0 γ 0 γ i (γ i E i ) νi T 0 0 • T λ 0 1 + K λ γ 0 γ i (γ i E i ) νi = 1 + K 0 γ 0 γ i (γ i E i ) νi 1 + K γ 0 γ i (γ i E i ) νi • 1 + α λ K γ 0 γ i (γ i E i ) νi 1 + K λ γ 0 γ i (γ i E i ) νi = 1 + K 0 Z 1 + KZ • 1 + α λ KZ 1 + K λ Z (25)
where

Z = γ0 γ i (γ i E i ) νi .
Using Taylor expansion on equations 22 and 23:

K 0 = K 1 + λ R λ 0 1 + R + o λ R 2 (26)
and

K λ = α λ K 1 + λ R λ 0 1 + R + o λ R 2 . ( 27
)
Then, incorporating 26 and 27 in 25:

R = 1 + KZ 1 + α λ + KZ λ R λ 0 1 + R λ 0 + α λ K 2 Z 2 1 + λ R λ 0 1 + R λ 0 1 + KZ (1 + α λ ) + α λ KZ λ R λ 0 1 + R λ 0 + α λ K 2 Z 2 1 + λ R λ 0 1 + R λ 0 . (28)
However,

α λ KZ λ R λ 0 1 + R λ 0 = 1 -λ KZ λ R λ 0 1 + R λ 0 = KZ λ R λ 0 1 + R λ 0 + o λ R λ 0 2 , (29)
and thus,

R = 1 + o λ R λ 0 2 . ( 30
)
This result shows that the error introduced by using the approximation is of the order o λ R λ 0 2 .

Limitations

Results show that the error on the isotopic ratio deriving from the use of the approximation K 0 = K and K λ = αK is in the order of ( R) 2 . With usually less (or much less) than 0.1, the assumption on thermodynamic constant can be safely in most natural applications. The ultra-major isotope has already been presented (see subsection Simplification based on major isotope hypothesis). For less discriminating isotopic ratios, we stressed that the choice of reference isotope was arbitrary; therefore, when one isotope has a clearly higher relative abundance than the other, the smartest choice is to choose it as a reference (with resulting R < 1, the error is further reduced). This still holds true for isotopes of similar concentrations, with R 1. In all cases, it is probably safer to check whether R is indeed small enough for the acceptable error. Finally, only one application might lead to unacceptable errors: when the most abundant isotope changes over time or space in the same system (e.g. R 0.2 to 0.8), and for larger values of , although the authors cannot imagine natural cases falling into these conditions.

Extension to double substitution

We next extend this methodology to address the case of double substitution. We consider here the following reaction:

E 0 + E 1 + N c-2 i=2 ν i E i Ē (31)
where E 0 and E 1 have two possible isotopes (0 and λ).

We suppose that:

Rλ = Ēλ Ē0 = Ēλ0 Ē00 = Ēλλ Ē0λ (32)
and

R λ = Ē λ Ē 0 = Ē0λ Ē00 = Ēλλ Ēλ0 . (33)
As for the case of one element bringing isotopes in the reaction, activity coefficients and E i terms will be removed from the calculation of the ratio as they will be divided by themselves.

The equilibrium constant can be written as:

K = Ē00 + Ē0λ + Ēλ0 + Ēλλ E 0 0 + E λ 0 E 0 1 + E λ 1 . (34)
We can then derive the specific equilibrium constants:

K = K 00 1 + R λ + Rλ + R λ Rλ 1 + R λ 0 1 + R λ 1 ( 35
)
where

R λ 0 = E λ 0 E 0 0 and R λ 1 = E λ 1 E 0 1 .
Thus,

K 00 = K 1 + R λ 0 1 + R λ 1 1 + α λ 0 R λ 0 + α λ 1 R λ 1 + α λ 0 α λ 1 R λ 0 R λ 1 . (36)
Using the same methodology, it appears that:

K λ0 = α λ 0 • K 1 + R λ 0 1 + R λ 1 1 + α λ 0 R λ 0 + α λ 1 R λ 1 + α λ 0 α λ 1 R λ 0 R λ 1 . (37)
We will now focus on E 0 and its isotopes, assuming that the E 1 isotopologue involved in the reaction is

E 0 1 .
We can write the following reactions:

E 0 0 + E 0 1 Ē00 (38)
E λ 0 + E 0 1 Ēλ0 . ( 39
)
Considering T 0 0 the total of species E 0 0 and T λ 0 the total of species E λ 0 :

E 0 0 = T 0 0 1 + K 00 E 0 1 (40)
E λ 0 = T λ 0 1 + K λ0 E 0 1 . (41)
As for the case of one species bringing isotopes in the system, we track the relative difference in the isotopic ratio when using K and α λ 0 for approximation relative to the isotopic ratio when using constants

K 00 and K λ0 . R = R λ 0 Rλ 0 = E λ 0 E 0 0 Êλ 0 Ê0 0 = 1 + K 00 E 0 1 1 + K λ0 E λ 0 • 1 + α λ 0 KE 0 1 1 + KE 0 1 = 1 + E 0 1 K 1 + R λ 0 1 + R λ 1 1 + α λ 0 R λ 0 + α λ 1 R λ 1 + α λ 0 α λ 1 R λ 0 R λ 1 1 + α λ 0 E 0 1 K 1 + R λ 0 1 + R λ 1 1 + α λ 0 R λ 0 + α λ 1 R λ 1 + α λ 0 α λ 1 R λ 0 R λ 1 • 1 + α λ 0 KE 0 1 1 + KE 0 1 . (42)
Noting that

1 + α λ 0 R λ 0 + α λ 1 R λ 1 + α λ 0 α λ 1 R λ 0 R λ 1 1 + R λ 0 1 + R λ 1 = 1 - λ 0 R λ 0 1 + R λ 1 + λ 1 R λ 1 1 + R λ 0 1 + R λ 0 1 + R λ 1 + o λ 0 λ 1 R λ 0 R λ 1 (43) it follows that R = 1 + o λ 0 λ 1 R λ 0 R λ 1 . ( 44
)
This result is comparable to equation 30 and shows that the error induced by using the bulk thermodynamic constant and the fractionation factor instead of isotopes specific thermodynamic constants is negligible. This assumption can then be safely used.

Calculation: Developments for kinetic rate constant calculation

Next we turn to a comparable calculation for kinetically regulated reactions. The formalism for kinetic reactions is consistent with that developed and used by [START_REF] Lasaga | Chemical kinetics of water-rock interactions[END_REF], [START_REF] Steefel | Modeling coupled chemical and isotopic equilibration rates[END_REF], [START_REF] Depaolo | Surface kinetic model for isotopic and trace element fractionation during precipitation of calcite from aqueous solutions[END_REF], [START_REF] Druhan | Calcium isotope fractionation in groundwater: Molecular scale processes influencing field scale behavior[END_REF].

Starting with the formation of a species Ē with two isotopes comprising species E 0 (E 0 0 and E λ 0 )

E 0 + Nc-1 i=1 ν i E i - → Ē ( 45
)
we can defined the rate + V of both reactions involving the different isotopologues

+ V 0 = + k 0 a 0 0 N c-1 i=1 a νi i ( 46
)
+ V λ = + k λ a λ 0 N c-1 i=1 a νi i . (47)
The global rate of the bulk reaction is then defined by

+ V = + ka 0 • N c-1 i=1 a νi i = + V 0 + + V λ (48)
Similarly to the thermodynamic developments, we investigate here the ratio of the exact global rate over the approximated rate: R =

+ V + V = k 0 a 0 0 + k λ a λ 0 k a 0 0 + α λ a λ 0 = k 0 a 0 0 + α λ a λ 0 k a 0 0 + α λ a λ 0 = k 0 k . (49)
From equation 48 it follows that:

ka 0 = k 0 a 0 0 + α λ k 0 a λ 0 = k 0 a 0 0 + a λ 0 -λ a λ 0 ( 50
)
and thus

k 0 k = a 0 0 + a λ 0 a 0 0 + a λ 0 -λ a λ 0 = 1 1 - a λ 0 a 0 0 + a λ 0 . (51)
Using Taylor expansion:

k 0 k = 1 + a λ 0 a 0 + o λ a λ 0 a 0 2 = 1 + λ R λ 1 + R λ + o λ R λ 1 + R λ 2 ( 52
)
where R stands for the isotopic ratio. We note that in the case of a kinetic reaction forming or dissolving a mineral, one may wish to develop the expressions shown above as a 'forward' reaction and track overall mineral growth or dissolution as the net rate arising from the difference between this reaction and a 'backward' reaction proceeding in the opposite direction. The same development can be done for the backward reaction and thus these unidirectional pathways can be combined into a net rate and applied to the common class of fractionating reactions associated with mineral growth and solubilization.

Discussion

The above derivations extend the validity of equations with isotopic fractionation which can be found in certain reactive transport codes, without the need for the major isotope assumption in the case of thermodynamic equilibrium. These mathematical developments are not exhaustive, but they expand the capacity to handle isotopic fractionation modeling for a large variety of chemical species. These developments show that no a priori assumption is required for the relative isotopic abundances. Instead, the ratios may be defined relative to an arbitrary choice of reference isotope. The possible deviation from the exact solution is introduced through assumptions regarding the (small) magnitude of enrichment factors. In natural systems these distributions are usually narrow, and thus this assumption improves upon the prior approach in which the presence of a major isotope in a given element was required. When neglecting terms in ( R) 2 , the equations correspond to the ones obtained using the major isotope hypothesis and the classical approach.

The derivations for kinetically controlled reactions such as mineral precipitation and dissolution show that the error introduced when using k and α λ k instead of k 0 and k λ depends on R, which changes over time during the progress of the reaction. This issue has been considered in prior studies, largely in the context of organic contaminants [START_REF] Hunkeler | Modeling chlorine isotope trends during sequential transformation of chlorinated ethenes[END_REF], [START_REF] Van Breukelen | Modeling 3d-csia data: Carbon, chlorine, and hydrogen isotope fractionation during reductive dechlorination of tce to ethene[END_REF]. A solution was offered using a probability term which is a function of the abundance of the heavy and light isotopes through the course of the reaction, the number of atoms possibly isotopically substituted in the molecule and number of reactive positions among them (described as the 'isotopologues approach' in [START_REF] Hunkeler | Modeling chlorine isotope trends during sequential transformation of chlorinated ethenes[END_REF]). The reaction rates resulting from this approach are distinct from the bulk rate, such that the light isotope reacts slightly faster, and the heavy isotope slightly slower than what would be the bulk rate. The average rate weighted by the probability of the presence of a heavy isotope then reproduces the bulk rate.

Here, in absence of specific experimentally determined rate constants for the isotopologues, we offer a generalized demonstration that assuming the rate of one isotopologue is comparable to the bulk rate constants and deriving the other rate constants using the fractionation factor is a reasonable approximation when dealing with elements which present a clear major isotope (isotopic ratio R is then small). The error introduced by this approximation may become problematic when dealing with elements with isotopes of comparable relative abundance, and thus the use of a correction term as detailed in [START_REF] Hunkeler | Modeling chlorine isotope trends during sequential transformation of chlorinated ethenes[END_REF] and [START_REF] Van Breukelen | Modeling 3d-csia data: Carbon, chlorine, and hydrogen isotope fractionation during reductive dechlorination of tce to ethene[END_REF] could become warranted. We provide a framework to evaluate the magnitude of these errors in application to kinetic rate expressions, and extend consideration of these assumptions to applications beyond organic contaminants, such as mineral precipitation and dissolution. These developments are provided with the intent of raising awareness among the reactive transport community about the potential errors introduced when modeling fractionation of isotopes lacking a clear major abundance, which have thus far largely been limited to the specific context of aqueous chlorinated compounds.

As analytical capabilities increasingly enable determination of the specific location of substituted atoms within a molecular structure ( [START_REF] Corso | High-precision position-specific isotope analysis[END_REF], [START_REF] Hoffman | Position-specific carbon stable isotope ratios by proton nmr spectroscopy[END_REF], [START_REF] Remaud | Position-specific isotope analysis by isotopic nmr spectrometry: New insights on environmental pollution studies[END_REF]), it remains to be shown whether such an approach will facilitate treatment of the complexity of partitioning coefficients in natural systems.

Application case study

As an illustrative example of the utility of the previous sections, we provide an application of geochemical modeling including equilibrium isotopic fractionation. Two simulations were performed with the CHESS speciation code [START_REF] Van Der Lee | Module-oriented modeling of reactive transport with hytec[END_REF], [START_REF] Van Der Lee | Thermodynamic and mathematical concepts of chess[END_REF] and are inspired by the distribution of zinc isotopes (no major isotope) within Mississippi Valley Type (MVT) ore deposits. In this system, zinc isotopic signatures have been shown

to differentiate depending upon the history of the deposit. MVT ore formation is characterized by the combination of hydrothermal systems and major crustal tectonic events [START_REF] Leach | Mississippi valley-type lead-zinc deposits through geological time: implications from recent age-dating research[END_REF], [START_REF] Anderson | Ore deposit models-6. mississippi valley-type lead-zinc deposits[END_REF]. The precipitation of the ore is related to fluid migration and more specifically to dense basinal brines with temperatures between 75 and 200C, typically in carbonated platforms [START_REF] Leach | Mississippi valley-type lead-zinc deposits through geological time: implications from recent age-dating research[END_REF]. MVT ores gather a large variety of deposits and some aspects are still debated, including the origin of the ore fluids, the timescale of fluid migration and the processes of precipitation. Despite this uncertainty, prior studies have shown that a plausible condition for ore formation involves mixing between hot fluids containing elevated metal concentrations and cooler fluids with high reduced sulfur content [START_REF] Leach | Precipitation of lead-zinc ores in the mississippi valley-type deposit at trèves, cévennes region of southern france[END_REF]. This overly simplified illustrative case can be linked to the second part of section 3 about double substitution, namely 2 isotopes of a single Zn atom in Zn 2+ and ZnS compounds, with comparable abundances for the two Zn isotopes. The two simulations illustrate equilibrium and kinetic fractionation respectively.

Equilibrium fractionation. We provide a simplified simulation of a progressive mixing of two (sulfide bearing and zinc bearing) solutions leading to sphalerite precipitation with associated zinc isotopes fractionation.

The simulation consists of a sequential addition of a constant volume (0.002 L per sample, 100 samples) of solution 2 into solution 1 (table 1). The cumulative volume of the mixed solution is denoted as ζ. We illustrate here that isotopic enrichment at equilibrium can be simulated through the definition of appropriate thermodynamic constants. In this model, we focus on the isotopic composition of the mineral sphalerite.

To this purpose, sphalerite was considered as a solid solution. This approach allows us to define different endmembers of pure 64 Zn or 66 Zn sphalerite and all intermediate compositions. As the ratio 66 Zn

64 Zn is close to 0.5, using the major isotope hypothesis to assume that 64 K Sphalerite =K Sphalerite seems hazardous. But the developments detailed in this study showed that this assumption can indeed be performed safely.

Initial concentrations of 64 Zn 2+ and 66 Zn 2+ were calculated in order to produce a cumulative concentration of 1 × 10 -2 mol/L with the following parameters : a constant fractionation factor α eq =0.9997

(corresponding to the fractionation during sphalerite precipitation from ZnSO 4 ) was selected [START_REF] Veeramani | Zinc isotope fractionation as an indicator of geochemical attenuation processes[END_REF]. We refer here to an equilibrium fractionation factor in order to illustrate the previous developments, even if using an example of mineral precipitation. The same example is used later to illustrate the kinetic developments (with a TST rate law, which is usually used for kinetically-mediated mineral growth). Due to the lack of data, no dependence of the fractionation factor to the temperature is taken into account here. Initial δ 66 Zn was set at 0 (arbitrarily chosen, δ 66 Zn of terrestrial geological materials are within a range of -0.80 and 2.5 [START_REF] Ducher | Equilibrium zinc isotope fractionation in zn-bearing minerals from first-principles calculations[END_REF]) and the standard isotopic ratio 66 Zn 64 Zn was set to 0.56502 [START_REF] Maréchal | Precise analysis of copper and zinc isotopic compositions by plasma-source mass spectrometry[END_REF].

Zn 2+ + H 2 S ZnS + 2H + (53)
All thermodynamic constants K for sphalerite precipitation [START_REF] Lichtner | Continuum formulation of multicomponent-multiphase reactive transport[END_REF] were calculated using equations [START_REF] Bouchez | Modeling novel stable isotope ratios in the weathering zone[END_REF] to [START_REF] Cloquet | Variation in the isotopic composition of zinc in the natural environment and the use of zinc isotopes in biogeosciences: a review[END_REF] at appropriate temperatures (table 2) following the methodology described in this study.

During the course of the simulation, pH decreases from 5.87 to 5.29 (figure 1a) and temperature increases from 60 to 75 • C from pure solution 1 to a mixing ratio (figure 1b). Zn 2+ decreases of 2 × 10 -3 mol in stoichiometric agreement with sphalerite precipitation (figure 1c). Evolution of δ 66 Zn in both zinc and sphalerite shows that the implemented thermodynamic constants accurately reproduce the reported enrichment of -0.3 during precipitation of sphalerite (figure 1d). Kinetic fractionation. We also provide an example of fractionation modeling during kinetically controlled reaction. The same system was used with the difference that precipitation of sphalerite is kinetically constrained. The simulation was created without comparison to actual experimental results: its main purpose is to illustrate appropriate behavior isotopic fractionation during kinetically controlled reactions where a major isotope is lacking. The law and parameters determined by [START_REF] Acero | Sphalerite dissolution kinetics in acidic environment[END_REF] were used to derive the kinetic rate constant for sphalerite precipitation at 70 • C: +5.04848 × 10 -12 mol/m 2 /s in precipitation for 64 Zn. Then we calculated the rate constant for 66 Zn using the same fractionation factor than for the equilibrium precipitation case (0.9997), leading to a rate of +5.046965456 × 10 -12 mol/m 2 /s. The model simulates an initial solution oversaturated in zinc at 70 • C and advances through time-dependent sphalerite formation over a duration of 1 hour (figure 2). At the end of the simulation, the solution is still over-saturated with respect to sphalerite.

The simulation results show that sphalerite precipitates faster at the beginning, when the conditions are the farthest from equilibrium (figure 2a). The precipitation rate decreases as the conditions get closer to equilibrium. The δ 66 Zn values show that sphalerite isotopic composition reaches the initial composition of aqueous Zn 2+ (i.e 0 ) rapidly (as sphalerite precipitates rapidly in early time) and remains below the values of aqueous δ 66 Zn (figure 2b). In the aqueous phase, δ 66 Zn increases rapidly due to the fast precipitation rate, and this kinetic fractionation persists as the solution is still over-satured with respect to sphalerite. By the end of the simulation, the rate has become reduced due to convergence to equilibrium conditions. The evolution of δ 66 Zn in the aqueous phase gives information about the δ 66 Zn of the incremental precipitating Zn (with a -0.3 enrichment) but the δ 66 Zn in the solid phase is an accumulated averaged, such that the sphalerite isotopic signature evolves mainly at the beginning of the simulation.

It is worthwhile to recall that in these CHESS simulations, as in most speciation and RT codes, the two isotopes of Zn are entered as unique 'species' with individual constants. The kinetic reaction driving sphalerite formation couples these two 'species' through the solid solution, and thus one may nominally sum up the concentration of Zn 64 and Zn 66 at any point in the simulation to retrieve a 'bulk' Zn concentration. In this example, as with many applications of isotope-enabled RT, there are no isotope-specific rate constant available that would uniquely constrain these constants. Rather, it is necessary to use the bulk rate constant and infer the associated values for the remaining isotopes. As a result, it is not possible here to determine the error (and thus the difference in isotopic ratio) introduced by the use of the bulk rate and the fractionation factor rather than isotopes specific rate constant as detailed by [START_REF] Hunkeler | Modeling chlorine isotope trends during sequential transformation of chlorinated ethenes[END_REF], [START_REF] Van Breukelen | Modeling 3d-csia data: Carbon, chlorine, and hydrogen isotope fractionation during reductive dechlorination of tce to ethene[END_REF]. As a point of reference, we illustrate the mass-weighted rate of Zn precipitated when using unique 'species' Zn 64 and Zn 66 with constants k and αk relative to that using a single 'bulk' Zn species (no isotopes) with rate constant k. The sum of the initial mass of Zn 64 and Zn 66 is equal to the initial mass of the 'bulk' Zn simulation, and as described by [START_REF] Hunkeler | Modeling chlorine isotope trends during sequential transformation of chlorinated ethenes[END_REF], [START_REF] Van Breukelen | Modeling 3d-csia data: Carbon, chlorine, and hydrogen isotope fractionation during reductive dechlorination of tce to ethene[END_REF] there is a slight (not exceeding 5%) difference in absolute value. It is vital to appreciate that this slight difference arises because one isotope of Zn was associated with the 'bulk' rate constant k while the other was associated with a slightly different constant αk. The alternative is to assign unique values of rate constant to both isotopes, where neither are precisely equivalent to k, but in doing we both (1) lack uniquely determined values for these constants and (2) return to a dependence on the isotope ratio R (equation 52)

which evolves over the progression of the reaction.

Conclusions

The mathematical equations defining the implementation of equilibrium and kinetic isotopic fractionation in RT codes are detailed in this paper. The classical approach is based on the hypothesis of an existing major isotope for each element of interest leading to the assumption that the constant of the isotopologue containing the major isotope is equal to the "classical" constant. In this study, we considered generic cases regarding the abundance of isotopes and the potential errors introduced when no isotope can be defined as major. We provided mathematical developments based on small deviations around enrichment factors, which are legitimate in virtually all natural systems. The calculations lead to thermodynamic constants for isotopologues depending on known thermodynamic constants and fractionation factors. Interestingly, the constants are the same as when considering major isotopes although the results stem from small deviations of R instead of negligible concentration(s) of the minor isotope(s).

This study demonstrates that no hypothesis regarding the isotopes abundances is needed to assume that the thermodynamics of one isotopologue can be defined as equal to the general thermodynamics of the studied reaction, extending the scope of validity for this common assumption. The difference between the thermodynamic constants assigned to each isotopologue is the key for accurate modeling of isotopic fractionation.

However, we show that the same methodology has to be taken with caution when applied to kinetically driven reactions such as mineral precipitation and dissolution. As the rate constant depends on the isotopic ratio, which evolves with the progression of the reaction, an error is introduced when using the bulk rate constant and fractionation factor only when determining the isotopologue rate constants. This issue has already been highlighted in other studies investigating the question of modeling isotopic fractionation during organic compound degradation but has not been detailed in the case of minerals precipitation and dissolution reactions, which are of primary importance in a large number of reactive transport applications. The introduced error can be neglected when one isotopes of the element of interest can be considered as major, leading to small isotopic ratios. However, when modeling isotopic fractionation of elements presenting isotopes of comparable relative abundance, modelers should be aware of the introduced error depending on the methodology employed.

Figure 1 :

 1 Figure 1: Simulation results during mixing: a. pH, b. temperature evolution, c. zinc and sphalerite evolution, d. δ 66 Zn evolution in zinc and sphalerite

Figure 2 :

 2 Figure 2: Simulation results during mixing: a.zinc and sphalerite evolution, b. δ 66 Zn evolution in zinc and sphalerite, c. Global rates of the bulk reaction and the isotopes dependent reaction

Table 1 :

 1 Simulation parameters (balance means that it was adjusted to ensure solution electroneutrality)

		solution 1	solution 2
	volume (L)	1	0.2
	temperature (C)	60	150
	dolomite (mol/L)	1	0
	64 Zn 2+ (mol/L)	6.33081e-03	0
	66 Zn 2+ (mol/L)	3.66919e-03	0
	Cl -(mol/L)	balance	0
	H 2 S (mol/L)	0	1e-02
	Na + (mol/L)	0	balance
	CO 2 fugacity (atm)	1	60
	pH	equilibrium (5.87)	4.2

Table 2 :

 2 Thermodynamic constants of sphalerite precipitation

	Temperature (C) log K Sphalerite 64 Zn log K Sphalerite 66 Zn
	0	4.632700000	4.632569692
	25	4.452300000	4.452169692
	60	4.185100000	4.184969692
	100	3.876800000	3.876669692
	150	3.51200000	3.5111869692
	200	3.190400000	3.190269692
	250	2.928600000	2.928469692
	300	2.740400000	2.740269692

Acknowledgements

The authors sincerely thank the reviewer for the very relevant remark made about the mathematical developments. This comment allowed a significant improvement of the work presented in this study.

Funding

This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.