A Sepetys

R Guirlet

J Rosato

Y Marandet

Doppler line shape model for the W I 5d 5 6s -5d 5 6p transition at 400.9 nm in tokamaks

Keywords: Tungsten, sputtering, fusion, line-shape, Zeeman

Tungsten spectroscopy is known to be a useful tool for quantifying W fluxes in tokamaks. We aim to analyze how background plasma parameters influence the W spectral-line shape measured in the divertor and scrape-off layer and what information could be retrieved from such an analysis. To that end, a Monte Carlo model is developed to simulate the line shape of the 400.9 nm W I line emitted by physically sputtered W ions from plasma facing components. The influence of background plasma parameters (ion and electron temperature) and geometrical effects (detector positioning with respect to the eroded surface) are investigated. The effect of Zeeman splitting on the spectral line-shape is also explored. A set of requirements for the experimental set-up needed to measure the analyzed spectral line features with line peak shifts of the order of 10 -2 Å and peak widths in the range from 10 -1 Å to 10 -2 Å are proposed. These requirements are shown to be experimentally accessible, through the installation of specific high-resolution spectrometers not typically used on magnetic fusion devices.

Introduction

For a long time a favored material for plasma facing components in fusion reactors was carbon. However due to its high hydrogen retention and high erosion rate it was deemed unfit for new generation fusion reactors. This has propelled a research into application of high Z materials in fusion reactors, due to their lower erosion and capability to handle high heat loads.

Tungsten was found to have the desired thermo-mechanical properties as well as low tritium retention. For these reasons ITER will be equipped with a W divertor and an increasing number of tokamaks are already running with partial or complete W walls (ASDEX-Upgrade [START_REF] Kallenbach | Overview of ASDEX upgrade results[END_REF], JET [START_REF] Matthews | Plasma operation with an all metal first-wall: Comparison of an ITER-like wall with a carbon wall in JET[END_REF] and now WEST [START_REF] Bucalossi | The WEST project: Testing ITER divertor high heat flux component technology in a steady state tokamak environment[END_REF]). However, tungsten is a high Z material and can lead to high radiation losses and eventually to a complete loss of confinement, if its content within the plasma is not controlled [START_REF] Ashby | A study of the effect of impurity radiation from the peripheral plasma of a tokamak reactor[END_REF][START_REF] Philipps | Comparison of tokamak behaviour with tungsten and low-Z plasma facing materials[END_REF].

One aspect of W content control is the characterisation of its sources within the tokamak. The primary tungsten source is expected to be a neutral particle flux from the plasma facing components originating from their physical erosion. The local neutral W flux can be inferred spectroscopically from a W I line and the so called S/XB factor [START_REF] Brezinsek | Spectroscopic determination of inverse photon efficiencies of W atoms in the scrapeoff layer of TEXTOR[END_REF][START_REF] Harder | ELM-resolved divertor erosion in the JET ITERlike wall[END_REF]. It is also well known that parameters such as the density and temperature of the emitting population can usually be inferred from the spectral line shape [START_REF] Hans | Principles of Plasma Spectroscopy[END_REF][START_REF] Stambulchik | Spectral Line Shapes in Plasmas II[END_REF][START_REF] Rosato | The TCV team, and The MST1 team. Stark broadening of Balmer lines with low and moderate quantum number in dense divertor plasmas[END_REF][START_REF] Guirlet | Experimental and theoretical evaluation of the He II 6560.1 å line contribution to the deuterium Dα spectral line shape in Tore Supra ergodic divertor plasmas[END_REF].

It is expected that the shape of tungsten low ionization stage spectral lines can help to characterize W sources by providing information about the neutral W velocity distribution. This comes from the Doppler effect as the wavelength for the transition changes due to motion of the source atoms. The resulting Doppler line profile should in principle reflect key aspects of sputtering physics (e.g. energy and angular distribution of sputtered atoms) and be sensitive to the energy of ions responsible for sputtering. However, as tungsten atoms have a large atomic mass, the velocities of sputtered particles are low, which results in a small Doppler shift. It is thus expected that accessing this information will require specific very high-resolution spectrometers, and one of the aims of this work is to quantify these requirements. To study this wavelength change a Monte Carlo based model was developed to estimate the Doppler profiles of neutral W lines. The model includes the physical sputtering process of the target material, the Doppler shift of the photons emitted by the sputtered atoms. Obtaining line profiles further requires accounting for the Zeeman effect, which is line dependent. In this work we apply the model to the case of the most commonly observed W I line at 400.9 nm [START_REF] Brezinsek | Spectroscopic determination of inverse photon efficiencies of W atoms in the scrapeoff layer of TEXTOR[END_REF][START_REF] Beigman | Tungsten spectroscopy for the measurement of W-fluxes from plasma facing components[END_REF], since this allows to discuss the main parameters determining the line shape and the instrumental requirements to analyse these line profiles. Applying the analysis to several W I lines simultaneously would allow to check the consistency of the measurements with the calculated Zeeman profiles and identify possible parasitic line contributions to the W I lines, and is thus highly recommended. Furthermore, the analysis could also be coupled to W II line measurements in order to further strengthen the results and evidence the degree of thermalization of W II ions.

The outline of the paper is as follows: the model for the Doppler line shape is discussed in section 2, and the results are presented in section 3 for a set of relevant projectiles and plasma conditions. The Zeeman effect resulting from the confining magnetic field is discussed in section 4. The results are then summarized in section 5, where an assessment of the information potentially retrievable from a W I line shape analysis is presented, together with the experimental requirements.

W atom source and emission

The model for the sputtered ion velocity distribution to be presented in this section is schematically represented in Figure 1. The W source process is assumed to be the sputtering of W atoms by an ion flux accelerated in the magnetized sheath and hitting a planar W surface representing a divertor plate or a limiter.

We aim to estimate the sputtered tungsten atom velocity distribution arising from the incoming ion flux characteristics and sputtering properties of the W target. From the velocity distribution thus derived and geometrical features of the experimental setup we calculate the spectral line shape. The analysis discussed here is aimed at neutral W spectral line emission which can only occur at a very narrow layer within the scrape off layer (SOL) due to the millimetric W ionisation length within SOL. This allows to disregard the line-of-sight (LoS) integration effect through the plasma volume as no neutral W sources exist there. In combination with the typical centimetric size of the target surface monitored by a LoS, this also allows one to consider that the magnetic field, and thus also its angle with respect to the LoS direction, are constant over the emission zone. The use of a Monte Carlo method for this model gives a lot of freedom regarding changes of the various parameters, and choices of different scenarios.

The input variables for the model are the background plasma parameters n e , T e , T i , detector geometry (the angle of the detector with respect to the surface and the distance between them), source dimensions, incident species and their charge.

Incident particles

Under tokamak conditions various impurities are generally present, depending of the wall materials and seeded species. Therefore, in this work we use O as a representative impurity, which will always be present to some extent. When applying our results to a particular experiment, this assumption could be refined on the basis of measured impurity concentrations. For instance, in ITER, one would expect beryllium to play a significant role while in a full W machine O could be responsible for the bulk of the erosion. The determination of the dominant ionisation stages of oxygen is substantially more complicated. This would require a combination of spectroscopic investigation and plasma edge modeling to estimate the concentrations of differently charged oxygen at the SOL. Note that self sputtering by W ions is yet another particular case which can be treated with the model, depending on the situation.

The incident ion flux is assumed to be thermalised and to have a Maxwellian energy distribution, shifted by the sound speed in the parallel direction at the entrance of the sheath. At the sheath entrance the ion and electron temperatures are assumed to be equal T i = T e . It should be noted here that while the ion temperature is often not measured, it has been experimentally observed that T i can be several times higher than T e [START_REF] Kocan | Edge ion-to-electron temperature ratio in the Tore Supra tokamak[END_REF][START_REF] Huber | Spectroscopic measurements of the ion temperature profile in front of a limiter in TEXTOR-94[END_REF], especially in the far scrape-off layer (SOL). The influence of the T i /T e ratio on the generated results is thus analyzed in Section 3.

This incident ion flux is considered to be accelerated in the sheath (Fig. 1, step 1), and the typical incident ion energy can be expressed as [START_REF] Kirschner | Modelling of the material transport and layer formation in the divertor of JET: Comparison of ITER-like wall with full carbon wall conditions[END_REF]:

E im = 2T i + 3T e Z ( 1
)
Here Z is the ion charge. In most experimental cases the ion flux composition onto the plasma facing components is not measured. As these are expected to play a role in the W I line shape, the mass and charge of the incident ions are considered to be free parameters in our model.

Sputtering

All the sputtered particles in the model are generated by physical sputtering (Fig. 1, step 2). The sputtering yield is considered to follow the revised Bohdansky equation [START_REF] Garcia-Rosales | Revised formulae for sputtering data[END_REF] and its dependence on the particle incidence angle is taken into account. However, the target is considered to be a completely smooth W surface.

The sputtered particles have a Thompson energy distribution [START_REF]Atomic and Plasma Material Interaction Data for Fusion[END_REF][START_REF] Thompson | The energy spectrum of ejected atoms during the high energy sputtering of gold[END_REF]:

