Quentin Bollaert

Mathieu Chassé

Hebatalla Elnaggar

Amélie Juhin

Alexandra Courtin

Laurence Galoisy

Cécile Quantin

Marius Retegan

Delphine Vantelon

Georges Calas

Bollaert Quentin

Chassé Mathieu

Elnaggar Hebatalla

Juhin Amélie

Courtin- Nomade Alexandra

Galoisy Laurence

Quantin Cécile

Retegan Marius

Vantelon Delphine

Calas Georges

Niobium speciation in minerals revealed by L2,3-edges XANES spectroscopy

Keywords: Niobium speciation in minerals revealed by L2, 3-edges XANES spectroscopy Niobium, XANES, multiplet, 10Dq, local structure, speciation

INTRODUCTION

Niobium (Nb) is a typical critical element (European Commission 2020). It is intensively used by industry to produce superalloys, superconducting magnets and catalysts [START_REF] Schulz | Niobium and Tantalum p[END_REF]). In the frame of the environmental transition, Nb-based materials are also valued for their large range of innovative properties that result from structural defects [START_REF] Ismael | A review and recent advances in solar-to-hydrogen energy conversion based on photocatalytic water splitting over doped-TiO2 nanoparticles[END_REF][START_REF] Ma | Niobium/tantalum-based materials: Synthesis and applications in electrochemical energy storage[END_REF].

This multitude of applications is expected to increase the demand in the coming years. This economic importance is combined with the scarcity of Nb deposits. The low concentration of Nb on Earth, ca. 12 ppm in the upper continental crust [START_REF] Rudnick | Composition of the Continental Crust[END_REF], and its geochemical properties limit the incorporation of Nb into most common rock-forming minerals. Under oxygen fugacities relevant to lithospheric and surficial environments, Nb occurs only as Nb 5+ [START_REF] Cartier | Evidence for Nb 2+ and Ta 3+ in silicate melts under highly reducing conditions: A XANES study[END_REF] with the electronic configuration [Kr] 4d 0 . Its high valence relative to its ionic radius (0.64 Å in octahedral coordination, i.e. [START_REF] Shannon | Revised effective ionic radii and systematic studies of interatomic distances in halides and chaleogenides[END_REF] makes Nb one of the high-field-strength elements (HFSE) along with Zr, Sn, Hf, Ta, Pb, Th and U. Niobium concentrates in silicic melts during partial melting and fractional crystallization until the late precipitation of Nb minerals of the pyrochlore and columbite supergroups [START_REF] Linnen | Granite-Related Rare-Element Deposits and Experimental Constraints on Ta-Nb-W-Sn-Zr-Hf Mineralization[END_REF]. During weathering, Nb is among the most immobile elements [START_REF] Hill | Geochemical evolution of a palaeolaterite: the Interbasaltic Formation, Northern Ireland[END_REF] due to the resistance of primary Nb minerals to alteration. It is consequently used for mass-balance calculations [START_REF] Kurtz | Refractory element mobility in volcanic soils[END_REF] and environmental studies [START_REF] Nabi | Concentrations and size distribution of TiO2 and Ag engineered particles in five wastewater treatment plants in the United States[END_REF] although the presence of organic ligands, which can increase Nb solubility, challenges this approach [START_REF] Friis | Niobium Is highly mobile as a polyoxometalate ion during natural weathering[END_REF]. Such geochemical behavior is basic to the formation of world-class lateritic Nb deposits, as a consequence of a significant residual enrichment during lateritic weathering of primary deposits, typically by an order of magnitude relative to the parent rock [START_REF] Mitchell | Primary and secondary niobium mineral deposits associated with carbonatites[END_REF].

Niobium is octahedrally coordinated with oxygen ligands in ca. 85 % of Nb minerals [START_REF] Waroquiers | Statistical analysis of coordination environments in oxides[END_REF]. The most abundant phases are multicomponent oxides with the general formula AB2O6 (columbite-group minerals: e.g. columbite, euxenite, aeschynite) and A2-site, a square antiprism, is occupied by rare-earth elements (REE), Ca, U, Th, Fe, while the octahedral B-site is occupied by Ti, Nb, Ta and W and forms double chains in a zigzag pattern by sharing at least two edges with neighboring octahedra [START_REF] Ewing | The crystal chemistry of complex niobium and tantalum oxides. IV. The metamict state: Discussion[END_REF]. Due to strong interionic repulsion, most NbO6 octahedra are distorted [START_REF] Kunz | Out-of-center distortions around octahedrally coordinated d 0 transition metals[END_REF]. Metamictization is common due to the presence of radioactive elements [START_REF] Zietlow | Thermal annealing of natural, radiation-damaged pyrochlore[END_REF]. It can form local structures made of corner-sharing octahedra at the expense of the edge-sharing linkages observed in the crystalline oxides [START_REF] Ewing | The crystal chemistry of complex niobium and tantalum oxides. IV. The metamict state: Discussion[END_REF]. Unambiguous identification of Nb minerals is needed to understand geological and geochemical processes concentrating Nb such as magmatic differentiation [START_REF] Walter | Pyrochlore as a monitor for magmatic and hydrothermal processes in carbonatites from the Kaiserstuhl volcanic complex (SW Germany)[END_REF][START_REF] Lee | Pyrochlore chemistry from the Sokli phoscorite-carbonatite complex, Finland: Implications for the genesis of phoscorite and carbonatite association[END_REF], metasomatic events [START_REF] Wu | Hydrothermal alteration of pyrochlore group minerals from the Miaoya carbonatite complex, central China and its implications for Nb mineralization[END_REF]) and supergene weathering [START_REF] Lottermoser | Compositional variation in pyrochlores from the Mt Weld carbonatite laterite, Western Australia[END_REF][START_REF] Wall | Pyrochlore from weathered carbonatite at Lueshe, Zaire[END_REF]. However, the similar crystal chemistry and metamictization of Nb minerals make the analysis of Nb speciation using X-ray diffraction analyses challenging.

Niobium is also found as a minor or trace element in minerals. Hydrothermal events can result in the formation of amorphous Fe oxides containing low amounts of Nb [START_REF] Wu | Hydrothermal alteration of pyrochlore group minerals from the Miaoya carbonatite complex, central China and its implications for Nb mineralization[END_REF]. In metasomatized mantle peridotites, rutile may host more than 95 % of the Nb budget [START_REF] Kalfoun | HFSE residence and Nb/Ta ratios in metasomatised, rutile-bearing mantle peridotites[END_REF]. Niobium can also be enriched in Fe and Ti oxides that occur as veinlets or degradation products in laterites [START_REF] Lottermoser | Compositional variation in pyrochlores from the Mt Weld carbonatite laterite, Western Australia[END_REF][START_REF] Giovannini | Mineralogy and geochemistry of laterites from the Morro dos Seis Lagos Nb (Ti, REE) deposit (Amazonas, Brazil)[END_REF]. These oxides highlight the capacity of Nb 5+ to substitute for cations such as Ti 4+ , Sn 4+ , W 6+ , Ce 4+ and Fe 3+ in various minerals [START_REF] Černý | Mineralogy of niobium and tantalum: Crystal chemical relationships, paragenetic aspects and their economic implications[END_REF][START_REF] Giovannini | Mineralogy and geochemistry of laterites from the Morro dos Seis Lagos Nb (Ti, REE) deposit (Amazonas, Brazil)[END_REF], which also hinders the determination of Nb speciation in such geological contexts.

