Shuman Jia

Loïc Cadour

Hubert Cochet

Maxime Sermesant

STACOM-SLAWT Challenge: Left Atrial Wall Segmentation and Thickness Measurement Using Region Growing and Marker-Controlled Geodesic Active Contour

Keywords: atrial fibrillation, left atrial wall thickness, 3-dimensional image segmentation, cardiac computed tomography (CT), region growing, geodesic active contour

Analyzing the structure of the left atrium can provide precious insights into the pathology of atrial fibrillation, eventually resulting in optimization of treatment plans. In this paper, an interactive and patient-specific method is presented to segment the left atrial endocardium 3 , the left atrial epicardium and measure the left atrial wall thickness from cardiac computed tomography images. A region growing algorithm was adapted to segment the left atrial endocardium, whereas the left atrial epicardium was segmented indirectly: a marker-controlled geodesic active contour model was defined on its surrounding environment. The results of the left atrial wall thickness were then mapped onto meshes generated from the endocardium segmentation. We tested our pipeline on 10 datasets as a part of the STACOM 2016 Left Atrial Wall Segmentation Challenge and we compared our method with manual segmentation. Aimed at facilitating the segmentation of the left atrial thin-wall structure, this pipeline is partially implemented in MUSIC software for clinical use. The expertise of clinicians can be added through the choice of specific parameters for each patient, although this remains optional owing to the robustness of the approach.

Introduction

Atrial fibrillation (AF) is the most common type of cardiac arrhythmia, characterized by uncoordinated electrical activation and disorganized contraction of the atria. Around 2% to 3% of the population in Europe and North America, as of 2014, were affected, and its prevalence rate is increasing worldwide [START_REF] Zoni-Berisso | Epidemiology of atrial fibrillation: European perspective[END_REF]. This epidemic, with or without symptoms, is likely to be associated with life-threatening consequences, including heart failure, heart attack and stroke.

Atrial ablation, an effective treatment for AF, may be recommended for drug refractory patients. This invasive procedure establishes transmural lesions to block the arrhythmia while avoiding impairing extra-cardiac tissues. Thereby the radio-frequency power dose, delivered for ablation during the procedure, relies on local myocardial thickness, which demonstrates variations by region and by subject [START_REF] Siew | Anatomy of the left atrium: implications for radiofrequency ablation of atrial fibrillation[END_REF][START_REF] Sánchez-Quintana | Anatomic relations between the esophagus and left atrium and relevance for ablation of atrial fibrillation[END_REF]. Reduction of left atrial (LA) wall thickness also appears in AF patients when compared with controls, according to study by autopsy [START_REF] Hassink | Morphology of atrial myocardium in human pulmonary veins: a postmortem analysis in patients with and without atrial fibrillation[END_REF].

Model-based segmentation and region growing approach have been tested on the segmentation of the left atrium [START_REF] Tobon-Gomez | Left atrial segmentation challenge: a unified benchmarking framework[END_REF][START_REF] Kutra | Automatic multi-modelbased segmentation of the left atrium in cardiac mri scans[END_REF]. Furthermore, several methods have been developed to measure automatically the LA wall thickness from computed tomography (CT) images. In [START_REF] Thomas A Dewland | Use of computed tomography to identify atrial fibrillation associated differences in left atrial wall thickness and density[END_REF], a segmentation of the wall on four regions, inter-atrial septum, below right inferior pulmonary vein, appendage and anterior wall was performed whereas [START_REF] Inoue | Left atrial wall segmentation from ct for radiofrequency catheter ablation planning[END_REF][START_REF] Inoue | Patient-specific left atrial wall-thickness measurement and visualization for radiofrequency ablation[END_REF] chose to build a pipeline based on multiregion segmentation method. However, this task remains challenging, because the LA wall is heterogeneous, possibly consisting of fat granules and fibrosis, especially in patients of myocardial diseases or persistent AF. The accuracy of the segmentation is also difficult to validate, as few manual segmentations or reliable ground truths are available.

In this paper, we present an interactive pipeline to segment the LA wall and measure wall thickness from cardiac CT images. The expertise of clinicians can be introduced through the choice of parameters. Compared to previous works, we propose an inverse way to segment the LA thin-wall structure. The segmentation was divided into 2 parts, LA endocardium and LA epicardium. The former was segmented with high accuracy using region growing combined with patientspecific intensity value threshold. The latter was segmented indirectly from its surrounding environment, to address the fuzzy boundaries problem. Algorithms employed were in 3-dimensional (3D) space.

The rest of the paper is organized as follows: we summarize our method using a flowchart in section 2, along with detailed illustration of region growing and marker-controlled geodesic active contour algorithms. Results and evaluation are shown in section 3, and a discussion of the limitations of this study in section 4.

Materials and Methods

Data Acquisition and Pre-processing

The method was applied to a database consisting of 10 3D CT image datasets, provided by the STACOM Left Atrial Wall Thickness Challenge.

Coronary Computed Tomography Angiography (CCTA) was performed on a Philips 256 iCT scanner. All patients were injected with an intravenous contrast agent. The scans were ECG-gated in a single breath hold. The images were acquired at 0.5 mm in-plane resolution with a slice thickness of 1 mm, and reconstructed to a 0.8 to 1.0 mm slice thickness with a 0.4 mm slice increment and a 250 mm field of view. The image matrix was kept at 512 × 512 matrix, constructed with a sharp reconstruction kernel [START_REF]STACOM Left Atrial Wall Thickness Challenge[END_REF]. The pixel values were in Hounsfield units (HU).

The images were manually cropped around the region of interest, namely the left atrium. The intensity values outside the range from -500 HU to 500 HU were set to nearer bound, so as to obtain a better visualization of blood pool, muscles and fat within the region that was meant to be segmented.

Methodology

We summarize our method into following steps, as shows the flowchart in Fig. 1: 1. Drawing polygons on axial slices of CT images to isolate the left atrium4 . 2. Sampling intensity values of blood and muscles.

3. Segmenting the LA endocardium using region growing according to intensity values statistics.

4. Thresholding neighboring tissues (fat, other blood pools) of the LA wall. 5. Segmenting the LA epicardium using geodesic active contour based on markers of neighboring tissues.

6. Calculating distance map of the LA endocardium segmentation. 7. Mapping the distance onto meshes of the LA endocardium.

Pre-Processing Segmentation Measurement

Crop and Display

Step 1

Step 2 and 3

Region Growing Endocardium Segmentation

Excluding Markers Epicardium Segmentation

Step 5

Mesh Mapped with Wall Thickness

Step 6 and 7 Currently most steps are implemented in MUSIC software (multimodality software for specific imaging in cardiology, developed by IHU Liryc, University of Bordeaux and Inria Sophia Antipolis) [START_REF] Cochet | Cardiac arrythmias: multimodal assessment integrating body surface ecg mapping into cardiac imaging[END_REF].

Marker-Controlled Geodesic Active Contour

Region Growing

We chose one pixel inside the left atrium, a connected component, as seed. The seed was raised into a gowning region. Neighboring pixels next to the growing region whose intensity values fell within a certain range were considered part of the left atrium. The region evolved iteratively until no more pixel was assigned.

To determine whether a pixel belonged to the left atrium, statistical analysis was performed on intensity values of tissues samples. Given the histogram of muscles samples and blood samples, we computed the mean and standard deviation of intensity values, for muscles I(µ, σ 2 ) and for blood I ′ (µ ′ , σ ′2 ). Then we set a threshold to distinguish these two classes as

th = µσ ′ + µ ′ σ σ + σ ′ . ( 1
)
Neighbor pixels of the growing region with intensity values upper than the threshold th were assigned to the left atrium. Accordingly, spatial constraint and user-provided intensity information are combined in the algorithm.

To further negate the impact of noise in original CT images, alternating sequential filters (ASFs) were used to post-process the LA endocardium segmentation. This step is compulsory for CT images with low quality. Alternating sequential filters are composed of pairs of idempotent morphological filters with structuring elements of increasing sizes. They extract the geometrical characteristics and minimize the distortion of objects [START_REF] Jean | Image Analysis and Mathematical Morphology[END_REF][START_REF] Pei | An efficient class of alternating sequential filters in morphology[END_REF]. In our experiment, we used pairs of closing and opening with structuring elements of increasing radius from 0.25 mm, 0.50 mm to 0.75 mm sequentially.

Marker-Controlled Geodesic Active Contour (MCGAC)

Geodesic active contour algorithm [START_REF] Caselles | Geodesic active contours[END_REF], represents contour as zero crossing of a level-set function. An initial contour propagates to touch the shape boundaries, as the level-set function progresses in time until

∂u ∂t = 0, ( 2
)
where u is the level-set function under following conditions:

u| t=0 = u 0 (3a) ∂u ∂t = g(c + κ)|∇u| + ∇u∇g, (3b)
where u 0 is the initial level-set function; c is a constant to provide a steady velocity; κ is related to the curvature of the level-set function as κ = div( ∇u |∇u| );