P (E em ) = 2E b • E em (E em + E b ) • (E i /E thr ) 2 (E i /E thr -1) 2 .
(

) 2
Here E b = 8.68 eV is the binding energy for W atoms, while E em is the energy of the emitted particles and E i is the energy of the incident ions. E thr (threshold energy) is the minimum energy that an ion must transfer to the target atom to extract it from the surface and defines the distribution peak position. The emitted particles have a maximum possible emission energy which is proportional to the incident particle energy and can be expressed as (for E i > E thr ):

E max = E b • E i E thr -1 (3)
The Thompson energy distribution for W sputtering at three different oxygen ion impact energies is presented in Figure . 2.

Emission

The energy distribution of the sputtered particles is a combination of the Thompson distribution and the Maxwellian distribution of the incident particles. The resulting distribution is presented in Figure 3 and shows features of both distributions: a steep rise at low energy from the Thompson distribution and a continuous tail (in contrast to the sharp cut off seen in Figure 2). The values of T e = 5, 10 and 20 eV were chosen as possible temperatures near the tokamak wall during normal operation [START_REF]Atomic and Plasma Material Interaction Data for Fusion[END_REF][START_REF] Sun | Study of near scrape-off layer (SOL) temperature and density gradient lengths with Thomson scattering[END_REF][START_REF] Tokar | Modelling the transport of deuterium and tritium neutral particles in a divertor plasma[END_REF]. Temperatures below T e = 5 eV are not considered since the sputtering yield becomes very low at these conditions.

The combination of the Maxwellian distribution of incident particle energy and the Thompson energy distribution of sputtered particles can generate distribution functions with relatively wide flat top areas. This effect is more noticeable for simulations at low temperatures (T i , T e ), where the cut-off of the Thompson distribution is located close to its peak.

Photon emission occurs randomly along the trajectory of sputtered atoms and within the W ionisation length (Fig. 1, steps 4 and 5). In order to reduce the computational time we consider that all the sputtered particles will get excited and undergo spontaneous decay in our specified spectral line, since the average number of excitation per ionization event is independent on the direction in which the atom is emitted. This is equivalent to setting the S/XB factor to 1. Absolute line intensities can be obtained either by using the measured S/XB value or by implementing the excitation process in the Monte Carlo process, which gives the same results if the atomic data are consistent.

It should be noted here that the model additionally relies on the facts I) that volume recombination of W + is negligible in the conditions where sputtering does occur and II) that the line considered is optically thin (which can be checked using e.g. Eq. ( 3) in Ref. [START_REF] Hey | Optical spectroscopy of tokamak plasmas[END_REF], assuming that the population of the lower state of the transition is equal to the total W density, which is a conservative assumption). Furthermore, we assume there are no charge exchange processes involving neutral W. To our knowledge, there are no data available in the literature on such processes.

Angular dependence

The angular distribution of the sputtered particles must also be taken into account. When incident particle energies are above 500 eV the angular distribution of the sputtered particles is a cosine distribution. This impact energy of 500 eV would correspond to plasma temperatures of 35 eV for particles with a charge Z = 4. However, when the impact energy is less than 500 eV, the emitted particles have a more complex angular distribution. The latter is plotted in Figure 4, for W erosion by O + ions. This is the situation expected during the inter-ELM phases in H-mode operation, with plasma temperatures close to the divertor targets below 30 eV and down to a few eV's. At these temperatures, W sputtering is dominated by impurities for which the energy transfer to heavy W atoms in the lattice is more favorable. Table 1 shows impact energies of oxygen ions of increasing charge as a function of plasma temperatures, and indicates the corresponding range of sputtering yields. With increasing charge, the sheath acceleration becomes larger and fully determines the impact energy when 3ZT e 2T i . Also, note that W self sputtering threshold for singly ionized ions becomes significant at plasma edge temperatures of 10 eV [START_REF] Eckstein | Sputtering of tungsten and molybdenum[END_REF]. It is thus fully relevant to discuss the angular distribution for low energy regime, i.e. E i < 500 eV, because significant sputtering also occurs in this regime.

It has been reported consistently in various sources that in this energy range the emitted particles have a butterfly-like spatial distribution [START_REF] Eksaeva | ERO modelling of tungsten erosion in the linear plasma device PSI-2[END_REF][START_REF] Marenkov | Angular and velocity distributions of tungsten sputtered by low energy argon ions[END_REF][START_REF] Nishijima | Sputtering properties of tungsten 'fuzzy' surfaces[END_REF]. This has to be taken into account in this study, as it will affect the velocity distribution function of W atoms along the observation direction. Assuming we have a detector looking along the normal to the surface, there will be a smaller number of emitted particles with velocity pointing towards the detector compared to a cosine distribution (because of the empty range in between the butterfly lobes). The particles with emission angles within the Table 1. Energies (in eV) of incident oxygen ions with different charges after their acceleration in the plasma sheath. Different shaded cells: sputtering yield above 10 -4 at/ion and 10 -2 at/ion butterfly lobes would have smaller velocity projections towards the line of sight and would generate photons with a smaller blue shift. Note that the situation would be different in another observation direction. Such a distribution can be represented as [START_REF] Yu | Angular distribution of atoms during the magnetron sputtering of polycrystalline targets[END_REF]:

T i , eV T e , eV O + O 2+ O 3+ O 4+ O 5+ O 6+ O 7+ O 8+
f (x) = A cos n (θ) + B cos m (θ) (4)
The exponent n increases with increasing atomic number and varies depending on the target material and incident species. The proposed estimation for the value of n is:

n = 0.3 Z E b 0.5 1 + M 1 2M 2 3 ( 5
)
where Z is the atomic number of the target material, E b the binding energy in eV, M 1 and M 2 masses for bombarding (M 1 ) and target (M 2 ) atoms. As the remaining parameters (A, B and m) were not given explicitly by the authors, we have fitted them on the angular distributions provided by the authors [START_REF] Eksaeva | ERO modelling of tungsten erosion in the linear plasma device PSI-2[END_REF][START_REF] Marenkov | Angular and velocity distributions of tungsten sputtered by low energy argon ions[END_REF][START_REF] Nishijima | Sputtering properties of tungsten 'fuzzy' surfaces[END_REF]. The polynomial fits obtained agree well with the experimental and numerical results reported by several sources [START_REF] Eksaeva | ERO modelling of tungsten erosion in the linear plasma device PSI-2[END_REF][START_REF] Nishijima | Sputtering properties of tungsten 'fuzzy' surfaces[END_REF][START_REF] Yu | Angular distribution of atoms during the magnetron sputtering of polycrystalline targets[END_REF][START_REF] Yamamura | Angular distributions of sputtered atoms for low-energy heavy ions, medium ions and light ions[END_REF].

Such an angular distribution arises because at low energies sputtering is considered to behave in the "knock off" manner. In this case all the collisions are considered to happen only at the few topmost surface layers. Following the minimum energy-loss process these particles have a preferential angular distribution which is analysed in more detail in [START_REF] Yamamura | Angular distributions of sputtered atoms for low-energy heavy ions, medium ions and light ions[END_REF][START_REF] Yamamura | Few collisions approach for threshold sputtering[END_REF].

For incident particle energies above 500 eV and up to the ˜10 keV range (e.g. during ELMs) the impact atoms can transfer enough energy to a target atom to displace other atoms within the target. This can be considered as a linear cascade of collisions. In such a case the angular distribution is close to a cosine law.

For example, the maximum incident particle energy during a type I ELM event is considered to be proportional to the electron pedestal temperature and equal to E impact,max = 4.23 • T ped e,max [START_REF] Guillemaut | Experimental estimation of tungsten impurity sputtering due to type I ELMs in JET-ITER-like wall using pedestal electron cyclotron emission and target Langmuir probe measurements[END_REF]. In the ITER case the estimated T ped e,max is about 5 keV, which means E impact,max ≈ 20 keV. The corresponding angular distribution of the sputtered atoms is thus well described by a cosine law.

On the contrary, during the inter-ELM phase in JET experiments the recorded T e by Langmuir probes did not exceed 30 eV [START_REF] Guillemaut | Experimental estimation of tungsten impurity sputtering due to type I ELMs in JET-ITER-like wall using pedestal electron cyclotron emission and target Langmuir probe measurements[END_REF]. Assuming T e = T i and Z = 1 would result in impact energies of E impact,max ≈ 150 eV. This indicates that in such a scenario even considering Z = 3 the angular distribution of the sputtered particles should still be considered non-cosine.

These geometrical considerations for the sputtered particle angles should contribute to the changes of the spectral line during inter-and intra-ELM phases of tokamak operation, on top of the effects resulting from very different incident energies. These two effects can either reinforce each other or partially cancel each other depending on the conditions and the observation direction.

