Romain Gojon
email: romain.gojon@ec-lyon.fr.

Christophe Bogey
email: christophe.bogey@ec-lyon.fr.

Olivier Marsden
email: olivier.marsden@ec-lyon.fr.

Large-eddy simulation of supersonic planar jets impinging on a flat plate at an angle of 60 to 90 degrees

Six planar supersonic jets are computed by compressible large-eddy simulations (LES) using low-dissipation schemes. At the exit of a nozzle of height h, they are ideally expanded and have an exit velocity ue, yielding a Mach number of Me = 1.28, and a Reynolds number of Reh = ueh/ν = 5 × 10 4 . Four jets impinging on a flat plate at distances L from the nozzle lips ranging from 3.94h up to 9.1h, with an angle θ between the jet direction and the plate of 90 degrees, are first considered. Two other jets with L = 5.5h and θ = 60 and θ = 75 degrees are also examined. In this way, the effects of both the nozzle-to-plate distance and the angle of impact on the flow and acoustic fields of the jet are studied. Mean velocity flows and snapshots of density, pressure and vorticity are shown. The pressure fields are also described by computing sound pressure levels and using Fourier decomposition. Several tones are obtained in certain cases and their corresponding Strouhal numbers and symmetric or antisymmetric natures are in agreement with experimental data and theoretical models. They are due to an aeroacoustic feedback mechanism occurring between the nozzle lips and the flat plate. This mechanism generates an hydrodynamicacoustic standing wave revealed by using Fourier decomposition.

I. Introduction

Jets impinging on a flat plate have been studied experimentally by many authors over the past sixty years. In some cases, very intense tones have been observed in the acoustic field. In his pioneering work, in 1953, Powell [START_REF] Powell | On edge tones and associated phenomena[END_REF] suggested that they are due to a feedback mechanism between hydrodynamic instabilities propagating downstream from the nozzle lip to the plate and acoustic waves propagating upstream from the plate to the nozzle lips.

Round subsonic jets impinging on a flat plate normally for an exit Mach number ranging from 0.3 up to 0.9 have been studied notably by Ho and Nosseir 2 and Nosseir and Ho 3 in the eighties. The near pressure field was described using microphones and a simplified model predicting the frequencies of the feedback mechanism based on cross-correlations between the microphones was proposed. Round supersonic jets impinging on a flat plate normally have also been investigated experimentally by Henderson and Powell, [START_REF] Henderson | Experiments concerning tones produced by an axisymmetric choked jet impinging on flat plates[END_REF] Krothapalli et al. [START_REF] Krothapalli | Flow field and noise characteristics of a supersonic impinging jet[END_REF] and Henderson et al. [START_REF] Henderson | An experimental study of the oscillatory flow structure of tone-producing supersonic impinging jets[END_REF] and more recently by Risbord and Soria 7 and Buchmann et al. [START_REF] Buchmann | High spatial resolution imaging of a supersonic underexpanded jet impinging on a flat plate[END_REF] using high-speed optical measurements. In some cases, a feedback mechanism similar to those in subsonic jets was found. This mechanism appears very often when the jet is ideally expanded but only for some nozzle-to-plate distances when the jet is imperfectly expanded. Henderson et al. [START_REF] Henderson | An experimental study of the oscillatory flow structure of tone-producing supersonic impinging jets[END_REF] suggested that in the latter case, the feedback establishes only when a Mach disk forms upstream from the plate.

Planar jets impinging on a flat plate normally also produce intense tone frequencies. However, compared to round jets, tones are visible at lower flow velocities, as noted by Arthurs and Ziada. [START_REF] Arthurs | Self-excited oscillations of a high-speed impinging planar jet[END_REF] Planar supersonic jets impinging on a flat plate normally has been studied by Norum [START_REF] Norum | Supersonic rectangular jet impingement noise experiments[END_REF] and Tam and Norum [START_REF] Tam | Impingement tones of large aspect ratio supersonic rectangular jets[END_REF] in the nineties. These authors obtained, for certain jets, acoustic spectra with two emerging tone frequencies. The lower tone frequency is due to a varicose (symmetric) mode of the jet, and the upper tone frequency is due to a sinuous (antisymmetric) mode. Similar tone frequencies have been found in a 2-D simulation by Hourigan et al. [START_REF] Hourigan | The feedback loop in impinging two-dimensional high-subsonic and supersonic jets[END_REF] In order to understand the varicose or sinuous behaviour of the jet, Tam and Norum [START_REF] Tam | Impingement tones of large aspect ratio supersonic rectangular jets[END_REF] proposed that the upstream propagating waves of the feedback mechanism are neutral acoustic wave modes of the jet. Thus, both the wave number k and the angular frequency ω of these modes are real. They found an allowable frequency range for each upstream propagating neutral acoustic wave mode of the jet flow using a vortex sheet jet model and noted that the two emerging tones observed in the experiment of Norum 10 fall into the first symmetric and the first antisymmetric allowable bands given by the model. Later, several tone frequencies for supersonic rectangular impinging jets have also been detected experimentally by Thurow et al. [START_REF] Thurow | Structure of a supersonic impinging rectangular jet via real-time optical diagnostics[END_REF] using a real-time flow visualization technique. A coupling between the organization of the hydrodynamic structures in the shear layers of the jets and the degree of resonance of the feedback mechanism was pointed out.

For non-normal impingement angle, the intensity of the feedback mechanism decreases because of the loss of symmetry. Indeed, Norum [START_REF] Norum | Supersonic rectangular jet impingement noise experiments[END_REF] conducted an experiment for a supersonic ideally expanded planar jet impinging on a inclined flat plate and noted that the amplitude of the strongest impingement tone decreases when the angle of impact deviates from 90 degrees. Besides, to understand the lift-off phase of a space launcher, Nonomura et al. 14,[START_REF] Nonomura | Aeroacoustic waves generated from a supersonic jet impinging on an inclined flat plate[END_REF] performed the 2-D simulation of a jet impinging on a inclined plate and the LES of a 3-D supersonic ideally expanded round jet impinging on a inclined plate. They obtained three types of waves in the acoustic field, namely the Mach waves generated in the jet shear layers, the Mach waves from the shear layer of the supersonic flow on the inclined flat plate, and the waves coming from the impingement region. These three types of waves were analysed by Tsutsumi et al. [START_REF] Tsutsumi | Analysis of Acoustic Wave from Supersonic Jets Impinging to an Inclined Flat Plate[END_REF] using a snapshot proper orthogonal decomposition with Fourier transformation. Moreover, the transient regime of supersonic jets was investigated by Dargaud et al. [START_REF] Dargaud | Numerical Study of Solid-Rocket Motor Ignition Overpressure Wave Including Infrared Radiation[END_REF] by conducting the simulation of the ignition overpressure wave generated by horizontal rocket motors interacting with the ground.

In the present work, LES of six planar jets are carried out in order to investigate the feedback mechanism between the nozzle lips and the flat plate. In particular, their spectral and hydrodynamic properties are studied and compared to empirical and theoretical models. The paper is organized as follows. The main characteristics of the jets and the simulations are presented in section II. Mean fields, snapshots and acoustical results are shown in section III and compared with experimental data. The feedback mechanism is then studied in section IV by applying Fourier decomposition to the pressure field. Finally, concluding remarks are given in section V.

II. Parameters

II.A. Jets parameters

Large-eddy simulations are performed for four jets impinging normally on a flat plate for nozzle-to-plate distances L/h of 3.94, 5.5, 8.27, and 9.1. For the nozzle-to-plate distance L/h = 5.5, two other simulations are carried out for the angles of impact θ = 60 and 75 degrees between the jet axis and the flat plate. The different cases are presented in table 1. They are referred to as JetL3, JetL5, JetL8, JetL9, JetL5-60d and JetL5-75d, respectively. The jets originate from a planar nozzle of height h, and width l = 3.25h in the spanwise direction. The lip thickness is e = 0.5h, and periodic conditions are imposed in the spanwise direction, allowing planar jets to be investigated. The ejection conditions of the jets and the nozzle-to-plate distances are similar to the parameters in the experimental study of Thurow et al. [START_REF] Thurow | Structure of a supersonic impinging rectangular jet via real-time optical diagnostics[END_REF] Thus, the jets are ideally expanded, and have a Mach number of M e = 1.28. Their Reynolds number is Re h = u e h/ν = 5 × 10 4 . A Blasius mean velocity profile is imposed at the nozzle exit with a boundary-layer thickness of δ BL = 0.075h.

II.B. Numerical parameters

In the four cases where the angle between the jet axis and the flat plate is 90 degrees, the unsteady compressible Navier-Stokes equations are solved on a Cartesian mesh for which x, y and z denote the longitudinal, lateral and spanwise directions, respectively. An explicit six-stage Runge-Kutta algorithm is used for time integration, and low-dispersion eleven-point explicit finite differences are used for spatial derivation. [START_REF] Bogey | A family of low dispersive and low dissipative explicit schemes for flow and noise computations[END_REF][START_REF] Berland | High-order, low dispersive and low dissipative explicit schemes for multiple-scale and boundary problems[END_REF] At the end of each time step, a relaxation filtering is applied to the flow conservative variables in order to remove grid-to-grid oscillations and to dissipate subgrid-scale turbulent energy. The radiation conditions of Tam and Dong 20 are implemented at the inflow and lateral boundaries of the computational domain. A sponge zone combining grid stretching and Laplacian filtering is also employed to damp the turbulent fluctuations before they reach the lateral boundaries. Adiabatic conditions are imposed to the nozzle wall and the flat plate. Examples of 3-D jets simulated using the present solver can be found in Bogey and Bailly [START_REF] Bogey | Large eddy simulations of transitional round jets: influence of the Reynolds number on flow development and energy dissipation[END_REF][START_REF] Bogey | Turbulence and energy budget in a self-preserving round jet: direct evaluation using large eddy simulation[END_REF] for transitional round jets at a Mach number M e = 0.9 and for self-similar turbulent round jets at a Mach number M e = 0.9. A shock-capturing filtering is applied in order to avoid Gibbs oscillations near shocks. It consists in applying a conservative second-order filter optimised in Fourier space at a magnitude determined each time step using a shock sensor. [START_REF] Bogey | A shock-capturing methodology based on adaptative spatial filtering for high-order non-linear computations[END_REF] This method was successfully used by de Cacqueray et al. [START_REF] De Cacqueray | Investigation of a High-Mach-Number Overexpanded Jet Using Large-Eddy Simulation[END_REF] for the LES of an overexpanded jet at a fully developed Mach number M j = 3.3. Finally, vortical disturbances not correlated in the spanwise direction [START_REF] Bogey | Large-eddy simulation of the flow and acoustic fields of a Reynolds number 10 5 subsonic jet with tripped exit boundary layers[END_REF] are added in the boundary layer in the nozzle near x = -0.5h and y = 0.45h in order to generate velocity fluctuations at the nozzle exit. The strength α of the forcing is settled to 0.02 to obtain turbulent intensities between 5% and 10% at the nozzle exit, as it will be reported in the next section. This value is similar to those used used in Bogey et al. In the two cases where the angle between the jet axis and the flat plate is 60 and 75 degrees, two Cartesian meshes are used as presented in figure 1. The mesh containing the nozzle is referred to as the primary grid and the mesh close to the flat plate as the secondary grid. The two meshes are identical in the spanwise direction z, allowing the flow variables to be transferred from one grid to another using a 2-D interpolation in the (x, y) plane. For the interpolation, a high-order Lagrangian 2-D interpolation based on a 10-point stencil is used in each direction. The numerical methods used are identical to those of the four LES with θ = 90 degrees.