g is an edge detector function of image, which is strictly decreasing towards the outside of the region to be segmented and has values near 0 at the boundary to stop the evolution of curves. In equation 3b, three terms affect the way that a level-set function evolves:

• Advection term: gc|∇u|, related to progression speed of curves.

• Curvature term: gκ|∇u|, related to curvature of curves.

• Propagation term: ∇u∇g, related to expansion of curves.

The algorithm is implemented in Insight Segmentation and Registration Toolkit (ITK), which uses three parameters (advection, curvature and propagation scaling, noted here as S a , S c and S p ) to assign weighting ratios to each term and adjust their influence in the evolution process of the level-set function. We chose manually S a = 10 × S c = S p = 5. The initial level-set function u 0 was computed as a distance map of the LA endocardium segmentation. The edge detector function g was simply set to 0 on excluded regions and 1 on the rest of the image.

As for excluded regions, neighboring tissues other than muscles outside the LA wall were marked as excluded for the segmentation. Since CT images display relative linear attenuation values of tissues, and correspond to Hounsfield units (HU) scale, they may be used to characterize tissues [START_REF] Huda | Review of Radiologic Physics[END_REF]. Based on approximate HU values for air (-1000), fat (-100 to -50), muscle (10 to 50) and blood (300 to 400), we set 0 HU as upper threshold for fat (and air) and th, as defined in equation 1, as lower threshold for blood. Besides, fat granules of size smaller than 1 mm 3 were regarded as inside the LA wall, thus not excluded. Margins were kept to avoid misleading markers resulting from noise or density variation. At last, all pixels further than 6 mm away from the left atrium were also excluded, according to studies of the LA wall thickness on excised hearts [START_REF] Siew | Anatomy of the left atrium: implications for radiofrequency ablation of atrial fibrillation[END_REF][START_REF] Sánchez-Quintana | Anatomic relations between the esophagus and left atrium and relevance for ablation of atrial fibrillation[END_REF].

The markers of neighboring tissues were introduced to address the fuzzy boundaries problem, as the intensity gradient between the LA wall and its neighboring tissues is small and noise may blur original CT images. Then we adapted geodesic active contour model based on excluding markers, because the LA wall may be attached to muscles belonging to other organs such as esophagus, aorta, lung etc., and hence cannot be segmented only based on intensity information. Geodesic active contour, on the other hand, involves curvature constraint of shape to solve this problem.

Results

Parameters

Parameter setting mentioned in Section 2 can be altered depending on cases, but in this study we followed the same rule for all 10 cases to test the robustness of the algorithm with respect to different inputs.

Intensity value thresholds th used during region growing process were shown in Table 1.

Left Atrial Wall Segmentation

Here we present preliminary results of the LA wall segmentation, and comparison with manual segmentation, as shown in Fig. 2. The proposed method closely approximated the wall thickness, as compared with manual segmentation. Differences lay in the segmentation of the anterior LA wall, which is not fully involved in the manual segmentation provided. The intersecting surfaces inside pulmonary veins and on mitral valve are presented in our segmentation with rather thinner wall, which should not have been included.

Remark. We segmented the LA appendage and pulmonary veins connected to the left atrium as well, which is different from the manual segmentation provided.

Wall Thickness Measurement

We computed the LA wall thickness using the nearest distance from each point on the LA endocardium to the LA epicardium. Meshes of the LA endocardium were generated using The Visualization Toolkit (VTK).

The LA wall thickness varies by subject as well as by region, as shown in Fig. 3. From the mean and range of the wall thickness, we can see that dataset #8 has a thicker LA wall compared with other datasets, whereas dataset #5 and #6 have a thinner LA wall. The results for dataset #8 are shown as an example. Alternating sequential filters may be used to post-process the segmentation results. The LA epicardium segmentation, using marker-controlled geodesic active contour, takes into account the heterogeneity of the LA wall, and therefore can have small granules excluded. Whether or not to smooth the contour of the LA epicardium, depends on how the LA wall is defined.

Conclusion

We proposed a new method to segment the LA wall, making use of patientspecific intensity value information and surrounding environment of the LA wall. Despite good match of wall thickness with manual segmentation, the accuracy of the segmentation is still hard to validate as few reliable ground truths are available.

Although the study reached its aims, there exist some limitations to be mentioned. The proposed MCGAC method, applying uniform parameters to all regions of the LA wall, cannot achieve a region-wise constraint on curvature, which could lead to incorrectly assigned pixels. Potential future work to improve the accuracy of the segmentation may include: further eliminating the influence of noise; region-wise segmentation of the LA wall, which takes into consideration different surrounding environments of the LA sub-regions; parameters testing. Besides, alternative methods may be suitable as well to define how to measure the wall thickness, such as using a symmetric minimum distance between the epicardium and the endocardium [START_REF] Varela | 3d high-resolution atrial wall thickness maps using black-blood psir[END_REF], measuring the thickness in the surface normal direction of the endocardium. Combination of different approaches may be able to ameliorate the results.

Fig. 1 .

 1 Fig. 1. The flowchart of our method.

Fig. 2 .

 2 Fig. 2. Axial slices of CT images, overlapped with manual segmentation of the left atrial wall in green, our segmentation in red, intersection of the two in chartreuse. (a) Dataset #2; (b) dataset #3; (c) dataset #5; (d) example slice zoomed in.

Fig. 3 .

 3 Fig. 3. (a) Mean and interquartile range of the left atrial wall thickness for 10 datasets tested, with maximum whisker length specified as 1.5 times the interquartile range; (b) left atrial wall thickness map for dataset #8, posterior view; (c) left atrial wall thickness map for dataset #8, anterior view.

Table 1 .

 1 Intensity value thresholds used for 10 datasets tested. Threshold th 224.3 HU 419.5 HU 485.7 HU 226.6 HU 131.6 HU Threshold th 136 HU 448.2 HU 358.4 HU 345.8 HU 267.1 HU

	Dataset	#1	#2	#3	#4	#5
	Dataset	#6	#7	#8	#9	#10

We refer to the segmentation of the region inside the left atrial endocardium as the segmentation of the left atrial endocardium, the same for epicardium.

For step 1, polygons are drawn manually on a dozen of slices for each case and then interpolated automatically on the rest. Images are shown in axial planes.

Acknowledgments. Part of the research was funded by the Agence Nationale de la Recherche (ANR)/ERA CoSysMed SysAFib and ANR MIGAT projects. The authors would like to thank Marc-Michel Rohé and Hervé Delingette for their constructive feedback and Inria dtk team (Nicolas Schnitzler and Thibaud Kloczko) for helping in implementing the pipeline.


Kévin Graaff

Vincent Viel

Benoit Carlier

Candide Lissak

Gilles Arnaud-Fassetta

Monique Fort

Malika Madelin

Gravel sediment tracing from small torrents to trunk channels using pit-tags method: A case study from the upper Guil catchment (Queyras, French Alps)

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

Gravel sediment tracing from small torrents to trunk channels using pit-tags method:

A case study from the upper Guil catchment (Queyras, French Alps)

Kévin Gra (1,3), Vincent Viel (1,3), Benoit Carlier (1,3), Candide Lissak (2,3), Gilles Arnaud-Fassetta (1,3), Monique Fort (1,3), Malika Madelin (1,3)

(1) University of Paris-Diderot (Paris 7), Sorbonne Paris Cité, UMR 8586 PRODIG, (2) University of Caen Normandie, LETG Caen Geophen, UMR 6554 LETG, (3) ANR 12 SENV-0004 SAMCO

Contacts :

-GRAFF Kévin -UMR 8586 PRODIG kevin.gra @etu.univ-paris-diderot.fr -VIEL Vincent -UMR 8586 PRODIG vincent.viel@univ-paris-diderot.fr EGU 2016 Session GM9.1/HS11.15/NH1.24 Mountain streams : processes, channel morphology and knickpoint dynamics (co-organized)

1 Background and scientific context

In mountainous areas, especially in large catchments with torrential tributaries, the production and sediment transport signi cantly increase ood impacts in the valley bottoms. The quanti cation and characterisation of sedimentary transfers are therefore major challenges to provide better ood risk management. As a part of SAMCO (ANR 12 SENV-0004 SAMCO) project, for mountain hazard assessment in a context of global changes, we tried to improve the knowledge of these hydromorphological systems at both spatial and temporal scales, by identifying sediment supply and sediment dynamics from torrential tributaries to the trunk channel. A sediment budget was used as a tool for quantifying erosion, transport and deposition processes. This research is focused on a part of the sedimentary budget to understand the torrential processes and to de ne debris-ow's contribution. It is localized on the Upper Guil catchment, prone to catastrophic ash oods in active channel and alluvial fans. It is characterised by huge sediment transport coming from small torrential tributaries. In this study, we intend to highlight sediment dynamic between small torrential channels and its connection with gravel-bed streams, and characterize the sediment volume mobilized during a large debris-ow event.