Here, we report the first analysis of Nb L2,3-edges XANES spectra using a selected set of Nb minerals (fluorcalciopyrochlore, hydropyrochlore, niocalite, Nb perovskite, columbite) and synthetic oxides (hematite, goethite, rutile, anatase, cerianite, lueshite and Nb2O5). A few studies have been published on synthetic Nb oxides and fluorides [START_REF] Sugiura | Niobium LIIl and LII X-ray absorption-edge spectra of Nb2O5 and NH4NbF6[END_REF][START_REF] Kubouchi | Direct observation of fractional change of niobium ionic species in a solution by means of X-ray absorption fine structure spectroscopy: XAFS spectroscopy of solutions[END_REF][START_REF] Wang | Atomic layer deposition of lithium niobium oxides as potential solid-state electrolytes for lithiumion batteries[END_REF]. The objective of this work is twofold. (i) Relate spectral properties to the local structure around Nb in crystalline references. Crystal-field and ligand-field multiplet (CFM and LFM) calculations provide the first data on crystal-field parameters and Nb-O bond covalence in minerals. This CFM and LFM approach has proven useful to model Fe bioavailability and mobility in terrestrial and subsurface environments [START_REF] Bhattacharyya | Ligands representing important functional groups of natural organic matter facilitate Fe redox transformations and resulting binding environments[END_REF]) and V speciation in fresnoite [START_REF] Höche | Vanadium L2,3 XANES experiments and first-principles multielectron calculations: Impact of second-nearest neighboring cations on vanadium-bearing fresnoites[END_REF]. (ii) Investigate the processes by which Nb substitutes for Fe, Ti and Ce in Nb-substituted oxides. These materials serve as models for model major Nb-bearing oxides found in geological contexts and used in industry. Overall, this first study on Nb speciation in mineralogical contexts demonstrates how the sensitivity of Nb L2,3-edges XANES spectroscopy can be used to better characterize and understand Nb mineralogy. Deciphering Nb speciation will improve our understanding of Nb enrichment in primary and lateritic deposits and may help to diversify the range of primary sources and to upgrade the efficiency of ore processing [START_REF] Černý | Mineralogy of niobium and tantalum: Crystal chemical relationships, paragenetic aspects and their economic implications[END_REF].

MATERIALS AND METHODS

Mineral synthesis procedures

Niobium-substituted Fe oxides (1 wt% Nb). The synthesis of Nb-bearing goethite (α-FeOOH) has been adapted from the hydrothermal synthesis of Sc-bearing Fe oxides at 70 °C [START_REF] Chassé | Scandium speciation in a worldclass lateritic deposit[END_REF]. Nb-bearing goethite was converted into Nb-bearing hematite (α-Fe2O3) by dehydroxylation at 300 °C for 2 h. The solid was washed using H2SO4 in order to remove potential adsorbed Nb species and poorly-crystallized material.

Niobium-substituted Ti oxides (5 wt% Nb). The synthesis of Nb-bearing anatase (TiO2) has been carried out using a sol-gel method at moderate temperature [START_REF] Dufour | Do TiO2 nanoparticles really taste better when cooked in a microwave oven[END_REF]. Niobiumbearing rutile (5 wt%) was obtained by heating Nb-bearing anatase (5 wt%) at 700 °C for 3 hours.

Niobium-bearing anatase was also synthesized with 1 and 10 wt% Nb following this procedure.

Niobium-substituted cerianite (1 wt% Nb). Niobium-bearing cerianite was synthesized using a high temperature solid-state reaction at 1,300 °C by homogenizing CeO2, Nb2O5 and Nd2O3 powders [START_REF] Olegário | Synthesis and characterization of Fe 3+ doped cerium-praseodymium oxide pigments[END_REF].

Lueshite. This reference material was obtained by following a hydrothermal synthesis procedure at 200 °C over 24 h based on the experimental protocol of [START_REF] Gouget | High ionic conductivity in oxygen-deficient Ti-substituted sodium niobates and the key role of structural features[END_REF].

Analytical methods

Powder X-ray diffraction (XRD) analyses. The successful completion of the syntheses was verified by X-ray diffraction analysis. X-ray diffraction was performed using a PANALYTICAL X'pert Pro MPD diffractometer. Measurements were carried out in Bragg-Brentano geometry using a Co Kα anode in order to minimize the X-ray absorption of Fe in hematite and goethite. Data were recorded with an X'Celerator detector between 3 °2θ and 90 °2θ with 0.017 ° steps. The incident beam mask was fixed at 20 mm and the Soller slits at 0.04 rad.

Total counting time of measurement was 3 hours per sample. The crystal structure of synthetic compounds was checked using the International Centre for Diffraction Data (ICDD) references (PDF-2 database). No secondary niobian phase was detected on the XRD patterns of Nb-bearing Fe, Ti and Ce oxides confirming that the Nb ions were incorporated into the crystal lattice and substituted for Fe, Ce and Ti (Supplemental Figs. S1 andS2).

Single crystal XRD. Single-crystal X-ray diffraction data were acquired on natural Nb reference materials on an Agilent Diffraction Xcalibur-S diffractometer equipped with a Sapphire CCD-detector with Mo K-α1 radiation (λ = 0.71073 Å, graphite monochromator) at 293 K. Data reduction, cell refinement and space group determination were performed using the CrysAlisPro software.

Electron microprobe (EMP) analyses. Electron microprobe (EMP) analyses were used to determine the crystal chemistry of the Nb reference materials sourced from the mineralogical collections. Analyses were conducted on a CAMECA SX Five electron probe at the CAMPARIS facility (Sorbonne Université, Paris). After checking the homogeneity of the samples, two to seven spots were probed on each sample at 15 kV and 40 nA for the major elements and at 15 kV and 299 nA for minor and trace elements. Wavelength-dispersive spectroscopy analyses were performed using the following standards: albite for Na, diopside for Mg, Si and Ca, orthoclase for Al and K, MnTiO3 for Mn and Ti, hematite for Fe, BaSO4 for Ba, topaz for F, LiNbO3 for Nb, metallic Ta for Ta, SrSi for Sr, zircon for Zr, galena for Pb, monazite for Th, uraninite for U, and allanite-(Y) for Y, La, Ce, Nd.

X-ray absorption spectra acquisition. Niobium L2,3-edges XANES spectra were recorded on the LUCIA beamline of the SOLEIL synchrotron radiation facility (Gif-sur-Yvette, France) operating with a storage ring current of 450 mA and energy of 2.75 GeV [START_REF] Vantelon | The LUCIA beamline at SOLEIL[END_REF]. A double crystal Si(111) monochromator crystal was used with an energy resolution of 0.25 eV at 2400 eV [START_REF] Schaefers | KMC-1: A high resolution and high flux soft Xray beamline at BESSY[END_REF]). The monochromator was calibrated at the energy of the Nb L3-edge using Nb 2 O 5 . Due to the high Nb contents in Nb minerals, the spectra were recorded on cellulosediluted pellets before being mounted on a holder. Pure pellets were used for synthetic Nbsubstituted samples. Data were obtained between 2300 and 2550 eV at room temperature, under vacuum, in fluorescence mode. Energy steps were (2, 0.2, 1, 0.2, 1) eV for energy ranges of (2300-2350), (2350-2400), (2400-2455), (2455-2490) and (2490-2550) eV, respectively, with a 1 s integration time. Two measurements per sample were sufficient to get adequate signal-to-noise ratios and to check the reproducibility of the edge features.

Data processing and simulations

X-ray absorption data processing. Processing of the data was achieved using the Larch software [START_REF] Newville | Larch: An analysis package for XAFS and related spectroscopies[END_REF]. The MBACK normalization algorithm was applied to prevent distortion in the XANES data and preserve the branching ratio [START_REF] Weng | A method for normalization of X-ray absorption spectra[END_REF]. This method applies a single smooth background function over the entire L2,3-edges data range rather than normalizing each edge independently. The normalization function is a first order Legendre polynomial (Supplemental Fig. S3) and the edge energy chosen as the input parameter of the MBACK function implemented in Larch is 2371 eV. The spectra of Nb-substituted hematite, Nb-substituted goethite, lueshite, and niocalite were trimmed beyond 2475 eV because the S K-edge hides the high-energy features of the Nb L2-edge.

Niobium L 2,3 -edges XANES features were fitted using three pseudo-Voigt functions in order to quantitatively compare and retrieve differences in the three features between the spectra.

The continuum background was removed using an error function with fixed inflection points (2373.0 and 2467.0 eV) and large width (ca. 8 eV) to reproduce the slow onset of the electronic continuum [START_REF] Brotton | Valence state fossils in Proterozoic stromatolites by L-edge X-ray absorption spectroscopy[END_REF]. No agreement or convention exists regarding the parameters of the background functions [START_REF] Henderson | X-ray Absorption Near-Edge Structure (XANES) Spectroscopy[END_REF]) but similar procedures have been applied to Fe L2,3edges [START_REF] Bourdelle | Quantification of the ferric/ferrous iron ratio in silicates by scanning transmission X-ray microscopy at the Fe L2,3edges[END_REF] andZr L2,3-edges [START_REF] Jollivet | An enhanced resolution of the structural environment of zirconium in borosilicate glasses[END_REF]). An error function was preferred over the commonly-used arctangent because its width is related to the instrumental resolution rather than to the lifetime of the excited state, which is not known a priori. Because the best mathematic fits were obtained for a relative weight of Lorentzian to Gaussian components between 0.8 and 1.0 depending on the spectra, this parameter was fixed to 0.9 in all fits for comparative purposes.