The simulations are carried out using an OpenMP-based in-house solver, and a total of 200, 000 iterations are computed in each case for the steady state. The simulation time is thus equal to 500h/u e . The Cartesian meshes contain between 184 and 263 million points for the primary grids and between 38 and 46 million points for the secondary grids, as noted in table 2 and in table 3. For JetL3, JetL5, JetL8 and JetL9, the variations of the axial and lateral mesh spacings are represented in figure 2. The minimal axial mesh spacing, near the nozzle lip and the flat plate, is equal to ∆x = 0.00375h and the maximal axial mesh spacing, between the nozzle and the plate, is equal to ∆x = 0.015h. The lateral mesh spacings are equal to ∆y = 0.00375h at y = ±h/2 and to ∆y = 0.03h for 2.5h ≤ y ≤ 8h. The spanwise spacings are equal to ∆z = 0.015h. This discretization allows acoustic waves with Strouhal numbers up to St h = f h/u e = 5.6 to be well propagated, where f is the frequency. The grids are stretched at rates lower than 1% outside the sponge zones in order to preserve numerical accuracy. For JetL5-60d and JetL5-75d, the primary meshes are designed with a minimal axial mesh spacing near the nozzle lip equal to ∆x = 0.00375h and a maximal axial mesh spacing downstream from the nozzle equal to ∆x = 0.015h. The lateral mesh spacings are similar to those for JetL3. The secondary meshes have a minimal axial mesh spacing equal to ∆x = 0.00375h close to the plate and a maximal axial mesh spacing equal to ∆x = 0.015h upstream from the plate. The lateral mesh spacings are similar to those for jetL3 and are designed to minimize the difference in mesh size between the two grid as illustrated in figure 1(b).

III. Results

III.A. Mean flow fields

Mean velocity fields obtained in the (x, y) plane are shown in figure 3 for JetL3, JetL5, JetL8 and JetL9 and in figure 4 for JetL5, JetL5-75d and JetL5-60d. As expected, very little variations (around 4% of the jet exit velocity) of the mean axial velocity are noticed near the axis, indicating that the jets are almost ideally expanded. In all cases, a stagnation point is visible on the flat plate at y = 0 and boundary layers appear on the plate on each side of the jet. For JetL3, JetL5, JetL8 and JetL9, the maximum velocity and the thickness of the boundary layer on the flat plate at y = 2h are given in table 4. The maximum velocity of the boundary layer decreases and the boundary-layer thickness increases with the nozzle-to-plate distance. For JetL5, JetL5-75d and JetL5-60d, the developments of the top and bottom boundary layers are different. The maximum velocity and the thickness of the boundary layer at y = ±2 are given in table 5. When the angle of impact is not normal, the top boundary layer is thicker and its maximum velocity increases. Moreover, shock cells are visible in the upper boundary layer for JetL5-60d in figure 4(c). On the contrary, the maximum velocity and the thickness of the bottom boundary layer decrease when the angle of impact deviates from 90 degrees.

III.B. Flow snapshots

Snapshots of the vorticity norm obtained for JetL3, JetL5, JetL8 and JetL9 in the planes (x, y) and (x, z) along the lipline, are presented in figure 5. For JetL3, in figure 5(a), the two mixing layers exhibit large-scale structures, of typical size 0.5h, in the (x, y) plane, which do not seem to be organized symmetrically with respect to the jet axis. In the spanwise direction, both small and large structures can be seen, namely tubes of vorticity along the entire spanwise extent and small vortical structures of typical size 0.05h. Besides, the spanwise correlation of the vortical structures appears strong, near the nozzle but also close to the flat plate. For jetL5, JetL8 and JetL9, as the nozzle-to-plate distance increases, the asymmetric organisation of the large-scale structures in the two shear layers is still visible but its spanwise correlation seems to decrease in the downstream direction. These observations are in agreement with the experimental results provided by Thurow et al. [START_REF] Thurow | Structure of a supersonic impinging rectangular jet via real-time optical diagnostics[END_REF] for a rectangular supersonic jet impinging on a flat plate normally. Very large organized structures in the shear layers and an asymmetric pattern were indeed noted by the authors using real-time flow visualisation techniques. Snapshots of the vorticity norm obtained for JetL5, JetL5-75d and JetL5-60d in the planes (x, y) are presented in figure 6. The asymmetric organisation of the large-scale structures in the two shear layers, visible for the case JetL5 in figure 6(a), is less apparent as the angle of impact deviates from 90 degrees. Indeed, the pattern may still be noticed in figure 6(b) for JetL5-75d but can not for JetL5-60d in figure 6(c). Moreover, for JetL5-60d, the shear layers seem not to interact between each other. In order to visualize the flow and acoustic fields of the jets, snapshots in the (x, y) plane of the density and the fluctuating pressure are provided in figure 7 for the four cases where the angle between the jet axis and the flat plate is 90 degrees and in figure 8 for the cases JetL5, JetL5-75d and JetL5-60d. Isolated snapshots need to be considered with some cautions but they can provide qualitative results. For all cases, large-scale structures in the shear layers of the jet and sound waves are both observed. In figure 7, the sound waves seem to come from the region of impact of the jet and to have a higher amplitude for JetL3 and JetL5 than for JetL8 and JetL9. In figure 8, sound waves also appear to be generated from the region of impact, with amplitudes decreasing significantly as the angle of impact deviates from 90 degrees. Finally, there seems to be a correlation between the degree of organization of structures in the shear layers, noted in figures 5 and 6, and the amplitude of sound waves, observed figures 7 and 8. A similar trend was found experimentally by Thurow et al.

y = 2h y = -2h U max (m.s -1 ) δ 99 /h U max (m.s -1 ) δ 99 /h

III.C. velocity fluctuations

The variations between the nozzle lips and the flat plate of the peak rms values of axial and transverse velocity fluctuations u ′ and v ′ obtained for JetL3, JetL5, JetL8 and JetL9 are shown in figure 9. For JetL3 and JetL5, the maxima of the peak rms values of axial velocity are found in figure 9(a) around x = 0.5h, and are equal to 31.8% and 32.7% of the jet exit velocity, respectively. For JetL8 and JetL9, the maxima of the peak rms values of axial velocity are located farther downstream, around x = 0.75h, and are equal to 26.6% and 27.3% of the jet exit velocity, respectively. For the four cases, other high peak rms values of axial velocity are also reached just before the plate. They decrease with the nozzle-to-plate distance, as they are equal to 28.5%, 27.3%, 24.4% and 24% of the jet exit velocity for JetL3, JetL5, JetL8 and JetL9. For the peak rms transverse velocity values in figure 9(b), it can be noted that higher values are observed for JetL3 and JetL5 than for JetL8 and JetL9. The nominal values, reached at x = 1.7h in all cases, are around 21% of the jet exit velocity for JetL3 and JetL5, but 17.5% for JetL8 and JetL9.

x/h < v ′ v ′ > 1/2 /u e

III.D. Convection velocity

The evolution of the large-scale structures in the shear layers and more precisely their convection velocity are of primary concern to understand the feedback mechanism between the nozzle and the flat plate. For two pressure-matched parallel streams with equal specific heats, the convection Mach number M c was defined by Papamoschou and Roshko [START_REF] Papamoschou | The compressible turbulent shear layer: an experimental study[END_REF] as

M c = u 1 -u 2 a 1 + a 2 = u 1 -u c a 1 = u c -u 2 a 2 (1)
where u 1 and u 2 are the high-speed and low-speed free stream velocities, a 1 and a 2 are the speeds of sound, and u c is the theoretical isentropic convection velocity. For the present planar jets, u 1 is the jet velocity at the nozzle exit u e , u 2 is the velocity outside of the jet, hence u 2 = 0, a 1 is the speed of sound inside the jet a e = √ γrT e and a 2 is the speed of sound outside of the jet a 0 = √ γrT amb where T e and T amb are the jet and the ambient temperatures. Thus, the convective velocity of the planar jet can be approximated by

u c = u e a e /a 0 + 1 (2)
For the present jets, a convection velocity of u c /u e = 0.57 is thus obtained. The convection velocity u c of turbulent structures along the lipline y = -h/2 was also calculated from the cross-correlations of axial velocity fluctuations u ′ . In practice, for a reference point, cross correlations with twenty regularlyspaced neighbours are computed and the time lag between the reference point and its neighbours is fitted by a straight line. The derivative of the line gives then the local convective velocity. Results are shown in figure 10(a) for JetL3, JetL5, JetL8 and JetL9 and in figure 10(b) for JetL5, JetL5-75d and JetL5-60d. The algorithm does not permit to have the convection velocity close to the nozzle lips and close the flat plate. In figure 10(a), for all cases, the convection velocities are of about 0.50u e at x = h. Farther downstream, they increase slightly until it reaches the location 2h upstream from the flat plate. A rapid grow just before the plate is then visible. The maximal convection velocities between x = h and the location 2h upstream from the flat plate are equal to 0.62u e , 0.64u e , 0.67u e and 0.68u e for JetL3, JetL5, JetL8 and JetL9, respectively. In figure 10(b), the convection velocities for JetL5 and JetL5-75d show the same variations with a maximal value of 0.64u e located around x = 3h whereas the convection velocity of JetL5-60d reaches a peak value of 0.69u e at x = 2.5h. Overall, the mean convection velocity of the large-scale structures along the lipline is of about 0.60u e . This result is in good agreement with the theoretical value of u c /u e = 0.57 and with the experimental results of Panda et al. [START_REF] Panda | Underexpanded screeching jets from circular, rectangular and elliptic nozzles[END_REF] who measured a convection velocity of 0.60u e for an ideally expanded rectangular supersonic jet with an exit Mach number of M e = 1.3.

