Maxime Julien
email: julien.m.aa@m.titech.ac.jp

Alexis Gilbert

Keita Yamada

Richard J Robins

Patrick Höhener

Naohiro Yoshida

Gérald S Remaud

Al. -Talanta

Expanded uncertainty associated with determination of isotope enrichment factors: comparison of two-points calculation and Rayleigh-plot

Keywords: Isotope Effects, isotope enrichment factor, expanded uncertainty, Rayleigh, twopoints

The Enrichment Factor () is a common way to express Isotope Effects (IEs) associated with a phenomenon. Many studies determine  using a Rayleigh-plot, which needs multiple data points. More recent articles describe an alternative method using the Rayleigh equation that allows the determination of  using only one experimental point, but this method is often subject to controversy. However, a calculation method using two points (one experimental point and one at t 0 ) should lead to the same results because the calculation is derived from the Rayleigh equation. But, it is frequently asked "what is the valid domain of use of this two-points calculation?" The primary aim of the present work is a systematic comparison of results obtained with these two methodologies and the determination of the conditions required for the valid calculation of . In order to evaluate the efficiency of the two approaches, the expanded uncertainty (U) associated with determining  has been calculated using experimental data from three published articles. The second objective of the present work is to describe how to determine the expanded uncertainty (U) associated with determining . Comparative methodologies using both Rayleigh-plot and two-points calculation are detailed and it is clearly demonstrated that calculation of  using a single data point can give the same result as a Rayleighplot provided one strict condition is respected: that the experimental value is measured at a small fraction of unreacted substrate (f < 30%). This study will help stable isotope users to present their results in a more rigorous expression:  ± U and therefore to define better the significance of an experimental results prior interpretation.

Capsule: Enrichment factor can be determined through two different methods and the calculation of associated expanded uncertainty allows checking its significance.

INTRODUCTION

Stable isotope analyses is widely used in different fields such as geochemistry [START_REF] Galimov | Isotope organic geochemistry[END_REF], biology [START_REF] Hayes | Fractionation of Carbon and Hydrogen Isotopes in Biosynthetic Processes[END_REF] or environmental sciences [START_REF] Elsner | Current challenges in compound-specific stable isotope analysis of environmental organic contaminants[END_REF]. The isotope ratio of light elements (C, H, O, N, S or Cl) is known to change through processes such as (bio)chemical reactions or physical processes due to preferential selection for the light or heavy isotope. This physical phenomenon, called the Isotope Effect (IE), can occur at equilibrium where heavy isotopes will be accumulated in one of the components of the system. Many processes such as liquid-vapor transformation can be associated with an Equilibrium Isotope Effect (EIE) [START_REF] Julien | Fractionation in position-specific isotope composition during vaporization of environmental pollutants measured with isotope ratio monitoring by 13 C nuclear magnetic resonance spectrometry[END_REF]. IE can also be associated with a (bio)chemical reaction [START_REF] Deniro | Mechanism of carbon isotope fractionation associated with lipid synthesis[END_REF] where the presence of heavy isotopes can influence the reaction rate constant.

Classically, this Kinetic Isotope Effect (KIE) is defined as the ratio of rate constants between two isotopes of a given element KIE = ( light k/ heavy k) and depends on the mechanism of the reaction/process and on environmental parameters such as temperature or pH [START_REF] Elsner | A New Concept Linking Observable Stable Isotope Fractionation to Transformation Pathways of Organic Pollutants[END_REF]. As a consequence, kinetic processes can be associated with an enrichment (KIE<1) or depletion (KIE>1) of heavy isotope of the reaction product(s) compared to the substrate(s). For a more detailed explanation, we recommended to read the guide "Practice and Principles of Isotopic Measurements in Organic Geochemistry" by J. M. Hayes where these concepts are explained in a very pedagogical way [START_REF] Hayes | Practice and principles of isotopic measurements in organic geochemistry[END_REF]. IE is not directly observable but it induces an alteration of isotopic abundance in the studied compounds: substrate(s) and/or product(s). This change in the isotope ratio, called isotopic fractionation, can be determined using isotope analysis and is often expressed as a fractionation factor () calculated from the equation described by Bigeleisen and Wolfsberg [START_REF] Bigeleisen | Theoretical and Experimental Aspects of Isotope Effects in Chemical Kinetics[END_REF]. In the present article, isotopic fractionation has been expressed as an enrichment factor () [START_REF] Mariotti | Experimental determination of nitrogen kinetic isotope fractionation: Some principles; illustration for the denitrification and nitrification processes[END_REF] instead of  (where  = ( -1)*1000). This notation has been chosen, because the sign of  directly informs on the IE associated with the studied process: when  is negative, the IE is considered as normal (light isotopes are preferentially used during the process) and when is positive, the IE is inverse (transformation is faster using heavy isotopes). Enrichment factors can easily be determined using a Rayleigh-plot where  is the slope of the trend line [START_REF] Elsner | Insight into Methyl tert-Butyl Ether (MTBE) Stable Isotope Fractionation from Abiotic Reference Experiments[END_REF][START_REF] Höhener | Rayleigh equation for evolution of stable isotope ratios in contaminant decay chains[END_REF].

Therefore, knowing the sign and amplitude of the IE allows conclusions to be drawn on the process and/or mechanism responsible for the production, degradation, or transformation of a given molecule.  (expressed in ‰) can be defined as follows:

 = ( heavy k/ light k -1)*1000 (1)
Note that in geo-and environmental sciences, in contrast to classical (bio)chemistry, the ratio of the heavy to light isotope (and not vice versa) is commonly used.

However, what is usually measured during a given process is not directly the relative rate of the reaction but rather the isotope ratio of the starting material and the product at two or more points during the progress of the process/reaction. The link between isotope ratios, extent of progress of the reaction f and  for a kinetic process is given by the Rayleigh equation [START_REF] Mariotti | Experimental determination of nitrogen kinetic isotope fractionation: Some principles; illustration for the denitrification and nitrification processes[END_REF][START_REF] Kendall | Isotope tracers in cachment hydrology[END_REF][START_REF] Aelion | Environmental Isotopes in Biodegradation and Bioremediation[END_REF][START_REF] Hunkeler | Quantification of Isotope Fractionation in Experiments with Deuterium-Labeled Substrate[END_REF]:

R ⋍ R 0 × 𝑓 Ɛ/1000 (2)
where, R 0 is the isotope ratio at t 0 , R the isotope ratio at t f and f is the unreacted fraction of substrate at t f . Equation 2 is a good approximation for all elements where the heavy isotope is much rarer than the light isotope. In cases of near equal abundance of isotopes, such as for B or when working with highly enriched compounds, equation 2 must be replaced by that proposed by Hunkeler et al. in 2002 [14].

The most common method for determining  is to use a Rayleigh-plot: the isotopic composition of the remaining reactant is measured at different fractions of unreacted substrate (f) and ln(R/R 0 ) is plotted as a function of ln(f) [START_REF] Elsner | Insight into Methyl tert-Butyl Ether (MTBE) Stable Isotope Fractionation from Abiotic Reference Experiments[END_REF][START_REF] Cretnik | Chlorine Isotope Effects from Isotope Ratio Mass Spectrometry Suggest Intramolecular C-Cl Bond Competition in Trichloroethene (TCE) Reductive Dehalogenation[END_REF]. The enrichment factor  can then be obtained from the slope of a linear regression fitted to all data points. The variability of  is then expressed from the standard deviation of the slope.

Most studies involving the determination of IEs have exploited this approach. In many situations, such as field sampling, measurement conditions do not allow isotope analysis at different times or different degrees of f to be determined [START_REF] Mcnevin | Measurement of (carbon) kinetic isotope effect by Rayleigh fractionation using membrane inlet mass spectrometry for CO 2 -consuming reactions[END_REF].

Nonetheless, an enrichment factor can still be obtained by using equation 2 above (Rayleigh equation). This calculation method only requires isotopic measurements on the starting substrate (at t = 0, f 0 = 1) and on one data point (at t = t 1 , f 1 < 1): this calculation method is called two-points calculation in this article. The question is then:

how to express the variability of the results?

The goal of the present article is to compare these two approaches, namely

Rayleigh-plot and two-points calculation, both in terms of trueness and precision [START_REF] Coplen | Guidelines and recommended terms for expression of stableisotope-ratio and gas-ratio measurement results[END_REF].

Note that these two  determination approaches are derived from the Rayleigh distillation equation, so they should lead to the same . We asked three fundamental questions: "what is the domain of use of the Rayleigh equation?", "how to express the uncertainty of the results from the two-points calculation mode?" and "what is the significance of the results based on the uncertainty determination in both  determination modes?". Three articles have been selected, with the authors' consent, to compare these two methods for the determination of the enrichment factor. Using data from these studies, we first compare  values obtained with these alternative approaches in order to evaluate under which conditions they give consistent results.

In a second part, we detail how to determine the limit of significance of the calculated . Many studies only use the standard deviation (SD) of the isotopic measurement to indicate the significant threshold of the determination of . But this method does not

M Julien et al. -Talanta 7 / 22
include other sources of uncertainty such as the measurement of f, signal size or linearity [START_REF] Sherwood Lollar | An Approach for Assessing Total Instrumental Uncertainty in Compound-Specific Carbon Isotope Analysis: Implications for Environmental Remediation Studies[END_REF]. Calculation of the expanded uncertainty (named U) associated with the determination of  is explained in the section below and a more detailed U calculation is described in the 'guide of expression of uncertainty in measurement' GUM [START_REF] Jcgm/Wg | Evaluation of measurement data-Guide to the expression of uncertainty in measurement. Evaluation of measurement data-Guide to the expression of uncertainty in measurement 100[END_REF].

Calculations of both  and its associated U will help in determining IEs in a more robust way and will provide a correct style to express the accuracy of the results.

MATERIALS AND METHODS

Data collection

Experimental data from three different articles have been collated, with the authors' Both carbon and chlorine isotope ratios were monitored throughout the biotransformation using irm-GC/MS and f was determined by gas chromatography.

Article C, Yamada et al. 2014 [START_REF] Yamada | Evaluation of commercially available reagents as a reference material for intramolecular carbon isotopic measurements of acetic acid[END_REF], deals with position-specific isotope analysis of acetic acid using on-line pyrolysis combined with GC-C-irm-MS (Py-irm-GC/MS) [START_REF] Corso | High-precision position-specific isotope analysis[END_REF].

The authors determined position-specific isotope effects (PSIEs) associated with the evaporation of acetic acid. For this experiment, 14 vials containing pure acetic acid were placed in a hermetically-sealed chamber containing sodium hydroxide (NaOH)

to remove evaporated acetic acid. Vials are weighted at t 0 and at t f in order to determine f and PSIA were performed on each sample at variable f.

Determination of the enrichment factor ()

The two methods for  determination are described in this section.