The simulated lineshapes at identical plasma parameters but with different sputtered atom angular distributions are presented in Figure 5, assuming normal observation. It shows that while the angular distribution of sputtered particles does influence the expected line shape, the change is only significant at low temperature plasma conditions (inter-ELM phase). In these conditions (T e = T i = 5 eV) the FWHM of the predicted line shape reduces from 0.024 Å to 0.01 Å when considering a butterfly angular distribution for sputtered particles. This is in line with simple considerations on the projection of velocities on the observation direction. At T e = 10 eV the line shape still exhibits a dependence on the type of angular distribution but both the peak position and the line width remain practically unchanged. The profiles become almost identical at higher plasma temperatures as incident particles have higher energies and the angular distribution becomes more cosine-like shaped. Finally, it should be noted that this comparison assumes that the material surface is flat, so that the normal to the surface is identical on the whole spatial extension from which W atoms originate. In reality the surfaces have some level of roughness, which will contribute to blur the differences between cosine and butterfly distributions. However, at low energies the lobes of the butterfly distributions are quite narrow and it would take large fluctuations of the surface normal to make the effect irrelevant.

Results

Simulations were performed with different plasma temperatures, detector orientations and different incident species and charges, both in the inter-and intra-ELM phases. The Zeeman effect is dealt with in the subsequent section.

Normal detector

The predicted line shape for the most commonly observed W I line so far at 400.9 nm is presented in Figure 6. Here the x axis represents the wavelength difference with respect to the unshifted line position (0 being the unshifted 400.9 nm wavelength). Note that for the Doppler profile, only the wavelength scale is specific to this line -the shape would be identical for other lines. The results show that the spectral line should be entirely blue shifted due to the geometry of the simulation (red shifted photons would be emitted by sputtered particles with velocities away from the detector therefore directed back into the material). The width of the line also increases as expected with increasing plasma temperature and is in the range of 10 -1 -10 -2 Å while the peak shift is of the order of 10 -2 Å. For comparison, the typical thermal Doppler broadening for kT ∼ E B ∼ 10 eV is indeed of the order of 3 • 10 -2 Å ( ∆λ λ = kT mc 2 ). These values are small compared to the common spectral resolution of the typical tokamak spectrometers and indicates that the linewidth changes should be less challenging to resolve than the peak position shift.

It is evident from the simulated line shape that it depends on the plasma temperature and type of incident ion, showing noticeable differences between sputtering by O + and O 4+ except for the high temperature cases, for which the distribution tends towards a limit distribution. When the ion charge increases keeping plasma parameters identical, the sheath acceleration 3ZT e increases, and thus the impinging ion energy increases. This shift of energy range is clearly evidenced on Fig. 6: for the low Z sputtering case, corresponding to lower impact energies, the general line shape is quite asymmetric at all temperature conditions: a steep rise on the red side of the peak and a long decaying tail on the blue side. However, for the same plasma conditions but with Z = 4, the line shape approaches the (projected) Thompson distribution (Eq. 2) for lower temperatures. The peak position becomes essentially insensitive to the incident particle energy, and only the tail of the line shapes gets longer with increasing energy, an effect difficult to evidence experimentally. This is most evident when comparing the 100 eV and 1 keV cases. This analysis shows that while at high impact energies the line shape is largely insensitive to the exact plasmas conditions, it does provides a direct check on the energy distribution of sputtered particles.

Tilted detector

Subsequently, simulations were performed with the detector at an angle with respect to the surface. This allows one to mimic common experimental situations where the line of sight is tilted with respect to the surface normal, and investigate to what extent measuring the line shape along two line of sights with different inclinations would provide a check on the angular distribution of emitted atoms. Throughout this section the line of sight was placed at 60 • to the surface normal. For comparison the simulations were performed with cosine and butterfly angular distributions of the emitted particles. The line shapes obtained with such configurations at various plasma conditions are presented in Figure 7.

As expected, the red wing of the spectral line appears when the line of sight is at an angle. This arises due to the fact that in tilted configurations it is possible for the particles to have velocity projections directed away from the detector.

Additionally, it is evident that accounting for the butterfly angular distribution for sputtered particles produces wider and more asymmetrical peaks. The strong asymmetry at low plasma temperatures results from the superposition of two distributions with different widths, one being blue-shifted and emitted by the lobe aligned with the observation direction and the other being red-shifted and emitted by the other lobe. The latter is narrower due to the asymmetry of the two lobes with respect to the observation direction. With increasing incident particle energy, the width of the red shifted component decreases as the angular distribution becomes more cosine- like and the maximum achievable red shift becomes smaller. This FWHM (Full Width at Half Maximum) increase could slightly reduce the difficulty of resolving such spectral lines. This analysis clearly shows that simultaneous observations along two different directions would provide information on the angular distribution of emitted atoms.

Inter-ELM phase

During normal tokamak operation the temperature at the plasma edge has to be as low as possible to minimise the physical erosion and the related plasma contamination and plasma facing component damage. As this is the desired operational regime, the simulations for temperatures relevant for divertor operation in the inter-ELM phases were performed at low temperature values. They are presented in Figure 8. where sputtering was assumed to be only by O 4+ ions. As the sputtering yield becomes smaller at lower particle energies, the results are normalised so the maximum of the distribution is 1, to better illustrate the change in the line shape.

The results reinforce the previous conclusion that the line shape is sensitive to the ion impact energy and that the line shape width changes might be more clearly visible than the peak shift itself. However such changes would be measurable only with a very high resolution spectrometer, tracking for instance the position of the line centroid and width during a divertor temperature scan. Performing such a line shape measurement in these inter-ELM/low temperature conditions is challenging, and would require a spectrometer with resolving power λ/∆λ of approximately 400000, to resolve spectral features of 10 -2 Å needed for the T e = 5 eV case. In other words, we assume that the FWHM of the apparatus function should not be larger than the line shape width in order to be able to retrieve information. This also gives an indication about the shifts which would be measurable. At even lower temperatures the peak shift is of the order of 10 -3 Å, but sputtering is then weak and the possibility to measure the line with a reasonable signal-to-noise ratio comes into question. Spectrometers routinely used in tokamak plasma spectroscopy are, to the best of our knowledge, not able to provide the required spectral resolution (common Czerny-Turner spectrometers can for instance provide a maximum spectral resolution of about 0.02 Å, while the width of the spectral line at T e = 2 eV conditions is of the order of 0.001 Å). However, it is possible to achieve larger spectral resolution (up to R = 10 7 according to manufacturers) with for example Littrow configuration spectrometers. Littrow configuration spectrometers can reach very high spectral resolutions, but operate with very small apertures. Another possibility would be to use a scanning Fabry-Perot interferometer with appropriate filters. Such instruments have been long used to analyse astronomical plasmas [START_REF] Epinat | Des galaxies proches aux galaxies lointaines: études cinématique et dynamique[END_REF], and recently attempts to measure ion temperature by the use of interferometry in weakly magnetized plasmas were performed [START_REF] Désangles | Forçage à grande échelle d'une colonne de plasma faiblement magnétisée: influence d'une cathode émissive de grande taille[END_REF]. More recent developments are combining the advantages of Fabry-Perot cavities with those of grating spectrometers, making use of virtually-imaged phased array (VIPA) [START_REF] Shirasaki | Large angular dispersion by a virtually imaged phased array and its application to a wavelength demultiplexer[END_REF], with compact R = 10 6 spectrometers commercially available. Such measurements might in any case require optimising the throughput of the system and using long acquisition times (possibly making use of long steady state discharges in tokamaks such as WEST [START_REF] Bucalossi | The WEST project: Testing ITER divertor high heat flux component technology in a steady state tokamak environment[END_REF] or EAST), accumulating data only in between -or during -the ELMs.

Intra-ELM phase

It is known that T e is very low during the ELM events (≈ 20 -30 eV or less [START_REF] Guillemaut | Experimental estimation of tungsten impurity sputtering due to type I ELMs in JET-ITER-like wall using pedestal electron cyclotron emission and target Langmuir probe measurements[END_REF]). Since the typical incident ion energy is 2T i + 3T e Z and T i T e , the analysis of the spectral line shape could yield indirect measurements of the plasma ion temperature T i .

The calculated line shape for plasma temperatures relevant to the intra-ELM phase is presented in Figure 9. In this case it is assumed that T e = T i = 1 keV. Since our previous argument states that the influence of T e can be neglected, it can be considered that the presented lines correspond to T i values of 2.5 keV for O + and D + and 7 keV for O 4+ . It is evident that at high temperatures the line shape does not depend on the temperature and only the tail broadens (because the cut-off energy of the Thompson distribution increases with the incident particle energy). In fact, from previous results it is evident that the line width varies insignificantly with plasma temperatures above 100 eV (Figure 6). While this shows that it is not possible to infer T i during the ELM event from the line shape, it should make the profile measured during ELMs markedly different from that obtained in the inter-ELM phase. This is a good check on the sensitivity of the diagnostic, but also on the underlying assumptions of the model, by providing direct access to the energy distribution of sputtered atoms on a wide energy range (which is not possible at lower energies, for obvious energy conservation reasons).

Sensitivity to T i

The initial assumption T e = T i at the sheath entrance may often underestimate the value of T i , for which measurements are rather scarce. As mentioned earlier T i is often observed to be several times higher than T e , especially in the region far from the separatrix [START_REF] Kocan | Edge ion-to-electron temperature ratio in the Tore Supra tokamak[END_REF]. This can lead to significant sputtering even rather far from the separatrix. The comparison between the line shapes estimated by our model considering T i = T e and T i = 3T e are presented in Figure 10. This result shows that as soon as the electron temperature is larger than a few eV, the increase in T i does not change the results significantly since we already are in the high impact energy regime (above 100 eV), where only the tail of the Doppler line shape is affected. In fact, for the conditions shown on Figure 10, we have impact energies of the order of 150 eV.