2 Study area 0 200 m

N

The upper Guil catchment (Queyras -317 km²) is characterised by a lithology composed of 90-95% of shale (Fig. 2-a), supplying the river bed in huge ne and coarse sediment. Functionning in high magnitude/low frequency (Fort et al., 2002), the Guil River is subjected to ash-ood events. The volumes of sediment mobilized during the 100-year ood in1957 have been estimed at 330.000 m 3 (Koulinski, 2002). Sediment load in the main channel of the Guil River results from erosion in small torrential tributaries, as the Peynin catchment (Fig. 2-b). Operating in debris ows, the volumes deposited during June 2000 (RI-30 year ood event) have been estimed at 12.000 m 3 aggraded on the Peyronnelle fan (Fig. 2-c) (Einhorn, 2003).

) Estimation of debris ow volumes

Grain-size representativity of the Peyronnelle has been tested in order to estimate the size of the bedload transported (Gra , 2015). Pit-tags (p.t.) method has been used to estimate spatial dynamics in the Peyronnelle catchment during a debris ow event ( g.3-a). In the Peyronnelle stream, 320 p.t. were installed in 10 injection sites along the thalweg (Fig. 3-b). Each monitored particle was measured (3 axes-Long/Intermediate/Short) (Fig. 3-a). Pit tags have been spreading around 30% in the thalweg, and 70% in left and right banks. Five injection forms were determined and recorded to trace source of receivers displacement.

-A Di erential Global Positionning System (DGPS) has been used to locate and map each monitored particules. Precision is estimated from 0,5 m for horizontal dimension to 1,5 m for vertical dimension.

-A meteorological station was implemented in 2014 on the right side of the Peynin (Fig. 2-b). Rainfall were estimed every hour then been computed for everyday.

Three equations have been used to estimate volumes mobilized during the debris ow of August 2015. The rst equation used (1) is estimated from eld mapping and observations on alluvial fan. Equations (2) and (3) are theoretical estimations of the mean and maximum volume mobilisable by the torrential catchment. From these equations we will estimate the volume mobilized in the Peynin during a debris ow.

Where V f is the volume dropped on the Peyronnelle fan (in m 3 ), S is the area of debris ow (in m 2 ), H is the height of debris ow (in m).

Where V m is the average of debris ow mobilisable (in m 3 ), A d is surface of catchment (in m 2 ), S f the slope of the fan (in %), IG is a geological index (in Remaître, 2006).

Where V max is the maximum debris ow mobilisable (in m 3 ), L c is the lenght of the torrent (in m) (in Remaître, 2006). Volumes mobilized in the Peyronnelle is estimed at 24.250 m 3 (margin of error: 17%). Volumes mobilized on the Peynin have been estimed at 6.000 m 3 (margin of error: 25%). Sediment volumes mobilized during a debris is signi ciant (+30.000 m 3 ) but its contribution in the Peynin thalweg is only 20% of total sedimentary volume.

b) Trajectories and grain-size pattern c) Mobilized volumes of debris ow

During the debris ow event (Fig. 4-b), 101 monitored particules have been funded (e.g. 31%). They were gathered in 5 di erent groups : 35% of the total p.t. found in the Peynin, 16% in the downstream of the Peyronnelle, 12% in old debris ow channel, 25% in the main channel of the Peyronnelle and 12% in apex of alluvial fan (Fig. 4-d).

Average distances covered by monitored sediments are in between 58 m and 384 m (Fig. 4-d). The maximum distance observed is 720 m. Only 9% of the total p.t. found have not moved. The grain size of the found Pit Tags, and the average geometry, were three times higher than size set by Muir method (Fig. 4-c). Overnight on August 8 and 9, 2015, an important rainfall event (23 mm) has been recorded at the weather station in the Peynin catchment, due to intense storm cells (Fig. 4-a). Debris particules have been found in all sub-catchment, located on the left bank of the Peynin River (Four à Chaux, Trois Arbres, Peyronnelle, Roche Rousse -Fig. 2-b).

Rainfall recorded in the upper Guil was not su cient to switch to ash-ood regime in the Guil and tributaries. These data have to be used to discuss lag time, processes and thresholds. These three elements are needed to observe significant sediments fluxes. Sediment supply has also to be precised as it is a critical importance in the triggering of debris-flow events.

Figure 2

 2 Fig.2-a : Inset map of France with location of Guil River and geological structure (Lissak et al., 2015)

  Fig.3-b : Injections site of Pit tags

Figure 4 -

 4 Figure 4-d : Path covered by Pit tags (101 pt found)

Figure

  Figure 4-c : Grain size box plot

  Figure 4-a : Meteorological data of the Peynin in July -August

  Figure 3-a : Pit tag TIRIS 134 Hz -Monitored particule

Acknowledgements -This research was funded through the project SAMCO : "Adaptation of the societies to mountain risks in a context of global changes" (2013-2017) of the Agence Nationale de la Recherche (ANR-12-SENV-0004).


David R Lageson
email: lageson@montana.edu

Monique Fort

Roshan Bhattarai

Raj

Mary Hubbard

DAMAGE FROM THE APRIL-MAY 2015 GORKHA EARTHQUAKE SEQUENCE IN THE SOLUKHUMBU DISTRICT (EVEREST REGION), NEPAL

ABSTRACT: Rapid assessments of landslides (including other mass movements of rock, snow and ice) as well as human impacts were conducted by many organizations immediately following the 25 April 2015 M7.8 Gorkha earthquake and its aftershock sequence. In particular, a NASA-led satellite mapping campaign identified over 4,300 coseismic and post-seismic landslides throughout Nepal, largely concentrated in the "down-dropped northern tectonic block" well north of Kathmandu in the Greater Himalaya (Kargel et al., 2016). The Solukhumbu District lies at the eastern margin of this zone, near the epicenter of the 12 May 2015 M7.3 aftershock. Given the enormous cultural and economic value of the Khumbu region to Nepal, we conducted two post-monsoon ground assessments of geomorphic and structural damage in the fall of 2015, and a third ground assessment of rebuilding efforts in May, 2016. While landslides and other mass movements were not as extensive as those to the west, numerous landslides and slumps were identified, mostly on steep, convex-up, east-facing slopes where the dip of metamorphic foliation (Namche orthogneiss and related rocks) is towards the Dudh Kosi River channel. Also, building damage was more extensive where villages were located on terraces of unconsolidated sand and gravel perched above the Dudh Koshi River (Nakchun Village). The area upstream from Phakding experienced several co-seismic mass movements on steep slopes that had clearly failed before (previous temblors?), as evidenced by old vegetated slump blocks and steep, unstable hummocky terrain. Particularly hard-hit was the Sherpa village of Thame northwest of Namche Bazaar, constructed on unconsolidated, water-saturated glacial outwash between two lateral moraines at the mouth of the Thame Khola Valley. Aside from Thame, one is struck by the randomness of structural damage in many Sherpa villages, reflecting micro-ground conditions beneath buildings and the style and/or age of construction (drystacked fieldstone versus the use of mud or cement mortar, the presence of wire mesh or gabion bands/spacers within stone walls, etc.). Fresh valleyflank slumps and landslides were markedly less obvious further north (towards Everest) where wide glacial valleys exist, although evidence of large rock falls and other mass movements in the past is clearly present.

Valley profile convexity:

Earthquake-triggered mass movements (past & recent): Traditional and new construction methods: Spectrum of structural damage:

Rebuilding the Khumbu:

The photo above shows the iconic Sherpa village of Namche Bazaar with Kongde Ri (6,187 m, 20,299 ft) in the back-ground. One of the most dramatic aspects of the geomorphology of the Greater Himalaya in general and the Khumbu specifically is the spectacular convexity of deeply incised river valleys. Valley profiles are steepest (nearly vertical) right above river channels, but gradually "lay-back" at higher elevations. This gives the impression that erosion is barely keeping up with tectonic uplift. Also note extreme gradients of tributary "channels" (most are avalanche chutes) above the Bhote Koshi River cut into bedrock below Namche Bazaar. This boulder (photo above) lies along the main trekking trail that passes through Ghat; the many painted Mani stones at Ghat are debris from a massive paleolandslide. Unconsolidated, hummocky and often water-saturated deposits like this form a poor substrate for field-stone buildings when earthquakes strike. Other villages that are constructed on unconsolidated landslide debris, perched high above adjacent rivers, include Khumjung and Namche Bazaar.

Courtesy Mike Stickney Montana Bureau of Mines and Geology

The high Sherpa village of Pheriche is located on a glacial ground moraine surrounded by high lateral moraines. Ground shaking of water saturated glacial sediment resulted in significant damage (inset photo).