The most common method is the determination of  using a Rayleigh-plot, in which ln(f) is plotted on the x-axis and ln(R/R 0 ) on the y-axis [START_REF] Elsner | Insight into Methyl tert-Butyl Ether (MTBE) Stable Isotope Fractionation from Abiotic Reference Experiments[END_REF][START_REF] Höhener | Rayleigh equation for evolution of stable isotope ratios in contaminant decay chains[END_REF]. In this graph, the enrichment factor () corresponds to the slope of the trend line and the associated The second method is the determination of  with the two-points calculation and has been performed using the following relationship [START_REF] Jeannottat | Chlorine and Carbon Isotopes Fractionation during Volatilization and Diffusive Transport of Trichloroethene in the Unsaturated Zone[END_REF], in which the isotope ratios R (as in equation 2) are expressed as isotopic composition :

Ɛ = ln ( δ A X t𝑓 +1000 δ A X t0 +1000 ) ( 1000 ln 𝑓 ) (3)
where  A X t0 and  A X tf are respectively the isotopic composition of the starting reactant (initial state) and the remaining reactant at t f . This relationship can be approximated as described in the following equation [START_REF] Julien | Fractionation in position-specific isotope composition during vaporization of environmental pollutants measured with isotope ratio monitoring by 13 C nuclear magnetic resonance spectrometry[END_REF]:

Ɛ = Δδ A X ln 𝑓 (4)
where  A X stands for the isotopic fractionation between the reactant at t 0 and at time f ( A X =  A X tf - A X t0 ). These two equations give very similar results [START_REF] Mariotti | Experimental determination of nitrogen kinetic isotope fractionation: Some principles; illustration for the denitrification and nitrification processes[END_REF] and the simplified equation 4 is employed here for the determination of  with results from the selected articles.

Determination of expanded uncertainty (U)

In order to determine the significance threshold of the enrichment factor, the expanded uncertainty needs to be calculated. This article is the first which offers a fully-described expanded uncertainty calculation method within this context.

The determination of the expanded uncertainty associated with the enrichment factor determined using the Rayleigh-plot can be directly calculated using the function "LINEST" in Microsoft Excel TM , as described by Elsner et al. [START_REF] Elsner | Insight into Methyl tert-Butyl Ether (MTBE) Stable Isotope Fractionation from Abiotic Reference Experiments[END_REF]. This function calculates both the slope and the standard deviation of the trend line (STDV slope).

The expanded uncertainty can thus be calculated as follows:

U = k × STDV slope (5)
where k = 2 for a coverage factor of 95% (t value in Student table).

In the case where  is determined with a single data point, U can be calculated from the uncertainty estimation described in the 'guide of expression of uncertainty in measurement' GUM [START_REF] Jcgm/Wg | Evaluation of measurement data-Guide to the expression of uncertainty in measurement. Evaluation of measurement data-Guide to the expression of uncertainty in measurement 100[END_REF]. This evaluation criterion has been chosen because all parameters that can affect the final result are taken into account by U calculation, unlike in most articles in which authors only use the SD of the isotope analysis method as described by the instrument manufacturer. Parameters that need to be considered in the calculation of expanded uncertainty can be assembled in an Ishikawa diagram (also called fishbone diagram or cause-and-effect diagram, see Figure 1). Two contributions have to be taken into account, (i) the isotopic measurement of both the starting reactant (at t 0 ) and the remaining reactant (at t f ) along with the standard deviation of the measurements and (ii) the determination of f (and its associated variability) which can be performed by weighting, chromatography or spectroscopy (like NMR or Mass Spectrometry).

The expanded uncertainty can be calculated using the following equation:

U = k × Ɛ × √ ( u(δ A X) Δδ A X ) 2 + ( u(p) p ) 2 (6)
where k = 2 for a coverage factor of 95% and p is the quantity measured during f determination (a mass, an area). Then, u( A X) is the uncertainty associated with isotopic measurement at t 0 and t f obtained by the propagation uncertainty law:

u(δ A X) = √ SD² + SD² = √ 2 × SD² (7) M Julien et al. -Talanta 11 / 22
where SD is the standard deviation of the isotopic measurement. The analyte at the initial and the final states is measured using strictly the same protocol, i.e. there is no linearity problem. The next step is the calculation of u(p), namely the precision of the analytical method used to determine f, which is expressed as an uncertainty using a rectangular distribution (also called continuous uniform distribution) [START_REF] Jcgm/Wg | Evaluation of measurement data-Guide to the expression of uncertainty in measurement. Evaluation of measurement data-Guide to the expression of uncertainty in measurement 100[END_REF] and calculated using the following equation:

u(p) = √ ( v √3 ) 2 × k ( 8
)
where v is the standard deviation of the analytical method used to determine f (balance, spectrometric technique), p is the quantity measured (mass, area) and k = 2 for a coverage factor of 95%.

In the general equation for U calculation (equation 6), the

negligible but needs to be taken into account in order to ensure that every step of the uncertainty calculation is performed in a correct way. The precision of f determination needs to be taken into account during  determination, because both the measurement method and the sampling can be associated with a systematic error [START_REF] Buchner | Experimental Determination of Isotope Enrichment Factors -Bias from Mass Removal by Repetitive Sampling[END_REF].

Once the enrichment factor and its associated expanded uncertainty are calculated, the final result is expressed as  ± U. This expanded uncertainty value is slightly overestimated but covers all sources of variability. In these conditions, the observed isotope effect is considered significant when ││ > │U│.

RESULTS AND DISCUSSION

Comparison between the two approaches

M Julien et al. -Talanta 12 / 22
In this study,  and the associated expanded uncertainty are calculated for each data point of the selected articles and the values are compared with those obtained using Rayleigh-plots. In this section, the studied element is specified in enrichment factor expression in order to make the reading easier (i.e.  C in the case of 13 C enrichment factor).

Figure 2 shows the Rayleigh-plots obtained for the datasets used in this study. In most cases, the Rayleigh-plot trend-line fits very well with the experimental points, with root mean square values (R 2 ) higher than 0.98, which clearly demonstrates the linearity of the experimental data (Figure 2). The only exception is the trend line obtained from the monitoring of  13 C of the carboxyl group of acetic acid during evaporation (Figure 2c) for which R² is about zero. This value is not surprising given that the slope of the curve is 0 (implying that there is no significant isotope fractionation associated with evaporation on this carbon-atom position).

For each dataset,  was determined using the slope of the trend line. Values for  obtained from each data point using the two-points calculation (equation 3) are presented in Table 1. In most cases, two-points calculation method gives very similar results to those from Rayleigh-plots, but only for small values of f (which corresponds to the fraction of unreacted substrate at t f ). As an example, results from article A (Figure 2a) show a very high  H of -65.7‰ obtained from Rayleigh-plot and twopoints calculations give values comprised between -63.0 and -69.3‰ when f is lower than 20%. The same phenomenon is observed in the case of 13 C with an  C of -1.5‰ according to Rayleigh-plot and very different values obtained with equation 3 when f > 24%. Rayleigh-plots drawn using experimental data from article B (Figure 2b) show an  Cl value of -5.0‰ and an  C of -19.3‰. These data also demonstrate the trueness of  calculated using equation 3, with very similar values obtained when f < 34.8% in as significant when │U/│< 1). These data mean that results obtained with the twopoints calculation are significant on a large range of transformation yields, even if the expanded uncertainty associated with the enrichment factor increases with the amount of remaining substrate. This graph demonstrates that the results obtained are significant when 1f > 0.7 (f < 30%). The relationship between the transformation yield (1f) and the gap between  values obtained with the two methods ((- Rayleigh )/ Rayleigh ) clearly demonstrates the trueness of  values determined with the two-points calculation (see Figure 3c), considering that the Rayleigh-plot allows the calculation of a true . Indeed, values obtained using the two-points method have a deviation of 10% relative to the true value ( Rayleigh ) when f < 20% (1f > 0.8) and a deviation of 20% when f < 30% (1f > 0.7), which means that the two-points method allows the calculation of true values of  when the selected data point has been measured at high transformation yield.

These results demonstrate that the two-points calculation method is efficient when the amount of remaining substrate is below 30% (which corresponds to a reaction yield of 70%).

CONCLUSIONS

Considering these observations, the determination of  using a single experimental point can be justified in conditions wherein the yield of the studied transformation is above 70% (f < 30%). Within this range, the  values calculated are similar to those obtained from Rayleigh-plots, and so is the associated uncertainty. Furthermore, expanded uncertainty calculations are shown as an excellent and rigorous way to check if the determined values of  are significant and, in particular, are a good tool to assess the validity of the two-points calculation.

EXCEL spreadsheet for calculation of enrichment factor and uncertainty by the twopoint method.

Tables Table 1: Values of  calculated using each experimental point from (a) Gray et al.

M

  Julien et al. -Talanta 5 / 22

M

  Julien et al. -Talanta 9 / 22 root mean square (R²) gives a first indication of the quality of the linearity. Note that the regression is not forced through the origin following a recommendation by Scott et al. in 2004 [24].

M

  Julien et al. -Talanta 10 / 22

22 Figure captions Figure 1 :

 22captions1 Figure captions

Figure 2 :Figure 3 :

 23 Figure 2: Rayleigh-plots drawn from experimental data of (a) Gray et al. 2002, (b) Cretnik et al. 2014 and (c) Yamada et al. 2014.

  part of article A fromGray et al. 2002 [20] deals with the determination of both deuterium and carbon IEs associated with the aerobic biodegradation of methyl tertbutyl ether (MTBE) by bacteria strains discovered in a military base located in California. Bioremediation experiments have been performed in vitro using sediments and groundwater extracted from the contaminated site. Stable isotope analysis was performed by Gas Chromatography -Combustionisotope ratio monitoring by Mass Spectrometry (irm-GC/MS)[START_REF] Brenna | High-precision continuous-flow isotope ratio mass spectrometry[END_REF] and the amount of unreacted substrate was measured by gas chromatography for each experimental point.

	[22]
	These articles have been chosen in order to explore diverse isotope fractionation
	conditions (evaporation, biodegradation) and also to give examples with different
	elements ( 2 H, 13 C, and 37 Cl). Moreover, this selection includes an article using
	position-specific isotope analysis (PSIA) in order to confirm that these calculation
	methods are applicable in such a case. In the present work, the articles are called

consent, to compare  values obtained with Rayleigh-plot and the two-points calculation and illustrate the advantage of performing the expanded uncertainty calculation. These are: Article A: Gray J. R. et al., 2002, Environ. Sci. Technol., 36, 1931-1938 [20] Article B: Cretnik S. et al., 2014, Molecules, 19, 6450-6473 [15] Article C: Yamada K. et al., 2014, Rapid Commun. Mass Spectrom., 28, 1821-1828 article A, B or C to make the text easier to read. Note that in all three articles a Rayleigh-plot was used to calculate the isotope fractionation factor. M Julien et al. -Talanta 8 / 22 A In article B, Cretnik et al. 2014 [15] report a study of the in vitro anaerobic biodegradation of tetrachloroethene (PCE) by Desulfitobacterium sp. strain Viet1.

Acknowledgment

