M Brodholt

Mezouar

Martha G Pamato
email: martha.pamato@unipd.it

Y Li

D Antonangeli

F Miozzi

G Morard

I G Wood

L Vočadlo

J P Brodholt

M Mezouar

Equation of State of hcp Fe-C-Si Alloys and the Effect of C Incorporation Mechanism on the Density of hcp Fe alloys at 300 K

Equation of State of hcp

Fe-C-Si Alloys and the Effect of C Incorporation Mechanism on the Density of hcp Fe alloys at 300 K 2

Introduction

The only direct investigation of the inaccessible Earth's deep interior is provided by seismological studies of sound wave velocities through the Earth. The comparison between seismological observations and elastic properties of candidate materials at the relevant conditions present in the deep Earth, allows to relate geophysical data with the chemistry and temperature of the Earth's interior. In the Earth's mantle, mineralogical and compositional models can match the observed velocities to approximately 1% (e.g. [START_REF] Irifune | Sound velocities of majorite garnet and the composition of the mantle transition region[END_REF][START_REF] Murakami | A perovskitic lower mantle inferred from high-pressure, high-temperature sound velocity data[END_REF]Bass & Zhang, 2015;[START_REF] Pamato | Single crystal elasticity of majoritic garnets: Stagnant slabs and thermal anomalies at the base of the transition zone[END_REF][START_REF] Kurnosov | Evidence for a Fe3+-rich pyrolitic lower mantle from (Al,Fe)-bearing bridgmanite elasticity data[END_REF]. Conversely, for the inner core, sound wave velocities of pure Fe or Fe-Ni alloys in the hcp structure usually obtained from mineral physics are much higher (10-30 %) than those from seismic observations [START_REF] Vočadlo | Ab initio calculations of the elasticity of iron and iron alloys at inner core conditions: Evidence for a partially molten inner core?[END_REF][START_REF] Belonoshko | Origin of the Low Rigidity of the Earth's Inner Core[END_REF]Martorell, Brodholt, et al., 2013). Numerous arguments were proposed over the years to explain these discrepancies, including, for example, large compositional effects, pervasive partial melting, anelasticity or the stabilization of bcc Fe [START_REF] Vočadlo | Ab initio calculations of the elasticity of iron and iron alloys at inner core conditions: Evidence for a partially molten inner core?[END_REF][START_REF] Belonoshko | Stability of the body-centred-cubic phase of iron in the Earth's inner core[END_REF][START_REF] Belonoshko | Origin of the Low Rigidity of the Earth's Inner Core[END_REF][START_REF] Belonoshko | Stabilization of body-centred cubic iron under inner-core conditions[END_REF][START_REF] Belonoshko | Low viscosity of the Earth's inner core[END_REF][START_REF] Antonangeli | Elastic anisotropy in textured hcp-iron to 112 GPa from sound wave propagation measurements[END_REF].

Theoretical calculations have shown that the sound velocities observed in the inner core could be matched by hcp Fe due to a strong non linear softening of the elastic properties at 360 GPa above 7000 K (Martorell, Vočadlo, et al., 2013). However, even though the velocities of the inner core were matched by pure Fe, the density remained too high (by 2-3%) to agree with geophysical observations, and so lighter elements are still needed to explain this density mismatch.

Furthermore, the suggested pre-melting softening is generally suppressed in iron alloys (e.g. [START_REF] Martorell | The elastic properties of hcp-Fe 1-x Si x at Earth's inner-core conditions[END_REF]. With the exception of Fe 7 C 3 [START_REF] Chen | Hidden carbon in Earth's inner core revealed by shear softening in dense Fe 7 C 3[END_REF]Prescher et al., 2015), all the iron compounds examined so far exhibit higher velocities than those found in the core. For the case of Fe 7 C 3, a strong non linear softening has been also observed [START_REF] Li | Thermoelasticity of Fe 7 C 3 under inner core conditions[END_REF], but, while the calculated velocities agree with seismic data, the density is far too low (by ~ 8%).