III.E. Acoustic results

Sound pressure levels (SPL) obtained for JetL3, JetL5, JetL8 and JetL9 at x = 0 and y = 1.5h are displayed in figure 11 as a function of the Strouhal number St h = f h/u e . Several tones emerge, especially for JetL3 and JetL5 for which about ten tones are visible. A large number of tones were also observed experimentally by Norum 10 for a supersonic rectangular jet impinging on a flat plate normally. The frequency of the tones, whose level are about 5 dB higher than the broadband noise, are given in table 6.

In the sound spectrum of JetL3 shown in figure 11(a), the tone frequencies are not the harmonics of one fundamental frequency. Instead, they can be given by the linear combination of the tone frequencies at Strouhal numbers St 1 = 0.115 and St 2 = 0.255, namely the experimental data of Norum. 10 Several tone frequencies were almost always obtained, and for some experimental conditions, two tone frequencies and their linear combinations were observed. Moreover, as detailed earlier, the jets parameters are similar to those of the experimental study of Thurow et al. [START_REF] Thurow | Structure of a supersonic impinging rectangular jet via real-time optical diagnostics[END_REF] and they found for L = 3.94 the seven tone frequencies observed in the present study. Eight other tone frequencies were also obtained by Thurow et al. [START_REF] Thurow | Structure of a supersonic impinging rectangular jet via real-time optical diagnostics[END_REF] and can be due to the rectangular nozzle in the experiment compared to the ideal planar present jets. For JetL5 in figure 11 in each case. However, for nozzle distances of 8.27h and 9.1h, the experimental jet of Thurow can not be seen as planar any more and differences are expected. It can also be noted that in the two latter cases,

St 3 = St 1 + St 2 , St 4 = 2St 1 + St 2 , St 5 = St 1 + 2St 2 , St 6 = 2St 1 + 2St 2 and St 7 = 3St 1 + 2St 2 . A
St 4 = St 1 + St 2 .
Finally, the case JetL5 is the most resonant of our four cases where the angle of impact is θ = 90 degrees with a maximum sound pressure level of 188 dB/St. It was also the most resonant case in the experiment of Thurow et al. [START_REF] Thurow | Structure of a supersonic impinging rectangular jet via real-time optical diagnostics[END_REF] Sound pressure levels acquired for JetL5, JetL5-75d and JetL5-60d at x = 0 and y = 1.5h are shown in figure 12. As earlier, several tone frequencies are noted and those whose levels are about 5 dB higher than the broadband noise are given in table 7. For JetL5-75d, in figure 7(b), two emerging tone frequencies are visible at St 1 = 0.19 and St 2 = 0.25. The first one is the same one as for the case JetL5 in figure 12(a) and the latter one do not exist in JetL5. For JetL5-60d, in figure 12(c), there is no clear tone frequency. However, the maximum sound pressure level seems to be found around St = 0.25.

III.F. Skewness and kurtosis factors

In order to investigate the non-linearity of the acoustic waves generated by the impingement, the properties of the pressure has been examined at x = 0 and y = 8.5h. Figure 13(a) displays the fluctuating pressure at this point for JetL5. Weak shock waves and N-shaped waves with sharp compressions associated with gradual expansions are visible. These shapes are similar to those obtained by de Cacqueray and Bogey [START_REF] De Cacqueray | Noise of an overexpanded Mach 3.3 jet: non-linear propagation effects and correlations with flow[END_REF] in the acoustic field of an overexpanded jet 24 at M j = 3.3. Moreover, for some times, the amplitude exceeds 0.15P amb . This is the case around t = 250h/u e . A zoom centered around this time is presented in figure13(b). Two strong shock waves are found at t 1 = 242h/u e and t 2 = 253h/u e . In order to quantify this observation, figure 13(c) shows the probability density function (PDF) of the fluctuating pressure normalized by the standard deviation as well as the skewness and kurtosis factors. A skewness factor of S = 0.59 and a kurtosis factor of K = 3.99 are obtained, indicating a strong non-linearity and intermittency in the signal, respectively. This behaviour with strong intermittent shock waves seems to correspond to that of the crackle noise first defined by Ffowcs-Williams et al., [START_REF] Ffowcs-Williams | Crackle: An annoying component of jet noise[END_REF] who obtained this specific noise when the skewness factor exceeds 0.4. Moreover, the evolution visible in figure 13(b) agrees well with the results of Baars and Tinney. [START_REF] Baars | Quantifying crackle-inducing acoustic shock-structures emitted by a fully-expanded Mach 3 jet[END_REF] In an experimental Mach 3 jet, they found out that strong shock waves, responsible for the crackle noise, are often present in groups of multiple shocks. For the six jets studied, the maxima of sound pressure levels shown in figures 11 and 12, and the skewness and the kurtosis factors at x = 0 and y = 8.5h, are reported in table 8. A significant correlation appears between the maxima in the sound spectra of figures 11 and 12, and the skewness factors of the fluctuating pressure at x = 0 and y = 8.5h. Moreover, the most resonant case, JetL5, show the strongest intermittency of the fluctuating pressure with a kurtosis factor of 3.99.

IV. Aeroacoustic feedback mechanism

IV.A. Tone frequencies

In order to explain the tone frequencies found experimentally, Powell 1 was the first to suggest that a feedback mechanism establishes between the nozzle lips and the flat plate. This mechanism consists of two steps. First, in the shear layer, a coherent structure is convected downstream from the nozzle to the plate. The structure impinges on the plate, and generates an acoustic wave propagating upstream towards the nozzle. This wave is then reflected by the nozzle lip, excites the shear layer, and leads to the formation of a new coherent structure. The fundamental period T 0 of this feedback loop corresponds to the sum of the time necessary for a coherent structure to travel downstream from the nozzle to the plate and the time of propagation of an acoustic wave from the plate to the nozzle, yielding

T 0 = L 0 1 u c (x) dx + L c 0 = N + p f (3)
where u c (x) is the local convection velocity of the large-scale structures in the shear layers, c 0 is the speed of sound outside of the jet, p is a phase lag at the nozzle exit, and the mode number N indicates the number of time the feedback mechanism occurs during the fundamental period T 0 . Powell 1 argued that the phase lag p is not necessarily null because the reflection of the acoustic wave on the nozzle lips and the creation of a coherent structure are not simultaneous phenomena. In experiments, Krothapalli et al. [START_REF] Krothapalli | Flow field and noise characteristics of a supersonic impinging jet[END_REF] found for instance p = 0 for subsonic round jets and p = -0.4 for supersonic round jets. Later, Ho and Nosseir 2 and Nosseir and Ho 3 investigated the near field pressure of jets impinging on a flat plate normally, and proposed an empirical model in order to predict the frequencies of the feedback mechanism, this model writes

N f = L u c + L c 0 (4)
where u c is the mean convective velocity of the large-scale structures in the shear layers. Here, the phase lag p is considered null. The Strouhal number of the tone frequencies obtained for the present jets in figures 11 and 12 are shown in figure 14 as a function of the nozzle-to-plate distance L/h. Only the fundamental tone frequencies, which are not harmonics of other tone frequencies and which are called source tone frequencies in the following, are represented. The tone frequencies observed in the experiments of Thurow et al., [START_REF] Thurow | Structure of a supersonic impinging rectangular jet via real-time optical diagnostics[END_REF] as well as the tone frequencies predicted by relation ( 4) with a mean value of convective velocity u c = 0.60u e , are also plotted. As mentioned previously, a good agreement is found between the simulation and the experimental data. Moreover, the tone frequencies seem to be well predicted by the empirical model. for JetL5-75d, and × experimentally by Thurow et al. [START_REF] Thurow | Structure of a supersonic impinging rectangular jet via real-time optical diagnostics[END_REF] The grey lines show the values predicted by the model (4) using uc = 0.60ue.

For JetL3, with L/h = 3.94, the two source tone frequencies, at St 1 = 0.115 and St 2 = 0.255, correspond to the first mode and the third mode predicted by the model. For JetL5, with L/h = 5.5, the three source tone frequencies, at St 1 = 0.12, St 2 = 0.19 and St 3 = 0.33, can be associated with the second mode, the third mode and the fifth mode. For JetL8 and JetL9, with L/h = 8.27 and L/h = 9.1, the three source tone frequencies can be connected with the second mode, the fourth mode and the fifth mode.

The influence of the angle of impact can be explored by comparing JetL5, JetL5-75d and JetL5-60d. For JetL5, the dominant source tone St 2 = 0.19 is linked to the third mode of the model (4). For JetL5-75d, the two observale tones, at St 1 = 0.19 and St 2 = 0.25, correspond to the third and the fourth modes of the model. Finally, for JetL5-60d, no tones are represented in the figure 14, because they are not clearly visible in the spectra. The maximum sound pressure level is however reached around St = 0.25. Thus, when the angle of impact deviates from 90 degrees, there seems to be an evolution from the third mode to the fourth mode predicted by (4). A similar jump is observed in figure 14 when the nozzle-to-plate distance increases from 5.5h to 5.8h in the experiment of Thurow et al. [START_REF] Thurow | Structure of a supersonic impinging rectangular jet via real-time optical diagnostics[END_REF]

IV.B. Fourier decomposition of the pressure field

For each jet, the pressure field in the (x, y) plane has been recorded every 50 th time steps. The results are arranged in a M × N matrix P all :

P all =      P 1 1 P 2 1 . P N 1 P 1 2 P 2 2 . P N 2 . . . . P 1 M P 2 M . P N M      (5)
where N is the number of samplings, and M = n x × n y is the total number of points. The pressure field obtained at a given time is thus provided by one column of P all . A Fast Fourier Transform is applied to each row of the matrix P all . In this way, for a given frequency, amplitude and phase fields can be displayed. For the four jets with θ = 90 degrees, the dominant tone frequencies are considered. The amplitude fields are shown in the top views of figure 15, and the phase fields are displayed in the bottom views. The latter fields enable to identify the nature of the modes. The dominant source tone frequencies in JetL3 and JetL5 correspond to antisymmetric modes because there is a 180 degree phase shift visible with respect to the jet axis in figures 15(e) and 15(f). On the contrary, the dominant source tone frequencies in JetL8 and JetL9 are associated with a symmetric mode, as the phase shift is null with respect to the jet axis in figures 15(g) and 15(h). For planar high subsonic jets impinging on a flat plate, Arthurs and Ziada 9 also found that the source tone frequencies obtained in the spectra are generated by both symmetric and antisymmetric modes.

Besides, these results agree well with the theoretical developments of Tam and Norum 11 who proposed that the upstream propagating waves of the feedback mechanism are due to neutral acoustic wave modes of the jet. Such modes are characterized by a wave number k and an angular frequency ω both real. For a supersonic jet, the existence of acoustic modes propagating upstream was demonstrated by Tam et Hu. [START_REF] Tam | On the three families of instability waves of high-speed jets[END_REF] Such waves propagate upstream along the shear layer of the jet. Tam and Norum 11 used a vortex sheet jet model to predict these modes and found an allowable frequency range for each upstream propagating neutral acoustic wave mode of the jet flow. The two dominant tone frequencies of JetL3 and JetL5, which are associated to symmetric modes, fall into the first symmetric allowable band given by this model. Moreover, the two dominant tone frequencies of JetL8 and JetL9, linked to antisymmetric modes, take place into the first antisymmetric allowable band given by this model. The experimental results of Norum [START_REF] Norum | Supersonic rectangular jet impingement noise experiments[END_REF] were also in agreement with this model.

The amplitude fields of the dominant tone frequencies for JetL3, JetL5, JetL8 and JetL9, represented in the top of figure 15, all exhibit a cell structure between the nozzle lips and the flat plate. By considering the two semi-cells near the nozzle lips and the plate as one cell, the cell structure contain three cells for JetL3 and JetL5 in figures 15(a) and 15(b) and two cells for JetL8 and JetL9 in figures 15(c) and 15(d). Thus, the number of cells correspond to the mode number predicted by the model of Ho and Nosseir. 2 Indeed, the third mode was identified for the dominant source tone frequencies of the two first jets, and the second mode for the two others. Moreover, for the present jets, the amplitude fields of all tone frequencies have been examined and the number of visible cells corresponds in all cases to the mode number given by the model of Ho and Nosseir. Information on the sound sources are also provided by the amplitude an phase fields of figure 15. For a better description of these fields, let α be the angle at the impingement region between the upstream direction and the waves propagating from the flat plate, illustrated in figure 15(a).

For the dominant tone frequency of JetL3, in figures 15(a) and 15(e), three acoustic contributions can be seen by analysing the intense spots and the destructive interferences in the amplitude field and by looking at the phase contours in the phase field. Indeed, the phase contours in the near acoustic field appear to spread concentrically from three regions. The first visible acoustic contribution can be seen in the shear layers of the jet and in the near field on each side of the nozzle, for 0 < α < 30 degrees and correspond to the upstream propagating neutral acoustic wave mode of the jet. The second contribution is visible for 30 < α < 50 degrees in figure 15(e) as the phase contours in this region seem to spread from a point around x = 2.5h and y = 1.5h. This position coincides with an intense spot in the magnitude field in figure 15(a). Finally, the third acoustic contribution is visible for α > 50 degrees in figure 15(e) and seems to propagate from a point on the plate at y = 3h. Again, this position matches an intense spot in the magnitude field displayed in figure 15(a). Finally, destructive interferences can be seen for α = 30 and α = 50 degrees between the three acoustic contributions in figure 15(e) with a strong reduction of the amplitude in these directions. The same reasoning can be applied to the dominant source tone frequency of JetL5, in figures 15(b) and 15(f). In this way, two acoustic contributions can be observed. The first one is still due to the upstream propagating neutral acoustic wave mode of the jet and can be seen for 0 < α < 30 degrees and the second one seems to originate from a point located on the plate at y = 3.5h. Destructive interferences can also be seen between the two acoustic contributions for α = 40 degrees. For the dominant tone frequencies of JetL8 and JetL9, two acoustic contributions seem to emerge. One still due to the upstream propagating neutral acoustic wave mode of the jet and another propagating from a point located close to the plate. However, the large wavelengths of these modes compared to the computational domain makes it difficult to look at the phase contours or to see destructive interferences.

The amplitude and phase fields obtained for the dominant source tone frequency of JetL5 and the two tone frequencies of JetL5-75d are displayed in figure 16. For JetL5-75d, the amplitude fields in figures 16(b) and 16(c) exhibit a cell structure, as in figure 16(a) for the dominant tone frequency of JetL5. There are three and four cells visible in figures 16(b) and 16(c), respectively, corresponding to the third and the fourth modes of the model of Ho and Nosseir. 2 This support the hypothesis of a switch from the third mode to the fourth mode, mentioned previously, as the angle of impact deviates from 90 degrees.

IV.C. Hydrodynamic-acoustic standing wave

The cell structures noted in the amplitude fields of figures 15 and 16 may be due to the presence of an hydrodynamic-acoustic standing wave. Indeed, when a feedback mechanism establishes in a flow, such waves can be obtained. In the experiment of Panda et al. [START_REF] Panda | Underexpanded screeching jets from circular, rectangular and elliptic nozzles[END_REF] on a screeched supersonic jet, an hydrodynamic-acoustic standing wave was for instance observed, corresponding to the sum of hydrodynamic coherent structures in the shear layers propagating downstream and acoustic waves propagating upstream.

A simple model can be built to determine the wave number of the hydrodynamic-acoustic standing waves in the present jets. Let ω be the angular frequency of the feedback mechanism and k p et k a the wave numbers of the hydrodynamic coherent structures and of the acoustic waves, respectively. In the shear layer, the resulting fluctuating pressure can be written as

p ′ (x, t) = A a sin(k a x + ωt + φ a ) + A p sin(k p x -ωt + φ p ) (6)
where A a , A p , φ a and φ p are the amplitudes and phase shifts of the acoustic waves and of the hydrodynamic structures. Then, the quadratic mean value of the fluctuating pressure is

p ′ rms (x) = 1 T T 0 p ′ (x, t) 2 dt = A 2 a 2 + A 2 p 2 -A a A p cos ((k a + k p )x + φ a + φ p ) (7)
where T = 2π/ω is the time period. The wave number k sw of the hydrodynamic-acoustic standing wave is thus equal to k p + k a . By noting that L sw = 2π/k sw is the wavelength of the hydrodynamic-acoustic standing wave, and that k p = ω/u c and k a = ω/c 0 , one obtains

1 L sw = f c 0 + f u c (8)
where f = 2π/ω, c 0 is the speed of sound outside of the jet and u c is the mean convective velocity of large-scale structures in the shear layers.

As seen in figures 15 and 16, there is a whole number N of cells between the nozzle lips and the flat plate. Considering that the cell structure obtained in the amplitude fields in figures 15 and 16 corresponds to an hydrodynamic-acoustic standing wave, the nozzle-to-plate distance L is a multiple of the wavelength, yielding L sw = L/N . Then, the relation (8) becomes

N f = L c 0 + L u c (9)
The relation ( 9) is identical to the relation (4) given by the model of Ho and Nosseir, 2 indicating that the number of cells of the standing wave is also the mode number. Thus the aeroacoustic feedback mechanism occurring between the nozzle lips and the flat plate corresponds to the formation of an hydrodynamic-acoustic standing wave as defined by Panda et al. [START_REF] Panda | Underexpanded screeching jets from circular, rectangular and elliptic nozzles[END_REF]

V. Conclusion

In this paper, the hydrodynamic and acoustic properties of six planar supersonic jets computed by compressible large-eddy simulations (LES) using low-dissipation schemes are presented. The jets are ideally expanded, and have a Mach number of M e = 1.28, and a Reynolds number of Re h = 5 × 10 4 . Four jets impinging on a flat plate normally, and two other impinging jets with angles of impact of θ = 60 and θ = 75 degrees are considered. Mean velocity flows and snapshots of density, fluctuating pressure and vorticity are described. The results agree well with experimental data. The development of the shear layers is investigated by computing rms values of velocity fluctuations and convective velocities of large-scale structures. The mean convective velocities are in line with the theoretical model proposed by Papamoschou and Roshko [START_REF] Papamoschou | The compressible turbulent shear layer: an experimental study[END_REF] and with the experimental values found by Panda et al. [START_REF] Panda | Underexpanded screeching jets from circular, rectangular and elliptic nozzles[END_REF] The near pressure fields are then detailed. First, sound pressure levels are computed. They revealed several tones emerging, which agree well with the measurements of Thurow et al., [START_REF] Thurow | Structure of a supersonic impinging rectangular jet via real-time optical diagnostics[END_REF] as well as with the empirical model of Ho and Nosseir. 2 Finally, the skewness and kurtosis factors of the fluctuating pressure are computed, and properties typical of crackle noise are observed for the most resonant case. A Fourier decomposition is also applied to the near pressure fields in the (x, y) plane. The amplitude and the phase fields of the dominant source tone frequencies for each case are considered. For the cases of normal impact, the antisymmetric or symmetric natures of the modes are identified. The results are in agreement with the theoretical developments of Tam and Norum. [START_REF] Tam | Impingement tones of large aspect ratio supersonic rectangular jets[END_REF] Besides, an hydrodynamic-acoustic standing wave is found for each source tone frequency, as observed by Panda et al. [START_REF] Panda | Underexpanded screeching jets from circular, rectangular and elliptic nozzles[END_REF] on a screeched supersonic jet. The number of cells in the structure corresponds to the mode number of the model given by Ho and Nosseir. [START_REF] Ho | Dynamics of an impinging jet. Part 1. The feedback phenomenon[END_REF] As the angle of impact deviates from 90 degrees, a jump from one mode to a higher mode and a reduction of the intensity of the feedback mechanism are observed.