Crystal-field and ligand-field multiplet calculations. Simulations of the Nb L2,3-edges spectra were performed using the Crispy software as implemented in the Quanty library [START_REF] Haverkort | Quanty for core level spectroscopy -excitons, resonances and band excitations in time and frequency domain[END_REF]Retegan 2019). Two approaches were tested to reproduce the spectral features. The first is a CFM calculation which takes into account the crystal-field effects acting on Nb. Ligand-field multiplet calculations offer a more comprehensive model by including the ligand-to-metal charge transfer, i.e. by mixing in the ground state the 4d 0 , 4d 1 L 1 , 4d 2 L 2 , 4d 3 L 3 configurations through a hybridization Hamiltonian (more details on this notation are provided in the Supplemental Material).

Crystal-and ligand-field multiplet calculations both require the specific input parameters.

To account for the lifetime broadening, Lorentzian functions with a full width at half-maximum of 1.66 eV and 1.87 eV for the L3 and L2-edges were used to simulate the features [START_REF] Krause | Natural widths of atomic K and L levels, Kα X-ray lines and several KLL Auger lines[END_REF]. The experimental Gaussian broadening was set to ca. 1.5 eV to best reproduce the experimental spectra. The simulations of the experimental spectra were carried out by finding the best scaling factors for the Slater-Condon integrals that best account for direct and exchange interactions between the 2p-4d orbitals (F 2 pd = 2.680 eV, G 1 pd = 1.627 eV and G 3 pd= 0.946 eV for the free ion). The reduction of these values modifies the relative intensity of the features. The values of the 4d and 2p spin-orbit coupling parameters were kept close to the atomic values of 0.124 eV and 63.97 eV. The crystal-field splitting parameter (10Dq) was chosen to best reproduce the energy difference between the two features, which corresponds to the splitting of the final state of the d orbitals. The extra parameters required for LFM calculations are the hybridization strength Vmix (Vt2g, Veg), the configuration-averaged energies of the charge transfer configuration (E(4d 1 L 1 )-E(4d 0 )) ∆, the core-hole potential Udc and the average Coulomb interaction energies of a corresponding electron pair Udd. Detailed explanations of these parameters and their influence on the spectral shapes are given in the Supplemental Material (Supplemental Fig. S7).

RESULTS

Characteristics of niobium L 2,3 -edges XANES spectra

The Nb L2,3-edges XANES spectra (Fig. 1) are characterized by a split edge feature with two peaks (A and B features) followed by a peak of lower intensity, ca. 10 eV after the edge (C feature), as reported for other Nb-bearing materials [START_REF] Sugiura | Niobium LIIl and LII X-ray absorption-edge spectra of Nb2O5 and NH4NbF6[END_REF][START_REF] Wang | Atomic layer deposition of lithium niobium oxides as potential solid-state electrolytes for lithiumion batteries[END_REF]. The energy of the L3-edge features (A3, B3, C3) varies around ca. 2372, 2376 and 2387 eV, respectively (Fig. 1a).

The energy of the L2-edge features (A2, B2, C2) varies around ca. 2466, 270 and 2481 eV, respectively (Fig. 1b). At the L3-edge, the energy position and width of the A3 feature are stable, whereas that of B3 varies among the reference materials. Niocalite, Nb-substituted hematite and Nbsubstituted goethite spectra show narrower B3 features than the other reference materials. In contrast, the wide B3 feature of the Nb perovskite spectrum shows a shoulder at high energy (B '3 at 2379'3 at -2380 eV, Fig 2b) eV, Fig 2b). The intensity ratio between the A3 and B3 features (A3/B3) varies but remains larger than 1 except in the Nb-substituted cerianite L3-edge spectrum. In all compounds, C3 is of lower intensity than A3 and B3 with limited variation in position and width. Similar trends are observed at the L2-edge.

Analysis of the niobium L 2,3 -edges XANES spectra Variation of the edge splitting. The ∆E parameter, corresponding to the energy splitting between the A and B features, ranges from ca. 2.6 to 4.0 eV at the L 3 -edge and from ca. 2.3 to 4.0 eV at the L 2 -edge (Table 2). The Nb pervoskite, lueshite, fluorcalciopyrochlore and hydropyrochlore spectra have the largest ∆E values (ranging from ca. 3.6 to 4.0 eV). Most spectra have ∆E values between ca. 3.1 and 3.5 eV except for the Nb-substituted cerianite spectrum (ca.

2.6 and 2.3 eV at L 3 -and L 2 -edges, respectively). Overall, the energy splitting at the L3-edge follows the same trend (R² = 0.94) as the energy splitting at the L2-edge (Supplemental Fig. S6a).

Variation of the intensity ratio. The intensity ratio (I

A /I B ) varies between ca. 0.8 and 2.1. It is higher at the L2-edge (1.43) than at the L3-edge (1.07). In most compounds, it ranges between ca.

1.0-1.2 and 1.4-1.6 at the L3-and L2-edge, respectively (Table 2). Intensity ratios of Nb-substituted cerianite and Nb2O5 spectra are below 1 at the L3-edge. There is a strong positive correlation of this parameter between the two edges (R 2 = 0.95, Supplemental Fig. S6b).

Variation of the bandwidth. The ratio of the full width at half maximum of feature B with respect to feature A (FWHMB/A), ranges between ca. 1.0 and 2.2 for both edges. The spectra of Nbsubstituted goethite and hematite (Fig. 2c) have the lowest FWHMB/A values (ca. 1.1 and 1.4 at the L3-and L2-edge, respectively). Most spectra exhibit intermediate FWHMB/A values (ca. 1.4-1.7 and ca. 1.6-1.9 at the L3-and L2-edge, respectively). The Nb perovskite spectrum (Fig. 2b) stands out by its high FWHMB/A value (ca. 2.1 and 2.4 at the L3-and L2-edges, respectively). Niobium-substituted anatase spectra (Fig. 2d) are also characterized by high FWHMB/A, especially at the L2-edge (ca. 1.8-1.9 eV). On average, the width ratios at the L 3 -edge (1.54) are smaller than those at the L 2edge (1.73). The FWHM B/A parameter is correlated between the two edges (R² = 0.93), after excluding the lueshite spectra, which deviate from the general trend (Supplemental Fig. S6c).

Multiplet calculations

Crystal-field multiplet calculations. The spectral features of Nb-substituted cerianite and hematite are well reproduced by considering Nb in a regular polyhedron (Oh) (Figs. 3a and3b) with crystal field values of 2.8 and 3.7 eV, respectively (Table 3). The spectra of hydropyrochlore and Nb perovskite are best modeled with a 10Dq value of 3.9 eV, which reproduces the splitting between the features (Figs. 3c and3d). Nonetheless, the B' features of low intensity at 2379-2380 eV do not appear. The good reproduction of the intensity ratios of the spectra for Nbsubstituted cerianite, hematite and hydropyrochlore results from the reduction of the Slater-Condon integrals to 50 %, an indication of the equivalent mixing between ionic and covalent character in the Nb-O bonds. The high intensity ratios of Nb perovskite spectra (Table 2) are reproduced with the reduction of the Slater-Condon integrals to 35 %.

Ligand-field multiplet calculations. We calibrated our calculations using Nb-substituted hematite, which exhibits the simplest spectral features. Charge-transfer parameters (Udc, Udd, ∆, Vmix) were chosen (Table 4) in agreement with the calculations done on Ti 4+ (3d 0 )-compounds [START_REF] Okada | Theory of core level X-ray photoemission and photoabsorption in Ti compounds[END_REF]. The intensity of the features of Nb-substituted hematite spectra is improved and the C feature at ca. 2385 eV is reproduced when charge transfers are considered (Fig. 4). We also attempted to reproduce the B' features in the spectra of hydropyrochlore and Nb perovskite by investigating the influence of the hybridization strength (Vmix) on the spectral shapes.

Using the same charge-transfer parameters, the LFM calculation fails to explain the coexistence of both B' and C features in hydropyrochlore and Nb perovskite (Fig. 5).

DISCUSSION

L 2,3 -edges XANES spectroscopy mirrors the Nb electronic environment Spectral energy splitting, interatomic distances and site symmetry. Niobium L2,3-edges correspond to transitions to empty 4d states, making them sensitive to crystal-field effects. In the frame of crystal-field theory and given that O 2-is the ligand of Nb 5+ in all of the reference materials, ∆E is correlated with the 10Dq crystal field parameter, mirroring Nb 5+ site geometry and ligand distances.

