Groundwater sampling, purge or no purge that is the question

Boisson A , Guérin V., Vergnes J-P., Quiot F. BRGM, 3 Av Claude Guillemin, 45060, Orléans, France INERIS (ARDEVIE), Domaine du Petit Arbois, BP 33, 13 545 Aix-en-Provence Cedex 04

Summary

There are uncertainties about to what extend the water in the well is representative of the water in the aquifer. So usually, for groundwater sampling a purging step is considered. This purging step used to consider a number of volumes purged (from 3 to 5 times the well's volume) or more recently the stabilization of the physicochemical measures before the sampling itself. This purging step is nowadays done with flow rates less important that they used to be, due to the emergence of the low flow procedure (Puls and Barcelona, 1996) and the availability of new pumps. This purge step presents several constraints [START_REF] Robin | Field Evaluation of Well Purging Procedures[END_REF], among which the time (and thus cost) and the treatment/ elimination of the water purged. In France, a new standard (AFNOR, NF X 31-615) will be published in 2017. It states that in some conditions the purge can be avoid. Moreover some questions remain on the spatial representativeness of the water sampled regarding the protocol used for the purge step (flow rate, position of the pump in the well screen interval). Considering also that in many wells there are very small differences in loads within the aquifer that can generate ambient fluxes in wells [START_REF] Zinn | Effects of intraborehole flow on groundwater age distribution[END_REF] which may have an impact on the results of the deployment of passive samplers compared to classic sampling.

The representativeness of the sample taken is of paramount importance because it has a direct influence on the interpretation of the quality of the water and / or the pollution present. It is therefore necessary to know precisely which part of the aquifer is considered and how to proceed (purge or not). A project financed by the French Ministry of Environment between two public research bodies (BRGM and INERIS) has the objective to state the influence of the purge on the sampling representativeness.

To answer those questions, three means will be used: literature review to know the position in other countries, on-site measurements campaigns for different hydrogeological context of the influence of some parameter (no purge, purge at different flow rate, at different positions in the well) and modelling of flows in and near the well for different scenarios (including the modelling of real cases) stating origin of the pumped water. This study will allow to evaluate the interest and the consequences of the purging processes during the realization of qualitative samples of water and will allow the production of guidelines for different cases.

Context

The realization of purges on several well's volumes induces different constraints [START_REF] Robin | Field Evaluation of Well Purging Procedures[END_REF]:

 High pumping times and/or high flow rates in the case of large diameter boreholes and long screened ones,  Dewatering of the well in cases of poorly permeable medium,  Degassing in the case of volatile compounds,  Lack of assurance regarding the effective remobilization of stagnant water,  Management of contaminated water in the event of major pollution.

In addition to the above-mentioned technical constraints, well structure, aquifer structure, pollution storage processes and pollution type may have an impact on the representativeness of the samples.

Thus there is a need to better understood consequences of the purge on the representativeness of the water sampled and the compulsory nature of this step.

State of the art on the purge procedure

In the literature, according to authors and standards, for purging procedure, purged volumes equal to 3 to 5 well volumes (US EPA, 2013) or 2 to 10 (AFNOR, 2000) are commonly used. However, this approach is considered too generic and many authors suggest that purging to stabilize one or more relevant criteria is preferable (eg [START_REF] Daughney | Field Trial to Support Development of a Nationally Standardised Protocol for State of the Environment Groundwater Sampling in New Zealand[END_REF]. In the new French Standard (NFX31-615, 2016), it is the stability of the physico-chemical parameters that is favored. This is defined by acceptable deviations on three measurements at 3 to 5 min minimum interval. This implies a purge time of at least 10 minutes. This purging must also meet requirements of flow rate that do not induce an important drawdown. The aim is both to ensure a renewal of the drilling water, but also to ensure the representativeness of the water withdrawn and to avoid significant vertical flows. The stability of physicochemical parameters may be difficult to obtain in the case of qualitative heterogeneity of the groundwater. In this case, a volume equivalent to 5 well volumes (Vp) can be retained as a condition for stopping the purge. The objective is to limit the volume purged and to maintain a punctual aspect of the sample.

The consequences of the purge could be different due to several considerations:  Linked to the well structure as the length of the screened slot. Very small load differences within the aquifer can generate ambient fluxes in wells [START_REF] Elci | Implications of observed and simulated ambient flow in monitoring wells[END_REF][START_REF] Zinn | Effects of intraborehole flow on groundwater age distribution[END_REF], a long screened slot will tend to increase these fluxes.  Pump's position along the casing. Although at long times the position of the pump has no longer any influence on the water composition, the purge times generally used in the sampling of industrial sites may not be sufficient to attain this equilibrium.  Influence of water table and seasonal variations. The hydrodynamic conditions within an aquifer are not always stable over time, mainly in the case of free table groundwaters, and local or regional changes in flows are frequent.  Dynamics of flows under ambient conditions and pumping conditions:

o One of the strong constraints on the representativeness of the samples is the difference of dynamics of the well under natural conditions and under pumping conditions. Under natural conditions, the flows are controlled by the loads with hydraulic loads flowing from the highest to the lowest. This can generate ambient flows within the well. o Under pumping conditions, the flows arriving at the pump (lowest load) are then weighted by the transmissivities of the various hydrogeological units and the heterogeneities within them. In this case, the sample is an average sample weighted by the flows. Thus, pumping will always tend to preferentially sample the most transmissive zones of the aquifer (eg [START_REF] Hutchins | Ground Water Sampling Bias Observed in Shallow, Conventional Wells[END_REF]Acree 2000, Mayo, 2010). and is therefore not representative as such. An element stored in an area with little transmissivity or a high sorption capacity will be systematically undervalued.

It is particularly the case in fractured aquifer or in aquifer with or clayed or peaty layers.

 Influence of ambient vertical flows on samples. The presence of vertical natural ambient fluxes within boreholes is regularly observed in heterogeneous aquifers (eg [START_REF] Alazard | Investigation of recharge dynamics and flow paths in a fractured crystalline aquifer in semi-arid India using borehole logs: implications for managed aquifer recharge[END_REF]. On the other hand, they are often neglected in the context of supposedly homogeneous aquifers.

However, various studies have shown the ubiquity of these flows [START_REF] Elci | Implications of observed and simulated ambient flow in monitoring wells[END_REF] even in cases where the aquifers are supposed to be homogeneous. The creation of wells induces the creation of a well can induce an ambient flow even in a homogeneous medium [START_REF] Zinn | Effects of intraborehole flow on groundwater age distribution[END_REF].

These considerations are particularly important in the context of the use of passive samplers or low flow purges.

The distribution of the pollution plays also a strong role. For example in a highly stratified environment with highly heterogeneous layer chemistry, [START_REF] Mayo | Ambient well-bore mixing, aquifer cross-contamination, pumping stress, and water quality from long-screened wells: What is sampled and what is not?[END_REF] shows the difficulty in obtaining a representative sample. In this particular case, a volume corresponding to 1600 times the volume of the well was necessary in order to obtain a complete purge allowing a representative sample to be obtained.

No purge is possible in the case of a well that is in a zone transmissive enough so that a renewal of the water in the screened casing occurs naturally. This hypothesis makes it possible to get rid of the purge. In this case, a sampling a double sealing sampler may be sufficient to obtain a good representativeness of the quality of the aquifer at a point. This renewal can be estimated by analytical solutions or numerical modeling with some simplifications. Experimental evaluation of this turnover rate can be carried out by injecting solutes [START_REF] Pitrak | A Simple Borehole Dilution Technique in Measuring Horizontal Ground Water Flow[END_REF] or a controlled temperature [START_REF] Read | Characterizing groundwater flow and heat transport in fractured rock using fiber-optic distributed temperature sensing[END_REF] into the well and evaluating the dilution at different points.

Field results:

Test on a calcareous fractured media, long well In this well, physico-chemistry of the water changes around one Vp and concentrations increases with the purge volume (Figure 2). There are huge differences in concentrations between the water obtained with a discrete interval sampler and by pumping.

This has been attributed to functioning of the well. In ambient condition, for high water level the water in the well comes from a fracture before 20 m and when a pumping is active most of the water comes from a deep fracture around 41 m (see Figure 1). For low water level, only the deeper fracture is active.

Test on an alluvial porous media, short well

On an alluvial aquifer with an acid tar contamination due to coke production, measurements were carried on a well screened from 8 m to 14 m for 2 hydrogeological conditions with the same procedure as the former case study but with purging depth: 6; 8,5 and 11.

The effect of the purging rate seems to be connected to the hydro-geological context and the history of the contamination. The monitoring well studied is located on the fringe of the plume, so increasing the purging rate modified the ratio of non-polluted water pumped over polluted water pumped and furthermore variation of water flow direction with hydro-geological context, explain time variation in PAHs concentrations, even though purging and sampling protocols are unchanged.

For low water table condition, increasing the purging rate decrease the concentrations in PAH, and for each purging rate, a decrease of the concentration is observed. For each of these dynamic samples, logs of electrical conductivity of the water were made during the purge (Figure 5). The sampling with the discrete interval sampler gives result consistent with the log of the conductivity at ambient condition. The evolution of the conductivity during pumping at 0.02 m3/h with pakers at 30 m, shows a slow increase and then a fast relative stabilization to reach the conductivity of 5 mS / cm. An increase to this value was not expected. If the fluxes were purely horizontal, a value of around 3.5 mS / cm should have been measured at rest by logging at this depth. This value higher than expected could be due to a default in the confinement by the paker, to a piston effect when introducing and then remove of the apparatus and or due to water bypass due to gravel pack. The sampling at 30 m, with a pumping rate of 0.15 m3/h, has a conductivity value of 13.7 mS / cm, which is significantly higher than that observed under ambient conditions (logs or pneumatic sampler) or dynamic at low flow rates with shutters. The increase in the stress induces an increase in the conductivity due to a greater stress from the deep zones to the higher conductivities. For a high pumping rate (1 m 3 /h), the two samples have relatively similar values of 21 mS / cm at 38 m and 20.2 mS / cm at 30 m. This confirms the tendency for the system to be homogenized in the long term or for large pumping flows in relation to the productivity of the aquifer under consideration. This confirms the tendency for the system to be homogenized in the case of long periods or large pumping flows in relation to the productivity of the aquifer under consideration. However, despite the pumping of nearly 9 volumes of wells at this flow, the homogenization of the system is not complete and therefore the purge of the well is not complete. 1) shows that the transmissivity of the deep zone is greater than that of the top. However, this estimation of the mixture is dependent on the hypothesis of the existence of only 2 poles and the ability to identify the conductivity of each of these poles. It is therefore only an approximation. When mixing between multiple poles the uncertainty increases strongly and without fine qualification of the local permeabilities, it is not possible to obtain a single solution. However, this rapid calculation gives an order of magnitude of the distribution of the flows and permits a more thorough interpretation of the sampling carried out.

The experiment carried out on this well made it possible to illustrate the dynamics in the casing depending on the purging and sampling procedure in the case of stratified pollutions in a transmissive and heterogeneous aquifer. The samples and their representativeness are always dependent on the well and the objectives sought, however, different conclusions can be generalized:

 Regardless of the flow rate, if pumping takes place in the screened interval, the water in the casing is never stressed (if a low drawdown is maintained as recommended),  The use of a discrete interval sampler allows sampling at different level in the well,  A dynamic sampling is always a mixture of waters of different levels weighted by the transmissivity,  This weighting can be reduced by the use of hydraulic pakers if the filtering mass does not allow connections between the levels isolated by the pakers,  The proportions of mixing between the different levels vary according to the pump depth and the transmissivity ratios of the different parts of the aquifer intersected until the purge is complete,  The position of the pump no longer influences the sampling once the asymptotic regime is reached,  3-5 well volumes may be insufficient to reach the asymptotic regime.

Modelling

The numerical modeling presented in this section aims at illustrating and confirming the conclusions of the experiments carried out some case studies and on the other hand to have consequences of purging step for theoretical cases. These different cases are illustrated in Figure 6. /s, respectively the domain. The 2D mesh consists of 50 columns and 60 rows. Four simulations are made by changing the position of the pump: in the casing, then at the top, middle and bottom of the screened interval. Simulations are carried out in transient conditions for hydrodynamics and concentration transport over a period of 3 hours with pumping in the well of 0.5 m 3 /h using the MARTHE code [START_REF] Thiéry | Code de calcul MARTHE -Modélisation 3D des écoulements dans les hydrosystèmes -Notice d'utilisation de la version 7[END_REF].

For example for the case f) (similar to the on site test) Figure 7 shows the origin of the pumped water over time. The modeling confirms that for long pumping times, once the asymptotic regime is reached, the origin of the water is the same and the position of the pump no longer matters.

Conclusion and perspectives

In the case of a transmissive aquifer, the logs in ambient conditions seem to be the most representative of the structure of the pollution and are more instructive on this point than the dynamic samples.

The dynamic samples (after obtaining the asymptotic regime) are more representative than the logs of the masses transported in the aquifer due to their weighting by the transmissivities. These two last points show that it is necessary to choose the sampling method (dynamic or static) according to the objectives sought (assessment of the situation or evaluation of mass transport). Thus sample must be representative of the situation, but this representativeness depends on the objectives of the sampling. Sampling for water quality analysis is carried out to meet objectives in contexts that may be different and therefore may require adaptation of sampling procedures.

In the case of site monitoring, establishment of a local geochemical background, the samples may be based on a standard or improved procedure, preferring stabilization of the water level and physicochemical characteristics of water extracted than to a purge volume. In that case, what we are looking for is an evaluation of the risk when the groundwater will be pumped and used. In that case the chemistry of the water extracted is representative of the global chemistry of the transmissive zones of the aquifer.

In the context of pollution's characterization, the development of temperature or conductivity profiles of the water column and the consideration of local hydrodynamics may be necessary to define the most suited purging and sampling procedure. The objective is then to have information on areas where there is little transmissivity where pollutants are likely to be stored. Thus a single sampling method cannot meet all the objectives. Because local hydrogeology has a strong influence on groundwater sampling, a sample is always a more or less important mixture of water from different portions of the screen column; it may be advisable to apply different methods depending on the objective and limit the size of the screen slot to limit mixing of water while purging and sampling.

In some context it is shown that purge has no effect on water quality: aquifer of homogeneous quality or well where the water naturally renews itself or it renews due to industrial pumping.

The possibility of not purging before sampling could be accepted if previously demonstrated:  Either by comparing the concentrations measured before and after purging, and / or by the absence of variation on physicochemical measurements in water column before and after purging and sampling,  Or by hydrodynamic proofs: measurement of velocities with a flowmeter without stress (pumping or injection), followed by the disappearance of tracers (fluorescein type or salt) after its injection into the well..

More modelling and on site tests on different hydrogeological contexts will be conducted before to achieve the production of guidelines for different cases.

  On a fractured media site with an acid tar contamination due to coke production, measurements were carried on a well screened from 15 m to 54 m for 2 hydrogeological conditions:  In ambient conditions: o Electrical conductivity log under ambient conditions before purge, o Electrical conductivity log under ambient conditions after purge, o Discrete interval sampling,  In dynamic conditions o Pump sampling MP1 at 4 flow rates from 0

Figure 1 :

 1 Figure 1 : Conductivity (mS/cm) before and after purge : purge at 32m deep

Figure 3 :Figure 4 :

 34 Figure 3 : PAH concentrations for 2 purging rates at 3 depths

Figure 5 :

 5 Figure 5 : Sampling of water for pumping at 0,02; 0,15 and 1 m3 / h with two pump positions at 30 and 38 m and and electrical conductivity in ambiant condition and after purging steps

Figure 6 :

 6 Figure 6 : Simulated configurations -radial model

Figure 7 :

 7 Figure 7: Time evolution of the origin of the pumped water in the well expressed in percentage

Table 1 :

 1 Estimation of flows by evaluation of mixtures

			Mixing with 2 poles
	Pole	EC (MS/cm)	Proportion	EC of the mixture (mS/cm)	Measured value (mS/Cm)
	Pole n 1: upper part of the aquifer Pole n 1: upper part of the aquifer	9.5 23	0.18 0.82	20.57	20.5
	This estimation (Table

Acknowledgement

This project was possible thanks to the financial support of the MEEM (French Ministry of the Environment, Energy and the Sea) and of Ademe (French Environment and Energy management Agency).


Soafara Andrianome

Laurent Hugueville

René De Seze

Brahim Selmaoui
email: brahim.selmaoui@ineris.fr

Increasing levels of saliva alpha amylase in electrohypersensitive (EHS) patients

Keywords: IEI-EMF, EMF, alpha amylase, markers, symptoms, electromagnetic fields, environmental illness, saliva

Increasing levels of saliva alpha amylase

Introduction

A rapid development of technologies has increased the number of electromagnetic field (EMFs) sources. These include extremely low frequencies (ELF) and radiofrequencies (RF).

While the debate is still focused on any potential public health impact due to their EMF emissions, some individuals report having symptoms that they relate to EMF exposure [START_REF] Bergqvist | Possible health implication of subjective symptoms and electromagnetic fields[END_REF]. This sensitivity to EMF has been termed "electromagnetic hypersensitivity" (EHS). These non-specific symptoms include dermatological symptoms (redness, tingling and burning sensations), as well as neurasthenic and vegetative symptoms (fatigue, tiredness, concentration difficulties, dizziness, nausea, heart palpitation and digestive disturbances) [START_REF] Eltiti | Development and evaluation of the electromagnetic hypersensitivity questionnaire[END_REF]). These symptoms are not part of any recognized syndrome.

Therefore, and due to the fact that these symptoms are medically unexplained, the world health organization termed this sensitivity as Idiopathic Environmental Intolerance (IEI) attributed to electromagnetic fields (EMF) [START_REF] Mild | Electromagnetic hypersensitivity[END_REF].

Despite the lack of scientific basis to objectively link EHS symptoms to EMF exposure, these symptoms and associated suffering certainly remain real and can vary widely in term of severity.

Whatever the cause, EHS impacts on the quality of life of sufferers. So far, there are no clear diagnostic criteria, and no medical diagnosis has been validated. Therefore, physicians find themselves unarmed when facing patients.

The last decade has also been marked with rare involvement of EHS participants to research in this field. Many studies have suffered from this difficulty of recruitment, and EHS enrollment remains as particular challenge.

In view of the above and despite this challenge of volunteer recruitment, our work aimed to look for biochemical markers in a population suffering from electrohypersensitivity in the hope of bringing new elements that could help diagnosis. In this regard, some variables of the endocrine

J U S T A C C E P T E D

and immune systems and of sympathetic activity will be analyzed and compared to a non electrohypersensitive population (control). Our hypothesis is that electrohypersensitive patients might have a different biological profile than the control.

Methods

Participants

After approval from the regional French ethics committee "Comité de protection des personnes (CPP) Nord Ouest CHU Amiens, France (CPP, 2014/8)", 30 adults reporting EHS across France and 25 control individuals were recruited for this study. Individuals declaring being electrohypersensitive participated in the study. EHS patients (23 females and 7 males) and controls (21 females and 4 males) were aged between 22 and 63 (mean age: 47 ± 9) and 22 to 66 (mean age: 46 ± 10) respectively. Controls and the EHS group were matched for age, sex and body mass Index (BMI). They were recruited by advertising in self-help blogs and hospital listservs or by traditional word of mouth.

The EHS group was selected according to the inclusion criteria from participants who completed a questionnaire of hypersensitivity conducted by the research group before testing.

Inclusion criteria for the EHS group were an attribution of experienced symptom(s) to an electromagnetic source and an absence of acute psychiatric disorders, such as acute depressive or paranoid psychosis. Exclusion criteria for participants were the presence of chronic illness and the consumption of antibiotics, anti-inflammatory agents or any other supplements that may interfere with the biological assessments.

Participants were instructed to abstain from consuming alcohol and coffee for 24 h before and during each experimental session. They were also advised not to brush their teeth and not to exercise for at least an hour before sample collection.

J U S T A C C E P T E D

Before entering the study, all participants gave written, informed consent to participate in the study. Used methods were carried out in accordance with the approved guidelines. EHS participants were from different regions of France and were reimbursed for their travel fees, and all participants received financial compensation for their participation. From the questionnaire filled out by the participants, we extracted the duration of their sensitivity.

Experimental protocol

The first part of the experiment took place at the participant's home. Saliva samples were collected with salivettes© tubes and a cotton plug (Sarstedt Nümbrecht, Germany) previously sent to the participants. Instructions were given on how and when to collect the saliva. Samples were self-collected on the day prior the experiment: one sample before bedtime and another one in the morning, just after waking up.

On the day of the experiment, participants came to the laboratory around 09:00 and spent their waking time in the laboratory. Participants were seated in a comfortable armchair inside an electrically shielded room with very low background fields. EMF levels were measured by a personal exposimeter, Satimo EME Spy 140 model. Newspapers, magazines and books were provided. From 10h00 to 11h30 and from 14h00 to 16h30, saliva was collected every 30 min.

The last saliva collection (16:30) was performed immediately after recording the electrodermal activity in response to auditory stimuli (results not shown).

In parallel, total urine was collected during the night from 20:00 to wake up time as the night fraction and from wake-up time to 16h30 as the day-time fraction. More details of the study protocol are available in [START_REF] Andrianome | Disturbed sleep in individuals with Idiopathic environmental intolerance attributed to electromagnetic fields (IEI-EMF): Melatonin assessment as a biological marker[END_REF].

J U S T A C C E P T E D

Sample collection, immunoassay and enzymatic protocols

After rinsing their mouth with water thoroughly, participants collected saliva for over 3 min with a cotton plugged salivette. Salivettes were immediately kept in an icebox. After transportation to the laboratory, the tubes were centrifuged (2 min, 3000 rpm) and stored without preservatives at -20°C until analysis. Urine fractions were aliquoted in appropriate tubes and stored at -20°C until analysis.

The manufacturer's recommended protocols for saliva were used without modification for saliva cortisol, Immunoglobuline A (sIgA), C-reactive protein (CRP) (Salimetrics, State College, PA) and saliva alpha amylase (sAA) (IBL, Hamburg, Germany). For urine, protocols of neopterin (IBL, Hamburg, Germany), creatinine (Enzo Life Sciences, Framingdate, NY) and cortisol (Arbor assay, MI) were applied.

Analyses were performed on all 12-points for saliva collected for cortisol, IgA and sAA while sCRP analyses were performed only on two points: 10:00 and 16:30.

All samples were tested in duplicate, and the average of duplicates was used as the final value.

Statistical analysis

Data were analyzed with SPSS statistical software, IBM SPSS Statistics for Windows, Version 20.0 (Armonk, NY: IBM Corp), and plotted by using Graphpad Prism, Version 5 (San Diego, CA). Descriptive results of continuous variables were expressed in mean ± SD. We applied a two factors ANOVA for assessing time and group effects.

U Mann Whitney or Kruskall Wallis non-parametric tests were used to compare subgroups of EHS in terms of the median duration. Values were considered significantly different when p < 0.05. Correlation tests were performed by the Spearman correlation.

J U S T A C C E P T E D

Results

Results from EMF recording have showed that levels of EMF for radiofrequency were below the detection limit. The required number of EHS participants for our study was achieved with 30 individuals self-reporting EHS who fully agreed to participate to our study. This number can be considered as satisfactory given the difficulty of patient's enrollment.

Due to aberrant results for biochemical analysis from one EHS male, all data from this participant were not included in statistical analysis, leading to 29 EHS participants. The mean duration of EHS was 7.51 years (range: 1-19 years). Groups did not differ in terms of gender proportion (84 % of female participants in control group versus 79.31 % in the EHS group, p=0.498), BMI (22.80 ± 2.60 kg/m² in the control group, and 22.35 ± 3.33 kg/m² in the EHS group, p=0.268) and age (46 ± 10 years in the control group, and 47 ± 9 years in the EHS group, p=0.999). No significant difference in proportion of smoker were found (24 % in the control group, and 13.8 % in the EHS group, p=0.307). In the control group, 24 participants were having an activity (96%), and 24 participants in the EHS group (82.75%), (p=0.123).

Salivary biomarkers

The cortisol saliva concentrations in both groups were not significantly different. The cortisol secretion profiles were similar (p > 0.05) with the lowest concentration at bedtime (control: mean: 0.04 µg/dL; EHS mean: 0.06 µg/dL) and the highest at waking time (control mean: 0.33 µg/dL; EHS mean: 0.29 µg/dL) (Fig. 1). sAA activity was significantly different (p < 0.0001) between the control and EHS group. The EHS group showed a significantly higher level of sAA (Fig. 2a).

ANOVA analysis showed that the IgA level in saliva was not significantly different between both groups but showed a significant time variation (p < 0. 0001) (Fig. 2b). A peak of sIgA was reached in the control group at awakening (106.53 ± 44.04 µg/ml) with a significant difference when compared to all the other time points. In the EHS group, no difference was disclosed between awakening sIgA levels and those of the other points (p > 0.05).

Figure 2. Salivary concentration of alpha amylase (sAA) (a) and immunoglobuline A (IgA) (b) in the control (black line) and EHS group (grey line) during experimentation (from bedtime to 16.30). Data expressed as mean ± SEM -Significant difference between groups for sAA (p < 0.0001) -Significant effect of time for sAA (p = 0.0001) and IgA (p < 0.0001) sCRP levels showed a significant increase in the morning, when compared to the afternoon levels (p < 0.0001). However, no significant difference could be detected between the EHS group and the control group (Fig. 3).

Saliva correlation with selected clinical parameters and EHS subgroup analysis

For the EHS group, sAA levels were different between the male and female participants. Male participants expressed higher levels than female participants. sIgA and sCRP were negatively correlated with age at r = -0.135, p = 0.017 and r = -0.36, p = 0.007 respectively.

The mean of sIgA showed a tendency to be positively correlated with EHS duration (r = 0.141, p = 0.010); sAA was also slightly but positively correlated with EHS duration (r = 0.254, p < 0.001).

J U S T A C C E P T E D

According to the median of duration of EHS, which was seven years, there is no subgroup difference in sCRP and sCortisol concentrations (Table 1). sAA activity is significantly higher in patients reporting EHS ≥ 7 years for saliva collected at 10:00. Likewise, sIgA shows a difference between the control and EHS groups at bedtime.

Urinary biomarkers

Urinary cortisol and neopterin concentrations, adjusted to creatinine, did not show any significant difference between the control and EHS groups. Urine neopterin to creatinine concentrations were significantly higher for the morning fraction than for the night fraction for both groups (p < 0.05) (Fig. 4).

Urine biomarkers and EHS subgroup analysis

According to the median duration of EHS, which is seven years, there is no difference in urinary cortisol levels between participants being EHS less than seven years and those being EHS more than seven years (Table 2). However, in the morning fraction of urine, a significant difference was found for neopterin between the control and EHS group and between the EHS subgroups.

J U S T A C C E P T E D

Discussion

This study was designed in an attempt to find out some biological markers that may differentiate EHS individuals from the rest of the population (non-EHS) and explain the occurrence of the atypical symptoms they complain about. In this regard, some markers of the immune, circadian and autonomic nervous systems were explored.

Among the markers analyzed, cortisol is considered as a reliable marker of the endogenous circadian oscillator [START_REF] Haus | Chronobiology in laboratory medicine[END_REF]) and as a primary biological indicator for stress reactions for the hypothalamic-pituitary-adrenal (HPA) axis [START_REF] Kirschbaum | Salivary cortisol in psychobiological research: an overview[END_REF][START_REF] U S T A C C E P T E D Kirschbaum | Salivary cortisol in psychoneuroendocrine research: recent developments and applications[END_REF]. Our present work showed no significant differences in saliva cortisol levels between a group suffering from EHS and a control group. In light of the above, our results are in agreement with the literature and show a peak in the morning immediately after wake up and a trough around bedtime [START_REF] Krieger | Rhythms of ACTH and corticosteroid secretion in health and disease, and their experimental modification[END_REF][START_REF] Selmaoui | Reproducibility of the circadian rhythms of serum cortisol and melatonin in healthy subjects: a study of three different 24-h cycles over six weeks[END_REF]. Neither the cortisol peak in the morning nor the lowest values of the day (bedtime) were different between the two groups. This may suggest that the HPA axis is not affected in EHS patients.