At lower impact energies, the effect of the T i /T e ratio is expected to be stronger, since it contributes in determining the impact energy. However, as mentioned earlier, ions reaching the surface are considered to have an energy distribution with mean energy E = 2T i + 3T e Z. One can rewrite this energy as E = T e ( 2T i Te + 3Z), which shows that T i controls the incident energy only when Z 2 3

T i

Te . This leads to the condition Z 2 3

T i
Te ≈ 2 if we assume that T i Te = 3. It shows that as soon as Z ≥ 2, the incident particle energy is defined by the incident ion charge rather than by the T i value at the sheath entrance.

Zeeman effect

It is well known that in tokamak environment the spectral lines exhibit a Zeeman splitting due to the presence of a magnetic field. Since the aim of this study is to investigate the possibility of inferring plasma parameters from the spectral line shape, this effect can not be overlooked.

To the best of our knowledge, as of this date there is no reported in-depth analysis of the 5d 5 ( 6 S)6s 7 S 3 -5d 5 ( 6 S)6p 7 P 0 4 W I 400.9 nm spectral line Zeeman splitting, even though Zeeman patterns have been published [START_REF] Brezinsek | Spectroscopic determination of inverse photon efficiencies of W atoms in the scrapeoff layer of TEXTOR[END_REF][START_REF] Beigman | Tungsten spectroscopy for the measurement of W-fluxes from plasma facing components[END_REF]. We provide in the Appendix the calculation of the positions and relative intensities of the Zeeman components for this line. An example of the resulting Zeeman pattern is shown in Figure 11. Here the Zeeman structure is well resolved and despite some overlapping of the central and lateral components the typical triplet shape starts to appear. The validity of the weak field approximation is marginal for this field strength (B = 2.5 T) but treating the case where fine structure is of the same order (Paschen-Back effect) or negligible compared to the Zeeman splitting (strong field) requires more information about the states involved than just the spectroscopic terms, and this is outside the scope of the current work.

The convolution of the Zeeman profile (at B = 2.5 T and with a line of sight perpendicular to the direction of the B field (θ = 90 • )) with the calculated Doppler line is shown in Figure 11. Since the width of the modelled line is lower than the gap between the Zeeman components the Zeeman profile can be clearly seen. However, in the higher plasma temperature cases (T e = T i > 5 eV) the modelled line shape width becomes wider than the spacing between the individual Zeeman components and the contributions from the different components start to overlap. The expected line-shape in such a case is presented in Figure 12. The magnetic field strength at the magnetic axis of the ITER plasmas is expected to be above 5 T. A similar magnetic field is expected on the ITER plasma divertor legs. In such a case the Zeeman splitting should be wider than proposed in this study and such overlapping should be minimal (but the weak field approximation used here is no longer valid).

As a final point, we discuss the question of the natural width of the line. The latter can be calculated from the value of the Einstein coefficient for the transition, namely A ki = 1.63 × 10 7 s -1 [START_REF] Kramida | NIST Atomic Spectra Database (ver. 5.6.1)[END_REF] for the line under discussion. The full width at half maximum (FWHM) of the Lorentzian resulting from the natural broadening of the line is then given (in angular frequency) by ∆ω = A ki [START_REF] Brezinsek | Spectroscopic determination of inverse photon efficiencies of W atoms in the scrapeoff layer of TEXTOR[END_REF] This gives the following value for the FWHM in wavelength, ∆λ = λ 2 0 2πc ∆ω 1.4 × 10 -5 Å (7) where λ 0 is the wavelength of the line center, and c the speed of light. This rough calculation (neglecting the fact that the line is split in several components) shows that the natural width is small enough to be ignored in our analysis.

Summary

A Monte Carlo model for simulating the W I spectral line shape was developed. The parameters controlling the line shape of the W I spectral line have been investigated, and the possibility to extract information on these parameters from a line shape analysis has been discussed. The results also apply to other W I lines, except for the Zeeman effect which is line specific. The simulations show that the modelled Doppler line shape depends on the ion impact energy (itself a function of T e , T i and the charge of the incident particles). For impact energies above ∼ 100 eV, the line shape essentially converges to the Thompson distribution projected on the observation direction, and thus provides a check on the model assumptions instead of information on the plasma conditions. This is relevant for the intra-ELM phase, while during the inter-ELM phase a line shape analysis would help characterizing ion impact energies, and thus plasma conditions in conjunction with other diagnostics, such as probe measurements for T e , Z eff ... The angular distribution of the sputtered atoms plays an important role in the analysis of the Doppler profiles, and simultaneous measurements at different viewing angles are shown to provide strong checks on model assumptions.

In the inter-ELM phase the predicted Doppler peak shift for the 4009 Å line is of the order of 10 -2 Å. For the higher wavelengths this result is expected to be slightly higher due to the linear dependence between the shift and the central wavelength of the line. While the Doppler shift is likely to be difficult to resolve, the line shape changes are shown to be measurable using appropriate commercially available high resolution spectrometers. Such measurements are likely to require averaging over a long acquisition time, the limits being set by the duration of achievable stationary phases in discharges, and proper optimization of the optical throughput of the system. During the intra-ELM phase the predicted line shape practically does not depend on the incident particle energy. According to our model, Doppler profiles should differ significantly during and in between ELMs, which could provide a simple check on the ability to resolve the profile.

Interpretation of the Doppler line shape is likely to require an accurate Zeeman pattern model when the Doppler shift is comparable to the magnetic field and/or the fine structure splitting of the line, so that merging of different components is observed. The measurements of the lineshape should include several W I lines in order to verify whether line shape changes are consistent with Zeeman pattern calculations, and allow identifying possible parasitic line contributions. Additionally, the measurement of several spectral lines simultaneously would help choosing the optimal spectral line in terms of sensitivity and Zeeman pattern -the commonly used 400.9 nm line does have a complex Zeeman structure, while some others will have more triplet like features for instance, which could facilitate the analysis. a vanishing magnetic field, the magnetic sublevels are degenerate, with a degeneracy 2J + 1, corresponding to |M J | ≤ J. ⊥ is the polarization vector, which is in the plane perpendicular to the observation direction k, and d is the electric dipole operator. The double sum on the upper and lower states is reduced to a single sum by the selection rule ∆M J = 0, ±1, yielding the following result for the intensity of the line

I σ ± (∆ω) = 3 M j =-3 C σ M J I 0 ∆ω ± ω Z 4 (M J -7) , (A.4)
where the Larmor frequency ω Z = eB 2me , ∆ω = ω -ω 0 with ω 0 the frequency of the unperturbed line, and the relative intensity of components, after evaluation of the relevant 3j-symbols for the transitions [START_REF] Robert D Cowan | The theory of atomic structure and spectra[END_REF], is given by

C π M J = 16

Figure 1 .

 1 Figure 1. Schematic representation of the model (see text for description)

Figure 2 .

 2 Figure 2. Thompson energy distribution for W sputtering at 3 different oxygen ion impact energies.

Figure 3 .

 3 Figure 3. Sputtered particle energies, accounting for the energy distribution of impinging ions.

Figure 4 .

 4 Figure 4. Angular distribution of sputtered W particles in polar coordinates for various incident particle energies

Figure 5 .

 5 Figure 5. Simulated line-shape with different sputtered particle angular distribution at various plasma temperatures (T e = T i , sputtering by O 4+ ions)

Figure 6 .

 6 Figure 6. Calculated line shape for normal detector placement when sputtering is generated by O + and O 4+ ions at various plasma temperatures T e = T i

Figure 7 .

 7 Figure 7. Simulated line shape for detector placement at a 60 • angle to the surface normal at various plasma conditions (T e = T i , sputtering by O 4+ ions, with different sputtered particle angular distributions)

Figure 8 .

 8 Figure 8. Predicted line shape for low temperature plasma conditions, when sputtering by O 4+ ions

Figure 9 .

 9 Figure 9. Calculated line shape when sputtering by O 4+ (O V), O + (O II) and D + (D II) ions at T e = T i = 1 keV

Figure 10 .

 10 Figure 10. Calculated line shape assuming different T i /T e ratios, for O 4+ incident ions.

Figure 11 .Figure 12 .

 1112 Figure 11. Convolved Zeeman -Doppler lineshape, at T e = T i = 5 eV, for Z = 4, θ = 90 • and B = 2.5 T

Figure 13 .

 13 Figure 13. Convolved Zeeman -Doppler lineshape, at T e = T i = 1 keV, for Z = 4, θ = 90 • and B = 2.5T

2 θ 3 M

 23 × (I σ + + I σ -) (A.2)where I π and I σ± correspond to different polarisations. Explicitly, we haveI π (∆ω) =

Acknowledgments

Discussions with Olivier Peyrusse are gratefully acknowledged. This work has been carried out within the framework of the EUROfusion Consortium and French Research Federation for Fusion Studies and has received funding from the Euratom research and training programme 2014-2018 and 2019-2020 under grant agreement No 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission.