Coalesced coarse talus/scree/alluvial cones at the base of Taboche Peak (6367 m; 20,889 ft.), overlain by a massive rock avalanche deposit that has been truncated at its toe by the Lobuche Khola River (which drains directly from the Khumbu Glacier at the foot of Mount Everest). Although this deposit was not triggered by the 25 April 2015 Gorkha earthquake, ancient rock avalanche deposits of this magnitude are good analogues for earthquake hazard analysis in the high alpine valleys of the Greater Himalaya. Note the yak pastures and huts below the slide. A major GLOF hazard lies just behind the terminus of Chola Glacier at the base of Arakam Tse.

Taboche (6367 m) Arakam Tse (6423 m)

Yak pasture with a fence constructed of dry-stacked "field stones" (no mortar) was toppled during the 2015 Gorkha earthquake sequence. This yak pasture complex is located on the high NE lateral moraine (Khumbu Glacier) between Pheriche and Dughla. Unsupported top-heavy construction, unstable substrate, and possible amplification of ground shaking on the moraine crest all contributed to failure. Older, traditional buildings in the Khumbu were constructed of uncut "field stone" with no mortar or external plaster, just small stones placed between larger stone blocks to fill gaps; some buildings used dried mud mortar as an external "plaster," as shown in the diagram and photo above (Lukla).

"X-cracking" of external structural walls was common damage from the Gorkha earthquake, often with cracks originating from corners of windows and door jams. Large, through-going sub-vertical cracks were also common, along with toppling of stones at the top of walls and from triangular spaces directly under roofs. The Buddhist stupa above, just west of Namche Bazaar, shows a dramatic example of x-cracking at a corner.

Traditional construction methods throughout many parts of the Himalaya have utilized hand-cut stone blocks with wood spacers, called "timber bands" or "tie beams." Timber bands tend to isolate smaller portions of a stone wall relative to seismic shaking, thus preventing the entire wall from undergoing catastrophic failure. Timber bands, carved into decorative snakes, are clearly seen in Kathmandu on Hanuman Dhoka Palace (yellow arrows in photo below) in Durbar Square. The old palace experienced far less damage from the Gorkha earthquake than did the newer addition to the palace (white building to the left).

This is typical of structural damage in the Khumbu region. Failure occurred at the top of external rock walls (unreinforced masonry), leaving inside wood paneling exposed but relatively intact. It is likely that this rock wall will be rebuilt in the same style as before the earthquake, due to lack of resources.

Buildings constructed in the aftermath of the 2015 Gorkha earthquake show great variety in the degree to which materials are used to mitigate future seismic damage, depending on available resources. The spectrum includes: 1) hand-cut rectangular stone blocks with dried mud mortar and external mud plaster 2) Hand-cut stone with cement mortar (rare) 3) Hand-cut stone with wire-mesh wrap (Gabion bands), with or without mud or cement mortar 4) Hand-cut stone with wood spacers, as described above 5) Metal sheet siding on upper levels of multi-story structures, instead of stacked rock blocks or field stone; this is expensive but provides for a safer dwelling during seismic ground shaking

External metal siding on upper levels

Example of insufficient timber banding in the wall in a small structure in Phakding; the wood boards do not go through the entire wall and did not support the wall during seismic shaking, resulting in catastrophic failure.

Wire-mesh gabion walls are used in conjunction with expensive steelframe construction at the Khumbu Climbing Center in the Sherpa village of Phortse, resulting in a reasonably strong structure that withstood seismic shaking in 2015.

Catastrophic failure of a traditional Khumbu building (stacked field stone -URM).

Catastrophic failure of a traditional Khumbu building (stacked field stone with external plaster) near Thame. The entire structure is on the verge of collapse (red circle). An elderly man lives here with a young yak, residing in a small tent on the opposite side of the building.

ABOVE: Catastrophic failure of roof and walls of a traditional Khumbu building (stacked field stone with external plaster) near Thame. BELOW: Thame experienced massive damage during the Gorkha earthquake and its aftershock sequence; this famous Sherpa village was constructed on water-saturated ground moraine and outwash between high lateral moraines. Many people still live in make-shift shelters of tarps, sheet metal, or whatever else can be salvaged.

  View looking south down the Khumbu Valley from a vantage point just north of Namche Bazaar showing dramatic valley profile convexity. Near-vertical valley flanks near the Dudh Koshi River channel increase the risk from earthquake-and monsoon-induced mass movements, particularly on perched benches of unconsolidated glacial or alluvial sediment. View looking southward over the village of upper Pangboche at steep east-dipping metamorphic foliation in the Namche orthogneiss which controls valley profile convexity in this area. Valley profile convexity is most apparent in distant profiles. Valley profile convexity dramatically lessens north of Phortse where the effects of alpine glaciation have modified the landscape. Nevertheless, extreme hazards associated with ground shaking exist, including rock falls, glacial outburst floods (GLOFS), avalanches from perched glaciers, and failure of steep glacial moraine deposits, to name a few. This photo shows the Everest massif, with Nuptse, Everest and Lhotse (left-to-right) in the distance, with the Imja Khola River in the channel. "The convexity of the walls of such valleys is the outcome of the increasing erosive power of the rivers, this in turn being due to the steady increase in the height of the mountains." Quote from Arnold Heim and Augusto Gansser, The Throne of the Gods: An account of the first Swiss expedition to the Himalayas (1939) A. Gansser, 1936 Large earthquake-induced landslide on the west flank of the Dudh Koshi River immediately north of Phakding above the Dudh Koshi River. Workers have started rebuilding the trail across the slide-path (Oct. 2015), despite obvious unstable perched boulders above. The unconsolidated and heterogeneous nature of the coarsely-stratified sediment is obvious. ABOVE: Lower portion of the landslide scar (see photo above) on the west flank of the Dudh Koshi River immediately north of Phakding above the Dudh Koshi River. Heterogeneity of the sediment is a mirror image of that higher up. BELOW: Twin landslides from the 25 April 2015 earthquake in perched unconsolidated sediment of mixed sources along the Dudh Koshi River between Ghat and Phakding.Evidence of past compound landslides and slump deposits, some massive, are common along the flanks of the Dudh Koshi River north of Lukla. The bench on which Lukla is built is riddled with enormous angular boulders adjacent to the Hillary landing strip; this large landslide likely dammed the Dudh Koshi River(Götz et al., 2015).

  Ladder crossing of crevasses in the Khumbu Icefall above Everest Base Camp. Shown in red is the break-away zone (between Pumori and Lingtren) and run-out path (air blast zone) of the massive ice/snow avalanche on 25 April 2015 that swept across Base Camp and killed over 20 people. Photo taken in April 2012 during the MSU/Nat Geo/TNF Everest Education Expedition. BELOW: Prayer flags on the bridge over the Dudh Koshi River below the Sherpa village of Phortse.


E Herbert

P.-P Cortet

F Daviaud

B Dubrulle

Eckhaus-like instability of large scale coherent structures in a fully turbulent von K árm án flow

The notion of instability of a turbulent flow is introduced in the case of a von Kármán flow thanks to the monitoring of the spatio-temporal spectrum of the velocity fluctuations, combined with projection onto suitable Beltrami modes. It is shown that the large scale coherent fluctuations of the flow obey a sequence of Eckhaus instabilities when the Reynolds number Re is varied from 10 2 to 10 6 . This sequence results in modulations of increasing azimuthal wavenumber. The basic state is the laminar or time-averaged flow at an arbitrary Re, which is axi-symmetric, i.e., with a 0 azimuthal wavenumber. Increasing Re leads to non-axisymmetric modulations with increasing azimuthal wavenumber from 1 to 3. These modulations are found to rotate in the azimuthal direction. However, no clear rotation frequency can be established until Re ≈ 4 × 10 3 . Above, they become periodic with an increasing frequency. We finally show that these modulations are connected with the coherent structures of the mixing shear layer. The implication of these findings for the turbulence parametrization is discussed. Especially, they may explain why simple eddy viscosity models are able to capture complex turbulent flow dynamics. C 2014 AIP Publishing LLC.

[http://dx.doi.org/10.1063/1.4855018]

I. INTRODUCTION

In classical phenomenology, turbulence arises after a sequence of symmetry breakings, successive instabilities, or bifurcations, which however progressively restore the system symmetries in a statistical sense. [START_REF] Frisch | Turbulence: The Legacy of A.N. Kolmogorov[END_REF] The study of these instabilities traditionally proceeds from (linear or nonlinear) perturbations of the so-called "basic state," the stationary laminar solution of the Navier-Stokes (NS) equation at low Reynolds number. At finite Reynolds number there is no general well-defined criterion to discriminate between the turbulent or laminar nature of the flow: one cannot define clearly a critical threshold beyond which the flow is turbulent and below which it is laminar. The consensus view is that the flow is turbulent when the Reynolds number is large enough and when a well established spatio-temporal energy spectrum is observed with broadband power laws. In such a case, a possible statistical equivalent of the laminar "basic state" can be defined using the (statistically or time) averaged flow. This flow is stationary by construction, but differs from a usual basic state in the sense that it is solution of the ensemble averaged Navier-Stokes equation instead of a solution of the plain Navier-Stokes equation. A natural question then arises: what happens once the statistically stationary turbulent state is reached? Is this the end of the story or can new instabilities of the averaged turbulent flow develop? Experimentally, the answer seems to be positive, since spontaneous bifurcations and flow reversals in fully developed turbulence have already been observed: (i) in a wake flow, there is a mean pattern transition at a critical Reynolds number. This corresponds to the so-called "drag-crisis" (see Ref. 2). This leads to a dramatic decrease of the mean drag of a sphere or a cylinder in a turbulent flow for a critical value of the Reynolds number a) eric.herbert@univ-paris-diderot.fr 1070-6631/2014/26(1)/015103/13/$30.00

C 2014 AIP Publishing LLC 26, 015103-1

Re ∼ 10 5 . The wake becomes narrower as the mean flow pattern changes. (ii) Spontaneous flow reversals have been observed in thermally driven (Rayleigh-Bénard) convection, both experimentally [START_REF] Sugiyama | Flow reversals in thermally driven turbulence[END_REF] and numerically. [START_REF] Van Der Poel | Flow states in two-dimensional Rayleigh-Bénard convection as a function of aspect-ratio and Rayleigh number[END_REF][START_REF] Podvin | Proper orthogonal decomposition investigation of turbulent Rayleigh-Bénard convection in a rectangular cavity[END_REF] In a cylindrical Rayleigh-Bénard geometry, the reversals of large scale circulation were reported to be induced by reorientation along the azimuthal direction. [START_REF] Mishra | Dynamics of reorientations and reversals of large-scale flow in Rayleigh-Bénard convection[END_REF] However, it is difficult to experimentally characterize a bifurcation from a non-stationary or periodic flow. (iii) In a magnetohydrodynamic turbulence, the large scale magnetic field spontaneously generated in a liquid metal at Re > 10 6 was shown [START_REF] Berhanu | Magnetic field reversals in an experimental turbulent dynamo[END_REF] to undergo spontaneous reversals, with a dynamics governed by a few magnetic modes despite the strongly turbulent background. This observation led to a simple model of geodynamos with a few modes undergoing chaotic dynamics. [START_REF] Gallet | Reversals of a large scale field generated over a turbulent background[END_REF] An open problem common to all these situations is to find suitable tools to study and characterize these types of instabilities.

The purpose of this paper is to present some answers to these questions in the specific case of a von Kármán flow. This flow is generated by two counter-rotating impellers in a cylindrical vessel. At low Reynolds number, the laminar flow is axi-symmetric and divided into two toric recirculation cells separated by an azimuthal shear layer. Its transition to turbulence for various counter-rotations and aspect ratio has been extensively studied theoretically, numerically, and experimentally with smooth or rough (fitted with blades) impellers (see, e.g., Refs. 9-16 and references therein). The study in Ref. 13 was performed in the exact-counter-rotating case. It was reported that a first bifurcation from the stationary and axi-symmetric laminar state occurs at Re = 175, yielding a stationary flow with an azimuthal modulation, i.e., for which the axisymmetry is broken. This stationary state persists up to typically Re ∼ 300 where time-dependence arises. The transition to turbulence further proceeds through plain, modulated, or chaotic traveling waves, until Re ∼ 10 4 , where a "fully developed turbulent state" seems to have been reached. In that state, the energy spectrum is broad and the dimensionless dissipation does not depend on Reynolds number anymore. [START_REF] Ravelet | Supercritical transition to turbulence in an inertially driven von Kármán closed flow[END_REF] In that turbulent state, the flow once time averaged regains the structure of the basic laminar flow, made of two shearing toric recirculation cells. In the present paper, we have worked using the same experimental setup. We demonstrate that the flow is organized into large scale coherent structures. These structures are subject to a sequence of Eckhaus type instabilities leading to an increasing azimuthal wavenumber as the Reynolds number varies from 10 2 to 10 6 , similar to the sequence of instabilities at much lower Re. Since this sequence of instabilities occurs on coherent structures in a turbulent flow rather than a laminar "basic state," we are facing a new paradigm that both require new tools of investigations and open interesting questions about turbulence parametrization. The tools are presented in Sec. II, the instability is described in Sec. III and some theoretical consequences are discussed in Sec. IV.

II. TECHNICAL BACKGROUND AND TOOLS

A. Experimental setup

We have worked with a von Kármán flow generated by two coaxial and counter-rotating impellers in a cylindrical vessel. More details about the experimental setup can be found in Ref. 17. The cylinder radius and height are R = 100 mm and H = 180 mm, respectively. The impellers consist of 185 mm diameter disks fitted with 16 curved blades of 20 mm height. The impellers are driven by two independent motors. The motor frequencies are, respectively, set to f 1 and f 2 . In the present work, the impellers are rotating with the convex face of the blades going forwards, contrarily to Ref. 13. The working fluid is either pure water, a water-glycerol mixture (26%-74% in weight, respectively), or pure glycerol. The resulting accessible Reynolds numbers, Re = π (f 1 + f 2 )R 2 ν -1 with ν the kinematic viscosity, vary from 10 2 to 10 6 . This setup is invariant under R π rotations around any radial axis passing through the center of the cylinder. Measurements are done thanks to a Stereoscopic Particle Image Velocimetry (S-PIV) system. The S-PIV provides time series of the 3 components of the velocity (radial v x (x, z, t), vertical v z (x, z, t), and azimuthal v y (x, z, t)) on a 63 × 58 points grid in a meridian plane [x, z] composed of two planes dephased by π in the cylindrical coordinates. x = rcos (ϕ 0 ), with the azimuthal angle ϕ 0 = [0, π ] and r the radial distance from the impellers axis. z is the vertical distance from the center of the cylinder, y is the out-of-plane direction, and t is the time. The azimuthal component is then written in the cylindrical coordinates using v ϕ (x, z, t) = x/|x| v y (x, z, t). The size of the recorded Phys.

B. Experimental perturbation analysis

For low Reynolds number Re ≤ 500, the flow is steady and laminar and only few fluctuations are present. When θ = (f 1f 2 )/(f 1 + f 2 ) = 0, the instantaneous flow is composed of two toric recirculation cells separated by an azimuthal shear layer located at z = 0 reflecting the R π -symmetry of the system. Increasing the Reynolds number, one expects to reach fully developed turbulence around Re = 10 4 as observed in Ref. 13, with increasing difference between the instantaneous flow and the time-averaged flow. Above Re = 500 and at θ = 0, the axisymmetry is broken by fluctuations of the instantaneous flow. However, this symmetry is restored for the time-averaged flow (cf. Fig. 1), at any Reynolds number. By construction, this average flow is time-independent and will be the equivalent of a "basic state." As the Reynolds number is increased, this "basic state" is subject to increasingly stronger fluctuations. [START_REF] Monchaux | Properties of steady states in turbulent axisymmetric flows[END_REF] The question we want to address here is: are these fluctuations purely disorganized as in thermal noise or can we identify a pattern in them when varying the Reynolds number, [START_REF] He | Thermodynamical versus log-Poisson distribution in turbulence[END_REF] in a way similar to a sequence of bifurcations of ordinary instabilities? To answer this question we need two steps: (i) to separate the fluctuations into "organized" motions and "thermal" noise and (ii) to build a suitable tool to study the temporal behavior of these organized motions. These steps are described below. For step (i), we use comparison with and projection onto Beltrami modes. For step (ii), we use spatio-temporal spectra, as used in wave turbulence. [START_REF] Cobelli | Space-time resolved wave turbulence in a vibrating plate[END_REF]

Projection onto Beltrami modes

The measured velocity fields time series are analyzed thanks to a projection on the basis of Beltrami modes as detailed in Ref. 21. Beltrami modes have been introduced as a general spectral decomposition basis, [START_REF] Constantin | The Beltrami spectrum for incompressible fluid flows[END_REF] and any velocity field in a cylindrical geometry can be decomposed as a superposition of such modes. Given a vector field, v, lying in a cylinder of radius R = 1 and height

H = 2h, with usual cylindrical components v = (v r , v ϕ , v z ), we first switch to the new components V = (V + , V -, v z ) such that V ± = (v r ± iv ϕ )/2.
The field V can further be decomposed over Beltrami modes B nmks as where k 0 = π /h, n and m are integer numbers, and (N, M, P) are the number of considered modes in the (r, ϕ, z) directions, respectively. The D nmk are complex amplitudes which conjugates verify D * nmk = D n,-m,-k and the Beltrami modes are given by