The ∆E values range between 3.3 and 4.0 eV for octahedrally-coordinated Nb 5+ . The ∆E variability among octahedral reference materials is correlated with the average interatomic distances (R 2 = 0.84), as expected from the relationship between 10Dq and R -5 (Fig. 6) in a point-charge model [START_REF] Burns | Mineralogical applications of crystal field theory[END_REF] polyhedra; i.e. [START_REF] Waring | Temperature-pressure phase relationships in niobium pentoxide[END_REF]) structures (Supplemental Fig. S2c). These structural differences may explain a lower splitting as compared to the other octahedrally-coordinated Nb reference materials.

The ∆E value of the cerianite spectrum in which Nb is in a cubic site, is 6/9-7/9 times smaller compared to octahedral reference materials (Table 2). Crystal-field splitting values are lower than expected in cubic coordination (10Dq c is 8/9 of 10Dq o ). Though not detected by XRD (Supplemental Fig. S2b), the occurrence of trace amounts of tetrahedrally-coordinated fergusonite-(Ce) could result from the low solubility of Nb in cerianite [START_REF] Kolodiazhnyi | Unusual lattice evolution and magnetochemistry of Nb doped CeO2[END_REF]. In fergusonite, the crystal-field splitting parameter of tetrahedral Nb 5+ is even lower (10Dq t is 4/9 of 10Dq o ) explaining the small ∆E value of cubic Nb-substituted cerianite spectra.

Comparison of the splitting of Nb 5+ XANES features with those of other 4d 0 elements indicates an increase of the splitting with increasing cation charge. The splitting of Zr 4+ , Nb 5+ and Mo 6+ coordinated to oxygen ligands in weakly distorted octahedra is 3.2 eV [START_REF] Galoisy | Evidence for 6-coordinated zirconium in inactive nuclear waste glasses[END_REF]), 3.5-3.9 eV (this study) and 4.5 eV [START_REF] Bare | Local site symmetry of dispersed molybdenum oxide catalysts: XANES at the Mo L2,3-edges[END_REF], respectively. This trend highlights the decrease of the ionic radius with increasing cation charge ( [6] Zr 4+ = 0.72 Å, [6] Nb 5+ = 0.64 Å,

[6]

Mo 6+ = 0.59 Å, i.e. Shannon 1976) demonstrating the high sensitivity of L2,3-edges XANES spectroscopy to the speciation of 4d 0 elements.

Spectral broadening, charge transfer and site distortion. Consideration of the C features

is required in the fitting procedure in order to avoid an overestimation of the broadening of the B features. Unlike their assignment to 2p → 5s transitions [START_REF] Sugiura | Niobium LIIl and LII X-ray absorption-edge spectra of Nb2O5 and NH4NbF6[END_REF], our LFM calculations show that this feature results from hybridization between the 4d orbitals of Nb 5+ and the 2p orbitals of O 2-, in agreement with the work of [START_REF] Okada | Theory of core level X-ray photoemission and photoabsorption in Ti compounds[END_REF] on Ti compounds.

Spectra of Nb-bearing minerals all display a high-energy shoulder of varying intensity at 2378 eV. The high degree of distortion of NbO6 octahedra in columbite and niocalite (Table 1) is not evidenced from the fit of their XANES spectra, which are characterized by B' shoulders of limited intensity (Fig. 1). Furthermore, spectra of Nb perovskite and pyrochlore display the most intense B' shoulder but these reference materials are composed of regular octahedra. The intensity of the shoulder is therefore not related to the distortion of the NbO6 octahedra. Spectra of niobian perovksite, which exhibit the most intense C feature, is also the one with the largest B' shoulder.

We therefore ascribe the B' shoulder to orbital hybridization. The failure of ligand-field multiplet calculations to reproduce both B' and C features (Fig. 5) may be due to the differences between the Ti 4+ charge-transfer parameters used here [START_REF] Okada | Theory of core level X-ray photoemission and photoabsorption in Ti compounds[END_REF] and those of Nb 5+ , which are unknown.

The broadening of the XANES features for Ti 4+ and Zr 4+ (De Groot et al. 1992; Ikeno et al. 2013; Jollivet et al. 2013) indicates a distortion of the cation site. In contrast, the broadening of the Nb L2,3-XANES features, probed by the FWHMB/A parameter, is shaped by charge-transfer effects and there is no correlation with the degree of distortion of the NbO6 octahedra. The specifically high contribution of charge-transfer effects to Nb L2,3-XANES spectra is related to the covalent character of the Nb-O bonds [START_REF] Ikeno | Ab initio charge transfer multiplet calculations on the L2,3 XANES and ELNES of 3d transition metal oxides[END_REF].

Spectra of Nb-substituted anatase and rutile display the poorest fit quality (χ² = 0.2-0.5) due to poor fitting of the wide B features. The absence of a B' shoulder, suggesting the limited effect of charge transfer on the spectra, reflects the distortion of the Nb sites (Fig. 2d). On the contrary, only the spectra of Nb-substituted hematite and goethite exhibit A and B features of similar widths, suggesting a symmetrical Nb site rather than the expected low symmetry of Fe 3+ and the high distortion of FeO6 octahedra. This hypothesis is supported by the good reproduction of the spectra of Nb-substituted hematite in Oh symmetry using LFM calculations. The similarities between the spectra of Nb-substituted hematite and goethite suggest similar Nb environments despite the different geometry of the Fe-octahedra in these minerals.

Spectral intensity variations and bond covalence. The intensity ratios (IA/IB) exhibit low standard deviations in our dataset (σ = 0.13 and 0.20 at the L3 and L2-edges, respectively). The low intensity ratio in the L3-edge spectrum of Nb-substituted cerianite (0.75) stands out from the other references (1.10 on average for octahedral references). This can be explained by approximating the intensity ratio to the number of unoccupied states, i.e., 2/3 or 3/2 in regular cubic and octahedral coordination, respectively [START_REF] Ogasawara | Relativistic cluster calculation of ligand-field multiplet effects on cation L2,3 X-ray-absorption edges of SrTiO3, NiO, and CaF2[END_REF]. The intensity ratio values retrieved from the fitting procedure are not equal to the statistical values due to ligand-field multiplet effects [START_REF] Ogasawara | Relativistic cluster calculation of ligand-field multiplet effects on cation L2,3 X-ray-absorption edges of SrTiO3, NiO, and CaF2[END_REF]). However, the intensity ratio at the L2-edge of the octahedral reference materials is close to the statistical value of 3/2 (1.46 on average) because this edge is less affected by multiplet effects [START_REF] De Groot | Differences between L3 and L2 X-ray absorption spectra of transition metal compounds[END_REF].

The intensity of the multiplet effects can be retrieved from our simulations through the reduction factor of the Slater-Condon integrals. The values of the reduction factor range from 0.35 to 0.54, indicating a major contribution from the covalent nature of the Nb-O bonds. These values are remarkably low compared to those around 0.7-0.8 observed in most oxides of 3d elements (De Groot 2005). However, the modeling of other Nb compounds, such as CsNbF 6 , requires an even greater reduction of these parameters down to 40 % [START_REF] De Groot | Differences between L3 and L2 X-ray absorption spectra of transition metal compounds[END_REF]. Therefore, the reduction of the Slater-Condon parameters mirrors the covalent nature of the bonds formed by 4d elements relative to 3d elements [START_REF] Singh | Covalency and chemical bonding in transition metal complexes: An ab initio based ligand field perspective[END_REF].

Why combine L 2 and L 3 -edge XANES spectroscopy? As the spectral features are correlated between the L 2 and L 3 -edge (Supplemental Fig. S5), the in-depth analysis of only one edge is sufficient to determine the coordination number of Nb 5+ and characterize the the local Nb environment. However, the simulation of XANES spectra using CFM and LFM calculations cannot be carried out on a single edge. The lower sensitivity of the L 2 -edge to atomic-multiplets effects involves several pairs of solutions of the10Dq value and the reduction of Slater integrals that satisfactorily reproduce the L 2 -edge spectra. Only the joint modeling of both edges constrains the reduction of the Slater integrals, because of their higher sensitivity to atomic-multiplet effects. The choice for the acquisition and analysis of one or combined edges depends on the precision of the crystallographic information one wants to extract from the spectra.

Mechanisms of Nb substitution in oxides