Appendix: Zeeman pattern for the 5d 5 ( 6 S)6s 7 S 3 -5d 5 ( 6 S)6p 7 P 0 4 transition at 400.9nm

In the following, we calculate the Zeeman pattern of the 5d 5 ( 6 S)6s 7 S 3 -5d 5 ( 6 S)6p 7 P 0 4 transition in the weak field approximation, namely when the Zeeman splitting is small compared to the fine structure splitting, as is the case here (the other lines of the multiplet are respectively at 407.43 nm and 428.46 nm). Here θ is the angle between the magnetic field and the observation direction, defined by a unit vector k. The intensity of the line can be written as

where |α are the upper levels and |β the lower levels of the transition, ω αβ = (E α -E β )/ and δ is the Dirac delta function. Here the upper level has a total angular momentum number J = 4, and the lower level is such that J = 3. For


F Bourgin

M-H Ramos

G Thirel

V Andréassian

Investigating the interactions between data assimilation and post-processing in hydrological ensemble forecasting

Keywords: hydrological ensemble forecasting, data assimilation, post-processing, ensemble dressing, uncertainty propagation

We investigate how data assimilation and post-processing contribute, either separately or together, to the skill of a hydrological ensemble forecasting system. Based on a large catchment set, we compare four forecasting options: without data assimilation and post-processing, without data assimilation but with post-processing, with data assimilation but without post-processing, and with both data assimilation and post-processing. Our results clearly indicate that both strategies have complementary effects. Data assimilation has mainly a very positive effect on forecast accuracy. Its impact however decreases with increasing lead time. Post-processing, by accounting specifically for hydrological uncertainty, has a very positive and longer lasting effect on forecast reliability. As a consequence, the use of both techniques is

Introduction

Addressing uncertainties in hydrological ensemble forecasting

Developing and improving operational hydrological ensemble forecasting systems is a critical step toward better decision-making and risk management.

The skill of operational hydrological ensemble forecasting systems is limited by two main sources of uncertainty [START_REF] Krzysztofowicz | Bayesian theory of probabilistic forecasting via deterministic hydrologic model[END_REF]: meteorological uncertainty and hydrological uncertainty. From a pragmatic point of view, the need to properly account for these two main sources of uncertainty arises because (i) a hydrological forecaster has no choice but to rely on uncertain meteorological forecasts; (ii) even with accurate inputs, hydrological forecasts will remain uncertain due to our limited knowledge of initial conditions and the inherent limitations of the forecast model used.

Meteorological uncertainty is commonly addressed by propagating an ensemble (or multi-scenario) input of weather forecasts. For instance, several operational and pre-operational flood forecasting systems across the globe have been set up to be forced by ensemble numerical weather predictions (see [START_REF] Cloke | Ensemble flood forecasting: A review[END_REF], for a review). Addressing the hydrological uncertainty issue is less common, although a general framework of probabilistic forecasting that includes a hydrological post-processing method has been introduced fifteen years ago by [START_REF] Krzysztofowicz | Bayesian theory of probabilistic forecasting via deterministic hydrologic model[END_REF]. Since then, a number of other hydrological uncertainty processors have been proposed [START_REF] Montanari | A stochastic approach for assessing the uncertainty of rainfall-runoff simulations[END_REF][START_REF] Montanari | Estimating the uncertainty of hydrological forecasts: A statistical approach[END_REF][START_REF] Solomatine | A novel method to estimate model uncertainty using machine learning techniques[END_REF][START_REF] Coccia | Recent developments in predictive uncertainty assessment based on the model conditional processor approach[END_REF][START_REF] Morawietz | Systematic evaluation of autoregressive error models as post-processors for a probabilistic streamflow forecast system[END_REF][START_REF] Weerts | Estimation of predictive hydrological uncertainty using quantile regression: examples from the National Flood Forecasting System (England and Wales)[END_REF][START_REF] Ewen | Prediction intervals for rainfall-runoff models: raw error method and split-sample validation[END_REF][START_REF] Pianosi | Dynamic modeling of predictive uncertainty by regression on absolute errors[END_REF][START_REF] Smith | Adaptive correction of deterministic models to produce probabilistic forecasts[END_REF][START_REF] Van Steenbergen | A non-parametric data-based approach for probabilistic flood forecasting in support of uncertainty communication[END_REF][START_REF] Yan | Characterizing the uncertainty in river stage forecasts conditional on point forecast values[END_REF], but their use is not widespread for operational ensemble forecasting.

Although generally dealt with separately, statistical post-processing and data assimilation (also called real-time model updating in the engineering community) can be intrinsically related in the hydrological forecasting framework. Both represent techniques that may be used in a forecasting system to improve the quality of the forecasts (i.e., to provide more accurate and reliable forecasts) and to, ultimately, enhance the usefulness of the forecasts in decision-making. Since forecasting deals with an uncertain future, these techniques aim to bring additional information to the forecast procedure and take into account the various uncertainty sources (or at least the major uncertainty sources) affecting the forecasting chain. This is usually achieved by merging information from model and observations. While data assimilation and post-processing share a general goal, the techniques applied may differ in the practice of hydrological forecasting.

These differences usually draw the separation between what is defined as data assimilation and what is defined as post-processing in a modelling framework.

The definitions used in this study are the following: we use the term "postprocessing" when using the hydrological uncertainty processor (Section 2.4), whose primary purpose is to dress deterministic forecasts with uncertainty based on distributions of past model errors and, this way, build probabilistic forecasts. "Data assimilation" refers to techniques applied to perform the updating of the system before it issues a deterministic forecast. Here it concerns the state updating of the hydrological model and a model error correction applied to its output (Section 2.3).

The fact that data assimilation has the potential to improve real-time streamflow forecasting is widely accepted (see Liu et al., 2012, for a review).