V = N n=1 M m=-M Pk 0 k=-Pk 0 (+) s=(-) D nmk B nmks , ( 1
B nmks = ⎛ ⎜ ⎝ B nmks + B nmks - B nmk z ⎞ ⎟ ⎠ = 1 2 ⎛ ⎜ ⎝ (λ s -k) J m+1 (μ nm r ) (λ s + k) J m-1 (μ nm r ) -2iμ nm J m (μ nm r ) ⎞ ⎟ ⎠ exp(imϕ + ikz). ( 2
)
Here, λ s = s μ 2 nm + k 2 with s = ± and J m is the Bessel function of order m. The orthogonality condition is ensured by the fact the μ nm coefficients are the nth root of J m . Equation ( 1) then corresponds to a decomposition into solenoidal Beltrami waves with polarization given by the sign s of λ s . From now on, we will omit the s superscript for simplicity. We define the synthetic velocity 2(a). We use the boundary condition explained in Ref. 21. When m > 0 the velocity fields depend on the azimuthal angle ϕ and the axi-symmetry is lost. However, for each meridian plane at a fixed [ϕ i , ϕ i + π ] the velocity field shows mirror (anti-)symmetry with respect to the r = 0-axis. Even-numbered m velocity fields are found to be symmetric while odd-numbered are found to be anti-symmetric. Introducing the space average = 1 0 rdr π -π dϕ h -h dz and using the classical properties of the Bessel functions, the orthogonality of B nmk writes

field b nmk = (b r , b ϕ , b z ) corresponding to B nmk = (B + , B -, B z ) using the real parts of, respectively, b r = B -+ B + , b ϕ = i(B --B + ),
B nmk • B n m k * = B nmk + B n m k * + + B nmk -B n m k * - + B nmk z B n m k * z , = B nmk • B nmk * δ mm δ nn δ kk . ( 3
)
This means that for any field satisfying the decomposition (1), we have

V • B nmk * = D nmk B nmk • B nmk * , ( 4
)
which provides a simple way to find the projection of the vector V on B nmk by spatial average over the fluid volume. In the experiment, we however have access to velocity measurements in one plane only, corresponding to ϕ = 0 and ϕ = π . This precludes the exact instantaneous projection onto given Beltrami modes. We remedy to these problems by two methods: (i) conditional time averages at a given phase mϕ. This is explained in Sec. III E. (ii) Direct comparison with "synthetic measurements." This is described in Sec. III B. These projections, limited to low order Beltrami modes, allow to reconstruct synthetic velocity field time series in the full 3D space, accounting for the full spatiotemporal evolution of the "organized" motions of the turbulent flow.

Spatio temporal spectrum

In traditional instability analysis, it is customary to draw space-time diagrams to detect wavelike pattern, see Refs. 23 and 24. In a turbulent system with a wide range of scales, a natural generalization is to resort to spatio-temporal spectra. This approach is routinely used in the field of wave turbulence (see Ref. 25 for a recent review) and led to computation of dispersion relations, see Ref. 20. To obtain these spectra from our experimental PIV measurements, we proceed as follows. From time series of the azimuthal velocity v ϕ (x, z, t) in the meridian plane we compute the time-averaged velocity field vϕ and its fluctuations v ϕ (x, z, t) defined by v ϕ = v ϕvϕ . From the corresponding time series, the full space-time (2D in space and 1D in time) power spectrum E(k x , k z , f) is computed, where k = K/2π with K the wavevector and f = ω/2π is the frequency. This is done in two steps. First, the instantaneous spatial (2D) Fourier Transform ṽ 2D (k x , k z , t) for each time t of the time series is computed. Then, the temporal Fourier Transform of ṽ 2D using a window size composed of n time steps is computed. Finally, ṽ 3D (k x , k z , f ) is obtained that leads to the spatio-temporal power spectrum in 3 dimensions E(k x , k z , f):

E(k x , k z , f ) = | ṽ 3D | 2 = dx dz dtv ϕ (x, z, t) e i(ωt+k x x+k z z) 2 . ( 5
)
Spatial (horizontal and vertical) and temporal resolutions of the spectra are comprised in the ranges (1/(2R),1/(2δ R )), (1/H,1/(2δ Z )), and (f s /n,f s /2), respectively.

Synthetic time series

Direct comparison of the experimental spectra with synthetic ones, based on a synthetic velocity fields time series, is then performed. The latter were obtained with single Beltrami modes of a given (cylindrical) geometry. For simplicity, we have restricted our comparison with modes with largescale poloidal structure n = k/k 0 = 1 and varying azimuthal structure (variable m). As soon as m > 0, the axisymmetry is broken. Like in any rotation-breaking instability, it is natural to expect that the corresponding perturbation mode will rotate with a frequency f r (that is a priori small with respect to the impellers rotation rate (f 1 + f 2 )/2). The corresponding Beltrami modes will then have the structure:

B 1mk 0 ∼ exp(imϕ + ikz -i f r t). ( 6
)
An artificial time series of such a mode is build in the meridional plane [x, z]. These synthetic data are then processed like the experimental signal to get synthetic space-time power spectra. Since the perturbation Beltrami mode is periodic (with m the number of period) in its azimuthal direction, a peak of energy is located at a frequency f = mf r as can be seen in Fig. 5(h). This peak is well resolved by setting the sampling frequency to a much higher value, typically 50 × mf r . In the sequel, we will refer to the spectra obtained using the experimental (resp., synthetic) data by E E (resp., E S ).

III. RESULTS

A. Isotropy: Evidence for the Eckhaus instability

A first characterization of the flow perturbation independently of its dynamics can be found by considering the spatial spectrum E(k x , k z ), obtained through the integration of the 3D spectrum E(k x , k z , f) over all temporal frequencies. An example for an experimental flow at Re ≈ 10 6 is provided in Fig. 3. The energy is mainly concentrated along the x-axis at k z ≈ 0 inside an ellipsoidal region, pointing out the anisotropic structure of the flow. For useful interpretation, we have compared this isotropy measurements with synthetic data. In Fig. 3, we also show an example obtained with n = k/k 0 = 1, m = 3, and f r = 0.2 Hz synthetic velocity fields time-series. For this field, the overall shape of the spatial spectrum is close to the experimental one. Both spatial spectrum E E (k x , k z ) and E S (k x , k z ) in Fig. 3 display 2 distinct maxima. They are found around k z = 0 for two values of k x symmetric with respect to k x = 0. As we can see in Fig. 2 using the synthetic velocity fields or in Fig. 9 using the projection of the experimental fields onto the Beltrami modes, the x-axis traces the Phys. Fluids 26, 015103 (2014) parity of m. This parity is reproduced through the spectral analysis in this way: in the even-numbered (resp., odd-numbered) case, the symmetry (resp., anti-symmetry) of the flow with respect to r = 0 leads to a fundamental wavelength of roughly the size (resp., half of the size) of the vessel. The wavenumbers k x ±1/(2R) corresponding to the maxima of energy in the symmetric case approach the resolution limit of the Fourier transform and practically appear as k x 0 (see Figs. 5(e) and 5(g)). On the contrary, the anti-symmetric case leads to two distinguishable symmetric peaks at k x ≈ ±1/R (see Figs. the maximum of the 1D-spatial spectrum shifts back and forth between k x = 0 and two symmetric values k x = ±7 m -1 . Therefore, the kind of transitions observed in Fig. 4 is suggestive of successive azimuthal changes in parity (even-odd-even-odd). This means that the flow is composed of a large scale coherent structure which is experiencing an Eckhaus type instability in its azimuthal direction starting from the axi-symmetric time-averaged flow which is a m = 0 mode. To bring some clues to the upper modes where m > 0 we resort to a different analysis of the spatial spectrum. The contour line at a fixed iso-energy level was fitted using a centered ellipsoidal curve, (x/w) 2 + (z/ h) 2 = 1, with the width w and the height h as fitting parameters. As a result w/ h is a measure of the anisotropy of the flow. In order to obtain error bars, the ratio w/ h was computed for various iso-energy levels corresponding to different fractions of the total energy of the spectra, from 65 to 90% (resp., 75% to 99%) for E E (resp., for E S ). The mean value A = w/ h and the standard deviation are finally computed to evaluate both the anisotropy and the error. Using the synthetic data, varying m leads to different values of A which is moreover found not to depend on f r . Specifically, A = 2.3 ± 0.3, for m = 0. For m = 1, A = 5.2 ± 2.4 reach its maximum. Increasing m, A tends to decrease, yielding for m = 2, 3, and 4 to A = 2.9 ± 0.8, 2.4 ± 0.3, and 2.3 ± 0.3, respectively. We have then studied the experimental evolution of the anisotropy ratio with the Reynolds number. This is shown in Fig. 4. One sees the error bar can be large in particular below Re = 5000. Consequently, this approach is not precise enough to determine at which Reynolds number the transitions occurred. However, we clearly distinguish different regimes. The anisotropy starts from a value of about 1.1 at very low Reynolds number where the amplitude of the fluctuations reach the resolution limit of our PIV. Increasing the Reynolds number, A increases up to 2.7 ± 0.7 corresponding to the laminar axi-symmetric (m = 0) case where A = 2.3 or to the time-averaged flow at any Reynolds number. The amplitude of the anisotropy then reaches a maximum (4.2 ± 1) at Re = 800. This last value is compatible with the odd-numbered m = 1 mode. At approximately Re = 5000, A reaches a minimum at 2.3 ± 0.1. Which seems not compatible with even-numbered m = 2 mode but could be a m = 4 mode. Finally, for Reynolds number larger than 10 4 , A shows a plateau around A = 2.5, which value is compatible with the odd-numbered m = 3 mode.

k z [m -1 ] k x [m -1 ] k z [m -1 ] k x [m -1 ]
k z [m -1 ] k x [m -1 ] k z [m -1 ] k x [m -1 ]

B. Temporal behaviour

A symmetry breaking generally coincides with a propagative instability. [START_REF] Coullet | Instabilities of one-dimensional cellular patterns[END_REF] To check this hypothesis and try to measure the rotating frequency f v of the azimuthal structure of the turbulent flow, we used the 2D-spatio-temporal spectrum E E (k x , f) which is computed through the integration over k z of the 3D-spectrum E(k x , k z , f). E E (k x , f) is then compared to the spectrum E S (k x , f) obtained using a synthetic Beltrami flow for different values of m. This is done in Fig. 5 using different Reynolds numbers. The rotation frequency f r of the synthetic flows is adjusted to get the best resemblance with the experimental field which is guessed to rotate at a f v frequency. At Re = 314, the experimental spectrum is found very close to E S (k x , f) with m = 0. This result can be extended to the time-averaged flow spectra (not shown) at any Reynolds number since it is similar to the laminar instantaneous flow. Increasing the Reynolds number to Re = 1473 (resp., Re = 7200) in a region where one (resp., two) Eckhaus bifurcation(s) has occurred, we observe a strong similarity with corresponding E S (k x , f) with m = 1 (resp., m = 2) with a small but non-zero rotation frequency f r = 0.01 Hz. For the last two cases, the rotation frequency f v of the azimuthal structure of the flow shall be non-zero. However, it reveals to be smaller than our spectral resolution, typically 0.01 Hz, and therefore remains unmeasurable. Increasing the Reynolds number the energy is found to remain mainly concentrated at a 0 frequency up to Re = 4 × 10 4 . In the vicinity of this last Re the energy becomes distributed on a wider range of frequencies, up to 1.5 Hz. However, no particular frequency is selected. Increasing further the Reynolds number, peaks at non-zero frequencies come up, as shown for Re = 10 6 . These peaks are the signature of a fixed rotation frequency at f v = 0.7 Hz. In this case, the experimental spectrum is regained by E S (k x , f) with m = 3 and f r = 0.2 Hz. Overall, our results are compatible with a rotating frequency bifurcating from a small but non-zero value around Re = 4 × 10 4 .

C. Rotation direction

At large enough Reynolds number, we have seen that a m = 3 perturbation of the "basic state" of the flow is established with a non-zero rotating frequency. In the perfectly symmetric case there is, however, no general argument imposing the direction of rotation of the pattern. Experimentally, this should be the case at θ = 0. However, our setup was not perfectly symmetric what led to observe the same direction of rotation for all experiments. A slight shift of θ was then introduced to compensate this lack of symmetry. Using this adjustment we observed for each experiment one or the other direction of rotation with an equal probability. The direction of rotation is best seen by focusing on the vertical direction, i.e., looking at the 2D-spatio-temporal spectrum E E (k z , f), as shown in Fig. 6. One observes either an upward (k z > 0 when f > 0) or a downward (k z < 0 when f > 0) tilt of the spatio-temporal spectrum. This can be modelled by considering a synthetic Beltrami spectrum corresponding to m = 3 with two opposite rotation directions, leading to an upward or downward tilt of the spectrum (also shown Fig. 6).

D. Link with the shear layer dynamics

It is tempting to associate the Eckhaus azimuthal instability evidenced in Sec. III A to a Kelvin-Helmoltz instability of the shear layer. Indeed, this instability was already pointed out as the physical mechanism for a stability loss of the axi-symmetric steady-state followed by successive flows with non-zero azimuthal wavenumber. This was numerically observed for small Reynolds number (Re < 500) in a geometry [START_REF] Nore | The 1[ratio]2 mode interaction in exactly counter-rotating von Kármán swirling flow[END_REF] close to ours. Here, we report a similar mechanism but at largely higher Reynolds numbers. It was also observed that the Kelvin-Helmoltz instability results in a series of large scale structures taking the shape of radial vortices with well defined azimuthal wavenumber that can be seen, e.g., through bubble seeding and white light. [START_REF] Ravelet | Supercritical transition to turbulence in an inertially driven von Kármán closed flow[END_REF][START_REF] Cortet | Normalized kinetic energy as a hydrodynamical global quantity for inhomogeneous anisotropic turbulence[END_REF] The possibility of changes in m of this pattern revealing an Eckhaus instability analog to that existing in small dynamical systems was already mentioned in Ref. 17. To sustain this possibility, we have monitored the 2D-spatial spectrum at a fixed Reynolds number, as a function of θ . It was already reported (see, for example, Ref. 17) a shifting of the shear layer toward the slower impeller with increasing |θ |. When |θ | > 0.1 the shear layer was found to be completely absorbed. If the azimuthal structures we observed are related to the shear layer, the structure of the 2D-spatial spectrum should also follow θ . This is indeed the case, as can be seen in Fig. 7 both from the isotropy and from the 1D-spatial spectrum E E (k x , θ ). With increasing |θ | the two bands of k x merge into one, at k x = 0. At |θ | = 0.1 the isotropy measurement and the peak at k x = 0 of the spectrum drops suddenly. Further increasing θ , both the isotropy and the spatial spectrum (not shown) are found to be independent of θ . This is indicative of a transition into a m = 0 structure, with no azimuthal modulations anymore. This thus proves that the azimuthal modulation is indeed strongly associated with the shear layer.

E. Visualization of the coherent structure

Finally, the existence of organized motions in the turbulent flow means v can be written as the sum of an in-phase term and a fluctuation term: v = v c + v , with v the fluctuations and v c the coherent structure velocity field. A visualization of the spatial structure of one azimuthal mode m can then be obtained using conditional average of the velocity field v . The condition is given by suitable projection onto the corresponding Beltrami mode, see Eq. ( 4). Due to the geometry of our PIV device, we however only have access to the velocity field in a meridional plane, i.e., for ϕ 0 = [0, π ]. Since the azimuthal modes are rotating we may however replace the azimuthal integration in the scalar product by a suitable time integration. This is made possible for high enough Reynolds number where the rotation frequency of the coherent structure is resolved. We used conditional average as follows: we first compute for each instantaneous experimental field v the quantity: in two planes, dephased by π . One indeed sees large scale structures looking like pairs of azimuthal vortices. By reproducing this conditional average using different phase ϕ 0 for B 1mk 0 at a fixed m, we obtained the complete movie (see Fig. 9) of the rotation of this pair of vortices (see the Multimedia view). The physical mechanism leading to this vortices can be the Kelvin-Helmoltz instability laying in the shear layer, a secondary or related instability or the centrifugal instability. It is interesting to compare this result with the m = 3 velocity field shown in Fig. 2(d) using the same phase ϕ 0 as in Fig. 9(a). Both the anti-symmetric nature of the velocity field with respect to the r = 0 vertical axis and v(r = 0) = 0 are recovered.

a m (t) = 1 -1 x dx h/R -h/R dz v (t) • B 1mk 0 (x, ϕ 0 , z), (7)

IV. SOME THEORETICAL CONSEQUENCES

In 1970, Kraichnan wrote a paper entitled "Instability in fully developed turbulence," [START_REF] Kraichnan | Instability in fully developed turbulence[END_REF] in which he focuses on the propagation of errors from one range of scale to another. This problem is related to predictability of turbulence. In the present work, we consider another aspect of the same problem related to the stability of the unsteady coherent structures of a turbulent flow (in other words, to the dynamic of the perturbation of this given coherent structures). We have shown that using spatio-temporal spectra, and comparison or projection onto suitable Beltrami modes, it is possible to identify coherent perturbations to the time-averaged flow. This is liable to an Eckhaus instability analog to what is observed in systems in the laminar limit, with a few degrees of freedom. In classical phenomenology of turbulence the number of degrees of freedom scales like Re 9/4 . Given our range of Reynolds number (from 10 2 to 10 6 ), we thus can expect our system to be described by a number of degrees of freedom in the range 7 × 10 4 to 3 × 10 13 . Our results thus prove that a huge number of degrees of freedom are irrelevant to describe the instability of our turbulent flow. This sets a number of interesting questions that we list below.

In the present paper, we have characterized experimentally the instability of the basic state. Is there any analytical or numerical way to explain our observations? Classical instability analysis for a given velocity field v(x, t) starts from stationary solutions of the Navier-Stokes equation:

∂ j v j = 0, (8)
∂ t v i + v j ∂ j v i = - 1 ρ ∂ i p + f i + ν∂ k ∂ k v i .
In our case, we consider the instability of a time-average basic state v(x, t) that obeys the equation:

∂ j v j = 0, (9)
∂ t v + v j ∂ j vi = - 1 ρ ∂ i p + fi + ν∂ k ∂ k vi + ∂ j R i j ,
where ρ R i j = ρ ( vi v j -v i v j ) is the Reynolds stress tensor. This Reynolds stress represents the influence of all the degrees of freedom of the flow onto its average, and can, in general, only be computed via full solution of the NS equation. Therefore, the problem of instability of a mean turbulent flow cannot be tackled analytically or is too demanding numerically, unless a prescription (parametrization) of the Reynolds stress is provided. In the case of the plane Couette turbulent flow, for example, this was attempted by Tuckerman et al. [START_REF] Tuckerman | Instability of Uniform Turbulent Plane Couette Flow: Spectra, Probability Distribution Functions and Kappa-Omega Closure Model[END_REF] via the K-closure model. They calculate steady 1D solution profiles of the K-model and their linear stability to 3D perturbations, but find no correspondence between this analysis and the onset of turbulent-laminar bands in experiments and simulations. In the same way, Legras and Villone 30 address the problem of Kolmogorov flow instability when molecular viscosity is replaced by a Smagorinsky parametrization for small-scale turbulence. Such a parametrization represents the motion at scales smaller than the large scale Kolmogorov flow. They claim that it may provide hints on large-scale instabilities at large Re and, hopefully, on the character of such instabilities, but it has never been checked. In our case, it was observed that the sequence of Eckhaus instability resembles the sequence of instability observed at much lower Reynolds number. [START_REF] Nore | The 1[ratio]2 mode interaction in exactly counter-rotating von Kármán swirling flow[END_REF] This suggests that the replacement of the Reynolds stress by an eddyviscosity might be a suitable parametrization by increasing the dissipation, and therefore decreasing the effective Reynolds number (and the effective number of degrees of freedom). Theoretically, we note that the special geometry of our experiment and its parity symmetry properties make it plausible that the first term in the expansion of the Reynolds stress as a function of the velocity gradient is indeed quadratic, resulting in a non-isotropic eddy-viscosity tensor (in other words, there is probably no anisotropic kinetic alpha (AKA) term in the Reynolds stress expansion). [START_REF] Dubrulle | Eddy viscosity of parity-invariant flow[END_REF] It is, however, not clear that the correct parameters of the instability (threshold, wavenumbers) can be captured with a simple eddy-viscosity model, since this procedure fails in the case of a plane Couette flow. [START_REF] Tuckerman | Instability of Uniform Turbulent Plane Couette Flow: Spectra, Probability Distribution Functions and Kappa-Omega Closure Model[END_REF] An interesting alternative would be to derive directly the mean state as the critical points of a suitable Arnold functional, making the problem liable to traditional tolls of bifurcation theory. Such strategy was followed in Ref. 32, resulting in a full description of a spontaneous parity-breaking turbulent bifurcation. The generalization of this model to non-axisymmetric mean state is however an open problem. Summarizing, we have shown that the fluctuations of a turbulent von Kármán flow obey a sequence of Eckhaus-like instabilities with varying Reynolds number that is similar to the sequence of Eckhaus instability observed in the laminar flow, at much lower Reynolds number and that may be explained through a suitable parametrization by increasing the dissipation, and therefore decreasing the effective Reynolds number (and the effective number of degrees of freedom). If this is indeed true, it may have interesting implications for other fields, especially astrophysics and geophysics where Reynolds number are huge and our interest is mainly in the dynamics of the large-scale mean flow. For example, current climate models that have presently very low effective Reynolds number may be able to capture efficiently certain types of instabilities of the mean (climatic) turbulent state. They may, therefore, be more predictable that Kraichnan thought, in a very different meaning.

(Received 21

 21 November 2013; accepted 9 December 2013; published online 9 January 2014)

Fluids 26 6 FIG. 1 .

 2661 FIG. 1. Left, time-averaged velocity field v, measured in a meridian plane using S-PIV, at Re = 10 6 and θ= (f 1f 2 )/(f 1 + f 2 ) = 0. x = rcos ϕ 0 , ϕ 0 = [0, π ]. Right, corresponding instantaneous fluctuating component of the velocity field v = vv.The color maps the azimuthal velocity v ϕ (normalized by R(f 1 + f 2 )/2) whereas the arrows map the (v x ,v z ) velocities.

FIG. 2 .

 2 FIG. 2. (a) (resp., (b), (c), and (d)) synthetic velocity field b nmk obtained from corresponding Beltrami modes B nmk in a meridional plane with n = k/k 0 = 1 and m = 0 (resp., m = 1, 2, and 3). Same color and arrow codes as in Fig. 1.

  and b z = B z and a fixed sign of λ. Some examples of the synthetic velocity fields obtained for k/k 0 = n = 1 and m = 0, 1, 2, 3 in the two planes ϕ 0 = [0, π ] are provided in Fig. 2. The recirculation cells, the shear layer, and the axi-symmetry of the experimental time-averaged flow shown in Fig. 1 are well recovered by the ϕ-invariant m = 0 mode shown in Fig.

FIG. 3 .

 3 FIG. 3. Left, experimental spatial power spectrum E E (k x , k z ) at θ = 0 and Re ≈ 10 6 , from time series from which the flow shown Fig. 1(b) is extracted. White line is an ellipsoidal fit of a fixed energy level contour line. w/ h = 2.5, with w the width and h the height of the ellipse. Right, synthetic spatial power spectrum E S (k x , k z ) with n = k/k 0 = 1 and m = 3, corresponding to the flow shown Fig. 2(d) rotating at f r = 0.2 Hz and for which w/ h = 2.3. Colors are log scaled.

  5(f) and 5(h)). Experimentally, these different patterns are reproduced using the Reynolds number as a driving parameter. Figures5(a) and 5(c) at Re = 714 and 7200 show the symmetric (even) case whereas Figs. 5(b) and 5(d) at Re = 1473 and 10 6 show the anti-symmetric (odd) case. Consequently, the flow is experiencing successive transitions in parity. To study more precisely the evolution of this pattern with the Reynolds number, one can concentrate on the 1Dspatial spectrum E E (k x ), obtained through integration of E E (k x , k z ) over all k z . This was done for all values of Re resulting in Fig. 4. One observes clear transitions at Re = 700, 3200, and 8800, where

FIG. 4 .

 4 FIG. 4. Top, anisotropy factor A = w/ h computed from experimental power spectra as a function of the Reynolds number (see details in the text and illustration in Fig. 3). Marker shapes (•, ×, +) represent the different working fluids. Bottom, experimental E E (k x ) at θ = 0 as a function of the Reynolds number. Each spectrum is normalized by its maximum. Colors are log scaled. Vertical white lines mark strong transitions in the E E (k x ) spectrum which highlight parity changes in the structure of the experimental flow. Black and white circles, see Fig.5.

FIG. 5 .

 5 FIG. 5. Top, experimental spatio-temporal spectrum E E (k x , f) at θ = 0 for Re = 314 (a), 1473 (b), 7200 (c), and 10 6 (d), corresponding to circles in Fig. 4. Colors are log scaled in the range [-1.5, 0]. Bottom, synthetic spatio-temporal spectrum E S (k x , f), using n = k/k 0 = 1 and f r = 0.01 Hz for m = 0 (e), 1 (f), 2 (g), and f r = 0.2 Hz for m = 3 (h). Colors are log scaled in the range [-2, 0] for ((e)-(g)) and [-6, 0] for (h). Note the frequency axis range change for (d) and (h).

FIG. 6 . 1 A

 61 FIG. 6. Left, experimental spatio-temporal spectrum E E (k z , f) for (a) θ = 0, Re ≈ 10 6 and (c) θ = -8 × 10 -3 , Re ≈ 8 × 10 5 . Right, synthetic spatio temporal spectrum E S (k z , f) with n = k/k 0 = 1, m = 3 and using (b) f v = 0.1 Hz and (d) f v = -0.1 Hz. Colors are log scaled.

FIG. 8 . 1 FIG. 9 .

 819 FIG. 8. (a) a m (t) = xdx dzv • B 1mk 0 with m = 3 and Re = 10 6 at a fixed angle ϕ (corresponding to the angle of Fig. 2(c)) as a function of time. a m (t) was computed using v (t) from the time-series used in Fig. 5(d) and was finally low-pass filtered at 4 Hz. Horizontal lines, see Fig. 9. (b) Power spectrum E am ( f ) of a m (t) (+ solid line) for the data shown in (a) and for other Reynolds numbers: Re/1000 = 41 (thin red/gray line), 21 (thin black line), 12 (thick black line), 7 (•). The arrow marks the main peak evolution with decreasing Re.

ACKNOWLEDGMENTS

E. Herbert acknowledges the support of the program DSM-Energie and is grateful to V. Padilla and C. Wiertel-Gasquet for techniques and B. Saint-Michel for interesting discussions.