Niobium-substituted Ti oxides. The discrepancy between the average Ti-O distance in the host structure and the observed Nb-O distances points out specific mechanisms of substitution in Nb-substituted TiO2 oxides (Fig. 6). Estimation of the Nb-O distances in Ti oxides by interpolation of the correlation between ∆E and R -5 gives values of ca. 2.03 Å for Nb-substituted anatase and ca. 2.01 Å and for Nb-substituted rutile. Such an increase of the Nb-O distances is expected from the larger ionic radius of Nb with respect to Ti ( [6] Nb 5+ = 0.64 Å and [6] Ti 4+ = 0.605 Å, i.e. [START_REF] Shannon | Revised effective ionic radii and systematic studies of interatomic distances in halides and chaleogenides[END_REF]. Moreover, the symmetry of the Nb 5+ site in TiO2 is reduced to a lower symmetrical group as evidenced by the broadening of the B3 and B2 features (Fig. 2d and Table 2). Increasing Nb substitution leads to lower ∆E and higher FWHMB/A (Table 2), indicating an increased distortion of the Nb sites and increased Nb-O distances. Such distortion of Nb sites and an increase in Nb-O distances with increasing amounts of Nb were also evidenced using EXAFS (1.98 Å and 1.99 Å for Ti0.97Nb0.03O2 and Ti0.92Nb0.08O2, respectively, i.e. [START_REF] Gardecka | On the nature of niobium substitution in niobium doped titania thin films by AACVD and its impact on electrical and optical properties[END_REF]. The mechanism allowing the incorporation of Nb 5+ in rutile and anatase is the formation of one Ti 4+ vacancy for every four Nb 5+ [START_REF] Sheppard | Effect of niobium on the structure of titanium dioxide thin films[END_REF]) following conditions of equilibrium and electroneutrality [START_REF] Sheppard | Defect chemistry and charge transport in niobium-doped titanium dioxide[END_REF] as illustrated by the Kröger-Vink notation:

2[𝑁𝑏 2 𝑂 5 ] ⇔ 4[𝑁𝑏 𝑇𝑖 • ] + [𝑉 𝑇𝑖 ′′′′ ] + 10[𝑂 𝑂 𝑥 ](1) 4[𝑉 𝑇𝑖 ′′′′ ] = [𝑁𝑏 𝑇𝑖 • ](2)
In geological environments, the substitution of Nb 5+ for Ti 4+ in Ti oxides is also expected to be accommodated by interstitial Fe 3+ in oxidized environments following the substitution 2Ti 4+ = Fe 3+ + 2Nb 5+ [START_REF] Ruck | Geochemical association of Fe and Nb in synthetic and natural cassiterites and rutiles[END_REF]) and by Fe 2+ in reducing environments, according to 3Ti 4+ = Fe 2+ + 2Nb 5+ [START_REF] Giovannini | Mineralogy and geochemistry of the Morro dos Seis Lagos siderite carbonatite, Amazonas, Brazil[END_REF].

Recent EXAFS analyses have highlighted the splitting of the first shell into two Nb-O contributions at 1.88 and 2.01 Å in Nb-substituted anatase [START_REF] Ribeiro | Influence of Nb-doping on the local structure and thermoelectric properties of transparent TiO2:Nb thin films[END_REF]. The increase of Nb-O distances up to ca. 2.01 Å in the latter study and this work could reflect the presence of Nb2O5 clusters in Nb-substituted anatase and rutile [START_REF] Atanacio | Niobium segregation in niobium-doped titanium dioxide (rutile)[END_REF], as the values of Nb-O distances in Nb2O5 are also 2.01 Å (Table 1). Even if the formation of Nb2O5 clusters cannot be excluded, the differences between the spectra of Nb2O5 and Nb-substituted TiO2 (Table 2) indicate that Nb clustering is not an important process.

Niobium-substituted Fe oxides. In hematite and goethite, Fe 3+ has a similar ionic radius ( [6] Fe 3+ = 0.645 Å, i.e. [START_REF] Shannon | Revised effective ionic radii and systematic studies of interatomic distances in halides and chaleogenides[END_REF] to Nb 5+ ( [6] Nb 5+ = 0.64 Å). The substitution of Nb 5+ for Fe 3+

in goethite has been shown to be isomorphic [START_REF] Oliveira | Catalytic properties of goethite prepared in the presence of Nb on oxidation reactions in water: Computational and experimental studies[END_REF]. Despite the difference between the Fe sites in hematite and goethite (Table 1), the similar Nb XANES spectra also argue for an isomorphic substitution in hematite. Charge balance may be maintained by substituting 5Fe 3+ for 3Nb 5+ , which induces cation vacancies. Our study shows that Nb-substituted Fe oxides can be reproduced by CFM and LFM calculations in a regular Oh site symmetry even though the initial Fe site has a C3 symmetry (Figs. 3b and4b). The regularity of the NbO6 octahedra is also supported by the absence of broadening of features B3 and B2 (Fig. 2c). The alignment of their experimental ∆E values on the trend with respect to theoretical Fe-O distances (Fig. 6) indicates that the average Nb-O distances are not significantly modified with respect to the average Fe-O distances.

However, EXAFS data on Nb-substituted hematite indicate low values of Nb-O distances

(1.87 and 2.05 Å) compared to the Fe-O distances in hematite (1.92 and 2.11 Å, i.e. [START_REF] Silva | Nb-containing hematites Fe2-xNbxO3: The role of Nb 5+ on the reactivity in presence of the H2O2 or ultraviolet light[END_REF]). The differences between our results and those of [START_REF] Silva | Nb-containing hematites Fe2-xNbxO3: The role of Nb 5+ on the reactivity in presence of the H2O2 or ultraviolet light[END_REF] Niobium-substituted cerianite. Charge and size differences (0.97 and 0.64 Å, for Ce 4+ and Nb 5+ respectively, i.e. Shannon 1976) limit the substitution of Nb in cerianite. The substitution of Nb 5+ for Ce 4+ was explained by the displacement of Nb 5+ cations from ideal cubic symmetry that leads to oxygen vacancies [START_REF] Hiley | Pair distribution function analysis of structural disorder by Nb 5+ inclusion in ceria: Evidence for enhanced oxygen storage capacity from under-coordinated oxide[END_REF]). It has also been suggested that half of the Ce ions become trivalent to account for charge balance. The other half of the charge compensation would be accommodated by interstitial oxygens [START_REF] Kolodiazhnyi | Unusual lattice evolution and magnetochemistry of Nb doped CeO2[END_REF] according to:

[𝑁𝑏 𝐶𝑒 • ] = [𝐶𝑒 𝐶𝑒 ′ ] + 2[𝑂 𝑖 ′′ ] + [𝑂 𝑖 ′ ]
Our simulation of the Nb-substituted spectrum in Oh point symmetry (Fig. 3a) indicates that a Nb 5+ -Nd 3+ co-doping eases the charge balance equilibrium (2[Ce 4+ ] = [Nb 5+ ] + [Nd 3+ ]) and prevents significant distortion of the Nb site by limiting the formation of interstitial oxygen and vacancies.

IMPLICATIONS Deciphering the origin of the properties of technological Nb-substituted oxides

The synthetic Nb-bearing oxides selected in this study are considered important in meeting technological challenges associated with climate change mitigation or environmental decontamination (e.g., [START_REF] Ismael | A review and recent advances in solar-to-hydrogen energy conversion based on photocatalytic water splitting over doped-TiO2 nanoparticles[END_REF][START_REF] Ma | Niobium/tantalum-based materials: Synthesis and applications in electrochemical energy storage[END_REF]). Among other applications, they are used for dye-sensitized solar cells [START_REF] Lü | Improvedperformance dye-sensitized solar cells using Nb-doped TiO2 electrodes: Efficient electron injection and transfer[END_REF], Na-ion battery anodes [START_REF] Usui | Nb-doped rutile TiO2: A potential anode material for Na-ion battery[END_REF], direct solar fuel conversion [START_REF] Fu | Nb-doped hematite nanorods for efficient solar water splitting: Electronic structure evolution versus morphology alteration[END_REF]) and catalysts for the removal of exhaust-gas nitrogen oxide [START_REF] Lian | Nb-doped VOx/CeO2 catalyst for NH3-SCR of NOx at low temperatures[END_REF]. These enhanced physico-chemical properties result from the substitution of Nb 5+ for a cation, leading to oxide mobility and transport properties due to charge unbalancing, oxygen vacancies or interstitials. The properties of these Nb materials may be rationalized through a knowledge of the mechanisms of Nb insertion, which can be unveiled by Nb L2,3-edges XANES spectroscopy.

Upgrading ore processing through quantitative analysis of Nb speciation