In contrast to probabilistic and ensemble-based data assimilation methods (e.g., [START_REF] Weerts | Particle filtering and ensemble Kalman filtering for state updating with hydrological conceptual rainfall-runoff models[END_REF][START_REF] Salamon | Disentangling uncertainties in distributed hydrological modeling using multiplicative error models and sequential data assimilation[END_REF][START_REF] Moradkhani | Evolution of ensemble data assimilation for uncertainty quantification using the particle filter-Markov chain Monte Carlo method[END_REF][START_REF] Vrugt | Hydrologic data assimilation using particle Markov chain Monte Carlo simulation: Theory, concepts and applications[END_REF], deterministic updating schemes are designed to improve forecasts without producing probabilistic outputs. They may be easier to implement, mainly operationally, but at the price of leaving the uncertainty quantification issue unanswered. In these cases, the use of statistical post-processing methods together with data assimilation procedures provides a way to reduce and quantify the predictive uncertainty in the hydrological forecasts.

Integrating uncertainties in hydrological ensemble forecasting

"Ensemble dressing" is an intuitive and operationally-appealing method that allows integration of uncertainties from hydrological modelling and meteorological (ensemble) forcing. The main difference with other ensemblebased post-processors (e.g., [START_REF] Wang | Improvement of ensemble reliability with a new dressing kernel[END_REF][START_REF] Fortin | Probabilistic forecasting from ensemble prediction systems: Improving upon the best-member method by using a different weight and dressing kernel for each member[END_REF][START_REF] Brown | A Nonparametric Postprocessor for Bias Correction of Hydrometeorological and Hydrologic Ensemble Forecasts[END_REF][START_REF] Boucher | Hydro-economic assessment of hydrological forecasting systems[END_REF][START_REF] Brown | Evaluation of a nonparametric post-processor for bias correction and uncertainty estimation of hydrologic predictions[END_REF] is that, for ensemble dressing, hydrological modelling errors are assessed separately, and later combined with ensemble forecasts. Distributions of modelling errors are obtained from long time series of simulated and observed data (i.e., learning from the past), and then applied to ensemble forecasts to obtain the total predictive distribution.

In recent studies, the use of ensemble dressing has been implemented and tested to improve the skill of hydrological ensemble forecasting systems.

For instance, [START_REF] Reggiani | Uncertainty assessment via Bayesian revision of ensemble streamflow predictions in the operational river Rhine forecasting system[END_REF] present a Bayesian ensemble uncertainty processor for medium-range ensemble flow forecasts in the Rhine river basin. [START_REF] Hopson | A 1-10-Day Ensemble Forecasting Scheme for the Major River Basins of Bangladesh: Forecasting Severe Floods of 2003-07[END_REF] use an uncertainty processor based on the knearest neighbours (k-NN) resampling method to dress probabilistic mediumrange forecasts for two large basins in Bangladesh. [START_REF] Zalachori | Statistical processing of forecasts for hydrological ensemble prediction: a comparative study of different bias correction strategies[END_REF] compare different strategies based on pre-and post-processing methods to remove biases in a streamflow ensemble prediction system developed for reservoir inflow management in French catchments, while [START_REF] Pagano | Ensemble dressing for hydrological applications[END_REF] present a hydrological application of ensemble dressing for 128 catchments in Australia.

The studies mentioned above have in common the fact that they focus on post-processors for operational applications and on the overall evaluation of the quality of post-processed forecasts. Like in the studies that develop and test data assimilation techniques, most of the forecast assessment is on the benefits (in terms of quality) that post-processors or data assimilation may bring to forecast quality (accuracy, reliability, sharpness, etc.) at fixed forecast lead times. Little is known about the interactions between these two components of a forecasting system and the impacts of implementing both post-processing and data assimilation on the performance of the forecasts along the forecast lead times.

Aim and scope of the study

This study aims to shed light on the interactions between data assimilation and post-processing in hydrological ensemble forecasting. We address the following questions:

1. How does data assimilation impact hydrological ensemble forecasts? 2. How does post-processing impact hydrological ensemble forecasts? 3. How does data assimilation interact with post-processing to improve the quality and skill of hydrological ensemble forecasts over the forecast lead times?

We address these questions with the help of a large set of catchments, making it possible to draw more general and robust conclusions.

Data and methods

Data set

A set of 202 unregulated catchments spread over France was used (Figure 1). The catchments represent various hydrological conditions, given the variability in climate, topography, and geology in France. This set includes fast responding Mediterranean basins with intense precipitation as well as larger, groundwater-dominated basins. Some characteristics of the data set are given in Table 1. Catchments were selected to have limited snow influence, since no snowmelt module was used in the hydrological modelling (Section 2.3).

Potential evapotranspiration (PE), precipitation, and discharge data were available at hourly time steps over the 1997-2006 period. Temperature inputs originate from the SAFRAN reanalysis [START_REF] Vidal | A 50-year high-resolution atmospheric reanalysis over France with the Safran system[END_REF]. PE was estimated using a temperature-based formula [START_REF] Oudin | Which potential evapotranspiration input for a lumped rainfall-runoff model? Part 2 -Towards a simple and efficient potential evapotranspiration model for rainfall-runoff modelling[END_REF]. Precipitation data come from a reanalysis dataset recently produced by Météo-France based on weather radar and rain gauge network [START_REF] Tabary | A 10-year (1997-2006) reanalysis of Quantitative Precipitation Estimation over France: methodology and first results[END_REF]. River discharge data were extracted from the HYDRO national archive (www.hydro.eaufrance.fr).

! !! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! !! ! ! ! ! ! ! ! ! ± 100

PEARP, the Météo-France ensemble forecast

A short-range meteorological ensemble prediction system, the Météo-France PEARP EPS [START_REF] Nicolau | Short-range ensemble forecasting[END_REF], was used to produce hydrological ensemble forecasts. The PEARP EPS runs once a day at 18:00 UTC; it has 11 members, a 60 h forecast range, and a 0.25 • (ca. 25 km in France) grid resolution. A spatial disaggregation to an 8 km x 8 km grid, which includes bias correction, was applied to the PEARP forecasts. Bias correction was applied to precipitation forecasts using a multiplying factor obtained from a comparison between the mean of the PEARP ensemble and the Météo-France SAFRAN reanalysis over a complete year (March 2005to March 2006). Details can be found in [START_REF] Thirel | On the Impact of Short-Range Meteorological Forecasts for Ensemble Streamflow Predictions[END_REF] PEARP forecasts were already used at the daily time step in recent hydrological studies [START_REF] Thirel | On the Impact of Short-Range Meteorological Forecasts for Ensemble Streamflow Predictions[END_REF][START_REF] Randrianasolo | Comparing the scores of hydrological ensemble forecasts issued by two different hydrological models[END_REF]. Overall, they showed good quality over France at this time step. The quality for shortterm forecasting at hourly time steps (with either raw and post-processed forecasts) is first assessed here.

The GRP rainfall-runoff forecasting model

The GRP model is a continuous, lumped storage-type model designed for flood forecasting. Its structure was derived from the GR4J model [START_REF] Perrin | Improvement of a parsimonious model for streamflow simulation[END_REF] and is composed of a production function and a routing function.

The production function consists of a non-linear soil moisture accounting (SMA) reservoir and a volume adjustment coefficient. The routing function includes a unit hydrograph and a non-linear routing store. The GRP model uses catchment areal rainfall and PE as inputs; it is parsimonious with three parameters to be calibrated against observed data: one in the production function (the volume adjustment coefficient) and two for the routing function (the base time of the unit hydrograph and the total capacity of the routing store). In this study, the three free parameters were calibrated for each catchment by minimizing the root mean square errors (RMSE) during the first five years of available data (1997)(1998)(1999)(2000)(2001).

Importantly, the hourly version of the GRP model uses together two data assimilation procedures for flood forecasting. The first exploits the last available observed discharge to directly update the routing store state, and the second exploits the last relative error to correct the model output with a multiplicative coefficient. More details about the forecasting model GRP and the two assimilation procedures can be found in [START_REF] Berthet | How crucial is it to account for the antecedent moisture conditions in flood forecasting? Comparison of event-based and continuous approaches on 178 catchments[END_REF].

Hydrological uncertainty processor

We used a hydrological uncertainty processor (HUP) to evaluate the conditional errors of the hydrological model. Only hydrological uncertainty is considered by the HUP here since the model is run with observed weather data. The meteorological uncertainty is subsequently considered through the joint use of the HUP with the PEARP forecasts, as described in Section 2.5.

The HUP used here is a data-based and non-parametric method that was applied by [START_REF] Andréassian | What is really undermining hydrologic science today?[END_REF] to assess model simulation uncertainties and compute empirical uncertainty bounds to flow simulations. Here it is applied to produce probabilistic flow forecasts. The basic idea is to estimate empirical quantiles of relative errors stratified by different flow groups. The HUP is trained during the period used for calibrating the parameters of the hydrological model (1997)(1998)(1999)(2000)(2001). Note that it is possible that this approach yield optimistic uncertainty estimates, since errors are usually larger on an independent period than during the calibration period. Since forecast error characteristics vary with forecast range when data assimilation is used, the HUP is trained at several lead times separately.

For each catchment, the HUP is trained as described below:

Step 1. The hydrological model is run with observed weather data as input and the time series of relative errors is evaluated:

Q f ct /Q obs , where (Q f ct ,
Q obs ) are the pairs of discharge forecasts and observations.

Step 2. The time series is stratified into 20 groups according to the magnitude of the Q f ct . The limits of each group are fixed so that each group contains the same number of values.

Step 3. Within each group, an empirical distribution of relative errors is defined and 99 quantiles are estimated (corresponding to the percentiles 1%, 2%, ... 98%, 99%).

Application of the HUP for another forecast period is described by the last step:

Step 4. Once defined during the training period, the empirical quantiles of relative errors can be applied to any forecast discharge at a certain lead time. The limits of each group are the same as those obtained during the training period. Note that when data assimilation is not used, the empirical quantiles of relative errors are the same whatever the forecast lead time is. Given a discharge forecast Q f ct , we first determine the flow group Q f ct belongs to; then Q f ct is multiplied by the 99 quantiles of relative errors; the 99 values obtained describe the predictive distribution at the considered time step and for a given forecast horizon. In cases of extrapolation (i.e., when the forecast discharge is out of the range of the flow groups defined during the training phase of the HUP), values of relative errors from the nearest flow groups (i.e., the lowest or the highest flow groups) are used.

Preliminary studies carried out to compare this approach to other similar post-processing approaches suggest that it can yield similar results in terms of forecast performance, while being simpler in its application.

Ensemble dressing method: an integrator of the meteorological and hydrological uncertainties

The ensemble dressing method is used as an integrator of the meteoro- Application and evaluation of the ensemble dressing method for the ensemble forecasts is done over an independent period, the 17-month period from August 2005 to December 2006.

Experiments

The hydrological ensemble forecast system combines meteorological and streamflow data from observation networks, the Météo-France PEARP ensemble forecast, the GRP rainfall-runoff model with its two data assimilation functions, the hydrological uncertainty processor (HUP) and the ensemble dressing method. Hereafter we will use the term "post-processing" to describe the joint use of the HUP and the ensemble dressing method, while the term "data assimilation" will refer to the two updating techniques used in the GRP model.

In order to assess the benefits of data assimilation and post-processing, considered together or separately, different configurations of the forecasting chain were analysed. Our experiments comprise a chain without data assimilation and post-processing (NoDA-NoPP), without data assimilation but with post-processing (NoDA-PP), with data assimilation but without postprocessing (DA-NoPP), and with both data assimilation and post-processing (DA-PP). The characteristics of the experiments and the acronyms used are given in Table 2.

In particular, the NoDA-NoPP experiment corresponds to the situation where the hydrological model is run in simulation mode, i.e., without using recent streamflow observations for data assimilation, and is then driven by the PEARP ensemble forecast when the forecast is issued. When data assimilation is used, the state of the routing reservoir of the hydrological model is first updated based on the last observed discharge, and the second procedure is then applied separately at each streamflow ensemble member. This structured analysis allows us to identify the influence of data assimilation and post-processing separately to assess the benefits of both components when used together in the forecasting chain.

Without data assimilation With data assimilation

Without post-processing NoDA-NoPP DA-NoPP

With post-processing NoDA-PP DA-PP

Forecast evaluation methods

The evaluation of the performance of probabilistic forecasts should reflect the different facets of probabilistic forecasts. In this study, the forecasts obtained from the four experiments set up (Table 2) were evaluated with both deterministic and probabilistic scores. We aimed to assess the influence of data assimilation and post-processing on the following characteristics of ensemble forecasts: accuracy of the ensemble mean, overall sharpness and reliability of the whole ensemble, and overall forecast quality of the ensemble.

More specifically, we evaluated the accuracy of the ensemble mean values with the relative bias (BIAS) and the normalized root-mean-square error (NRMSE). To assess the overall reliability of the forecasts, we used the Probability Integral Transform (PIT) diagram (see e.g., [START_REF] Laio | Verification tools for probabilistic forecasts of continuous hydrological variables[END_REF][START_REF] Thyer | Critical evaluation of parameter consistency and predictive uncertainty in hydrological modeling: A case study using Bayesian total error analysis[END_REF] and an index that quantifies deviation from the ideal case, the alpha score [START_REF] Renard | Understanding predictive uncertainty in hydrologic modeling: The challenge of identifying input and structural errors[END_REF]. The overall sharpness of the forecasts was measured with an index based on the interquartile range that we called normalized mean interquartile range (NMIQR). Finally, we assessed the overall forecast quality of the whole ensemble with the mean Continuous Rank Probability Skill Score (mean CRPSS). The mean CRPSS is computed with the unconditional streamflow climatology as the reference. These scores are presented in more details in A. Nonetheless, for all lead times, forecast accuracy is best when DA and PP are used together, which indicates the benefits of the combined use of data assimilation and post-processing. 2) and lead times 6 h, 12 h, 24 h and 48 h. Boxplots (5th, 25th, 50th, 75th and 95th percentiles) synthesize the variety of scores over the 202 catchments of the data set.

Results and discussion

Forecast accuracy

Reliability

Figure 3 presents the PIT diagrams obtained for each of the 202 catchments, when considering 24 h ahead ensemble forecasts. Since similar figures were obtained for the other lead times (not shown).

From Figures 3a and3b, it can be seen that most of the curves are almost horizontal straight lines, while they would follow the bisector (black lines in the graphs) in the ideal case of reliable ensemble predictions. Figures 3a and3b clearly reveal that the raw ensembles are lacking reliability for all of the catchments. The impact of post-processing on reliability is apparent when looking at the results on Figures 3c and3d: the curves of the ensemble streamflow forecasts with post-processing follow the ideal situation much more closely than the curves shown in Figures 3a and3b The PIT diagrams convey a visual evaluation of the overall reliability of probabilistic forecasts. To quantify it, we used the alpha score, a reliability index that measures the deviation of the PIT curves from the ideal situation. These results suggest that for the 202 catchments studied the spread obtained by propagating solely the precipitation ensembles into the hydrological model is too small to properly reflect the range of forecast errors. The deterministic data assimilation strategy used here is effective in improving the reliability of the ensemble forecasts, but it is not self-sufficient to correct the under dispersion of the streamflow ensemble forecasts as revealed by the PIT diagrams in Figure 3 and the alpha scores in Figure 4. This is a strong indication that the hydrological uncertainty issue should be specifically addressed in order to improve the overall reliability of hydrological ensemble forecasts.

Sharpness

Sharpness is a desirable characteristic of any probabilistic forecast. The sharper the forecast, the less uncertain it is, and thus the more information is conveyed. The four experiments we used made it possible to investigate how meteorological and hydrological uncertainties interact and affect sharpness.

Figure 5 shows the distributions of a sharpness index, the normalized mean interquartile range (NMIQR), over 202 catchments. 2) and for lead times 6 h, 12 h, 24 h and 48 h.

Boxplots (5th, 25th, 50th, 75th and 95th percentiles) synthesize the variety of scores over the 202 catchments of the data set. Perfect score is 1.0.

It can be seen that the ensemble spreads of three experiments, NoDA-NoPP, DA-NoPP and DA-PP, increase significantly with increasing lead time, while it is more stable over lead times for the experiment NoDA-PP. For NoDa-NoPP and DA-NoPP, the median value of NMIQR over the 202 catchments raises in a very close behaviour for both experiments, from around 0.05 for 6 h ahead forecasts to 0.13 for 48 h ahead forecasts. For the experiment DA-PP, the increase in the median values is much more important: from 2) and for lead times 6 h,