Many recent computational and experimental studies addressing core compositions point to the need of two or more distinct light elements (e.g. Badro et al. 2104;Morard, Nakajima, et al., 2017;[START_REF] Li | The elastic properties of hcp-Fe alloys under the conditions of the Earth's inner core[END_REF][START_REF] Edmund | Velocity-Density Systematics of Fe-5wt%Si: Constraints on Si Content in the Earth's Inner Core[END_REF]. However, the identity of this light element mixture is still largely debated. Very recently [START_REF] Li | The elastic properties of hcp-Fe alloys under the conditions of the Earth's inner core[END_REF], reported ab initio molecular dynamics (AIMD) calculations which reveal the observed inner-core density and sound velocities can be simultaneously matched by a range of ternary and quaternary Fe alloys, all containing carbon. This constrains possible solid solution's space in a multiple light elements model explaining all of the major seismic constraints. However, even when limiting the discussion to mineral physics properties, several aspects still remain uncertain, including the different incorporation mechanisms of C in Fe and its effect on seismic properties. From density functional theory (DFT) calculations, [START_REF] Huang | Carbon in iron phases under high pressure[END_REF] suggested that C takes the substitutional form in hcp Fe at 0 K above approximately 300 GPa, as confirmed by [START_REF] Li | Carbon Partitioning Between the Earth's Inner and Outer Core[END_REF] who, using ab initio molecular dynamics, also found that C takes the substitutional form at 360 GPa and 6500 K.

However, at lower pressures, interstitial C and other complex defect structures may be stabilized [START_REF] Huang | Carbon in iron phases under high pressure[END_REF][START_REF] Caracas | The influence of carbon on the seismic properties of solid iron[END_REF][START_REF] Li | Carbon Partitioning Between the Earth's Inner and Outer Core[END_REF]. In fact, the thermodynamic model of [START_REF] Fei | Experimental study and thermodynamic calculations of phase relations in the Fe-C system at high pressure[END_REF], based on experiments at pressures <20 GPa, and more recently the work by [START_REF] Yang | Effect of Carbon on the Volume of Solid Iron at High Pressure: Implications for Carbon Substitution in Iron Structures and Carbon Content in the Earth's Inner Core[END_REF] indicate that the interstitial form is preferred up to inner-core pressures. The change of incorporation mechanism may be related to the variance of solubility. [START_REF] Lord | Melting in the Fe-C system to 70 GPa[END_REF], for example, found a negligible C solubility in hcp Fe approaching 50 GPa. However, more recently, [START_REF] Mashino | Melting experiments on the Fe-C binary system up to 255 GPa: Constraints on the carbon content in the Earth's core[END_REF] performed melting experiments up to 255 GPa and predicted a solubility of 1 wt.% in hcp Fe at inner-core boundary conditions. A question remains as to whether different incorporation mechanisms of C in iron have an effect on the seismic properties of an Fe-C alloy; does it matter if the carbon is on an interstitial, substitutional or other site within the hcp crystal structure?

Multiple lines of evidence, from geochemical and cosmochemical arguments (e.g. [START_REF] Allègre | The chemical composition of the Earth[END_REF], to core formation models based on metal-silicate partitioning (e.g. [START_REF] Siebert | Terrestrial accretion under oxidizing conditions[END_REF][START_REF] Fischer | High pressure metal-silicate partitioning of Ni, Co, V, Cr, Si, and O[END_REF], to isotopic considerations (e.g. [START_REF] Fitoussi | Si isotope systematics of meteorites and terrestrial peridotites: implications for Mg/Si fractionation in the solar nebula and for Si in the Earth's core[END_REF][START_REF] Shahar | Experimentally determined Si isotope fractionation between silicate and Fe metal and implications for Earth's core formation[END_REF], support silicon to be present in the core. Due to the almost equal Si partitioning between liquid and hcp iron [START_REF] Alfè | Composition and temperature of the earth's core constrained by combining ab initio calculations and seismic data[END_REF], the quantity of Si present in the inner core is not significantly different from that expected for the liquid core. However, many recent studies point out that Si is very unlikely to be the only light element in the core [START_REF] Badro | A seismologically consistent compositional model of Earth's core[END_REF]Morard, Andrault, et al., 2017) and specifically in the inner core [START_REF] Martorell | The elastic properties of hcp-Fe 1-x Si x at Earth's inner-core conditions[END_REF][START_REF] Antonangeli | Sound velocities and density measurements of solid hcp-Fe and hcp-Fe-Si (9 wt.%) alloy at high pressure: Constraints on the Si abundance in the Earth's inner core[END_REF][START_REF] Edmund | Velocity-Density Systematics of Fe-5wt%Si: Constraints on Si Content in the Earth's Inner Core[END_REF]. An alloy in the ternary Fe-C-Si system thus seems to provide the simplest composition possibly accounting for the bulk of geophysical and geochemical evidences of the inner core. In particular, very recent AIMD calculations [START_REF] Li | The elastic properties of hcp-Fe alloys under the conditions of the Earth's inner core[END_REF] indicate that an hcp Fe 60 C 2 Si 2 alloy is expected to have density, compressional and shear sound velocities matching those of PREM at 360 GPa for simulation temperatures between 6000 and 7000 K.

To confirm this prediction, both high-pressure experiments and further simulations need to be performed to clarify the effects of C on the seismic properties of Fe and Fe-Si alloys over a wider range of pressures and temperatures. Thus, here we address the effect of the C incorporation mechanism on density by combining synchrotron X-ray diffraction (XRD) measurements on an hcp Fe-C-Si alloy (with 4 at % C and 3 at % Si at 300 K and up to ∼150 GPa), with DFT calculations. In particular, we examined the different incorporation mechanisms of carbon in the structure of the alloy as a function of pressure. Our results are then used to discuss the density of a Fe-C-Si alloy at inner core conditions.

Materials and Methods

Experiments

The starting material, with a nominal composition of Fe 93 at%, C 4 at%, Si 3 at%, was synthesised by an ultra-rapid quench method at the ICMPE (Institut de Chimie et des Matériaux de Paris-Est, Paris, France) using nominally pure Fe (Neyco F-12734), Si (Wacker) and Fe 3 C (Neyco FC-76857/1). Homogeneous Fe-alloys, in the form of ribbons, approximately 25 µm thick and of 20-30 mm width, were synthesized following the procedure described in [START_REF] Morard | Melting of Fe-Ni-Si and Fe-Ni-S alloys at megabar pressures: Implications for the coremantle boundary temperature[END_REF].

The chemical composition of the synthesized alloys was measured by means of a Cameca SX100 electron microprobe analyser equipped with a cold finger and an O 2 flux at Camparis center-Sorbonne Université. The resulting composition for the alloy is 3 (±0.2) at% Si, 4 (±0.4) at% C and 93 (±0.06) at% Fe (Fe 93 C 4 Si 3 ). In addition, scanning electron microscopy (SEM) analyses were performed on the sample using a Zeiss Ultra 55 scanning electron microscope installed at the IMPMC (Institut de Minéralogie, de Physique des Matériaux et de Cosmochimie), confirming homogeneity at the micron scale [START_REF] Miozzi | Eutectic melting of Fe-3 at% Si-4 at% C up to 200 GPa and implications for the Earth's core[END_REF].

Samples were prepared for the high-pressure experiments by crushing parts of a ribbon specimen and compressing it between two diamonds (~ 600 μm culet), in order to obtain the desired thickness (usually between 5 and 10 microns). The sample was then loaded in a rhenium gasket of 20 μm thickness, with a hole of 60 μm, with 100/300 μm bevelled anvils. High pressures were generated using a Le Toullec-type (membrane driven) diamond anvil cell. Tungsten carbide seats designed for diamonds with conical support were used [START_REF] Boehler | New anvil designs in diamond-cells[END_REF]. The sample was loaded alongside Mo as a pressure calibrant, with Ne acting as the pressure transmitting medium; -the equation of state parameters for the Mo pressure calibrant were taken to be: V 0 = 31.14 Å 3 [START_REF] Litasov | Thermal equation of state and thermodynamic properties of molybdenum at high pressures[END_REF][START_REF] Ross | High temperature X-ray metallography. I. A new debye-scherrer camera for use at very high temperatures II. A new parafocusing camera III. Applications to the study of chromium, hafnium, molybdenum, rhodium, ruthenium and tungsten[END_REF], K 0 = 260(1) GPa and K' = 4.19(5).

Angle-dispersive XRD measurements were carried out at the ID27 beamline of the European Synchrotron Radiation Facility in Grenoble [START_REF] Mezouar | Development of a new state-of-the-art beamline optimized for monochromatic singlecrystal and powder X-ray diffraction under extreme conditions at the ESRF[END_REF]. In the experiments, a monochromatic beam of 33 keV (λ= 0.3738 Å) was focused on an area < 3 × 3 μm 2 (full width at half maximum). Data collection times were typically between 10 and 30 s; diffraction patterns were collected on a MAR CCD detector at pressures between 48 and 150 GPa at 300 K.

2-D diffraction images were integrated using the software DIOPTAS [START_REF] Prescher | DIOPTAS: A program for reduction of twodimensional X-ray diffraction data and data exploration[END_REF], using calibration parameters (sample detector distance and orientation) derived from the diffraction pattern of a cerium dioxide (CeO 2 ) standard. Diffraction patterns were indexed with PD Indexer (http://pmsl.planet.sci.kobe148u.ac.jp/~seto/?page_id=20). The sample started to transform from the bcc structure stable at ambient pressure to the hcp structure at ~19 GPa. All bcc reflections were absent by ~22 GPa. At least five sample diffraction lines are clearly observed for the hcp phase over the entire pressure range. However, the diffraction patterns are, to some extent, non-ideal in terms of randomness of the distribution (both 101 and 102 reflections showed intensity variation around the powder ring). Also, the 002 reflection, while present over the entire pressure range, is very weak, due to a small preferential orientation, likely induced by sample thinning before loading. Therefore, it was not possible to perform full profile data refinements. The unit-cell parameters were thus calculated from the d-spacings of the reflections 100, 101, 102, and 110 using the software Unit Cell by [START_REF] Holland | UNITCELL: A nonlinear least-squares program for cell-parameter refinement implementing regression and deletion diagnostics[END_REF].

DFT Calculations

DFT calculations were performed by using a 4 × 4 × 2 hcp supercell, the VASP code [START_REF] Kresse | From ultrasoft pseudopotentials to the projector augmentedwave method[END_REF] and PAW potential [START_REF] Blöchl | Projector augmented-wave method[END_REF]. The valence states 3p 6 3d 6 4s 2 , 3s 2 3p 2 and 2s 2 2p 2 were adopted for Fe, Si and C, respectively. We employed the generalised gradient approximation (GGA); the Perdew-Wang [START_REF] Perdew | Accurate and simple analytic representation of the electrongas correlation energy[END_REF] exchange-correlation functional was used. We adopted the ordered symmetrical structure for the alloy by placing solid solution atoms with the maximum sum of distances between them to achieve homogeneous dispersion. The octahedral interstitial site was adopted for the interstitial carbon atoms as the tetrahedral interstitial site was reported to be less favourable [START_REF] Caracas | The influence of carbon on the seismic properties of solid iron[END_REF]. A plane-wave cutoff energy of 400 eV and 12 irreducible k-points were used to converge the results.

Results and Discussion

Equations of state

The diffraction pattern of the investigated alloy collected at 140 GPa (the upper end of the pressure range investigated) and indexed as hexagonal closed packed (hcp) iron, is shown in Figure 1.

Measured unit-cell volumes of hcp Fe 93 C 4 Si 3 alloy are presented in Figure 2 and show a smooth, continuous trend with increasing pressure, indicating the absence of any phase transition within the pressure range investigated.

The P-V data, reported in Table S1, were fitted using a third-order Birch-Murnaghan equation of state (BM3 EoS; [START_REF] Birch | Finite elastic strain of cubic crystals[END_REF] with the software EosFit7 [START_REF] Angel | EosFit7c and a Fortran module (library) for equation of state calculations[END_REF]. All the parameters were refined resulting in the following EoS parameters: V 0 = 22.7(2) Å 3 , K 0 = 148(18) GPa, and 1). The EoS parameters are usually defined at 0 GPa (e.g., V 0 , K 0 and K 0 ʹ).

K 0 ʹ = 5.8(4) (Table
However, as iron alloys in the hcp-structure are not stable and cannot be quenched at room pressure and temperature, the EoS parameters can be defined at any other pressure in order to compare different EoS in the pressure stability field of hcp-phases. The parameters refined relative to a pressure scale from 50 GPa are: V P=50 = 18. 80(1) Å 3 , K P=50 = 394(6) GPa, andKʹ P=50 = 4.8(2).

The measured unit-cell volumes of Fe 93 C 4 Si 3 are larger than those of pure hcp Fe, Fe-Si, and Fe-C alloys reported in the literature (Figure 2a). If compared to the pure hcp Fe EoS as determined by [START_REF] Dewaele | Quasihydrostatic equation of state of Iron above 2 Mbar[END_REF], the P-V trend established here shows a different curvature, with the difference in the volumes increasing at high pressure, thus highlighting a difference in compressibility. On the other hand, differences are less pronounced when compared to recent EoS of pure hcp Fe reported by [START_REF] Sakai | Equation of state of pure iron and Fe0.9Ni0.1 alloy up to 3 Mbar[END_REF], [START_REF] Yang | Effect of Carbon on the Volume of Solid Iron at High Pressure: Implications for Carbon Substitution in Iron Structures and Carbon Content in the Earth's Inner Core[END_REF] and [START_REF] Miozzi | A new reference for the thermal equation of state of iron[END_REF]. Specifically, the new data on the Fe 93 C 4 Si 3 alloy determined here, show a compressibility curve comparable to the more recent results for pure hcp Fe, but with higher volumes.

Comparison of Fe-C-Si alloys

A comprehensive comparison of the equations of state parameters K 0 and K 0 ʹ of the different alloys considered here is illustrated in the inset in Figure 2a. The data published by [START_REF] Dewaele | Quasihydrostatic equation of state of Iron above 2 Mbar[END_REF], [START_REF] Sakai | Equation of state of pure iron and Fe0.9Ni0.1 alloy up to 3 Mbar[END_REF], [START_REF] Yang | Effect of Carbon on the Volume of Solid Iron at High Pressure: Implications for Carbon Substitution in Iron Structures and Carbon Content in the Earth's Inner Core[END_REF] and [START_REF] Miozzi | A new reference for the thermal equation of state of iron[END_REF] for pure Fe and by [START_REF] Edmund | Velocity-Density Systematics of Fe-5wt%Si: Constraints on Si Content in the Earth's Inner Core[END_REF] and [START_REF] Yang | Effect of Carbon on the Volume of Solid Iron at High Pressure: Implications for Carbon Substitution in Iron Structures and Carbon Content in the Earth's Inner Core[END_REF] for Fe-Si and Fe-C alloys, respectively, were reevaluated in order to obtain the variance-covariance matrixes [START_REF] Allègre | The chemical composition of the Earth[END_REF] and to construct the confidence ellipses (Table 1, Inset Figure 2a). The confidence ellipses of our dataset are larger than those reported in the literature. This might be due to the different pressure and volume uncertainties, how these are taken into account in the analysis, and to the number and range of measured pressure points in the different studies. Here, both pressure and volume uncertainties were weighted, whereas most of the studies from the literature report an equallyweighted fit. When performing an equally-weighted fit, errors are smaller and thus also the correlation between K 0 and K 0 ʹ is smaller. In addition, given the trade-off between V 0 , K 0 and K 0 ʹ, it is important to measure the volume of the hcp structure at the lowest pressure in order to constrain the refinement of the room temperature EoS (see e.g., [START_REF] Miozzi | A new reference for the thermal equation of state of iron[END_REF]). In our study, due to the lack of low-pressure points (we report P-V data from 48 GPa, see Table S1)

the EoS refinement provides a less robust constraint of the reference volume at P = 0 (i.e. V 0 ).

The ellipse representing our data and the one of [START_REF] Yang | Effect of Carbon on the Volume of Solid Iron at High Pressure: Implications for Carbon Substitution in Iron Structures and Carbon Content in the Earth's Inner Core[END_REF] for pure Fe intersect at the 68.3% confidence level. Our data are also in good agreement with the recent EoS of Fe reported by [START_REF] Miozzi | A new reference for the thermal equation of state of iron[END_REF], which extends the data set to lower pressures, and also, within experimental error, with that of [START_REF] Sakai | Equation of state of pure iron and Fe0.9Ni0.1 alloy up to 3 Mbar[END_REF]. The confidence ellipses are in fact parallel and the error bars intersect. All this indicates that the value of the isothermal room pressure bulk modulus obtained in this study for Fe 93 C 4 Si 3 is very similar to the values obtained for pure Fe from [START_REF] Sakai | Equation of state of pure iron and Fe0.9Ni0.1 alloy up to 3 Mbar[END_REF], [START_REF] Yang | Effect of Carbon on the Volume of Solid Iron at High Pressure: Implications for Carbon Substitution in Iron Structures and Carbon Content in the Earth's Inner Core[END_REF] and [START_REF] Miozzi | A new reference for the thermal equation of state of iron[END_REF]. A critical review of the EoS parameters of pure Fe is beyond the aim of this paper. It is important to note, however, that the confidence ellipse built with the K 0 and K 0 ʹ values of [START_REF] Miozzi | A new reference for the thermal equation of state of iron[END_REF] intersects only with that of [START_REF] Sakai | Equation of state of pure iron and Fe0.9Ni0.1 alloy up to 3 Mbar[END_REF] for pure Fe, while, that of [START_REF] Dewaele | Quasihydrostatic equation of state of Iron above 2 Mbar[END_REF] does not intersect with our results and also does not intersect with other results in the literature for pure Fe, as a direct consequence of the different compressibility already mentioned. [START_REF] Edmund | Velocity-Density Systematics of Fe-5wt%Si: Constraints on Si Content in the Earth's Inner Core[END_REF], recently suggested that hcp Fe can contain up to 5 wt% Si (approximately compression curve and the confidence ellipses, in fact, the data by [START_REF] Edmund | Velocity-Density Systematics of Fe-5wt%Si: Constraints on Si Content in the Earth's Inner Core[END_REF] are in good agreement with the data from [START_REF] Dewaele | Quasihydrostatic equation of state of Iron above 2 Mbar[END_REF] for hcp Fe. [START_REF] Yang | Effect of Carbon on the Volume of Solid Iron at High Pressure: Implications for Carbon Substitution in Iron Structures and Carbon Content in the Earth's Inner Core[END_REF] reported that the incorporation of C has a direct effect on the volume of hcp Fe. In particular, their fitting results showed that addition of C to Fe leads to larger unit cell volumes at the same pressure, with the magnitude of the increase remaining constant up to ~135 GPa, the highest investigated pressure, with compression curves almost parallel. The authors fitted a BM3 EoS with K 0 ′ fixed to 4.79 (value taken from [START_REF] Fei | Thermal equation of state of hcp-iron: Constraint on the density deficit of Earth's solid inner core[END_REF] for pure Fe) to systematically compare the EoS parameters for pure Fe, Fe-0.31 wt% C (~1.4 at% C), and Fe-1.37 wt% C (~6.1 at% C).

However, when the data of [START_REF] Yang | Effect of Carbon on the Volume of Solid Iron at High Pressure: Implications for Carbon Substitution in Iron Structures and Carbon Content in the Earth's Inner Core[END_REF] are fitted without constraining K 0 ′, the resulting magnitude of the volume increase is not the same. Furthermore, it can be seen that extrapolation of the EoS established for Fe-6.1 at% C and Fe-1.4 at% C crosses above 120 GPa, with the volumes of the first becoming smaller than the second and approaching the volume of pure Fe (Figure 2b). Even when limiting comparison to the directly investigated pressure range, we note that Fe-6.1 at% C has smaller volumes than the sample studied here, which contains less C (4 at% C), and volumes similar to those of Fe-Si alloy with 8 at% Si [START_REF] Edmund | Velocity-Density Systematics of Fe-5wt%Si: Constraints on Si Content in the Earth's Inner Core[END_REF]) and pure Fe by [START_REF] Dewaele | Quasihydrostatic equation of state of Iron above 2 Mbar[END_REF]. Moreover, as shown in the inset in Figure 2a, the confidence ellipses representing our data and those of [START_REF] Yang | Effect of Carbon on the Volume of Solid Iron at High Pressure: Implications for Carbon Substitution in Iron Structures and Carbon Content in the Earth's Inner Core[END_REF] on pure Fe and Fe-1.4 at % C intersect, indicating a shift only in V 0 without changing the compressibility, while that of Fe-6.1 at. %C is closer to those of Fe-8 at% Si [START_REF] Edmund | Velocity-Density Systematics of Fe-5wt%Si: Constraints on Si Content in the Earth's Inner Core[END_REF]) and pure Fe according to [START_REF] Dewaele | Quasihydrostatic equation of state of Iron above 2 Mbar[END_REF].

Effect of the incorporation mechanism of carbon on the compressibility of hcp

Fe

While Si is incorporated via substitution in iron without introducing major changes in the volume [START_REF] Edmund | Velocity-Density Systematics of Fe-5wt%Si: Constraints on Si Content in the Earth's Inner Core[END_REF], the strong effect of small amounts of C on the volume of hcp Fe may well be related to its interstitial position in the hcp structure [START_REF] Caracas | The influence of carbon on the seismic properties of solid iron[END_REF][START_REF] Yang | Effect of Carbon on the Volume of Solid Iron at High Pressure: Implications for Carbon Substitution in Iron Structures and Carbon Content in the Earth's Inner Core[END_REF].

Impurities, such as C, in a solid can be present as substitutional or interstitial atoms. In addition, impurities can develop defect clusters, possibly modifing the physical and mechanical properties of the material [START_REF] Smidt | PROPERTY CHANGES RESULTING FROM IMPURITY-DEFECT INTERACTIONS IN IRON AND PRESSURE VESSEL STEEL ALLOYS[END_REF][START_REF] Li | Impurity effects on the grain boundary cohesion in copper[END_REF][START_REF] Li | Carbon Partitioning Between the Earth's Inner and Outer Core[END_REF] including its volume and its compressibility. Therefore, it is essential to address the different incorporation mechanisms of carbon in the structure of the alloy at different pressures, to evaluate their relevance in defining the properties of the Earth's core. To this end, in our ab initio calculations, we considered both binary Fe-C and ternary Fe-C-Si alloys. At first we relaxed the structures and calculated enthalpies for the binary substitutional Fe 62 C 2 , interstitial Fe 64 C 2 and dimer Fe 63 C 2 (two C at one atomic site) alloys (see the structural illustration in Figure S1). Then, we also calculated volumes and fitted the EoS of pure Fe, Fe 60 C 2 Si 2 via substitution, Fe 62 C 2 Si 2 via interstitial carbon and Fe 61 C 2 Si 2 with a dimer (Table S2 andTable 1, respectively). The composition of the calculated alloys corresponds to Fe-~3 at% C and ~3 at% Si.

Calculations using a 16 Fe atoms supercell [START_REF] Huang | Carbon in iron phases under high pressure[END_REF]) support a transition from interstitial to substitutional C at pressure of ~ 300 GPa at 0 K. At high temperature, calculations based on the quasiharmonic approximation (QHA) [START_REF] Li | Carbon Partitioning Between the Earth's Inner and Outer Core[END_REF] showed that the substitutional alloy is more stable than the interstitial one at 360 GPa and 6500 K, confirming that the temperature effect does not significantly vary the relative stability (at least at inner core pressures). They also showed that the inclusion of C as a dimer is more stable than both substitutional and interstitial alloying at 360 GPa and 0 K, but entropy does not favour the dimer and it becomes less stable than substitutional C at 6500 K [START_REF] Li | Carbon Partitioning Between the Earth's Inner and Outer Core[END_REF]. The enthalpies of C in different forms calculated at 0 K in this study are shown as a function of pressure in Figure 3.

The enthalpy calculations show that interstitial C is more stable than substitutional C below ~350 GPa. The difference in transition pressure with [START_REF] Huang | Carbon in iron phases under high pressure[END_REF] may be due to the size effect and concentration dependence. The dimer is shown to be more stable than both the substitutional and interstitial forms at all examined pressures from 60 to 400 GPa. However, on qualitative grounds, considering the temperature effect and entropy cost in forming the dimer, this could become unfavourable at high temperatures. Quantitative high-temperature calculations will be necessary to determine at what temperature the dimer structure becomes unstable, but, in any case, it definitely has its own stability field over a large pressure range extending to a certain temperature. Overall, these theoretical studies show that the interstitial alloy is more stable at low pressures, while the substitutional defect is more stable at 360 GPa. This also implies that a change of C incorporation mechanism is likely under inner-core conditions.

The calculated volumes of the alloys considered, together with those that we measured experimentally, at different pressures are displayed in Figure 4.

Substitutional C decreases the volume of the hcp unit cell relative to that of pure Fe, while interstitial C increases it. The dimer also increases the volume of hcp Fe but to a much lesser degree than the interstitial defect. Calculated and experimentally determined volumes are in excellent agreement at pressures above ~60 GPa, suggesting that C in the experimental sample is interstitial.

At lower pressures (below 60 GPa) the calculated volumes are smaller than those determined experimentally. This difference is attributed to the approximations in the theoretical approach.

DFT may fail to capture the correct electron correlation at low pressures, becoming more reliable at high pressures as a consequence of the increased delocalization of the electrons This mismatch between calculations and experiments at low pressures has also been observed for hcp Fe (Figure 5).

Actually, in the case of pure Fe the difference in vlumes is much more pronounced, with the experimental and computational results starting to agree closely only at ~200 GPa.

To better understand the effect of C and the incorporation mechanism, Figure 6 shows the pressure dependence of the volume difference ΔV=V-V Fe (where V is the unit-cell volume of the alloy and V Fe is the reference unit-cell volume for hcp Fe) for hcp Fe and Fe-Si-C alloys. The experimental volumes for Fe 93 C 4 Si 3 are compared to the different experimental V Fe available in the literature, whereas the calculated volumes are compared to hcp Fe from our DFT simulations.

Regardless of the reference V Fe used, we can unambiguously confirm that C in the experimental sample is not substitutional as the ΔV for substitutional C is negative. In the pressure range of the experiments (50-150 GPa), the experimental volume differences are close, although slightly higher, to those reported for the dimer configuration, whereas at pressures above 200 GPa, all (except when using the EoS of [START_REF] Miozzi | A new reference for the thermal equation of state of iron[END_REF] approach the ΔV expected for interstitial C. Once considering the different EoS of hcp Fe and the differences among the experiments (pressure scales, hydrostaticity, etc) and possible problems with the calculations at low pressures, this analysis seems to suggest that C in the experimental sample is interstitial. The more stable dimer might be present in some domains at lower pressures. Possibly the sample synthesis process at ambient pressure and high temperature favored interstitial carbon. Dimer formation in the quenched samples might then be entropically unfavorable or prevented by kinetic barriers and slow diffusion of vacancies. Previous computational studies [START_REF] Huang | Carbon in iron phases under high pressure[END_REF][START_REF] Caracas | The influence of carbon on the seismic properties of solid iron[END_REF][START_REF] Li | Carbon Partitioning Between the Earth's Inner and Outer Core[END_REF] suggest that C preferentially occupies the octahedral interstitial site at low pressure. The presence of C slightly expands the octahedron containing the impurity, resulting in an increase in unit cell volume (Figure 2 and4). According to [START_REF] Caracas | The influence of carbon on the seismic properties of solid iron[END_REF], the volume increase is almost linear with the amount of C and is approximately independent of pressure within the stability field of the hcp phase (see their Figure 2). In agreement with [START_REF] Caracas | The influence of carbon on the seismic properties of solid iron[END_REF], the volume increase is almost linear with C content, from 1.4 at% [START_REF] Yang | Effect of Carbon on the Volume of Solid Iron at High Pressure: Implications for Carbon Substitution in Iron Structures and Carbon Content in the Earth's Inner Core[END_REF]) to 4 at% C (this study) and almost independent of pressure. Contrarily to the other datasets, Fe-6.1 at% C by [START_REF] Yang | Effect of Carbon on the Volume of Solid Iron at High Pressure: Implications for Carbon Substitution in Iron Structures and Carbon Content in the Earth's Inner Core[END_REF] seems not to follow the same trend; above ~80 GPa the volume ratio decreases as a function of pressure and is much lower than that expected on the basis of the other data for lower C contents, and below that calculated by [START_REF] Caracas | The influence of carbon on the seismic properties of solid iron[END_REF] for a similar composition (Figure S2). These data by [START_REF] Yang | Effect of Carbon on the Volume of Solid Iron at High Pressure: Implications for Carbon Substitution in Iron Structures and Carbon Content in the Earth's Inner Core[END_REF] also disagree with the EoS predicted for a similar composition (Fe-6 at% C) by [START_REF] Huang | Carbon in iron phases under high pressure[END_REF] and, therefore, it appears that Fe-6.1at% C by [START_REF] Yang | Effect of Carbon on the Volume of Solid Iron at High Pressure: Implications for Carbon Substitution in Iron Structures and Carbon Content in the Earth's Inner Core[END_REF] may not be an interstitial solid solution.

To summarise, even a small amount of C in the alloy has a strong effect on the volume of hcp Fe.

Calculated and experimentally determined volumes are in excellent agreement at pressures above ~60 GPa, suggesting that C in the experimental sample is interstitial. While substitutional C decreases the volume of the hcp unit cell, interstitial C increases it. The inclusion of C as a dimer, predicted by calculations to be the most stable, also increases the volume of hcp Fe, but to a lesser degree than the interstitial defect.

Effect of the incorporation mechanism of carbon on the density

The effect of C on the density of hcp Fe is critical to estimation of the effect of C on the physical properties of the Earth's inner core. A common practice in mineral physics is to calculate the density of the alloy assuming a substitutional solid solution, which is valid, for example, for an Fe-Si alloys as Si only enters in the hcp structure as a substitutional impurity. In this way, the density of a substitutional alloy 𝐀 𝟏-𝒙 𝐁 𝒙 𝐒 (s denotes substitutional) is calculated as follows:

𝝆 = 𝑵 • [(𝟏 -𝒙) • 𝐌 𝐀 + 𝒙 • 𝐌 𝐁 ] 𝑽 ⁄ (eq. 1)
where N is the number of atoms in unit cell, V is the volume of unit cell. 𝐌 𝐀 , 𝐌 𝐁 are the molar mass for A, B and x is the molar fraction.

However, this is not directly applicable for other alloys like Fe-C, as C can have other solid solution mechanisms, such as incorporation interstitially. Following the established approach for fcc Fe-C alloys reported for example in [START_REF] Waseda | X-Ray Diffraction Crystallography[END_REF], the density of an interstitial solid solution 𝐀 𝟏-𝒙 𝐁 𝒙 𝐢 , (i denotes interstitial) in this case, is calculated as follows:

𝝆 = 𝑵 • [𝐌 𝐀 + 𝒙 𝟏-𝒙 𝐌 𝐁 ] 𝑽 ⁄ (eq. 2)
where the number N of atoms in unit cell in hcp Fe is two while in fcc is four.

The density in the case of the dimer defect, with two C atoms sitting at one lattice site, can be computed as a combination of one substitutional and one interstitial C.

In this study, the density of Fe-C-Si alloy considering Si as substitutional and C as interstitial 𝐀 𝟏-𝒙-𝒚 𝐁 𝒙 𝐒 𝐂 𝒚 𝐈 , is calculated as follow:

𝝆 = 𝑵 • [( 𝟏-𝒙-𝒚 𝟏-𝒚 ) • 𝐌 𝐀 + 𝒙 𝟏-𝒚 • 𝐌 𝐁 + 𝒚 𝟏-𝒚 𝐌 𝐂 ] 𝑽 ⁄ (eq. 3)
while if both Si and C are substitutional, the density for 𝐀 𝟏-𝒙-𝒚 𝐁 𝒙 𝐒 𝐂 𝒚 𝐒 is simply:

𝝆 = 𝑵 • [(𝟏 -𝒙 -𝒚) • 𝐌 𝐀 + 𝒙 • 𝐌 𝐁 + 𝐲 • 𝐌 𝐂 ] 𝑽 ⁄ (eq. 4)
Further details of density calculations of the Fe-C-Si alloy considering Si as substitutional and C both as interstitial and as substitutional are reported in the supplementary material. The derived density evolution for hcp Fe-C-Si with pressure is shown in Figure 7.

The density of the alloy experimentally determined in this study, with 4 at% C, is ~6.5 % smaller than that of pure Fe when C is assumed to be substitutional, whereas it is only is ~2.5 % smaller when C is assumed to be interstitial, roughly independent of pressure. Density calculations assuming the interstitial incorporation mechanism are in remarkable agreement with DFT calculations also showing a density reduction with added interstitial C in the Fe structure [START_REF] Caracas | The influence of carbon on the seismic properties of solid iron[END_REF][START_REF] Li | The elastic properties of hcp-Fe alloys under the conditions of the Earth's inner core[END_REF]. [START_REF] Caracas | The influence of carbon on the seismic properties of solid iron[END_REF] indicated that 1 wt.% (~4.5 at%) C would decrease hcp Fe density by approximately 2.8%. [START_REF] Yang | Effect of Carbon on the Volume of Solid Iron at High Pressure: Implications for Carbon Substitution in Iron Structures and Carbon Content in the Earth's Inner Core[END_REF] reported the density of Fe-1.37 wt % C (~6.1 at% C) to be 5.6% smaller than the of pure hcp Fe. These authors state that interstitial C in the hcp structure enhanced the density difference between pure Fe and the alloy. However, although the authors clearly assume C as interstitial, the densities they report (see their Figure 2) were calculated assuming substitutional solid solution (see our Figure 7). When the densities are calculated starting from the reported volumes using the equation for interstitial carbon, C, the resulting values are anomalously large, even larger than those of pure Fe (see Figure 7). As already discussed, it appears that the Fe-1.37 wt. % C studied by [START_REF] Yang | Effect of Carbon on the Volume of Solid Iron at High Pressure: Implications for Carbon Substitution in Iron Structures and Carbon Content in the Earth's Inner Core[END_REF] might not be an interstitial solid solution.

The correct interpretation of the incorporation mechanism is essential for evaluating the identity and amount of light element(s) in the Earth's core. As reported above, the incorporation mechanism affects both the volume and the mass of the unit cell. Substitutional C lowers the volume of the hcp unit cell while the interstitial and, to a lesser degree, the dimer mechanisms increase it (Figure 4 and 6). At the same time, the interstitial mechanism increases the mass while the substitutional mechanism reduces it. The two effects act to cancel each other out when calculating density. In fact, as shown by our DFT simulations, in which both the changes in unitcell volume and unit-cell contents are taken into account in an internally consistent way, the density difference between a substitutional, interstitial and dimer alloy is almost negligible (Figure 8). This is different in the case of the experiments, in which the volume is directly measured, without any prior assumption of C distribution. However, assumptions of the effective mass of the unit-cell have to be made to derive density from the measured volumes. As illustrated in Figure 8, the way this is done might lead to very different densities with very dissimilar conclusions concerning the amount of light elements required to account for PREM values. In particular, if the volume measured in the experiments is due to an interstitial mechanism but is interpreted as being substitutional, the resulting density will be significantly lower, and a much smaller concentration of light element would be required to match PREM. These results indicate that it is vital to know the exact incorporation mechanism of the alloy and its effect on the elastic properties with changing pressure and temperature, to correctly interpret the properties of the inner core. To this end, the combination of experiments and theoretical studies provide unique capabilities.

We have carried out volume and density measurements on an hcp Fe-C-Si alloy with 4 at% C and 3 at% Si up to ∼150 GPa, using a combination of synchrotron powder diffraction and DFT simulations, carried out at the same pressure conditions, to facilitate interpretation of the experiments and to validate the theoretical approach. In particular hcp Fe-C-Si alloys with a different incorporation mechanism -substitutional, interstitial, and dimerwere examined to precisely address the effect of C on the physical properties of the alloy.

Enthalpy calculations suggest that interstitial C is more stable than substitutional C below ~350

GPa. While the dimer, with two C atoms replacing an Fe atom, is more stable than both the interstitial and substitutional C, it is destabilized by temperature and, at ambient temperature, it can be kinetically inhibited by a diffusion barrier. Indeed, we find that C in the experimentally investigated hcp Fe-C-Si alloy sample takes the interstitial form, as supported by the very good agreement at pressure above ~60 GPa between the experimentally determined compression curve and that obtained from calculations for the interstitial alloy. Calculated values for density in the literature are often determined assuming a substitutional solid solution. However, this is not applicable for an interstitial-type solid solution, such as that reported here. Assuming the wrong incorporation mechanism when estimating density from measured volumes can lead to incorrect density determination, resulting in a more enhanced density contrast between the alloy and pure Fe. In turn, this can affect density-velocities systematics (Birch's law) resulting in erroneous interpretation of inner core composition. Finally, the agreement between our experiments and calculations, supports previous work [START_REF] Li | The elastic properties of hcp-Fe alloys under the conditions of the Earth's inner core[END_REF] suggesting Si and C combined can account fore core densities and velocities. [START_REF] Edmund | Velocity-Density Systematics of Fe-5wt%Si: Constraints on Si Content in the Earth's Inner Core[END_REF][START_REF] Yang | Effect of Carbon on the Volume of Solid Iron at High Pressure: Implications for Carbon Substitution in Iron Structures and Carbon Content in the Earth's Inner Core[END_REF] and reference EoS for pure Fe [START_REF] Dewaele | Quasihydrostatic equation of state of Iron above 2 Mbar[END_REF][START_REF] Sakai | Equation of state of pure iron and Fe0.9Ni0.1 alloy up to 3 Mbar[END_REF][START_REF] Yang | Effect of Carbon on the Volume of Solid Iron at High Pressure: Implications for Carbon Substitution in Iron Structures and Carbon Content in the Earth's Inner Core[END_REF][START_REF] Miozzi | A new reference for the thermal equation of state of iron[END_REF] [START_REF] Dewaele | Quasihydrostatic equation of state of Iron above 2 Mbar[END_REF][START_REF] Sakai | Equation of state of pure iron and Fe0.9Ni0.1 alloy up to 3 Mbar[END_REF][START_REF] Yang | Effect of Carbon on the Volume of Solid Iron at High Pressure: Implications for Carbon Substitution in Iron Structures and Carbon Content in the Earth's Inner Core[END_REF][START_REF] Miozzi | A new reference for the thermal equation of state of iron[END_REF] S1 andS2, reporting unit cell parameters and density of hcp Fe-Si-C from experiments and calculations, respectively; Table S3, summarizing the results from enthalpy calculations; Table S4 (uploaded separately) reporting the resulting stress matrices.

Text S1.

As reported in the main text, we denote the alloy with both substitutional and interstitial solute elements as ! !!!!! ! ! ! ! ! ! , where A is the solid solvent, B and C are the substitutional and interstitial solutes, respectively. x and y are the molar fraction. The density of the alloy can be calculated by

! = ! • ! -! -! ! -! • ! ! + ! ! -! • ! ! + ! ! -! ! ! !
where N is the number of atoms in unit cell, and it is two for the hcp structure. V is the volume of unit cell. ! ! , ! ! and ! ! are the molar mass for A, B and C respectively.

In the case of Fe-4 at% C and 3 at% Si investigated in this study, the density at P= 48.8 GPa, V = 18.816 Å 3 (see Table S1) is calculated:

! = ! • ! -!.
!" -!. !" ! -!. !" • !!. !"# + !. !" ! -!. !" • !". !"# + !. !" ! -!. !" !". !"" !!. !!" ! = 9.792 !/!" !

Where ! = 18.816×10 !!" ×6.022×10 !" = 11.331 !" !

The density calculations for other alloys can be derived from the equation above. If both B and C are substitutional, the density for ! !!!!! ! ! ! ! ! ! is then:

! = ! • ! -! -! • ! ! + ! • ! ! + ! • ! ! !
In the case of Fe-4 at% C and 3 at% Si investigated in this study, the density at P= 48.8 GPa, V = 18.816 Å 3 (see Table S1) is calculated:

3 ! = ! • ! -!.
!" -!. !" • !!. !"# + !. !" • !". !"# + !. !" • !". !"" !!. !!" ! = 9.400 !/!" ! If both B and C are interstitial, the density for

! !!!!! ! ! ! ! ! ! is then ! = ! • ! ! + ! ! -! -! ! ! + ! ! -! -! ! ! !
The dimer defect with two C atoms sitting at one lattice site is actually also a combination of one substitutional and one interstitial C. P(GPa) V(Å 3 ) a (Å) c(Å) ρ(g/cm 3 ) int ρ(g/cm 3 ) sub. ρ(g/cm 3 ) dimer 48.6( 4 Table S3. Enthalpy calculations for the different alloys

Figure Captions Figure 1 .

 Captions1 Figure Captions

Figure 2 .

 2 Figure 2. a) Unit-cell volume of Fe-C-Si alloy as a function of pressure. Error bars are smaller

  are shown for comparison. Inset: confidence ellipses centred on each data point, representing the 68.3% confidence level. The reported 1 error bars for K 0 and K 0 ʹ were derived from least-squares refinements. b) Unit-cell volumes of Fe-C alloys showing effect of C in the unit cell volume of hcp Fe.

Figure 3 .

 3 Figure 3. Calculated enthalpy of carbon (referenced to the isolated atom) in the form of

Figure 4 .

 4 Figure 4. Calculated and measured unit-cell volumes of Fe-C-Si alloys with pressure showing

Figure 5 .

 5 Figure 5. Comparison between calculations and experiments for hcp Fe and Fe-C-Si alloys as a

Figure 6 .

 6 Figure 6. Volume variation, ΔV=V-V Fe (where V is the unit-cell volume of the alloy and V Fe is

  were taken from the literature. The calculated volumes of Fe 62 C 2 Si 2 (C interstitial), Fe 61 C 2 Si 2 (dimer) and Fe 60 C 2 Si 2 (substitutional) are normalized to hcp Fe from our DFT simulations.

Figure 7 .

 7 Figure 7. Density evolution with pressure for different hcp Fe-C-Si alloys. Solid lines represent

Figure 8 .

 8 Figure 8. Density variation with pressure for different hcp Fe-C-Si alloys calculated at 300 K

  Figure 1.

Figure 2 .

 2 Figure 2.

  Figure 3.

Figure 8 .

 8 Figure 8.

Figure S1 .Figure S2 .

 S1S2 Figure S1. Illustration of the different configurations of C in hcp Fe. (a) substitutional; (b) octahedral interstitial site; (c) 2 C atoms occupying one atomic site (dimer defect)

Table 1 .

 1 EoS parameters of hcp Fe alloys

	V 0 (Å 3 )	K 0 (GPa)	K'	composition	Reference
				experiments
	22.7(2)	148(18)	5.8(4)	Fe-4 at% C-3 at% Si	This study
	22.524(62)	172.4(6.0) 4.64(14)	Fe-8 at% Si	Edmund et al. (2019)
	22.5(1)	170(10)	4.7(2)	Fe-8 at% Si	Edmund et al. (2019) (our fit)
	22.37(4)	168.9(4.8) 5.18(14)	Fe-1.4 at% C	Yang et al. (2019)
	22.40(6)	164(8)	5.4(3)	Fe-1.4 at% C	Yang et al. (2019) (our fit)
	22.37(6)	182.2(4.2) 4.68(17)	Fe-6.1 at% C	Yang et al. (2019)
	22.44(5)	174(6)	4.88(18)	Fe-6.1 at% C	Yang et al. (2019) (our fit)
	22.26(6)	169.7(5.2) 5.19(16)	Fe	Yang et al. (2019)
	22.32(5)	163(5)	5.36(5)	Fe	Yang et al. (2019) (our fit)
	22.468(24)	165(fixed)	4.97(4)	Fe	Dewaele et al. (2006)
	22.31(8)	178(7)	4.75(12)	Fe	Dewaele et al. (2006) (our fit)
	22.80(2)	129(6)	6.2(2)	Fe	Miozzi, Matas, et al. (2020)
	22.80(7)	130(6)	6.2(2)	Fe	Miozzi, Matas, et al. (2020) (our fit)
	22.468 (fixed) 155.3(2.2) 5.37(14)	Fe	Sakai et al. (2014)
	22.5(2)	151(15)	5.5(4)	Fe	Sakai et al. (2014) (our fit)
				calculations
	20.40	281	4.4	Fe 60 C 2 s Si 2	This study
	20.85	305	4.3	Fe 62 C 2 i Si 2	This study
	20.39	313	4.3	Fe 61 C 2 d Si 2	This study
	21.46	291	4.34	Fe-6 at% C i	Huang et al. (2005)
	20.29	312	4.3	Fe	This study
	s substitutional C, i interstitial C, d dimer C

Table S2 .

 S2 Volumes and densities of hcp Fe and Fe-C-Si alloys from DFT calculations.

	)		)			)						a)	m 3 )
		17.61	10.52	59.98	17.53			60.01		18.0	10.18		17.69
	60	7	8	3	5	10.153	1		53	3	60.33	7	10.224
		17.02	10.89	80.00	16.87			80.00		17.4	10.53		17.10
	80	9 16.52 2 P (GPa)	1 11.22 5	3 99.98 7 Fe 62 C 2	7 16.39 5 Enthalpy (eV/atom) 10.549 10.859 Fe 64 C 2	1 100.0 1 Fe 63 C 2	42 16.9 19	9 10.86 5	79.702 100.254	4 16.59 5	10.578 10.903
		16.08 60		-3.68	15.93	-5.39		119.9 -5.85		16.4	11.16		16.15
		6 15.69 80	11.53 11.81	120 -3.05	7 15.53	11.171 -4.59	9 140.0 -5.22		63 16.0	6 11.44	119.975	1 15.76	11.481
		8	5	140 -2.48	7	11.459 -3.77	1 -4.53		60	7	140.389	0	11.737
		15.34	12.08	159.9 -1.92	15.18	-3.14		-4.10		15.6	11.71		15.41
		9 15.03	4 12.33	8 179.9 -1.20	3 14.86	11.726 -2.39	160 -3.40		98 15.3	1	159.567	5 15.09	11.987
		2	8	9 -0.79	2	11.98 -1.77	180 -2.93		71	11.96	179.979	5	12.205
		14.74 3	12.58	200 -0.27	14.57 0	12.22 -1.14	200 -2.26		15.0 72	12.19 7	197.562	14.82 5	12.87
		14.00	13.24	259.9 0.37	13.82	-0.47		260.0 -1.85		14.3	12.84		14.05
		8	1	9 2.00	8	12.875 1.36	1 -0.35		09	7	260.014	8	11.203
		13.59 6	13.64 1	300.0 1 3.01	13.41 6	13.271 2.77	299.9 9 0.93		13.8 85	13.24	299.51	13.62 6	13.279
		13.23	14.01	3.81	13.05	3.90		339.9 1.66		13.5	13.60		13.27
		1	8	340 4.18	3	13.639 4.24	9 2.13		09	8	339.157	8	13.627
				5.11		5.27		3.11
		)	18.816(11)	2.388(3)	3.812(15)	9.792	9.400	9.592
													8