The choice of Nb ore-processing strategies depends on the nature of the Nb minerals, their degree of alteration and their crystal-chemical variability. It is essential to correctly characterize the Nb ores to optimize ore processing [START_REF] Černý | Mineralogy of niobium and tantalum: Crystal chemical relationships, paragenetic aspects and their economic implications[END_REF]. Niobium recovery from pyrochlore using amine collectors is currently being used in the three ore processing plants worldwide [START_REF] Guimaraes | The complexity of the niobium deposits in the alkalineultramafic intrusions Catalao I and II -Brazil[END_REF]. Nevertheless, this process induces a significant loss of columbite and Nb-rutile between feed and concentrate [START_REF] Gibson | Challenges in Niobium Flotation[END_REF]. Moreover, high concentration of Fe in pyrochlore, which is a typical feature of hydropyrochlore, reduces its floatability and causes Nb loss during ore processing [START_REF] Chehreh Chelgani | Study of pyrochlore matrix composition effects on froth flotation by SEM-EDX[END_REF]. Niobian rutile has also been identified in the world-class deposit of Araxá but its significance in the overall Nb budget remains elusive [START_REF] Neumann | Comprehensive mineralogical and technological characterisation of the Araxá (SE Brazil) complex REE (Nb-P) ore, and the fate of its processing[END_REF]. Niobium XANES spectroscopy is a useful tool to identify Nb speciation in ore minerals to optimize ore-processing strategies.

Unravelling Nb geochemical behavior in natural systems