Mean CRPSS

The analysis of the impacts of data assimilation and post-processing on two important characteristics of probabilistic forecasts, reliability and sharpness, showed that post-processing was necessary to improve reliability, but at the cost of lower sharpness, i.e., greater ensemble spread and uncertainty, even if sharpness could be improved with the application of a data assimilation procedure. We now turn our attention to the mean CRPSS, a probabilistic score that provides an assessment of the overall quality of ensemble forecasts.

Figure 6 shows the distributions of the mean CRPSS over 202 catchments.

We note that performance decreases with increasing lead time for the two ex- Furthermore, the comparison with the reference climatology shows that data assimilation alone is sufficient to generate skillfull forecasts for more than 95% of the catchments for lead times up to 24 h, but post-processing (DA-PP) is necessary to achieve forecasts that have better overall performance than climatology at 48 h.

These results show the general added value of data assimilation and postprocessing to the overall quality of ensemble forecasts. When evaluating the overall quality of ensemble forecasts with the CRPSS, the benefits in terms of reliability overcome the loss of sharpness that results from accounting for hydrological uncertainty. The streamflow ensemble forecasts that explicitly account for both sources of uncertainty, meteorological and hydrological uncertainties, through post-processing, while reducing as much as possible hydrological uncertainty, here through data assimilation, are the most skillfull forecasts. 2) and for lead times 6 h, 12 h, 24 h and 48 h. Boxplots (5th, 25th, 50th, 75th and 95th percentiles) synthesize the variety of scores over the 202 catchments of the data set. Perfect score is 1.0.

Summary and conclusions

We investigated the relative contributions of data assimilation and postprocessing to the skill of hydrological ensemble forecasts. The study assessed the benefits of data assimilation and post-processing with the help of four configurations of a short-range hydrological ensemble forecasting system:

without data assimilation and post-processing (NoDA-NoPP), without data assimilation but with post-processing (NoDA-PP), with data assimilation but without post-processing (DA-NoPP), and with both data assimilation and post-processing (DA-PP).

We applied deterministic and probabilistic scores to streamflow forecasts of a large catchment set which brought into light the main general conclusions listed below:

• We verify the well-known fact that short-range hydrological forecasts benefit from data assimilation. Data assimilation has a strong impact on improving the quality of the ensemble mean, and a much lesser effect on the variability of the ensemble members (i.e., their spread).

• The benefits of a simple yet efficient hydrological uncertainty processor to improve the reliability and the overall quality of the short-range hydrological ensemble forecasts were demonstrated. Post-processing has a strong impact on forecast reliability.

• The benefits of the combined use of data assimilation and post-processing were demonstrated: both contribute to achieve reliable and sharp forecasts, with impacts acting differently according to the target lead time.

The stronger impact on forecast reliability comes from the use of postprocessing. Adding data assimilation to the system helps in improving sharpness and reliability at all lead times, with higher gains in performance at shorter lead times.

We acknowledge some limitations. It was only possible to evaluate the forecasting chain over a 17-month period of ensemble forecasts, since this was the common period between observations and forecasts we had available.

Furthermore, PEARP ensembles are ran only once a day, which limits the number of hourly evaluation pairs. For these reasons, it was not possible to evaluate flows over specific flooding thresholds . However, with increasing data archives, we expect that such an issue will be treated in future work.

Our study considered only one data assimilation technique (state updating with error output correction) and one post-processing method (ensemble dressing with hydrological errors) together with one rainfall-runoff model forecasting (GRP model). There are several other techniques and models in the literature that could also be tested using the methodology presented here. For instance, a comparison between different configurations of the method used, or different hydrological uncertainty processors, including methods that take into account the autocorrelation of errors (e.g., Schoups

and Vrugt, 2010) could be investigated. Besides, while a bias correction was applied to the PEARP forecasts, a more sophisticated pre-processor (see e.g., [START_REF] Verkade | Post-processing ecmwf precipitation and temperature ensemble reforecasts for operational hydrologic forecasting at various spatial scales[END_REF] could be used to further investigate how meteorological and hydrological biases interact and contribute to the quality of the final hydrological ensemble.

Also, the effectiveness of a data assimilation technique or a post-processing method (and hence the choice of the procedures to operate in a forecasting system) is affected by different sources of uncertainties present in a flow forecasting system, including the forcing data, initial conditions, parameter uncertainty and model structural uncertainty. In our study, we followed the works of [START_REF] Krzysztofowicz | Bayesian theory of probabilistic forecasting via deterministic hydrologic model[END_REF] and focused on a decomposition of the total uncertainty into meteorological and hydrological uncertainty. Observational or parameter uncertainties were thus not explicitly considered. Additional sources of uncertainty may however affect the performance of data assimilation techniques and post-processors, as well as the way they interact in the forecasting system. Further investigations would be necessary to better assess the extent to which this may affect forecast quality.

Although our findings may be related to the configuration used, they are based on common techniques and on the study of a large set of catchments, which helps in giving robustness and generality to the results obtained. The study also shows that, for a given system configuration, it is interesting to analyse how data assimilation and/or post-processing techniques set up to improve forecast quality affect the attributes of the forecasts and interact to provide overall good forecasts. The aim of a forecaster may then be to achieve a good combination of hydrological model, data assimilation and post-processing procedures that results in an overall good quality of his/her operational system (eventually over specific space and time scales of interest), rather than to search for the best data assimilation technique or post-processor available, without taking into account how they will interact between them and with the probabilistic forecasting system as a whole.

Despite those limitations, our results strongly suggest that data assimilation and post-processing techniques based on hydrological uncertainty processors should be more widely tested to foster their implementation in preoperational and operational hydrological ensemble forecasting systems and their use in real-time probabilistic forecasting. The use of both strategies is highly recommended since they have complementary effects: data assimila-tion has a very positive effect on forecast accuracy, and thus helps reduce hydrological uncertainty, but its impact diminishes with lead time, while post-processing, by accounting for hydrological uncertainty, has a very positive and longer lasting effect on forecast reliability.

The mean Continuous Rank Probability Skill Score (CRPSS) is a skill score based on the CRPS. Skill scores (SS) are used to assess the relative quality of two forecasting systems. They are generally defined as:

SS = 1 - Score A Score B ( 7
)
where Score A and Score B are the scores of the forecasting system A and B respectively. The forecasting system B is usually termed the reference forecast.

Climatology is commonly used as a reference. To compute the mean CRPSS with the unconditional climatology as the reference, an unconditional streamflow ensemble forecast is first obtained from the empirical distribution of all observed discharges over the evaluation period, and then used for all forecast occasions. 2). Each line represents one of the 202 catchments of the data set.