In further work, the underlying physics of the feedback mechanism will be explored. This will be done by extracting the most energetic modes by proper orthogonal decomposition and by using dynamic mode decomposition. Moreover, the theoretical developments of Tam and Norum [START_REF] Tam | Impingement tones of large aspect ratio supersonic rectangular jets[END_REF] will be applied to the present jets.

5 ×

 5 10 4 9.1 90 degrees JetL5-60d 1.28 5 × 10 4 5.5 60 degrees JetL5-75d 1.28 5 × 10 4 5.5 75 degrees

25

 25

Figure 1 .

 1 Figure 1. Representation of the two Cartesian meshes used in the case JetL5-15d: (a) meshes sketch of the two meshes, with the primary grid in light grey, the secondary grid in dark grey and the nozzle and flat plate in black and (b) representation of the grids, every fifteen points is shown.

Figure 2 .

 2 Figure 2. Representation of (a) the lateral mesh spacings, and (b) the axial mesh spacings: JetL3, . JetL5, JetL8, and JetL9.

Figure 3 .Table 4 .

 34 Figure 3. Mean velocity field in the (x, y) plane for (a) JetL3, (b) JetL5, (c) JetL8 and (d) JetL9. The colour scale ranges from 0 to 400 m.s -1 . The nozzle is represented in black.

Figure 4 .

 4 Figure 4. Mean velocity field in the (x, y) plane for (a) JetL5, (b) JetL5-75d and (c) JetL5-60d. The colour scale ranges from 0 to 400 m.s -1 . The nozzle and the flat plate are represented in black.

Figure 5 .

 5 Figure 5. Snapshots in the planes (x, y) and (x, z) along the bottom lipline of the vorticity norm |ω| for (a) JetL3, (b) JetL5, (c) JetL8 and (d) JetL9. The colour scale ranges up to the level of 10ue/h. The nozzle is represented in black.

Figure 6 .

 6 Figure 6. Snapshots in the (x, y) plane of the vorticity norm |ω| for (a) JetL5, (b) JetL5-75d and (c) JetL5-60d. The colour scale ranges up to the level of 10ue/h. The nozzle and the plate are represented in black.

Figure 7 .

 7 Figure 7. Snapshots in the (x, y) plane of the density in the jet axis and close to the flat plate and of the pressure fluctuations for (a) JetL3, (b) JetL5, (c) JetL8 and (d) JetL9. The colour scale ranges from 1 to 2 kg.m -3 for the density and from -7500 to 7500 Pa for the fluctuating pressure. The nozzle and the flat plate are represented in black.

13

 13

Figure 8 .

 8 Figure 8. Snapshots in the (x, y) plane of the density in the jet axis and close to the flat plate and of the pressure fluctuations for (a) JetL5, (b) JetL5-75d and (c) JetL5-60d. The colour scale ranges from 1 to 2 kg.m -3 for the density and from -7500 to 7500 Pa for the fluctuating pressure. The nozzle is represented in black.

Figure 9 .

 9 Figure 9. Variations of peak rms values of velocity fluctuations (a) u ′ and (b) v ′ for JetL3, . JetL5, JetL8, and JetL9.

Figure 10 .

 10 Figure 10. Convection velocity of the large-scale structures in the shear layers for (a) JetL3, . JetL5, JetL8, and JetL9 and (b) JetL5, JetL5-75d, and . JetL5-60d.

Figure 11 .

 11 Figure 11. Sound Pressure Levels (SPL) at x = 0 and y = 1.5h as a function of the Strouhal number for (a) JetL3, (b) JetL5, (c) JetL8 and (d) JetL9.

5 = 2 +

 52 (b), a fundamental tone frequency at St = 0.19 and its first six harmonics can be seen. Two other tone frequencies St 1 = 0.12 and St 3 = 0.33 are visible and St St St 3 . The agreement is excellent with the experimental data of Thurow et al.[START_REF] Thurow | Structure of a supersonic impinging rectangular jet via real-time optical diagnostics[END_REF] who obtained a fundamental tone frequency occurring at St exp = 0.20. For JetL8 and JetL9 in figures 11(c) and 11(d), four tone frequencies are visible and only the third one St 3 is observed experimentally by Thurow et al.[START_REF] Thurow | Structure of a supersonic impinging rectangular jet via real-time optical diagnostics[END_REF]

Figure 12 .

 12 Figure 12. Sound Pressure Levels (SPL) at x = 0 and y = 1.5h as a function of the Strouhal number for (a) JetL5, (b) JetL5-75d and (c) JetL5-60d.

Figure 13 .

 13 Figure 13. Fluctuating pressure at x = 0 and y = 8.5h for JetL5: (a) time variations over the entire simulation, (b) variations over a short period and (c) probability density function ( probability density function of a Gaussian distribution).

Figure 14 .

 14 Figure 14. Strouhal number of the tone frequencies obtained: • for JetL3, JetL5, JetL8 and JetL9,for JetL5-75d, and × experimentally by Thurow et al.[START_REF] Thurow | Structure of a supersonic impinging rectangular jet via real-time optical diagnostics[END_REF] The grey lines show the values predicted by the model (4) using uc = 0.60ue.

Figure 15 .

 15 Figure 15. Amplitude (top) and phase (bottom) fields obtained for the dominant tone frequency in (a) and (e) JetL3 (Sth = 0.255), (b) and (f ) JetL5 (Sth = 0.19), (c) and (g) JetL8 (Sth = 0.092) and (d) and (h) JetL9 (Sth = 0.085).

2

 2

Figure 16 .

 16 Figure 16. Amplitude (top) and phase (bottom) obtained for the pressure fields for (a) and (d) JetL5 at Sth = 0.19, JetL5-75d (b) and (e) at Sth = 0.19 and (c) and (f ) at Sth = 0.25.

Table 1 .

 1 Jet parameters: Mach number Me, Reynolds number Reh = ueh/ν, nozzle-to-plate distance L, and angle θ between the jet direction and the flat plate.

Table 2 .

 2 Parameters of the primary grids containing the jets.

	n x	n y	n z number of points
	JetL3 JetL5 JetL8 JetL9 JetL5-60d 1196 801 219 799 1051 219 903 1051 219 1087 1051 219 1142 1051 219 JetL5-75d 966 901 219	184 × 10 6 208 × 10 6 250 × 10 6 263 × 10 6 210 × 10 6 191 × 10 6
	n x	n y	n z number of points
	JetL5-60d 181 1151 219 JetL5-75d 181 961 219	46 × 10 6 38 × 10 6

Table 3 .

 3 Parameters of the secondary grids close to the plate, used for the cases with inclined plates.

Table 5 .

 5 Maximum velocity and thickness of the boundary layer on the flat plate at y = ±2h.

Table 6 .

 6 Strouhal numbers emerging in the spectra of figure11. The Strouhal number of the main tone of each case appear in bold.

		St 1	St 2	St 3	S 4	St 5	St 6	St 7
	JetL3 0.115 0.255 0.37 0.485 0.625 0.74 0.855
	JetL5	0.12	0.19 0.33 0.38	0.52 0.57 0.76
	JetL8 0.092 0.165 0.21 0.255 JetL9 0.085 0.145 0.19 0.23	--	--	--

Table 7 .

 7 Strouhal numbers emerging in spectra of figure12. The Strouhal number of the main tone of each case appear bold. The Strouhal number of the maximum sound pressure level for JetL5-60d appears in grey.

		St 1	St 2	St 3	S 4	St 5	St 6	St 7
	JetL5	0.12 0.19 0.33 0.38 0.52 0.57 0.76
	JetL5-75d 0.19 0.25 JetL5-60d 0.25 -	--	--	--	--	--

Table 8

 8

		JetL3 JetL5 JetL8 JetL9 JetL5-75d JetL5-60d
	SP L max (dB/S t )	185	188	182	183	175	160
	Skewness	0.39	0.59	0.27	0.29	0.25	0.21
	Kurtosis	3.16	3.99	3.36	3.25	3.17	3.42

. Maxima levels in the sound spectra of figures 11 and 12, skewness and kurtosis factors of the fluctuating pressure at x = 0 and y = 8.5h.

Acknowledgments

This work was performed using HPC resources of P2CHPD (Pôle de Calcul Hautes Performances Dédiés) and IDRIS (Institut du Développement et des Ressources en Informatique Scientifique) under the allocation 2014-2a0204 made by GENCI (Grand Equipement National de Calcul Intensif). This work was performed within the framework of the Labex CeLyA of Université de Lyon, within the program "Investissements d'Avenir" (ANR-10-LABX-0060/ ANR-11-IDEX-0007) operated by the French National Research Agency (ANR).


Claude Cosandey
email: cosandey@cnrs-bellevue.fr

Claude Martin

Lionel Savina

Jean-François F Didon-Lescot

The impact of interception losses on the water balance in forested mountains range

Keywords: Forest hydrology, interception, forest evaporation. 1.INTRODUCTION

come   L'archive ouverte pluridisciplinaire

The hydrological impact of forest is one of the most controversial issues when it comes to studying the consequences of human activities on water resources (Calder, 1979(Calder, , 1985;;Morton, 1984Morton, , 1985)). As early as 1982 -and the study is still a reference -Bosch and Hewlett concluded their study "A review of catchment experiments to determine the effect of vegetation changes on water yield and evaporation" by emphasizing the complex nature of the results. Only one thing was certain: "No experiment in deliberately reducing cover caused reductions in yield, nor have any deliberate increases in cover caused increases in yield". Obviously, then, if the forest has an impact on annual runoff, it is to reduce it. The reduction is all the greater when conditions include a large water deficit, abundant water reserves, and also frequent but light rainfall (Cosandey and Robinson, 2000).

Two mechanisms account for the reduction:

• One, the forest affords a larger surface for interception during rainfall and enough roughness to favour high air turbulence and, therefore, more efficient use of advective energy and a higher rate of evaporation. On the other hand, the evaporation can limit plant transpiration, and it is difficult to evaluate the increase in overall evapotranspiration that results from the direct evaporation of intercepted water. It is not possible to know the value for the increase from direct measurements of the intercepted water that does not reach the soil, since a partial compensation can occur with a decrease in plant transpiration (the available energy is used more quickly in order to evaporate the water that is more easily available on the surfaces of leaves or branches). • Two, because the forest is more deeply rooted, it has a greater potential soil moisture reserve (which secures its water supply when evaporation exceeds rainfall). A greater reserve allows more evaporation during the summer months, and therefore lower flows in autumn that start later in the season due to the larger amount of water taken up by the soil.