Recent findings in hydrothermal contexts [START_REF] Andersson | Mineralogy, paragenesis, and mineral chemistry of REEs in the Olserum-Djupedal REE-phosphate mineralization, SE Sweden[END_REF][START_REF] Dietzel | Hydrothermal processes in the Fen alkaline-carbonatite complex, southern Norway[END_REF]) and nearsurface conditions [START_REF] Wall | Pyrochlore from weathered carbonatite at Lueshe, Zaire[END_REF][START_REF] Chakrabarty | Rinkite, cerianite-(Ce), and hingganite-(Ce) in syenite gneisses from the Sushina Hill Complex, India: Occurrence, compositional data and petrogenetic significance[END_REF][START_REF] Giovannini | Mineralogy and geochemistry of laterites from the Morro dos Seis Lagos Nb (Ti, REE) deposit (Amazonas, Brazil)[END_REF] dynamics in supergene environments as has been done for rare-earth elements [START_REF] Janots | Ce(III) and Ce(IV) (re)distribution and fractionation in a laterite profile from Madagascar: Insights from in situ XANES spectroscopy at the Ce LIII-edge[END_REF][START_REF] Chassé | Australian laterites reveal mechanisms governing scandium dynamics in the critical zone[END_REF][START_REF] Borst | Adsorption of rare earth elements in regolith-hosted clay deposits[END_REF].

  could reflect the high amounts of Nb incorporated in the Fe oxides (5 and 10 wt% Nb) or the synthesis protocol. The incorporation of Nb in Fe oxides may result in a symmetrization of the Nb site due to the relaxation of the crystal structure around Nb arising from the induced cation vacancies.

  demonstrate the alteration of Nb minerals and question the presumed immobility of Nb (MacLean and Barrett 1993). The demonstrated sensitivity of Nb L2,3-edges XANES spectroscopy opens up the possibility to track Nb mineral speciation in weathered Nb-enriched samples. It can also be used to determine the mechanisms of alteration of primary Nb minerals and reveal the mechanisms of Nb incorporation in newly-formed phases. A combination of bulk and micro-XANES analyses can thus be used to enlighten the molecular-scale and geochemical factors that induce contrasted Nb

  . Nonetheless, Nb-substituted Ti oxides deviate from this trend, showing

	distinct Nb-O and average Ti-O distances.
	Among octahedral Nb reference materials, the spectra of Nb2O5 show the lowest ∆E values
	(3.06 and 3.28 eV at the L2-and L3-edge, respectively), slightly higher than previously reported (2.6
	and 2.8 eV; i.e Sugiura et al. 1988). Nb2O5 powders may differ from each other in the nature and
	proportion of the polymorphs they contain (Pilarek et al. 2017). Unfortunately, the material
	investigated by Sugiura was not characterized by XRD. Our Nb2O5 reference consists of several
	polymorphs in equivalent proportions, of which two crystallize in monoclinic (ca. 90 % of NbO6
	octahedra and 10 % NbO5 polyhedra: Andersson 1967) and orthorhombic (NbO5 and NbO7

ACKNOWLEDGMENTS

We are grateful to Jean-Claude Boulliard and Éloïse Gaillou for the supply of rare mineral species from the Sorbonne Université and École Nationale Supérieure des Mines de Paris mineralogy (MINES ParisTech) collections. Alain Demourgues and Guillaume Gouget are thanked for providing the niobate samples and commenting the manuscript. We thank Sophie Cassaignon and Tsou Hsi Camille Chan Chang for the synthesis of Nb-substituted anatase. We thank Benoît Baptiste, Ludovic Delbes, Michel Fialin and Nicolas Rividi for experimental support during XRD and EPMA analyses. We acknowledge SOLEIL for provision of synchrotron radiation facilities and thank the staff of the LUCIA beamline for their help in the measurement of Nb L 2, 3 -edges (Proposal No. 20191239).

CFM parameters

Nb


Camilla Fiorini
email: camilla.fiorini@lecnam.net

Pierre-Marie Boulvard

Long Li

Étienne Mémin

A two-step numerical scheme in time for surface quasi geostrophic equations under location uncertainty

In this work we consider the surface quasi-geostrophic (SQG) system under location uncertainty (LU) and propose a Milstein-type scheme for these equations, which is then used in a multi-step method. The SQG system considered here consists of one stochastic partial differential equation, which models the stochastic transport of the buoyancy, and a linear operator linking the velocity and the buoyancy. In the LU setting, the Euler-Maruyama scheme converges with weak order 1 and strong order 0.5. Our aim is to develop higher order schemes in time, based on a Milstein-type scheme in a multi-step framework. First we compared different kinds of Milstein schemes. The scheme with the best performance is then included in the two-step scheme. Finally, we show how our two-step scheme decreases the error in comparison to other multi-step schemes.

Introduction

The main aim of the modelling under location uncertainty (LU) consists in simulating on coarse meshes an enriched system mimicking a high resolution deterministic chaotic dynamics. Such LU models allow one to recover phenomena such as backscattering, dissipation and reorganisation on very coarse meshes. Furthermore, it provides a natural framework for uncertainty quantification analysis [START_REF] Resseguier | Quantifying truncationrelated uncertainties in unsteady fluid dynamics reduced order models[END_REF]. The LU framework, first introduced in [START_REF] Mémin | Fluid flow dynamics under location uncertainty[END_REF], is based on the decomposition of the Lagrangian velocity into two components: a large-scale smooth component and a small-scale fast oscillating one. This decomposition leads to a stochastic transport operator, and one can, in turn, develop the stochastic version of classical fluid-dynamics systems derived from the Navier-Stokes equations. SQG in particular consists of one stochastic partial differential equation (SPDE), which models the stochastic transport of the buoyancy, and a linear operator relating the velocity and the buoyancy:

     db t = 1 2 ∇ • (a∇b t )dt -v * • ∇b t dt -∇b t • σdB t , b t = N (-∆) 1/2 ψ, u = ∇ ⊥ ψ, (1)
where b t is the buoyancy at time t, u the large-scale smooth velocity, N a constant depending on the vertical oscillation frequency of the buoyancy and a Coriolis parameter, B a Wiener process, ψ the stream function and v * = u - 1 2 ∇ • a + σ∇ • σ is a corrected velocity associated with the effect of the noise inhomogeneity on the advected variables. The spatial correlations of the noise are given through an integral kernel operator σ (here assumed deterministic and symmetric for sake of simplicity), and the variance matrix, a, given by the matrix kernel of the operator σσ provides a local measure of the noise strength. For more details on the derivation of this system, see [START_REF] Li | Stochastic modeling and numerical simulation of ocean dynamics[END_REF][START_REF] Resseguier | New trends in ensemble forecast strategy: uncertainty quantification for coarse-grid computational fluid dynamics[END_REF]. In the rest of this work we will mainly focus on the first equation, and the last two will be condensed in u = H(b). Concerning the modelling of the noise, we use the equivalent convenient spectral definition:

σdB t = m ϕ m dβ m t ,
where β m = β m (t) are independent one-dimensional standard Brownian motions and ϕ m = [ϕ m x , ϕ m y ] T (x) are basis functions. The number of terms involved in the sum is in theory infinite, but in numerical application a truncation is considered. In the definition of the numerical schemes we will thus assume that it is a finite sum. For the computation of the basis functions, two strategies are possible: an offline strategy, where they are defined from the eigenfunctions of an empirical covariance tensor built from high-resolution data as described in [START_REF] Li | Stochastic modeling and numerical simulation of ocean dynamics[END_REF][START_REF] Resseguier | New trends in ensemble forecast strategy: uncertainty quantification for coarse-grid computational fluid dynamics[END_REF]; an online strategies where the functions are updated during the simulation and in this case they are a function of the buoyancy b. With this representation, the variance tensor reads:

a = m ϕ m (ϕ m ) T .

Numerical schemes

In this section we derive a two-step numerical scheme in time for the SQG system under LU (SQG-LU). We compare this scheme to other multi-step schemes for the SPDE, in particular the ones developed in [START_REF] Cotter | Numerically modeling stochastic lie transport in fluid dynamics[END_REF] and [START_REF] Cotter | Modelling uncertainty using stochastic transport noise in a 2-layer quasi-geostrophic model[END_REF], and show how our scheme improves the precision. Concerning discretisation in space, standard spectral methods are used: the linear terms are treated in the Fourier space, whilst the nonlinear terms are discretised in the physical space. The derivation of the time scheme consists of two steps : first, we derive a class of Milstein schemes for SQG-LU and we empirically verify their convergence, then a two-step scheme is proposed.

Derivation of a Milstein scheme

To design the Milstein schemes, we consider the integral form of the SPDE in (1), namely

b t = b t0 + t t0 1 2 ∇ • (a∇b s ) -v * • ∇b s ds - t t0 m ∇b s • ϕ m dβ m s , (2)
and we can define the following functions:

f (b t , t) = 1 2 ∇ • (a∇b t ) -v * • ∇b t and g m (b t , t) = -∇b t • ϕ m . (3)
We can now use the functional extension of the Itô formula [START_REF] Cont | A functional extension of the Itō formula[END_REF] for both f and g to write their differential forms:

f (b t , t) = f (b t0 , t 0 ) + t t0 ∂f ∂s (b s , s)ds + t t0 ∂f ∂b (b s , s)db s + 1 2 t t0 ∂ 2 f ∂b 2 (b s , s)d b, b s (4)
g m (b t , t) = g m (b t0 , t 0 ) + t t0 ∂g ∂s m (b s , s)ds + t t0 ∂g ∂b m (b s , s)db s + 1 2 t t0 ∂ 2 g ∂b 2 m (b s , s)d b, b s (5)
We remark that, since the basis ϕ m is constant in time then so is a and the functions f and g m do not depend explicitly on time, therefore ∂f /∂t = ∂g m /∂t = 0.

Concerning the first derivatives with respect to b, it has to be interpreted as a Fréchet derivative. The Fréchet derivative of an operator F is the bounded linear operator DF (x) which satisfies the following relation:

lim h →0 F (x + h) -F (x) -DF (x)h h = 0, (6)
which implies that for a linear operator DF (x)h = F (h). We start for g and use the fact that ∇ is a linear operator:

∂g ∂b (b)b = -∇b • ϕ m -∇b • ∂ϕ ∂b m . (7)
If the basis is computed offline, ϕ m does not depend on b and therefore the second term in ( 7) is zero. If the basis is computed online and ϕ m does depend on b, we can rewrite the second term of the sum by components and, using the chain rule, one has:

∇b • ∂ϕ ∂b m = ∂b ∂x ∂ϕ m x ∂b + ∂b ∂y ∂ϕ m y ∂b = ∇ • ϕ m . ( 8
)
For the second term of f , i.e. v * • ∇b, the same considerations are valid. To compute the derivative of the first term of f , we remark that it is a composition and product of three operators, two of which are linear. We can define:

F 1 (h) = 1 2 ∇ • h, F 2 (b) = a(b), F 3 (b) = ∇b. (9)
Using the chain rule and the linearity of F 1 and F 3 one has:

D F 1 F 2 (b)F 3 (b) b = DF 1 F 2 (b)F 3 (b) DF 2 (b)F 3 (b) + F 2 (b)DF 3 (b) b = F 1 F 3 (b)DF 2 (b)b + F 2 (b)F 3 (b) = 1 2 ∇ • ∂a ∂b ∇b + a∇b . (10)
Finally, with the same considerations used above, we remark that we can write (∂a/∂b)∇b = ∇ • a. Therefore:

∂f ∂b (b)b = f (b) + 1 2 ∇ • ∇ • a -∇ • v * , ∂g ∂b m (b)b = g m (b) -∇ • ϕ m . ( 11
)
As for the Itô covariation bracket, one has:

b, b t = • t0 m g m (b s , s)dβ m s , • t0 k g k (b τ , τ )dβ k τ t = t t0 m g m (b s , s) 2 ds
We now suppose to be in either one of the following cases:

the basis functions ϕ m (and therefore a) do not depend on b and

∇ • v * = 0, -the basis functions ϕ m depend on b but are such that ∇ • v * = ∇ • ∇ • a = ∇ • σ = ∇ • ϕ m = 0.
It can be noticed that the first case corresponds to a noise defined from external high-resolution data (and thus that does not depend on the solution) while the second case boils down to impose an incompressibility condition constraint on the large scale component, ∇ • u = 0, that is indeed often considered in practice with particular scaling of the noise [START_REF] Bauer | Deciphering the role of small-scale inhomogeneity on geophysical flow structuration: a stochastic approach[END_REF][START_REF] Brecht | Rotating shallow water flow under location uncertainty with a structure-preserving discretization[END_REF]. With these assumptions, we have then:

∂f ∂b = ∂ 2 f ∂b 2 = f, ∂g ∂b m = ∂ 2 g m ∂b 2 = g m . ( 12
)
We can now replace all these expressions into (4) and ( 5), and then (4) and ( 5) into [START_REF] Brecht | Rotating shallow water flow under location uncertainty with a structure-preserving discretization[END_REF]. Keeping only the terms of order one or lower, we obtain:

b t = b t0 + f (b t0 )∆t + m g m (b t0 )∆β m + t t0 s t0 m,k g m (g k (b τ ))dβ k τ dβ m s , (13)
where ∆t = t -t 0 and ∆β m = β m t -β m t0 . We define the following quantities:

G m,k := g m (g k (b t0 )), I m,k := t t0 s t0 dβ k τ dβ m s ,
then the double iterated Itô integral in ( 13) can be approximated as follows :

m,k G m,k I m,k = m,k G m,k I m,k + I k,m 2 + G m,k I m,k -I k,m 2 .
The first symmetric term can be computed analytically from Itô integration by part formulae, I m,k + I k,m = ∆β m ∆β k -δ m,k ∆t, however the second antisymmetric term (I m,k -I k,m )/2 =: A m,k t0,t cannot and it is known as the Lévy area.

Lévy area simulation In this subsection, we briefly introduce the methods we used to simulate the Lévy area. More details can be found in [START_REF] Flint | Pathwise approximation of sdes by coupling piecewise abelian rough paths[END_REF][START_REF] Foster | An optimal polynomial approximation of Brownian motion[END_REF], where these methods were proposed. The first method to simulate the Lévy area will be referred to as the weak approximation in the rest of this work: in this method, we simulate a random variable that has the same moments as the Lévy area.

The second method, which will be referred to as the conditional method, is a recursive method: the time interval (t 0 , t) is recursively split into two subintervals of the same length, and the two following relations are used:

A m,k t0,t = A m,k t0,u + A m,k u,t + 1 2 (β m u -β m t0 )(β k t -β k u ) -(β k u -β k t0 )(β m t -β m u ) (14) E[A t0,t |B t -B t0 ] = 0.
For more details on these two methods, see [START_REF] Foster | Lévy area simulation[END_REF]. Finally, we consider a third approach, where we neglect the Lévy area. We remark that this approach is exact if G m,k = G k,m , which is not the case here.

Multi-step schemes

We next propose a two-step scheme in which the Milstein method is used as the prediction step and the Euler method is adopted as the correction step, it reads:

           b * t = b t0 + f (b t0 , u t0 )∆t + m g m (b t0 )∆β m + m,k G m,k S m,k t0,t + Ãm,k t0,t u * t = H(b * t ) b t = 1 2 b t0 + 1 2 b * t + f (b * t , u * t )∆t + m g m (b * t )∆β m (15)
where S m,k t0,t := (∆β m ∆β k -δ m,k ∆t)/2 and Ãm,k t0,t is one of the approximations of the Lévy area described in the previous subsection. This scheme will be referred to as SRK2-EM (EM stands for Euler-Milstein not for Euler-Maruyama) in the rest of the paper.

In the next section, we first analyse the results of the Milstein schemes with the different Lévy area approximations in order to select the best one. Then, we compare our multi-step scheme to two other multi-step schemes developed in [START_REF] Cotter | Numerically modeling stochastic lie transport in fluid dynamics[END_REF] and [START_REF] Cotter | Modelling uncertainty using stochastic transport noise in a 2-layer quasi-geostrophic model[END_REF]. We briefly recall them here. The first one, is based on third-order Runge-Kutta scheme with a strong stability scheme (SSPRK3) [START_REF] Cotter | Numerically modeling stochastic lie transport in fluid dynamics[END_REF], is:

                       b (1) = b t0 + f s (b t0 , u t0 )∆t + m g m (b t0 )∆β m u (1) = H(b (1) ) b (2) = 3 4 b t0 + 1 4 b (1) + f s (b (1) , u (1) )∆t + m g m (b (1) )∆β m u (2) = H(b (2) ) b t = 1 3 b t0 + 2 3 b (2) + f s (b (2) , u (2) )∆t + m g m (b (2) )∆β m (16)
where f s = f -∇ • (a∇b)/2 denotes the modified drift under Stratonovich integral. The second one, relies on Euler-Heun method [START_REF] Cotter | Modelling uncertainty using stochastic transport noise in a 2-layer quasi-geostrophic model[END_REF] equally for Stratonovich integral, reads: 1) , u (1) )∆t + m g m (b (1) )∆β m (17)

         b (1) = b t0 + f s (b t0 , u t0 )∆t + m g m (b t0 )∆β m u (1) = H(b (1) ) b t = 1 2 b t0 + 1 2 b (1) + f s (b (

Numerical results

In this section we show some numerical results. First, the effect of the different approximations of the Lévy area is studied on the Milstein scheme. Then, the multi-step scheme is assessed and compared to the ones already proposed in the literature. We focus on two variations of one specific test case plotted in Figure 1: the initial condition (left) consists of two warm elliptical anticyclones on the bottom of the domain and two cold elliptical cyclones on the top. After one day under moderate noise (centre), the four structures have rotated of approximately 45 o . After one day under strong noise (right) the nonlinearity of the dynamic is more noticeable. One can find all the configuration details used for these simulations in Chapter 6 of [START_REF] Li | Stochastic modeling and numerical simulation of ocean dynamics[END_REF] for the moderate noise configuration. For the strong noise, all the basis functions ϕ m are multiplied by a factor 10. We will use the following abbreviations for the different numerical schemes -Euler: Euler-Maruyama scheme.

-Milstein-0: Milstein scheme without the Lévy area.

-Milstein-weak: Milstein scheme with the weak approximation of the Lévy area. -Milstein-cond-n: Milstein scheme with the conditional approximation of the Lévy area. Here n stands for the number of times the interval is recursively split (cf. 14).

-SRK2-EM: scheme [START_REF] Schmitz | Measure order of convergence without an exact solution, Euler vs Milstein scheme[END_REF] with Ãm,k t0,t = 0. -SSPRK3: scheme (16).

-Heun: scheme (17).

In Figures 23one can see the difference among the Euler-Maruyama scheme and all the Milstein schemes proposed. In Figure 2 we plot for each scheme for a period of 30 day the root mean squared error (RMSE), defined as:

RMSE = 1 |Ω| E b h -b 2 L 2 (Ω) 1/2 , ( 18
)
where Ω denotes the spatial domain, b h is the numerical solution of stochastic system (1), and b stands for the reference solution downsampled from a highresolution deterministic simulation (recall that the aim of the stochastic setting is to reproduce on coarse grid high-resolution deterministic simulations). The downsampling procedure consists of a first low-pass filtering performed in the Fourier domain and a subsequent subsampling operation. The expectations are estimated from 30 realizations. These results are obtained with a ∆t twice as small for the Euler scheme with respect to the other schemes. One can observe that Milstein-0 performs slightly better than the other Milstein schemes. In Figure 3, we show the rate of strong convergence γ of all the schemes discussed, under weak and strong noise. Since the exact solution is unknown, we use the following method [START_REF] Schmitz | Measure order of convergence without an exact solution, Euler vs Milstein scheme[END_REF] to estimate γ, for a sufficiently small ∆t:

γ log 2 e 1 e 2
, with e

i := E b h T, ∆t 2 i-1 -b h T, ∆t 2 i 2 L 2 (Ω) 1/2
, where b h (T, ∆t) is the numerical solution at the final time T obtained with a time step ∆t. It is important to underline that in order for this method to work, the Brownian trajectories must be fixed. We applied this method for time steps 30, 60, 120, 240, hence obtaining two estimates for γ. Is is important to remark that the value of the time steps is given in seconds and the time-scale of the studied phenomenon is of the order of one day. For reference, the CFL condition for this problem at the initial time would give a time step around 300s. The smallest time step we considered to obtain this estimate is ten times smaller than this. As one can see from Figure 3, under weak noise all the one-step schemes provide almost identical results and all the multi-step schemes are very similar.

It is hard to distinguish among the different numerical schemes proposed. In particular, for the considered span of time steps, the error of the Euler scheme under moderate noise displays a linear trend and the prevailing convergence order in this case is one. The reason of that is explained in Appendix.

Under strong noise, it is easier to see the differences among the schemes. Milstein-weak is a slight improvement on the Euler-Maruyama, but its rate of convergence is far from 1. Milstein-0 has the highest rate of convergence among all the schemes.

In conclusion, Milstein-0 seem to perform better than the other Milstein schemes. Furthermore, it is less computationally demanding. For these reasons, we built our two-step scheme based on Milstein-0.

In Figure 3 we also compare the multi-step schemes mentioned above: they all have a similar behaviour, with a rate of convergence 0.5 ≤ γ ≤ 1, but a much smaller error when compared to the one-step schemes. In particular, the two-step scheme proposed in this work (SRK2-EM in the Figures) yields the smallest error of all for this test case. The SRK2-EM schemes also yields the smallest RMSE (cf. Figure 2).

Conclusion and perspectives

The Milstein schemes analysed in this work improve the numerical results, in particular when used in a multi-step framework. The Lévy area does not seem to play a key role in these test cases, which allows us to drastically reduce the computational costs. It must be pointed out that under weak noise, all the schemes tested provide very similar results. Some ongoing and future work include the understanding of the (non) importance of the Lévy area and whether this is related to the test case, the equations, or other factors.

Appendix: convergence of Euler-Maruyama scheme under moderate noise

To study the behaviour of our system under moderate noise, we use the formalism of [START_REF] Milstein | Numerical methods in the weak sense for stochastic differential equations with small noise[END_REF]; in particular, we write our system in the following generic form: Let Y δ • be an Euler-Maruyama integration scheme for X • with integration step δ. Then we may prove in a similar fashion to theorem 4.5.4 in [START_REF] Kloeden | Numerical solution of stochastic differential equations[END_REF] that:

dX t =
1. E[X t ] 2 ≤ C, ∀t ∈ [0, T ] 2. E |X t+δ -Y δ t+δ | X t+δ = x ≤ K(x)(δ + √ √ δ + •(δ)).
Using this and the Lipschitziannity of the coefficients in (19), we may prove a result, to some extent similar to theorem 2.1 in [START_REF] Milstein | Numerical methods in the weak sense for stochastic differential equations with small noise[END_REF], namely that

E sup t0≤t≤T |X t -Y δ t | X t0 = x ≤ K (x)(δ + √ √ δ + •(δ)). (20)
In light of this estimate, we may interpret the convergence rate displayed in Figure 3 as a case where δ is not small enough when compared to so that √ √ δ does not necessarily prevail over δ which is evidenced by the linear rate of convergence.

Fig. 1 .

 1 Fig. 1. Euler-Maruyama simulation of system (1) on a 128 × 128 spatial grid.

Fig. 2 .

 2 Fig. 2. RMSE (normalised by the amplitude of buoyancy B0 = 10 -3 m/s 2 ) of different schemes during 30 days of simulation under moderate noise.

Fig. 3 .

 3 Fig. 3. Convergence of different schemes under weak and strong noise. Order 1 in dotted black, order 0.5 in dashed black.

  a(x, t)dt + b(x, t)dW t + 2 c(x, t)dt, t ∈ [0, T ] (19)with a, b, c, being jointly L 2 -measurable in (x, t), Lipschitz, bounded lineargrowth functions in x.