1

Highlights

• Data assimilation and post-processing impact hydrologic ensemble forecasts' skill.

• Data assimilation has a strong impact on forecast accuracy.

• Post-processing has a strong impact on forecast reliability.

• The combined benefits of data assimilation and post-processing were demonstrated.

• We recommend the use of both data assimilation and post-processing in forecasting.

Figure 1 :

 1 Figure 1: Locations of the 202 French catchments used in this study (dots correspond to the gauging stations, and blue color is catchment areas).

  logical and hydrological uncertainties. It consists in two steps. Firstly, each time an ensemble PEARP forecast is available, the hydrological model is run with the ensemble forecast and the HUP is applied, according to Step 4 of Section 2.4, to each of the 11 members of the ensemble for each lead time considered. Secondly, the 11 × 99 values obtained at each lead time are pooled together and an empirical cumulative distribution is estimated. From this distribution, 99 quantiles are retained as the members of the dressed ensemble.

Figure 2

 2 Figure2shows the distributions of the two deterministic scores used to assess forecast accuracy: the relative bias (BIAS) and the normalized rootmean-square error (NRMSE). Each score is computed for lead times 6 h, 12 h, 24 h and 48 h and for all 202 catchments. The distribution of the 202 values is summarized with boxplots.We note that forecast accuracy decreases with increasing lead time for the four experiments. For NoDA experiments (NoDA-NoPP and NoDA-PP), the loss of performance is quite limited: it is only related to the decreasing performance of the PEARP ensemble precipitation forecasts. For DA experiments (DA-NoPP and Da-PP), the decrease is stronger and the performances converge toward those of NoDA experiments: the effects of the two DA procedures used in the GRP forecasting model vanish with larger horizons; the decrease in performance of the hydrological model is then added to the losses in performance of the PEARP ensemble precipitation forecasts. Figure2also reveals that post-processing does not significantly impact forecast accuracy, whether or not DA is used. DA has a much stronger impact on the ensemble mean values than post-processing, especially for shorter lead times and, to a lower extent, for larger lead times. The two DA procedures used in the GRP forecasting model have been designed to improve the performance of deter-

Figure 2 :

 2 Figure 2: Distributions of two deterministic scores, the relative bias (BIAS) and the normalized root-mean-square error (NRMSE), for ensemble streamflow forecasts from the four experiments (see Table2) and lead times 6 h, 12 h, 24 h and 48 h.Boxplots (5th, 25th,

  (ensemble streamflow forecasts without post-processing). It means that the overall reliability of the ensembles is clearly improved with post-processing and this for both cases, with and without DA. A comparison of solely Figures3c and 3dconfirms also the positive impact of data assimilation on the reliability of the ensembles: the PIT curves of the dressed ensembles are substantially closer to the diagonal (perfect reliability) when DA is applied.

Figure 4

 4 Figure 4 presents the distributions of the alpha scores obtained for each experiment over the 202 catchments. Results in Figure 4 confirm the visual evaluation obtained with the PIT diagrams: the two experiments that do not account for hydrological uncertainty (NoDA-NoPP and DA-NoPP) lack

Figure 4 :

 4 Figure 4: Distributions of the alpha score reliability index for streamflow ensemble forecasts from the four experiments (see Table2) and for lead times 6 h, 12 h, 24 h and 48 h.

Figure 5 :

 5 Figure 5: Distributions of the normalized mean interquartile range (NMIQR) for streamflow ensemble forecasts from the four experiments (see Table2) and for lead times 6 h,

  periments with data assimilation: median values of the CRPSS are equal to 0.84 (DA-NoPP) and 0.87 (DA-PP) for 6 h range forecasts, and equal to 0.45 (DA-NoPP) and 0.57 (DA-PP) for 48 h range forecasts. Mean CRPSS values of the two experiments without data assimilation decrease only slightly but are much lower than values obtained with data assimilation (median values around 0.10 for NoDA-NoPP and around 0.45 for NoDA-PP). This is especially true for shorter lead times and, to a lower extent, for larger lead times.

Figure 6 :

 6 Figure 6: Distributions of the mean CRPSS for streamflow ensemble forecasts from the four experiments (see Table2) and for lead times 6 h, 12 h, 24 h and 48 h. Boxplots

Figure 3 .

 3 Figure 3. PIT diagrams of the 24 h ahead streamflow ensemble forecasts from the four a-d) experiments (see Table2). Each line represents one of the 202 catchments of the data set.

Table 1 :

 1 . PEARP forecasts were available

	Percentiles
	0.05 0.25 0.50 0.75 0.95

Characteristics of the 202 catchments. P-precipitation, PE-potential evapotranspiration, Q-discharge. over the 2005-2009 period, but only the period matching the observed data could be used here, i.e. from August 2005 to December 2006.

Table 2 :

 2 Acronyms used for the different experiments used in this study.

h, 24 h and 48 h. Boxplots (5th, 25th, 50th, 75th and 95th percentiles) synthesize the variety of scores over the 202 catchments of the data set. Perfect score is 0.

Acknowledgments

The authors thank Météo-France for providing the meteorological data and Banque HYDRO for the hydrological data. The financial support of SCHAPI to the first author is also gratefully acknowledged.

The authors thank Dr. Thomas Pagano, two anonymous reviewers, and the Guest Editor Hamid Moradkhani for their critical and constructive evaluation of the manuscript, which helped improving its quality.

A. Evaluation scores

The evaluation scores used in this article are defined and briefly described below. For more details, the reader may refer to [START_REF] Wilks | Statistical methods in the atmospheric sciences[END_REF].

A.1. Relative Bias

The relative bias (BIAS) is defined as the ratio between the mean of deterministic forecasts and the mean of observations, BIAS =

where

is the kth of N pairs of deterministic forecasts and observations.

Values higher (lower) than 1 indicate an overall overestimation (underestimation) of the observed values.

A.2. Normalized root-mean-square error

The root-mean-square error (RMSE) is a widely used measure of accuracy for point forecasts,

where

is the kth of N pairs of forecasts and observations.

The lower the RMSE, the better. For a perfect deterministic forecast, RMSE=0.

The normalized root-mean-square error (NRMSE) is obtained by dividing the RMSE by the mean runoff. The use of a non-dimensional score facilitates the comparison of the results obtained over different catchments.

A.3. PIT diagram and alpha score

The Probability Integral Transform (PIT) diagram is a graphical tool used to assess the reliability of probabilistic forecasts [START_REF] Gneiting | Probabilistic forecasts, calibration and sharpness[END_REF][START_REF] Laio | Verification tools for probabilistic forecasts of continuous hydrological variables[END_REF].

A.4. Normalized mean interquartile range

To assess the sharpness of probabilistic forecasts, we defined the mean interquartile range (MIQR) as the mean of the interquartile range of forecasts over the evaluation data. The interquartile range, defined as the range between the upper quartile (75th percentile) and the lower quartile (25th percentile) of a distribution, is a robust measure of the spread of a distribution.

MIQR is computed as

where

)

is the kth of N pairs of quartiles of the forecasts.

Similarly to the NRMSE, we divided the MIQR by the mean runoff to obtain a non-dimensional score.

A.5. Mean CRPS and mean CRPSS

For a forecast-observation evaluation pair, the Continuous Rank Probability Score (CRPS) (e.g., [START_REF] Matheson | Scoring rules for continu-ous probability distributions[END_REF][START_REF] Gneiting | Probabilistic forecasts, calibration and sharpness[END_REF] measures the quadratic distance between two cumulative distribution functions, the cumulative predictive distribution F (x) and a Heaviside function based on the observed value 1{Q obs x} ) :

The mean CRPS, CRPS, is the average value of the CRPS over the N pairs of evaluation data: