Francois Lacan

Kazuyo Tachikawa

Catherine Jeandel

Neodymium isotopic composition of the oceans: a compilation of seawater data

Keywords: neodymium isotopes, GEOTRACES, seawater, neodymium concentrations, Rare Earth Elements

A global compilation of the neodymium isotopic composition of seawater is presented. With 880 data points, it confirms the gradual  Nd increase for intermediate/deep water masses from the northwest North Atlantic, via the Austral and Indian oceans, to the Pacific. This confirms the usefulness of this tracer for studying large scale oceanic circulation. The compilation stresses the need for documenting the oceans south of 30°S, from which less than 4% of the data are derived. The associated neodymium concentrations display heterogeneous vertical gradients amongst major oceanic basins. In addition to particle remineralization along the global thermohaline circulation, the database suggests that basin size differences also contribute to that heterogeneity.

INTRODUCTION

The neodymium isotopic composition, hereafter abbreviated Nd IC, is expressed as , where CHUR stand for CHondritic Uniform Reservoir and represents a present day average Earth value ( 143 Nd/ 144 Nd) CHUR =0.512638 [START_REF] Jacobsen | Sm-Nd isotopic evolution of chondrites[END_REF]. The Nd IC is heterogeneous on the continents, ranging from  Nd =-56 in old granitic cratons to +12 in recent mid-oceanic ridge basalts (GEOROC database: http://georoc.mphmainz.gwdg.de/georoc/), reflecting both the heterogeneous ages and Sm/Nd ratios of the geological formations. It is therefore used as a tracer of the origin of the terrigeneous matter, such as aerosols [START_REF] Grousset | Antarctic (Dome C) ice-core dust at 18 k.y. B.P.: Isotopic constrains on origins[END_REF] and river particles [START_REF] Mearns | A samarium-neodymium isotopic survey of modern river sediments from Northern Britain[END_REF]. In the ocean, the Nd is of lithogenic origin. Its concentration is of the order of 10 -12 g/g (10 pmol/kg). It is predominantly found in the dissolved form (90 to 95 %, [START_REF] Sholkovitz | Ocean particle chemistry: The fractionation of rare earth elements between suspended particles and seawater[END_REF][START_REF] Jeandel | Exchange of Nd and its isotopes between seawater small and large particles in the Sargasso Sea[END_REF][START_REF] Alibo | Rare earth elements in seawater: particle association, shalenormalization, and Ce oxidation[END_REF]). The oceanic Nd IC distribution is heterogeneous, controlled both by Nd sources and sinks, and by the oceanic circulation. Besides river and aeolian inputs, Boundary Exchange -dissolution of Nd from margin sediments associated to removal of dissolved Nd from seawater -seems to dominate Nd sources and sinks [START_REF] Tachikawa | Neodymium budget in the modern ocean and paleoceanographic implications[END_REF]Lacan and Jeandel, 2005b;[START_REF] Arsouze | Reconstructing the Nd oceanic cycle using a coupled dynamical -biogeochemical model[END_REF]. The major process involved in Nd removal is considered to be adsorption onto organic and Fe-Mn oxide coatings of sedimenting particles [START_REF] Sholkovitz | Ocean particle chemistry: The fractionation of rare earth elements between suspended particles and seawater[END_REF][START_REF] Haley | Rare earth elements in pore waters of marine sediments[END_REF]. The oceanic residence time of Nd is of the order of 300 to 1000 years [START_REF] Tachikawa | Neodymium budget in the modern ocean and paleoceanographic implications[END_REF][START_REF] Arsouze | Reconstructing the Nd oceanic cycle using a coupled dynamical -biogeochemical model[END_REF], which is shorter than the mean oceanic mixing time. These features have led to use the Nd IC as a tracer of oceanic circulation for the modern and past ocean [e.g. [START_REF] Piepgras | Rare earth element transport in the western North Atlantic inferred from isotopic observations[END_REF][START_REF] Rutberg | Reduced North Atlantic deep Water flux to the glacial Southern Ocean inferred from neodymium isotope ratios[END_REF][START_REF] Piotrowski | Oscillating glacial northern and southern deep water formation from combined neodymium and carbon isotopes[END_REF] and dissolved/particles interactions [e.g. [START_REF] Tachikawa | Distribution of rare earth elements and neodymium isotopes in settling particulate material of the tropical Atlantic Ocean (EUMELI site)[END_REF]. Although the usefulness of Nd IC in oceanography has been recognized since the early 80s, global modeling of the Nd IC of the oceans has became possible only for a few years because of the improvement of our understanding of the Nd oceanic cycle and of an increase of the global oceanic database [e.g. [START_REF] Arsouze | Modeling the neodymium isotopic composition with a global ocean circulation model[END_REF]. Modeling is a powerful approach to examine the sensitivity of Nd IC to the various processes that are still difficult to quantify. Model-data comparison is one of the most efficient approaches that allow better constraining Nd behavior in the marine environment. Besides, Nd IC application to paleoceanographic studies necessitates the validation that marine archives faithfully record seawater isotopic signatures.

In such a context, databases of the Nd IC of i) the lithogenic inputs to the ocean and ii) the seawater itself appeared useful. A database of the lithogenic inputs to the ocean was published a few years ago [START_REF] Jeandel | Isotopic Nd compositions and concentrations of the lithogenic inputs into the ocean: A compilation, with an emphasis on the margins[END_REF]. We have maintained an up-to-date and freely available database of seawater Nd IC measurements on the LEGOS marine geochemistry group website (http://www.legos.obs-mip.fr/recherches/equipes/geomar/Nd_DataBase). In the present paper, we publish this database as of September 1 rst 2011, we briefly discuss the methodology and intercalibration issues and finally we detail the most significant features of this global data compilation. We believe that this compilation could be useful, notably it could be used as a reference for modeling and paleoceanographic studies.

Nd IC has been identified as a "key parameter" of the GEOTRACES program (Scientific Committee on Oceanographic Research, [GEOTRACES Planning Group, 2006]). The number of seawater Nd IC data should therefore significantly increase in the coming years. These new and forthcoming data, acquired in the framework of GEOTRACES labeled cruises, will be archived in the GEOTRACES International Data Assembly Centre (http://www.bodc.ac.uk/geotraces/).

METHODOLOGY

All seawater Nd IC data published up to September 2011 to our knowledge are reported here [START_REF] Piepgras | The isotopic composition of Nd in different ocean masses[END_REF][START_REF] Piepgras | Neodymium isotopic variations in seawater[END_REF][START_REF] Piepgras | Isotopic composition of neodymium in waters from the Drake Passage[END_REF][START_REF] Piepgras | Influence of the Mediterranean outflow on the isotopic composition of neodymium in waters of the North Atlantic[END_REF][START_REF] Stordal | Neodymium isotopic study of Baffin Bay water: sources of REE from very old terranes[END_REF][START_REF] Piepgras | Rare earth element transport in the western North Atlantic inferred from isotopic observations[END_REF][START_REF] Piepgras | The isotopic composition of neodymium in the North Pacific[END_REF][START_REF] Spivack | Neodymium isotopic composition of the Mediterranean outflow and the eastern North Atlantic[END_REF][START_REF] Bertram | The geochemical balance of the rare earth elements and Nd isotopes in the oceans[END_REF][START_REF] Jeandel | Concentration and isotopic composition of neodymium in the South Atlantic Ocean[END_REF][START_REF] Henry | Particulate and dissolved Nd in the Western Mediterranean Sea: sources, fates and budget[END_REF][START_REF] Shimizu | Cerium and neodymium ratios and REE patterns in seawater from the North Pacific Ocean[END_REF][START_REF] Jeandel | Concentrations and Isotopic compositions of Nd in the Eastern Indian Ocean and Indonesian Straits[END_REF][START_REF] Tachikawa | A new approach to Nd residence time in the ocean: the role of atmospheric inputs[END_REF][START_REF] Amakawa | Nd isotopic and REE pattern in the surface waters of the eastern Indian Ocean and its adjacent seas[END_REF][START_REF] Lacan | Tracing Papua New Guinea imprint on the central Equatorial Pacific Ocean using neodymium isotopic compositions and Rare Earth Element patterns[END_REF][START_REF] Lacan | Masses d'eau des Mers Nordiques et de l'Atlantique Subarctique tracées par les isotopes du néodyme[END_REF]Amakawa et al., 2004aAmakawa et al., , 2004b;;Lacan and Jeandel, 2004a, 2004b, 2004c;[START_REF] Tachikawa | Neodymium isotopes in the Mediterranean Sea: Comparison between seawater and sediment signals[END_REF][START_REF] Vance | The use of foraminifera as a record of the past neodymium isotope composition of seawater[END_REF][START_REF] Dahlqvist | The concentration and isotopic composition of diffusible Nd in fresh and marine waters[END_REF][START_REF] Lacan | Masses d'eau des Mers Nordiques et de l'Atlantique Subarctique tracées par les isotopes du néodyme[END_REF]Jeandel, 2005a, 2005b;Tazoe et al., 2007aTazoe et al., , 2007b;;[START_REF] Andersson | Neodymium isotopes in seawater from the Barents Sea and Fram Strait Arctic-Atlantic gateways[END_REF][START_REF] Amakawa | Nd isotopic composition in the central North Pacific[END_REF][START_REF] Godfrey | Hafnium and neodymium isotope variations in NE Atlantic seawater[END_REF][START_REF] Porcelli | The distribution of neodymium isotopes in Arctic Ocean basins[END_REF][START_REF] Rickli | The hafnium-neodymium isotopic composition of Atlantic seawater[END_REF]Zimmermann et al., 2009aZimmermann et al., , 2009b;;[START_REF] Rickli | Hafnium and neodymium isotopes in surface waters of the eastern Atlantic Ocean: Implications for sources and inputs of trace metals to the ocean[END_REF][START_REF] Copard | Nd isotopic composition of water masses and dilution of the Mediterranean outflow along the southwest European margin[END_REF][START_REF] Tazoe | Coupled isotopic systematics of surface cerium and neodymium in the Pacific Ocean[END_REF]. When Nd concentration is also available, it is also reported. In contrast, studies reporting seawater Nd concentration only (without Nd IC) are not considered in this study. Seawater is defined here as waters with salinities >25. Both unfiltered and filtered sampled are included (authors usually assume that unfiltered samples well represent dissolved Nd, as particulate Nd usually accounts for only ~5% of the total; cf. introduction).

Since the early work of G.J. Wasserburg and his collaborators, various protocols have been used to measure seawater Nd IC. They are briefly summarized. They include four main steps: the sampling (and sometimes filtration), the preconcentration, the purification and the analysis of Nd isotopic ratios. Sampling has been performed using standard Niskin bottles, trace metal clean Go-Flo bottles, or various surface water sampling devices, such as inlets from the ship hulls, towed fishes or snorkels. The preconcentration step has been performed by Fe or Mn REE co-precipitation, REE complexation by a HDEHP-H 2 MEHP mixture adsorbed on a C 18 cartridge, or REE adsorption on MnO 2 fibers. One study used Nd diffusion through a layer of poly acrylamide hydrogel (DGT), in situ, allowing sampling and preconcentration to be performed simultaneously. MnO 2 fibers were sometimes also deployed in situ, also allowing combining these two steps. The purification has been performed by various ion exchange chromatographies. Finally Nd isotope ratios have been determined by multi-collector thermal ionization and inductively coupled plasma mass spectrometry.

Thanks to the recent intercalibration conducted on a world scale in the framework of the GEOTRACES program, most of these analytical procedures have been validated and a detailed review of those currently in use is given in [van de Flierdt et al., submitted].

The intercalibration exercise revealed that the inter-laboratory reproducibility for dissolved  Nd measurement is 0.5 to 0.6  Nd units (2 standard deviation)[van de Flierdt et al., submitted], that is well below the natural variation range in the ocean (cf. below). Most of the laboratories having reported seawater Nd IC data took part to the Geotraces intercalibration exercise. Nevertheless, the pioneer group of G.J Wasserburg and his co-workers could not take part to this recent exercise. Therefore, we compared our own data with those published by this pioneer group. Because surface waters are more susceptible of temporal variations than intermediate and deep waters, the following comparisons are restricted to the latters. The signature of the core of the North Atlantic Deep Water first characterized by  Nd =-12.6 by [START_REF] Piepgras | Rare earth element transport in the western North Atlantic inferred from isotopic observations[END_REF] (TTO/TAS St. 63, 8°N, 2910m), was then characterized by  Nd =-12.3 by [START_REF] Jeandel | Concentration and isotopic composition of neodymium in the South Atlantic Ocean[END_REF] (SAVE St. 302, 33°S, 2763m). This consistency was later underlined on the basin scale [START_REF] Von Blanckenburg | Tracing Past Ocean Circulation?[END_REF] and further confirmed by the GEOTRACES intercalibration data,  Nd =-12.14 [van de Flierdt et al., submitted] (BATS, 32°N, 2000m). The data published by [START_REF] Spivack | Neodymium isotopic composition of the Mediterranean outflow and the eastern North Atlantic[END_REF] in the North-eastern Atlantic for the North Atlantic Central Water and in the Gibraltar Strait for the Mediterranean Water compare also well with the values published by [START_REF] Tachikawa | A new approach to Nd residence time in the ocean: the role of atmospheric inputs[END_REF][START_REF] Tachikawa | Neodymium isotopes in the Mediterranean Sea: Comparison between seawater and sediment signals[END_REF] (cf. Table 1 and Fig. 1). Finally, an excellent agreement was also found in the Labrador Sea Water, the Iceland Scotland Overflow Water, the Norwegian Sea Deep Water and the Greenland Sea Deep Water, when comparing Piepgras and Wasserburg data [1987] with Lacan and Jeandel data [2004a, 2004b, 2004c, 2005a], despite the fact the latter samples were taken 16 to 18 years later than the former (cf. Table 1 and Fig. 1). The good consistency between all these data made us confident to include them in the present database. , 1987]; b [Lacan and Jeandel, 2005a]; c [Lacan and Jeandel, 2004b]; d [Lacan and Jeandel, 2004a]; e [START_REF] Spivack | Neodymium isotopic composition of the Mediterranean outflow and the eastern North Atlantic[END_REF]; f [START_REF] Tachikawa | A new approach to Nd residence time in the ocean: the role of atmospheric inputs[END_REF]; g [START_REF] Tachikawa | Neodymium isotopes in the Mediterranean Sea: Comparison between seawater and sediment signals[END_REF] [START_REF] Schlitzer | Ocean Data View[END_REF].

RESULT AND COMMENTS

The database is presented in Table 2 (supplementary material). Figure 2 shows the geographical distribution of the data. The database includes a total of 880  Nd data points, from 280 different stations, and 686 Nd concentrations. Out of these 280 stations, only 110 (i.e. 39%) include more than 3  Nd data points (Fig. 2b). About two thirds (66%) of the  Nd data are located in the top 1000 m of the water column (Fig. 3). The global mean Nd IC is  Nd =-8.8, with values ranging from -26.6 to +0.4. The latitudinal distribution is very uneven, with 84% of the  Nd data in the Northern hemisphere, and very few data south of 30°S (less than 4%). Most data are found in the North Atlantic and North Pacific oceans, whereas the South Atlantic, the Indian, the Austral and the South Pacific oceans are very scarcely documented (see also Supplementary Fig. 1).

Figure 4 displays the Nd IC values averaged between the surface and 400m depth and between 800m depth and the bottom. Although slightly different, the two distributions display the same trends on the first order. Values are minimal in the northwest Atlantic, especially in the Baffin Bay, intermediate in the Austral and Indian oceans, and maximal in the Pacific. This is also visible in Fig. 3 and Supplementary Fig. 1. This global pattern is in good agreement with the Nd isotopic signature of the continents close to the coasts, non-radiogenic (i.e. relatively depleted in Sm borne 143 Nd, therefore with relatively low  Nd values) around the Atlantic Ocean and radiogenic (i.e. relatively enriched in Sm borne 143 Nd, therefore with relatively high  Nd values) around the Pacific [Jeandel et al., 2007]. It confirms what was previously known about the global Nd isotope cycle [e.g. [START_REF] Piepgras | The isotopic composition of Nd in different ocean masses[END_REF][START_REF] Goldstein | Long lived Isotopic Tracers in Oceanography, Paleoceanography, and Ice sheet dynamics[END_REF]. Neodymium concentrations are displayed in Fig. 3 and Supplementary Fig. 1. Surface values are very variable. A few data points (13; less than 2% of the total number of data) display extremely high values (up to 152 pmol/kg), which likely reflect external inputs in coastal areas (e.g. Fjords). Excluding those, as a general rule, Nd concentration increases with depth. This increase contrasts with the relative homogeneity of the Nd isotopic composition (Fig. 3 and Supplementary Fig. 1). The Nd concentration vertical gradient is maximum in the Pacific, less pronounced in the Austral and Indian, much more moderated in the Atlantic, and not the dominant feature in the Mediterranean and the Arctic. The Nd concentration vertical gradient increase along the global thermohaline circulation has been attributed to the cumulative effect of particle remineralization as deep water age (as for nutrients) [START_REF] Siddall | Towards explaining the Nd paradox using reversible scavenging in an ocean general circulation model[END_REF][START_REF] Arsouze | Reconstructing the Nd oceanic cycle using a coupled dynamical -biogeochemical model[END_REF]. The present compilation also suggests a possible impact of the size of each basin (i.e. ratio of margin length to basin surface area), itself being related to the intensity of external inputs to each basin. In other words, the vast Pacific Ocean is less subjected to external inputs than smaller basins (e.g. the Mediterranean and the Arctic), which likely induces the lower concentrations observed in the upper part of the water column in the Pacific than in closer basins.

Figure 6 shows the temporal distribution of the  Nd data. It shows a relatively stable and limited number of data produced in the 80s and 90s, followed by a sharp increase in the last decade. This illustrates the pioneer nature of the work conducted during these two first decades, followed a drastic increase in the interest of the oceanographic community for this tracer, which is now recognized as a core parameter of the international marine geochemical program GEOTRACES. This trend will therefore very likely go on in the framework of the GEOTRACES program. As the GEOTRACES International Data Assembly Centre will take over the management of the oceanic Nd isotope database, we will keep this pre-GEOTRACES database available on the internet. Any contribution from colleagues to correct or complement it would be appreciated. The geographical distribution of the data described above suggests that, in the present context, documenting the southern oceans (in particular south of 30°S) is a high priority (cf. Fig. 2b). In that respect, major improvements are expected in the years to come, resulting from GEOTRACES and the 2007-2008 International Polar Year cruises (see the GEOTRACES International Data Assembly Centre web site for maps showing the location of these cruises).

Fig. 1 :

 1 Fig. 1: Locations of the stations used for comparing the data published by Wasserburg's group (black circles), with those published by the authors of the present work (white circles). Where 2 stations are too close to each other to be represented by 2 distinct symbols, a black and white symbol is used. Cf. Tab. 1 for  Nd data. Figure made with Ocean Data View[START_REF] Schlitzer | Ocean Data View[END_REF].

  Fig. 3. All  Nd data as a function of depth (left), and corresponding Nd concentration when documented (right). Colors indicate locations, Blue: Atlantic and Arctic; Green: Mediterranean; Yellow: Austral and Indian; Red: Pacific. The Nd concentration scale is restricted to values lower than 65 pmol/kg for clarity, excluding 13 shallow near shore data (< 2% of the total) from the North Atlantic. Figure made with Ocean Data View[START_REF] Schlitzer | Ocean Data View[END_REF].

Fig. 4 .

 4 Fig. 4.  Nd averaged between the surface and 400m depth (top) and between 800m depth and the bottom (bottom). Figure made with Ocean Data View [Schlitzer, 2009].

Figure 5 Fig. 5 .

 55 Figure 5 displays a statistical distribution of  Nd values in distinct oceanic regions. This is an updated version of a figure published in 2005 [Lacan and Jeandel, 2005b], with 61% more data points, and some more details. Overall, this distribution confirms the preceding one. It includes a new basin, the Arctic Ocean, which was almost not documented in 2005. The Arctic distribution is strikingly similar to that of the Nordic Seas. Again, globally, this figure illustrates the gradual  Nd increase from the northwest North Atlantic, via the Austral and Indian oceans, to the Pacific, confirming the usefulness of this tracer for studying large scale oceanic circulation.

  Fig. 6. Number of  Nd data as a function of cruise and publication years. Note that cruise years are not always reported in the publications.

Table 1 .

 1 Comparison of data from G.J Wasserburg and his co-workers with data produced at LEGOS (Toulouse, France)

	Cruise, station	Depth. (m) Salinity  Nd 2 Conc. (pmol/kg)	Water mass	Ref .
	Hudson 83-036 St. 11; 52.08°N -47.02°E Signature/GINS St. 6; 50.20°N, -45.68°E	1500 1650	34.872 -14.4 0.4 34.86 -13.9 0.4	18.1 17.4	LSW	a b
	Hudson 83-036 St. 11; 52.08°N -47.02°E	2500 3000	34.932 -13.8 0.5 34.933 -13.4 0.4	16.7 17.3	NEADW	a a
	Signature/GINS St. 6; 50.20°N, -45.68°E	2499	34.90 -13.2 0.4	18.3		b
	TTO/NAS Station 142; 61.35°N -8.02°E Signature/GINS St. 23; 60.50°N, -5.00°E	750 988	34.905 -7.7 0.6 34.90 -7.3 0.4	21.4 24.6	ISOW	a c
	TTO/NAS Station 144; 67.67°N -3.28°E Signature/GINS St. 26; 69.03°N, 7.95°E	3750 2972	34.909 -9.5 0.5 34.90 -9.8 0.4	16.3 27.3	NSDW	a c
	TTO/NAS Station 149; 76.88°N 1.03°E	2800	34.901 -10.7 0.4	16.8		a
	Signature/GINS St. 30; 76.74°N, -2.33°E	2512 (unfiltered)	34.90 -10.8 0.4	19.6	GSDW	d
	TTO-TAS station 80; 27.83°N, -30°53°E	389	35.942 -11.0 0.8	13.9		e
	EUMELI Station E3O; 21°N, -31°E	250 500	36.47 -10.5 0.2 35.63 -10.5 0.2	18.6 18.2	NACW	f f
	Station MED-15; 36.07°N, -6°E	250 400	-38.517 -10.1 0.7 -9.9 0.5	28.0 32.4	MW	e e
	Station 5, 35.9°N, -5.65°E	350	38.48 -9.4 0.2	-		g
	Cf. Fig 1 for station locations.

References: a [Piepgras and Wasserburg

Acknowledgements: The CNRS (French National Center for Scientific Research) is thanked for supporting this study. We thank two anonymous reviewers for their constructive comments.

FIGURE CAPTIONS

Supplementary Fig. 1. All  Nd data as a function of depth, and corresponding Nd concentration when documented, for different oceanic areas. The Nd concentration scale is restricted to values lower than 65 pmol/kg for clarity, excluding 13 shallow near shore data (< 2% of the total) from the North Atlantic. Figure made with Ocean Data View [START_REF] Schlitzer | Ocean Data View[END_REF].


M F Azam

Patrick Wagnon
email: patrick@lgge.obs.ujf-grenoble.fr

Alagappan L Ramanathan

Christian Vincent

Parmanand Sharma

Yves Arnaud

Anurag Linda

Jose George Pottakkal

Pierre Chevallier

Virendra Bahadur Singh

Mohd Farooq

Etienne Berthier

From balance to imbalance: a shift in the dynamic behaviour of Chhota Shigri

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Although Himalayan glaciers have important social and economic impacts (e.g. Barnett and others, 2005), they have never been monitored on a long-term basis and little is known about recent glacier trends or their contribution to local and regional water supplies, even giving rise to a controversial statement in the IPCC 4 th assessment report saying that "the likelihood of them disappearing by the year 2035 or perhaps sooner is very high if the Earth keeps warming at the current rate" (Cogley and others, 2010). A general negative mass balance of mountain glaciers on a global level is clearly revealed from recent research (e.g. [START_REF] Cogley | Geodetic and direct mass-balance measurements: comparison and joint analysis[END_REF]Zemp and others, 2009) but the effect of global warming in the Himalaya is still under debate (Yadav and others, 2004;[START_REF] Roy | Analysis of trends in maximum and minimum temperature, diurnal temperature range, and cloud cover over India[END_REF]. Though temperate glacial mass-balance change is one of the best indicators of climate change [START_REF] Oerlemans | Glaciers and climate change[END_REF]Vincent and others, 2004;Ohmura and others, 2007) the paucity of mass-balance data in Himalaya makes it difficult to obtain a coherent picture of regional climate change impacts in this region. In the Indian Himalaya the first massbalance study started on Gara Glacier (Himachal Pradesh) in September 1974 (Raina and others, 1977) and ended in 1983 (Dobhal and others, 2008). According to [START_REF] Dyurgerov | Glaciers and the changing Earth System: a 2004 snapshot[END_REF], eight glaciers in the Indian Himalaya were surveyed for mass balance at least for one year during the 1980s. Unfortunately each study was restricted to short periods, not more than one decade (Dobhal and others, 2008). Remote sensing studies were also attempted in this part of Himalaya, but these studies either deal with only surface area changes (e.g. Kulkarni and others, 2007;Bhambri and others, 2011) or cover short periods [START_REF] Kulkarni | Mass balance of Himalayan glaciers using AAR and ELA methods[END_REF]Berthier and others, 2007).

The present study is based on mass balance and surface ice flow velocity measurements conducted on Chhota Shigri Glacier, Himachal Pradesh, India between 2002 and 2010, and on a comparison to data collected in 1987/1989. In the Indian Himalaya, this is one of the longest continuous field mass-balance dataset. Moreover, in October 2009, a Ground Penetrating Radar (GPR) survey was also conducted to measure ice thickness. Eight years of mass-balance measurements, surface ice velocities and ice thickness data provide an opportunity to study the behaviour of this glacier. The main objectives of this paper are (i) to present the recent mass balance of Chhota Shigri Glacier, (ii) to determine the ice fluxes at five cross sections from thickness and ice velocities, and (iii) to compare these data with the ice fluxes inferred from cumulative surface mass balance upstream of the same cross sections. These results give insights into the mass-balance trend of the glacier over the last two to three decades, and allow us to assess whether it is in equilibrium with the climate of the 21 st century.

Site description and methodology

Site description

Chhota Shigri Glacier (32.2 o N and 77.5 o E) is a valley-type glacier located in the Chandra-Bhaga River basin of Lahaul and Spiti Valley, Pir Panjal Range, Western Himalaya. This glacier extends from 6263 m to ~4050 m a.s.l., is ~9 km long and covers 15.7 km 2 of area.

The snout of this glacier is easy to locate from one year to another because it is well defined, lying in a narrow valley and giving birth to a single proglacial stream. The main orientation of this glacier is north except for its tributaries which have a variety of orientations (Fig. 1). The lower ablation area (<4400 m a.s.l.) is partly covered by debris representing ~3.4% of the total surface area. This glacier is located in the monsoon-arid transition zone and is influenced by two atmospheric circulation systems: the Indian monsoon during summer (July-September) and the northern-hemisphere mid-latitude westerlies during winter (January-April) (Singh and others, 1997;[START_REF] Bookhagen | Topography, relief, and TRMM-derived rainfall variations along the Himalaya[END_REF]Gardelle and others, 2011).

Mass balance

The first series of mass-balance measurements was performed on Chhota Shigri Glacier between 1987and 1989(Nizampurkar and Rao, 1992;Dobhal and others, 1995;[START_REF] Kumar | Chhota Shigri Glacier: its kinematic effects over the valley environment, in the northwest Himalaya[END_REF].

The bedrock topography and surface ice velocity were also surveyed over the same period by gravimetric and stake displacement methods respectively (Dobhal and others, 1995;[START_REF] Kumar | Chhota Shigri Glacier: its kinematic effects over the valley environment, in the northwest Himalaya[END_REF]. We re-initiated the mass-balance observations in 2002. From that year, annual surface mass-balance measurements have been carried out continuously on Chhota Shigri Glacier at the end of September or beginning of October using the direct glaciological method [START_REF] Paterson | The physics of glaciers[END_REF]. Ablation was measured through a network of ~22 stakes distributed between 4300 and 5000 m a.s.l. (Fig. 1) whereas in the accumulation area, the net annual accumulation was obtained at six sites (by drilling cores or pits) between 5100 and 5550 m a.s.l. (Wagnon and others, 2007). In the accumulation area, the number of sampled sites is limited due to the difficulty in access and high elevation. The glacier-wide mass balance, B a is calculated according to:

B a = ∑ b i (s i /S) (1)
Where b i stands for the mass balance of the altitudinal range (m w.e. yr -1 ), i, of map area s i and S symbolizes the total glacier area. For each altitudinal range, b i is obtained from the corresponding stake readings or net accumulation measurements.

Surface velocity

Annual surface ice velocities were measured at the end of each ablation season (September-October) by determining the annual stake displacements (~22 stakes) using a Differential Global Positioning System (DGPS). These geodetic measurements were performed in kinematic mode relative to two fixed reference points outside the glacier on firm rocks. The accuracy in x (easting), y (northing) and z (elevation) of each stake position is estimated at ±0.2 m depending mainly on the size of the hole in which the stake has been set-up. Thus the surface ice velocities measured from stake displacements have an accuracy of ±0.3 m yr -1 .

Ice thickness

Ground penetrating radar measurements were conducted in October 2009 to determine ice thickness on five transverse cross sections (Fig. 1) between 4400 and 4900 m a.s.l. A pulse radar system (Icefield Instruments, Canada) based on the Narod transmitter [START_REF] Narod | Miniature high-power impulse transmitter for radio-echo sounding[END_REF] with separate transmitter and receiver, was used in this study with a frequency centered near 4.2 MHz and antenna length of 10 m. Transmitter and receiver were towed in snow sledges along the transverse profile, separated by a fixed distance of 20 meters, and used to record measurements every 10 m. The positions of the receiver and the transmitter are known through DGPS measurements, within an accuracy of ±0.1 m. The speed of electromagnetic wave propagation in ice has been assumed to be 167 m µs -1 [START_REF] Hubbard | Field Techniques in Glaciology and Glacial Geomorphology[END_REF]. The field measurements were performed in such a way as to obtain reflections from the glacier bed located more or less in the vertical plane with the measurement points at the glacier surface, allowing the determination of the glacier bed in two dimensions. The surface of the bedrock was constructed as an envelope of all ellipse functions, which give all the possible reflection positions between sending and receiving antennae. Ice thickness was measured along four transverse profiles (profiles 1-4) on the main glacier trunk and one (profile 5) on a western tributary (Fig. 1).

Data analysis and results

Glacier-wide mass balance and mass-balance profile

The annual glacier-wide mass balance and cumulative mass balance of Chhota Shigri Glacier between 2002 and 2010 are plotted in Figure 2. The glacier-wide mass balance was negative except for three years (2004/2005; 2008/2009 and 2009/2010) The quantitative uncertainty associated with the glaciological mass balance requires a distinction between the accumulation zone and the ablation zone. In the accumulation zone, the surface mass-balance measurements were obtained from shallow boreholes (auger). Therefore, they are based on core length and density determination. In the ablation zone, the measurements have been carried out from ablation stakes. The overall error (standard deviation) on point measurements are estimated at 0.30 m w.e. and 0.15 m w.e. in the accumulation zone and in the ablation zone, respectively. The overall error comes from a variance analysis (Thibert and others;2008) applied to all types of errors (ice/snow density, core length, stake height determination, liquid-water content of the snow, snow height). Although conducted on a glacier in the Alps, the analysis of Thibert and others (2008) can be generalized to other glaciers because it is based on measurement errors which are similar on every glacier when using the glaciological method.

However, only 6 sites are sampled in the accumulation zone (11.6 km²), and 22 sites in the ablation zone (4.1 km²). The uncaptured spatial variability of surface mass balance may cause systematic errors on the glacier-wide mass balance. In the accumulation zone, the spatial variability remains unknown and is probably very high as observed for other glaciers (e.g. Machguth and others, 2006). In the ablation zone, stakes set up at the same altitude show similar values except on the terminal tongue which is debris covered (0.54 km²). Consequently, the overall uncertainties on mass-balance profile have been assessed at 0.5 m w.e. in the accumulation zone, 0.25 m w.e. in the white ablation zone and 0.5 m w.e. in the debris covered area of the glacier. Moreover, the surface area estimation also causes systematic error. The uncertainty on the surface area calculated for each altitudinal range is estimated at 5%.

Combining these errors at different altitudinal ranges using Equation ( 1), the uncertainty on the annual glacier-wide mass balance is 0.4 m w.e. yr -1 . As revealed by other studies (e.g. [START_REF] Vincent | Influence of climate change over the 20 th century on four French glacier mass balances[END_REF]Thibert and others, 2008;Huss and others, 2009), this estimation confirms that the glaciological method needs to be calibrated by a volumetric method over a long period of monitoring (i.e. >5 years) in order to limit the systematic errors and to improve the accuracy of absolute values of mass balance. Note that the uncertainty of relative changes in mass balance from year to year is smaller than those inherent in annual mass balances, given that the influence of systematic errors can be reduced.

We also calculated the mass-balance profile between 2002 and 2010 (Fig. 3). For each altitudinal range, we computed the average of all available measurements. Figure 3 reveals that melting in the lowest part of the ablation area (below 4400 m a.s.l.) is reduced by about 1 m w.e. yr -1 irrespective to its altitude. This is due to the debris cover (~5 to 10 cm thick debris mixed with isolated rocks) which reduces the melting in this region (Mattson and others, 1993;Wagnon and others, 2007). Moreover the lower part of Chhota Shigri Glacier flows in a north-south oriented deep and narrow valley (Fig. 1), causing the glacier tongue to receive less solar radiation due to the shading effect of the steep valley slopes.

Ice thicknesses and cross section areas

Thanks to clear reflections, ice-bedrock interface was generally easy to determine on all profiles. Figure 4 provides an example of the radargram obtained at cross section 2. A radar wave velocity of 167 m µs -1 was used for calculations of ice thickness at all the profiles. The cross sections obtained from GPR measurements reveal a valley shape with maximum ice thickness greater than 250 m (Fig. 5). The centerline ice thickness increases from 124 m at 4400 m a.s.l.

(cross section 1 in Fig. 1) to 270 m at 4900 m a.s.l. (cross section 4), which confirms that the thicknesses obtained by gravimetric methods in 1989 (Dobhal and others, 1995), twice lower than the present results, were under-estimated as proposed by Wagnon and others (2007). The cross sectional areas are given in Table 1. The accuracy of the calculated ice thickness is determined, in part, by the accuracy of the measurement of the time delays and the antenna spacing. Additional errors may arise because the smooth envelope of the reflection ellipses is only a minimal profile for a deep valley-shape bed topography, with the result that the ellipse equation will be governed by arrivals from reflectors located toward the side and thus not directly beneath the points of observation. Further errors may be introduced by assuming that all reflection points lie in the plane of the profile rather than on an ellipsoid. No errors associated to radar wave velocity variations between snow and ice have been accounted for because all cross sections were surveyed in the ablation zone or slightly above (with the firn-ice transition depth at the surface or < 2 m deep). Hence, the radar wave velocity for ice (167 m µs -1 ) was used to calculate all ice depths. The overall uncertainty in ice thickness is estimated as 15 m. Given that the uncertainty in ice surface coordinates is low (±0.1 m), the uncertainty on cross section areas mainly arises from the uncertainty in ice thickness. The uncertainties in cross section areas are 16, 9, 10, 10 and 15% for the cross sections 1, 2, 3, 4 and 5 respectively.

Ice velocity

Annual surface ice velocities were also measured between 2002 and 2010. However, some data gaps exist due to discontinuous DGPS signal, or loss of stakes. The ice velocities from 2003/2004 were used in this study because they provided the most complete dataset (Fig. 6). The center line horizontal ice velocities at each cross section were calculated by linear interpolation method along the center line between the velocities measured immediately upstream and downstream of the cross section (ablation stakes visible on Fig. 1). Mean cross section velocities are required to compute the ice fluxes (see section below). A map of the surface ice velocity field has been derived by correlating 2.5-m SPOT5 images acquired on 13 November 2004 and 21 September 2005 (Berthier and others, 2005). Comparison of the satellite-derived velocities with 16 nearly simultaneous DGPS velocity measurements shows a mean difference of 0.2 m yr -1 and a standard deviation of 1.6 m yr -1 . The ratio between the center line horizontal velocity and the mean surface velocity (all extracted from the satellite-derived 2004/2005 velocity field) was found to be 0.80 and 0.78 for cross sections 2 and 3 respectively. Reliable velocity measurements could not be measured from SPOT5 imagery for other cross sections. Using the mean value of 0.79, the mean horizontal velocity has been calculated from the center line velocity for each gate cross section (Table 1).

Ice fluxes from kinematic method

The ice flux Q (m 3 of ice per year) through each cross section was calculated using the cross sectional area S c (m 2 ) and depth-averaged horizontal ice velocity U (m yr -1 ).

Q = U S c
(2)

The depth-averaged horizontal ice velocity was derived from the mean surface ice velocity calculated in the previous section. [START_REF] Nye | The flow of a glacier in a channel of rectangular, elliptic or parabolic cross section[END_REF] gives ratios of depth-averaged horizontal ice velocity to mean surface ice velocity varying from 0.8 (no sliding) to 1 (maximum sliding). Here, we assume a mean basal sliding, with a constant ratio of 0.9. The calculated ice fluxes and maximum depth at each cross section are given in Table 1. The flux through cross section 3 at 4750 m a.s.l. is higher than the flux through cross section 4 at 4900 m a.s.l. This is due to the ice influx from the western part of the glacier (flux through cross section 5) which contributes to cross section 3 and not to cross section 4 (Fig. 1).

The largest uncertainty on the depth-averaged horizontal ice velocity results from the ratio between the depth velocity and the surface flow velocity. The estimated factor 0.9 and unknown variations in the basal sliding lead to an uncertainty of roughly ±10% in the calculated flux, which lies within the range of uncertainty of the other variables as discussed by Huss and others (2007). Consequently, we can assess that depth-averaged horizontal ice velocity at each cross section is known with an accuracy of 1.0 to 3.0 m yr -1 depending on the cross sections.

Combining these errors on the cross sectional area and mean velocity, the uncertainties on the ice fluxes are 0.21, 0.69, 0.92, 0.89 and 0.38 x 10 6 m 3 yr -1 for the cross sections 1, 2, 3, 4 and 5

respectively. Here we have considered that the errors are systematic, so these uncertainties are probably over-estimated.

Ice fluxes obtained from surface mass balance

We also calculated ice fluxes using annual surface mass balance measured during 2002/2010. Although the dynamic changes are neglected here, this method allows us to estimate the ice fluxes for each section from mass-balance data according to the following equation:

i z z i s b Q   max 9 . 0 1 (3)
Where Q is the ice flux (converted into m 3 of ice per year using an ice density of 900 kg m -3 , hence the factor 1/0.9) at a given elevation, z, and b i is the annual mass balance of the altitudinal range i of map area s i . The altitudinal ranges taken into account in the calculation are located between z and the highest range of the glacier z max (highest altitude of the glacier area contributing ice to the cross section). We assume here that on each point of the glacier above this altitude, z, the surface elevation has remained unchanged from one year to the next.

The ice fluxes calculated from annual mass-balance data at the 5 cross sections each year are given in Table 3, while the average ice fluxes for the eight years are given in Figure 7. The uncertainties on ice fluxes resulting from surface mass balance are directly derived from the mass-balance uncertainties (see section 3.1) applied to areas contributing to each cross section.

Discussion

The first and main objective of this section is to discuss the mass-balance change of Chhota Shigri Glacier over the last two to three decades using not only direct mass-balance observations (over the last eight years) but also ice-flux analysis. The second goal is to give insights into the specific dynamics and the future retreat of this glacier that can be expected in relation to its recent surface mass balance (hereafter referred to as SMB).

Null to slightly positive mass balance during the 1990s inferred from ice fluxes

The ice fluxes obtained from the kinematic method using ice thickness and 2003/2004 ice velocities are much higher than the average fluxes derived from the 2002/2010 SMBs, the latter being often negative (Table 2). Thus in this section, to assess the mean state of the glacier corresponding to the ice fluxes obtained by the kinematic method, we compare these measured ice fluxes to theoretical ice fluxes calculated from SMB assuming the glacier to be in steady state. The glacier-wide mass balance obtained by the glaciological method is -0.67 m w.e. yr -1 over the 2002/2010 period. Consequently, the SMB needs to be increased by 0.67 m w.e. yr -1 , for the glacier to be in steady state with the present surface area. For each year (2002/2010), we calculated the theoretical ice flux from SMB at each cross section assuming the glacier was in steady state. For this purpose, every year, a theoretical SMB at each elevation has been calculated by substracting the overall annual specific SMB of the same year. This result is also supported by other observations. First, the ice velocities measured in 1987/1988 (Dobhal and others, 1995)

Glacier dynamics starting to adjust to 21 st century negative SMB

In theory, the response of ice fluxes to surface mass balance is immediate (Cuffey and Paterson, 2010, p. 468) but observations show a 1-5 year delay (Vincent and others, 2000;[START_REF] Span | Simulating annual glacier flow with a linear reservoir model[END_REF]Vincent and others, 2009). For instance, [START_REF] Span | Simulating annual glacier flow with a linear reservoir model[END_REF] found synchronous decrease in ice velocity between eight glaciers in the Alps, which are driven by the same mass-balance changes (Vincent and others, 2005). Consequently, the recent dynamic behaviour of Chhota Shigri Glacier should be affected by the negative mass balance since 1999.

However, the stake network on Chhota Shigri Glacier, originally designed for SMB measurements, is not best suited to accurately compare either the ice velocities or the thickness variations because the measurements have not been performed exactly at the same location every year and they are mainly restricted to the ablation area.

In spite of the above limitation, an attempt has been made to compare ice velocities and elevations from the available stake network. For this purpose, stakes measured at the beginning and at the end of the series have been selected on five short longitudinal cross sections (A, B, C, D and E in Fig. 1) along the center line of the glacier where the network is most dense. The elevations in 2003 and 2010 and the ice velocities in 2003/2004 and 2009/2010 have been reported on these longitudinal cross sections to deduce thickness and velocity changes in the ablation area (Fig. 8, Table 4). Although the accuracy of the results is affected by the distance between the point measurements, we can conclude that the part of the glacier below 4750 m a.s.l.

is in strong recession. First, the thickness has decreased annually by 0.7 to 1.1 m yr -1 over the last seven years. Second, the ice velocities have been reduced by 7 m yr -1 between 2003 and 2010 resulting in a 24 to 37% decrease in the ice fluxes since 2003. Despite an improvable monitoring network, it may be surmised that the ice fluxes have been affected by the negative glacier-wide mass balance during (at least) the last eight years, and the dynamics of this glacier are progressively adjusting to the negative SMB. Consequently, we expect an accelerated terminus retreat in the coming years. If the SMB remained equal to its 2002/2010 average value in the future, the terminus would retreat by 5.6 km to reach an altitude of 4870 m a.s.l. (altitude where the ice flux is equal to zero) (Table 2).

Conclusion

The Chhota Shigri Glacier experienced negative mass balance over the 2002/2010 period.

The glacier-wide mass balance of the glacier is estimated at -0.67 m w. yr -1 between 1999 and 2004 using satellite images, our analysis suggests that the glacier experienced a period of slightly positive or close to zero mass balance at the end of the 20 th century, before starting to shrink. As Chhota Shigri seems to be representative of other glaciers in the Pir Panjal Range (Berthier and others, 2007), it is possible that many Western Himalayan glaciers of northern India experienced growth during the last 10-12 years of the 20 th century, before starting to shrink at the beginning of the 21 st century. Table 2: Ice fluxes (in 10 6 m 3 ice yr -1 ), inferred at each cross section from annual mass-balance data.

Table and figure captions

Table 3: Ice fluxes (in 10 6 m 3 ice yr -1 ), obtained at every cross section, using steady state massbalance assumption for every surveyed year. Altitude (m a.s.l.)

Flux derived from the dynamic method Flux derived from the SMB method, with bn = 0 Flux derived from the SMB method, with bn = -0.67 m w.e.

  . It varies from a minimum value of -1.40 m w.e. in 2002/2003 to a maximum value of +0.33 m w.e. in 2009/2010. The cumulative mass balance of Chhota Shigri is -5.37 m w.e. between 2002 and 2010 while the glacier-wide mass balance averaged over the same period is -0.67 m w.e. yr -1 .

  For instance, year 2002/2003 was characterized by a negative annual glacier-wide SMB of -1.40 m w.e. so we calculated a new SMB profile by adding +1.40 m w.e. to the SMB at each elevation. In contrast year 2009/2010 was characterized by a positive annual glacier-wide SMB of +0.33 m w.e. so we calculated a new SMB profile by subtracting 0.33 m w.e. from the SMB at each elevation. The resulting ice fluxes are reported in Table 3, together with the mean ice flux at each cross section over the eight years and the corresponding standard deviations. These ice fluxes are close to the 2003/2004 ice fluxes obtained by the kinematic method (Fig. 7) indicating that the dynamic behaviour of the glacier in 2003/2004 is representative for steady-state conditions. This suggests that in the years preceding 2003/2004, the glacier-wide mass balance of this glacier has probably been close to zero and that, in 2003/2004, the ice fluxes had not adjusted to previous year negative SMB.

  are very close to the 2003/2004 values (Fig. 6) suggesting that the dynamic behaviour of this glacier did not change a lot between 1988 and 2004. Second, the terminus fluctuation measured between 1988 and 2010 show a moderate retreat equal to 155 m, equivalent to only 7 m yr -1 , in agreement with conditions not far from steady state. Given that Berthier and others (2007) observed a glacier-wide SMB of Chhota Shigri Glacier of approximately -1 m w.e. yr -1 during the period 1999 to 2004, the glacier is likely to have experienced a null to slightly positive mass balance between 1988 and the end of the 20 th century.

  e. yr -1 between 2002 and 2010, revealing strong unsteady-state conditions over this period. Conversely, ice fluxes calculated through 5 transverse cross sections by the kinematic method correspond to near steady-state conditions before 2004. Given that ice velocities measured in 2003/2004 are close to those measured in 1988, and that terminus has retreated only 155 m between 1988 and 2010, it seems that the dynamic change was moderate between 1988 and 2004. Therefore, considering that Berthier and others (2007) observed a negative glacier-wide mass balance of about -1 m w.e.

  Since 2003, ice velocities and elevation are decreasing in the ablation area. Our data suggest that the ice fluxes have diminished by 24 to 37% below 4750 m a.s.l. between 2003 and 2010. Even if we account for a 37% decrease in ice fluxes calculated from 2003/2004 ice velocities to obtain present ice fluxes values, it remains a very large imbalance with ice fluxes coming from glacier-wide mass balance of the last eight years. Thus the present dynamics (thickness and ice velocities) of this glacier are far from surface mass balance and climate conditions of the last eight years, even if it is progressively adjusting. Therefore the glacier is likely to undergo accelerated retreat in the near future.This glacier is almost free of debris and thus its mass-balance variations are closely related to climate changes. This glacier has the longest running series of mass balance measurements in the Himalaya range. In the future, the dynamic behaviour and the mass balance need more detailed investigations, although such field measurements are demanding due to the very high altitude. In order to investigate the annual thickness and the ice velocity changes, we recommend performing elevation and ice velocity measurements on ~12 cross sections including some in the accumulation zone. We also recommend measuring the ice velocities from a dense network of stakes to be set up on longitudinal center lines in order to compare the annual velocity changes at the same points. Finally, we recommend calibrating and checking the massbalance field measurements from a volumetric method (from photogrammetry or remote sensing techniques). Zemp, M., M. Hoezle and W. Haeberli. 2009. Six decades of glacier mass-balance observations: a review of worldwide monitoring network. Ann. Glaciol., 50(50),101-111.

Figure 1 :

 1 Figure 1: Map of Chhota Shigri Glacier with the measured transverse cross-sections (lines 1 to

Figure 2 :

 2 Figure 2: Cumulative (line) and annual glacier-wide mass balances (histograms) of Chhota Shigri

Figure 3 :

 3 Figure 3: The 2002/2010 average mass-balance profile and the hypsometry of Chhota Shigri

Figure 4 :

 4 Figure 4: Radargram of cross section 2: radar signals plotted side by side from west to east in

Figure 5 :

 5 Figure 5: Ice depth and surface topography of cross-sections 1-5. The horizontal and vertical

Figure 6 :

 6 Figure 6: Measured ice velocities plotted as a function of the distance from the 2010 terminus

Figure 7 :

 7 Figure 7: Ice fluxes at every cross section derived (i) from 2003/2004 ice velocities and section

Figure 8 :

 8 Figure 8: Elevation (dots) and surface ice velocity (triangles) between 2003 (continuous lines)

  Figure 6

  Figure7

Figure 8

 8 Figure 8

Table 1 :

 1 Calculated ice flux, mean surface ice velocity and maximum ice depth at each cross section. The mean surface horizontal ice velocities are from DGPS measurements performed in 2003/2004. The satellite-derived mean ice velocities are from the correlation of satellite images acquired on 13 November 2004 and 21 September 2005 (NA: Not available).

Table 4 :

 4 Thickness and surface velocity changes between 2003 and 2010 on 5 longitudinal cross

sections (NA: Not available).

Acknowledgements

This work has been supported by the IFCPAR/CEFIPRA under the project n°3900-W1 and by the French Service d'Observation GLACIOCLIM as well as the Department of Science and Technology (DST), and Space Application Centre, Government of India. The French National Research Agency through ANR-09-CEP-005-01/PAPRIKA provided DGPS devices to perform field measurements. We thank J. E. Sicart, J. P. Chazarin, our field assistant Mr. B. B. Adhikari and the porters who have been involved in successive field trips, sometimes in harsh conditions.

We would also like to thank Dr D. P. Dobhal for his kind cooperation in answering our queries regarding the earlier research on Chhota Shigri Glacier. E. Berthier acknowledges support from the French Space Agency (CNES) through the TOSCA and ISIS proposal #0507/786 and from the Programme National de Télédétection Spatiale (PNTS). We thank Jawaharlal Nehru University for providing all the facilities to carry out this work. Prof. K. A. Brugger and another anonymous reviewer have provided constructive suggestions and comments which helped to significantly improve the manuscript. They are greatly acknowledged here.


Rolando Biscay

Hélène Lescornel

Jean-Michel Loubes

Adaptive Covariance Estimation with model selection

Keywords: covariance estimation, model selection, adaptive procedure

published or not. The documents may come     Adaptive Covariance Estimation with model selection.

Introduction

Estimating the covariance function of stochastic processes is a fundamental issue in statistics with many applications, ranging from geostatistics, financial series or epidemiology for instance (we refer to [START_REF] Stein | Interpolation of spatial data[END_REF], [START_REF] Journel | Kriging in terms of projections[END_REF] or [START_REF] Cressie | Statistics for spatial data[END_REF] for general references). While parametric methods have been extensively studied in the statistical literature (see [START_REF] Cressie | Statistics for spatial data[END_REF] for a review), nonparametric procedures have only recently received attention, see for instance [START_REF] Elogne | Non parametric estimation of smooth stationary covariance functions by interpolation methods[END_REF][START_REF] Bigot | Nonparametric estimation of covariance functions by model selection[END_REF][START_REF] Bigot | Adaptive estimation of spectral densities via wavelet thresholding and information projection[END_REF][START_REF] Bigot | Group lasso estimation of high-dimensional covariance matrices[END_REF] and references therein.

In [START_REF] Bigot | Nonparametric estimation of covariance functions by model selection[END_REF], a model selection procedure is proposed to construct a non parametric estimator of the covariance function of a stochastic process under mild assumptions. However their method heavily relies on a prior knowledge of the variance. In this paper, we extend this procedure and propose a fully data driven penalty which leads to select the best covariance among a collection of models. This result constitutes a generalization to the matricial regression model of the selection methodology provided in [START_REF] Baraud | Model selection for regression on a fixed design[END_REF].

Consider a stochastic process (X (t)) t∈T taking its values in R and indexed by T ⊂ R d , d ∈ N. We assume that E [X (t)] = 0 ∀t ∈ T and we aim at estimating its covariance function σ (s, t) = E [X (s) X (t)] < ∞ for all t, s ∈ T . We assume we observe X i (t j ) where i ∈ {1 . . . n} and j ∈ {1 . . . p}. Note that the observation points t j are fixed and that the X i 's are independent copies of the process X. Set x i = (X i (t 1 ) , . . . , X i (t p )) ∀i ∈ {1 . . . n} and denote by Σ the covariance matrix of X at the observations points Σ =E x i x i = (σ (t j , t k )) 1≤j≤p,1≤k≤p .

Following the methodology presented in [START_REF] Bigot | Nonparametric estimation of covariance functions by model selection[END_REF], we approximate the process X by its projection onto some finite dimensional model. For this, consider a countable set of functions (g λ ) λ∈Λ which may be for instance a basis of L 2 (T ) and choose a collection of models M ⊂ P (Λ). For m ⊂ M, a finite number of indices, the process can be approximated by

X (t) ≈ λ∈m a λ g λ (t) .
Such an approximation leads to an estimator which depends on the collection of functions m, denoted by Σm . Our objective is to select in a data driven way, the best model, i.e. the one close to an oracle m 0 defined as the minimizer of the quadratic risk, namely

m 0 ∈ arg min m∈M R (m) = arg min m∈M E Σ -Σm 2 .
This result is achieved using a model selection procedure.

The paper falls into the following parts. The description of the statistical framework of the matrix regression is given in Section 2. Section 3 is devoted to the main statistical results. Namely we recall the results of the estimate given in [START_REF] Bigot | Nonparametric estimation of covariance functions by model selection[END_REF] and prove an oracle inequality with a fully data driven penalty. Section 4 states technical results which are used in all the paper, while the proofs are postponed to the Appendix.

Statistical model and notations

We consider an R-valued process X (t) indexed by T a subset of R d with expectation equal to 0. We are interested in its covariance function denoted by σ (s

, t) = E [X (s) X (t)].
We have at hand the observations x i = (X i (t 1 ) , . . . , X i (t p )) for 1 i n where X i are independent copies of the process and t j are deterministic points. We note Σ ∈ R p×p the covariance matrix of the vector x i .

Hence we observe

x i x i = Σ + U i , 1 i n (1)
where U i are i.i.d. error matrices with expectation 0. We denote by S the empirical covariance of the sample : S = 1 n n i=1 x i x i . We use the Frobenius norm defined by A 2 = Tr AA for all matrix A. Recall that for a given matrix A ∈ R p×q , vec(A) is the vector in ∈ R pq obtained by stacking the columns of A on top of one another. We denote by A -the reflexive generalized inverse of the matrix A, see for instance in [START_REF] Seber | A matrix handbook for statisticians[END_REF] or [START_REF] Engl | Regularization of inverse problems[END_REF].

The idea is to consider that we have a quite good approximation of the process in the following form

X (t) ≈ λ∈m a λ g λ (t) , (2)
where m is a finite subset of a countable set Λ , (a λ ) λ∈Λ are random coefficients in R and (g λ ) λ∈Λ are real valued functions. We will consider models m among a finite collection denoted by M .

We note G m ∈ R p×|m| where (G m ) jλ = g λ (t j ) and a m the random vector of R |m| with coefficients (a λ ) λ∈m .

Hence, we obtain the following approximations :

x = (X (t 1 ) , .., X (t p )) ≈ G m a m xx ≈ G m a m a m G m Σ ≈ G m E a m a m G m
Thus, this point of view leads us to approximate Σ by a matrix in the subset

S (G m ) = G m ΨG m /Ψ symmetric in R |m|×|m| ⊂ R p×p . (3)
Hence, for a model m, a natural estimator for Σ is given by the projection of S onto S (G m ). We can prove using standard algebra (see in [START_REF] Bigot | Nonparametric estimation of covariance functions by model selection[END_REF] for a general proof) that it has the following form :

Σ m = Π m SΠ m m ∈ M ∈ R p×p , (4)
where

Π m = G m G m G m -G m ∈ R p×p (5)
are orthogonal projection matrices. Set

D m = T r (Π m ⊗ Π m )
which is the dimension of S (G m ) assumed to be positive, and Σ m = Π m ΣΠ m the projection of Σ onto this subspace.

Hence we obtain the model selection procedure defined in [START_REF] Bigot | Nonparametric estimation of covariance functions by model selection[END_REF]. The estimation error for a model m ∈ M is given by

E Σ -Σ m 2 = Σ -Π m ΣΠ m 2 + δ 2 m D m n , (6)
where

δ 2 m = Tr ((Π m ⊗ Π m ) Φ) D m , Φ=V vec x 1 x 1 .
Given θ > 0, it is thus natural to define the penalized covariance estimator Σ = Σ m by

m = arg min m∈M 1 n n i=1 x i x i -Σ m 2 + pen (m) , where pen (m) = (1 + θ) δ 2 m D m n . (7)
The following result proved in [START_REF] Bigot | Nonparametric estimation of covariance functions by model selection[END_REF] states an oracle inequality for the estimator Σ.

Theorem 2.1. Let q > 0 be given such that there exists κ > 2 (1 + q) satisfying E x 1 x 1 κ < ∞. Then, for some constants K (θ) > 1 and C (θ, κ, q) > 0 we have that

E Σ -Σ 2q 1/q ≤ 2 (q -1 -1) + K (θ) inf m∈M Σ -Π m ΣΠ m 2 + δ 2 m D m n + ∆ κ n δ 2 sup ,
where

∆ q κ = C (θ, κ, q) E x 1 x 1 κ m∈M δ -κ m D -(κ/2-1-q) m and δ 2 sup = max δ 2 m : m ∈ M .
However the penalty defined here depends on the quantity δ m which is unknown in practice since it relies on the matrix Φ = V vec xx . Our objective is to study a covariance estimator built with a new penalty involving an estimator of Φ.

More precisely, we will replace pen(m) by an empirical version pen(m), where

pen (m) = (1 + θ) δ 2 m D m n , (8)
and

δ 2 m = Tr (Π m ⊗ Π m ) Φ D m ,
with Φ an estimator of Φ.

The objective is to generalize Theorem 2.1 and to construct a fully adaptive penalized procedure to estimate the covariance function.

Main result : adaptive penalized covariance estimation

Here we state the oracle inequality obtained for the new covariance estimator introduced previously.

Set

y i = vec x i x i , 1 i n,
which are vectors in R p 2 and denote by S vec = 1 n n i=1 y i their empirical mean. Consider the following constant C inf = inf m∈M Tr ((Π m ⊗ Π m ) Φ), and assume that the collection of models is chosen such that C inf > 0. Set

Φ = 1 n n i=1 y i y i -S vec S vec , δ 2 m = Tr (Π m ⊗ Π m ) Φ D m .
Given θ > 0, we consider the covariance estimator Σ = Σ m with

m = arg min m∈M 1 n n i=1 x i x i -Σ m 2 + pen (m) ,
where

pen (m) = (1 + θ) δ 2 m D m n . (9)
Theorem 3.1. Let 1 q > 0 be given such that there exists

β > max (2 (1 + 2q) , 3 + 2q) satisfying E xx β < ∞.
Then, for a constant C depending on θ, β and q, we have for n n(β, θ, C inf , Σ), and ∀κ ∈ ]2 (1 + 2q) ; min (β, 2β -4)[ :

E Σ -Σ 2q 1/q C inf m∈M Σ -Σ m 2 + δ 2 m D m n (10)
+ C n ∆ β E xx β 2 β + Σ 2 + δ 2 sup ∆ κ ( 11
)
where

∆ q β = c (θ, β, q) E xx β m∈M δ -β m D -β/2 m 1-2q κ ∆ q κ = C (θ, κ, q) E xx κ m∈M δ -κ m D -(κ/2-1-q) m and δ 2 sup = max δ 2 m : m ∈ M .
We have obtained in Theorem 3.1 an oracle inequality since the estimator Σ has the same quadratic risk as the "oracle" estimator except for an additive term of order O 1 n and a constant factor. Hence, the selection procedure is optimal in the sense that it behaves as if the true model were at hand.

The proof of this theorem is divided into two parts. First, as in the of Theorem 2.1 proved in [START_REF] Bigot | Nonparametric estimation of covariance functions by model selection[END_REF], we will consider a vectorized version of the model [START_REF] Baraud | Model selection for regression on a fixed design[END_REF]. In this technical part we will obtain an oracle inequality under some particular assumptions for a general penalty. In a second part, we will prove that our particular penalty verifies these assumptions by using properties of the estimator Φ.

Technical results

Vectorized model

Here we consider the vectorized version of model [START_REF] Baraud | Model selection for regression on a fixed design[END_REF]. In this case, we observe the following vectors in R p 2 :

y i = f i + ε i 1 i n. ( 12
)
Here y i corresponds to vec x i x i in the model (1), f i to vect (Σ) and ε i to vec (U i ). We set f = f 1 , . . . , f n , y = y 1 , . . . , y n and ε = ε 1 , . . . , ε n , which are vectors in

R np 2 .
We estimate f by an estimator of the form

f m = P m y m ∈ M,
where P m is the orthogonal projection onto a subspace S m of dimension D m . We note f m = P m f and we consider the empirical norm

f 2 n = 1 n n i=1 f i f i with the corresponding scalar product •, • n .
First we state the vectorized form of Theorem 2.1. Write

δ 2 m = Tr (P m (I n ⊗ Φ)) D m , δ 2 sup = max δ 2 m : m ∈ M .
Given θ > 0, define the penalized estimator f = f m , where

m = arg min m∈M y-f m 2 n + pen (m) , with pen (m) = (1 + θ) δ 2 m D m n .
Then, the proof of Theorem 2.1 relies on the following proposition proved in [START_REF] Bigot | Nonparametric estimation of covariance functions by model selection[END_REF]:

Proposition 4.1. : Let q > 0 be given such that there exists κ > 2 (1 + q) satisfying E ε 1 κ < ∞. Then, for some constants K (θ) > 1 and C (θ, κ, q) > 0 we have that

E f -f 2q n 1/q ≤ 2 (q -1 -1) + K (θ) inf m∈M f -P m f 2 n + δ 2 m D m n + ∆ κ n δ 2 sup , (13)
where

∆ q κ = C (θ, κ, q) E ε 1 κ m∈M δ -κ m D -(κ/2-1-q) m .
The new estimator Σ defined previously corresponds here to the estimator f = f m , where

m = arg min m∈M y-f m 2 n + pen (m) , with pen (m) = (1 + θ) δ 2 m D m n ,
and δ 2 m is some estimator of δ 2 m . Next Proposition gives an oracle inequality for this estimator under new assumptions on the model. As Proposition 4.1, it is inspired by the paper [START_REF] Baraud | Model selection for regression on a fixed design[END_REF]. Proposition 4.2. Let 1 q > 0 be given such that there exists κ

> 2 (1 + 2q) satisfying E ε 1 κ < ∞. For α ∈ ]0; 1[, set Ω = ∩ m∈M δ 2 m (1 -α) δ 2 m . Assume that A1. E δ 2 m δ 2 m .
A2. P (Ω c ) C (α) 1 n γ for some γ q 1-2q/κ . Then, for a constant C depending on κ, θ and q, and we have

E f -f 2q n 1/q C inf m∈M f -P m f 2 n + δ 2 m D m n (14) + C n ∆ κ E [ ε 1 κ ] 2 κ + f 2 n + δ 2 sup ∆ κ ( 15
)
where

∆ q κ = C (α) (1-2q κ ) with α = α (θ) is fixed in ]0; 1[ and
∆ q κ = C (θ, κ, q) E ε 1 κ m∈M δ -κ m D -(κ/2-1-q) m
Theorem 3.1 is thus a direct application of Proposition 4.2. Hence only remain to be checked the two assumptions A1 and A2.

Auxiliary concentration type lemmas

Here we state some propositions required in the proofs of the previous results.

To our knowledge, the first is due to von Bahr and Esseen in [START_REF] Bahr | Inequalities for the rth absolute moment of a sum of random variables, 1 ≤ r ≤ 2[END_REF].

Lemma 4.3. Let U 1 , . . . , U n independent centred variables with values in R. For any 1 κ 2 we have :

E n i=1 U i κ 8 n i=1 E [|U i | κ ]
The next proposition is proved in [START_REF] Bigot | Nonparametric estimation of covariance functions by model selection[END_REF]. . For all β ≥ 2 such that E ε 1 β < ∞ it holds that, for all x > 0,

P ζ 2 (ε) ≥ δ 2 Tr A + 2δ 2 Tr A ρ A x + δ 2 ρ A x ≤ C 2 (β) E ε 1 β Tr A δ β ρ A x β/2
, (16) where the constant C 2 (β) depends only on β.

Appendix

Proof of Proposition 4.2

This proof follows the guidelines of the proof of Theorem 6.1 in [START_REF] Baraud | Model selection for regression on a fixed design[END_REF]. The following lemma will be helpful for the proof of this proposition

Lemma 5.1. Choose η = η (θ) > 0 and α = α (θ) ∈ ]0; 1[ such that (1 + θ) (1 -α) (1 + 2η). Set H m (f ) = f -f 2 n -κ (θ) f -f m 2 n + Dm n δ 2 m + where κ (θ) = 2 + 4 η (1 + θ). Then, for m 0 minimizing m → f -f m 2 n + Dm n δ 2 m in m ∈ M E [H m 0 (f ) q 1 Ω ] ∆ q κ δ 2q sup 1 n q . ( 17
)
where ∆ κ was defined in Proposition 4.2.

Proof. Lemma 5.1 First, remark that on the set Ω, for all m ∈ M

pen (m) (1 -α) (1 + θ) δ 2 m D m n (1 + 2η) δ 2 m D m n . Set pen(m) = (1 + 2η) δ 2 m Dm
n , which corresponds to the penalty of Proposition 4.1. The proof of this lemma is based on the proof of Proposition 4.1 in [START_REF] Bigot | Nonparametric estimation of covariance functions by model selection[END_REF]. In fact, it is sufficient to prove that for each x > 0 and κ ≥ 2

P H (f ) 1 Ω ≥ 1 + 2 η x n δ 2 m ≤ c (κ, η) E ε 1 κ m∈M 1 δ κ m D m ∨ 1 (ηD m + x) κ/2 , ( 18
)
where we have set

H (f ) = f -f 2 n -2 + 4 η f -f m 0 2 n + pen (m 0 ) + . Indeed, for each m ∈ M, f -f m 0 2 n + pen (m 0 ) = f -f m 0 2 n + (1 + θ) δ 2 m 0 n D m 0 ≤ (1 + θ) f -f m 0 2 n + δ 2 m 0 n D m 0
then we get that for all q > 0,

H q (f ) 1 Ω ≥ H q m 0 (f ) 1 Ω (19)
Using the equality

E [H q (f ) 1 Ω ] = ∞ 0 qu q-1 P (H q (f ) 1 Ω > u) du
and following the proof of Propositon 4.1 in [START_REF] Bigot | Nonparametric estimation of covariance functions by model selection[END_REF] we obtain the upper bound (17) of Lemma 5.1. Now we turn to the proof of (18). For any g ∈ R np 2 we define the empirical quadratic loss function by γ n (g) = y -g 2 n . Using the definition of γ n we have that for all g ∈ R np 2 ,

f -g 2 n = γ n (g) + 2 g -y, ε n + ε 2 n and therefore f -f 2 n -f -P m 0 f 2 n = γ n f -γ n (P m 0 f ) + 2 f -P m 0 f, ε n . ( 20
)
Using the definition of f , we know that

γ n f + pen ( m) ≤ γ n (g) + pen (m 0 )
for all g ∈ S m 0 . Then

γ n f -γ n (P m 0 f ) ≤ pen (m 0 ) -pen ( m) . (21)
So we get from (20) and (21) that

f -f 2 n ≤ f -P m 0 f 2 n + pen (m 0 ) -pen ( m) + 2 f -P m 0 f, ε n + 2 P m f -f, ε n + 2 f -P m f, ε n . (22)
In the following we set for each m ∈ M,

B m = {g ∈ S m : g n ≤ 1} , G m = sup t∈B m g, ε n = P m ε n , u m = P m f -f P m f -f n if P m f -f n = 0 0 otherwise. Since f = P m f + P m ε, (22) gives f -f 2 n ≤ f -P m 0 f 2 n + pen (m 0 ) -pen ( m) +2 f -P m 0 f n | u m 0 , ε n | + 2 f -P m f n | u m , ε n | + 2G 2 m . (23)
Using repeatedly the following elementary inequality that holds for all positive numbers ν, x, z

2xz ≤ νx 2 + 1 ν z 2 (24)
we get for any m ∈ M

2 f -P m f n | u m , ε n | ≤ ν f -P m f 2 n + 1 ν | u m , ε n | 2 . ( 25
)
By Pythagora's Theorem we have

f -f 2 n = f -P m f 2 n + P m f -f 2 n = f -P m f 2 n + G 2 m . (26)
We derive from ( 23) and ( 25) that for any ν > 0

f -f 2 n ≤ f -P m 0 f 2 n + ν f -P m 0 f 2 n + 1 ν u m 0 , ε 2 n +ν f -P m f 2 n + 1 ν u m , ε 2 n + 2G 2 m + pen (m 0 ) -pen ( m) .
Now taking into account that by equation ( 26)

f -P m f 2 n = f -f 2 n -G 2 m the above inequality is equivalent to (1 -ν) f -f 2 n ≤ (1 + ν) f -P m 0 f 2 n + 1 ν u m 0 , ε 2 n + 1 ν u m , ε 2 n + (2 -ν) G 2 m + pen (m 0 ) -pen ( m) . (27)
We choose ν = 2 2+η ∈ ]0, 1[, but for sake of simplicity we keep using the notation ν. Let p 1 and p 2 be two functions depending on ν mapping M into R + . They will be specified as in [START_REF] Bigot | Nonparametric estimation of covariance functions by model selection[END_REF] to satisfy

pen (m ) ≥ (2 -ν) p 1 (m ) + 1 ν p 2 (m ) ∀(m ) ∈ M . ( 28
)
Remember that on Ω, pen(m) pen(m) ∀m ∈ M. Since 1 ν p 2 (m ) ≤ pen (m ) and 1 + ν ≤ 2, we get from ( 27) and (28) that on the set

Ω (1 -ν) f -f 2 n ≤ (1 + ν) f -P m 0 f 2 n + pen (m 0 ) + 1 ν p 2 (m 0 ) + (2 -ν) G 2 m -p 1 ( m) + 1 ν u m , ε 2 n -p 2 ( m) + 1 ν u m 0 , ε 2 n -p 2 (m 0 ) ≤ 2 f -P m 0 f 2 n + pen (m 0 ) + (2 -ν) G 2 m -p 1 ( m) + 1 ν u m , ε 2 n -p 2 ( m) + 1 ν u m 0 , ε 2 n -p 2 (m 0 ) . ( 29
)
As 2 1-ν = 2 + 4 η we obtain that

(1 -ν) H (f ) 1 Ω = (1 -ν) f -f 2 n -(1 -ν) 2 + 4 η f -P m 0 f 2 n + pen (m 0 ) + 1 Ω = (1 -ν) f -f 2 n -2 f -P m 0 f 2 n + pen (m 0 ) + 1 Ω ≤ (2 -ν) G 2 m -p 1 ( m) + 1 ν u m , ε 2 n -p 2 ( m) + 1 ν u m , ε 2 n -p 2 (m 0 ) + For any x > 0, P (1 -ν) H (f ) 1 Ω ≥ xδ 2 m n ≤ P ∃m ∈ M : (2 -ν) G 2 m -p 1 (m ) ≥ xδ 2 m 3n + P ∃m ∈ M : 1 ν u m , ε 2 n -p 2 (m ) ≥ xδ 2 m 3n ≤ m ∈M P (2 -ν) P m ε 2 n -p 1 (m ) ≥ xδ 2 m 3n + m ∈M P 1 ν u m , ε 2 n -p 2 (m ) ≥ xδ 2 m 3n := m ∈M P 1,m (x) + m ∈M P 2,m (x) . ( 30
)
From now on, the proof of Lemma 5.1 is exactly the same as the end of the proof of Proposition 4.1 in [START_REF] Bigot | Nonparametric estimation of covariance functions by model selection[END_REF] with L m = ν.

Proof. Proposition 4.2

We first provide an upper bound for E f -f 2q n 1 Ω , where the set Ω depends on α chose as in Lemma 5.1.

As q 1, we have (a + b) q a q + b q . Together with Lemma 5.1 we deduce that

E f -f 2q n 1 Ω ∆ q κ δ 2q sup 1 n q + E κ (θ) q f -f m 0 2 n + D m 0 n δ 2 m 0 q .
Using the convexity of x → x 1 q together with the Jensen inequality, we obtain

E f -f 2q n 1 Ω 1 q 2 1/q-1 ∆ κ δ 2 sup 1 n + 2 1/q-1 E κ (θ) f -f m 0 2 n + D m 0 n δ 2 m 0 ,
and by using the assumption A1 we have that

E f -f 2q n 1 Ω 1 q 2 1/q-1 ∆ κ δ 2 sup 1 n + 2 1/q-1 κ (θ) f -f m 0 2 n + D m 0 n δ 2 m 0 . (31)
Now we need to find an upper bound for the quantity E f -

f 2q n 1 Ω c . First, remark that f -f 2 n = f -P m y 2 n = f -P m f 2 n + P m (f -y) 2 n f -P m f 2 n + ε 2 n = f 2 n -P m f 2 n + ε 2 n And thus f -f 2 n f 2 n + ε 2 n .
So we have

E f -f 2q n 1 Ω c f 2q n P (Ω c ) + E ε 2q n 1 Ω c .
Using Hölder's inequality with κ 2q > 1 we obtain

E ε 2q n 1 Ω c E [ ε κ n ] 2q κ P (Ω c ) (1-2q κ ) . But E [ ε κ n ] = 1 n κ 2 E   n i=1 ε i 2 κ 2   ,
and as κ 2, we can use Minkowsky's inequality to obtain

E [ ε κ n ] 1 n κ 2 n i=1 (E [ ε i κ ]) 2 κ κ 2 = 1 n κ 2 n (E [ ε 1 κ ]) 2 κ κ 2 ,
that is

E [ ε κ n ] E [ ε 1 κ ] .
So we have

E f -f 2q n 1 Ω c E [ ε 1 κ ] 2q κ + f 2q n P (Ω c ) (1-2q κ ) ,
and with assumption A2

E f -f 2q n 1 Ω c E [ ε 1 κ ] 2q κ + f 2q n C(α) 1 n γ (1-2q κ ) .
As γ q 1-2q/κ , we deduce that

E f -f 2q n 1 Ω c 1 q 2 1 q -1 E [ ε 1 κ ] 2 κ + f 2 n C(α) 1- 2q κ q 1 n (32)
To conclude, we use again the convexity of x → x 1 q and the inequality (31) to get

E f -f 2q n 1 q 4 1 q -1 E [ ε 1 κ ] 2 κ + f 2 n C(α) 1- 2q κ q 1 n +4 1/q-1 ∆ κ δ 2 sup 1 n +4 1/q-1 κ (θ) f -f m 0 2 n + D m 0 n δ 2 m 0

Proof of Theorem 3.1

Recall that β > max (2 (1 + 2q) , 3 + 2q) and κ ∈ ]2 (1 + 2q) ; min (β, 2β -4)[. In order to use Proposition 4.2 , we need to prove the following inequalities :

A1. E δ 2 m δ 2 m A2. P (Ω c ) C (α) 1 n γ for γ q 1-2q/κ . First we prove A1. Remember that δ 2 m =
Tr((Πm⊗Πm) Φ) Dm

. By using the linearity of the trace and the equal-

ity E Φ = n-1 n Φ, we obtain that E δ 2 m = n-1 n δ 2
m which proves the result.

For the second, write

Ω c = ∪ m∈M δ 2 m (1 -α) δ 2 m .
We bound up the quantity

P δ 2 m (1 -α) δ
P δ 2 m (1 -α) δ 2 m 1 n γ C 2 (β)2 β+1 + C 1 (β) 1 α β 2 E xx β δ -β m D -β 2 m ,
for γ q 1-2q/κ . This Proposition concludes the proof of A2 with

C (α) = C 2 (β)2 β+1 + C 1 (β) 1 α β 2 E xx β m∈M δ -β m D -β 2 m .
Proof. Proposition 5.2

We start by dividing

P m = P δ 2 m (1 -α) δ 2 m
into two parts with one of them involving a sum of independent variables with expectation equal to 0.

P m = P Tr (Π m ⊗ Π m ) Φ (1 -α) Tr ((Π m ⊗ Π m ) Φ) P m = P Tr (Π m ⊗ Π m ) Φ -Φ + µµ + µµ -αTr ((Π m ⊗ Π m ) Φ) P m P Tr (Π m ⊗ Π m ) 1 n n i=1 y i y T i -Φ -µµ + µµ -S vec S vec αTr ((Π m ⊗ Π m ) Φ) P m P Tr (Π m ⊗ Π m ) 1 n n i=1 y i y T i -Φ -µµ α 2 Tr ((Π m ⊗ Π m ) Φ) +P Tr (Π m ⊗ Π m ) µµ -S vec S vec α 2 Tr ((Π m ⊗ Π m ) Φ) Set Q1 = P Tr (Π m ⊗ Π m ) 1 n n i=1 y i y T i -Φ -µµ α 2 Tr ((Π m ⊗ Π m ) Φ) and Q2 = P Tr (Π m ⊗ Π m ) µµ -S vec S vec α 2 Tr ((Π m ⊗ Π m ) Φ)
Study of Q1 First we use Markov's inequality to obtain

A1 2 β 2 E Tr (Π m ⊗ Π m ) 1 n n i=1 y i y T i -Φ -µµ β 2 (αTr ((Π m ⊗ Π m ) Φ)) β 2
.

We must consider the two following cases :

• If β 2 2, Rosenthal's inequality gives E   1 n n i=1 Tr (Π m ⊗ Π m ) y i y T i -Φ -µµ β 2   C β 2 1 n β 2 -1 E Tr (Π m ⊗ Π m ) y 1 y 1 -Φ -µµ β 2 +C β 2 1 n E Tr (Π m ⊗ Π m ) y 1 y 1 -Φ -µµ 2 β 4
.

As β 2 2, 1 n β 2 -1 1 n β 4
and we can use Jensen's inequality on the second term to obtain

E   1 n n i=1 Tr (Π m ⊗ Π m ) y i y T i -Φ -µµ β 2   C β 2 E Tr (Π m ⊗ Π m ) y 1 y 1 -Φ -µµ β 2 2 n β 4 . • If 1 β 2 2, we use Lemma 4.3 of subsection 4.2 to get E   1 n n i=1 Tr (Π m ⊗ Π m ) y i y T i -Φ -µµ β 2   8 n β 2 -1 E Tr (Π m ⊗ Π m ) y 1 y 1 -Φ -µµ β 2
In both cases, we can use the fact that x → x β 2 is a convex and increasing function to obtain

E Tr (Π m ⊗ Π m ) y 1 y 1 -Φ -µµ β 2 2 β 2 -1 E Tr (Π m ⊗ Π m ) y 1 y 1 β 2 + Tr (Π m ⊗ Π m ) Φ + µµ β 2 .
And by using the Jensen's inequality on the second term we have that

E Tr (Π m ⊗ Π m ) y 1 y 1 -Φ -µµ β 2 2 β 2 E Tr (Π m ⊗ Π m ) y 1 y 1 β 2
. Now consider the following lemma.

Lemma 5.3. If Ψ is symmetric non-negative definite, then Tr ((Π m ⊗ Π m ) Ψ) ∈ [0; Tr (Ψ)] (33)
From this fact we get that

Tr (Π m ⊗ Π m ) y 1 y 1 β 2 Tr y 1 y 1 β 2 = y 1 β = xx β .
In conclusion, we have

Q1 C 1 (β) E xx β α β 2 δ β m D β 2 m 1 n γ , ( 34
) with γ = min β 4 , β 2 -1 and C 1 (β) = 2C β 2 if β 4 where C β 2 is the constant in Rosenthal's inequality and C 1 (β) = 8 if 2 β 4. Remark that β 4 κ 4 and β 2 -1 κ 4 , so γ κ 4 .

Study of Q2

Recall that

Q2 = P Tr (Π m ⊗ Π m ) µµ -S vec S vec α 2 Tr ((Π m ⊗ Π m ) Φ) . Set B2 = Tr (Π m ⊗ Π m ) µµ -S vec S vec .
Using the properties of the trace, we can write

B2 = Tr (Π m ⊗ Π m ) µµ -Tr (Π m ⊗ Π m ) S vec S vec = Tr µ (Π m ⊗ Π m ) µ -Tr S vec (Π m ⊗ Π m ) S vec .
But Π m ⊗ Π m is an orthogonal projection matrix, then

In conclusion, with

A m = H m G n = 1 n      Π m ⊗ Π m . . . Π m ⊗ Π m Π m ⊗ Π m . . . Π m ⊗ Π m . . . . . . . . . Π m ⊗ Π m . . . Π m ⊗ Π m      ∈ R p 2 n×p 2 n , we have that A m (y -f ) = 1 n ⊗ ((Π m ⊗ Π m ) (S vec -µ)) .
Moreover, A m is an orthogonal projection matrix and we have the following equalities

A m (y -f ) 2 = n (Π m ⊗ Π m ) (S vec -µ) 2 = (y -f ) A m (y -f ) , Tr (A m ) = n n Tr (Π m ⊗ Π m ) = D m , Tr (A m (Id n ⊗ Φ)) = n n Tr ((Π m ⊗ Π m ) Φ) = Tr ((Π m ⊗ Π m ) Φ) .

Now we can use Proposition 4.4 with

A = A m , ε i = y i -µ, Tr (A m ) = D m , ρ (A m ) = 1, δ 2 = δ 2
m and β 2. This gives for all x > 0 .

P (y -f ) A m (y -f ) Tr ((Π m ⊗ Π m ) Φ) 1 + x Dm 2 C 2 (β) E[ y 1 -µ β ]D
Moreover E y 1 -µ β E ( y 1 + µ ) β ,
and by using properties of convexity we obtain

E y 1 -µ β 2 β-1 E y 1 β + µ β .
With the Jensen's inequality we get:

E y 1 -µ β 2 β E y 1 β .
In conclusion, with r =

Proposition 4 . 4 .

 44 Given N, k ∈ N, let A ∈ R N k×N k {0} be a non-negative definite and symmetric matrix and ε 1 , ..., ε N i.i.d random vectors in R k with E (ε 1 ) = 0 andV (ε 1 ) = Φ. Write ε = ε 1 , ..., ε N , ζ (ε) = ε Aε, and δ 2 = Tr( A(I N ⊗Φ)) Tr( A)

P 1 n 2 C 2 2 . 8 µn γ with γ q 1 - 2

 1222812 (Π m ⊗ Π m ) (S vec -µ) 2 Tr ((Π m ⊗ Π m ) Φ) 1 + x Dm (β) E[ y 1 -µ β ]DIn order to use this deviation bound to obtain the inequalitiesP (Π m ⊗ Π m ) (µ -S vec ) 2 α 4 Tr ((Π m ⊗ Π m ) Φ) C 1 n γ and P (Π m ⊗ Π m ) (µ -S vec ) α Tr ((Π m ⊗ Π m ) Φ) C 1 2q/κ ,we need to find x > 0 satisfying the three following facts Tr ((Π m ⊗ Π m ) Φ) and (37) hold for the choice x = D m n r with r < 1 and if n is large enough to have 1 obtain (38) with x = D m n r , we use the inequality D m n which gives

κ 2 +2β < 1 C 2 δ 2 m( 1 -α) δ 2 m 1 n 2 m 4 > 2+2q 4 q as q 1 .

 2122112441 we obtain for n n(κ, β, α, C inf , Σ)Q2 2 β+1 C 2 (β) (β) is the constant which appears in Proposition 4.4.In conclusion, combining (34) and (41)P κ/4 C 2 (β)2 β+1 + C 1 (β) for n n(κ, β, α, C inf , Σ) .To conclude, remark that κ 4 (1 -2q/κ) = κ-2qProof. Lemma 5.3Recall that Π m ⊗Π m is an orthogonal projection matrix. Hence there exists an orthogonal matrix P m such that P m (Π m ⊗ Π m ) P m = D, with D a diagonal matrix with D ii = 1 if i D m , and D ii = 0 otherwise. Then if Ψ is symmetric non-negative definite we have :Tr ((Π m ⊗ Π m ) Ψ) = Tr DP m ΨP m

  2 m in the following Proposition. Proposition 5.2. For all m ∈ M , α ∈]0; 1[ and n n(κ, β, α, C inf , Σ) we have for some constants C 1 (β), C 2 (β) :

Now we need to provide an upper bound for the quantities

For this we will use the deviation bound provided by Proposition 4.4 stated in subsection 4.2 .

Set

Now, if

we have Indeed, P m ΨP m is non-negative definite so all its diagonal entries are non-negative.