Additionally, the amount of cortisol in urine did not differ between both groups (p > 0.05), confirming the finding of the saliva analysis.

Besides, the salivary enzyme α-amylase has been proposed as a marker for the stress-induced activity of the sympathetic nervous system (SNS), or more precisely of the sympathetic adrenomedullary (SAM) system, which is also involved in the secretion of adrenaline and noradrenaline [START_REF] Chatterton Rt | Salivary alpha-amylase as a measure of endogenous adrenergic activity[END_REF][START_REF] Rohleder | Determinants of salivary alpha-amylase in humans and methodological considerations[END_REF] . Unlike cortisol, saliva alpha amylase levels were found to be significantly higher in the EHS group, in comparison to the control group. The question that arises here is whether this change in α-amylase activity was a result of outside exposure to EMF or simply a consequence of chronic stress.

J U S T A C C E P T E D

The question whether EHS patients are stressed by various environmental/ biological factors other than EMF is at the center of various research. Nowadays, there are only few studies dealing with the EHS population and which are not so far conclusive. Some few studies have tried to find a link with other diseases or other environmental factors such as study by [START_REF] Palmquist | Overlap in prevalence between various types of environmental intolerance[END_REF] who found an overlap between 4 types of environmental intolerances: exposure to odorous/pungent chemicals, certain buildings, everyday sounds and EMFs. While some other authors [START_REF] Rubin | Idiopathic environmental intolerance attributed to electromagnetic fields (formerly "electromagnetic hypersensitivity"): An updated systematic review of provocation studies[END_REF], proposed the presence of a "nocebo effects". In our population: no specific questionnaire about environmental stressing factors were presented. On the other hand, levels of cortisol do not indicate the presence of stress.

With regard to the biological stress system, one may wonder why only SAM has been activated

and not the HPA axis if there are stress-related changes. It is known that these two physiological systems are independent in terms of neuronal circuits (Ulrich-Lai & Herman 2009) and that they may react differently to stressors. Some results suggested that the responses of the HPA and SAM systems differ depending on the type of stressor [START_REF] Maruyama | Differences in salivary alpha-amylase and cortisol responsiveness following exposure to electrical stimulation versus the Trier Social Stress Tests[END_REF]. Moreover, salivary cortisol levels increase a few minutes after the initiation of stress [START_REF] Droste | Corticosterone levels in the brain show a distinct ultradian rhythm but a delayed response to forced swim stress[END_REF][START_REF] Takai | Effect of psychological stress on the salivary cortisol and amylase levels in healthy young adults[END_REF], meaning that the HPA axis is relatively slower to respond while sAA levels are characterized by a rapid increase after a stress action, followed by a rapid decrease to normal levels after the cessation of the stressor [START_REF] Speirs | The influence of sympathetic activity and isoprenaline on the secretion of amylase from the human parotid gland[END_REF][START_REF] Takai | Effect of psychological stress on the salivary cortisol and amylase levels in healthy young adults[END_REF]). However, in the present study, sAA levels were shown to be constantly elevated in EHS individuals, even during their stay (from 09h00 to 16h30) in the EMF shielded room in our laboratory. The role of EMF on the increase of AA is still to be proven. However, our results indicated that EHS population seems to have higher levels of sAA than control group. And the cause of this increase is to be yet determined.