M. Julien thanks the "Région Pays de la Loire" for funding his postdoctoral bursary through the project PLAISIR. A. Gilbert thanks the Grant-in-Aid for Young Scientist (B) (15K17774) and the Grant-in-Aid for Scientific Research (S) (23224013), MEXT, Japan, for financial support. The authors thank Prof. Barbara Sherwood-Lollar and Dr. Martin Elsner for agreeing to share their experimental data.

the case of 37 Cl and f < 31.1% for 13 C. Figure 2c demonstrates that the two-points calculation method is also efficient with PSIA data. Rayleigh-plot results from article C indicate that acetic acid evaporation is associated with a  C-Bulk of -1.0‰ mostly located on the methyl group ( C-Methyl = -2.0‰) which is confirmed by the comparable values obtained using the two-points calculation when f < 39.0%. Calculated  values also confirm the absence of significant isotope effect located on the carboxyl function.

According to data from Figure 2 and Table 1, the two methods seem to be exploitable to determine  from varied experimental data and elements. The limit of the twopoints calculation is the value of f, but this limit is not clearly established.

Determination of expanded uncertainty should help to fix the limit of the domain of use of the Rayleigh equation.

Expanded uncertainty

Expanded uncertainties associated with  determination using both Rayleigh-plot and two-points calculation are detailed in Table 1. When  is determined through Rayleigh-plot using experimental data from Article A, the expanded uncertainty associated with 13 C IE determination is very small (U = 0.1‰) compared to . That result suggests that  is significant. The U value is higher in the case of 2 H (3.7‰) but this result is in accordance with the higher  H of -65.7‰ detected through deuterium isotope analysis. These two examples demonstrate that an IE can be considered as significant even if its associated uncertainty is high; when ││ > │U│ (or │U/│< 1), the result is significant. Carbon-13 bulk enrichment factors obtained in the other articles are also significant with 1). The same conclusion can be made with results from 37 Cl experimental data of article B with  Cl = -5.0 ± 0.2‰. Enrichment factors calculated using the two-points approach are similar to those from Rayleigh-plots, as previously described by Jeannottat and Hunkeler [START_REF] Jeannottat | Chlorine and Carbon Isotopes Fractionation during Volatilization and Diffusive Transport of Trichloroethene in the Unsaturated Zone[END_REF], but U obtained with the two methods are very different and U seems to be increased with f (Table 1), which is coherent with the observations made above. In each study, when the process is not well advanced (f > 30%), the equation gives different  results than those obtained with the Rayleigh-plot. When less than 70% of substrate has reacted, the calculated  is not significant, as demonstrated in Figure 3b with │U/ Rayleigh │> 1.

For instance,  C-Bulk is -1.0 ± 0.1‰ (article C) according to the Rayleigh-plot and the same value is obtained from the two-points calculation with │U/ Rayleigh │ comprised between 0.3 and 0.7‰ provided f < 30% (Figure 3b). From all the data studied here, data points calculated at f > 30% give different values from those determined using Rayleigh-plot and this difference is confirmed by high expanded uncertainty values. The determination of │U/│ is an excellent tool to ensure that the uncertainties associated with the determination of IE do not exceed the measured  [START_REF] Singleton | High-Precision Simultaneous Determination of Multiple Small Kinetic Isotope Effects at Natural Abundance[END_REF].

Significance of the results

The results obtained from each article were used to draw three graphs (Figure 3) in order to evaluate the trueness of the  determined with the two-points calculation and the minimum of substrate transformation that must be used to ensure that results are significant.

The variation of the expanded uncertainty depending on the amount of remaining substrate is described in Figure 3a. This result confirms the need to set a limit of f value in order to use the two-points calculation in a correct way. The graph presented in Figure 3b shows the correlation between the reaction yield (1f) and the degree of significance of  calculated using the two-points method (the result is considered Enrichment factors can therefore be directly determined from experimental data using two different methods. The classical Rayleigh-plot approach is very efficient and the resulting  associated expanded uncertainty can be calculated using the function "LINEST" of Microsoft Excel TM . Equations that allow calculating enrichment factor with only one experimental point can be very useful when it is difficult to collect several samples during a transformation. When using this method, it is recommended to work with transformations of more than 70%.

For both methods, determining the expanded uncertainty associated with a result is highly recommended in order to draw any conclusion about the studied phenomenon. This calculation requires the determination of precision for both the isotopic measurement and the yield. In these conditions, the final result will take into account all uncertainty sources associated with both the experiments and the measurements.

The final IE result can be expressed as ± U and considered as significant when ││> │U│ (or │U/│< 1). In these conditions, the two-points calculation is an excellent way to determine the IE associated with a transformation and the use of the expanded uncertainty reinforces the validity of conclusions that can be made about the studied mechanism.

Such calculations can be complex to set up; so this article provides a spreadsheet in which both  and U can be automatically calculated after appropriate input: the initial isotopic composition (t = 0), the isotopic composition at t = f, the isotopic measurement precision (SD), the transformation yield (1 -f), the quantification measurement (mass or area) and the precision associated with f determination.

Supporting information available


Maxime Julien

Patrick Höhener

Richard J Robins

Julien Parinet

Gérald S Remaud
email: gerald.remaud@univ-nantes.fr

Position-Specific 13 C Fractionation during Liquid-Vapor Transition is Correlated to the Strength of Intermolecular Interaction in the Liquid Phase

Keywords: Position-specific isotope fractionation, isotope ratio monitoring 13 C NMR spectrometry, Volatilization, Distillation, Craig-Gordon model

 C fractionation observed during distillation was investigated. A range of molecules showing different intermolecular interactions in terms of mode and intensity were incorporated in the study. Although it had previously been suggested that during evaporation the diffusive 13 C isotope effect in the thin liquid layer interfaced with vapor is not position-specific, herein we show that this is not the case. In particular, the position-specific effect was demonstrated for a series of alcohols. Our hypothesis is that intermolecular interactions in the liquid phase are the source of position-specific 13 C fractionation observed on the molecule. A clear trend is observed between the 13 C isotope effect of the carbon bearing the heteroatom of chemical function and the relative permittivity, the solvent hydrogen-bond acidity and the solvent hydrogen-bond basicity, while only a weak trend was observed when using the 13 C content of the whole molecule.

Furthermore two families of products appeared when using the hydrogen-bond acidity parameter for the correlation by distinguishing H-acceptors and H-donors molecules to those H-acceptors only. This strongly reinforces the hypothesis of an important role of the 13 C positioned close to the interaction center.

Introduction

Stable isotopes were widely used in the past to study the nature of phase transfer processes like liquid-vapor transfer (References Jancso and van Hook, book of Wolfberg et al, etc.). [START_REF] Reichardt | Solvatochromism, thermochromism, piezochromism, halochromism, and chiro-solvatochromism of pyridinium N-phenoxide betaine dyes[END_REF] Past studies used fully substituted isotopic molecules for which large isotope effects are observed which could be better measured with past analytical techniques. Progress in analytical techniques has led to the development of Compound-Specific Isotope Analysis (CSIA) capable of monitoring stable isotope ratios at natural abundance ratios. With the help of CSIA, the mechanisms of phase transfer processes were furthermore investigated and elucidated (Hunkeler and Jeannottat, Kuder, Huang, etc.). However, it has been shown that CSIA based on isotope ratio monitoring by Mass Spectrometry (irm-MS) is not always sufficient to describe the [START_REF] Meyer | Interactions with Aromatic Rings in Chemical and Biological Recognition[END_REF] C fractionation during the whole evaporation phenomenon. Recently, we have reported the study of several liquid-vapor transition processes using isotope ratio monitoring by [START_REF] Meyer | Interactions with Aromatic Rings in Chemical and Biological Recognition[END_REF] C Nuclear Magnetic Resonance spectrometry (irm- [START_REF] Meyer | Interactions with Aromatic Rings in Chemical and Biological Recognition[END_REF] C NMR) in order to simulate position-specific fractionation during evaporation. [START_REF] Julien | Fractionation in position-specific isotope composition during vaporization of environmental pollutants measured with isotope ratio monitoring by 13 C nuclear magnetic resonance spectrometry[END_REF] We have shown that Position-Specific Isotope Analysis (PSIA) by irm- [START_REF] Meyer | Interactions with Aromatic Rings in Chemical and Biological Recognition[END_REF] C NMR can provide new insights into physical processes, including in circumstances in which CSIA may not detect any effect on the whole molecule. [START_REF] Julien | Insights into Mechanistic Models for Evaporation of Organic Liquids in the Environment Obtained by Position-Specific Carbon Isotope Analysis[END_REF] Furthermore, bulk [START_REF] Meyer | Interactions with Aromatic Rings in Chemical and Biological Recognition[END_REF] C composition ( [START_REF] Meyer | Interactions with Aromatic Rings in Chemical and Biological Recognition[END_REF] C bulk ) analysis (CSIA), obtained by irm-MS, only indicates which isotopologue reacts preferentially in the studied process, whereas PSIA shows which isotopomer is the most sensitive (see Figure 1 for a further definition).

Isotope fractionation is conveniently described by the isotopic fractionation factor  or the isotope effect  and is associated with discrimination between the behavior of isotopomers. [START_REF] Jeannottat | Chlorine and Carbon Isotopes Fractionation during Volatilization and Diffusive Transport of Trichloroethene in the Unsaturated Zone[END_REF] A previous study showed that the Craig-Gordon isotope model, originally derived for water evaporation, [START_REF] Craig | Deuterium and oxygen 18 variations in the ocean and the marine atmosphere[END_REF] is also valid for organic liquids as demonstrated by Kuder and collaborators, [START_REF] Kuder | Effects of volatilization on carbon and hydrogen isotope ratios of MTBE[END_REF] especially for the description of the volatilization of nonpolar organic liquids with air-side limitation of the volatilization rate. [START_REF] Julien | Insights into Mechanistic Models for Evaporation of Organic Liquids in the Environment Obtained by Position-Specific Carbon Isotope Analysis[END_REF] In the Craig-Gordon model in its original form, it is assumed that the diffusive effect in the stagnant liquid  diff-liq (thin liquid layer interfaced with vapor, see ref 4) is not position-specific, i.e. is independent of the position of the [START_REF] Meyer | Interactions with Aromatic Rings in Chemical and Biological Recognition[END_REF] C in the molecule.

However, among the compounds studied, ethanol and propan-1-ol both showed a positionspecific behavior. Working with the hypothesis that intermolecular interactions in the liquid phase were the source of position-specific [START_REF] Meyer | Interactions with Aromatic Rings in Chemical and Biological Recognition[END_REF] C fractionation observed on the molecule collected from the vapor phase, we have therefore generated a modified Craig-Gordon equation to take these position-specific interactions into account. [START_REF] Julien | Insights into Mechanistic Models for Evaporation of Organic Liquids in the Environment Obtained by Position-Specific Carbon Isotope Analysis[END_REF] The most widespread non-covalent molecular interaction is hydrogen bonding. The strength of the hydrogen bond depends on the properties of both the hydrogen-bond acceptor (HBA) and the hydrogen-bond donor (HBD). [START_REF] Laurence | The pKBHX Database: Toward a Better Understanding of Hydrogen-Bond Basicity for Medicinal Chemists[END_REF] Recently, a specific methodology has been proposed by which both solvent hydrogen-bond basicity and solvent hydrogen-bond acidity can be quantified, leading to the parameters  1 [START_REF] Laurence | Solvatomagnetic Comparison Method: A Proper Quantification of Solvent Hydrogen-Bond Basicity[END_REF] and ET(30), [START_REF] Cerón-Carrasco | Solvent polarity scales: determination of new ET(30) values for 84 organic solvents[END_REF] respectively.  1 is linked to the original scale of Kamlet and Taft 11 but is specific for the solvent scales. Thus, the scale of solvent hydrogen-bond basicity was established by the method that compares the 19 F NMR chemical shift of 4fluorophenol and 4-fluoroanisole in hydrogen-bond acceptor solvent. In addition to the empirical parameter describing the hydrogen-bond acidity  1 (from Kamlet-Taft), the quantitative solvent polarity parameter ET(30) was introduced. This is based on the negative solvachromism of 2,6diphenyl-4-(2,4,6-triphenylpyridinium-1-yl)phenolate (known as betaine 30). [START_REF] Reichardt | Solvatochromism, thermochromism, piezochromism, halochromism, and chiro-solvatochromism of pyridinium N-phenoxide betaine dyes[END_REF] Weaker interactions are encompassed under van der Waals, π-bonding [START_REF] Meyer | Interactions with Aromatic Rings in Chemical and Biological Recognition[END_REF] or halogen-bond [START_REF] Politzer | An overview of halogen bonding[END_REF] interactions that explain intermolecular associations between apolar molecules, molecules possessing an aromatic-ring, or halogen-containing compounds. The unit accounting for these interactions is the relative permittivity, [START_REF] Braslavsky | Glossary of terms used in photochemistry, 3rd edition (IUPAC Recommendations[END_REF] previously known as the dielectric constant diE .

11

The goal of the present work was to investigate the relationship between the strength of the intermolecular interaction in liquid and the position-specific [START_REF] Meyer | Interactions with Aromatic Rings in Chemical and Biological Recognition[END_REF] C fractionation observed during vaporization. We asked three fundamental questions: "Does the presence of [START_REF] Meyer | Interactions with Aromatic Rings in Chemical and Biological Recognition[END_REF] C in the molecule increase (or decrease) the interaction?", "Is this interaction stabilized in a specific [START_REF] Meyer | Interactions with Aromatic Rings in Chemical and Biological Recognition[END_REF] C isotopomer?", and "What is the influence of the proximity of [START_REF] Meyer | Interactions with Aromatic Rings in Chemical and Biological Recognition[END_REF] C to the interaction site?" To answer these questions conditions are required in which the isotope fractionation measured should reflect the transition liquid-vapor only with the influence on the liquid phase, but with no further isotope effect in the vapor phase. These conditions are essentially met in distillation in which no diffusive boundary layer is occurring when the compounds are boiling. Furthermore, this process is quick, simple, and known to induce significant isotope effects in 2 H and [START_REF] Meyer | Interactions with Aromatic Rings in Chemical and Biological Recognition[END_REF] C. [START_REF] Zhang | Determination of Multiple Equilibrium Isotopic Fractionation Factors at Natural Abundance in Liquid-Vapor Transitions of Organic Molecules[END_REF] Hence, this allows insight into isotope effects occurring on the liquid side. To express the isotope effect during the liquid-vapor transition we used the nomenclature already proposed for such transformation. 16a In the case of a small transformation of the substrate, such as the collection of a very small distillate volume (≤ 3%), any associated isotope effects can readily be detected in the distillate, but impact negligibly on the isotopomer composition of the remaining liquid. The corresponding fractionation factor  was defined in the classical literature as separation factor  (Jancso and Van Hook, 1974):

 = (N'/N) vap /(N'/N) liq (1)
where N' and N are the mole fractions of light and heavy isotopomers respectively in the system at equilibrium. The specific-isotope fractionation of site i (the isotopomer with [START_REF] Meyer | Interactions with Aromatic Rings in Chemical and Biological Recognition[END_REF] C at position isotopomer i) is then defined as the experimental fractionation factor  exp :

Distillation experiments

For each experiment, 100 mL of pure compound was introduced into a 1 L round bottom flask and distillation was carried out using a Cadiot distillation column equipped with a Teflon spinning band (see Murray and reference therein for its characteristics [START_REF] Murray | A modified spinning band column for low pressure fractionation[END_REF] ). This distillation equipment allows the reflux to be rigorously controlled. Once reflux conditions were reached, the distillate was collected up to approximately 3% maximum of the mass of the starting substrate.

The distillate was then submitted to [START_REF] Meyer | Interactions with Aromatic Rings in Chemical and Biological Recognition[END_REF] C-irm-EA/MS and irm- [START_REF] Meyer | Interactions with Aromatic Rings in Chemical and Biological Recognition[END_REF] C NMR.

In the present study, distillation has been chosen in order to study liquid-vapor equilibrium effects, because this method is quick and easy to perform which is perfect to study a large panel of chemical compounds in very controlled conditions. The inconvenient of this protocol is that the resulting isotope effects are magnified, especially using a spinning band distillation column which displays many theoretical plates and the capability of maintaining the system at equilibrium. Nevertheless, these overexpressed isotopic fractionations have the same direction (normal or inverse) as those observed in static phase equilibrium (Refs. Jeannottat+Hunkeler) or vapor pressure measurements (Jancso+Van Hook, 1974).

In the case of isotopic studies, the concept of "pure compound" needs to be further defined. The term "pure" means that there is only one chemical species present. However, it is in fact a mixture of all potential isotopologues and isotopomers of the considered compound. As an example, a sample of pure ethanol contains four types of molecules isotopically different in 13 C (see Figure 1). to the international reference (Vienna-Pee Dee Belemnite, V-PDB) using the relationship:

δ C VPDB 13 (‰) = ( R sample R standard
-1) × 1000. [START_REF] Muccio | Isotope ratio mass spectrometry[END_REF] The instrument was calibrated for δ [START_REF] Meyer | Interactions with Aromatic Rings in Chemical and Biological Recognition[END_REF] C using the international reference materials NBS-22 (δ [START_REF] Meyer | Interactions with Aromatic Rings in Chemical and Biological Recognition[END_REF] C PDB = -30.03 ‰), SUCROSE-C6 (δ [START_REF] Meyer | Interactions with Aromatic Rings in Chemical and Biological Recognition[END_REF] C PDB = -10.80 ‰), and IAEA-CH-7 PEF-1 (δ [START_REF] Meyer | Interactions with Aromatic Rings in Chemical and Biological Recognition[END_REF] C PDB = -32.15 ‰) (IAEA, Vienna, Austria) and instrumental deviation followed via a laboratory working standard of glutamic acid.

Isotope ratio monitoring [START_REF] Meyer | Interactions with Aromatic Rings in Chemical and Biological Recognition[END_REF] C NMR spectrometry (irm- [START_REF] Meyer | Interactions with Aromatic Rings in Chemical and Biological Recognition[END_REF] C NMR)

9 / 29
The sample preparation consisted of the successive addition in a 4 mL vial of the compound, the lock substance and the relaxation agent Cr(Acac) 3 (which acts to decrease the longitudinal relaxation time T1). [START_REF] Caytan | Precise and accurate quantitative 13 C NMR with reduced experimental time[END_REF] The respective amount of each was adapted according to (i) the T 1 values of the analyte, (ii) the reciprocal solubility with the deuterated solvents and/or the relaxation agent and (iii) the [START_REF] Meyer | Interactions with Aromatic Rings in Chemical and Biological Recognition[END_REF] C NMR spectrum: no peak overlapping. The exact preparation procedures were described previously. [START_REF] Julien | Fractionation in position-specific isotope composition during vaporization of environmental pollutants measured with isotope ratio monitoring by 13 C nuclear magnetic resonance spectrometry[END_REF] Quantitative 13 C NMR spectra were recorded using an AVANCE I 400 spectrometer (Bruker Biospin, Wissembourg, France), fitted with a 5 mm i.

Calculation of Ɛ bulk and Ɛ i

The isotope effect IE (‰) was obtained from the equation:

 = ( = -1) x 1000 (3)
where the isotope isotope effect (IE) is calculated as explained in detail in the work of Julien et al. [START_REF] Julien | Insights into Mechanistic Models for Evaporation of Organic Liquids in the Environment Obtained by Position-Specific Carbon Isotope Analysis[END_REF] where IE is equivalent to .

Here < is a normal isotope effect is observed while > indicates an inverse isotope effect.

Results and discussion

Position-specific [START_REF] Meyer | Interactions with Aromatic Rings in Chemical and Biological Recognition[END_REF] C analysis by NMR

10 / 29
PSIA may be performed by fragmenting the molecule in such a way that, after subsequent analysis of the resulting fragments by irm-MS, [START_REF] Meinschein | Intramolecular isotopic order in biologically produced acetic acid[END_REF] information on the isotopomer composition is retrieved. This fragmentation can conveniently be achieved on-line by pyrolysis prior to mass spectrometry analysis. [START_REF] Corso | High-precision position-specific isotope analysis[END_REF] However, this approach is not easily applicable to molecules containing more than 4 carbon atoms. In contrast, quantitative NMR gives access to all spectrally-resolved sites in the compound, i.e. relative titration of the [START_REF] Meyer | Interactions with Aromatic Rings in Chemical and Biological Recognition[END_REF] C isotopomers constituting the molecule when irm- [START_REF] Meyer | Interactions with Aromatic Rings in Chemical and Biological Recognition[END_REF] C NMR conditions are used. [START_REF] Bayle | Conditions to obtain precise and true measurements of the intramolecular 13 C distribution in organic molecules by isotopic 13 C nuclear magnetic resonance spectrometry[END_REF] It offers the advantage of directly analyzing the target molecule, hence avoiding the risk of technique-associated fractionation. The precision in irm-NMR depends primarily on the signal-to-noise ratio (SNR) and it is described in the present work by a standard deviation of 0.3‰ on the  13 C scale (‰). We have shown in a previous study that, under these conditions, the uncertainty on IE i is ±0.2‰ (95% confidence level) [START_REF] Jcgm/Wg | Evaluation of measurement data-Guide to the expression of uncertainty in measurement[END_REF] for the evaporation processes and the molecules used. [START_REF] Julien | Fractionation in position-specific isotope composition during vaporization of environmental pollutants measured with isotope ratio monitoring by 13 C nuclear magnetic resonance spectrometry[END_REF] This value is the threshold of significance that will be used for the interpretation of the results presented.

Bulk [START_REF] Meyer | Interactions with Aromatic Rings in Chemical and Biological Recognition[END_REF] C and position-specific fractionation upon distillation

Zhang et al. showed that distillation can be a useful tool to determine isotopic fractionation factors at equilibrium during liquid-vapor transitions. 16a For a given compound heated at reflux, an equilibrium between the liquid and the vapor is established. When a very small amount of distillate is collected (≤ 3%) it can be assumed that the isotope profile of the distillate represents the isotopic composition of the vapor at equilibrium with the liquid. These authors measured 2 H the position-specific isotope effect (by irm-2 H NMR) as well as isotope bulk effect (by irm-MS) on heteroatoms ( 13 C, [START_REF] Braslavsky | Glossary of terms used in photochemistry, 3rd edition (IUPAC Recommendations[END_REF] N and 18 O) in a series of organic molecules containing typical functional groups. They observed that the isotopic fractionation enrichments varied with (i) the element, (ii) the type of chemical function and (iii) the 2 H isotopomers ( 2 H position in the molecule). We previously confirmed that distillation under similar conditions also led to position-specific [START_REF] Meyer | Interactions with Aromatic Rings in Chemical and Biological Recognition[END_REF] C fractionation in ethanol: the distillate is enriched in [START_REF] Meyer | Interactions with Aromatic Rings in Chemical and Biological Recognition[END_REF] C mainly on the CH 2 moiety. 3, 16b In the present work we have extended this protocol to other organic compounds that have several specific chemical characteristics in terms of substituent, heteroatoms and chemical functions. The isotopomer profiles expressed as  13 C bulk ,  13 C i and  13 C (difference between the  13 C in distillate and  13 C in starting material) are summarized in Table 2. Two key observations can be made. First, in the majority of cases the distillate is enriched in [START_REF] Meyer | Interactions with Aromatic Rings in Chemical and Biological Recognition[END_REF] C. Second, a large variability of  [START_REF] Meyer | Interactions with Aromatic Rings in Chemical and Biological Recognition[END_REF] C i within and between molecules is seen. In a previous study (ref 4), [START_REF] Meyer | Interactions with Aromatic Rings in Chemical and Biological Recognition[END_REF] C isotopic fractionation associated with distillation has been determined analyzing the remaining substrate after distillation of 93 to 98% of pure compound. Under these conditions, it can be assumed that the detected isotopic fractionation corresponds to the isotope effect associated with the liquidvapor transformation without contribution of the diffusion. In the present case, the distillate has been analyzed in order to study the isotopic fractionation associated with distillation. The low distillation yields (around 3%) along with the presence of multiple theoretical plates both induce the detection of larger isotope effects than the analysis of the remaining substrate. However, this overestimation of the isotopic fractionation associated with distillation is not misleading for the interpretation because the intramolecular distribution of the fractionation is conserved. As an example, the analysis of the distillation tail of TCE showed an isotope effect of +0.8‰ mostly located on the carbon bearing two chlorine atoms (see ref 4) and this isotope effects is about +5.1‰ in the present study but still mostly located on the same carbon position (see Table 3).

The distillation protocol employed in this study doesn't allow obtaining true values of liquidvapor equilibrium isotope effect but it is still available to study its mechanism. Table 1 shows the numbering of the carbon positions, according to the shielding order of the peak in the 13 C NMR spectrum. This shielding (or deshielding) is associated with the electron density surrounding the given carbon atom and depends to some degree on the polarity of the neighboring heteroatoms (or chemical functions). The isotopomer 'C-1' which contains the [START_REF] Meyer | Interactions with Aromatic Rings in Chemical and Biological Recognition[END_REF] C directly linked to this heteroatom is clearly identified in alcohols. It is, however, less straightforward for other chemical species, such as ester, ketone, nitrile, etc. (see Table 1 for the attribution of C-1 in each molecule). In general, the [START_REF] Meyer | Interactions with Aromatic Rings in Chemical and Biological Recognition[END_REF] C ratio at the C-1 is seen to be the most effected by the liquid-vapor transition, i.e. it is the isotopomer containing [START_REF] Meyer | Interactions with Aromatic Rings in Chemical and Biological Recognition[END_REF] C at that position which distils first. Exceptions are notable for acetonitrile (acni), for which  13 C C-1 is depleted, and for n-heptane (hept) and toluene (tol) for which there is a homogeneous enrichment effect for each carbon. In the series of alcohols, a further effect is observed: the preference of the [START_REF] Meyer | Interactions with Aromatic Rings in Chemical and Biological Recognition[END_REF] C isotopomer for the vapor phase is diminished, stabilized or even increased according to the distance from the OH function along the carbon chain. Thus, in butan-1-ol (al-5) and pentan-1-ol (al-8)  [START_REF] Meyer | Interactions with Aromatic Rings in Chemical and Biological Recognition[END_REF] C upon distillation is large for the carbon-bearing OH (C-1): +5‰ and +4.1‰, lower for the inner carbon(s) of the chain:

+3.0/+2.6‰ and +2.7/+2.0‰ and slightly larger for the terminal carbon positions: +2.9 and +2.7‰ respectively (Table 2). This is clear evidence that the presence of [START_REF] Meyer | Interactions with Aromatic Rings in Chemical and Biological Recognition[END_REF] C close to the OH function influences strongly the facility for the molecule to be transformed into vapor.

Interestingly, a remote effect is associated with the terminal carbon.

Position-specific 13 C fractionation versus intermolecular interactions

The above discussion of the alcohol series addresses the issue of interactions within the liquid. As a first deduction, it appears that the proximity of [START_REF] Meyer | Interactions with Aromatic Rings in Chemical and Biological Recognition[END_REF] C to the OH function favors the distillation of the corresponding isotopomer. Therefore the question is: does  C-1 (the amplitude of the [START_REF] Meyer | Interactions with Aromatic Rings in Chemical and Biological Recognition[END_REF] C fractionation or the discrimination between the isotopomer C-1 and the other isotopomers)

represent the strength of the interaction between alcohol molecules in the liquid state? Among several parameters used for describing the potential level of intermolecular associations, the relative permittivity (diE) is the most commonly used. Nevertheless, there are some exceptions which deserve comment. In all cases, toluene (tol) and even more so n-heptane (hept) do not follow these trends: this is readily explained by the lack of polarity (no diE) and an inability to form H-bonds or, more generally, interactions involving a specific position in the molecule. Both n-heptane and toluene are only able to form van der Waals and π bonding interactions, respectively, that do not induce any significant position-specific isotope effect during distillation. Interestingly, the data for bromoethane (brom) do not fit with the 1 values either (Figure 2c) when it appears to have a behavior similar to alcohols concerning ET(30) trends (Figure 2b). Since the parameters 1 and ET(30) describe the potential interaction of the molecule as a solvent, [START_REF] Laurence | Solvatomagnetic Comparison Method: A Proper Quantification of Solvent Hydrogen-Bond Basicity[END_REF] it can be seen here that the molecule acts both as solute and solvent, i.e. in a self-association. Therefore, a molecule may be in a state to give or to receive H (or another interaction) but not to do both at the same time. Thus bromoethane shows a weak baseproperty, with a halogen-type association likely to be occurring.

As can be seen from Figure 2b, a clear double trend is observed. One encompasses the alcohols (al-x), but also chloroform (chlo), which are both H-acceptors and H-donors (a proton can readily be exchanged with deuterium [START_REF] Iwasaki | Hydrogen/deuterium exchange reaction between chloroform and water-d2 in two-liquidphase system[END_REF] ), and to some extent bromoethane (brom). The case of chloroform, with the isotope effect associated with its distillation correlating very well with those for alcohols (See Figure 2b), is interesting. This could be explained by (i) the high electron density in this molecule and (ii) the resulting presence of the acidic hydrogen. These two features of chloroform can interact forming a halogen bond and its strength seems to be dependent of the presence of a [START_REF] Meyer | Interactions with Aromatic Rings in Chemical and Biological Recognition[END_REF] C on the central carbon position. Trichloroethene (TCE) presents an isotope effect mostly located on the carbon bearing the acidic hydrogen that can form hydrogen bonds, but this compound better correlates with H-acceptors (Figure 2b). Of the remaining products studied, that are H-acceptors only, acetonitrile (acni), acetone (acet) and ethyl acetate (etac) all show good correlation with amphiprotic molecules (alcohols) (Figures 2 and3).

A change in the solvation shell of ethanol modifies the [START_REF] Meyer | Interactions with Aromatic Rings in Chemical and Biological Recognition[END_REF] C position-specific fractionation. As an illustration, the distillation of a dilute acidic aqueous ethanol solution (10% v/v with 1 M H 2 SO 4 ), using the same device, leads to the following isotope effect , measured on the distillate: bulk, +3.4‰, C-1, +6.3‰, C-2, +0.5‰. By increasing the H-bond network around the ethanol molecule, the  C-1 is accentuated, thus confirming that the isotopomer containing a 12 C at the C-1 position reinforces the intermolecular interaction. Interestingly, the reverse is found for the position C-2 for which it is the isotopomer with a 12 C on the CH 3 that distils first. The significance of such behaviors is beyond the scope of this paper, but it is clear that the 13 C position-specific fractionation is sensitive to the types of interaction and to the models of association: most probably a large network (complex) for alcohols and chloroform versus a binary association (simple) for the other compounds investigated.

Conclusion

The main outcomes of this study are two-fold. First, it is shown that position-specific fractionation in liquid-vapor transitions can occur in a range of molecules. When the molecular architecture within the liquid enables interactions, then the isotope effect of diffusion  diff-liq (ref 4) is position-specific. This is a novel finding, since up to now it has been considered that [START_REF] Meyer | Interactions with Aromatic Rings in Chemical and Biological Recognition[END_REF] C fractionation during diffusion shows a mass effect only, i.e. is not position-specific. These new data now need to be included in the model of fractionation during the process of volatilization according to the Craig-Gordon isotope model for organic liquids. [START_REF] Julien | Insights into Mechanistic Models for Evaporation of Organic Liquids in the Environment Obtained by Position-Specific Carbon Isotope Analysis[END_REF][START_REF] Kuder | Effects of volatilization on carbon and hydrogen isotope ratios of MTBE[END_REF] Second, [START_REF] Meyer | Interactions with Aromatic Rings in Chemical and Biological Recognition[END_REF] C position-specific fractionation appears to be exploitable as a marker of the strength of interactions in the liquid phase within the self-association framework. Fundamentally, the presence of a 13 C atom in a molecule according to its position decreases the interaction strength. How the presence of a [START_REF] Meyer | Interactions with Aromatic Rings in Chemical and Biological Recognition[END_REF] C acts to (de)stabilize the H-bond in alcohols now needs to be further investigated by theoretical modelling. Other intermolecular interactions such as the dimerization of carboxylic acids could also be studied using the same approach. Moreover, more complex systems such as VOCs diffusive transport in soil could be studied in the same way where the presence of heavy carbon isotopes should also have an influence on isotopic fractionation associated with diffusion.

8 / 29 13 C

 813 -irm-EA/MS Bulk[START_REF] Meyer | Interactions with Aromatic Rings in Chemical and Biological Recognition[END_REF] C abundance (δ[START_REF] Meyer | Interactions with Aromatic Rings in Chemical and Biological Recognition[END_REF] C bulk ) was determined by irm-MS using an Integra2 spectrometer (Sercon Instruments, Crewe, UK) linked to a Sercon elemental analyser (EA) (Sercon Instruments, Crewe, UK). A precise amount of each compound was weighted into tin capsules (2x5 mm, Thermo Fisher scientific) using a 10 -6 g precision balance (Ohaus Discovery DV215CD) to give approx. 0.4 mg of carbon for each compound. Great care was taken to ensure that there was no leakage from the capsule. First, during the weighing of the Volatile Organic Compound (VOC) introduced into the tin-capsule, the sealed capsule was left on the balance for a short delay to verify the stability of the mass. No change in mass indicated that the capsule was effectively sealed. Secondly, the percentage of carbon was checked by the operator by comparing this value with that usually obtained on the solid working reference. This gives a check in relation to the intensity of the signal of the ions obtained from the CO 2 generated from the different isotopologues (m/z = 44, 45, 46). Moreover, each sample is analyzed three times to ensure measuring true δ[START_REF] Meyer | Interactions with Aromatic Rings in Chemical and Biological Recognition[END_REF] C values and a repeatability study has also been performed in a previous study (see ref 3) which demonstrates the capability of these experimental conditions to measure the bulk isotopic composition of VOCs such as MTBE. The δ 13 C (‰) values were expressed relative

  This strongly reinforces the deduction of an important role of the 13 C isotopomer C-1 in the interaction between molecules.

Figure 1 .

 1 Figure 1. Definition of isotopologues and isotopomers: example of ethanol.

Figure 2 .

 2 Figure 2. Relationship between the isotope effect located on position C-1 (IE C-1 ) and (a) the relative permittivity (diE), (b) the hydrogen bond acidity (ET(30)), y=-0.31x+19.89, R 2 =0.853 for both H-acceptors and H-donors, y=-0.80x+34.69, R 2 =0.938 for H-acceptors only and (c) the hydrogen bond basicity (β 1 ), y=16.1x-6.37, R 2 =0.619.

Figure 3 .Figures Figure 1 .

 31 Figure 3. Relationship between the bulk isotope effect (IE Bulk ) and (a) the relative permittivity (diE), (b) the hydrogen bond acidity (ET(30)) y=-0.20x+13.39, R 2 =0.359 for both H-acceptors and H-donors, y=-0.36x+15.65, R 2 =0.643 for H-acceptors only and (c) the hydrogen bond

  Figure 2.

  d. 1 H/ 13 C dual +

probe, carefully tuned at the recording frequency of 100.61 MHz, or a Bruker AVANCE III 400 spectrometer fitted with a 5 mm i.d. BBFO probe tuned at the recording frequency of 100.62

MHz. The temperature of the probe was set to 303 K. Spectra were acquired without tube rotation. Isotope 13 C/

[START_REF] Reichardt | Solvatochromism, thermochromism, piezochromism, halochromism, and chiro-solvatochromism of pyridinium N-phenoxide betaine dyes[END_REF]

C ratios were calculated from processed spectra as described previously.

[START_REF] Julien | Fractionation in position-specific isotope composition during vaporization of environmental pollutants measured with isotope ratio monitoring by 13 C nuclear magnetic resonance spectrometry[END_REF][START_REF] Silvestre | Isotopic 13 C NMR spectrometry to assess counterfeiting of active pharmaceutical ingredients: Site-specific 13 C content of aspirin and paracetamol[END_REF]

  Table3reports the isotope effects ( bulk and  C-1 ), diE, ET(30) and 1 (see Introduction for a definition of these terms). Unfortunately, these parameters are not available in the literature for all the molecules studied, but sufficient values are available to establish trend lines. These lines have been drawn according to the visual distribution of data points, which means that Figures2 and 3do not present real correlations between isotope effects and the physical properties of compounds but, rather, lines drawn without

a priori to make the discussion easier. Trend lines have been drawn using Microsoft Excel so they have an equation and a root mean square (R²) which also are detailed in Figures captions in order to help the discussion the results. Figures

2 and 3

illustrate the potential trends between the isotope effects and diE, ET(30) and 1, respectively. A clear trend is observed between  C-1 and diE, ET(30) and a weaker for 1 (Figures 2), while a loose trend can even be discerned between  bulk and diE, ET(30) and 1 (Figures 3).

Table 1 .

 1 Molecular structure of each compound investigated with the carbon atoms numbered in decreasing chemical shift in the[START_REF] Meyer | Interactions with Aromatic Rings in Chemical and Biological Recognition[END_REF] C NMR spectrum. The code of each sample used in Fig.3 and4 and in the Tables is also defined.

	MTBE
	n -heptane

Table 2 .

 2 [START_REF] Meyer | Interactions with Aromatic Rings in Chemical and Biological Recognition[END_REF] C isotope composition  in ‰ and the difference between[START_REF] Meyer | Interactions with Aromatic Rings in Chemical and Biological Recognition[END_REF] C content of the starting material and of the distillate  in ‰ for bulk isotopologue and each isotopomer of the studied samples (see Table1for denomination and carbon numbering).

		Sample	Bulk	C-1	C-2	C-3	C-4	C-5
		δ 13 C initial -30.7	-30.7
	al-1	δ 13 C distillate -28.5	-28.5
			Δδ 13 C +2.2	+2.2
		δ 13 C initial -28.9	-32.2	-25.5
	al-2	δ 13 C distillate -24.8	-27.5	-22.0
			Δδ 13 C +4.2	+4.7	+3.6
		δ 13 C initial -32.6	-35.5	-28.1	-34.1
	al-3	δ 13 C distillate -29.2	-31.1	-25.1	-31.5
			Δδ 13 C +3.4	+4.5	+3.1	+2.6
		δ 13 C initial -26.9	-22.6	-29.1
	al-4	δ 13 C distillate -24.4	-17.6	-27.7
			Δδ 13 C +2.6	+5.0	+1.3
		δ 13 C initial -31.0	-42.7	-29.1	-24.9	-28.5
	al-5	δ 13 C distillate -28.3	-37.7	-26.1	-22.4	-25.6
			Δδ 13 C +2.8	+5.0	+3.0	+2.6	+2.9
		δ 13 C initial -25.4	-17.1	-28.2
	al-6	δ 13 C distillate -23.9	-15.0	-26.9
			Δδ 13 C +1.5	+2.1	+1.3
		δ 13 C initial -31.2	-20.6	-28.0	-28.5	-50.5
	al-7	δ 13 C distillate -28.8	-15.1	-26.3	-27.1	-48.3
			Δδ 13 C +2.4	+5.5	+1.7	+1.4	+2.1
		δ 13 C initial -28.3	-36.5	-23.9	-27.8	-26.2	-27.1
	al-8	δ 13 C distillate -25.5	-32.4	-21.3	-25.3	-24.2	-24.4
			Δδ 13 C +2.8	+4.1	+2.6	+2.5	+2.0	+2.7
		δ 13 C initial -24.4	-39.0	-17.2	-18.9	-26.9
	al-9	δ 13 C distillate -22.2	-33.9	-17.8	-17.4	-23.0
			Δδ 13 C +2.2	+5.1	-0.5	+1.5	+3.9
		δ 13 C initial -41.7
	chlo	δ 13 C distillate -33.4
			Δδ 13 C +8.4

Table 2 .

 2 (

	507	continued)

Table 3 .

 3 Isotope effect  in ‰ for bulk and isotopomer C-1, relative permittivity (dielectric constant) diE, solvent hydrogen-bond acidity ET(30) and hydrogen-bond basicity  1 of the studied samples (see Table1for denomination).

	code	IE bulk	IE C-1	ET(30)	β 1	dielectric constant
	al-1	+2.3	+2.3	55.4	0.54	33.0
	al-2	+4.3	+4.8	51.8	0.62	24.5
	al-3	+3.5	+4.6	50.5	0.65	20.3
	al-4	+2.6	+5.1	48.4	0.68	20.3
	al-5	+2.9	+5.2	49.7	0.67	17.5
	al-6	+1.5	+2.2	nd	0.73	12.5
	al-7	+2.5	+5.6	41.5	0.71	nd
	al-8	+2.9	+4.3	49.3	0.70	nd
	al-9	+2.3	+5.2	47.2	0.73	nd
	chlo	+8.7	+8.7	39.1	nd	4.8
	brom	+5.1	+8.6	37.8	0.12	9.4
	TCE	+5.1	+6.6	35.9	nd	nd
	acet	-0.2	+1.8	42.3	0.49	21.1
	2hep	+1.3	-0.2	nd	nd	nd
	3hep	+1.8	+1.1	nd	nd	nd
	4hep	+1.9	+7.1	36.0	nd	nd
	acni	-0.7	-2.5	45.6	0.37	36.6
	etac	+3.2	+4.1	nd	0.52	6.1
	MTBE	+2.6	+7.1	33.0	nd	nd
	TAME	+3.5	+6.7	nd	nd	nd
	tol	+3.2	+3.9	33.9	0.15	nd
	hept	+2.5	+2.1	30.9	0.00	nd

Acknowledgments

This work is funded by the French National Research Agency (ANR), project ISOTO-POL in the program CESA (N° 009 01). M. Julien thanks the Région Pays de la Loire for funding his postdoctoral bursary through the project PLAISIR. We thank Ms. Anne-Marie Schiphorst for help with irm-MS measurements and Jérôme Graton and Jean-Yves Lequestel for fruitful discussions on inter-molecular interactions.


E Beauvier

S Bodea

A Pocheau

Front propagation in a vortex lattice : dependence on boundary conditions and vortex depth †

We experimentally address the propagation of reaction-diffusion fronts in vortex lattices by combining, in a Hele-Shaw cell and at low Reynolds number, forced electroconvective flows and an autocatalytic reaction in solution. We consider both vortex chains and vortex arrays, the formers referring to mixed free/rigid boundary conditions for vortices and the latters to free boundary conditions. Varying the depth of the fluid layer, we observe no variation of the mean front velocities for vortex arrays and a noticeable variation for vortex chains. This questions the two-dimensional character of front propagation in low Reynolds number vortex lattices, as well as the mechanisms of this dependence.

Introduction

Stirred flows are well-known to enhance the transport properties of a medium by their ability in intimately mixing components. This feature extends to reacting flows, as the net reaction time of a medium is found to sharply decrease with the stirring rate. In particular, in regimes where reaction proceeds by the propagation of fronts separating reacted and unreacted zones, the mean front velocity is found to largely increase with the vortex intensity. Important consequences follow in the fields of combustion 1 , chemistry [START_REF] Field | Oscillations and travelling wave in chemical systems[END_REF][START_REF] Epstein | [END_REF] or even ecology [4][5][6] . However, the actual mechanisms of front velocity enhancement as well as the whole factors on which this velocity depends call for a better understanding in both turbulent [START_REF] Williams | Combustion Theory[END_REF][7][8][9][10][11][12] and laminar [13][14][15][16][17][18][19][20][21] contexts.

We focus attention here on laminar single scale flows made of a periodic planar array of vortices. Propagation of fronts in this context has been addressed by theories [13][14][15][16] and numerical simulations [14][15][16][17] by considering a two-dimensional array of vortices involving equal width and length. The same boundary conditions (hereafter called b.c.) were thus displayed on all their separatrices (actually the free b.c. provided by neighbor vortices) and the set of separatrices was found to form a square array. In comparison, apart from a study in Taylor vortex flows 21 , experiments [18][19][20] addressed a single chain of vortices confined in a channel. Then, different boundary conditions were applied on vortex separatrices depending on whether they were a channel wall (rigid b.c.) or a frontier with † Electronic Supplementary Information (ESI) available: Movies of front propagation in a channel accelerated 25 times. Movie "FreeBC" : large channel (5 vortices between walls) and free effective boundary conditions, U/V 0 = 35. Movie "RigidBC" : narrow channel (1 vortex between walls) and rigid effective boundary conditions, U/V 0 = 27. See DOI: 10.1039/b000000x/ a Aix Marseille Univ, CNRS, Centrale Marseille, IRPHE, Marseille, France a neighbor vortex (free b.c.). In addition, all experiments naturally took place in a fluid layer whose depth d was thus a natural parameter of them.

Usually, experiments involved a single vortex depth d. However, one of them varied this parameter and concluded that the mean front velocity V f not only depends on the vortex intensity but also on the vortex depth d 19,20 . As it was performed in a vortex chain, it involved mixed free/rigid boundary conditions on vortices. Its vortices were also longer than large, so that the set of separatrices was rectangular instead of square 19,20 . Accordingly, this experiment involved two kinds of asymmetry : (i) a boundary asymmetry yielding different b.c. on a vortex, (ii) a geometric asymmetry induced by the difference between vortex width and vortex length.

As theories or simulations considered a two-dimensional planar flow, they implicitly assumed a 2-dimensional character of front propagation in a vortex lattice [13][14][15][16][17] . However, the dependence of V f on d questions this assumption and raises the following interrogations regarding its possible origin : could it follow from the geometric asymmetry of vortices or from their boundary asymmetry and then, by what kind of mechanism ? In particular, could it refer to a modification of the structure of a still planar flow or to the occurrence of secondary flows breaking the planar flow assumption ?

To answer these questions, we perform here front propagation experiments in an array of vortices involving equal width and length and either similar or different boundary conditions on their separatrices. We then vary the channel depth so as to highlight the implications of geometrical symmetry and of boundary symmetry on the dependence of front propagation on channel depth. Altogether, this provides the first comparison of the effects of boundary conditions and of vortex depths on front propagation in a stirred medium.

1-8 | 1

The experimental set-up allowing these configurations to be achieved is reported in section 2. We then address qualitative and quantitative results on section 3 and discuss them in section 4. A conclusion about the study and its general issue is drawn in section 5.

Experiment

The experiment is conducted in a fluid enclosed in a channel with top, bottom and lateral boundaries and whose depth d is imposed by spacers (Fig. 1-a). Three depths have been considered : d = 2, 4 or 6 mm.

Flows are induced by electroconvection thanks to both a pattern of magnets placed below the channel and an electric current imposed by electrodes placed at the extremities of the fluid layer (Fig. 1-a). All the magnets are squares with dimension L = 20mm. They are placed one beside the others so as to form a square lattice. Their magnetization direction alternates between neighbors. It thus induces in the fluid layer a periodic magnetic field and thus periodic Laplace forces which give rise to alternating vortices. Electromagnetic forcing has been used for a while in many different contexts ranging from boundary layer effects in salty water 22,23 , mixing 24,25 , single scale flows [26][27][28] multi-scale flows 29 or quasi-two dimensional turbulent or chaotic flows [30][31][32] . Beyond hydrodynamics, it has also been used to study front propagation [18][19][20] in a laminar regime. In particular, in the Stokes regime in which our experiment fits, the current density, the magnitude of the Laplace force density and then the vortex intensity are all proportional, thereby providing a simple mean to tune the flow magnitude by the current intensity.

Front propagation is provided by an autocatalytic oxydoreduction reaction between chlorite and iodure ions in the fluid layer :

ClO - 2 + 4I -+ 4H + → Cl -+ 2I 2 + 2H 2 O (1) 5ClO - 2 + 2I 2 + 2H 2 O → 5Cl -+ 4IO - 3 + 4H + (2) IO - 3 + 5I -+ 6H + → 3I 2 + 3H 2 O (3)
The autocatalytic feature is provided by the fact that the consumption of iodide I -worked out in (1) is enhanced in (3) by a product of reaction, the iodate IO - 3 produced in (2) thanks to the iodine I 2 provided by the consumption of iodide itself in (1) and ( 3). This provides a positive feed-back of iodide consumption on itself possibly yielding a run-away and, in specific regimes, a reaction zone localized in a thin interface, the front, standing in between reacted and unreacted zones [START_REF] Saul | Oscillations and travelling Waves in chemical system[END_REF][START_REF] Hann | [END_REF] . The detection of these zones, and thus of the front as their frontier, is provided by a change of colour as reaction proceeds. For this, a Starch indicator introduced in the solution forms a blue complex with iodide and iodine which disappears after reaction following the depletion of iodide. The unreacted blue zone then transforms into an uncolored reacted zone. The front velocity may thus be deduced from the kinematics of their frontier, as recorded on a computer by a camera (Fig. 2). In a fluid at rest, this velocity is labelled V 0 and amounts in our experiment from 1.0 to 1.4mm/mn depending on the imposed current density 19,20 . For this, the temperature has been stabilized to 20 • C to minimize the variability of V 0 to less than few percents.

When the medium is stirred by the vortex lattice, front propagation may take advantage of advection. This is in particular so when the front is enrolled by vortices. Then, following the large value of the flow velocity compared to the proper front velocity V 0 (from 5 to 40 times greater here), the front velocity stays close to it. In particular, it appears that at large flow velocity, after having crossed a vortex separatrix, the front is advected on a streamline (Fig. 2) where its mean velocity provides a direct information on the vortex intensity U. Thanks to the linear relationship between flow velocity and current density, this determination may then be extended on the whole range of the experiment.

Considering a planar vortex flow V, its advective term V.∇V may be estimated to about U 2 /L in the bulk since the variation scale in the flow plane is then L. However, close to a lateral boundary, this scale decreases to about the channel depth d, yielding another estimate of about U 2 /d there. In both cases, the dissipative term ν V may be estimated to νU/d 2 for a near Poiseuille profile where ν ≈ 1 mm 2 /s denotes the water viscosity. Accordingly, the Reynolds number reads Re = U 2 d 2 /νUd ≡ Ud/ν close to a boundary and Re = Re d/L < Re in the bulk. As L = 20mm, one gets in our experimental range U ≤ 40mm/mn, 2 ≤ d ≤ 6, a maximal Reynolds number Re = 1.2 in the bulk and Re = 4 close to boundaries. Accordingly, the experiment may be considered to stand within or close to the Stokes regime.

The channel is 400 mm long and either 20, 80 or 100 mm large. In particular, its width could then include 1, 4 or 5 magnets. Placing the channel sides on magnet boundaries, this provided 1, 4 or 5 vortices in between them (Fig. 1-b,c). As these vortices display equal length and width, they satisfy the geometrical symmetry. When a single vortex fits in between the channel width, it undergoes rigid b.c. on the separatrices made by the channel sides and free b.c. on the remaining ones made by the frontier with neighbor vortices. This mixed case corresponds to asymmetric boundaries (Fig. 1-b). In contrast, when the channel width is large compared to the vortex size, the experiment shows that front propagation in the central vortices is favored, so that they act as a leading zone for the overall propagation (Fig. 2-a). As these central vortices are surrounded by other vortices, their b.c. are free on all their separatrices. This corresponds to symmetric boundaries (Fig. 1-c).

This way, the set-up is able to provide an experimental achievement of front propagation with geometric symmetry and either symmetric or asymmetric boundaries.

Fronts and Velocity

We report typical evolutions of fronts in the vortex lattice together with quantitative measurements of their mean front velocity V f along the channel axis, at various vortex intensity U and channel depths d. The two situations of free b.c. (symmetric vortices) and mixed free/rigid b.c. (asymmetric vortices), are considered to allow for their comparison. Here, the channel width englobes 4 or 5 vortices, so that 2 or 3 central vortices undergo free b.c. on all their separatrices (Fig. 1-c). A typical evolution of reaction fronts is reported on figure 2-a. It appears that the fronts display a periodic motion involving alternating phases between a slow crossing of separatrices and a quick advection in between. In addition, the contamination of a new line of vortices always occur on central vortices so that the mean front propagation is actually set by symmetric vortices. Figure 3-a reports, at various channel depths, the enhancement of the mean front velocity V f by the flow in term of nondimensional variables V f /V 0 , U/V 0 . It appears no definite dependence on the channel depth d, the velocity enhancement displaying the same linear trend in this velocity range. This denies relevance the channel depth d on the mean front propagation, in agreement with a two-dimensional view of this propagation issue. Figure 3-b reports, at various channel depths, the enhancement of the mean front velocity V f by the flow in term of non-dimensional variables V f /V 0 , U/V 0 . It shows evidence of a dependence on the channel depth d, the velocity enhancement displaying an increasing concavity with the channel depth. This reduction of velocity enhancement at large d yields in particular, at U/V 0 ≈ 26, a decrease of about 40% of the mean front velocity V f when increasing the channel depth from d = 2mm (V f /V 0 = 11.7) to d = 6mm (V f /V 0 = 6.8). This evidences a noticeable influence of the channel depth on front propagation, a fact which questions the two-dimensional assumption considered in theories and simulations [13][14][15][16][17] .

Discussion

Experimental data on velocity enhancement by vortices reveal an ambivalent behavior regarding the dependence on the channel depth : no dependence for free b.c. ; a noticeable dependence for mixed free/rigid b.c.. As these data have been obtained on vortices involving equal length and width, this indicates that their dependence on d has to be traced back not to geometrical asymmetry 19,20 but to the asymmetry between the boundary conditions applied at the vortex periphery.

For all vortices, rigid b.c. are in order at the top and bottom plates. However, for mixed free/rigid b.c., rigid b.c. are also undergone at two opposite separatrices of vortices. We discuss the implications of these features on front propagation and especially its evolution with the channel depth. The first feature is reminiscent of a Poiseuille flow and may be thought to be satisfied for electroconvective flows at the dominant order. The second feature follows from the local relationship between front velocity and flow, as studied on a Poiseuille flow 35,36 . It was then found that front velocity turns out to be the sum of the proper front velocity V 0 plus either the mean flow velocity on the channel depth for low enough d compared to the front thickness (i.e. η = dV 0 /2D < 1) or the maximal flow velocity in the opposite case 35,36 , D = 2.10 -3 mm 2 /s denoting the molecular diffusivity of the relevant species. In all cases, front velocity appears as the sum of V 0 and a velocity linearly related to U by a constant prefactor. In particular, in our experiment, as η is larger than 8, fronts locally advance at velocity V 0 + U(x, y).n where n denotes the front normal in the mid-depth plane z = d/2 and where the z-profile has been normalized so that f (1/2) = 1.

Altogether, both features may thus be thought as relevant for explaining the independence of mean front velocity on d.

Asymmetric vortices : mixed free and rigid b.c.

In comparison to the above case of symmetric vortices, the dependence of mean front velocity on d displayed when some vortex boundaries are rigid calls for either a structure of planar flows dependent on d or a breaking of the planar flows assumption.

Dependence of planar flows on d.

This issue has been explored in a previous study where the mean front velocity was argued to correspond to the quickest trajectory of the front head in the vortex field 19,20 . Then, an optimization model based on a simplified family of trajectories inspired by observations was worked out to derive the mean front velocity within a prescribed modeling of planar flows. In particular, these flows involved circular streamlines and boundary layer thicknesses of half the channel depth. The latter feature was motivated by the fact that, in a Hele-Shaw cell, the channel depth monitors the width of boundary layers, following the dominant dissipation implied by the presence of parallel plates. Its value was then taken to be d/2, yielding the quickest streamline to be located at this distance to the channel wall. This constraint thus led the flow field to depend on d. Its implication on the mean front velocity was deduced by selecting the trajectory that minimizes the crossing time through a vortex.

The predictions of the optimization model for geometrically symmetric vortices are reported on figure 3-b for the three channel depths. They are found to satisfactorily recover the experimental data, especially the increasing concavity of the enhancement curves with the channel depth. Beyond this agreement, it should be realized that the optimization model is based on a phenomenological family of front trajectories and on a crude flow modeling. In particular, the circular form of streamlines stands for an approximation of actual elliptical streamlines, which prevents the flow velocity to vary on them. In this sense, the streamline velocity is intended to reflect the average velocity on it instead of its local values. It thus corresponds to the nature of the velocity measured experimentally and is sufficient to evaluate the advection time of the front head on a streamline. However, this approximation also implies circular streamlines which are obviously inadequate around the square vortex separatrices, especially around the vortex corners. Finally, the scaling of the boundary layer thickness with the channel depth, although physically sound,

1-8 | 5
involves an approximate prefactor.

For these reasons, the optimization model must be taken as a toy model whose merit is to provide a kinematic origin for the dependence of front velocity on the channel depth : as the quickest front trajectory proves to belong to the boundary layer 19,20 , front propagation appears to be sensitive to its structure and thus to the channel depth which modifies its thickness. However, beyond this qualitative relevance, the optimization model appears too crude for allowing a quantitative comparison. In particular, as its major outcome is the sensitivity of front propagation to the structure of the boundary layer, it actually calls for an accurate flow modeling to provide a relevant quantitative comparison. In addition, as it relies on a restricted family of trajectories, its predicted mean flow velocities would nevertheless be approximate. Accordingly, given an accurate model of planar flow, a relevant quantitative comparison to experiment would call either for an extension of the optimization model to any possible trajectory or for numerical simulations.

Breaking of the planar flows

For rotating flows, a generic mechanism for the generation of secondary flows is provided by the Ekman pumping 37 . It relies on the fact that viscous effects at boundaries make velocity and centrifugal force shrink whereas pressure is maintained. This therefore induces a force disequilibrium in the fluid which results in setting the fluid into motion. Another possible origin for secondary flows is the occurrence of hydrodynamical instabilities but these are not expected at the low Reynolds numbers at which the experiment stands (Re = Ud/ν < 4 here ; see section 2). An important difference between these two mechanisms should be noticed : whereas the Ekman pumping must be present at all Reynolds numbers or at none, instabilities manifest themselves above some onset. Accordingly, if secondary flows are generated by Ekman pumping, their amplitude should continuously grow with the vortex intensity without onset.

In our configuration, this phenomenon underlying the Ekman pumping might first occur at the top and bottom rigid plates of the channel. It should thus induce inward flows towards the vortex centers giving rise, by incompressibility, to flows directed normally to the plates and finally to two spiral flows, following the top-bottom symmetry of the configuration. These secondary flows should however be equally present on symmetric and asymmetric vortices, i.e. for free or mixed free/rigid b.c., since both configurations involve top and bottom plates. No spirals have however been detected from front motion in our experimental range so that the corresponding secondary flows are presumably too faint to noticeably modify front dynamics. Anyway, they could not solely explain the specificity of the asymmetric vortex configuration, i.e. its noticeable dependence on the channel depth.

In contrast, the presence of rigid lateral b.c. stands as a specificity of asymmetric vortices from which specific inward flows could be generated. They might in particular modify the advection of the front head all along its travel inside the lateral boundary layer and possibly slightly deflect its trajectory towards the vortex center. This could result, at the end of the advection phase, in a slight increase of the front distance to the neighbor vortex separatrix. As the front must cross it by its sole propagation at its slow proper velocity V 0 , this effect, although faint, could result in a noticeable decrease of the mean front velocity. This mechanism is however in competition with the one induced by the increase of the boundary layer thickness of planar flows, since both yield an increase with d of the distance of the front head to the neighbor vortex separatrix. Determining the role of each thus calls for a quantitative evaluation of their respective effects.

Beyond the qualitative explanation of secondary flow generation, it should be noticed that a quantitative determination of their amplitude requires first an accurate modeling of the primary flows, i.e. of planar flows. Accordingly, the discrimination between the two above mechanisms calls for an accurate modeling of planar electroconvective flows in both symmetric (free b.c.) and asymmetric vortices (mixed free/rigid b.c.). Then, the most direct way to determine the implications of these flows on front propagation will be to perform numerical simulations. Comparison with experiments could then be able to conclude whether the planar flows for free b.c. as well as their modification by lateral walls allow to recover both the independence to the channel depth d for symmetric vortices as well as the dependence on it for asymmetric vortices.

Conclusion

In theories and simulations, front propagation in laminar flows is addressed on two-dimensional flows involving a prescribed structure [13][14][15][16][17] . This assumption actually denies a relevant influence of three-dimensional features of vortices as their depth or asymmetries induced by boundary conditions or by direct forcing. However, a previous experiment performed on a vortex chain enclosed in a channel revealed a variation of the mean front velocity with the channel depth 19,20 . This observation thus calls for determining to what extent front propagation in vortex arrays may be considered as 2-dimensional and to what kind of features it may be dependent.

To address this issue, we performed experiments so as to draw the first comparison of front propagation for different flow b.c. : the free b.c. considered in theories or simulations and the mixed free/rigid b.c. involved in the experiment which evidenced a dependence on the channel depth. For this, we considered both narrow and large channels as compared to the vortex width. Large channels provided effective free b.c., as they leave place for several vortices between their walls (vortex array). Narrow channels provided mixed free/rigid b.c. for vortices, as they englobe a single line of vortices between their walls (vortex chain). In both cases and in contrast with a previous experiment 19,20 , we considered only geometrically symmetric vortices involving equal width and length, as considered in theories and simulations.

The dependence of mean front velocity on channel depth d has been found to be bivalent : no dependence for vortex arrays (free b.c.) but a noticeable dependence for vortex chains (mixed free/rigid b.c.). In particular, in the latter case, the mean front velocity was found to increase as d was reduced, up to eventually reach at low d the values found on vortex arrays. Thus, in contrast with primary intuition, the most dissipative configuration regarding hydrodynamics (vortex chain at low cell depth) provided a similar efficiency regarding propagative transport as the less dissipative one (vortex array). This might provide a mean for restoring in confined systems a propagation efficiency similar to that displayed in unconfined ones. In particular, this property might apply to chemical engineering in stirred systems, possibly even down to microfluidic scales 38 . It would then stress that the net reaction efficiency in narrow geometries be monitored by their depth, up to eventually recover in thin configurations the efficiency involved in extended ones.

On the other hand, on vortex arrays, the independence of front velocity from the channel depth shows that the vertical profile of flows plays no significant role regarding propagation efficiency. This property might apply to natural situations which proceed by reaction-diffusion-advection within a set of quasi two-dimensional vortices, such as in plankton dynamics 4,5 or in ozone hole evolution 6 . It would then point out the planar structure of their vortices as the essential feature regarding the global dynamics.

Finally, the analysis of the possible origin of the dependence of front propagation on channel depth in vortex chains (mixed free/rigid b.c.) has stressed two possible mechanisms : a modification of planar flows by boundary conditions or a breaking of the planar flow assumption by secondary flows. In the first case, front propagation would still be considered as a two-dimensional issue, sensitive enough for being influenced by b.c. ; in the latter case, it could no longer be considered as two-dimensional since noticeably dependent on faint threedimensional features. Discrimination between these mechanisms calls for an accurate modeling of electroconvective planar flows together with simulations of front propagation, so as to provide relevant comparisons to experiments.

2 | 1 - 8 Fig. 1

 2181 Fig. 1 Sketch of the set-up showing the chemical solution stirred by Laplace forces induced by an electric current flowing over alternate magnets (a) and the resulting chain (b) or array (c) of counter-rotating vortices. The magnetic field and the electric current are labelled B and I and the channel depth and the vortex size d and L respectively. Vortex boundary conditions are mixed free/rigid (b) or free (c), depending on whether vortices make contact with the channel sides (b) or are surrounded by other vortices (c ; except for side vortices).

3. 1

 1 Symmetric vortices : free b.c.

3. 2

 2 Asymmetric vortices : mixed free/rigid b.c. Here, the channel width englobes only a single vortex which thus undergoes mixed free/rigid b.c (Fig.1-b). Regarding front motion, the typical evolution of reaction fronts reported on figure 2-b shows that the quick advection phase occurs near the channel sides. As in the previous symmetric situation, fronts undergo a periodic motion made of a succession of slow crossing of separatrices and quick advection in between, but the latter occurs close to a rigid b.c. this time.

Fig. 2

 2 Fig. 2 Front propagation in a vortex array (a) or in a vortex chain (b). Times regularly increases from left to right while fronts propagate upwards. The vortex size is L = 20mm, the channel depth is d = 2mm and the channel width is labelled l. (a) Free boundary conditions, except for the side vortices. l/L = 5, U/V 0 = 17.3. Time interval between pictures : 65s. (b) Rigid lateral boundary conditions. l/L = 1, U/V 0 = 14.6. Time interval between pictures : 58s.

4. 1

 1 Symmetric vortices : free b.c. The rigid b.c. on the top and bottom plates of the channel imply a near parabolic profile of the flow field on the channel depth direction z. Nevertheless, front propagation shows no dependence on its extent d. This feature is coherent with : 1. planar flows involving a common structure U independent of z and d and an intensity simply with z on the depth axis : V = f (z/d)U, the flow U belonging to the plane (x, y) and being only dependent on its variables : U ≡ U(x, y).

2 .

 2 a front propagation monitored by the planar flow U, independently of the z-profile f (z/d).

Fig. 3

 3 Fig. 3 Reduced front velocity V f /V 0 as a function of the reduced vortex intensity U/V 0 . (a) Free b.c. : a single linear relationship is displayed whatever the channel depth d. (b) Mixed free/rigid b.c. : concave relationships are displayed with a concavity increasing with the channel depth d. Lines correspond to predictions from an optimization model 19,20 .

| 1-8

Acknowledgments

This work has been carried out in the framework of the Labex MEC (ANR-10-LABX-0092) and of the A*MIDEX project (ANR-11-IDEX-0001-02), funded by the Investissements d'Avenir French Government program managed by the French National Research Agency (ANR).