If a classic sequence of rain events for an hydrologic year in temperate climate is considered (Cosandey and Robinson, 2000), the forest delays the resumption of flows at the start of the rainy season, because the soil moisture reserves have been depleted by the trees during the dry season. Of course, the delay is accompanied by a drop in runoff. The forest also reduces winter flows directly, due to the direct evaporation of part of the precipitation that is intercepted by the canopy and does not reach the soil. Such higher winter 2003 -Interdisciplinary approaches in small catchment hydrology : Monitoring and research, 9 th ERB Conference (Demänovská dolina, 2002), L. HOLKO et P. MIKLANED édit., TDH in Hydrology, n° 67, UNESCO, p. 79-84. evaporation depends a lot on the type of vegetation, of course, and a heath with broom growing on it is probably much more efficient than a broad leaved forest.

Evaporation during the summer period (when evaporation demand exceeds rainfall) is dependent on the available water (precipitation plus the soil water storage). Interception has no impact on precipitation, none on soil water storage, and so it does not change the values for evapotranspiration during the summer. The resulting interception losses, therefore, have an effect on flows only during the winter period and can be estimated only within that framework.

When a forest is cut, the estimation of evaporation before and after the cut should make it possible to estimate the differences between a forested catchment and an unforested one, and therefore the impact of interception losses on runoff.

2.EXPERIMENTAL SET-UP

Three small catchments on granite in the Mount Lozère Experimental Research Basin (ERB) have been monitored since 1981. These are the spruce-forested Latte (0.195 km 2 ), the beech-forested Sapine (0.54 km 2 ), and the grazed-grassland Cloutasses (0.81 km 2 ). The slopes there are moderately steepapproximately 12° (Latte), 18° (Sapine), and 10° (Cloutasses). The Latte catchment (Fig. 1) was the most closely observed catchment, as 80% of its surface was clear-cut of its spruce forest from 1987 to 1989. At between 1100 and 1500 metres in elevation, the Mount Lozère ERB has a mediterranean climate with mountain characteristics. The mean temperature at 1300 m is of 6.9°C. Mean annual precipitation is about 2000 mm, and ranges between 1100 and 3500 mm. Rainfall can be very violent, especially in autumn during "cévenols" events; for example on the Latte catchment, maximum intensities for 30 minutes reached 179 mm h -1 on August 28 1999 and 131 mm h -1 on September 22 1993. On the average, the soils and superficial deposits are from 60 cm (Sapine) to 70 cm (Latte and Cloutasses) thick. Filtration rates for the soils as determined in simulated rain conditions range from 78 to more than 123 mm h -1 under undisturbed vegetation and for well-protected soils (Cosandey et al., 1990).

3.METHODOLOGY

Estimation of "interceptions losses" (increase in global evapotranspiration resulting from interception) can be done only from water balance. Direct measurement of interception allows knowing the part of the rain that never reaches the soil. But the energy used for evaporated this intercepted water is not able for transpiration or other evaporation, and the global evaporation is not accrues by the total value of measured interception.

The method is based on the water balance equation calculated during the "excess water period", when rainfall exceeds atmospheric water demand (what is called "hydrologic winter", usually -but not alwaysfrom October to April).

It is known that during the hydrologic winter (defined by R>Pe), the equation may be expressed as: In the following study, the thin soils prevent the vegetation from developing a very deep root system, so there is very little difference in the soil moisture storage. This difference has been estimated at 20 mm (100 mm for the forest and 80 mm for heath/grassland). These reserves probably were somewhat depleted when the forest was cut, at least during the two year following the cut. Because these values are a bit arbitrary, they do constitute a source of uncertainty in the following developments. If the estimated interception losses into consideration, however, the magnitude of the uncertainty is negligible.

From this basic equation, you obtain for the duration of the hydrologic winter:

Ae = R -D -∆R
This calculation of Ae is done for the two basins. The problem, well known, is that the value of Ae is the residual term of the calculation, which includes errors in rainfall and runoff storage measurements. For this reason, not the values themselves, but only the ratio between the values are taking in account, according to the "comparative basins method". This method allows to reduce error concerning rainfall (which are quite the same for the two basins, and have no consequence on the difference of Ae estimation) and error concerning runoff, due to calibration curves, which remain the same for each basin respectively. Of course, the forest cut didn't change only evaporation (from the angle of evaporation of intercepted water). It also affected plant transpiration, although we are operating here on the hypothesis that plant transpiration can be ignored during the months of winter dormancy, especially given the low temperatures involved.

4.RESULTS

Differences in Ae, attributed to differences in interception losses were estimated over a period of twelve years for two catchments, (the heath/grassland Cloutasses, and the initially forested Latte, five years before the cur; seven years after the cut during the summer of 1987). The results are shown in the following table (Table 1).

Two preliminary remarks: 1) Values of Actual evapotranspiration are higher than values of Potential evapotranspiration from Turc formula (from local data). Turc formula is not the best for Pe estimation. But there is not data for Penman equation and theses Pe values are taken in account only for the limits of winter season. Even if Turc formula drives to lower values, it seems that interception losses are effective even on the grazed heath/grassland catchment. The presence of broom plants and clumps of trees can probably explains this. 2) The wide range of values is difficult to explain, even if they are clearly in relation with Etp values. It is clear that measurement error both on rainfall and runoff play a role; but as seen above consequences are minor on the differences. In order to determine the impact of the forest cut on interception losses, we can take a look at the evolution in the ratio between Ae in the two basins before and after the cut.

Figure 2 shows that after the cut, the ratio between evaporation in the two basins, with was positive before the cut became negative after.

5.DISCUSSION

Before the cut, the interception losses were higher in the spruce-forested Latte catchment than on the Cloutasses catchment. There was a large decrease in interception losses during the winter that followed the cut, explained by the reduction in evaporating surface area: more than a third of the trees were cut and, although debris from the cut was strewn over the ground, there was no plant colonisation to cover the deforested areas. In these conditions, it is not surprising that interception losses were lower than on the heath/grassland catchment.

The situation rapidly changed, however. By the following summer, an abundant herbaceous (great willow herb) and shrub (raspberry bush) vegetation had developed; nowhere was bare ground to be found. By the second post-cut winter, the deforested catchment's behaviour was similar to that of the heath/grassland catchment, even though cutting operations had resumed in an area where tree growth was less successful; it remained similar even though some very small trees were planted (very slow growth in the difficult climatic conditions, so their impact was still very limited).

6.CONCLUSION

This study was based on a certain number of hypotheses and rough calculations. In particular, it is clear that using water balances involves the risk of accumulating errors in data measurements (principally rainfall and discharge). Comparing neighbouring catchments is a way to reduce this type of risk, however, especially since it makes it easier to take any possible impact of climate variability into account. Beyond such uncertainty, we felt that it was of interest to give further thought to the hydrologic processes that bring about lower flows under forest cover. Once a higher level of evaporation during the winter period was identified on the forested catchment, and it was ascertained that it could be due only to the involvement of loss by interception, it became easier to understand why the overall increase in runoff that was observed after the cut concerned only minor floods. During the "cévenols" rainfall events, the amount of rain and the depth of runoff can be enormous -and the effects of interception losses become totally unnoticeable. The same cannot be said of a relatively light rainfall event. In that case, the depth of runoff measures only a few millimetres, and can therefore be heavily influenced in terms of relative value by a reduction in the effective rainfall due to interception of part of the rain by the plant cover.

REFERENCES:

Calder

Figure 1 -

 1 Figure 1 -Experimental pattern.

R

  evaporation (calculated with Turc formula; Turc, 1961): ∆R = ∆Ru + ∆Rh with ∆Ru soil moisture storage (function of soil characteristics and depth of root development) and ∆Rh groundwater storage (deduced from beginning and end-of-period base flow and the recession curve)

Figure 2 -

 2 Figure 2 -Ratio between Actual evapotranspiration values for the Latte and the Cloutasses Basins.

Table 1 -

 1 Actual evaporation in the two basins for the whole period (see annexe 1 for details):

	hydrologic	Ae	Ae Latte	Ae Lat -Ae	Ae Lat / Ae Clout
	winter	Cloutasses		Clout
	82/83	555	671	-116	1.21
	83/84	349	456	-107	1.31
	84/85	463	572	-109	1.24
	85/86	121	198	-77	1.64
	86/87	348	499	-151	1.43
	87/88	504	394	110	0.78
	88/89	300	283	17	0.94
	89/90	538	498	40	0.93
	90/91	226	198	28	0.88
	91/92	447	444	3	0.99
	92/93	457	475	18	1.04
	93/94	328	347	-19	1.06
	94/95	488	470	18	0.96
	95/96	347	276	71	0.80

  I.R., 1979, Do trees use more water than grass? Water Services, 83: 11-14 Calder I.R., 1985, What are the limits on Forest Evaporation? -Comment. J. Hydrol., 82: 179-192. Cosandey C., Boudjemline D., Roose E., Lelong F., 1990, Etude expérimentale du ruissellement sur des sols à végétation contrastée du Mont Lozère. Z. für géomorphologie, N.F. 34-1: 61-77. Cosandey C. et M. Robinson, 2000, Hydrologie continentale. Coll. U, Armand Colin, Paris, 360 p. Morton F. I., 1984, What are the limits of forest evaporation? J. Hydrol., 74: 373-398 Morton F.I., 1985, What are the limits of forest evaporation? -Reply. J. Hydrol., 82: 184-192 Turc L., 1961, Evaluation des besoins en eau d'irrigation. Formule climatique simplifiée et mise à jour. An. Agro., 12 (1): 13-49.

	Annexe 1 -Data for Ae calculation
			Latte Basin				Cloutasses Basin
	hydrol. winter	D (mm) R (mm) ∆R	Ae	hydrol. winter	D (mm) R (mm) ∆R	Ae
	82/83 1500 2269	98	671	82/83 1513 2184	116	555
	83/84	874	1490	160	456	83/84 1018 1490	123	349
	84/85 1455 2150	123	572	84/85 1452 2052	137	463
	85/86	840	1310	272	198	85/86 1049 1310	140	121
	86/87 1273 1890	118	499	86/87 1221 1738	169	348
	87/88 1959 2608	255	394	87/88 1925 2608	179	504
	88/89	624	1065	158	283	88/89	620	1065	145	300
	89/90	752	1384	134	498	89/90	720	1384	126	538
	90/91	950	1280	132	198	90/91	934	1276	116	226
	91/92	805	1380	131	444	91/92	826	1392	119	447
	92/93 1185 1798	138	475	92/93 1194 1769	118	457
	93/94 1690 2216	179	347	93/94 1754 2216	134	328
	94/95 1369 2012	173	470	94/95 1386 2012	138	488
	95/96 2797 3227	157	273	95/96 2760 3227	120	347