Moreover, a significant and positive correlation was shown between sAA levels and the duration of EHS. Indeed, the longer the duration of EHS, the higher the level of sAA. This may suggest

J U S T A C C E P T E D

that the longer they are EHS, the more they react. One may speculate that this could be due to a degraded system in patients presenting EHS for a long-time.

Besides, a study has shown that exposure to radiofrequency electromagnetic fields (RF-EMF) emitted by mobile phone base stations affected salivary α-amylase in healthy participants [START_REF] Augner | Effects of exposure to GSM mobile phone base station signals on salivary cortisol, alpha-amylase, and immunoglobulin A[END_REF]. In addition, a study of people living next to phone base stations, selfdeclaring base-station neighbors (distance to base station ≤ 100 meters), had higher levels of αamylase accompanied by a higher general strain [START_REF] Augner | Are people living next to mobile phone base stations more strained? Relationship of health concerns, self-estimated distance to base station, and psychological parameters[END_REF][START_REF] Augner | Effects of exposure to GSM mobile phone base station signals on salivary cortisol, alpha-amylase, and immunoglobulin A[END_REF].

This may suggest that people self-reporting EHS could share similarities with this population, but whether this is due to actual exposure to EMF or to other underlying factors is not clear. In our study, unfortunately, recordings of exposure field density before entering the study (outside the experimental lab) were not available.

Moreover, levels of IgA in saliva vary in response to physical and psychological stress through interactions with the autonomic nervous system [START_REF] Tsujita | Secretory IgA in saliva can be a useful stress marker[END_REF][START_REF] Bishop | Acute and chronic effects of exercise on markers of mucosal immunity[END_REF]. Saliva IgA concentrations were found to be similar between the EHS population and the matched control group. However, we could observe a significant variation in the pattern of IgA between morning and afternoon. This is in accordance with previous studies showing that IgA patterns follow a circadian rhythm [START_REF] Casale | Circadian rhythm of immunoglobulins in aged persons[END_REF]. Subgroup analysis showed that people declaring longer EHS had higher levels of IgA for bedtime collected saliva.

C reactive protein (CRP) is widely used as an inflammatory marker and could indicate bio inflammation of the body. High levels of serum CRP have been correlated to cardiovascular risks and other inflammatory or immune diseases [START_REF] Tracy | Relationship of C-reactive protein to risk of cardiovascular disease in the elderly. Results from the Cardiovascular Health Study and the Rural Health Promotion Project[END_REF][START_REF] Clos | C-reactive protein as a regulator of autoimmunity and inflammation[END_REF]. Recently, elevated levels of CRP and hsp70 have been detected in frequent users of mobile phones [START_REF] Balakrishnan | Hsp70 is an independent stress marker among frequent users of mobile phones[END_REF]. In this cohort, concentrations of saliva CRP in the EHS group were

J U S T A C C E P T E D
not different from those of the control group, indicating the absence of an acute inflammatory process. Furthermore, highly sensitive CRP was higher in the morning than in the evening. This finding is compatible with results in healthy young adults [START_REF] Izawa | The diurnal patterns of salivary interleukin-6 and C-reactive protein in healthy young adults[END_REF].

Neopterin in urine showed a light-dark variation in both groups, which corroborates previous findings [START_REF] Garcia-Gonzalez | Diurnal variations in serum neopterin levels are associated with the pineal hormone melatonin circadian rhythm in healthy human subjects[END_REF]. Surprisingly, subgroups analysis showed a variation according to the duration of the EHS effect. We found, in the day fraction of urine, a lower level of neopterin in patients reporting EHS for more than seven years, compared to those with shorter EHS. Little information is available on immune markers in EHS and on the immune effect of electromagnetic fields. Immune cells or cytokines were not modified under acute exposure to ELF magnetic fields [START_REF] Selmaoui | Acute exposure to 50 Hz magnetic field does not affect hematologic or immunologic functions in healthy young men: a circadian study[END_REF][START_REF] Selmaoui | Acute exposure to 50-Hz magnetic fields increases interleukin-6 in young healthy men[END_REF] or GSM 900 MHz exposure [START_REF] Tuschl | In vitro effects of GSM modulated radiofrequency fields on human immune cells[END_REF]. In an experiment with exposure to a GSM signal, neopterin was not modified in saliva [START_REF] Radon | No effects of pulsed radio frequency electromagnetic fields on melatonin, cortisol, and selected markers of the immune system in man[END_REF]. Other studies showed that acute exposure to ELF did not have an effect on blood concentrations of TSH, FSH, LH, testosterone, cortisol or neopterin [START_REF] Gamberale | Acute effects of ELF electromagnetic fields: a field study of linesmen working with 400 kV power lines[END_REF][START_REF] Selmaoui | Acute exposure to 50 Hz magnetic field does not affect hematologic or immunologic functions in healthy young men: a circadian study[END_REF].

In our study, we used low tech and bloodless methods to analyze markers in the EHS population.

Besides an accepted general method of collection (a cotton swab), all the stress-minimizing methods used here are shown to be appropriate in quantifying such markers.

Regarding the age of the EHS group, variability due to the age widths of participants was reduced by recruiting matched volunteers as far as possible. Besides this age characteristic, for a more powerful analysis and according to the subgroup exploration, a larger number of participants would be required for future studies.

Here, the classification of EHS subgroups according to the duration of their EHS can be a disputable point due to the declarative nature of this information. It could be understood that the A retrospective analysis of exposure and recording the lifetime information of participants could be an important issue in understanding EHS that could be coupled with such data. Until now, limited information about the mechanism of EHS is available, and few studies have explored the biological mechanism in EHS, though particular emphasis has been put on "provocation studies," which were not conclusive. In the absence of a known biological mechanism, the exploration of more markers involved in immunological activity could complete our study, without minimizing the complex interaction among symptoms, health concerns and risk perceptions in this particular population.

Conclusions

In summary, an exploration of selected and measurable markers of EHS was conducted in our study to understand this condition. Our finding indicates that markers of the circadian rhythm and immune system were not affected in EHS compared to non EHS individuals. Our results suggest different profiles of EHS according to the duration of the syndrome. Furthermore, the mechanism underlying the higher expression of α-amylase should be focused on. These findings could be new elements in the understanding, explanation and comprehension of the physiopathology of EHS.

Figure 1 .

 1 Figure 1. Saliva level of cortisol in control (black line) and EHS (grey line) groups, data

Figure 3 .

 3 Figure 3. Mean ± SEM of salivary CRP in the control (black bars) and EHS group (grey bars)

Figure 4 .

 4 Figure 4. Mean ± SEM of urine cortisol:creatinine ratio (a) and neopterine:creatinine ratio (b) in

  exact timing of the triggering the syndrome may be an estimate without any known diagnostic. It could correspond here to the development of major symptoms or the occurrence of a special event. Here again, the lack of an existing diagnosis of EHS poses a serious problem.

Figure 2 .Figure 3 .Figure 4 .

 234 Figure captions Figure 1. Saliva level of cortisol in control (black line) and EHS (grey line) groups, data represent mean ± SEM. Significant time effect (p < 0.05) with typical morning peak (waking point). No significant effect of group was observed (p > 0.05) Figure 2. Salivary concentration of α-amylase (a) and IgA (b) in the control (black line) and EHS group (grey line) during experimentation (from bedtime to 16.30). Data expressed as mean ± SEM -Significant difference between groups for α-amylase (p < 0.0001)

Table 1 .

 1 Demographic characteristics and levels of markers in control and subgroups of EHS based on EHS duration (median duration = 7 years) in saliva collected at night, wake up, 10:00, and 16:30 *sCRP only analyzed on 10:00 and 16:30 points; Data are quoted as mean ± SD value for

	continuous variables, or number;
	symbols: ‡ p < 0.05; § p < 0.01
		Control group	EHS	EHS	P-	Global P*
			duration <7	duration ≥ 7	values	values
			years (A)	years (B)	(A vs.
					B)
	n	25	17	12	-	-
	Mle	4	3	3	-	-
	Female	21	15	8	0.494 -
	Age, years	46.24	47.50	44.82	0.674 0.942
	BMI, kg/m²	22.80 ± 2.60	22.11 ± 2.54	21.91 ± 3.52	0.707 0.470
	sCRP
	*10:00 ng/mL	2.74 ± 0.90	2.56 ± 1.11	3.09 ± 0.88	0.173 0.419
	*16:30 ng/mL cortisol bedtime µg/dL awakening µg/dL 0.331 ± 0.178 1.75 ± 0.90 0.041 ± 0.041 10:00 µg/dL 0.187 ± 0.095 16:30 µg/dL 0.074 ± 0.058 sAA bedtime U/mL 161.26 ± 124.07 182.14 ± 125.07 240.10 ± 108.44 0.166 0.124 1.90 ± 0.73 2.10 ± 0.90 0.674 0.489 0.061 ± 0.047 0.049 ± 0.028 0.777 0.333 0.279 ± 0.194 0.299 ± 0.250 0.698 0.049 ‡ 0.152 ± 0.079 0.186 ± 0.074 0.227 0.298 0.095 ± 0.075 0.092 ± 0.052 0.744 0.584 awakening U/mL 75.12 ± 60.45 146.29 ± 140.12 146.82 ± 62.68 0.501 0.014 ‡ 10:00 U/mL 121.62 ± 129.84 185.67 ± 109.54 337.50 ± 163.42 0.043 ‡ 0.007 § 16:30 U/mL 119.20 ± 152.86 208.79 ± 142.71 226 ± 89.71 0.781 0.033 ‡ IgA bedtime µg/mL 47.74 ± 30.32 30.41 ± 37.86 57.46 ± 39.04 0.03 ‡ 0.034 ‡ awakening µg/mL 106.54 ± 44.04 55.70 ± 89.12 64.29 ± 82.56 0.913 0.003 § 10:00 µg/mL 26.20 ± 31.06 53.57 ± 33.13 33.38 ± 34.90 0.073 0.025 ‡ A C C E P T E D
	16:30 µg/mL J U 92.27 ± 64.62 S T	52.74 ± 36.67	77.90 ± 96.56	0.781 0.061

Table 2 .

 2 Biochemical urinary characteristics of the study population: controls (n =25) and the two subgroups (based on the median duration of EHS), data are quoted as mean ± SD

		Control	EHS duration < 7	EHS duration ≥ 7	P (A	Global P
		(n=25)	years (A) (n=17)	years (B) (n=12)	vs. B)
	Night fraction (20:00 to awakening)
	Cortisol:creatinin	54.90 ±	38.85 ± 47.63	34.18 ± 33.33	0.914 0.989
	(nmol/mmol)	12.14
	Neopterin:creatinin	223.82 ±	229.40 ± 322.17	107.15 ± 81.53	0.377 0.153
	(µmol/mol)	314.10
	Day fraction (awakening to 16:30)
	Cortisol:creatinin	33.51 ±	57.10 ± 10.33	52.47 ± 53.71	0.603 0.368
	(nmol/mmol)	33.64
	Neopterin:creatinin	354.82 ±	472.58 ± 438.35	138.68 ± 133.43	0.011 ‡ 0.037 ‡
	(µmol/mol)	519.18
	‡ p < 0.05; § p < 0.01

Acknowledgements

This project was financed by a grant from the French National Research Program for Environmental and Occupational Health of ANSES (Project n° 2013-2-02). Particular thanks to

J U S T A C C E P T E D

all volunteers for their participation. Thanks to Kelly Blazy and Maïté Hanot-Roy for their contribution on elisa analysis.

Declaration of interest statement

The authors report no declarations of interest.

J U S T A C C E P T E D


Sara Abada

Martin Petit

Amandine Lecocq

Guy Marlair
email: guy.marlair@ineris.fr

Valérie Sauvant-Moynot

François Huet

Combined experimental and modeling approaches of the thermal runaway of fresh and aged lithium-ion batteries

Keywords: Li-ion battery modeling, safety issues, thermal runaway, multiphysics modeling, battery ageing, electrochemical storage

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

Introduction

Lithium-ion batteries (LIB) are currently considered as the best available solution to store energy on board of new generation of electrical, plug-in, or hybrid electrical vehicles, or to develop stationary applications with alternative green energy sources (solar, wind, photovoltaic). However, their ageing, high cost, and safety aspects still represent critical issues for these storage systems.

To assess and mitigate these drawbacks, major research efforts have been carried out in the recent past. Thermal runaway has been identified as a major concern with LIBs on the full value chain, potentially leading to uncontrolled fires and explosions causing the failure of the entire battery pack [1]. Catastrophic hazardous events reported to originate from LIBs, albeit not so frequent, such as violent venting, smoke, fire, explosion [START_REF] Doughty | Interface, the Electrochemical Society[END_REF] still act as a clear restraint to LIB market development for high power/high energy applications, and still even in rare cases for consumer market appliances [START_REF] Duh | [END_REF]4].

The thermal runaway of LIBs has been described as a multistep process implying a series of exothermic reactions related to the decomposition of different components of the cell (electrodes, electrolyte, and binder), together with their interactions occurring above a critical temperature, which led to an increased temperature inside the cell [5].

The most consensual analysis of the thermal runaway of a LIB considers that the phenomenon proceeds in four main stages occurring successively according to the onset of exothermic decomposition reactions as follows:

1. The first degradation reaction leading to the initial subsequent increase in cell temperature is usually reported to concern the solid electrolyte interface (SEI) that forms a passivation film on the surface of the carbon negative electrode. The onset of this reaction is ranging between 90°C and 120°C depending on the chemistry of the cell [6].

2. When the SEI starts to decompose, the electrolyte reacts violently with the active material of the negative electrode, which accelerates the cell temperature rise [START_REF] Guo | Safety: thermal runaway in[END_REF].

3. At temperatures above 120°C [START_REF] Van Schalkwijk | Advances in Lithium-Ion Batteries[END_REF], the positive electrode materials start to decompose, leading to evolution of oxygen that can react with the electrolyte and generate additional heat According to thermal stability studies of commercial LIBs, the onset temperatures of the exothermic reactions at positive electrode materials rank in the following order: LiNiO2 (~180°C) < LiCoO2 (~200°C) < LiMn2O4 (~220°C) < LiFePO4 (> 240°C)

LiFePO4 is then considered as the most thermally stable material of positive electrodes [START_REF] Guo | Safety: thermal runaway in[END_REF]. 4. The electrolyte, most commonly a lithium salt dissolved in a mixture of organic carbonates, can decompose exothermically at elevated temperatures (> 200°C) depending on its composition (salt, solvents, additives) [START_REF] Binotto | Review Thermal Runaway Reactions Mechanisms, HELIOS (High Energy Lithium-Ion Storage Solutions) Project[END_REF].

At the cell level, the thermal runaway features globally depend on the chemistry of the cell components, as well as the shape and design of the cell [START_REF] Hallaj | [END_REF]. They also depend on the state of charge (SOC) of the cell. Cells of higher SOC have revealed a lower thermal stability in dedicated studies, meaning a lower onset temperature of the thermal runaway reactions [5] and a higher selfsustained rate of increase of the cell temperature [START_REF] Del Corso | High energy lithium ion storage solutions: comparative assessment of 4 chemistries of cathode for EV and PHEV applications[END_REF]. Some recent studies have investigated the effect of ageing on the thermal stability of cells by performing thermal stability tests under abuse conditions (overcharge, overheating…) to determine correlations between battery ageing and safety [START_REF] Fleischhammer | [END_REF][13][14][15]. For example, the influence of cyclic ageing on the thermal behavior of Li-ion cells with LixNi1/3Mn1/3Co1/3O2/LiyMn2O4-blend electrode positive materials was studied by Fleischhammer et al. [START_REF] Fleischhammer | [END_REF]. They showed that cells aged by high-rate cycling had about the same safety behavior as fresh cells. In contrast, a significant increase of self-heating was observed for cells showing plating of metallic lithium on the negative electrode due to low temperature charging. On the other hand, Röder et al. [13] studied the influence of calendar ageing (at full SOC and 60°C) on the thermal stability of a commercial Li-ion cell containing a mixture of LiMn2O4 and Li(Ni1/3Mn1/3Co1/3)O2 as positive electrode material and graphite as negative electrode. They showed that the self-heating onset temperature of the aged cells was much lower than that of the fresh cells. These experimental studies show that the ageing history of aged cells plays a strong role on their safety so that an open issue is the way battery ageing affects its sensitivity to thermal runaway, in particular for the graphite-LiFePO4 (LFP/C) technology claimed to be safer than the others. So far, battery safety studies have mainly consisted in experimental approaches based on safety and abuse tolerance tests as well as post-mortem battery analyses [16]. However, these approaches show significant shortcomings due to the complexity of the phenomena involved in battery operation. Moreover, the safety and abuse tolerance tests are destructive so that the cost in battery test samples becomes an issue. Modeling and simulation are then used to predict battery performance and safety and help the design of new batteries; especially as experimental databases are now available for model calibration and validation, at least at cell level.

The number of investigations on LIB safety has increased in the last decade, including safety studies implementing modeling approaches [17][18][19]. However, battery modeling in abuse conditions was exclusively aimed for some time at modeling single events leading to battery failures as short circuit, mechanical abuse, overcharge or fire propagation. Combined modeling of electrical, chemical, and thermal behaviors of LIBs in abuse operating conditions emerged only after 2012. NREL (National Refreshable Energy Laboratory), Battery design LLC, and Sandia National Laboratories works constitute a large step forward in battery modeling and simulation in abuse conditions at different scales (cell/module/pack). However, all models of thermal runaway in abuse conditions available at the cell level are presently based on the chemical degradation reactions considered in Kim et al.'s model [6], none of them taking into account the influence of cell ageing. In this context, the objective of the present work is the development and validation of accurate physical models of the thermal behavior of fresh and aged LIBs in conditions leading to thermal runaway, for further practical applications to sizing and design support of safer batteries.

A 3D thermal runaway model including ageing phenomena was developed and calibrated for cylindrical 26650 LFP/C cells (A123 Systems) using a battery test calorimeter (BTC) before validation with oven tests. The thermal safety of such LFP/C cells (fresh and aged at two different residual capacity) was then experimentally determined for comparison with expectations from the model to find out the influence of calendar ageing on battery safety.

Materials and methods

In order to build a multi-scale predictive model of thermal runaway adapted to fresh and aged LIBs, two complementary studies were carried out in parallel, a modeling work to develop several sub-models and achieve a coupled multi-physics and multidimensional model, and an experimental work to initially calibrate the sub-models and validate the model. As indicated above, A123 (2.3 Ah) cylindrical LFP/C cells were used. Their characteristics are presented in Table 1.

All the tests were carried out with fully charged cells (100% SOC), which is known to be the worst-case scenario [5].

Modeling approach

As mentioned above, thermal runaway models available at the cell scale are all based on the chemical degradation reactions introduced in Kim et al.'s model [6] devoted to lithium-cobaltoxide/carbon (LCO/C) technology. Reaction parameters were determined from accelerated rate calorimetry and differential scanning calorimetry tests, mainly performed by Hatchard and coworkers [20,21]. No effort has been made to validate Kim et al's model for the LFP/C technology and, to our knowledge, no degradation reaction parameters for this technology are available in the literature. Moreover, battery ageing phenomena have never been taken into account in thermal abuse models and their influence on thermal runaway has not been verified by simulations yet. To remedy to these shortcomings, the present model, based on Kim et al's work, is structured as an interactive assembly of 3 sub-models, as shown in Fig. 1.

Sub-model 1: 3D cell degradation chemical reactions

This sub-model considers the physicochemical mechanisms leading to thermal runaway to predict the chemical heat released by the exothermic reactions related to the degradation of the various cell components in extreme operating conditions, using Equation ( 1):

e pe ne abuse SEI Q Q Q Q Q          (1)
where the volumetric heat generation rates of the degradation reactions are denoted SEI Q  for the SEI, ne Q  for the negative electrode, pe Q  for the positive electrode, and e Q  for the electrolyte. Each volumetric heat generation rate is calculated according to the following relation:

i i i i R W H Q   ( 2
)
where Hi (J g -1 ) is the specific heat release, Wi is the specific active material content (g m -3 ) and Ri (s -1 ) is the reaction rate in cell component i. The reaction rates obey the following Arrhenius laws (Equations ( 3) to ( 6)) [6,22]:

