M.-C Firpo

D Constantinescu

Study of the interplay between magnetic shear and resonances using Hamiltonian models for the magnetic field lines

Keywords: numbers: 52.55.Fa, 05.45.+a, 52.25.Gj

The issue of magnetic confinement in magnetic fusion devices is addressed within a purely magnetic approach. Using some Hamiltonian models for the magnetic field lines, the dual impact of low magnetic shear is shown in a unified way. Away from resonances, it induces a drastic enhancement of magnetic confinement that favors robust internal transport barriers (ITBs) and stochastic transport reduction. When low-shear occurs for values of the winding of the magnetic field lines close to low-order rationals, the amplitude thresholds of the resonant modes that break internal transport barriers by allowing a radial stochastic transport of the magnetic field lines may be quite low. The approach can be applied to assess the robustness versus magnetic perturbations of general (almost) integrable magnetic steady states, including non-axisymmetric ones such as the important single helicity steady states. This analysis puts a constraint on the tolerable mode amplitudes compatible with ITBs and may be proposed as a possible explanation of diverse experimental and numerical signatures of their collapses.

I. INTRODUCTION

In the present hot plasma toroidal devices, it is widely recognized that the essential limitation that prevents ignition to be reached is due to an insufficient energy confinement resulting from small-scale turbulence. During the last two decades, some decisive progress has however been made to increase the energy confinement time through the widespread experimental achievement of internal transport barriers (ITBs) [START_REF] Wolf | Internal transport barriers in tokamak plasmas[END_REF][START_REF] Connor | A review of internal transport barriers physics for steady-state operation of tokamaks[END_REF]. In these regimes, there exists a region of the tokamak plasma core where the local ion and electron thermal diffusivities are substantially reduced, nearly down to their neoclassical values [START_REF] Connor | A review of internal transport barriers physics for steady-state operation of tokamaks[END_REF]. This is likely to be a signature of the regularity of the central magnetic surfaces. In this way, the emergence of ITBs results in a sort of virtuous circle for fusion with both high central temperature T and high energy confinement time τ E increasing the triple product n e T τ E , which implies a better fusion yield. The importance of the magnetic shear, through the form of safety factor profiles, in reaching these states by quenching the core plasma turbulent transport has been firmly established. Decisive ingredients were given notably by the experimental results obtained by Eriksson et al. [START_REF] Eriksson | Discharges in the jet tokamak where the safety factor profile is identified as the critical factor for triggering internal transport barriers[END_REF] who demonstrated experimentally that the q-profile could be used as a single control parameter to obtain ITBs. In these experiments, they could obtain ITBs for non-monotonic q-profiles but not with monotonously growing q, while other parameters such as the E × B flow shear were kept the same. Additional evidence was later given by Sauter et al. [START_REF] Sauter | Inductive current density perturbations to probe electron internal transport barriers in tokamaks[END_REF] through a dedicated experiment. They controlled the q-profile using inductive current to generate positive and negative current density perturbations in the plasma center, with negligible input power. In this way, they were able to demonstrate that the electron confinement could be modified significantly solely by perturbing the current density profile. Other experiments in various tokamaks showed that having an inner region with low magnetic shear was indeed sufficient to obtain ITBs [START_REF] Connor | A review of internal transport barriers physics for steady-state operation of tokamaks[END_REF][START_REF] Razumova | Formation of electron transport barriers under ecr control of the q(r) profile in the t-10 tokamak[END_REF].

In reverse field pinches (RFPs) also, improved performance related to magnetic chaos healing [START_REF] Escande | Chaos healing by separatrix disappearance and quasisingle helicity states of the reversed field pinch[END_REF] has been recently achieved experimentally, yet through non-axisymmetric states [START_REF] Piovesan | RFX mod team, and collaborators. Self-organized helical equilibria as a new paradigm for ohmically heated fusion plasmas[END_REF][START_REF] Puiatti | Helical equilibria and magnetic structures in the reversed field pinch and analogies to the tokamak and stellarator[END_REF][START_REF] Bonfiglio | [END_REF]. It was observed that, by increasing the plasma current, the RFP operated some change of magnetic topology through a spontaneous transition towards a singlehelicity equilibrium strongly resilient to magnetic chaos.

These various experimental results emphasize the importance of the magnetic structure as they show that some configurations -reversed-shear, or low-shear, q-profiles in tokamaks or some single-helicity equilibria in RFPs -enable a substantial improvement of the confinement. This is an invitation to consider a purely magnetic approach of the problem [START_REF] Nasi | Enhanced confinement with increased extent of the low magnetic shear region in tokamak plasmas[END_REF] that will be pursued here. This amounts to address the issue of the plasma confinement in terms of the magnetic confinement, namely in terms of the confinement of the magnetic field lines. The ensuing approach is thus obviously a simplified approximation of a more complex, self-consistent, reality involving a larger number of effective degrees of freedom, including in particular the electric field. It is a sort of zero order approach that relies on the essential principle of magnetic confinement devices which is to use the magnetic field to confine the charged particles composing the fusion fuel in a plasma state [START_REF] Morrison | Magnetic field lines, Hamiltonian dynamics, and nontwist systems[END_REF]. In this sense, ensuring the confinement of the magnetic field lines may be viewed as a prerequisite to particle and energy confinement. As it is well known, this purely magnetic approach may be addressed within an Hamiltonian formalism that results from the universal divergence-free nature of B [12]. An intermediate, first-order, step towards a fully realistic description of the physics of magnetic confinement devices would be to treat the test particle motion [START_REF] Abdullaev | Spatschek Mapping of drift surfaces in toroidal systems with chaotic magnetic fields[END_REF] instead of the magnetic field lines. This is left for a future work.

In the absence of resonances, that will be specified below, the dramatic improvement on magnetic confinement induced by a low magnetic shear condition will be recalled in Section II and its link to ITBs and core stochastic transport reduction discussed. The improved confinement observed in RFPs in single-helicity states will be shown to result from the same mechanism when moving to appropriate action-angle variables. However this optimistic scenario will be tempered since, in actual tokamaks, ITBs may eventually collapse. In Section III, it will be shown, using some Hamiltonian models for the magnetic field lines, that this low shear condition happens to be strongly deleterious to magnetic confinement when approaching resonances by enabling stochasticity for much lower magnetic perturbation amplitudes. This will be further analyzed in Section IV. A discussion of this phenomenon and of its possible consequences will be finally given in Section V.

II. BENEFITS OF A LOW MAGNETIC SHEAR A. Hamiltonian formulation for the magnetic field lines

It is well known that the equations of magnetic field lines can be cast in Hamiltonian form. Actually, the field lines of a magnetic, B(x), or any other divergence-free field are the trajectories of a Hamiltonian (See e.g. the review by Boozer [12] or Abdullaev's book [START_REF] Abdullaev | Construction of Mappings for Hamiltonian Systems and Their Applications[END_REF]).

In tokamaks, a dependance on the toroidal angle breaking the rotational invariance is analogous to the usual timedependence in Hamiltonian dynamical systems. Generically one is thus left with a one-and-a-half degrees of freedom Hamiltonian, that can be decomposed into an axisymmetric, integrable part and a perturbation term that may be Fourier decomposed in the poloidal and toroidal angles. For tokamaks, the Hamiltonian reads

H(ψ, θ, φ) = H 0 (ψ) + m,n H mn (ψ) cos(mθ -nφ + χ mn ), (1)
with H 0 (ψ) = ψ W (ψ ′ )dψ ′ and W (ψ) = 1/q(ψ) the winding profile, inverse of the safety factor profile. In the tokamak configuration, H is the poloidal magnetic flux and ψ the toroidal magnetic flux conjugated to the poloidal angle θ.

B. Some results on the benefits of low shear

In real magnetic fusion devices, the magnetic perturbations given by the wave amplitudes H mn (ψ) may be induced by a variety of phenomena. They may result from external effects, such as the ripple, or may result from intrinsic phenomena such as MHD activity. In order to address the magnetic confinement problem, one may turn to simplified models for the form of the perturbed magnetic field. Restricting to the single poloidal number m = 1, taking equal phases χ 1n and some special form of the radial perturbation H 1n (ψ) identical for all n's yields the continuous version of the tokamap [START_REF] Balescu | Tokamap : A hamiltonian twist map for magnetic field lines in toroidal geometry[END_REF]. Then, considering the limit of a large number of toroidal modes makes the associated Poincaré (toka)map appear naturally through the use of periodic δ distributions.

Using the symmetric tokamap framework, it was recently shown in Ref. [START_REF] Nasi | Enhanced confinement with increased extent of the low magnetic shear region in tokamak plasmas[END_REF] that the existence of a low magnetic shear region produced a drastic enhancement of the confinement of the magnetic field lines. More specifically, the existence of some ψ s such that H ′′ 0 (ψ s ) = 0, that amounts to W ′ (ψ s ) = 0, was shown to improve significantly the robustness w.r.t. magnetic perturbations of the regular core magnetic field lines inside the ψ = ψ s surface, acting then as a transport barrier, in agreement with previous results obtained in so-called non-KAM systems including mainly maps with reversed shear profiles (See e.g. [START_REF] Del Castillo-Negrete | Area preserving nontwist maps: periodic orbits and transition to chaos[END_REF][START_REF] Balescu | Hamiltonian nontwist map for magnetic field lines with locally reversed shear in toroidal geometry[END_REF][START_REF] Constantinescu | Transport barriers and diffusion phenomena for the magnetic field lines in tokamak[END_REF][START_REF] Wingen | Stochastic transport of magnetic field lines in the symmetric tokamap[END_REF][START_REF] Portela | Tokamak magnetic field lines described by simple maps[END_REF][START_REF] Rypina | Robust transport barriers resulting from strong kolmogorov-arnold-moser stability[END_REF]). Interestingly, Ref. [START_REF] Nasi | Enhanced confinement with increased extent of the low magnetic shear region in tokamak plasmas[END_REF] gave evidence that this effect does not necessitate q to be reversed shear: the q-profile can also be monotonous with an inflexion point, which agrees with the diversity of ITBs experimental scenarios. The beneficial effect of the radial extent of the low shear zone for safety profiles with q > 1 was demonstrated. Recently, these results were more rigorously grounded: Theorems of existence of invariant circles for non-twist maps were formulated for the symmetric tokamap [START_REF] Constantinescu | Integrability versus chaos in 1 1/2 hamiltonian systems with broad perturbation spectrum. applications to the study of some transport phenomena[END_REF] using a theorem about small shear systems formulated by Ortega [START_REF] Ortega | Invariant curves of mappings with averaged small twist[END_REF].

It is important to note that the domain of applicability of this favorable regime may be extended from axisymmetric equilibria to non-axisymmetric integrable states provided ψ is replaced by some suitable action variable. Let us consider the case where some magnetic configuration involves not only an axisymmetric component H 0 (ψ) but also some (m 0 , n 0 ) component whose amplitude A is much larger than that of the other (m, n) waves of order ε. The Hamiltonian reads then

H(ψ, θ, φ) = H 0 (ψ) + AH m0n0 (ψ) cos(m 0 θ -n 0 φ), (2)
up to O(ε) terms. The dominant part of the Hamiltonian, given in Eq. ( 2), may be put into an explicitly integrable form

H(Ψ, Θ) = H 0 (m 0 Ψ) -n 0 Ψ + AH m0n0 (m 0 Ψ) cos Θ (3)
through a canonical change of variables, with the generating function F 2 (Ψ, θ, φ) = (m 0 θn 0 φ)Ψ, that amounts to moving to the (m 0 , n 0 ) wave frame. Even if H ′′ 0 (ψ) does not vanish, it is then easy to recognize that, provided A is large enough, a shearless condition in the action variable

J ≡ 1 2π ΨdΘ (4)
associated to the integrable part of the Hamiltonian (2) may however be obtained. In this way, there can be nonaxisymmetric magnetic states that are much more resilient to magnetic chaos than what would be their axisymmetric components. An illustration of this phenomenon was given in Ref. [START_REF] Escande | Chaos healing by separatrix disappearance and quasisingle helicity states of the reversed field pinch[END_REF]. It was shown that the magnetic structure involving some small m = 1 magnetic island with a separatrix was less robust towards magnetic perturbations than a reconnected structure with a higher m = 1 amplitude. This transition in terms of the winding profile in the action variable J may be interpreted as follows: In the low amplitude case where a magnetic island is present, the separatrix is associated to a vanishing of W (J) within a cusp singularity, whereas, as the amplitude grows, the reconnection process eventually makes the separatrix disappears inducing locally the growth of W (J) and smoothing the singularity, previously associated to the separatrix, which causes the appearance of a local minimum of W (J). This later favorable process may be responsible of the unexpectedly good confinement of the non-axisymmetric self-organized states very recently achieved in a RFP [START_REF] Piovesan | RFX mod team, and collaborators. Self-organized helical equilibria as a new paradigm for ohmically heated fusion plasmas[END_REF] and in a compact toroidal [START_REF] Cothran | Observation of a helical self-organized state in a compact toroidal plasma[END_REF] plasmas. For the tokamak setting, this non-axisymmetric formulation of the low shear or zero-twist condition may open new perspectives to improve magnetic confinement.

III. DETRIMENTAL EFFECT OF LOW SHEAR IN THE VICINITY OF RESONANCES USING SOME HAMILTONIAN MODELS

Let us now go back to the notations introduced in Eq. (1). The above analysis did not take into account the possibility of resonances. In particular, the fact that the winding number W s ≡ W (ψ s ), associated to the action ψ s such that H ′′ 0 (ψ s ) = 0, could be in the vicinity of some phase velocity n/m corresponding to the (m, n) magnetic perturbation was not considered. For instance, works using the tokamap framework [START_REF] Nasi | Enhanced confinement with increased extent of the low magnetic shear region in tokamak plasmas[END_REF][START_REF] Balescu | Tokamap : A hamiltonian twist map for magnetic field lines in toroidal geometry[END_REF][START_REF] Misguich | Noble internal transport barriers and radial subdiffusion of toroidal magnetic lines[END_REF] ignored de facto this possibility by restricting to the m = 1 case while taking typically q-profiles above 1 (with the noticeable exception of the sawtooth-related studies of Ref. [START_REF] Weyssow | Hamiltonian mapping of magnetic reconnection during the crash stage of the sawtooth instability[END_REF]).

In order to illustrate the interplay between low shear and resonances, let us first consider a simplification of the system (1) similar to the continuous tokamap case but where only the m = 4 poloidal modes are retained, with

H (ψ, θ, φ) = ψ dx q(x) + K H(ψ) M n=-M cos (4θ -nφ) . (5)
For M ≥ 1, this represents a one-and-a-half degrees of freedom Hamiltonian in which the perturbation amounts to a superposition of plane waves with winding numbers equal to n/4. It is convenient, and physically relevant, to consider the large M case for which the derivation of the symplectic mapping associated to ( 5) is immediate through the use of the periodic-δ distribution. For computational purposes, we used a symmetric version [START_REF] Wingen | Stochastic transport of magnetic field lines in the symmetric tokamap[END_REF][START_REF] Abdullaev | Construction of Mappings for Hamiltonian Systems and Their Applications[END_REF] of the later mapping that yields a very accurate modeling of the continuous dynamics [START_REF] Razumova | Formation of electron transport barriers under ecr control of the q(r) profile in the t-10 tokamak[END_REF] with H(ψ) = 2ψ/(1 + ψ). The physical motivation to consider this specific form of magnetic perturbation was discussed in Ref. [START_REF] Balescu | Tokamap : A hamiltonian twist map for magnetic field lines in toroidal geometry[END_REF]. Basically, this form roughly mimics the radial (or ψ) form of MHD tearing modes. It gives a vanishing perturbation in the limit ψ → 0 and ensures that the polar axis cannot be crossed so that, under the map, ψ remains positive as it should be. It would have been possible to use fully realistic MHD profiles for each (m, n) mode obtained from experiments (as e.g. in the study of Ref. [START_REF] Igochine | Zvejnieks and the ASDEX Upgrade Team[END_REF]). Our point in the present Section is however to use a simpler framework in order to emphasize the specific role of W s .

In order to use W s as a single control parameter, we considered the family of winding profiles W p (ψ) = W s -c(ψ-ψ s ) 2 with fixed c and ψ s . These correspond to reversed shear q-profiles with zero shear at ψ s . To be specific, W s was allowed to vary while putting c = .5 and ψ s = 0.4 in the whole study. In order to quantify the impact of the specific value of the minimum of the q-profile, q min , (or of its inverse W s = 1/q min ), we defined by K c the value of the stochasticity parameter K above which there exists some initial condition with ψ < 0.1 whose evolution under (5) eventually escapes the physically accessible phase space by crossing ψ = 1. Practically, for each value of W s and a given value of K, up to one hundred initial conditions were chosen at random in the range 0 < ψ < 0.1 and allowed to evolve under the symmetric map obtained from Eq. ( 5) for up to 5000 iterations. If some ψ happened to cross ψ = 1, the process was relaunched for a smaller K or else, for a larger K, in a dichotomic way, so as to converge towards the percolation-like threshold K c value. The obtained results have been plotted in Fig. 1. Let us first briefly comment on the meaning of the indicator K c . Using the form of the magnetic perturbations as H(ψ) = 2ψ/(1 + ψ) ensures that core magnetic field lines remain regular whereas border field lines (at ψ = 1) are the first to feel chaos as K increases. Having in mind fusion applications, K c has a practical meaning in terms of a measure of the confinement of the magnetic field lines, as it gives a threshold above which magnetic field lines starting from the core eventually reach the border. It is also an indicator of the global quantitative measure of the stochasticity of the magnetic field lines since, in the present model, there cannot be a connection between the core and the border of the ψ-domain through a regular magnetic field trajectory. However, with respect to the magnetic confinement issue, K c may underestimate chaos, in the sense that there may exist non-negligible subsets of the physical central phase space domain in which chaotic field lines are chaotic already for values of K K c .

It is straightforward to realize that, for the profiles W p (ψ) or W l (ψ) ≡ W s -c(ψψ s ), the dynamics is invariant under a translation in W s of some multiple of 1/4. Effectively the curve giving K c as a function of W s is 1/4-periodic in W s up to statistical discrepancies. Some part of the curve obtained for the linear strong shear case W l (ψ) has been drawn for comparison for about one period in W s . Basically, the effect of low shear happens to be efficient for W s between the n/4 resonances by increasing significantly K c w.r.t. the high shear case. The Poincaré plot of Figure 2 corresponds to such a favorable situation. This reflects the fact that, in spite of a non-negligible perturbation of magnitude K, winding profiles having a low shear region can sustain almost regular core magnetic surfaces. This magnetic chaos healing is likely to explain the stochastic transport reduction in ITBs due to the sole effect of the q-profile reported in Refs. [START_REF] Eriksson | Discharges in the jet tokamak where the safety factor profile is identified as the critical factor for triggering internal transport barriers[END_REF][START_REF] Sauter | Inductive current density perturbations to probe electron internal transport barriers in tokamaks[END_REF]. However, in the vicinity of resonances, this tendency reverses: low shear induces abrupt drops of the magnetic confinement whereas quite higher threshold values of K are reached for the high-shear linear winding profile. In this case, the statement formulated by Rosenbluth, Sagdeev, Taylor and Zaslavski [START_REF] Rosenbluth | Zaslavski Destruction of magnetic surfaces by magnetic field irregularities[END_REF] that strong shear is beneficial to preserve the regularity of the magnetic surfaces applies.

For the reversed-shear winding profile W p (ψ), low minima of K c occur for values of W s slightly above the linear values of the resonances given by n/4. Figure 3 depicts what happens in this case. The abrupt degradation of magnetic confinement in the vicinity of the resonances coincides with the emergence of some local stochastic web pattern [START_REF] Zaslavsky | Weak Chaos and Quasi-Regular Patterns[END_REF][START_REF] Zaslavsky | [END_REF][START_REF] Chernikov | Minimal chaos and stochastic webs[END_REF] that enables a large diffusion in the action ψ. Under the perturbation, separatrices of the island chains, a characteristic phase-space pattern of the non-twist Hamiltonian systems, are destroyed and replaced with stochastic channels of finite width, visible on this Figure. This produces a path to some deleterious radial transport already at low values of the amplitude of the resonant mode. The merging of separatrices is the phenomenon through which global chaos enters the system and will be more closely examined in Section IV. In another physical context, Soskin et al. studied a close phenomenon of nonlinear resonance in the minimal case of a single exciting wave [START_REF] Soskin | Drastic facilitation of the onset of global chaos[END_REF]. Here, the consideration of a large mode number M , that is physically relevant to account for turbulence in toroidal fusion devices, should provide an additional effective randomization of motion in the vicinity of the separatrix.

All other things being equal, we then considered the Hamiltonian system for the magnetic field lines (5) obtained by just replacing m = 4 by m = 2 and computed its associated symmetric map. Figure 4 displays the obtained values of the threshold K c defined for two cases of winding profiles, the parabolic one W p (ψ) and the quadratic one W q (ψ) ≡ W s -c(ψψ s ) 4 , together with the curve previously obtained for m = 4 for the reversed-shear parabolic profile. This time K c is plotted as a function of the minimum of the safety factor profile q min ≡ 1/W s , a parameter more commonly used by the tokamak community. In the m = 2 case, the positive role on magnetic confinement induced by the existence of a low shear region away from resonances appears in full clarity. The threshold K c attains there very large values and the comparison with the m = 4 case is illuminating: the span in W s between two subsequent m = 2 resonances is twice larger than the one between two subsequent m = 4 resonances and K c is then allowed to grow to larger amplitudes between resonances. The effect of the amount of flatness around q min is investigated by comparing the results obtained for W p to that for the quadratic flattest profile W q in the m = 2 case. Away from resonances, it is clear that expanding the low shear zone improves the resilience of the magnetic structure towards perturbations as shown in Ref. [START_REF] Nasi | Enhanced confinement with increased extent of the low magnetic shear region in tokamak plasmas[END_REF]. However, in the vicinity of resonances, the effect reverses. A lower minima of K c is obtained for W q and it is reached closer to the resonance q = 2/n (with n = 1 in Fig. 4).

IV. LOCAL ANALYSIS

We wish to interpret the above results within the continuous Hamiltonian system. Let us consider then the Hamiltonian

H (ψ, θ, φ) = ψ W (x)dx + K H(ψ) M n=-M cos (m 0 θ -nφ) (6)
for some given poloidal mode m 0 and take W (ψ) = W s -c(ψψ s ) k , with c > 0 and W s = W (ψ s ). Let us proceed to a canonical change of variable with the generating function

F 2 (ψ ′ , θ, φ) = (m 0 θ -n 0 φ) ψ ′ + ψ s θ. (7)
This yields

∂F 2 ∂θ = ψ = m 0 ψ ′ + ψ s , (8)
∂F 2 ∂ψ ′ = θ ′ = m 0 θ -n 0 φ, (9)
and the new Hamiltonian reads, up to a constant,

H (ψ ′ , θ ′ , φ) = (m 0 W s -n 0 ) ψ ′ - c k + 1 (m 0 ψ ′ ) k+1 + K H (m 0 ψ ′ + ψ s ) M n=-M cos (θ ′ + n 0 φ -nφ) . (10)
Sufficiently close to the resonance n = n 0 , θ ′ + n 0 φnφ is slowly varying, so that one can average over the fast "time" φ. This yields the one degree of freedom Hamiltonian

H (ψ ′ , θ ′ , φ) = (m 0 W s -n 0 ) ψ ′ - c k + 1 (m 0 ψ ′ ) k+1 + K H (m 0 ψ ′ + ψ s ) cos θ ′ . (11)
The study of its topology as a function of W s and K should help to clarify the phase space location where the stochasticity, induced by the perturbation to [START_REF] Morrison | Magnetic field lines, Hamiltonian dynamics, and nontwist systems[END_REF], should first emerge as K increases from zero. Figure 5 shows the transition from an homoclinic topology with dimerized islands (plot e) to an heteroclinic topology with twin islands (plots g and h) as W s increases while being slightly above n 0 /m 0 = 1/2. The bifurcation between these two topologies [START_REF] Van Der Weele | The birth process of periodic orbits in non-twist maps[END_REF] occurs at some reconnection threshold W * s (see e.g. Ref. [START_REF] Howard | Stochasticity and reconnection in Hamiltonian systems[END_REF][START_REF] Del-Castillo-Negrete | Coherent structures and self-consistent transport in a mean field Hamiltonian model[END_REF]) that depends on the perturbation amplitude K and can be easily computed by equaling the energy values at the x-points obtained for θ = 0 and θ = π.

As it is well known, the nonintegrable perturbation coming from the other modes will gradually replace, as K increases, separatrices by stochastic layers. On the basis of Figure 5, one can argue that the largest extent of the stochastic domain shall take place in the case where W s is slightly above n 0 /m 0 , close to the situation corresponding to the plot e, in which case the stochastic layers coming from the two separatrices shall merge for a sufficiently large value of the perturbation amplitude. In this way, this analysis gives account of the fact that the smallest values of the critical thresholds take place, for the given ψ-form of the perturbation, not for W s = n 0 /m 0 but slightly larger values. The Reader is referred to [START_REF] Abdullaev | Construction of Mappings for Hamiltonian Systems and Their Applications[END_REF][START_REF] Morrison | Nontwist maps[END_REF] for extra reference studies on the transition to chaos in Hamiltonian systems with non-monotonic frequency.

V. DISCUSSION

The previous qualitative, yet systematic, study has emphasized the dual impact of low magnetic shear on magnetic confinement in toroidal devices. The main achievement of this approach has been to hopefully clarify the dual impact of low magnetic shear through a unified picture. In particular, it has been shown using some Hamiltonian models for the magnetic field lines that, away from any consideration on dynamics and stability, there exists potentially very low thresholds on magnetic perturbations above which safety factor profiles having a low shear region occurring for q-values close to q = m/n induce a loss of magnetic confinement. This may now be used both to predict and to interpret experimental observations.

For instance, assume that some reversed-shear q-profile associated to an ITB regime evolves in such a way that q min decreases and crosses e.g. q = 2. Then, if there exists some magnetohydrodynamic (MHD) mode with m = 2n, typically a (2,1) mode, having a sufficiently large amplitude, one predicts a collapse of this ITB due to an abrupt macroscopic degradation of magnetic confinement. This prediction does not need any consideration on the nature of the instability that produced this particular (2,1) mode. However, to make a quantitative statement on the amplitude threshold inducing the collapse one should have some knowledge on the radial form of the modes involved . This prediction does not either preclude the fact that ITBs may be more easy to trigger experimentally for rational values of q min for low MHD activity, a fact on which there does not seem to be a general agreement [START_REF] Sauter | Inductive current density perturbations to probe electron internal transport barriers in tokamaks[END_REF] anyway. There is some subtlety here: For instance, for the original symmetric tokamap, that corresponds to the case where only m = 1 (non-resonant) perturbations are retained, the most robust magnetic surfaces would actually be obtained for q min close to 2. It could however explain the fact that the good confinement phase obtained in these cases is always transient [START_REF] Tuccillo | Development on jet advanced tokamak operations for iter[END_REF] as various (m, n) modes should eventually be destabilized either linearly or nonlinearly through mode couplings. Moreover, as it is well known from the tokamak linear MHD theory, the low m and n MHD modes are typically the most dangerous ones, in the sense that they are associated to the largest linear growth rates. Consequently, one expects drops of magnetic confinement to take place in the first place for low-order rational values of q min , as the amplitude thresholds requested for the large-scale degradation of the magnetic confinement should be first attained for those low (m, n) modes.

As illustrated in Fig. 4, the negative impact of low shear should be either benign and transient in the case where the q-profile happens to rapidly cross low rational values of q min , or particularly deleterious, if the q-profile happens to be clamped and very flat about q min close to some q ≡ m/n, with a large enough amplitude in the (qn, n) modes. These ITBs collapses triggered by stochastization of the magnetic field lines should manifest through flattenings of the pressure and temperature profiles in the stochastic region.

Interestingly enough, these predictions agree with numerous experimental and numerical observations on ITBs dynamics. We should first note that experimental evidences indicate that the so-called Edge-Localized modes (ELMs) may not be responsible for ITBs collapses since these are observed to happen already before the ELMs occur [38]. Additionally, ITB do collapse while the q-profile is reversed-shear and there is a variety of experimental report of core confinement collapses in reversed-shear cases when q min reaches values close to low order rational values such as 3/2, 2 or 3 [38][START_REF] Takechi | 60U Team. Mhd instabilities leading to disruptions in low beta jt-60u reversed shear plasmas[END_REF][START_REF] Hogeweij | and contributors to the EFDA-JET workprogramme. Electron heated internal transport barriers in jet[END_REF]. The extensive experimental and numerical results of Ref. [START_REF] Maget | Nonlinear magnetohydrodynamic simulation of tore supra hollow current profile discharges[END_REF] consistently show that pressure crashes associated to core magnetic stochasticity are strongly correlated with the closeness of q min to 2 and with a sufficient amplitude in the (2, 1) mode. The proposed picture is also consistent with recent numerical results showing that an ITB collapse occurs for q min in the vicinity of 4/3 when turbulence involves a sufficiently large amplitude of the (4, 3) mode [START_REF] Tokunaga | Multi-scale transport simulation of toroidal momentum source profile effect on internal transport barrier collapse[END_REF]. Moreover, the above predictions should be most easily evidenced on electron ITBs, since due to their small Larmor radius, electrons provide a scan of the magnetic field structure. In this respect, it would be interesting to investigate further the connections with the experimental results of Ref. [START_REF] De Baar | Tokamak plasmas with dominant electron cyclotron heating; evidence for electron thermal transport barriers[END_REF] in which the disappearance of electron transport barriers are associated with q min crossing simple rational values.

Finally, resonances with m = n associated to the q = 1 surface are by no way special. Sawteeth crashes could then be interpreted as ITB collapses associated to q = 1. A complete modeling would however require an enlarged framework that is beyond the scope of the present article.

FIG. 1 :

 1 FIG.1:The critical value of the stochasticity parameter K, as defined in the text, as a function of the inverse of the minimum of the q-profile, Ws, for the parabolic family Wp(ψ) (dashed line) and as a function of the parameter Ws in the linear case W l (ψ) (plain line). The vertical dashed lines mark the locations of the resonances Ws = 1/4, 2/4 and 3/4.

FIG. 2 :

 2 FIG.2: Poincaré plot of the system (5) with the parabolic reversed-shear profile Wp(ψ) for K = 0.03 and Ws = 0.645.

FIG. 3 :FIG. 4 :

 34 FIG.3: Poincaré plot of the system (5) with the parabolic reversed-shear profile Wp(ψ) for K = 0.03 and Ws = 0.52.

1 dFIG. 5 :

 15 FIG. 5: Phase space of the Hamiltonian (11) for some values of Ws close to n0/m0 with m0 = 4 and n0 = 2, for c = 0.5 and K = 0.01.

Acknowledgments

Comments by A. Vasiliev on stochastic webs are gratefully acknowledged. MCF thanks L. Nasi and J.-M. Rax for useful discussions. This work was carried out within the framework the European Fusion Development Agreement and the French Research Federation for Fusion Studies. It is supported by the European Communities under the contract of Association between Euratom and CEA. The views and opinions expressed herein do not necessarily reflect those of the European Commission.


Sylvain Biancamaria
email: sylvain.biancamaria@legos.obs-mip.fr

Anny Cazenave

Nelly M Mognard

Frédéric William Llovel

Frédéric Frappart

Satellite-based high latitude snow volume trend, variability and contribution to sea level over 1989/2006

Keywords: SSM/I, snow volume, high latitude, climate indices, sea level

come

Introduction

High latitude regions are the most affected by current climate change (e. g. [START_REF] Trenberth | Observations: Surface and Atmospheric Climate Change[END_REF]. Different components of the Arctic hydrological cycle have experienced important modifications since the beginning of the 20 th century. For example river discharge has significantly increased [START_REF] Stocker | Water cycle shifts gear[END_REF], snow extent is decreasing [START_REF] Déry | Recent Northern Hemisphere snow cover extent trends and implications for the snow-albedo feedback[END_REF][START_REF] Brown | A multi-data set analysis of variability and change in Arctic spring snow cover extent, 1967-2008[END_REF] and the annual duration of the period with unfrozen soil conditions has increased [START_REF] Groisman | State of the ground: climatology and changes during the past 69 years over Northern Eurasia for a rarely used measure of snow cover and frozen land[END_REF]. Snow volume is a key variable to understand the evolution of the high latitude hydrological cycle. High latitude rivers discharge is mainly driven by the accumulated snow volume and the timing of its melting, leading to extremely important floods in spring [START_REF] Yang | Streamflow response to seasonal snow cover extent changes in large Siberian watersheds[END_REF]. So, snow cover extent and depth are among the Essential Climate Variables (ECV) of the Global Climate Observing System (GCOS; [START_REF] Sessa | Terrestrial Essential Climate Variables[END_REF] and observations of their temporal evolution are of critical importance.

Few in situ snow observations are available at high latitude, thus providing limited information on global and regional snow depth fields [START_REF] Brown | Northern Hemisphere snow cover variability and change, 1915-97[END_REF]. Remote sensing techniques complement the in situ data, but most analyses focused on snow extent change (see [START_REF] Trenberth | Observations: Surface and Atmospheric Climate Change[END_REF] for a review), which only partially characterize snow variability. So far, interannual to multidecadal changes in snow volume have been mainly estimated using hydrological model outputs (e.g. [START_REF] Milly | Contribution of climate-driven change in continental water storage to recent sea-level rise[END_REF]. In this study, we use satellite-based microwave observations to derive and analyze high latitude snow volume changes over 1989/2006. Correlations between snow volume and climate indices have been investigated and the Arctic snow contribution to the global mean sea level variation has been estimated.

Data analysis

Snow volume used in this study has been computed from Special Sensor Microwave/Imager (SSM/I) data. SSM/I measures the Earth brightness temperature for different microwave frequencies in both horizontal and vertical polarizations (19.35 GHz,37 GHz,85.5 GHz and 22.235 GHz). Since July 1987, this instrument has been operating on board the operational Defense Meteorological Satellite Program satellite series. Daily SSM/I data, mapped to the Equal Area SSM/I Earth Grid projection with a 25×25 km 2 resolution, are provided by the National Snow and Ice Data Center (Armstrong et al., 1994). A dynamic algorithm that takes into account the temporal evolution of the snow grain size, developed by [START_REF] Mognard | Northern Great Plains 1996/97 seasonal evolution of snowpack parameters from satellite passive-microwave measurements[END_REF], has been used to retrieve snow depth from SSM/I measurements. This algorithm has been first developed for the U.S. Northern Great Plains, then amended by [START_REF] Grippa | Siberia snow depth multiyear average derived from SSM/I data using a combined dynamic and static algorithm[END_REF] and applied over Western Siberia. Recently, a multiyear (1988/1995) averaged snow depth field computed using this algorithm has been validated over Siberia [START_REF] Boone | The impact of simulated soil temperatures on the estimation of snow depth over Siberia from SSM/I compared to a multi-model multi-year average[END_REF] and the whole high latitude regions, Greenland excluded [START_REF] Biancamaria | A satellite snow depth multi-year average derived from SSM/I for the high latitude regions[END_REF]. Interannual variability of these data has been validated over the Ob river basin, in Western Siberia, by comparison with discharge measurements at the estuary [START_REF] Grippa | Comparison between the interannual variability of snow parameters derived from SSM/I and the Ob river discharge[END_REF]. The latter study found a significant correlation between the snowmelt date and the discharge in May (correlation coefficient = -0.92) and between the winter snow depth and the discharge in June (correlation coefficient = 0.61).

The inputs for this algorithm are the difference between 19.35 GHz and 37 GHz brightness temperature in horizontal polarization from SSMI, the air/snow interface temperatures from the National Centers for Environmental Prediction global (NCEP) reanalysis [START_REF] Kalnay | The NCEP/NCAR 40-year reanalysis project[END_REF] and the snow/ground interface temperatures modeled by the "Interaction between Soil-Biosphere-Atmosphere" (ISBA) scheme forced by the Global Soil Wetness Project-Phase2 P3 precipitation field [START_REF] Boone | The impact of simulated soil temperatures on the estimation of snow depth over Siberia from SSM/I compared to a multi-model multi-year average[END_REF]. For the present study, we consider all continental surfaces above 50°N, Greenland excluded, composed of two sub-regions: Eurasia (0°E<longitude<191°E) and North America (191°E<longitude<360°E).

Regions below 50°N have not been taken into account as the snowpack is highly variable spatially and of low amplitude, hence difficult to observe with a 25x25 km 2 spatial resolution.

According to the snow climatology over North America from [START_REF] Brown | Gridded North American monthly snow depth and snow water equivalent for GCM evaluation[END_REF], between December and March, 80% of the total North American snow volume is found above 50°N.

This study focuses on monthly and yearly-averaged total snow volume (sum of all non-zero snow depth pixels multiplied by the pixel area). Yearly averages are centered on winter months (i.e. the yearly average for year n corresponds to the temporal average from October year n-1 to September year n).

Snow volume derived from SSM/I measurements for January, and temporally averaged over 1989/2006, has been compared to the 1979/1996 climatology for North America from [START_REF] Brown | Gridded North American monthly snow depth and snow water equivalent for GCM evaluation[END_REF] and to the global climatology from U.S. Air Force/Environmental Technical Applications Center (USAF/ETAC) [START_REF] Foster | Global snow depth multi-year average[END_REF].

The correlation coefficient between snow depth from SSM/I and [START_REF] Brown | Gridded North American monthly snow depth and snow water equivalent for GCM evaluation[END_REF] is 0.36.

However most of the differences between the two datasets are found over regions covered with tundra, according to the snow classification from [START_REF] Liston | Global Seasonal Snow Classification System[END_REF]. For tundracovered regions (42% of North America), the correlation coefficient is equal to 0.20, whereas over the remaining of North America the correlation coefficient is equal to 0.60 (the tundra covers 42% of North America). This result is consistent with the previous study by [START_REF] Biancamaria | A satellite snow depth multi-year average derived from SSM/I for the high latitude regions[END_REF], which found that the dynamic algorithm does not perform well over the Northern part of the continent (where the tundra is located) due to the presence of numerous lakes. The SSM/I data over these regions need to be processed using a specific algorithm, such as the one developed by [START_REF] Derksen | Development of a tundra-specific snow water equivalent retrieval algorithm for satellite passive microwave data[END_REF], to take into account the specificity of snow emissivity over lakes. Yet, this still remains an open issue, which requires further investigations since, the in situ snow gauges, used by [START_REF] Brown | Gridded North American monthly snow depth and snow water equivalent for GCM evaluation[END_REF] to compute a climatology by interpolation are very scarce above 50°N, especially for the Northern part of the continent. snow depth from SSM/I correlates better with the USAF/ETAC climatology (correlation coefficient equals to 0.62 over the whole study domain, 0.67 over Eurasia and 0.52 over North America ), These correlation coefficients are highly significant, as their pvalues (i.e. the probability to obtain these coefficients by random chance, whereas the variables are uncorrelated) are extremely small (lower than 0.001). The snow depth retrieved from SSM/I measurements are not directly compared to in situ measurements as the spatial resolution of SSM/I is too coarse to be compared to local measurements, and high latitude networks of in situ snow depth measurements are not dense enough to allow an estimation of the mean snow depth over a 25x25 km 2 area. Based on a statistical analysis over the U.S. Northern Great Plains, [START_REF] Chang | Analysis of ground-measured and passive microwave derived snow depth variations in mid-winter across the Northern Great Plains[END_REF] showed that error between one single in situ measurement and the mean snow depth over a 1°x1° region can be up to 20 cm (for a range of snow depth values between 1.5 cm and 45.4 cm).

The temporal variability of the SSM/I-derived snow volume is analyzed in the next section. It agrees well with previous published studies, giving high confidence in the quality of the snow volume time series estimated from SSM/I observations.

Snow volume temporal variability

Figure 1 shows the monthly anomalies of snow volume time series averaged over the study area from SSM/I, from an inversion of the Level-2 GFZ Gravity Recovery And Climate Experiment (GRACE) products [START_REF] Ramillien | Time variations of the land water storage from an inversion of 2 years of GRACE geoids[END_REF][START_REF] Frappart | Evolution of high-latitude snow mass derived from the GRACE gravimetry mission (2002-2004)[END_REF], from the European Centre for Medium-Range Weather Forecasts ReAnalysis (ERA)-interim product [START_REF] Uppala | Towards a climate data assimilation system: status update of ERA-Interim[END_REF], and outputs from two Land Surface Models (LSM) used by the Global Land Data Assimilation System (GLDAS) [START_REF] Rodell | The global land data assimilation system[END_REF]: MOSAIC [START_REF] Koster | Energy and water balance calculations in the Mosaic LSM[END_REF] and NOAH 2.7 [START_REF] Chen | Modeling of land-surface evaporation by four schemes and comparison with FIFE observations[END_REF], between January 2003 and June 2006. The mass variations measured by GRACE have been converted into snow volume assuming a constant snow volume density of 300 kg.m -3 . The time variations of the snow volume are dominated by the seasonal cycle. NOAH outputs present a better agreement with GRACE data both in terms of amplitude and timing. Snow volume from SSM/I is in better agreement with MOSAIC for the amplitude and with ERA-interim for the phase. However, all datasets agree well in phase and their amplitudes have the same order of magnitude, except for ERAinterim which seems to overestimate the amplitude. The GRACE land water and snow solutions used in this study, are based on the development of geopotential harmonic coefficients up to a degree 50, which correspond to a spatial resolution of 400 km (see Frappart et al. (in press) for more details about this dataset). This leads to smaller amplitudes and smoothes the temporal time series. Table 1 presents the annual snow volume trends over 2003/2006 computed from SSM/I-based snow depth, from the snow reservoir extracted from GRACE measurements and from the total GRACE signal over the whole study domain, Eurasia and North America. Two trends computed from GRACE data are shown: the first one has been directly computed from the GRACE (snow and total) time series and the second one has been corrected from the Post-Glacial Rebound (PGR) trend estimated by [START_REF] Paulson | Inference of mantle viscosity from GRACE and relative sea level data[END_REF], available on the GRACE Tellus website (http://grace.jpl.nasa.gov). The uncertainties on PGR trends are supposed to be around 20% [START_REF] Paulson | Inference of mantle viscosity from GRACE and relative sea level data[END_REF] and are given in mm.year -1 of equivalent water thickness, which have been converted into snow volume trends using a constant snow density of 300 kg.m -3 . The PGR trend (in equivalent snow volume per year) is equal to 495.6 km 3 .year -1 , 31.2 km 3 .year -1 and 464.4 km 3 .year -1 over the whole study domain, Eurasia and North America, respectively. SSM/I and PGR corrected GRACE trends are all negative, as previously observed at basin-scale for the GRACE data by Frappart et al. (in press), yet they are statistically significant only over North America. The differences between SSM/I and GRACE trends are likely caused by the sensitivity due to the small number of years available for the computation, the truncation of the GRACE data (which caused a loss of energy in the short spatial wavelengths) and PGR uncertainty.

Interannual total snow volume time series (seasonal signal removed) from SSM/I has been computed over 1989/2006 for the whole study domain (Figure 2a) and separately for Eurasia (Figure 2b) and for North America (Figure 2c). Table 2 presents ERA-interim data has not been used because of biases in the background forecast and in the assimilated observations that makes them unreliable for trend estimation [START_REF] Trenberth | Observations: Surface and Atmospheric Climate Change[END_REF]. Even if an effort has been undertaken to reduce the biases in ERA-interim [START_REF] Dee | Variational bias correction of satellite radiance data in the ERAinterim reanalysis[END_REF], trends computed with this dataset should still be used with caution. Outputs from GLDAS are not shown, due to the presence of an obvious bias in these datasets between the 1989/1999 and 2000/2006 time periods. For NOAH, the mean snow volume and standard deviation over 1989/1999 are equal to 6239 km 3 and 329 km 3 , respectively, whereas over 2000/2006 they are equal to 4760 km 3 and 183 km 3 , respectively. Thus, trends computed from these snow depth fields will mainly be the result of this bias. SSM/I snow volume over Eurasia displays a positive, but not statistically significant, trend of 11.3±9.3 km 3 .year -1 (p-value=0.25, Figure 2b), while over North America the trend is negative and statistically significant (-9.7±3.8 km 3 .year -1 , p-value=0.02, Figure 2c). Trends have been computed using the generalized linear model regression [START_REF] Dobson | An introduction to generalized linear models[END_REF], and the uncertainty given with each trend corresponds to the standard error on the estimation of the slope from the linear regression algorithm used. These uncertainties are high as snow volume time series have large variability. When the yearly snow volume is averaged over the whole study domain, the trend is positive, remains small, with a very large uncertainty and is not statistically significant (1.5±10.5 km 3 .year -1 , p-value=0.88, Figure 2a). Figure 3 2006) observed a significant increase in the number of days with unfrozen soil conditions between 1956 and 2004. Yet, this increase is most frequently due to a reduction of days with frost and ice on the ground rather than a snow cover retreat. Besides, these modifications tend to diminish during the last decade of the twentieth century. These observations agree with our results (a positive, but not significant, snow volume trend over Eurasia). [START_REF] Bulygina | Changes in snow cover over Northern Eurasia in the last few decades[END_REF] used 820 in situ stations over Russia between 1966 and 2007 to infer snow depth trend maps. They found that snow cover periods tend to be shorter, however the amount of snow fall tends to increase, leading to positive snow depth trends over Eurasia: maximum storage change increased from 0.2 cm.year -1 to 0.6/0.8 cm.year -1 (with maximum rates in Western Siberia). These trends have the same order of magnitude than trends shown in Figure 3. Based on in situ measurements, [START_REF] Kitaev | Distribution of snow cover over Northern Eurasia[END_REF] found a positive snow depth trend (0.09 cm.year - 1 ) over Eurasia (for latitudes above 40°N) for February between 1936/2000. This study also showed opposite trends between snow water equivalent in February for 1966/2000 over Eurasia and North America (the trends are equal to 0.743 mm/decade and -1.231 743 mm/decade, respectively). For February of the 1989/2006 time span, SSM/I-based snow depth trends are equal to 0.14 cm.year -1 (p-value=0.15) over Eurasia and -0.18 cm.year -1 (p-value=0.02) over North America. These results are in agreement with the results found by [START_REF] Kitaev | Distribution of snow cover over Northern Eurasia[END_REF].

Using different data sets (visible and microwave satellite observations, objective analyses of surface snow depth observations, reconstructed snow cover from daily temperature and precipitation, and proxy information derived from thaw dates), [START_REF] Brown | A multi-data set analysis of variability and change in Arctic spring snow cover extent, 1967-2008[END_REF] showed that snow cover extent in June and May respectively decreased by 46% and a 14% in the Arctic (latitude>60°N), during the 1967/2008 time period. This reduction is observed over both Eurasia and North America and 56% of snow cover extent variability for June and 49% for May is explained by air temperature. Using a less accurate data set, [START_REF] Déry | Recent Northern Hemisphere snow cover extent trends and implications for the snow-albedo feedback[END_REF] showed that snow cover extent also decreased in Eurasia and North America during winter time for the 1972/2006 time span. However, the observed decline is smaller than during spring. The snow volume variability computed in our study does not seem to be consistent with the trend in snow cover extent observed in these previous studies, especially over Eurasia. Yet, this was expected. In effect, [START_REF] Ge | Observed inconsistencies between snow extent and snow depth variability at regional/continental scales[END_REF] showed that, at continental/regional scales, high latitude snow extent and snow depth are largely unrelated.

The decreasing trend of snow volume over North America estimated from SSM/I is in agreement with the previous study from [START_REF] Dyer | Spatial variability and trends in snow depth over North America[END_REF]. Using interpolated in situ measurements during the 1960/2000 time span, they found that snow depth has a decreasing trend in January/February over the time period, which becomes even steeper around March, along with an earlier onset of spring thaw, which could explain the decreasing North American snow volume trend observed in our study.

It is extremely difficult to estimate the implications of the observed decreased in North

American snow volume on other snow related parameters, like glaciers mass, and an answer to this issue is far beyond the scope of this paper. For example, very recently, [START_REF] Berthier | Contribution of Alaskan glaciers to sea-level rise derived from satellite imagery[END_REF] confirmed that Alaskan glaciers are losing mass. However, glaciers mass loss is not only observed in North America, but is widely measured on all continents [START_REF] Kaser | Mass balance of glaciers and ice caps: Consensus estimates for 1961-2004[END_REF].

Therefore, it is hard to assess if there is any relation between the Alaskan glaciers mass loss and the decreasing snow volume in North America, and this issue will require further investigations.

Relationship between snow volume and climate indices

To investigate the causes of snow volume variability in North America and Eurasia, yearly mean SSM/I snow volume anomaly has been correlated to climate indices representing the dominant modes of atmospheric and ocean variability. The following climate indices have been considered:

-Arctic Oscillation (AO, leading mode from the Empirical Orthogonal Function analysis of monthly mean height anomalies at 1000-hPa, poleward of 20°N). A positive (negative) AO index corresponds to a lower (higher) than normal atmospheric pressure over the Artcic, which leads to stronger (weaker) westerly winds. Therefore, in positive AO phase , cold Arctic air is maintained in the Northern part of America (Arctic coast and Quebec), while the rest of America, Europe and Asia experiences a warmer then averaged winter, with more precipitation in Northern Europe. On the contrary, in negative AO phase, cold Arctic air reaches lower latitude (South Canada, US, Asia and Europe), whereas the Northern parts of America is warmer than during the positive AO phase [START_REF] Thompson | The Arctic Oscillation signature in the wintertime geopotential height and temperature fields[END_REF].

-Atlantic Multidecadal Oscillation (AMO, North Atlantic mean sea surface temperature anomaly north of the equator). AMO corresponds to cycles of warming and cooling of the North Atlantic Ocean with a period comprised between 50 and 80 years. This cycle affects the North Atlantic branch of the thermohaline circulation and therefore the whole oceanic system [START_REF] Kerr | A North Atlantic climate pacemaker for the centuries[END_REF]. A positive (negative) phase of the AMO leads to more (less) summer precipitations in Northern Europe and Alaska and less (more) summer precipitations in the U.S. and South Canada [START_REF] Enfield | The Atlantic Multidecadal Oscillation and its relationship to rainfall and river flows in the continental U.S[END_REF]. [START_REF] Knight | Climate impacts of the Atlantic Multidecadal Oscillation[END_REF] show, using a climate model, that positive AMO phase tends to strengthen broad cyclonic pressure anomalies over the Atlantic and Europe in winter, therefore increasing precipitations on these regions. During the time span of the study (1989/2006), the AMO has shifted from a negative to a positive phase around 1995.

-Pacific Decadal Oscillation (PDO, leading principal component of monthly sea surface temperatures in the North Pacific, poleward of 20°N). PDO is an El Niño-like pattern characteristic the North Pacific climate variability with interannual to interdecadal fluctuations. PDO influences mainly North America climate during winter time and is positively correlated with precipitation along the coasts and central Gulf of Alaska and negatively correlated over much of the interior of North America [START_REF] Mantua | A Pacific interdecadal climate oscillation with impacts on salmon production[END_REF].

-Pacific North American pattern (PNA, second component of the Northern Hemisphere extra tropical sea level pressure anomalies). During positive (negative) phase of the PNA, geopotential height anomalies are positive (negative) along the West coast of North America and negative (positive) in the mid-Pacific and Eastern US. Therefore, negative PNA phase is characterized by a strong East Asian jet stream, which is blocked during positive PNA phase.

The spatial scale of the PNA pattern is at its most extent during winter [START_REF] Wallace | Teleconnections in the geopotential height field during the Northern Hemisphere winter[END_REF]. 3 gives the correlation coefficients between the climate indices and the annual snow volumes over North America and Eurasia. The AO index is relatively well correlated with snow volume over North America (correlation=0.51 and p-value=0.03, Figure 4b) and anti-correlated with snow volume over Eurasia (correlation=-0.57 and p-value=0.01, Figure 4a), thus the climatic conditions represented by the AO index (which is the dominant mode of interannual variability in North Hemisphere) play a significant and opposite role over the two continents.

It is worth mentioning that SSM/I snow volume over Eurasia and North America are not correlated (correlation=-0.07 and p-value=0.80). AMO and PNA are anti-correlated with snow volume over North America (correlation=-0.59 with p-value=0.01 and correlation=-0.66 with p-value=0.003, respectively) and not, or only weakly, correlated with snow volume over respectively). On the contrary, PDO is more strongly linked with snow volume over Eurasia (correlation=0.49 with p-value=0.04) than over North America (correlation=-0.18 with p-value=0.47). Figure 5 presents snow depth linear regression maps based on each climate index (i.e. the correlation coefficient between snow depth and the climate index multiplied by the snow depth standard deviation for each pixel), which clearly show the locations where snow depth co-varies with each climate index. On these maps, regression coefficients are presented only for pixels which have a statistically significant (p-value<0.1) correlation coefficient. Regression map between snow depth and AO index (Figure 5a) clearly shows the opposite impact of AO index over snow depth between North American Arctic coast and the Eastern part of Siberia (with positive regression coefficients) and middle Siberia and Eastern Europe (with negative regression coefficients). There are few statistically significant regression coefficients between snow depth and AMO index over Eurasia (Figure 5b), some negative regression coefficients over Quebec and North American Arctic coast and some positive coefficients over the Rockies. Snow depth weakly co-varies with PDO index over North America (Figure 5c), contrarily to Southern middle Siberia and Northeastern Europe, where regression coefficients are positive. Finally, regression coefficients between snow depth and PNA index (Figure 5d) show that PNA is linked with snow depth especially in central North America (negative coefficients), in the Rockies (positive coefficients) and in Southern middle Siberia (positive coefficients). From these results, it seems that AO is the only mode affecting significantly snow volumes over both Eurasia and North America and could explain the different behavior of snow volume over these two continents.

Previous studies [START_REF] Cohen | Stratopshere-Troposphere coupling and links with Eurasian land surface variability[END_REF][START_REF] Orsolini | Role of Eurasian snow cover in wintertime circulation: Decadal simulations forced with satellite observations[END_REF] have shown that high snow cover in late autumn over Eurasia can create upward propagating planetary wave pulses in winter, which weaken the polar vortex and therefore lead to negative AO in late winter, explaining a negative correlation between snow cover extent and AO index. The same process also explains the negative correlation between snow volume over Eurasia and winter AO index observed in our study. [START_REF] Ge | North American snow depth and climate teleconnection patterns[END_REF] compared monthly AO, NAO (North Atlantic Oscillation, commonly seen as a regional manifestation of AO), PDO and PNA indices with monthly in situ interpolated snow depth over North America during the 1956/2000 time span. They found that snow depth has weak correlation with monthly AO and NAO indices, but strong anticorrelation with monthly PDO and PNA indices for the winter months (December to April). However, for PDO, only correlation coefficients for February and March are significant at 90% confidence level (for PNA all the correlations during all winter months are above the 90% confidence level). Our results seem to differ partially, as we found that snow volume significantly correlates with AO, significantly anti-correlates with PNA and has no significant linear relationship with PDO over North America. However, it should be noted that Ge and Gong ( 2009) studied a wider domain (with latitudes below 35°N) and found the highest relations between snow depth and PDO/PNA over interior central-western North America, which expands far below 50°N (the Southern limit of our study domain). [START_REF] Ghatak | North American temperature, snowfall, and snow-depth response to winter climate modes[END_REF] used the same North American in situ-based snow depth field than [START_REF] Ge | North American snow depth and climate teleconnection patterns[END_REF] and found that globally snow depths correlate negatively with both winter NAO and PNA. Yet, their study domain includes the whole North American continent. If only the latitudes above 50°N are considered (Figure 4 and 6 in [START_REF] Ghatak | North American temperature, snowfall, and snow-depth response to winter climate modes[END_REF], they showed that snow depth correlates positively with NAO and negatively with PNA, as found in our study (Table 3).

Their explanation of these correlations is the following: positive phase of winter NAO leads to higher air temperature anomalies over eastern North America, reducing snow volume; positive phase of the winter PNA leads to stronger East Asian jet stream and thus more snowfall in Northwest America but to less snowfall in Northeast America and near the Arctic coast (which is in agreement with the regression map between snow depth and PNA, Figure 5d).

To better examine inter-annual co-variability between snow volume and climate indices, the linear trend in all time series has been removed and the corresponding correlation coefficients have been computed (Table 3). Over Eurasia, the most statistically significant correlations are still obtained with the AO (correlation=-0.51 with p-value=0.03) and PDO (correlation=0.41 with p-value=0.09) indices. Over North America, the only statistically significant correlation is obtained with PNA index (correlation=-0.58 with p-value=0.01), which means that North American snow volume linear trend could be related to the AO and AMO indices (probably due to the shift from a negative to a positive AMO around 1995), whereas snow volume inter-annual variability is more linked with the PNA pattern.

The correlation coefficients between yearly mean snow volume anomaly and all possible linear combinations of two different climate indices have also been computed. For Eurasia, the best correlation coefficient (0.68, p-value=0.002) is obtained for a linear combination between AO and AMO (-156.AO-611.AMO, blue curve in Figure 2b), whereas for North America the best correlation coefficient (0.75, p-value<0.01) is obtained for a linear combination between PNA and AMO (-235.AMO-92.PNA, blue curve in Figure 2c). AO and PNA, when combined with AMO, influence respectively the most Eurasia and North America snow volume and represent regional atmospheric processes influencing the two continents. If the trend is removed from both climate indices and snow volume, the best correlation coefficient is still obtained with a linear combination of AO and AMO indices over Eurasia (-136.AO-940.AMO, correlation=0.68 and p-value=0.002). Surprisingly, linear combination of AO and PDO indices, which individually gives respectively the first (-0.51) and second (0.41) best correlations with snow volume over Eurasia, only corresponds to the second best correlation between a linear combination of two climate indices and snow volume . Over North America, the best correlation is obtained with a linear combination of , however the second best correlation is obtained with AMO and PNA indices (-165.AMO-87.PNA, correlation=0.61 and p-value=0.007).

Snow and sea level

High latitude snow has a large impact on river discharge and thus is the main source of fresh water input to the Arctic Ocean. Snow volume change (V snow ) presented in the previous sections can be used to estimate the snow contribution to the global mean sea level (SLV snow ), using equation (1).

snow ocean water snow snow V A SLV ⋅ ⋅ - = ρ ρ (1)
where ρ snow =300 kg.m -3 (snow density), ρ water =1000 kg.m -3 (liquid water density) and A ocean =3.6x10 8 km 2 (total oceanic domain).

Over 1989/2006, the snow volume trend from SSM/I converted into equivalent sea level is very small (-0.0013±0.0087 mm.yr -1 ) and not statistically significant. The snow volume trend over the altimetry time span (1993/2006) The mean sea level data have been corrected for steric effects (mainly thermal expansion of ocean waters) using the methodology developed by [START_REF] Llovel | Regional and interannual variability in sea level over 2002-2009 based on satellite altimetry, Argo float data and GRACE ocean mass[END_REF] and based on Argo data [START_REF] Guinehut | On the use of satellite altimeter data in Argo quality control[END_REF]. Mean and trend have been removed from this corrected sea level and the seasonal cycle (i.e. the sinusoid with a 365.25 day period which best fits the time series) has been least square adjusted. The sea level seasonal cycle has a maximum amplitude of 6.2 mm which occurs around 15 October. Snow volume converted into sea level (mean and trend removed) has a seasonal cycle with 4.1 mm maximum amplitude around 10 August. Its amplitude is smaller, but has the same order of magnitude than the global mean sea level seasonal cycle, yet breaks earlier. This phase lag could be due to the time taken by water from snow melt to be routed to the ocean by the river network. It could also be explained if water stored in other reservoirs like ground water and water vapor in the atmosphere is taken into account. To test this hypothesis, the seasonal cycle of the ground water has been approximated by a sinusoid with amplitude of 3 mm (in sea level equivalent) and a yearly maximum at the beginning of September [START_REF] Cazenave | Global ocean mass variation, continental hydrology and the mass balance of Antarctica ice sheet at seasonal time scale[END_REF]. Similarly, the atmospheric water vapor has been approximated by a sinusoid with amplitude of 2 mm and a yearly maximum at the beginning of December [START_REF] Cazenave | Global ocean mass variation, continental hydrology and the mass balance of Antarctica ice sheet at seasonal time scale[END_REF]. The sum of these three contributors (snow, ground water and water vapor) has an amplitude equal to 6.9 mm, which is very close to the global mean sea level seasonal cycle amplitude with a reduced phase lag (the maximum is around mid-September). This is a surprisingly good result given the large approximation used and clearly shows that Arctic snow variability is one of the main contributors to the global mean sea level seasonal cycle, as previously shown from model outputs by [START_REF] Chen | Seasonal global water mass budget and mean sea level variations[END_REF], [START_REF] Minster | Annual cycle in mean sea level from Topex-Poseidon and ERS-1: inference on the global hydrological cycle[END_REF], [START_REF] Cazenave | Global ocean mass variation, continental hydrology and the mass balance of Antarctica ice sheet at seasonal time scale[END_REF] and [START_REF] Milly | Contribution of climate-driven change in continental water storage to recent sea-level rise[END_REF].

Conclusions and perspectives

From passive microwave data acquired between 1989 and 2006, it has been possible to estimate the high latitude snow volume variability. Over Eurasia, the mean annual snow volume trend is positive (11.3±9.3 km 3 .year -1 ), yet not statistically significant. Over North America the snow volume trend is equal to -9.7±3.8 km 3 .year -1 and is statistically significant.

This difference between the two continents could be due to AO which correlates with North American snow volume (correlation = 0.51) and anti-correlates with Eurasian snow volume (correlation = -0.57). These differences are also linked with regional climatic conditions as snow volume anomaly over North America better (anti-)correlates with the PNA index (correlation=-0.66), and AMO index (correlation=-0.59). However, the correlation between AMO and North American snow volume is mainly due to the trend and not to the interannual variability. Moreover, snow volume over Eurasia correlates well with a linear combination of the AO and AMO indices (correlation = 0.68), whereas over North America it correlates with a linear combination of the PNA and AMO indices (correlation = 0.75).

Finally, this study shows that high latitude snow volume does not contribute to the global mean sea level trend observed by satellite altimetry, but is a main component of the global mean sea level seasonal cycle.

In the future, it will be interesting to compare the snow volume trends observed in this study and trends from other hydrologic parameters in order to better understand the interaction between snow and the whole North Hemisphere high latitude hydrological cycle.

  the mean, standard deviation and trend (with the p-value) of the SSM/I-based snow volume over 1989/2006.

  presents the regional distribution of snow depth trends over1989/2006 (only statistically significant trends, i.e. p-value<0.1, are shown). Over Eurasia, the highest positive trends are found over the Lena basin, the southern part of the Yenisey basin (between 60°N/70°N and 100°E/130°E) and Eastern Europe. Over North America, Quebec, Baffin Island and the Arctic Ocean coast show negative trends, whereas positive trends are found over the Rockies and Southern Alaska. Previous studies interpolated sparse in situ measurements to infer temporal evolution of snow cover.[START_REF] Groisman | State of the ground: climatology and changes during the past 69 years over Northern Eurasia for a rarely used measure of snow cover and frozen land[END_REF] analyzed 1811 in situ observations of the soil condition (classified as frozen or unfrozen) between 1956 and 2004, within 1197 1°x1° grid cells over the former Soviet Union (most of these grid cells containing only one in situ station). The in situ network used is very sparse above 55°N and East of the Ural Mountains.Groisman et al. (

  Data have been computed by the National Oceanic and Atmospheric Administration (NOAA)/Climate Prediction Center (CPC), the NOAA/Earth System Research Laboratory (ESRL) and the Joint Institute for the Study of the Atmosphere and Ocean (JISAO)/University of Washington (UW), and have been downloaded from http://ioc3.unesco.org/oopc/state_of_the_ocean/atm. Each index has been averaged from January through March for each year.

Figure 4

 4 Figure 4 presents annual snow volume time series over Eurasia and North America

Figure captions : Figure 1 .

 :1 Figure captions:

Figure 2 .

 2 Figure 2. Annual snow volume from SSM/I data (solid red curve) over the whole study

Figure 3 .

 3 Figure 3. Map of the annual snow depth trends over the 1989/2006 time span. Only

Figure 4 .

 4 Figure 4. Annual snow volume over Eurasia (red curves in a., b., c., d.) and over North

Figure 5 .

 5 Figure 5. Regression maps between annual snow volume and January to March average of

Figure

  Figure 1:

Table captions :

 captions Table 1. Annual snow volume trends over 2003/2006 computed from SSM/I snow volume, snow reservoir extracted from GRACE and from the total GRACE signal over the whole study domain, Eurasia and North America

Table 2 .

 2 Mean, standard deviation and trend of annual snow volume retrieved from SSM/I measurements during 1989/2006 period over the whole study domain, Eurasia and North America; p-values (p) for trends are indicated in parentheses.

Table 3 .

 3 Correlation coefficients from January to March average of the climate indices (AO, AMO, PDO and PNA) and annual snow volumes over North America and Eurasia (the pvalue of each correlation is indicated in parentheses). The correlations are computed both with and without trend in the time series of both snow volume and climate indices.

Table 1 :

 1 Annual snow volume trend for 2003/2006 (km 3 .year -1 )

	Tables:
		SSM/I		GRACE snow		GRACE total
				no PGR corr.	PGR corr. no PGR corr.	PGR corr.
	Whole domain -71.4±83.5 (p=0.48)		46.8±103.4 (p=0.70)	-448.8	179.0±234.3 (p=0.52)	-316.6
	Eurasia	-20.0±82.8 (p=0.83)		-91.3±57.5 (p=0.25)	-122.5	-237.2±231.4 (p=0.41)	-268.4
	North America -51.4±9.6 (p=0.03)		138.1±46.2 (p=0.10)	-326.3	416.3±52.0 (p=0.02)	-48.1
	Table 2:
				SSM/I annual snow volume
			Whole domain	Eurasia		North America
	1989/2006 mean (km 3 )	3713		2272		1441
	Std (km 3 )		218		189		95
	Trend (km 3 .year -1 )	1.5±10.5 (p=0.88)	11.3±9.3 (p=0.25)	-9.7±3.8 (p=0.02)
	Table 3 :
				AO	AMO	PDO	PNA
	Eurasian snow	with trend	-0.57 (p=0.01)	0.04 (p=0.87)	0.49 (p=0.04)	0.30 (p=0.22)
	depth	without trend	-0.51 (p=0.03)	-0.33 (p=0.18)	0.41 (p=0.09)	0.21 (p=0.39)
	North American	with trend	0.51 (p=0.0 3)	-0.59 (p=0.01)	-0.18 (p=0.47)	-0.66 (p=0.003)
	snow depth	without trend	0.25 (p=0.31)	-0.31 (p=0.20)	0.10 (p=0.71)	-0.58 (p=0.01)

Acknowledgments

The NSIDC is greatly thanked for processing and freely distributing SSM/I data.

The GLDAS data used in this study were acquired as part of the mission of NASA's Earth Science Division and archived and distributed by the Goddard Earth Sciences (GES) Data and Information Services Center (DISC). The Authors are thankful to J.-P. Vergnes for downloading and post-processing GLDAS data.

ECMWF ERA-Interim data used in this study have been obtained from the ECMWF data server.

The NOAA/CPC, NOAA/ESRL and JISAO/UW are acknowledged for letting the time series of the climate indices used in this study freely available to the community.

The authors are particularly grateful to two anonymous reviewers for their constructive comments and suggestions, which significantly improved the quality of the manuscript.

Three of the authors are supported by a CNES/Noveltis PhD grant (S. Biancamaria), a STAE foundation grant in the framework of the CYMENT project (F. Frappart) and a CNRS/Région Midi-Pyrénées PhD grant (W. Llovel).