P Bais

T Caldes
email: maite.caldes@cnrs-imn.fr

C Guillot-Deudon

Adèle Renaud

M Boujtita

S Jobic

A Lafond

Influence of the copper deficiency and anionic composition on band-energy diagram of bulk kesterite CZTSSe

Keywords: Cu2ZnSn(S, Se)4, photovoltaics, Cu-poor, flat band potentials, XPS

Highlights  The influence of the Cu-deficiency on band-energy diagram of CZTSSe was investigated  Cu-deficiency enhances optical absorption without modification of the optical gap  Cu off-stoichiometry influences notably the positioning of EF into the band gap  The VBM shifts towards a deeper energy level J o u r n a l P r e -p r o o f 2  Cu non-stoichiometry appears as a key point to design suitable CZTSSe PV absorbers

Introduction

Although copper zinc tin chalcogenides (CZTSSe) have been extensively studied as photovoltaic absorbers, the related thin film solar cells remain with unsatisfactory performances ( < 12 %). Short minority carrier lifetime and low open circuit voltage (Voc) have been identified as the key issues for CZTSSe materials [1][2]. Formally, cationic and anionic disorders, i.e. random distributions of Cu and Zn as well as S and Se at specific J o u r n a l P r e -p r o o f crystallographic sites of the kesterite-type structure [3][4], could impact these characteristics.

However, Cu-poor CZTSSe materials appear to have better performances than their stoichiometric derivatives [5]. In that context, using advanced crystallographic techniques, a classification of point defect clusters occurring in Cu-poor Cu2-2aZn1+aSn(SxSe1-x)4 series was proposed (e.g. 2CuCu → ZnCu +VCu corresponds to the A-Type defects, 2CuCu + SnSn → 2ZnCu +ZnSn to the B ones) [6][7][8]. It was shown that A-type defect could limits cationic disorder which could explain the best optoelectronic properties of Cu-poor CZTSSe compounds [4].

On the other hand, it is well known that the efficiency of a solar cell is strongly correlated to the band alignment between the absorber and the buffer materials, and specifically to the conduction band offsets (CBO) [START_REF] Crovetto | Solar Energy Materials and Solar Cells[END_REF]. For CZTSSe compounds, the dependence of the bandgap and the band edge position versus the anionic composition (S/Se) has been studied theoretically (first-principles calculations) and experimentally (both electrochemical and spectroscopic methods). In the specific case of a CZTSSe/CdS heterojunction, the offset has been directly determined by several spectroscopic techniques [START_REF] Udaka | [END_REF]. However, its magnitude and its type (i.e. "spike-like" or "cliff-like") are still subject to debate. One possible reason to account for this still-existing controversy is chemical interdiffusion at the interface region that would strongly impact CBO [START_REF] Crovetto | Solar Energy Materials and Solar Cells[END_REF]. At that time, most studies reported in the literature on CZTS concern solely thin films or nanocrystals where an accurate determination of the chemical composition remains naturally difficult. Unfortunately, stoichiometry control is a key parameter since CZTSSe single-phased compounds exist only in a very narrow region of the pseudo-ternary Cu2(S,Se)-Zn(S,Se)-Sn(S,Se)2 phase diagram [3,[11][12][13][14]. Thus, off-stoichiometric CZTS and secondary phases are easily formed during the growth and post-growth processes, impacting possibly CBO and consequently solar cell performances [5,[15][START_REF] Pal | Solar Energy Materials and Solar Cells[END_REF]. To better understand the influence of the anionic composition of CZTSSe on band edges (with possibility to fine-tune the CBO in J o u r n a l P r e -p r o o f ultimate device approach), we prepared Cu2Zn0.5Sn0.5(SxSe1-x)4 (0  x ≤1) compounds as bulk materials. In parallel, to get insight on the property change induced by a copper substoichiometry, two Cu-poor compounds, namely Cu1.70Zn1.15SnSSe3 and Cu1.70Zn1.15SnS2Se2, were prepared and their characteristics compared to those of related stoichiometric compounds. The energy diagrams of all compounds on an absolute scale were determined combining diffuse reflectance, electrochemical impedance spectroscopy (EIS) and X-ray photoelectron spectroscopy (XPS) measurements.

Experimental

The Cu2ZnSn(SxSe1-x)4 stoichiometric compounds (0 ≤ x ≤ 1) (i.e. Cu2ZnSnSe4 (S00), Cu2ZnSnS1Se3 (S25), Cu2ZnSnS2Se2 (S50) and Cu2ZnSnS4 (S100)) and the Cu-poor compounds Cu1.70Zn1.15SnS1Se3 (A25) and Cu1.70Zn1.15SnS2Se2 (A50) were elaborated from the constitutive chemical elements (elemental Cu, Zn, Sn, S and Se) via a ceramic route as previously reported [6]. The purity of the samples was checked by powder X-ray diffraction and by scanning electron microscopy (SEM) backscattered imaging. The exact chemical composition was determined by Energy-dispersive X-ray analyses (EDX).

The analyses were performed on polished sections of the samples embedded in epoxy using calibrated standards. Table 1 sums up the target and measured compositions as well as the labelling used hereafter. The UV-visible-NIR diffuse reflectance spectra of finely ground samples were recorded with a UV-vis/NIR spectrometer (Perkin Elmer Lambda 1050) equipped with an integration sphere and controlled using the UV Winlab 6 program. The reflectance measurements were made in the 1653-496 nm range (i.e. from 0.75 to 2.5 eV) with a 2 nm step. Barium carbonate was used as white standard. The absorption (K/S) data were then calculated from the raw reflectance data using the K/S = (1-R) 2 /(2R) Kubelka-Munk formula. Spectra have been normalized for easy comparison. X-ray photoemission measurements were performed on a Kratos AXIS ultra-spectrometer using a monochromatic Al Kα X-ray source (1486.6 eV). The X-ray source ran at 150 W and the analyzed area was 700x300 m 2 . The spectrometer was calibrated using Au 4f7/2 at 83.97 eV and Cu 2p3/2 peaks at 923.63 eV measured from sputter-cleaned Au and Cu films. The Ag 3d5/2 FWHM is 0.47 eV at 20 eV pass energy giving an instrumental resolution of 0.1 eV. The sound samples were ground and pressed into 6 mm diameter pellets and annealed at 350°C in sealed silica tubes in order to avoid grain surface oxidation. The tubes were opened in a glove box and mounted on double-sided carbon tape on an aluminium plate and transferred to the XPS spectrometer using a bag under argon to limit surface oxidation. The wide scans were J o u r n a l P r e -p r o o f systematically recorded with a pass energy of 160 eV. For all compounds, high resolution spectra of C1s, O1s, Cu2p, Zn2p, Sn3d, S2p/Se2p core levels and valence bands (VB) were acquired at 20 eV pass energy. The intensities of the C1s and O1s peaks were very low and similar for all compounds, in good agreement with a little contamination of the surfaces. In addition, the C1s peak was used as energy reference (284.8eV). Spectra analysis was done using CASA XPS program [START_REF] Fairley | [END_REF]. All the spectra were fitted using Voigt functions and a linear background. A constraint on the area ratio according to their multiplicity (i.e. S2p1/2 and S2p3/2) was systematically applied.

J o u r n a l P r e -p r o o f

Electrochemical impedance spectroscopy (EIS) measurements were performed in a threeelectrode electrochemical cell using platinum as counter electrode, a saturated calomel electrode (SCE) as reference electrode and CZTSSe pellets as working electrode. A LiClO4 (1 M) aqueous solution (5.4 -7.7 pH range) was used as electrolyte. Pellets were pressed under 0.6 MPa and sintered at 600°C for 48 h into a sealed silica tube. The more friable Cu-poor compounds were consolidated by spark plasma sintering (FAST/SPS, FCT System GmbH, model HP D10-SD) as reported in [START_REF] Yang | [END_REF]. The sample powder was pressed at 410°C for 5 min under uniaxial pressure of 60 MPa under an argon atmosphere in a high-density graphite die.

Each sintered pellet was connected to a cooper wire using a silver conductive paste, then mounted in a chemical epoxy resin and polished to give a smooth electrode surface. For Cupoor samples, a KCl (0.5 M) aqueous solution (6.3 -7.2 pH range) was used as electrolyte.

Impedance spectra were obtained under AC voltage (10 mV amplitude and 100 Hz-100 kHz frequency range) in a potential ranging from -0.8 to 0 V vs. SCE. In order to determine the flat-band potential (Vfb), EIS spectra were analyzed using the Mott-Schottky method [19] at high frequencies (1-10 kHz).

J o u r n a l P r e -p r o o f

Results and Discussion

Optical band gaps.

Figure 1a shows the Kubelka-Munk transformed reflectance spectra of the four stoichiometric CZTSSe compounds. Their chemical compositions are indicated in Table 1 (see Experimental section). The estimated absorption thresholds, taken at the inflexion point, continuously increase from 1.05 to 1.16, 1.28 and 1.48 eV for corresponding targeted x values (i.e. 0, 0.25, 0.5 and 1), according an empirical y(eV) = 0.43 x' + 1.05 formula where x' is the exact S/(S+Se) ratio determined by EDX. The linear variation of the gap is directly related to that of the unit cell parameters which follow a Vegard' law. Let us notice that these measured gap values are in a good agreement with those reported in the literature [20][21][22]. 1)). It is worth noting that the optical gaps remain unchanged for both materials but the absorption onset is much steeper for A50 than S50. To J o u r n a l P r e -p r o o f check whether this is related to a scattering effect depending on particle size, both powdered samples underwent several grindings and no modification of the reflectance spectra was observed. Thus, this strong difference in the reflecting behavior at low energy is associated to the change in chemical composition of the bulk materials. and more precisely to the Cu offstoichiometry known to limit naturally the Cu/Zn disorder in kesterite [3].

Band-Edge energies and Flat band potentials

In order to determine the position of the valence band maximum (VBM), XPS measurements were undertaken. For all collected spectra, the signal-to-noise ratio is high and makes stronger our forthcoming arguments. The valence band spectra of the stoichiometric and Cu-poor CZTSSe compounds are displayed in Figure 2a. Namely, going from Cu2ZnSnSe4 to Cu2ZnSnS4, the binding energy associated to the VBM varies from 0.1 eV to 0.3 eV with respect to the Fermi level set at zero eV. For A25 and A50, J o u r n a l P r e -p r o o f the valence band is clearly shifted towards higher binding energies (VBM at 0.6 eV and 0.5 eV, respectively) supporting a deeper valence band. It is also striking that a much higher density of states at the top of the VB (built upon Cu-3d orbitals) [23] is evidenced for all stoichiometric compounds compared to non-stoichiometric ones, in agreement with a higher concentration of Cu in the targeted composition. This last observation is fully supported by comparing the analyses of surface compositions of S25 and S50 samples to those of A25 and A50 ones (see Table 2). Table 2. Chemical composition of the surfaces of S25, A25, S50 and A50 samples determined by XPS measurements. Expected values from the EDX analyses are given in parentheses in italics. To do this, the cationic core level spectra (i.e. Cu2p, Zn2p, Sn3d) were normalized with respect to those of S2p/Se3p and superimposed for each cation (see Figure 2). Let us consider A25 and S25 samples. Results indicate an under-representation of copper coupled to an overrepresentation of zinc at the surfaces of A25 compared to S25. Although both compounds exhibit the S/Se ratios quite similar to those expected, their cation surface compositions are clearly distinguishable. Indeed, the already existing Cu deficit in the A25 bulk material is significantly reinforced at the surface. The Cu/(Zn+Sn) ratio shifts from 1.3 to 0.1 from S25 to A25. Simultaneously, the Zn/(Zn+Sn) ratio increases from 0.4 to 0.6 while the Sn/(Zn+Sn)

At. % ratios Cu/Zn+Sn Zn/Zn+Sn Sn/Zn+Sn S/Se

follows the exact opposite trend (from 0.6 to 0.4). This analysis is also valid for the A50 and In order to propose energy-level diagrams of all samples on an absolute energy scale (AVS), flat-band potentials (Vfb) were determined via EIS measurements to position their Fermi level.

The method is founded on the Mott-Schottky relationship: 1/Csc 2 = 2/εε0eA 2 N (V -Vfb -kT/e) where CSC is the capacitance of the space charge region in the semiconductor, At the interfacial surface area between the semiconductor electrode and the electrolyte, ε0 the vacuum permittivity, ε the relative permittivity of the semiconductor, N majority charge carrier density, k the Boltzmann constant, T the temperature and e the electron charge. It involves thus determining the capacitance at the semiconductor/electrolyte interface (C) as a function of potential (V). In fact, at the semiconductor-electrolyte interface two capacitances in series could be considered, that of the space charge region (CSC) and that of the Helmholtz layer (CH). Thus, the total capacitance is the addition of their reciprocals (C -2 = CSC -2 + CH -2 ).

For calculations, C -2 is assumed to be equal to CSC -2 since the capacity of the Helmholtz layer is notably higher than the space charge region capacity. Mott-Schottky plots (1/Csc 2 vs. V) for all stoichiometric CZTSSe compounds are shown in Figure 3a.

J o u r n a l P r e -p r o o f The relationship between 1/Csc 2 and V is linear in depletion region and the negative slope observed for all samples confirms the expected p-type behavior (see Figure 3a) [19].

For potentials higher than 0.5V, a second slope is observed on all plots due to a surface reactivity into the electrolyte. This reaction being very pronounced for the compound S00, its flat-band potential could not be determined by EIS. As summed up in Table 3, the S/Se ratio has no influence on the Fermi level position since the value of EF is the same for all compounds. As for Fermi level, the nature of the anion has also slight influence on the VBM. However, the CBM increases with the sulfur content (from -4.4 to -3.9 eV). These results are in very good agreement with those described in the literature for thin films and nanocrystals [21][22]24]. Moreover, since the conduction band of all stoichiometric compounds is lower than that of CdS, a "spike-like" offset has to be privileged over "cliff-like" one.

Conversely, the influence of the copper deficiency on the band edge parameters is notable.

Comparing the energy diagram of the stoichiometric compounds S25 and S50 with those of their corresponding Cu-poor compounds A25 and A50 (see Figure 5), a band energy shift towards a deeper level is experimentally observed, as predicted by density functional theory (DFT) calculations [25][26][27]. This is concomitant with a slight but significant decrease of the flat band potential, around 110-190 mV for A25 and A50, respectively (see Figure 3b).

Consequently, a notable difference of the Fermi level positioning with respect of VBM is observed between the Cu-poor compounds and the stoichiometric ones.

J o u r n a l P r e -p r o o f for A25 and A50) compounds were therefore estimated from the slope of Mott-Schottky plots.