SEI SEI SEI a, SEI SEI SEI ] exp[- = d d - = m c RT E A t c R ( 3
)
where cSEI is the normalized concentration of metastable species containing Li in the SEI layer, ASEI (s -1 ) is the frequency factor of the SEI decomposition reaction, Ea,SEI (J mol -1 ) its thermal activation energy, and mSEI the order of the reaction.

] exp[- ] exp[- d d - = d d = ne a, ne ref SEI, SEI ne ne SEI ne ne RT E c t t A t c t t R m  (4)
in which the increase in normalized thickness of the SEI, tSEI, is expressed as a function of the normalized concentration, cne, of Li intercalated in the negative electrode, the frequency factor, Ane (s -1 ), of the reaction, the activation energy, Ea,ne (J mol -1 ), and the order, mne, of the reaction. The SEI normalized thickness is defined as tSEI,0 (SEI/SEI,0), where SEI is the SEI thickness and the subscript 0 indicates the initial value. The tSEI,ref parameter allows the diffusion limitation of the SEI growth to be taken into account. It is sometimes considered equal to the initial dimensionless thickness of the SEI layer [23] but other authors use it as an independent tunable parameter [6]. In the present work, this parameter was set to 1 and the other parameters of the negative electrode degradation reaction were fitted consequently.

] exp[- ) 1 ( = d d = pe a, pe pe r pe p pe RT E A t R m m     ( 5
)
where  is the degree of conversion of the active material of the positive electrode, Ape (s -1 )

and Ea,pe (J mol -1 ) are the reaction kinetic parameters and mpe is the reaction order.

= d d - = m c RT E A t c R ( 6
)
where ce is the dimensionless concentration of electrolyte, Ae (s -1 ) and Ea,e (J/mol) are the kinetic parameters of the electrolyte decomposition reaction and me is the reaction order.

The kinetic parameters of the degradation reactions were extracted both from literature data and the experimental test carried out in this work.

Sub-model 2: 0D calendar ageing

Multiple and complex ageing mechanisms have been described for various chemistries of positive and negative electrodes. The mechanism considered in the present model reflects the growth of the SEI at the negative electrode since it has been identified as the main cause of calendar ageing of graphite based LIBs [24]. In this model, in which the thickness of the SEI layer is assumed uniform inside the cell (0D model), a static ageing approach has been followed to evaluate the thickness of the SEI layer corresponding to a specific cell state of health. This approach does not take into account the ageing conditions of the cell and is based on the assumption that the capacity loss measured was entirely due to the SEI growth.

The SEI growth mechanism leading to calendar ageing is limited to the solvent (S) reduction at the surface of the negative electrode active material [24], according to the following reaction:

P Li 2 e 2 S      (7)
The amount of SEI formed can be expressed using Faraday's law as a function of the cell capacity loss C in Coulomb:

F C n 2 1 SEI  (8)
so that the volume change of SEI can be derived from the change in capacity loss:

F C M V 2 1 d d SEI SEI SEI   (9)
where MSEI and SEI are the molar mass and density of the film, respectively. To estimate the increase in SEI thickness, dSEI, it is assumed that the radius, Rs,n, of the graphite particle modeling the negative electrode is larger than the thickness of the formed SEI [24]. In that case:

F C S M n 2 1 d d SEI SEI SEI    (10)
where the electroactive surface of the negative electrode, Sn, is given by the following equation:

n s, n n s, n 3 R A S    (11)
in which A is the geometric area of the negative electrode (assumed to be flat), n its thickness and s,n the volumic fraction of carbon [24].

Integration of Equation ( 10) between times 0 (initial state) and t (calendar ageing state)

gives the relationship between the thickness of the SEI layer, SEI, and the cell capacity loss, Closs, induced by calendar ageing:

F C S M t n 2 1 0 , ) ( loss SEI SEI SEI SEI      (12)

Sub-model 3: 3D-thermal model

This sub-model is based on the energy balance equation involving the different heat sources in the core of the battery cell:

2 2 z 2 2 r abuse p z T k r T k Q t T C            ( 13
)
where abuse Q  is provided by the sub-model of cell degradation chemical reactions (Equation ( 1)),

 is the density of the cell, Cp its mean thermal capacity, T its temperature, and kr, kz are the thermal conductivities in the radial and vertical directions of the cylindrical cell, respectively. The heat exchange between the cell and its external environment is integrated in the boundary conditions expressed by Equation ( 14):

) ( ) ( 4 ext 4 skin ext skin conv exchange n T T T T h Q n T k            (14)
where hconv, , , Tskin, and Text are the convective heat transfer coefficient, the emissivity of the cell surface, Stefan-Boltzmann's constant, and the skin and ambient temperatures, respectively.

The left-hand side term represents the conductive heat flux from the cell core in the direction n (n = r, z) while the right-hand side term is the sum of the convective and radiative heat transfer rates.

For simplification, this sum was written as heq (Tskin -Text) in this study, where heq is an equivalent heat transfer coefficient. With such assumption, the thermal exchanges were correctly described for the battery cell but this would not be the case for a battery module or pack.

In order to calibrate the thermal sub-model, the thermal properties of the cell under investigation (heat capacity Cp in J kg -1 K -1 , thermal anisotropic conductivity k in W m -1 K -1 , and heat transfer coefficient h in W m -2 K -1 ) were estimated by fitting experimental data obtained from thermal characterization tests in a BTC.

The sub-models described in sections 2.1 have been ultimately coupled to achieve a multiphysics predictive model of the thermal runaway of fresh and aged Li-ion cells. This model was implemented on Comsol Mutiphysics platform with fit-to-purpose experimental calibration and validation procedures for the cell under study, as described in the following section.

Experimental study

Several experimental tests were performed to calibrate and validate the thermal runaway model and sub-models, as indicated in Table 2.

Calendar ageing characterization

The initial state of the cells under investigation was characterized by their capacity and internal resistance before accelerated ageing. The actual capacity of each cell was determined through two charge/discharge cycles. The charge was performed successively at a constant current rate of 1 C until a 3.6 V cell voltage, and at a constant voltage of 3.6 V until the current was lower (resp. 30% ± 3%).

Experimental calibration of the degradation chemical reactions sub-model

The calibration of the degradation chemical reactions parameters was performed using a BTC to ensure pseudo-adiabatic experimental conditions. The BTC (HEL Company, Hertfordshire, UK) has a cylindrical test chamber (diameter 35 cm and height 32.5 cm) heated on all sides (bottom, top and lateral). The tested cell was positioned in the center of the BTC enclosure. Two type-K thermocouples were placed on the negative and positive electrode sides of the cell to monitor the cell surface temperature during the test. The cell was then covered with an aluminum foil to ensure good heat diffusion and finally surrounded by a heater wire of length 1 m and resistance 11 . The energy for heating the cell was provided by the heater wire while the energy for heating the climatic chamber was provided by the BTC, both heating procedures being controlled by the program operating the BTC. The cell voltage was also monitored during the test.

To determine the onset of self-heating, the BTC was operated with a MultiRamp program similar to the well-known HWS (heat, wait, search) mode. In this study, fully charged cells were first heated up to a temperature of 30°C before a calibration step of 90 min. Then, the temperature was increased in discrete steps of 5°C, each of them (under safe behavior of the cell) followed by a 30 min long wait step during which the cell was held at constant temperature to reach complete thermal equilibrium with the calorimetric system. If the cell temperature did not increase at a rate above a threshold value, typically 0.02°C /min, indicating no adverse thermal event, the temperature was increased by another heating step and the algorithm was repeated as needed until some exothermal reaction occurred. By contrast, if the cell temperature increased at a rate above the mentioned threshold value during the wait step, the BTC switched to the exothermic mode during which its temperature closely matched the cell temperature, thus maintaining the adiabatic state. This was followed by search step of 10 min during which the system tried to detect any selfheating phenomena on the cell surface. The experiments in the BTC were terminated by cooling the cell once it reached its upper temperature limit of 350°C. Thanks to the adiabatic conditions and the stepwise process allowing the observation of individual self-heating reactions, this procedure allowed the reaction rate parameters of the different exothermic reactions degrading the cell components during a thermal runaway to be determined.

The use of oven tests is very frequently reported in the literature, particularly to apply thermal abuse conditions on the investigated cells in a well-controlled environment, which allows experimental validation of thermal runaway models [6,21,22,[25][26][27]. In such tests, cells are heated to excessive temperatures in order to examine their thermal stability performance. Being closer to actual environmental conditions (no adiabaticity), oven tests were used in this work to validate the thermal runaway model. The tested cell was placed in the center of a thermostatically controlled oven (France Etuves, model XL 1000) and the air was circulated at full speed to obtain a uniform temperature around the cell. Up to five type-K thermocouples were fixed on the cell surface (two on the positive and negative electrodes and three on its lateral surface) and four additional thermocouples were installed to measure the oven temperature. The oven was first preheated to the target cell temperature of 30°C that was maintained during 60 min. Thereafter, the cell temperature was increased at a rate of 5°C / min until the activation of thermal runaway or the observation of a plateau value of 200°C stabilized during two hours as a proof of a complete reaction of the cell. The surface temperature of the cell was measured at a sample rate of 1 Hz with a data logger.

Influence of ageing on the cell degradation chemical reactions

As an example, Fig. 2 shows the evolution of the surface temperature of a fresh A123 cell and its voltage during the test in the BTC. Three critical temperatures can be identified [15]:

 the initial self-heating temperature, T1,, related to the SEI layer decomposition reaction: it is important to note that its detection depends on the accuracy and sensitivity of the testing device and also on the threshold detection value indicated in the BTC operating program (0.02°C / min for exothermic reaction in this work).

 the starting temperature, T2, of the melting process of the battery separator material: this temperature was detected by the endothermic reaction indicated by the hardly visible decrease of the recorded temperature at time 1,077 min. It is known that the melting point of commercial separator materials is expected to range between 130°C and 160°C depending on their composition. Above this critical temperature the separator loses its mechanical integrity and starts to induce internal -micro-short-circuits.

 the thermal runaway trigger temperature, T3: the accelerated increase in temperature was caused by the heat released by the exothermic reactions and by the occurrence of short circuits, as shown by the steep decrease of the cell voltage.

Fig. 2 shows five stages in the evolution of the cell temperature: (1) end of safe state, [START_REF] Doughty | Interface, the Electrochemical Society[END_REF] and ( 3) still reversible self-heating regions, (4) triggering of thermal runaway, and (5) cell cooling down. However, the peaks detected do not allow a clear identification of the corresponding degradation reactions of the cell components.

Several such experiments were performed in the BTC with fresh and aged A123 cells to investigate the influence of ageing on thermal runaway. Table 3 represents the three critical temperatures averaged for 3 fresh cells, 1 cell aged at 10% and 3 cells aged at 30%. The data were collected using the procedure described in Section 2.2.2. The results show that T1 increased with ageing, indicating that a higher temperature is needed to cause the onset of self-heating. In contrast, T3 decreased with increasing ageing, which means that the aged cells underwent thermal runaway at a lower temperature than fresh cells. As for T1, the critical temperature T2 related to the separator melting process increased with ageing, which might be caused by mechanical stress generated during the ageing process. However, the effect of ageing on the cell separator is currently not well known and would deserve more investigation. It is not possible to conclude on the positive effect of ageing on the thermal safety of the cells in spite of the delayed self-heating process. All the parameters examined here are important and have to be considered to improve the safety management of battery cells during their entire lifetime in various operating conditions [13].

Model calibration and validation

Model implementation and parameters setting

The thermal runaway model was developed for a cylindrical 26650 LFP/C cell. Since the battery cell case is made of metal with very good thermal conductivity and very small thickness (< 1 mm), its thermal influence was neglected and subsequently only the electrochemical part of the cell (positive electrode, negative electrode, and electrolyte) was considered in the modeling work. Moreover, the battery internal spiral wound structure was not introduced in the model so that the cell was considered as an homogenous system with equivalent thermal parameters for the different materials, as in the approach used by Chen et al. [28]. This approximation made it possible to reduce the number of input parameters needed to operate the model and to simplify its geometry. However, the anisotropy of the thermal conductivity was actually taken into account.

Typically, the thermal conductivity can be several tens of times less in the radial direction than in the axial direction. The equivalent thermal parameters given in Table 4 were determined by fitting experimental data obtained on a battery cell in quasi-adiabatic conditions in an accelerated rate calorimeter with a simple thermal model based on Equation ( 13).

The geometry of the cylindrical battery cell being simple, a regular mesh was used with 1,129 domain elements, 424 boundary elements, and 80 edge elements. The 3D thermal sub-model giving the evolution of the temperature gradient in the entire cell was implemented using the Heat Transfer in Solids module, in which the energy balance in the core of the cell (Equation ( 13))

corresponds to the differential form of the Fourier equation. Equations ( 3) to (6) accounting for the chemical degradation processes characterizing the thermal runaway were solved using the ODE (Ordinary Differential Equations) module to evaluate the temporal evolution of the concentrations of the species involved in the exothermic reactions. The equations were solved locally at each node of the mesh. Both chemical and thermal sub-models were coupled since the abuse heat source term given by the chemical sub-model was an input of the thermal sub-model.

Model calibration

The model was first calibrated for fresh A123 battery cells from experiments performed in the BTC to determine the degradation chemical reaction parameters (Table 5). For that, a lumped model, in which variables are considered uniform and only depend on time, was used to save calculation time without losing efficiency, as shown in a comparative study by Kim et al. [6]. This model allows the temperature of a LFP/C cell to be simulated in pseudo-adiabatic conditions representative of thermal abuse tests in a BTC in which the model input is the heat power provided by the heater wire surrounding the cell and the model output is the cell surface temperature defined as the average temperature of the two thermocouples.

Fig. 3 presents the simulated evolution of the surface temperature of the fresh cell shown in Fig. 2. The kinetic parameters (Table 5) have mostly been taken from literature [6,29]. However, the specific carbon content in the negative electrode, Wn,e, and the specific active material content in the positive electrode, Wp,e, were calculated from ante mortem analysis of A123 cells performed in the Simcal project [START_REF] Delaille | SIMCAL Project: calendar aging results obtained on a panel of 6 commercial Li-ion cells[END_REF] while the specific active material content in the electrolyte, We, was chosen as 20% of that of the whole cell. Large differences can be noted between the values of Wn,e

and Wp,e compared to literature values, while the We values are closer. Moreover, the specific heat released from the electrolyte, He, was fitted from the experimental data since the solvent composition of the electrolyte in the A123 cells identified in the Simcal project was different from that of the electrolyte in the literature [29]. The values of the activation energies ,Ea,SEI, Ea,ne, Ea,pe, Ea,e, were also fitted from the experimental data. Indeed, the exothermic degradation reactions of the different cell components (SEI, negative electrode, positive electrode and electrolyte) occur in cascade. It is then difficult to identify them separately using the BTC data, all the more so as the first two reactions related to the SEI and to the negative electrode are correlated through the Wn,e parameter. In addition, the degradations of the positive electrode and electrolyte solvent are coupled.

For the calendar ageing model, no parameter was fitted from experimental data. As shown in Table 6, most of the parameters were taken from literature. The dimensionless thickness, tSEI, of the SEI layer formed at the two levels of accelerated calendar ageing considered in this work (10% and 30% capacity loss) are also reported. These parameters, which were calculated from Equations ( 11) and ( 12), are the inputs of the cell degradation chemical reactions sub-model. In other words, the only effect of calendar ageing on thermal runaway was the change in value of the initial thickness of the SEI layer used in Equation ( 4), as described in Table 6.

Model validation

The calibrated thermal runaway model for the LFP / C chemistry was then validated for A123 cells with overheating tests carried out in the oven. Fig. 4 compares the experimental cell surface temperature measured during the oven test to the average surface temperature simulated with the thermal runaway 3D model. At first glance, it comes out that both temperature curves have a similar form even if the model better fits measured temperatures evolution at the beginning of the simulation (up to t  100 min), than during the later stages of the thermal runaway process.

Fig. 5 shows the corresponding evolution of the dimensionless average volume concentration and degree of conversion of the chemical species. It can be noted from the cSEI and cne variations that the chemical degradation of the SEI and negative electrode shown by, occurred at by occurred times between t = 90 and 100 min. Indeed, until that time the cell temperature prediction was in good agreement with the experimental data according to Fig. 4. By contrast, the major degradation of the electrolyte reflected in Fig. 5 by the steep ce decrease occurred beyond t = 100 min where the predicted cell temperature significantly differs from the temperature measured in the oven test.

The current limitation in the prediction capability of the model is mainly related to the poor estimation of the kinetic parameters associated to the complex degradation process of the electrolyte. However, a slight loss or airtight integrity of the oven door due to its deformation under thermal stresses was observed at a temperature around 170°C. This could be, at least partly, at the origin of the limited thermal dissipation that caused a slower increase in the measured temperature of the cell surface at t = 100 min (see Fig. 4). This thermal loss is obviously not taken into account in the model and results in the overestimation of the simulated peak temperature of the cell under the operating conditions.

Simulation results and discussion

Using this calibrated and validated modeling approach for the A123 cells, it is then possible with computer simulations to better understand the thermal runaway mechanism and to investigate the influence of battery ageing on this destructive phenomenon.

Thermal runaway onset temperature

Fig. 6 shows computer simulations of the average surface temperature of a fresh cell for different values of the external temperature, which is the parameter of major influence on thermal runaway. The initial temperature value was fixed at 20°C, and various abuse external temperatures (from 180°C to 250°C) were set as single increments at t = 0 minute. The heat exchange coefficient, heq, was set at 20 W K -1 m -2 as in the oven test (Fig. 4), which corresponds to forced convection conditions. At the beginning of the simulation, the cell average surface temperature increased at the same rate, independently of the external temperature. Then, for Text= 250°C the cell temperature rapidly reached a peak temperature of 335°C, indicating that thermal runaway occurred under these conditions, before decreasing, to stabilize at the external temperature of 250°C. For Text= 200°C the simulation shows a decrease in the cell heating rate between 180°C and 200°C, and a rapid increase from t = 40 minutes to reach 300°C in 10 minutes before cooling down In contrast, for Text = 180°C, the cell temperature increased gradually from t = 30 min up to 185°C and did not undergo thermal runaway

It can be concluded that the higher the external temperature, the higher the probability to initiate a thermal runaway, and the higher the severity of this thermal event. It has also been found that there is a critical temperature below which the cell does not go on thermal runaway.

Effect of heat transfer on thermal runaway

According to Fig. 6, below a critical external temperature, Text,c, ranging between 180°C and 200°C, no thermal runaway occurs. However, the result was shown for an equivalent heat transfer of rather high value (heq = 20 W K -1 m -2 ). It is then important to investigate the influence of this coefficient on the surface temperature of the cell as a function of the external temperature [21]. Fig. 7 shows the evolution of the average surface temperature of a fresh cell for two values of the external temperature, 180°C and 200°C, and several values of heq. When the external temperature is higher than Text,c (Fig. 7A), thermal runaway takes place for all values of the heat transfer coefficient since the cell reaches temperatures around 200°C at which violent degradation reactions already occurred, as that of the electrolyte. So, the external temperature is excessive for the cell even when decreasing the heat exchange, which only delays the time to reach 200°C and, therefore, thermal runaway. In contrast, when the external temperature is lower (180°C), no thermal runaway is expected according to Fig. 6 but the simulations in Fig. 7B show it does occur when the heat exchange is too low, due for example to a failure of the cooling system, to dissipate the heat outside the cell and prevent it from reaching abuse temperatures.

Effect of ageing on thermal runaway

The effect of accelerated calendar ageing on thermal runaway is analyzed through simulations of the cell temperature of cells aged at 10% and 30% of capacity loss exposed at an oven temperature Text = 180°C and a low heat exchange (heq = 5 W K -1 m -2 ). As for the fresh cell in Fig. 7B, the initial temperature of the cell is 20°C and the only difference is the thickness of the SEI layer corresponding to the ageing state (see Section 3.2 and Table 6). Fig. 8A shows a slower temperature increase for the aged cells from 60 to 140 min when the cell is mainly heated by the simulated oven up to a temperature of 180°C. From that time, the cell is heated by the degradation reactions and its temperature increases faster for the aged cells, before decreasing down to the temperature of the oven once the degradation reactions are complete. These results are consistent with those observed experimentally (Table 3).

To get a deeper understanding of the evolution of the cell temperature, Fig. 8B-E shows the effect of ageing on the average volume concentrations of the reactive species. It can be seen in Fig. 8B and 8C that the degradation reactions of the SEI and the negative electrode take place between 60 and 140 min for the fresh cell, causing a slight temperature increase at about 80 min.

As a consequence, after 140 min only two exothermal reactions happen (degradation of the positive electrode and electrolyte, see Fig. 8D and8E) instead of three for the aged cells since reactants are still present in the negative electrode at time 140 min (2min2.4% and 43.6% for the cells aged at 10% and 30%, respectively). Moreover, the power released by the aged cells is more important and their temperature increases faster and higher than for the fresh cell. The delay in degradation of the negative electrode can be explained by increased diffusional limitations of Li ions due to a thicker SEI layer after accelerated calendar ageing.

Conclusion

No model able to predict thermal runaway regardless its origin (electrical, thermal or mechanical abuse) has yet been published. No model either takes into account all physical phenomena involved in a Li-ion battery or considers all heat generation sources inside the battery.

Moreover, most of past studies deal with fresh materials or new cells only. Nevertheless, for using Li-ion batteries in growing high energy demanding applications, it is a key issue to guarantee safety for the entire lifetime of batteries and also for their second life targeting less demanding applications. The present work is a step forward in understanding the thermal behaviour of an aged LIB cell.

To investigate the mechanisms leading to thermal runaway, especially for aged Li-ion The primary effect of accelerated calendar ageing on the thermal stability of A123 cells was shown to be the onset of the first exothermic reaction at a higher temperature for the aged cells. In addition, the aged cells were revealed to undergo thermal runaway at a lower temperature than new cells. Computer simulations were able to reproduce these experimental results.

Moreover, the simulated evolution of the reactive species concentration provided evidence of the relationship between the diffusional limitations generated by the growth of the SEI layer induced by calendar ageing and the delay of the self-heating onset temperature. Besides calendar ageing considered in this work, it is important investigate the influence of other ageing processes on thermal stability, as for example Li-plating that is expected to be faced in high energy cells.

Other perspectives can be put forward to improve the thermal runaway model. Some events

Figure and table captions (Table list)

Table 1.

Characteristics of the cell under investigation (EC: Ethylene Carbonate, PC: Propylene Carbonate, EMC: Ethyl-Methyl Carbonate, NMP: N-Methyl Pyrrolidone, PP: PolyPropylene).

Table 2.

Models and experimental tests used in this work.

Table 3.

Evolution of the critical temperatures related to the thermal safety of A123 cells with ageing (from experimental data in the BTC).

Table 4.

Thermal parameters of the A123 cells.

Table 5. Table 6.

Parameters of the degradation

Parameters of the calendar ageing model for the A123 LFP/C cells: (a) from Ref. [24], (b) from

Ref. [6], (c) calculated. Ea,SEI 1.38×10 5 (a) 1.35×10 5 (b) Ea,ne 1.32×10 5 (a) 1.35×10 5 (b) Ea,pe 0.99×10 5 (a) 1.03×10 5 (c) Ea,e 2.70×10 5 (a) 2.74×10 5 (b) Reaction factors (s -1 ) ASEI 1.66×10 15 (b) 1.66×10 15 (b) Ane 2.50×10 13 (b) 2.50×10 13 (b) Ape 2×10 8 (c) 2×10 8(c) Ae 5.14×10 25 (b) 5.14×10 25 (b) Reaction orders (-) Wn,e 220 (d) 1700 (b) Wp,e 520.74 (d) 960 (c) We 334.68 (d) 500 (b) Table 6 Parameters of the calendar ageing model for the A123 LFP/C cells: (a) from Ref. [24], (b) from

(Figure list)

mSEI 1 (b) 1 (b) mne 1 (b) 1 (b) p pe m , r pe m 1 (b) 1 (b)
Ref. [6], (c) calculated.

Parameters A123

MSEI (kg mol -1 ) 0.162 (a) SEI (kg m -3 ) 1690 (a) s,n 0.58 (a) δn (m) 3.45×10 -5 (a)

A (m²) 0.18 (a) Rs,n (m) 5×10 -6 (a)

Closs (Ah) 10% capacity loss 0.23 (c) 30% capacity loss 0.69

  than 0.115 A (C/20). The discharge was performed at a constant current rate of 1 C until a 2 V cell voltage. The internal resistance was determined from the voltage relaxation after a discharge current step of 2 C rate during 1 s on cells charged at a 50% SOC. A 50 V -200 A Digatron test bench and a Biologic multipotentiostat were used in addition to climatic chambers thermoregulated at 25°C for the experimental electrical characterization of the cells. Accelerated calendar ageing of the cells under investigation was obtained by placing fully charged cells in a climatic chamber thermo-regulated at 60°C. Considering an end of first life for LIBs usually consisting in 20% capacity loss, two levels of cell capacity loss, 10% and 30% corresponding respectively to mid-range of first life and to beginning of second life, were targeted and obtained by means of periodic measurements of remaining capacity all along the accelerated calendar ageing process. Higher capacity losses were not considered in this work since ageing phenomena other than the SEI growth would require the development of new sub-models. The cells collected after 75 (resp. 340) days of calendar ageing showed a capacity loss of 10% ± 3%

  Some other limitations and shortcomings related to the validation process of the model should be mentioned in addition to the fact that the model considers only the contribution of the degradation chemical reactions of the cell components during thermal runaway. Under abuse testing conditions, the cylindrical cells were observed to vent once a critical internal pressure was reached. In such conditions, the further heating of the cell results in vaporization of the organic solvents used in the electrolyte. Additional degradation chemical reactions between the electrolyte and the cell components also lead to accelerated gas generation. Depending on the overpressure reached, cell venting can lead to ejection of the electrolyte or even the jellyroll contents. Thus, the heat dissipation and the reactive mass loss due to such ejection should be considered in a more robust thermal runaway model as shown by Coman et al.[23] in their lumped model of venting during thermal runaway for a cylindrical LCO cell. During the experimental study, separator melting was observed to lead most of the time to almost immediate strong internal short circuit, which is not taken into account in the present model, and then to the associated thermal runaway. Contribution of internal short circuits could be considered in an improved version of the model by introducing a specific electrochemical sub-model or localized hot spots in the thermal model with corresponding heat releases.

  batteries, it was necessary to adopt a strategy combining modeling and experimental work to achieve a multiphysics 3D model of thermal runaway of fresh and aged LIB cells. The model was structured into three coupled sub-models, a 3D thermal model, a 3D degradation chemical reaction model, and a 0D calendar ageing model based on a simplified SEI growth mechanism. This model was developed at the cell level and calibrated for the LFP/C chemistry. The thermal parameters were fitted from experimental data while the kinetics model parameters were taken from literature or fitted from calorimetric tests. Accelerated calendar ageing tests allowed the calibration of the ageing model at two different levels of capacity loss. The loss thermal runaway model including ageing was then validated for commercialized A123cylindrical cells by comparing simulations to cell responses measured in oven abuse tests.

  observed during the abuse tests, which would deserve due consideration, are not yet included in the current version of the model, such as the separator melting process and cell venting. Separator melting, which causes an internal short circuit that triggers the cell thermal runaway can be treated by considering the internal short circuit as a localized hot spot in a future refinement of the thermal runaway model. The development of a coupled electrochemical-thermal short circuit model would also deserve a dedicated research effort. With regard to the cell venting issue, it could be taken into account by introducing a pressure dependent relationship associated with the evaporation of the electrolyte solvent. Calibration of the model for application to other chemistries is obviously another direction for future work towards a more valuable tool to optimize the design of lithiumion batteries for safety management purposes.

  chemical reaction model for the A123 LFP/C cells: (a) fitted from our calibration work, (b) from Ref. [6], (c) from Ref. [29], (d) calculated.

Fig. 1 .

 1 Fig. 1. Thermal runaway coupled multiphysics model.

Fig. 2 .Fig. 3 .

 23 Fig. 2. Surface temperature and voltage of a fresh A123 battery cell measured during the calibration experiment in the BTC.Fig. 3. Comparison of the simulated average surface temperature of a fresh A123 battery cell with the experimental data in Fig. 3.

Fig. 4 .

 4 Fig. 4. Experimental and simulated evolution of the average surface temperature of a fresh A123 cell during the oven test.

Fig. 5 .

 5 Fig. 5. Simulated evolution of the dimensionless concentrations of the reactive species during the oven test.

Fig. 6 .

 6 Fig. 6. Simulation of the average surface temperature of a fresh A123 cell as a function of the external temperature (heq = 20 W K -1 m -2 ).

Fig. 7 .

 7 Fig. 7. Influence of the equivalent heat transfer coefficient on the average surface temperature of a fresh A123 cell: Text = 200°C (A), Text = 180°C (B).

Fig. 8 .

 8 Fig. 8. Simulated evolution of the cell average surface temperature (A) and dimensionless concentrations of the reactive species (B to E) during the oven test (T ext = 180°C, h eq = 5 W K -1 m -2 ).

5

 5 Parameters of the degradation chemical reaction model for the A123 LFP/C cells: (a) fitted from our calibration work, (b) from Ref. [6], (c) from Ref. [29], (d) calculated.

Fig. 1 .Fig. 2 .Fig. 3 .Fig. 4 .

 1234 Fig. 1. Thermal runaway coupled multiphysics model.

Fig. 5 .

 5 Fig. 5. Simulated evolution of the dimensionless concentrations of the reactive species during the oven test.

Fig. 6 .

 6 Fig. 6. Simulation of the average surface temperature of a fresh A123 cell as a function of the external temperature (heq = 20 W K -1 m -2 ).

Fig. 7 .

 7 Fig. 7. Influence of the equivalent heat transfer coefficient on the average surface temperature of a fresh A123 cell: Text = 200°C (A), Text = 180°C (B).

Fig. 8 .

 8 Fig. 8. Simulated evolution of the cell average surface temperature (A) and dimensionless concentrations of the reactive species (B to E) during the oven test (T ext = 180°C, h eq = 5 W K -1 m -2 ).

Table 1

 1 Characteristics of the cell under investigation (EC: Ethylene Carbonate, PC: Propylene Carbonate, EMC: Ethyl-Methyl Carbonate, NMP: N-Methyl Pyrrolidone, PP: PolyPropylene).

	Cell technology	A123
	Chemistry	LFP/C
	Packaging	Cylindrical
	Nominal capacity	2.3 Ah
	Nominal voltage	3.2 V
	Dimensions (D * h) 26 mm * 65 mm
	Weight	77 g
	Electrolyte	LiPF6 in EC-PC-EMC-NMP
	Separator	PP

Table 2

 2 Models and experimental tests used in this work.

	Models	Calendar ageing	Cell degradation	Thermal	Thermal
		sub-model	chemical reactions	sub-model	runaway
			sub-model		model
	Experimental	-Initial state cell	Battery Test	Characterization	Oven: external
	tests	characterization	Calorimetry	through fine	heating test
		-Accelerated	(for calibration)	tuning thermal	(for validation)
		calendar ageing		parameters
		-Regular check-up

Table 3 .

 3 Evolution of the critical temperatures related to the thermal safety of A123 cells with ageing (from experimental data in the BTC).

	Temperature	T1 (°C)	T2 (°C)	T3 (°C)
	Fresh	79.7	137.7	171.3
	Aged at 10%	83	145.7	154
	Aged at 30%	92	147	150.7

Table 4

 4 Thermal parameters of the A123 cells.