Acknowledgments

The authors thank Jeroen Jacobs (ESRF) for his assistance with the diamond anvil cell preparation. M.G.P. received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 796755. This project has received funding from NERC (grant NE/M015181/1 to L.V.) and from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (Grant agreement No. 72460 to D.A.). Data will be available at the repository of University of Padova (http://researchdata.cab.unipd.it/385/).

Contents of this file

Text S1 Figure S1 to S2 Tables S1 to S3 Tables S4 (File uploaded separately)

V (Å 3 ) ρ (g/cm 3 ) P (GP V (Å 3 ) ρ (g/c


F Errigo

L Chédot

F Morel

P Venet

A Sari

A Hijazi

R A Peña

Assessment of Aging and Performance Degradation of Supercapacitors Integrated into a Modular Multilevel Converter

The interest for modular multilevel converter (MMC) with energy storage systems (ESSs) has increased quickly over the last decades. Since an MMC has several hundreds of sub-modules (SMs) and ESS performances usually deteriorate quickly, the ability of the converter to provide the expected service during its lifetime must be investigated. In this work, supercapacitors are used as ESS and the reduction of their capacity is first calculated through lifetime simulation of the ESS using aging models. Then, simulations are used to analyze the available energy of the proposed converter considering dispersions for the initial parameters of ESSs and for the aging rate. Finally, the results are compared with different levels of redundancy and maintenance intervals to discuss the viability of the solution. For the depicted case study, results show that oversizing the ESSs without carrying any maintenance operations during the lifetime of the converter minimizes the total number of cells required.

Introduction

In the upcoming decades, most of the fossil-fuel power plants are planned to be replaced by renewable energy sources connected to the grid through power electronics. This development is particularly critical for the stability of power systems because the transmission grids were not developed to integrate resources with uncontrolled production, located far from the center of consumptions. One effect is the expected decrease of system inertia, which traditionally depends on the kinetic energy stored in the synchronous rotating machines. This implies that can result in dangerously low system frequencies which could lead further to a blackout [START_REF] Tielens | The relevance of inertia in power systems[END_REF]. This transition must inevitably be accompanied by strengthening the electrical grid.

Grid-scale Energy Storage Systems (ESSs) can be a promising solution to provide fast operating reserves. Because the required overall energy for these new short ancillary services is low compared to traditional bulk storage applications, storage technologies as supercapacitors have recently gained interests to meet these specifications [START_REF] Alvarez | Modular multilevel converter with short-time power intensive electrical energy storage capability[END_REF]. Usually, utility-scale ESSs are considered to be connected to the grid through dedicated converters. Recent works have studied the opportunity to integrate ESS within converter used for power transfer capabilities such as the Modular Multilevel Converter (MMC) for HVDC (High-Voltage Direct Current) transmission in order to provide ancillary services and, in this way, to extend the range of services they can offer [START_REF] Judge | Modular multilevel converter with partially rated energy storage with intended applications in frequency support and ancillary service provision[END_REF]- [START_REF] Zeng | A New Hybrid Modular Multilevel Converter with Integrated Energy Storage[END_REF].

Even if the energy storage units can be partially distributed among Sub-Modules (SMs) [START_REF] Errigo | Assessment of the Impact of Split Storage within Modular Multilevel Converter[END_REF], in this work, each SM of an MMC arm has its capacitor connected to a string of supercapacitors by means of a DC-DC converter to make an "energy storage sub-module" (ES-SM) (cf. Figure 1). Energy storage technologies are usually based on low-voltage elements while the voltage of a SM capacitor is of the order of one or few thousands of volts. Consequently, the number of components must be inherently high to sustain voltage as well as fulfill the minimum energy requirement.

In addition, the life expectancy of HVDC converters is in the range of decades. As it is often stated a 40-year converter operation is usually expected [START_REF] Wylie | Analysis of the effects from constant random and wear-out failures of sub-modules within a modular multi-level converter with varying maintenance periods[END_REF]. This time span is large with regard to the lifetime of electrochemical energy storage technologies. Thus, it is highly challenging for them to maintain high levels of availability throughout this period. In addition, an MMC must remain fully operational between two scheduled maintenance operations to avoid important penalties. Because a 1-GW MMC can have approximately 2400 SMs [START_REF] Peralta | Detailed and averaged models for a 401level MMC-HVDC system[END_REF], with a supercapacitor stack, the availability over time of such a solution is a key performance indicator to be considered. Furthermore, the Failure In Time (FIT), 1 FIT equaling 1 failure per billion hours, of the main components in a SM (IGBTs, capacitors, etc.) varies from hundred to few hundreds [START_REF] Wylie | Analysis of the effects from constant random and wear-out failures of sub-modules within a modular multi-level converter with varying maintenance periods[END_REF]. Thus, it can be assumed that their failure probability is low during the lifetime of a supercapacitor cell and they are not considered in the following study.

Most of the time, a capacitance decrease or an increase of the equivalent series resistance (ESR) [9], caused by the influence of temperature and voltage, is the main root of supercapacitors performance degradation leading to noncompliance with the expected service. Wear out with age can be a good representation of the outage of supercapacitors. Moreover, the characteristics of these components show natural dispersions. Such dispersion leads to unbalanced aging and uneven energy distribution among cells that can make the system prematurely unable to provide the expected service. This paper proposes a straightforward method that has never been done before for the assessment of the ability of an MMC with embedded supercapacitors to provide a given service over time. Thanks to physical models and accurate simulation, the lifetime of the storage unit of an ES-SM under a real mission profile is obtained. The inherent probabilistic properties of classical reliability model are avoided. Afterwards, Monte Carlo simulation is done to cope with the uncertainty of the manufacturing process of supercapacitors and the high number of cells to be considered due to the high modularity of an MMC. At the end, the aim of the proposed method is to evaluate the effectiveness of different maintenance intervals and redundancy schemes over the lifetime expectancy of the converter.

The content of the paper is organized as follows. The first section presents the main design parameters of an MMC with ESS and the case study. Secondly, the method to assess the available energy is described. Finally, the last section compares the obtained results, such as the number of supercapacitors required according to different maintenance policies during the lifetime of an HVDC project, and draws the conclusion of this work.

Modular multilevel converter with energy storage system for fast frequency response

The purpose of this part is to present the case study and the parameters for the multi-physical modelling and simulation related to the aging of the ESS. Facing future challenges of lacking inertia have introduced the need for short-term dynamic ancillary services from hundreds of milliseconds to tens of seconds, such as fast frequency support. It has been shown that ESS integrated into MMCs are well suited to bring this contribution since they can act quickly and since these service provisions require only a small part of the bulk power flow of a conventional unit, usually less than 50 MW [START_REF] Mosca | Mitigation of frequency stability issues in low inertia power systems using synchronous compensators and battery energy storage systems[END_REF].

Recently, some Transmission System Operators (TSO) have launched their own market or products because they required these new fast services from a grid point of view [START_REF]Enhanced frequency response[END_REF]- [START_REF] Eirgrid | DS3 System Services Implementation Project -Regulated Arrangements[END_REF]. By definition, they are used for the purpose to mitigate uncommon system frequency fluctuation in case of disturbance and one can ask of the interest to estimate the availability of such system which are rarely used under those circumstances. However, some specific products deal with the purpose to dynamically manage the normal second-by-second frequency within a tolerance band such as in the United Kingdom (UK) where the National Grid Electricity Transmission (NGET) recently introduced the "enhanced frequency response" (EFR) [START_REF]Enhanced frequency response[END_REF] or in the Nordic system with the standard Normal Frequency Containment Reserves (FCR-N) [START_REF] Entso-E | Nordic Balancing Philosophy[END_REF]. The meaning of these services is to provide a one-second real-time proportional response to deviations in the grid frequency to maintain the frequency in the normal range or to provide additional power during a short time until primary reserves reach their maximal set-point or secondary reserves start in case of severe disturbances.

According to the call for tender and public grid frequency data in the UK [START_REF]Enhanced frequency response[END_REF] [16], a one-year load profile mission for the ESS is developed. It is assumed, that the ESS has a maximum power of 50 MW with a storage capacity of 900 MWh. It has to provide a proportional response, in either direction, to a change in system frequency with a maximum support duration of 30 seconds when a range, set at 50 Hz plus or minus 0.05 Hz, is exceed. The output of the ESS can be varied from ± 9% within the frequency deadband to manage its State of Charge (SoC). Figure 2 shows these key features.

Then, a simulation of an MMC ES-SM with a storage function coupled with a multi-physical modelling of the supercapacitor stack has been implemented under MATLAB/Simulink [17][18]. It is carried out throughout the lifetime of the ESS with the previous load profile, which is repeated for every year. The current that will flow through the ESS can be deduced from the simulation to feed the aging model and track ESS lifetime. The obtained data is used to assess the performance of the proposed solution. The main parameters of an ESS per ES-SM and the converter specifications are given in Table I. Note hereinafter in this paper, an energy of 1 p.u. is the energy available in an arm of an MMC for the expected service. 460 cells per SM are needed to provide this energy at the beginning of the lifetime of the station. Due to aging, this level of energy cannot be sustained during all the operating time of the converter and a goal of this work is to determine the number of cells, 𝑁 𝑐𝑒𝑙𝑙 , to always ensure this energy content.

Energy storage remaining lifetime

Proposed method

In reliability studies, the bathtub curve (cf.

Figure 3) is widely used to estimate failure rate and thus the lack of ability of an equipment to perform its function [START_REF] Venet | Amélioration de la sûreté de fonctionnement des dispositifs de stockage d'énergie[END_REF]. By neglecting, early failures, due to manufacturing process, the main root of failures can be divided into two parts. Random failures which occurs throughout the useful lifetime of a product, with approximately a constant failure rate 𝜆 0 , and wear out failures that increase with age.

However, supercapacitors are electronic components whose FIT is negligible when the failure rate is constant (equal to 𝜆 0 ). Only wear out "failures" can be considered. In fact, applied to supercapacitors, the main failure criteria that are generally considered are a 20% capacitance loss or a 100% increase of its equivalent series resistance [9]. It can be clearly seen that these criteria influence primarily the usable energy and does not prevent the ESS to keep working but in a degraded manner. In that case, it does not conduct to an instantaneous loss of ability of the component. Meanwhile, reliability predictions aim at calculating the rate at which a component will fail, based on an estimation of 𝜆 0 from reliability standard, using probabilistic models for failure rate such as exponential or Weibull laws. One of the goals of this paper is to provide an estimation of the degradation of supercapacitor stack's performances integrated within an ES-SM of an MMC through damage analyses. Therefore, aging speed is first studied through electro-thermal simulation to determine ESS lifetime and then statistical analysis are carried out.

System modelling

In this purpose, the proposed method models the real behavior of components, by applying a specific mission profile and thus the real constraints, to observe the evolution of health indicators such as capacitance. Similarly, in this paper, the assumption that all the supercapacitors within a stack of one ES-SM has the same characteristics is made. Thus, they can be modelled as one equivalent capacitance. Eyring's law can be used to evaluate supercapacitor lifetime [START_REF] Kreczanik | Study of Supercapacitor Aging and Lifetime Estimation According to Voltage, Temperature, and RMS Current[END_REF]. It is considered the main End of Life (EoL) criteria corresponds to a doubling of the initial equivalent series resistance or if the capacitance evolves until 80% of its initial value. However, this model is valid when a supercapacitor undergoes the same constraints all along its lifetime. This is the reason why, in this work, the modelling of a supercapacitor stack of one ES-SM is based on a dynamic lifetime model from [START_REF] Kreczanik | Study of Supercapacitor Aging and Lifetime Estimation According to Voltage, Temperature, and RMS Current[END_REF]. This latter takes into account that constraints are variable over time and lead to a non-constant degradation rate of the component and so the estimation of the lifetime can change. It considers not only voltage and Figure 3 Bathtub curve temperature as the main parameters of influence but also current. The dynamic lifetime of a device 𝜏 𝑑 is estimated according to the cell voltage, its temperature and the RMS current as shown in equation ( 1):

𝜏 𝑑 (𝑡) = 𝜏 𝑜 1 𝑇 ∫ exp ( 𝑣(𝑡) 𝑣 𝑜 + 𝜃(𝑡) 𝜃 𝑜 + 𝑖 𝑟𝑚𝑠 𝑖 𝑟𝑚𝑠 𝑜 ) 𝑡 𝑡-𝑇 𝑑𝑡 (1)
where 𝑣(𝑡) and 𝜃(𝑡) are respectively the instantaneous value of the cell voltage (V) and temperature (°C) while 𝑖 𝑟𝑚𝑠 is the value of the RMS current (A) during the time window considered. 𝜏 𝑜 is the theoretical lifetime for a cell at 0 V and 0 °C while 𝑣 0 , 𝜃 0 and 𝑖 𝑟𝑚𝑠 0 are constant whose any decrease of 10 °C, 0.2 V or 30 A double the lifetime of the component. These parameters are given as follows [START_REF] Kreczanik | Study of Supercapacitor Aging and Lifetime Estimation According to Voltage, Temperature, and RMS Current[END_REF]:

𝜏 0 = 439 500 𝑦𝑒𝑎𝑟𝑠 𝜃 0 = 10 ln(2) 𝑣 0 = 0.2 ln(2) 𝑖 𝑟𝑚𝑠 0 = 30 ln(2)
During the simulation, the electrical parameters are provided by the electrical model of the ESS whereas the thermal effects are supplied through a thermal model that analyzes heat transfers and thermal behavior of a supercapacitor cell with the ambient temperature [START_REF] Hijazi | Thermal Network Model of Supercapacitors Stack[END_REF]. As reminder, in the literature, the lifetime of a supercapacitor is defined as the time required for a 20% drop of the capacitance. It does not imply a failure of the component. In this work, this is only used to obtain a model of the evolution of the capacitance over time and we do not refrain from using them beyond this duration.

Simulation results

In Figure 4, the obtained results are presented. It shows the evolution of the dynamic lifetime 𝜏 𝑑 over time. At the beginning of the simulation (t=0), when the ESS experiments a calendar aging, without current constraints, the lifetime is evaluated at 8.45 years. Then, it is clearly observed that estimation is almost constant, with only a slight decrease, when the mission profile is applied. This implies that the degradation rate of the component is not affected by the variation of the thermal and electrical parameters. Under those circumstances, supercapacitors undergo a calendar aging. Consequently, 𝜏 𝑑 can be assimilated as a static lifetime and thus the time at the equivalent capacitance of the supercapacitor stack of one ES-SM decreases by 20%, or its ESR doubled, can be determined. It is when the dynamic lifetime estimation equals the simulation time (~ 8.36 years) as described in Figure 4.

In addition, during the simulation, the temperature of the energy storage device remains the same regardless the variation in power. This is because the grid frequency excursions are most of the time within the intended deadband during one year as seen in Figure 2.(c). It means that the power that has to be supplied is relatively low (cf. Figure 2.(a)) and as a consequence the current that flows through the ES-SM is small with regard to the nominal current that the supercapacitor cell can supply (or absorb) at full power without increasing the device temperature beyond a given value, in that case 40°C [9]. This hypothesis is confirmed during the simulation as shown in Figure 5.(b) where the current is always under the nominal current (1 p.u.). In the same manner, the voltage of the ES-SM is always kept close to the voltage reference. This voltage is calculated to make use of of the maximal stored energy and handle the same amount of available energy for charging as discharging as seen in Figure 5.(a). Throughout this paper, the energy ratio is defined as the available energy from an ES-SM for the expected service and 1 p.u. is the energy content for an arm according to the initial sizing (cf. Table I). Based on the previous results, it can be assumed in that situation an ES-SM undergoes calendar aging. Thus, the evolution of the capacitance of an ES-SM over time can be supposed to follow a square root dependency law characteristic of a floating aging that traduced the formation of a solid layer at the interface electrode/electrolyte [START_REF] German | Ageing law for supercapacitors floating ageing[END_REF], as shown in equation ( 5) with 𝐴 𝑐 a proportional coefficient.

𝐶 𝑠𝑐 (𝑡) = 𝐶 𝑠𝑐0 -𝐴 𝑐 √𝑡 [START_REF] Zeng | A New Hybrid Modular Multilevel Converter with Integrated Energy Storage[END_REF] Finally, according to the knowledge of the time at which the capacitance drops by 20% thanks to the simulation, it is easy to extrapolate the value of this coefficient with equation (5).

Monte Carlo simulation for reliability evaluation

Probabilistic method

Generally, reliability modelling uses statistical techniques based on reliability standards such as US Military Handbook 217 F. Thus, probabilistic analysis is made assuming a constant failure rate or by introducing a time varying failure rate to define the expected lifetime of a device [START_REF] Wylie | Analysis of the effects from constant random and wear-out failures of sub-modules within a modular multi-level converter with varying maintenance periods[END_REF][START_REF] Wang | Reliability Model of MMC Considering Periodic Preventive Maintenance[END_REF]. These techniques are often criticized because of the differences between field results and predictions [START_REF] Venet | Amélioration de la sûreté de fonctionnement des dispositifs de stockage d'énergie[END_REF]. Moreover, these mathematical models can quickly become complex when it models periodic maintenance and the provision of redundancy, where SMs are at different levels of aging. In this paper, the wear out of components under a load profile is used to calculate performance degradation before carrying out well-known Monte Carlo methods to an MMC.

The above-mentioned simulations focus only on one ES-SM. Therefore, it is necessary to evolve from component level to converter level since an MMC contains six arms comprising hundreds of ES-SMs. as shown in Figure 1. It is assumed that all the ES-SMs are assembled identically to describe an arm with redundancy, which means they all share load during normal operation. It implies an active redundant mode. Similarly, an ESS in a SM is just a series connection of supercapacitors. Otherwise, the total stored energy of an ES-SM can be easily estimated by observing the degradation rate of a stack and subsequently the available energy in an arm to assess its lack of ability to complete its intended function. Because in average, the six arms of an MMC behave identically, it is possible to restrict the study to only one converter arm.

Note that supercapacitors have natural dispersion in capacitance or ESR due to their manufacturing but also in reason of aging. Similarly, the kinetic of degradation of the capacity also expect discrepancies. In addition, a stack of supercapacitors of one ES-SM can be replaced during maintenance operations. As a consequence performance degradation of ES-SMs within an MMC is not uniform over time because they experience different levels of aging. To define the availability of a storage function within an MMC according to these dispersions, a sequential Monte Carlo simulation has been implemented. This flexible approach allows to consider the time progress and take into account repairable systems. By random sampling of ES-SM characteristics with a given probability distribution, a high number of independent lifetime scenarios can be generated and played. From these simulations, performance indices are measured and compared. The method is as follows.

A normal probability distribution is used to randomly select the initial capacitance assuming 95% of production samples have a value around ±20% the typical capacitance value of an ES-SM as shown in Figure 6.(a). The minimum number of supercapacitors, 𝑁 𝑐𝑒𝑙𝑙 , per ES-SM to provide the expected service is defined at the initialization stage. As shown, in Table I, one row of 460 cells of 310 F is necessary. It leads to an equivalent supercapacitors of 0.67 F. The same process is applied to the degradation speed of the capacitance (ie. slope deviation) around the average value of the coefficient 𝐴 𝑐 of (5), found with the simulation described in the previous paragraph, at the beginning of each scenario (cf. Figure 6.(b)). Then, the capacity of all ES-SMs is observed during all the life expectancy of the converter. So, the evolution of the total stored energy from the storage function can be estimated over time and thus the energy available for the service. Afterward, the availability function of the ESS of the MMC can be drawn depending on the desired energy threshold. As a reminder, in this study, the availability corresponds to the ability to provide a minimal level of energy in case of disturbance to sustain the grid. However, because aging was not considered during the design phase, it is expected that the energy requirement will not be ensured due to performance degradation. To make sure that the available energy remains above the desired threshold, redundant cells are placed in series within an ES-SM. Thus, the scenario is restarted as long as the minimum number of cells to be inserted, allowing the desired availability of the storage function, is not found. It is considered that increasing the number of cells in series allows to increase the total stored energy per stack by a factor 𝑁 𝑐𝑒𝑙𝑙 if they are all charged at its nominal voltage. Note that the topology of the interface converter and the cell configurations are not looked at for these redundant designs.

Scheduled events, such as different preventive maintenance policies, are also included in the model. As an example, at each maintenance interval 𝑇 𝑚 (years), all the stacks of supercapacitors of the ES-SMs that reach a threshold capacitance value 𝐶 𝑚𝑖𝑛 , which corresponds to a certain decrease in the rated capacitance, are replaced by a new one while the others are still operational but with a lower capacity than when they were commissioned. Note that if an ES-SM overcomes the threshold capacitance before a maintenance interval, it is still operational even if its energy content is highly reduced.

The combination of which will provide an indication of the effects from oversizing the energy storage systems within a SM and varying the maintenance interval to maintain proper operation of the energy storage function throughout all the HVDC project lifetime. To conclude, the Monte Carlo estimation is reproduced many times for each scenario and the results analyzed for the average availability of the case study.

Simulation results

The results generated from the Monte Carlo simulations are presented in this part. The estimation was reproduced 100 times for each case. Firstly, the influence of oversizing the ES-SMs to sustain different energy levels at the end of the lifetime of the converter without any maintenance operations is described. Hence, Figure 7.(a) depicts the evolution of the energy available in one arm from the ESSs over time while Table II shows that as the energy threshold increases, the number of cells per ES-SM. Obviously, the greater the initial energy available, the more redundant cells needed. It can be observed the initial sizing, that allows 1 p.u. of available energy at the commissioning, leads to less than 0.6 p.u. at the end of the lifetime of the converter. Otherwise, 818 cells per ES-SM are needed to ensure an energy level greater than 1 p.u., which implies to oversize the ESS by 80%. Finally, it should be noted in that case, each increment of 0.1 p.u. at 40 years requires a constant number of extra cells. However, this solution can face some difficulties. In fact, it is assumed supercapacitors are still used even if what is called "EoL criteria" in the literature is reached, which stands for a 20% drop in their initial capacitance. As observed in Figure 7.(a), most of the cells would lose approximately 45% of their rated value during the time span of the project which is much more than what is considered in most aging studies in the literature. Therefore, it should be clarified experimentally whether or not supercapacitors can be used with such a low capacity without leading to important concerns from a reliability point of view. In fact, the aging of supercapacitors is accompanied by a loss of performance but also with gases development. The main risk is an overpressure inside the cell that can guide to a cracking of the component or in the worst case to a blast.

In a second step, the same process is reproduced considering a maintenance interval 𝑇 𝑚 . This latter is varied from 2 to 10 years with a step time of 2 years. Indeed, an HVDC converter station scheduled outage is regularly planned for maintenance activities. All the ES-SMs, for which the capacity fade under the threshold capacitance value 𝐶 𝑚𝑖𝑛 at each maintenance operation, are replaced. In Figure 7.(b), that shows the evolution of the energy available from the energy storage units over time, it was decided to replace it when they reached a milestone of 70% of the initial capacitance. Similarly, an arm must keep a target availability for the expected service of 1 p.u. during the lifetime of the converter, whatever the maintenance interval, that correspond to the minimum energy needed as explained previously. Figure 7.(b) does not follows a similar shape as in Figure 7.(a) because an increase of energy is brought at each maintenance intervals. Thus, the availability is improved. Furthermore, the number of cells per ES-SM to keep an energy level for the expected service greater than 1 p.u. is reduced compared to the same case without maintenance activities.

Clearly in that case, the available energy is a combination of the total number of supercapacitors per ES-SMs and the number of time this stack has to be changed, which implies an optimum configuration. For this purpose, the relationship between the number of supercapacitors in a stack per ES-SM and the capacitance for which they are changed according to different maintenance intervals is presented in Figure 8.(a). Figure 8.(b) illustrates in the same manner the number of stacks that should be changed while Figure 8.(c) describes the total thousands of supercapacitors required during all the lifecycle of the converter.

Less frequent maintenance needs more cells per ES-SM but less stacks of supercapacitors to be replaced. This latter is explained by the fact that each stack in all ES-SMs are oversized with regard to the capacitance threshold. Similarly, decreasing the capacitance threshold allows less replacement. These choices lead to an important gain with less supercapacitors needed during the lifetime of the converter (cf. Figure 8.(c)), even if it is clearly necessary to oversize the system (cf. Figure 8.(a)). At the same time, a regular change of the stack of an ES-SM is not necessarily the best choice. Moreover, such maintenance activities are usually expensive particularly for offshore converters. In order to highlight this assessment and compare the different strategies, the total number of cells required to have at least 1 p.u. of available energy all the time, and so a full availability, with any maintenance is plotted in dashed line on Figure 8.(c). Undoubtedly, this strategy to regularly have maintenance suffers of the quick decrease of capacitance at the beginning of the lifetime of a supercapacitor, due to the square root dependency, that leads to often change cells. As shown in Figure 7.(a) with no maintenance operations, a 20% decrease in capacity is expected during the first 8.4 years but only another loss of 25% is noticed during the last 31.6 years. There is an obvious trade-off that concerns the capacitance threshold for which supercapacitors have to be changed against the risk of failures due to significant aging. Nonetheless, it seems that oversizing the energy storage system without carrying any maintenance operations or the bare minimum seems a pertinent solution.

Conclusion

This paper proposes a general methodology with a practical application that considers the physics of aging mechanisms to analyze performance degradation of ESSs integrated into an MMC during the lifetime of an HVDC project. A multi-physical simulation was conducted to estimate the expected lifetime of one supercapacitor stack under a high resolution mission profile. Then, a lifetime model of a supercapacitor stack to track the evolution of the capacitance as a function of time was defined. Afterwards, Monte Carlo simulations were performed to evaluate the availability of a storage service provided by an MMC under different redundancy and maintenance interval schemes.

Results show that increasing the number of planned maintenance lead to an important increase of the number of supercapacitor stacks required during the project. It was also found that oversizing the ESS at the beginning of the lifetime of the converter per ES-SM allows to minimize maintenance activities, which are expensive tasks, and to keep a good level of energy to fulfill requirements. However, it is done at a cost of using a supercapacitor beyond its normal utilization range. It will be necessary to evaluate if it does not lead to a high probability of failure. In future works also, cost assessment indicators should be equally proposed to evaluate the cost effectiveness of both strategies to define an optimum with regard to the total cost associated. Finally, the proposed methodology can potentially be extended to other energy storage technologies such as batteries.

Figure 1

 1 Figure 1 Topology of a modular multilevel converter with integrated energy storage systems

Figure 2 (

 2 Figure 2 (a) ESS power response to grid frequency deviation (b) 5 hours extract from the one year second-by-second frequency load profile from the UK system (c) Frequency duration curve of the associated mission profile

Figure 4

 4 Figure 4 Evolution of the dynamic lifetime 𝜏 𝑑

Figure 5

 5 Figure 5 Simulation results for the 1st year with: (a) Voltage at the terminal of the supercapacitor stack of one ES-SM (b) Current that flows through it

Figure 6 (

 6 Figure 6 (a) ES-SM capacitance distribution (b) Example of a random selection of five ES-SMs with slope and initial capacitance dispersion

Figure 7 (

 7 Figure 7 (a) Available arm energy for the expected service for different energy requirements at the end of life of the converter without maintenance (b) Available arm energy for the expected service for different maintenance intervals with a full availability

Figure 8 (

 8 Figure 8 (a) Number of supercapacitors in a stack per ES-SM against capacitance threshold according to different maintenance intervals (b) Number of stacks of supercapacitors replaced against capacitance threshold according to different maintenance intervals (c) Total thousands number of cells capacitance threshold according to different maintenance intervals

Table I : Grid parameters and MMC specifications for the simulation

 I Power rating, 𝑃 𝑚𝑚𝑐 1 GW DC bus voltage, 𝑣 𝑑𝑐 640 kV AC line to ground voltage, 𝑣 𝑎𝑐 192 kV Number of ES-SM per arm, 𝑁 400 Sub-module nominal voltage, 𝑣 𝑐 𝑠𝑚 1600 V Energy needed per ES-SM for the expected service, 𝐸 𝑒𝑠𝑠 375 kJ Maximum ESS power per ES-SM, 𝑃 𝑒𝑠𝑠 20.83 kW Minimum number of supercapacitors, 𝑁 𝑐𝑒𝑙𝑙 , per ES-SM 460

	Supercapacitor cell reference	Maxwell Standard Series BCAP0310
	Cell capacitance, 𝐶	310 F
	Cell Equivalent Series Resistance, 𝐸𝑆𝑅	2.2 mΩ
	Ambient temperature	55°C

Table II : Minimum number of cells per ES-SMs to ensure different energy requirements for the expected service at the end of the lifetime of the converter without maintenance

 II

			Energy requirements at 40 years
	Energy (p.u.)	Initial sizing	0.8	0.9	1
	Number of cells 𝑁 𝑐𝑒𝑙𝑙𝑠	460	655	736	818

Acknowledgements

This work was supported by a grant overseen by the French National Research Agency (ANR) as part of the "Investissements d'Avenir" Program ANE-ITE-002-01.