Although these values are overestimated because of the underestimated of the interfacial surface area between the porous semiconductor electrodes and the electrolyte, it appears indeed that NA of the Cu-poor compounds are two or three orders of magnitude lower than those of stoichiometric ones. These features could induce a significant change in the electrical conductivity of these compounds.

J o u r n a l P r e -p r o o f

Conclusion

The synthesis of CZTSSe compounds with a well-determined chemical composition made it possible to study precisely the influence of the anionic composition as well as that of the non-stoichiometry of copper, on the position of the bands edges and the Fermi level.

The influence of the S/Se ratio had already been described in the literature while that of the Cu deficiency had been studied mainly by DFT calculations. Our results allowed to experimentally corroborate these theoretical studies.

While the anionic composition mainly influences the optical gap and the CBM, the Cudeficiency enhances optical absorption without modification of the optical gap. Moreover, the positioning of the Fermi level with respect to the VBM is notably modified and this latter shifts towards a deeper energy level. This feature could induce a significant change in the electrical conductivity of these compounds These results confirm the key importance of controlling precisely the chemical composition of CZTSSe compounds concerning to their use as absorbers in photovoltaic devices.

Figure 1 :

 1 Figure 1: (a) Kubelka-Munk transformed reflectance spectra of stoichiometric CZTSSe

Figure 2 :

 2 Figure 2: (a) XPS spectra of stoichiometric and Cu-poor CZTSSe materials, (b) XPS

J o u r

  n a l P r e -p r o o f S50 compounds excepting for the S/Se ratio which is lower than expected, indicating an enrichment of the surface in Se.The strong Cu-depletion observed in both Cu-poor samples at their surface, might be due to the occurrence of A-type substitution defects (i.e. 2CuCu → ZnCu +VCu), already observed in these bulk materials. Furthermore, B-type substitution defects (i.e. 2CuCu +SnSn → 2ZnCu +ZnSn ) could explain the opposite evolution of the Zn/(Zn+Sn) and Sn/(Zn+Sn) ratios. This change in composition from the core to the shell should also drastically modify the optoelectronic characteristics of surfaces compared to bulk, and surely impact significantly the overall photovoltaic performances of a CZTSSe-based solar cell.

Figure 3 .

 3 Figure 3. Mott-Schottky plots for (a) stoichiometric and (b) off stoichiometric CZTSSe samples

  The flat-band potentials versus SCE (VSCE = 0.241 V vs NHE) were extrapolated to CSC -2 = 0 and then calculated versus RHE (reversible hydrogen electrode) on an absolute vacuum scale using formulae (1) et (2): J o u r n a l P r e -p r o o f Vfb versus RHE = Vfb versus SCE + 0.059 pH + VSCE (1) Efb versus AVS = -4.5 eV -Vfb versus RHE (2)

  , VBM and EF values, an energy diagram sketched in Figure4can be proposed for all compounds. Because CdS is commonly used as buffer material in CZTS solar cells, its energy diagram is also shown for comparison. For S00, in order to minimize surface reactivity, the VBM was deduced from a shorter cyclic voltammetry analysis, in which one of the anodic peaks observed could be co-related to the electron transfer mediated through valence band edge. The VBM value thus obtained ( -5.5 eV), matched well with the one reported previously in the literature.J o u r n a l P r e -p r o o f

Figure 4 :

 4 Figure 4: Schematic representation of the band structures of CZTSSe compounds on an

Figure 5 :

 5 Figure 5: Schematic representation of the band structures of the stoichiometric compounds

Table 1 :

 1 Labelling and measured chemical compositions of Cu2ZnSn(SxSe1-x)4 and Cu-

	poor samples
	sample targeted composition	EDX Composition*	x (EDX)
	S00	Cu2ZnSnSe4	Cu1.98(1)Zn1.02(2)Sn0.996(6)Se4.00(2)	0
	S25	Cu2ZnSnS1Se3	Cu1.99(2)Zn0.99(1)Sn1.019(7)S1.05(1)Se2.947(7) 0.263(2)
	A25	Cu1.70Zn1.15SnS1Se3	Cu1.81(2)Zn1.18(4)Sn1.006(9)S1.08(1)Se2.923(9) 0.268(3)
	S50	Cu2ZnSnS2Se2	Cu1.98(2)Zn1.01(2)Sn1.005(7)S2.09(2)Se1.910(7) 0.523(5)
	A50	Cu1.70Zn1.15SnS2Se2	Cu1.83(4)Zn1.15(2)Sn1.01(1)S2.11(3)Se1.99(1)	0.528 (2)
	S100	Cu2ZnSnS4	Cu1.99(2)Zn0.99(2)Sn1.02(1)S4.000(9)	1
	*Estimated standard deviations are given in parentheses

Table 3 :

 3 Characteristic of energy diagrams of studied powdered CZTSSe samples

	sample	S00 S25 A25 S50 A50 S100
	Eg (eV)	1.05 1.16 1.16 1.28 1.28 1.48
	Vfb RHE (V) -	0.57 0.47 0.51 0.32 0.54
	EF (eV)	-

Acknowledgements

The authors thank the Region Pays de la Loire and the Synchrotron SOLEIL for the cofinancing of P.B. Ph-D. thesis.

J o u r n a l P r e -p r o o f

Declaration of interests

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Author statement All authors have contributed in multiple roles.