Carsten Riedel
email: carsten.riedel@awi.de

Vera Schlindwein

Did the 1999 earthquake swarm on Gakkel Ridge open a volcanic conduit? A detailed teleseismic data analysis

Keywords: Gakkel Ridge, HYPOSAT, Relocation, Station residuals, Spatiotemporal organization, 1999 swarm, Polar research

In 1999, a seismic swarm of 237 teleseismically recorded events marked a submarine eruption along the Arctic Gakkel Ridge, later on also analyzed by sonar, bathymetric, hydrothermal, and local seismic studies. We relocated the swarm with the global location algorithm HY-POSAT and analyzed the waveforms of the stations closest to the events by cross-correlation. We find event locations scattered around 85 • 35 N and 85 • E at the southern rift wall and inside the rift valley of the Gakkel Ridge. Waveforms of three highly correlating events indicate a volumetric moment tensor component and highly precise referenced double-difference arrival times lead us to believe that they occur at the same geographical position and mark the conduit located further southeast close to a chain of recently imaged volcanic cones. This result is supported by station residual anomalies in the direction of the potential conduit. Seismicity is focused at the crust-mantle boundary at 16-20 km depth, but

Introduction

The Arctic Gakkel Ridge (Fig. 1) is one of the slowest-spreading zones of the Earth's lithosphere [START_REF] Edmonds | Discovery of abundant hydrothermal venting on the ultraslow-spreading Gakkel Ridge in the Arctic Ocean[END_REF]) with a full spreading rate lower than 1.2 cm/year. Along the rift axis, volcanic and amagmatic units alternate [START_REF] Michael | Magmatic and amagmatic seafloor generation at the ultra-slow Gakkel Ridge, Arctic Ocean[END_REF]. Details of the active rifting processes are still unknown due to the remote location of Gakkel Ridge and its cover with Arctic sea ice which only few research vessels can tackle. However, in 1999, a sequence of earthquakes of body wave magnitude 3.3 (International Seismological Centre [ISC]) and higher was recorded worldwide [START_REF] Müller | Seismic evidence for volcanic activity discovered in Central Arctic[END_REF][START_REF] Tolstoy | Seismic character of volcanic activity at the ultraslow-spreading Gakkel Ridge[END_REF], which is so far the largest teleseismically recorded earthquake swarm at any mid-ocean ridge (MOR) volcanic system and occurred in the area surrounding 85 • 40 N, 85 • E (Fig. 1). Global centroid moment tensor solutions with large isotropic components Fig. 1 The ISC locations of the 237 events mentioned in the text have been plotted as dark gray circles on a bathymetric map of Gakkel Ridge in the survey area on the left side of the figure . A scale for the bathymetry where black is related to topographic lows and white is related to topographic highs is displayed at the top. On the right side the position of the swarm area is marked as a white box on a plot of the high latitudinal areas of the northern hemisphere. Below, a histogram of magnitudes is presented. The size of the circles in the location graph is proportional to the square values of magnitude (scale see below magnitude distribution) to ease the view of high-magnitude events.

A small white triangle (volcano) marks the position of the crater of Oden volcano which was discovered by [START_REF] Sohn | Explosive volcanism on the ultraslow-spreading Gakkel Ridge, Arctic Ocean[END_REF]. All figures including this one have been plotted using Generic Mapping Tools [START_REF] Wessel | New, improved version of generic mapping tools released[END_REF] and high b values indicated a typical volcanic character [START_REF] Müller | Seismic evidence for volcanic activity discovered in Central Arctic[END_REF] and sonar images from the USS Hawkbill only 4 months after the main activity showed areas of high backscatter interpreted as fresh lava flows [START_REF] Edwards | Evidence of recent volcanic activity on the ultraslow-spreading Gakkel Ridge[END_REF].

During the joint AMORE 2001 cruise on the USCGC Healy and RV Polarstern rock sampling, bathymetric surveys and hydrothermal studies discovered that the seismically active area represents a major volcanic complex [START_REF] Michael | Magmatic and amagmatic seafloor generation at the ultra-slow Gakkel Ridge, Arctic Ocean[END_REF]. A giant hydrothermal plume in the water column indicated sustained volcanic activity [START_REF] Edmonds | Discovery of abundant hydrothermal venting on the ultraslow-spreading Gakkel Ridge in the Arctic Ocean[END_REF][START_REF] Baker | Hydrothermal venting in magma deserts: the ultraslow-spreading Gakkel and Southwest Indian Ridges[END_REF]. A local seismic experiment with seismic arrays on ice floes recorded a swarm of sea bottom explosive signals [START_REF] Schlindwein | Seismoacoustic evidence for volcanic activity on the ultraslowspreading Gakkel Ridge, Arctic Ocean[END_REF][START_REF] Schlindwein | The 1999-2001 explosive volcanic episode on Gakkel Ridge, Arctic Ocean: reconstruction of the volcanic events from multidisciplinary observations[END_REF] even 2 years after the major earthquake swarm. In 2007, the AGAVE cruise revealed thick blankets of volcaniclastic sediments in the area which contain Limu o Pele, i.e., very thin glassy crystallized rims of gas bubbles in magma and, thus, clear evidence for considerable explosive volcanism at water depths of about 4,000 m and a chain of small volcanic peaks with Oden volcano (white triangle in Fig. 1) at its easternmost location [START_REF] Sohn | Explosive volcanism on the ultraslow-spreading Gakkel Ridge, Arctic Ocean[END_REF].

While MOR volcanism is generally considered to be effusive, the relevance of explosive eruptions at MORs has only recently been recognized [START_REF] Clague | Widespread strombolian eruptions of mid-ocean ridge basalt[END_REF]. Understanding the seismovolcanic episode at the 85 • E volcanic complex at Gakkel Ridge may, therefore, yield valuable new information on deep submarine eruption processes. In this light, we intensively analyze the teleseismically recorded earthquake swarm of 1999, relocating it with an advanced classical location tool, HYPOSAT, using various velocity models and exploring new ways of describing the dataset in terms of its spatiotemporal organization.

Dataset

We used all seismic events between January and September 1999 registered in the ISC online bulletins (we downloaded them on the 6 September 2007 from http://www.isc.ac.uk) in an area from 85 • N to 87 • N and 70 • E to 95 • E. This dataset is not completely identical to previous studies [START_REF] Müller | Seismic evidence for volcanic activity discovered in Central Arctic[END_REF][START_REF] Tolstoy | Seismic character of volcanic activity at the ultraslow-spreading Gakkel Ridge[END_REF], which used a preliminary dataset available at that time. Our dataset contains 237 events in contrast to 252 events of [START_REF] Tolstoy | Seismic character of volcanic activity at the ultraslow-spreading Gakkel Ridge[END_REF] with body wave magnitudes (m b ) between 3.3 and 5 (Fig. 1).

The original dataset was located using ISCloc in the ISC standard procedures as outlined in [START_REF] Storchak | Results of locating the IASPEI GT(0-5) reference events using the standard ISC procedures[END_REF] using the Jeffreys-Bullen travel time tables and a complicated weighting scheme that delivers unrealistically small errors. The root mean square (rms) travel time residual for the original dataset is ∼1.14 s with only four of these events not fixed in depth.

For further analysis, we also inspected publicly available waveforms of the nearest stations that recorded the events (see Fig. ). Two stations from each geographical sector will be exemplarily used in the figures. Particularly, station KIV located at only 42 • E and 43 • N was only used because the avail-ability of data for some of the decisive events from stations farther east is scarce. Thus, we have relatively huge azimuthal gaps for waveform analysis and did not try a reanalysis of moment tensor solutions.

Methods and results

HYPOSAT locations of the swarm events

After intercomparison with several other location tools such as HYPO71 [START_REF] Lee | HYPO71 (revised): a computer program for determining hypocenter, magnitude and first motion pattern of local earthquakes[END_REF], HYPOCENTER [START_REF] Lienert | A computer program for locating earthquakes both locally and globally[END_REF], and ISCloc [START_REF] Luckett | ISC loc software[END_REF], we used HYPOSAT [START_REF] Schweitzer | HYPOSAT-an enhanced routine to locate seismic events[END_REF] as one of the advanced classical global location tools for relocation of the swarm:

-It is classical in the way it uses the Geiger method for locating earthquakes. -It allows the selection of confidence levels for error ellipses that we chose as 95%, the error is handled as rms values as well as distance errors in all spatial directions. -It documents the search for location solutions far away from solutions already found in previous iterations and thus documents the robustness of the hypocenter solution (better than HYPOCENTER). -It uses an arbitrary regional velocity model and can also use a global velocity model in farther distances from the source location (better than ISCloc). In our case, we used the IASPEI-91-model (IASP91; [START_REF] Kennett | Traveltimes for global earthquake location and phase identification[END_REF][START_REF] Kennett | Constraints on seismic velocities in the Earth from travel times[END_REF] or the preliminary reference earth model (PREM; [START_REF] Dziewonski | Preliminary reference Earth model[END_REF] which has been shown to work better in marine environments (Reves-Sohn, personal communication) as the global model. We varied the regional models according to Fig. 3 later in the processing. -It can process all the information available in the ISC dataset, i.e., slowness values, amplitudes, and arrival times (better than HYPO71 or HYPOELLIPSE).

The original dataset from ISC was mainly weighted according to distance to the source and for the North American sector, thus KBS, which is very similar to SPITS, and ZAL, which does not show sufficient data, were not displayed). Note the marked divergence of station residuals from 0 s in Spitsbergen and Eastern Siberia/Alaska. All stations used for waveform correlation are shown except KBS which is located too close to SPITS array travel time residual, although heterogeneities in the subsurface can disturb such a trivial setup and the regional velocity model is unknown. Velocity models from nearby locations investigated by [START_REF] Jokat | Geophysical characteristics of the ultra-slow spreading Gakkel Ridge, Arctic Ocean[END_REF] deviate strongly from IASP91 and PREM. To unfold the effects of location algorithm, weighting, local heterogeneities, and regional velocity model, we relocated the events using the IASP-91 model and no phase weighting in HYPOSAT with a starting depth of 10 km. [START_REF] Schweitzer | How can the ISC location procedures be improved?[END_REF] describes a similar procedure to improve ISC locations. The rms error of the original swarm data in HY-POSAT we obtained was ∼1.46 s (HYPOSAT baseline quality). Naturally, our aim was not only to achieve better than this relatively high values but to achieve better values than the original IS-Cloc dataset of ∼1.14 s, which naturally delivers a Velocity in km/s Fig. 3 Compressional velocity (vp) over depth plots which have been used in the hypocenter inversion to optimize the regional vp model are displayed. Where vp assesses an upper mantle value of 7.8 km/s or higher, the area of the plot is underlain by gray. In blue, five arbitrarily chosen models with various textures are shown that have been used for creating the OLD dataset described in the text. The models in red exhibit the lowest rms time residuals and are interpreted as the optimum velocity models. They adopt mantle velocities at 16 or 20 km depth, respectively. The continuous red line shows the slightly better model with a Moho at 16 km, which has been used in further elaborations in the text better value since only P arrivals are evaluated for location.

In the following, we tried to optimize the regional velocity model. Seismic events usually occur in the solid earth. In the teleseismic case, they are mostly registered on land. Thus, although the survey area around the Gakkel Ridge is entirely submarine with an average of around 4,000 m water depth, inserting a water layer into our velocity models would be erroneous because it implies that a water layer is modeled both at source and receiver location. Thus, the shallowest part of the velocity model should be an average image of the subseafloor structures.

We used classical P velocities for MOR environments [START_REF] Kearey | Global tectonics[END_REF] and local seismic velocities deduced from unreversed re-fraction profiles by [START_REF] Jokat | Geophysical characteristics of the ultra-slow spreading Gakkel Ridge, Arctic Ocean[END_REF] which showed an upper crustal velocity of around 5.6 km/s and a lower crustal velocity of 7.4 km/s in a similar geologic setting further west at Gakkel Ridge. We varied the depth of the classical layers of the MOR and the Gakkel Ridge layers and came up with a set of 124 different velocity models (62 different regional models, see Fig. 3, multiplied by two global models). Five velocity models from this set with and without sediment layers (i.e., velocities below 4 km/s) and with crustal thicknesses varying from 2 to 30 km (blue in Fig. 3) were used to create an offsetlimited input dataset (OLD) from the whole ISC bulletin readings.

The near-offset arrivals appear more sensitive to variations of the regional velocity model and are best used to optimize the regional velocity model. We checked the nearest arrivals of our event database up to a distance of 30 • for robust solution in location. Any event that delivered extremely large rms errors and/or could not find a better hypocenter solution in comparison to the HYPOSAT baseline quality and/or was not robust in space (i.e., located at a completely different spot when varying the velocity model) was deleted from the test dataset, so that only 87 events made up the OLD to optimize our velocity model.

So, we tested all 124 regional P velocity models for those 87 OLD events. The velocity model in depths higher than those shown in Fig. 3 was filled by either of the two global models mentioned above (IASP91 and PREM). Also, arrivals further away than the maximum range of validity of the regional model (variable maxdis in HY-POSAT, see Schweitzer 2001) were only modeled by IASP91 or PREM.

On top of varying the velocity model, we varied maxdis from 5 • to 30 • and the starting location depth from 5 to 50 km. All tested parameters change the location solution in a nonlinear manner. However, there is a trend that a maxdis of 20 • and a starting location depth of 10 km generally deliver the best results.

The best models from our selection provide a best average rms time residual of ∼0.90 s for the OLD and correspond to a sediment-stripped oceanic-type crust of 16 to 20 km thickness (both red in Fig. 3). We will restrict our further analysis to the slightly better model with 16 km crustal thickness. When using the complete dataset with this optimized velocity model, we decrease the rms error increases to ∼1.23 s, which is still 0.09 s higher than the ISCloc value, but 0.23 s better than the HYPOSAT baseline quality. We could not cross-check our model with gravity data because the nearest detailed surveys were conducted at 60 • E or 95 • E (Coakley and Cochran 1998), respectively.

In the original records, the epicenter locations cluster around 85 • 37 N and 84 • 40 E and the swarm is elongated along the rift axis of Gakkel Ridge. The focus of seismicity in the complete relocated set of data is situated at 85 • 35 N and 85 • E, i.e., slightly southwest of the original records and thus shifting the swarm location from the northern flank to the rift valley slightly north of the southern flank of Gakkel Ridge (most relocations in Fig. 4 show some easterly component), the elongation remains similar. The lateral error is uniformly distributed for the best events with a circular error ellipse of around 2-3 km radius for 95% confidence; for the other events, an elongation of the error ellipse along the rift valley is observed with azimuths of 110 • to 150 • (Fig. 5). Most events with a high number of reported phases in the ISC catalog (>50 phases) plot in the central valley of Gakkel Ridge in a narrow longitude band of ±1 • along 86 • E (Fig. 5).

The semimajor axis of the overall error ellipse in [START_REF] Tolstoy | Seismic character of volcanic activity at the ultraslow-spreading Gakkel Ridge[END_REF] is equal to the HYPOSAT baseline quality of 16 km. By using our best model, we relocate the earthquakes as accurately as possible, but a semimajor axis of the overall error ellipse of 10 km remains. This uncertainty is m below sea surface 2), error ellipses have been drawn (circles with ∼3 km radius for events 1 and 2 and an ellipse with semimajor axis of 9.5 km for event 3). The high-magnitude events which provide moment tensor solution are shown by upper hemisphere projections of the moment tensor. All event epicenters and moment tensor solutions are marked in colors according to their phase of occurrence (blue initial phase A, red potential vent opening phase B, green late phase C) discussed in the text. The close location in space of a blue moment tensor and one of event 2 is random and bears no significance for the rupture process encountered in event 2. The size of the event stars displays the number of phases used in the final HYPOSAT location. All events marked in red and nearly all events with more than 100 phases used in the final inversion are located in the central rift valley. A small white triangle (volcano) close to events 1 and 3 marks the position of the crater of Oden volcano which was discovered by [START_REF] Sohn | Explosive volcanism on the ultraslow-spreading Gakkel Ridge, Arctic Ocean[END_REF] unavoidable as the nearest stations to the earthquake swarm are 1,000 km away. Our geological interpretation bears this location accuracy in mind.

The hypocentral depth of the majority, i.e., 71 of 87 (81%), events lies between 10 and 60 km using the best velocity model and the OLD (Fig. 6). Using the complete dataset, the maximum depth of the majority of events decreases to 45 km (187 of 237 events, i.e., 78%). The average depth error is 32.4 km at a confidence level of 95%; however, for the OLD, it is only 20 km; for intracrustal events, it decreases to 10 km; and for events with more than 30 registered phases (115 out of 237), it drops to 5.5 km. Most events in both cases are situated slightly below the seismic Moho at 16 km depth or close to the surface. Even some of the OLD events in the upper mantle exhibit depth errors smaller than 30 km and cannot be located within the crust. Using a starting location at 20 km depth produced worse results than at 10 km, indicating that the inversion has a tendency to increase the depth of events, a smearing along the raypath that is observed in all kinds of nonlinear seismic inversion processes such as local earthquake tomography. We are unable to quantify the extents of smearing, but the documented errors of HYPOSAT indicate they are not of the order to lift the events from the upper mantle into the crust.

Station residuals

In the following, we studied the travel time residual analysis of single stations. We estimated the typical station residual by plotting a histogram for all events located in a band of 1 • latitude surrounding the focus of seismicity only allowing events with rms residual better than 1 s and a depth smaller than 45 km and binned the values into 0.5 s residual stripes (for examples, see Fig. 6).

Most frequently, we obtain a unimodal distribution and searched for the peak. Peaks, generally, scatter around 0 s residual (light blue color on Fig. 2). However, there is a trend toward big delays (>2 s) in eastern Siberia (stations TIXI and SEY) and westernmost Alaska (station ANM) and strong negative values in Spitsbergen (stations KBS/SPITS). This trend is so strong in the direction of the northernmost Atlantic that the Spitsbergen stations KBS and SPITS usually are not considered in the final solution of HYPOSAT (the number of residuals considered for the analysis of residuals at SPITS is only five of 20 plotted in Fig. 6). There have been reports about time shifts in Spitsbergen stations [START_REF] Gibbons | On the identification and documentation of timing errors: an example at the KBS station, Spitsbergen[END_REF]), however, they did not occur during the observation period and, thus, bear no relevance to our observation. Only two arrivals were picked for Jan Mayen island and no arrival at all was documented in Iceland (Vogfjord, personal communication). The heterogeneity of structures in the Barents Sea [START_REF] Ritzmann | A three-dimensional geophysical model of the crust in the Barents Sea region: model construction and basement characterization[END_REF]) might play a minor role in producing this delay. The shal-low structures in the Barents Sea are not sampled directly by seismic rays arriving at Spitsbergen, but could interfere for some events.

Last but not the least, we observe that the phases notoriously arrive earlier in the northern North Atlantic direction than predicted by theoretical models. Hence, they must propagate faster and traverse a large area of higher seismic velocities than incorporated in the regional velocity model. The opposite is true for the stations in westernmost Alaska.

Cross-correlation analysis

We reassessed the quality of the location dataset by inspecting the waveforms of the nearest stations with publicly available waveforms. This procedure helped to exclude picking errors by the reporting agencies, which we eradicated from the set of data, and allowed us to find spatially clustered events with coherent waveforms. We calculated the coherency value as in [START_REF] Maurer | Microearthquake cluster detection based on waveform similarities, with an application to the western Swiss Alps[END_REF]: similar waveforms are grouped by calculating the coherency value as ratio between cross-correlation and the square root of the autocorrelations of the single events.

In our case, we used a time window of 16 s around the first arrival (mostly the normal P wave) to include the decisive parts of the single waveforms in our analysis. The correlation window was shifted sample by sample up to a maximum of 8 s in both directions to calculate the maximum coherency in all three components as described in [START_REF] Hensch | Hypocenter migration of fluid-induced earthquake swarms in the Tjörnes Fracture Zone (North Iceland)[END_REF] and this value was thresholded at a value of 0.6 to 0.7 to group similar waveforms depending on the station. Events with a coherency higher than the threshold are considered as similar. During this initial correlation step, many well-correlated events were detected for single-station datasets. An example of the results The stations are sorted according to the filter width, which depends on sample rate, distance to the earthquake source, and quality of the data is given for station SPB4 of the SPITS array in Fig. 7, which shows good correlation even at high frequencies due to the relative proximity to the source area. Sample rates of the data used for the various stations are shown in Table 1.

Within each group, a master event was defined as the event with the least rms residual in the relocated OLD dataset and its delay was set to zero. The arrivals of other events were shifted accordingly by using the delay calculated from the shift of the time series toward each coherency maximum. We relocated the entire dataset using the improved phases but obtained no significantly different results as the percentage of corrected phases of coherent events in the total dataset is small.

If similar events can be detected on more than one station, [START_REF] Maurer | Microearthquake cluster detection based on waveform similarities, with an application to the western Swiss Alps[END_REF] pro-pose that the source process of such events has to be almost identical. Two earthquakes only produce similar waveforms at two seismic receivers if the distance between the earthquakes does not exceed a quarter of the dominating seismic wavelength and the source process is the same. A simple binary logical test was performed to find matching pairs of events. Two events were considered similar if they produced at least at two stations coherent waveforms. Only three of the 237 events (see Tables 2 and3) are so similar that we actually find matching patterns (Fig. 8). All of them occur between 5 and 9 April 1999 and HYPOSAT places them in the central rift valley of Gakkel Ridge with events 1 and 3 being located extremely close to Oden volcano [START_REF] Sohn | Explosive volcanism on the ultraslow-spreading Gakkel Ridge, Arctic Ocean[END_REF]) and to each other (see Fig. 5). The significance of the time shifts introduced by the cross-correlation can be explored in the following chapter about double-difference location where raw and corrected values of interevent times are explored (see Table 3).

We further analyzed the first motion of the P phases of events 1 to 3 in our set of equiazimuthally distributed waveform records. All three events yield a compressive first motion on all analyzed vertical components with the exception of the station INK where it is unclear whether the positive first motion can be confirmed. This strongly indicates a deviation from the classical P wave radiation pattern of a double-couple source and indicates a volumetric expansion component in the source mechanism of these events.

Double-difference location

The double-difference algorithm of [START_REF] Waldhauser | A doubledifference earthquake location algorithm: method and application to the Northern Hayward fault, California[END_REF] has been widely used to calculate relative locations of earthquakes in a Table 2 The locations of the three clustered events with similar waveforms (see Fig. 8) marked in Fig. 5 From these values, referenced double-difference times (i.e., interevent differences) have been calculated (raw) and corrected by delay times calculated from crosscorrelation (corrected) between events 1 and 2 and events 1 and 3. The stations closely analyzed in the figures have been printed in italics limited source area greatly improving the location quality. This, however, requires many well-connected coherent events with similar waveforms on several stations. Since only three of our events fulfill the similarity criterion, a relocation of the entire swarm with this method did not promise considerable improvement of the location quality. We still analyze relative travel times of the three coherent events to assess their relative location.

The double difference as defined by [START_REF] Waldhauser | A doubledifference earthquake location algorithm: method and application to the Northern Hayward fault, California[END_REF] using travel time for the first event t 1 and travel time for the second event t 2 is:

dt DD = t 1 -t 2 .
We can also express this in terms of an approximately constant time arrival t S 0 at a close station S 0 as:

t * 1 = t 1 -t S 0 ,1 and t * 2 = t 2 -t S 0 ,2 .
Assuming the events are close in space and, thus, travel time differences to the reference station are equal, we can calculate a referenced doubledifference time with respect to S 0 :

dt RDD (1, 2) = t * 1 -t * 2 = t 1 -t S 0 ,1 -t 2 -t S 0 ,2
which we calculated in Table 3. If t S 0 ,1 equals t S 0 ,2 , then the referenced double difference (RDD) will equal the real double difference. So, in Table 3, we show the referenced dt RDD for "raw" travel time values from the bulletin and for "corrected" travel times including delay times calculated from the cross-correlation. Only stations with arrivals from all three events and within 3,600 km from the source were used, with the exception of KIV that was used because of the scarcity of data in the Siberian sector (see above).

The RDD with respect to SPITS between events 1 and 2 is nearly 0 s for all the stations listed. Locating only by these stations should provide us with the same location for the two events. This not being the case shows that the emphasis during the location procedure is not on the few stations situated in close proximity to the source.

Between events 1 and 3, which are located very close in space to each other using HYPOSAT, an average dt RDD of ∼1 s remains. So, the overall HYPOSAT location was actually produced by phases further away from the source. On top, there is neither an obvious relation between distance to the source and dt RDD nor between azimuth and dt RDD , so the average dt RDD of 1 s reflects the minimum error of location. Since the rays travel mostly within the upper mantle, this is equivalent to around 8 km lateral error. The lateral error of HYPOSAT at 95% confidence level is thus an underestimation of around 25%.

All three events as laid out must reflect the same process to assume similar waveforms, all events are located between 12.6 to 16.4 km depth, all of them are located in the central rift valley, events 1 and 2 are indistinguishable up to 3,600 km offset from the source and events 1 and 3 plot nearly on the same spot close to Oden volcano [START_REF] Sohn | Explosive volcanism on the ultraslow-spreading Gakkel Ridge, Arctic Ocean[END_REF]. Thus, either this volcano or another vent aligned with Oden's chain of volcanoes and located further to the SE is the most probable source region for the repetitive events.

Magnitudes

The occurrence of high magnitudes is the most unusual feature about the earthquake swarm and is the basic reason why it was detected worldwide. Typical oceanic earthquake swarms at volcanic ridge areas normally reach their highest magnitude somewhere around a maximum magnitude of 5 or higher [START_REF] Fox | Hydroacoustic detection of volcanic activity on the Gorda Ridge, February-March 1996[END_REF][START_REF] Dziak | Evidence of a recent magma dike intrusion at the slow spreading Lucky Strike segment, Mid-Atlantic Ridge[END_REF]); however, the amount of such earthquakes is unique at Gakkel Ridge. Making a crude but empirically supported estimate [START_REF] Utsu | Relationship between magnitude scales[END_REF] by equating body wave magnitude with moment magnitude scale following [START_REF] Kanamori | The energy release in great earthquakes[END_REF]:

M 0 = 10 1.5 mb +16.1 ,
estimating a typical stress drop sigma of 1,000 bar for the average magnitude 4 earthquake and assuming a small earthquake with a circular fault shape [START_REF] Ruff | State of stress within the earth[END_REF], i.e.: sigma = 7 16 M 0 r 3 , we expect a median fault radius of 176 m for each of the events. Just using the well-located OLD events, a total of 30 km of crust in vertical and/or horizontal direction could be cracked by the events if they imaged, for example, a tectonic process or propagation of fluid-filled cracks [START_REF] Dahm | Numerical simulation of the propagation path and arrest of fluid-filled fractures in the earth[END_REF]. For such a process to occur, seismicity should cluster densely, marking one fault plane or the tip of the propagating dike. Instead, we observe that, even after relocation, seismicity remains distributed over a large area of 60 × 30 km throughout the swarm (Fig. 5), discouraging a single dyke ascent as the driving force. However, a pure tectonic origin of the swarm is also unlikely. It should produce aftershock sequences with decaying magnitudes which is not observed. Highmagnitude events (with body wave magnitude m b > 5; National Earthquake Information Center [NEIC]) occur throughout the entire swarm period (Table 4); however only 1 day, 1 February 1999 (14 days after the beginning of the swarm) is marked by three such events and shows that a large amount of energy is released in the initial phase of the swarm. The first and the third event on this day show four aftershocks, i.e., further teleseismically registered events within a radius of 10 km with decaying magnitudes. This tectonic process is overprinted by additional seismicity scattered in the area. The largest event occurs at the end of March shortly before the three similar events mentioned above. The largemagnitude events in April and later are single events without any related low-magnitude seismicity in time and space. Most of these late high-magnitude events are centered slightly north of the focus of seismicity on the northern rift flank.

Except for one event (26 May 1999), all large events have a moment tensor solution. We show the calculated percentage of double-couple component following [START_REF] Jost | A students guide to and review of moment tensors[END_REF] in Table 4. There is a first low value in mid-March 1999 and a dramatic permanent decrease to values around 80% and lower at the end of April 1999, i.e., shortly after the occurrence of the three similar events mentioned above.

Table 4 All events with higher body wave magnitude than 5, their hypocentral depth for the best location solution, and the percentage of the double-couple component in the moment tensor solution from one of three sources (NEIC: National Earthquake Information Center/USGS, Harvard: Harvard centroid moment tensor solution, GCMT: http://www.globalcmt.org) Only 11 of the 13 moment tensors are located in the region displayed in Fig. 5 3.6 Spatiotemporal organization of the swarm Figures 5 and9 show the evolution of the swarm in space and time. To be independent of single outliers in focal depth and, thus, to get a representative view of depth of hypocenters over time, we calculated the median depth of five subsequent events for the entire swarm by shifting the averaging window by three events such that an overlap of two events exists for all neighboring groups. These results are also plotted in Fig. 9. The swarm sequence starts at shallow depth. Hypocenters range across the entire crust and upper mantle in ascending and descending phases before focusing in the lower part of the crust and the uppermost part of the mantle between 10 and 20 km at day 40 and, finally, the crust-mantle boundary at day 60. In this initial phase (phase A, blue in Figs. 5 and9) of the swarm, five out of 13 seismic events with magnitude above 5 occur. These events are predominantly normal faulting earthquakes following rift valley bounding faults at the base of the crust in 14 to 19 km depth (accepting a dip of 60 • , the events would situate 8 to 10 km inside from the morphologically high rift walls-all events fulfill this criterion) with high percentage of double-couple mechanism. Laterally, the hypocenters are distributed in a 60 × 30-km long band around the median valley and show the highest depth variation, which indicates brittle behavior of the entire lithosphere/crust typical for tectonic events.

After day 70, seismicity enters a second phase (phase B). Hypocenters now follow a migration path of decreasing depth in the crust (bars marked in red and events surrounding it on Fig. 9) from the crust-mantle boundary to values below 10 km until about day 100. Apart from the odd events occurring further west, seismicity in phase B concentrates in the rift valley close to the crater of Oden volcano [START_REF] Sohn | Explosive volcanism on the ultraslow-spreading Gakkel Ridge, Arctic Ocean[END_REF]) and around the three events with very similar waveforms (encircled stars in Fig. 5) which occur at the end of this phase.

Subsequently, hypocenters can range down to 30 km before they essentially focus in a 5-km thick region above the crust-mantle boundary (Fig. 9). During the entire phase, the epicen-

Depth in km

Julian day in 1999 Fig. 9 Hypocenter depth of the events in the study area is plotted over time displayed as Julian days in the complete year 1999. No events were registered after Julian day 267, 24 September. The single hypocenter depths are marked as black stars, which are sized according to the number of phases used in the final hypocenter inversion (see scale in upper right corner). Colored horizontal bars represent the median and 80% confidence level of a group of five followup events and the time range the group is spanning. In time periods when the 80% confidence level spans the entire depth range, no horizontal bar is drawn. An overlap of two events occurs from one bar to the next. The bars are painted in colors according to their phase of occurrence as in Fig. 7a (blue initial phase A, red and gray potential vent opening phase B, green late phase C). The seismic Moho is marked by a line and the upper mantle is shaded in gray to visualize hypocenter depth in comparison to the underlying velocity model. The three events with similar waveforms are encircled ters scatter again in a larger area (Fig. 5), but the best located events occur in the uppermost mantle from 16 to 20 km depth. High-magnitude events in phase C expose the highest percentage of nondouble-couple components. Without the nondouble-couple component disappearing from the high-magnitude event, seismicity suddenly ceases with an intracrustal event on 24 September 1999 (Julian day 267) marking the end of the swarm (Fig. 9).

Interpretation and discussion

The seismic swarm forms a NW-SE-oriented elongated hypocenter cluster extending a minimum of 60 km along the rift valley of eastern Gakkel Ridge. The focus of activity is at 85 • 35 N and 85 • E where a prominent chain of volcanic cones including Oden volcano was detected by [START_REF] Sohn | Explosive volcanism on the ultraslow-spreading Gakkel Ridge, Arctic Ocean[END_REF]. Seismicity concentrates at this central volcanic ridge in the rift valley, the bounding rift valley faults to the north and south and extends over the southern rift flank. The swarm is untypical for a purely tectonic event sequence. The earthquake magnitudes mostly do not follow an Omori law because typical aftershock sequences are the exception and seismic energy release occurs in a wide range of magnitudes not only related to previous earthquakes in the series.

The widely scattered lateral distribution of events even after relocation supports a seismic activation of a large area rather than a failure of distinct faults. Concentration around a narrow area occurs after 70 days, when also a constant decrease of hypocentral depth is observed. This is typical for fluid-driven seismicity. Late high-magnitude events showing nondouble-couple moment tensors support such a scenario. On the other hand, the swarm differs from typical MOR dike intrusion events as it contains highmagnitude earthquakes and shows an unusually long duration of 8 months rather than a few weeks. [START_REF] Tolstoy | Seismic character of volcanic activity at the ultraslow-spreading Gakkel Ridge[END_REF] analyzed the original ISC locations and grouped large portions of events over time so that they found a lateral migration of seismicity over time. A migration of seismicity over time would be typical for fissure eruptions as, e.g., in the classical case study at Krafla in Iceland [START_REF] Einarsson | Seismological evidence for lateral magma intrusion during the July 1978 deflation of the Krafla volcano in NE-Iceland[END_REF]. Instead, after a detailed analysis of the error of location, seismicity remains scattered over a large area throughout most of the swarm (Fig. 5).

We, therefore, propose a different interpretation of the swarm dividing it into three phases with different characteristics:

Phase A

We believe that the initial phase of the swarm represents a tectonic activation of the ridge segment. This phase includes five of the high-magnitude events (m b > 5) of the swarm. It activates the entire crust and upper mantle in the early part, focusing to lower crustal levels after day 40. Tectonic seismicity from brittle failure is typical for upper crustal and mid-crustal material. Occasional shock-aftershock sequences and predominant double-couple source mechanisms support this idea. Not well-correlating waveforms of the seismic events throughout this phase indicate that the crustal structure is very heterogeneous and that events occur in a large subsurface volume. The spatially wide distribution of epicenters from HYPOSAT (60 × 30 × 45 km) or even the original ISCloc results are reasonable images of the distribution of seismicity (blue on Figs. 5 and9).

Phase B

After day 60, the b value increases [START_REF] Müller | Seismic evidence for volcanic activity discovered in Central Arctic[END_REF], indicating an increasing proportion of smaller events and less large events, thus imaging the decrease to lower effective stress in the subsurface. Since fluids in magma chambers [START_REF] Daines | An ultraslowspreading class of ocean ridge[END_REF] and the increase of pore pressure through interaction of magma with groundwater [START_REF] Wiemer | Variations in the magnitudefrequency distribution in two volcanic regions, Mt St Helens (Washington) and Mt. Spurr (Alaska)[END_REF] pose as buffers for effective stress, high b values are frequently produced by volcanic seismicity (red on Figs. 5 and9).

Hypocentral depth starts to decrease after 70 days and the phase culminates in three events with highly similar waveforms at seismic stations in the entire Arctic region. These must occur at the same geographic location and image the same process. First arrival analysis of these three events implies that compressive stress is propagating from the earthquake source in all direc-tions and indicates a strong volumetric component in the source as an explosion or expanding dike.

Phase B occurs just after the largest magnitude event and just before the percentage of doublecouple components of high magnitudes drop permanently. We speculate that, when the largest event occurs in March 1999, stress drop opens the pathway for material from the upper mantle to penetrate the crust. In the following, the lateral distribution of events can narrow down so much that the waveforms of different events actually mirror the same process and the same location. A typical geological process that focuses energy into a single narrow pathway is volcanism.

Thus, we speculate that magma accumulates below the crust and only when the amount of buoyancy inside this accumulating body is high enough, stress intensity at the top reaches a threshold value. A dike can form and penetrates the crust. [START_REF] Pedersen | Controlling factors on earthquake swarms associated with magmatic intrusions; constraints from Iceland[END_REF] suggest that, in cases of high magmatic stress rate, i.e., in voluminous and rapid intrusions, high seismic energy release can be triggered.

Typical volcanic vents in the survey area and within the seismicity during phase B are the cones of a volcanic chain imaged by [START_REF] Sohn | Explosive volcanism on the ultraslow-spreading Gakkel Ridge, Arctic Ocean[END_REF], which are related to explosive volcanism. [START_REF] Sohn | Explosive volcanism on the ultraslow-spreading Gakkel Ridge, Arctic Ocean[END_REF] suggest deep gas accumulation and its rapid explosive discharge following models of [START_REF] Head | Deep pyroclastic eruptions, theory and predicted landforms and deposits[END_REF] for these vents. Thus, it is likely that the seismic swarm with its high energy release could occur due to a rapid and contemporary magmatic intrusion that later delivers the magmatic ground mass for the discovered young Limu o Pele on the rim of these vents.

Seismicity during phase B rises (within 25 days) from the Moho to mid-crustal depths with an ascent velocity of 0.4 km/day or 0.46 cm/s (Fig. 9). We speculate that the shallower ascent phase to the seafloor produced only events of low magnitude and could not be imaged by globally monitorable seismicity. The red bars of median hypocentral depth in Fig. 9 assume a convex shape, indicating a sharp increase in ascent velocity which is typical for volcanic eruptions [START_REF] Battaglia | Pre-eruptive migration of earthquakes at the Piton de la Fournaise volcano (Réunion Island)[END_REF][START_REF] Rivalta | Acceleration of buoyancydriven fractures and magmatic dikes beneath the free surface[END_REF].

The fastest ascent occurs at the end from 12.5 to 8 km depth in around 1 day, i.e., an increase of ascent velocity to 7.5 cm/s, which is typical for daily swarms at slow-spreading ridges, such as offshore North Iceland [START_REF] Hensch | Hypocenter migration of fluid-induced earthquake swarms in the Tjörnes Fracture Zone (North Iceland)[END_REF]), but on the lower end of rise velocities of mantle-derived magmatic dykes between 1 cm/s and 10 m/s globally [START_REF] Rubin | Propagation of magma-filled cracks[END_REF].

Phase C

During mid-April or around day 110, the rate of seismicity decreases drastically [START_REF] Tolstoy | Seismic character of volcanic activity at the ultraslow-spreading Gakkel Ridge[END_REF]. We speculate that the area now adjusts to a changed stress field after the major phase of the volcanic eruption showing seismic events again scattered over a larger area, but mostly within a narrow band in depth. While [START_REF] Müller | Seismic evidence for volcanic activity discovered in Central Arctic[END_REF] attribute the nondouble-couple events of this phase to potential caldera-collapse phenomena, [START_REF] Tolstoy | Seismic character of volcanic activity at the ultraslow-spreading Gakkel Ridge[END_REF] point out that the creation of nondouble-couple focal mechanism occurs by distortion of focal mechanisms due to spatial heterogeneity within the source volume, for example, due to reduced seismic velocities in an area of partially molten or hot material, both of which are not excluded by our studies. The latter model would explain why earthquakes along the valley bounding faults show larger nondoublecouple events after a volcanic event. Additionally, it would fit observations that large delays for teleseismic arrivals occur in East Siberia and Alaska. They could well be the large-scale effect of the spatial change in the source volume.

However, in contrast to [START_REF] Tolstoy | Seismic character of volcanic activity at the ultraslow-spreading Gakkel Ridge[END_REF], we are confident that magnitude and nondoublecouple component are not linked with the increasing double-couple character for larger events (Table 4). Though seismicity stops suddenly at the end of September 2009 and the deep crust has potentially adjusted to the new stress situation, it seems likely that a strong eruption would leave remaining shallow heat centers or even open effusive vents for later years. Only locally installed seismometers would pick up the remaining seismic energy within the shallow crust or small seafloor explosions which have been described by [START_REF] Schlindwein | Seismoacoustic evidence for volcanic activity on the ultraslowspreading Gakkel Ridge, Arctic Ocean[END_REF].

Conclusions

The 1999 seismic swarm along Gakkel Ridge was analyzed for locations by optimizing the subsurface velocity models, calculating coherency of the waveforms of arrivals within the entire Arctic region, and verifying the double-difference behavior of earthquakes exhibiting similar waveforms. This lead to a list of overall station residuals during the entire swarm. Magnitude trends within the swarm, first arrival amplitudes, and nondoublecouple components of the moment tensor added to the following interpretation.

We find three events with almost equal waveforms and volumetric source component as crown witnesses for a volcanotectonic origin of the swarm and speculate that the location of the volcanic vent at the seafloor lies within the central rift valley at the southwestern edge or slightly outside a volcanic chain recently imaged from R/V Oden. We propose that opening of the vent occurred between 5 and 9 April 1999 when the source area narrows down laterally and an ascent of hypocenters can be observed through upper crustal material previously shattered by a tectonic earthquake sequence that occurred in the initial phase of the swarm and culminated in mid-March 1999.

Most of the seismicity is focused on the crustmantle boundary that we imaged at 16 to 20 km depth and the neighboring layers in the lower crust and upper mantle and that we suggest to be the source area for magmatic material involved in the potential eruption. Our velocity model of the source area though sparsely sampled in the upper few kilometers shows a rather crystalline crust with no large velocity gradients and thus indicating next to no sedimentation. Heterogeneity of the source area, probably related to the eruption and post-eruptive phenomena, leads to high station residuals in the direction of Spitsbergen, Eastern Siberia, and Western Alaska.

  2), taking care to realize an equiazimuthal coverage. This set of waveforms includes data from the North Atlantic sector (Kings Bay [KBS], Spitsbergen array [SPITS], Daneborg-Greenland [DAG], North Norway [ARCES], and Lovozero-European Russia [LVZ]), the North American sector (Alaska array [ILAR], Inuvik-Canada [INK], and Yellowknife array-Canada [YKA]), and the Siberian sector (Kislovodsk [KIV], Norilsk [NRIS], Peleduy [PDY], Tiksi [TIXI], Yakutsk [YAK], and Zalesovo [ZAL]

Fig. 2

 2 Fig. 2 Average station residuals constructed from histograms as in Fig. 6 for the entire Arctic region north of 60 • N. Each circular line is a 10 • latitude separation and the label surrounding the map indicates longitudes. Three sectors are discussed in the text and the following figures, which are separated by dashed lines that meet where the swarm activity is centered. Two stations of each sector will be exemplarily displayed (SPITS and DAG for the North Atlantic sector, KIV, situated outside the map-see

Fig. 4 Fig. 5

 45 Fig.4The epicenters of the seismic events of the swarm are plotted on top of the bathymetry. Original positions of the ISC location are marked by dark gray circles as in Fig.1, whereas new HYPOSAT locations are marked by black stars as in the following figures. Here, they have been outlined in white to help the eyes of the reader to distinguish between the two locations. The same events are

Fig. 6

 6 Fig. 6 Histograms of event time residuals within the six stations plotted in and mentioned in the caption of Fig. 2, i.e., SPITS and DAG for the North Atlantic sector, INK and YKA for the North American sector, and KIV and

Fig

  Fig. 7 Three groups of similar events were found for the station SPB4 of the Spitsbergen array SPITS during the station-wise cross-correlation procedure described in the text. The waveforms of the vertical component are plotted in 16 s time windows normalized to the highest amplitude

Fig. 8

 8 Fig.8When correlating waveforms at different stations, only three similar events can be detected. The vertical component of the event waveforms for the six stations also used in Fig.6are displayed in the figure in 16 s time windows normalized to the maximum amplitude. The picked first

Table 1

 1 Stations used in the waveform analysis in this study with sample rate of downloaded data and bandpass filter limits used for processing

	Station	Sample	Filter lower Filter upper
	(array if applicable) rate (Hz) limit (Hz)	limit (Hz)
	SPB4 (SPITS)	40	1	20
	INK	40	1	20
	PDY	40	1	20
	NRIS	40	1	15
	ARC2 (ARCES)	40	1	15
	YKW3 (YKA)	100	1	10
	YAK	40	1	10
	LVZ	20	1	10
	DAG	20	1	10
	IL31 (ILAR)	20	1	10
	KIV	40	1	8
	TIXI	40	1	5
	ZAL	40	1	5
	ULM	40	1	5

Table 3

 3 Arrival times in day (DD), hour (HH), minute (MIN), and seconds (SEC) and travel time difference with respect to the arrival of the P wave at SPITS are displayed for events 1-3 for those stations in a maximum distance of 3,600 km distance from the source (Offset) which registered all three events and for KIV out of consistency reasons laid out in the text

							location in latitude, longitude, and depth, the correspond-
						with the	ing axes of the error ellipses, and finally, their body wave
	final origin time as date and time in GMT, HYPOSAT	magnitude
	Event Date	Time	Latitude Longitude Error latitude Error	Depth Error depth Magnitude
				( • N)	( • E)	(km)	longitude		(km)	m b
	1	05/04/99 08:08:36 85.59	85.92	3	3	12.6	1.4	4.5
	2	06/04/99 10:05:24 85.51	86.84	3	3	14.6	1.2	4.4
	3	09/04/99 01:05:58 85.52	85.4	4	6	16.4	1.6	4.5

Acknowledgements

We thank the ISC, IRIS, NORSAR, and the Geological Survey of Canada for providing the used location and waveform data via their online access tools. We also want to thank Maya Tolstoy for providing their original dataset and Simone Cesca for providing some insight into moment tensor analysis. Last but not the least, we gratefully acknowledge the constructive feedback from H. Soosalu and an anonymous reviewer. This research was funded via the Emmy Noether programme by the Deutsche Forschungsgemeinschaft (Schl 853/1-1).


L Momayez

P Dupont

B Popescu

O Lottin

H Peerhossaini
email: hassan.peerhossaini@univ-nantes.fr

GENETIC ALGORITHM BASED CORRELATIONS FOR HEAT TRANSFER CALCULATION ON CONCAVE SURFACES

Laboratory experiments were carried out on a boundary layer over a concave wall subject to Görtler instability in the presence of forced wavelengths.

Measurements

of the local heat transfer along the concave wall permit calculation of the local Stanton number St for various axial positions x, nominal velocity n U , the wavelength λ and strength w d of the forced initial perturbations. The relation among these variables was constructed by a phenomenological argument in the form of classical power law correlations:

INTRODUCTION

Boundary layers with concave streamlines arise in several situations of practical importance such as the pressure side of gas turbine blade, combustion chamber, diffusers etc. In laminar regimes such boundary layers are susceptible to hydrodynamic centrifugal instability, generally called Görtler instability after the original work of Görtler [START_REF] Görtler | On the three-dimensional instability of laminar boundary layers on concave walls[END_REF]. At the start of the instability, the boundary layer is still a laminar boundary layer over which an array of streamwise vortices (called Görtler vortices) is superposed, as is shown on figure 1. Görtler vortices in their nonlinear stage of development drastically modify momentum and heat transfer on the wall.

The control parameter θ G , called the Görtler number, is defined as

2 / 1 R U G θ υ θ = ∞ θ
, where ∞ U , θ, and R are freestream velocity, momentum thickness and radius of curvature of the concave wall respectively, and υ is the kinematic viscosity of the fluid.

The Görtler number can be recast as

2 / 3 ) ( G αθ Λ = λ θ , where λ Λ
is the wavelength parameter, defined as

2 / 3 R R U λ υ = Λ ∞ λ
, and the vortex wavelength

α π = λ 2 .
The parameter λ Λ is constant in an experiment. with a given wavelength.

Centrifugal instabilities have a great influence on the heat transfer over curved surfaces such as gas turbine blades [2]& [START_REF] Peerhossaini | On the effects of streamwise vortices on wall heat transfer[END_REF]. The heat transfer is enhanced over a concave surface while it is inhibited over a convex wall. In the special case of gas turbine blades Görtler vortices play an important role. In fact, for calculation of thermal stresses in the turbine blade one needs to know the amount of heat transferred from the combustion products to the blade. Heat transfer correlations

existing in the open literature do not take into account the extra thermal load caused by the streamwise vortices. Therefore, they underestimate the convective heat transfer coefficient which is used as a boundary condition for the calculation of blade thermal stress loading. This coefficient is also crucial for the design of blade internal cooling. In the blade film cooling Görtler vortices can be harmful.

Duchêne [START_REF] Duchêne | Etude numérique du refroidissement par film: phénomènes de base et application à la protection thermique des aubes de turbines[END_REF] showed that the Görtler vortices exist on the pressure side of the gas turbine blades and strongly interact with the longitudinal vortices generated by the film cooling jets. In this special case he demonstrated that the Görtler instability is responsible for a slight reduction in the cooling efficiency.

In the present work we measure detailed heat transfer rates for different In previous work, Momayez et al. [START_REF] Momayez | Effects of vortex organization on heat transfer enhancement by Görtler instability[END_REF] have demonstrated the difficulties in reducing the number of governing parameters for heat transfer on a concave wall.

Those complex situations have not been expressed in terms of correlations. The problem arises either from the number of parameters [START_REF] Duchêne | Etude numérique du refroidissement par film: phénomènes de base et application à la protection thermique des aubes de turbines[END_REF] or from the unknown mechanisms that influence the heat transfer.

In this paper we derive and optimize correlations for the convective heat transfer coefficient in a boundary layer over a concave wall by using Genetic Algorithm (GA) methods.

EXPERIMENTAL SET UP AND METHODS

Experiments were carried out in the boundary layer on a concave-convex model mounted in a laminar open-loop wind tunnel described in [START_REF] Momayez | Effects of vortex organization on heat transfer enhancement by Görtler instability[END_REF]. The nominal freestream velocity could vary between 1.5 and 10 m/s with constant turbulence intensity of 0.7%.

The concave-convex model (Figure 2) has four main parts:

-the leading edge in the shape of a thick laminar airfoil (NACA-0025)

-the concave part (radius of curvature 65 cm) in which the measurements were

EFFECTS OF GÖRTLER INSTABILITY ON THE WALL HEAT TRANSFER

The cross-sectional heat advection brought about by Görtler vortices in an unstable boundary layer over a concave wall makes the heat transfer very different from that in a plain laminar boundary layer. Görtler vortices at their onset are very weak (Floryan [8], Peerhossaini and Wesferid [START_REF] Peerhossaini | On the inner structure of Görtler vortices[END_REF]). According to the linear theory, the boundary layer becomes unstable at some critical Görtler number (for instance, at θ G = 0.46 in Floryan and Saric [START_REF] Floryan | Stability of Görtler vortices in boundary layers[END_REF]), and its linear amplification (half of its growth in terms of axial distance) occurs before any disturbance can be detected.

The first effect of Görtler vortices on wall heat transfer thus appears in the nonlinear state of the vortices. This fact is evident in Figure 3 • zone AB: here heat transfer on the concave wall deviates gradually from the flat plate by the appearance of a "plateau" with roughly constant Stanton number. This heat-transfer intensification is related to the growth of the Görtler vortices under the effect of centrifugal instability (3D steady-laminar flow).

• zone BC: the heat transfer coefficient gradually reaches values close to or above turbulent boundary layer values on a flat plate. This rapid increase in heat transfer coefficient (St) is due to the secondary instability of the Görtler vortices (Toé et al. [START_REF] Toé | Etude expérimentale de l'instabilité de Görtler : Instabilité secondaire et effects des tourbillons de Görtler sur les phénomènes de transfer thermique[END_REF]). Indeed, once the Görtler vortices reach sufficient strength, subsequent longitudinal velocity gradients cause inflection points where shear instability develops (3D unsteady flow); this secondary instability grows rapidly and induces a premature (compared to the flat plate) boundary-layer transition to turbulence.

• zone CD: heat transfer ceases to increase and follows the flat-plate turbulent curve. It has been observed (Ajakh et al. [START_REF] Ajakh | Influence of forced perturbations in the stagnation region on Görtler instability[END_REF]) that turbulent spots are present in the boundary layer at these Görtler numbers. As Kestoras and Simon [START_REF] Kestoras | Combined effects of concave curvature and high free-stream turbulence intensity on boundary layer heat momentum transport[END_REF] and also Toé et al. [START_REF] Toé | Heat transfer enhancement by Görtler instability[END_REF] have noted, the heat-transfer level of the turbulent boundary layer is U , λ , R, X, ) can better correlate the results. For this purpose and given the very large number of the data collected in 44 series of measurement we have used the Genetic Algorithm method to search such correlations.

GENETIC ALGORITHM

Genetic algorithm is a branch of evolutionary computing, a swiftly growing area of reproduction intelligence, which is a search-and-optimization technique based on natural selection. Genetic algorithms (GAs) are motivated by Darwin's theory of evolution. Simply put, problems are solved by an evolutionary process that yields a best (fittest) solution (survivor) -in other words, the solution evolves. The GA is a computer simulation of such evolution in which the user provides the environment (function) in which the population must evolve.

Evolutionary computing was introduced in the 1960s by Rechenberg in his book

Evolution strategies 1 and was later developed by other researchers. GAs were created by John Holland [START_REF] Holland | Genetic algorithms[END_REF] and developed by him, his students and colleagues, culminating in Holland's 1975 book, Adaptation in Natural and Artificial Systems.

In 1992 John Koza [START_REF] Koza | Genetic Programming" III: Darwinian Invention and Problem Solving[END_REF][18] [START_REF] Koza | Evolving inventions[END_REF] used GA to develop programs to carry out certain tasks, calling his method 'genetic programming' (GP). If we are solving a problem, we are usually looking for a solution that is better than all others. The space of all reasonable solutions (the set of solutions among which the desired solution resides) is called the search space (also state space). Each point in the search space represents one possible solution. Each possible solution is distinguished by its value (or fitness) for the problem. With GA we look for the best solution, represented by one point in the search space, among a number of possible solutions. Looking for a solution is then equivalent to looking for some extreme value (minimum or maximum) in the search space. At times the search space may be well defined, but usually only a few points in it are known. The process of finding solutions in GA generates other points (possible solutions) as evolution proceeds.

In the field of acoustics, Sato et al. [START_REF] Sato | Applying genetic algorithms to the optimum design of a concert hall[END_REF] used GA to design a concert hall with optimal acoustic properties, maximizing the sound quality for the audience, the conductor, and the on-stage musicians. Porto et al. [START_REF] Porto | Alternative neural network training methods[END_REF] used evolutionary programming to train neural networks to distinguish between sonar reflections from different types of objects: man-made metal spheres, seamounts, fish and plant life, and random background noise. Tang et al. [START_REF] Tang | Genetic algorithms and their applications[END_REF] survey the uses of genetic algorithms in acoustics and signal processing.

In aerospace engineering, Obayashi et al. [START_REF] Obayashi | Multiobjective evolutionary computation for supersonic wing-shape optimization[END_REF] used a multiple-objective GA to design the wing shape for a supersonic aircraft. Williams et al [START_REF] Williams | Average and maximum revisit time trade studies for satellite constellations using a multiobjective genetic algorithm[END_REF] applied genetic algorithms to spacing satellite orbits in order to minimize coverage black-outs.

Keane and Brown [START_REF] Keane | The design of a satellite boom with enhanced vibration performance using genetic algorithm techniques[END_REF] used a GA to design a load-bearing truss or boom to be assembled in orbit and used for satellites, space stations and other aerospace construction projects.

The GA used in this study is written in C. A MATLAB ® .Implementation of the algorithm in MATLAB simplifies its use and also gives access (from the command line) to some of its main parameters in order to increase performance in particular cases. It is designed to minimize/maximize complex functions and is based on real number encoding. It can handle various characteristics (discontinuity, multimodality, and the like) without being optimized for a particular type of function and returns good results with the most common benchmark functions and with new ones as well [START_REF] Keane | The design of a satellite boom with enhanced vibration performance using genetic algorithm techniques[END_REF] [27] [START_REF] Michaelewicz | Genetic Algorithms + Data Structures = Evolution Programs[END_REF]. The algorithm integrates classic and new features, such as the ability to shape the sampling bias and a new sampling type for the BLX-crossover operator [START_REF] Tsutsui | Multi-parent recombination in genetic algorithms with search space boundary extension by mirroring[END_REF][30][31] [START_REF] Tsutsui | Search space boundary extension method in realcoded genetic algorithms[END_REF]. We recommend this GA for solving problems on which traditional methods fail to return the global minimum/maximum and especially when the function characteristics are not known.

CHOICE OF THE CONTROL PARAMETER

In a first correlation among the different problem parameters we assume that the wall temperature depends on , , , ( )

U x R λ ∞ , that is: T ( ) R , , x , U λ
d c b d c b R x a R x R G a St = = θ λ θ λ θ θ θ

Re

Coefficients a, b, c, and d are calculated for different zones (OA, AB, BC, CD) and shown in Table 1 (to save space we have shown the coefficients only for this correlation which, as we see later, fits the best the experimental results). Values of mean errors are shown in Table 2. In the third try for correlation we attempted to take into account all dimensional parameters in a minimum number of dimensionless numbers:

c b n R U a St θ υ λ =
. Again we have calculated by GA the values of a, b, c, for the different zones; the mean errors are shown in Table 2.

RESULTS

Before presenting the results we recall that the range of parameters covered in the experiments and in simulations reported in this work is given in Table 3. In Table 1 the relative magnitudes of coefficients (b, c and d) reveal the sensitivity of the heat transfer to the Reynolds number, curvature and wavelength respectively. In zone OA, b is larger than c and d, which is consistent with a nearly 2D steady laminar flow. Otherwise, in zone AB, c is larger than b and d, implying that curvature is the dominant mechanism of heat transfer sensitivity. This is in agreement with a strong 3D steady laminar flow which dominates in this zone. In zone BC, we again have a strong Reynolds-number effect with sensitivity to the wavelength. This result is in agreement with a 3D unsteady flow with transition governed by the secondary Görtler instability. Finally in zone CD, the dominant parameter is the Reynolds number with sensitivity to curvature since in this zone the flow is fully turbulent therefore, the organized vortex structure is already disappeared and the notion of wavelength has lost its significance. However, the centrifugal effects due to the curvature are still active.

Zones in which each set of constants is valid are defined in terms of the value of the non-dimensional Görtler number. Therefore, in applying the correlations these sets of constants should also be valid in the same range of Görtler number defined based on the parameters of the problem in hand..

FIGURE CAPTIONS

LIST OF TABLES

H

  derived by genetic algorithm methods. The calculated Stanton number is in good agreement with the experimental results. In this paper we propose correlations for the calculation of heat transfer on concave surfaces (such as gas turbine blade pressure side) where the flow is complex due to the presence of streamwise vortices.Key words: Görtler instability, concave boundary layer, thermal boundary layer,wall heat transfer, heat transfer intensification, genetic algorithm, correlations, height of the counter wall, H=0.15 (m)

  Görtler vortex wavelengths and two upstream perturbation strengths. The initial disturbances are generated by means of fixed manipulation grids made up of a series of forcing wires of diameter 0.18 mm or 0.80 mm and with 5, 10, 20, 30, 40, 50, and 60 mm spanwise wavelengths. The grid is placed vertically 4 mm upstream of the leading edge of the concave wall. The principal role of the perturbation (manipulator) grid is to fix the distance between vortices and their spanwise positions. The concave wall is heated with a constant heat flux layered resistive sheets. Wall temperature was measured with one hundred and ninety chromel-alumel-type thermocouples of 80 µm diameter embedded in the wall. Thermocouples are placed in such a way that the longitudinal evolution of the wall temperature can be followed as well as its spanwise variation at several longitudinal positions. The vortex wavelength of the measurements carried out in this work covers a very large range in λ Λ , between 58 to 3841, corresponding to regions in the damped zone of the linear theory, passing through the most amplified Görtler vortex ( ≈ Λ 210) and extending into the weak amplification ( ≈ Λ 3841) zone. Data used in the present work was obtained from more than 44 experimental configurations with constant wall radius of curvature.

made-

  the convex part (radius of curvature 15 cm) -the trailing edge, a flat plate that can rotate around the center of curvature of the convex section In order to minimize the longitudinal pressure gradient, a Plexiglass counter wall (opposite to the concave wall) was fixed on top of the model. The origin of the curvilinear axial coordinate x is fixed at the leading edge and the concave wall starts at x = 9 cm. Görtler vortices are generated as the result of the amplification by the centrifugal instability of upstream perturbations entering the concave boundary layer. In this study the spanwise position of the Görtler vortices was fixed by forcing predetermined wavelengths upstream of the leading edge. This was done by a manipulation (perturbation) grid made of variable-diameter wires with different wavelengths (spanwise distances λ ) placed vertically 4 mm upstream of the leading edge. In order to analyze the response of the boundary layer to the upstream perturbations, we study the reference case ( ∞ = λ , no forcing) and then the cases corresponding to increasing numbers of wires: λ = 6, 5, 4, 3, 2, 1 and 0.5 cm, with a small wire diameter of w d = 0.18 mm (weak perturbations) and a large wire diameter w d = 0.8 mm (strong perturbation). The model surface was covered with a thin (130 µm) resistance film made of a 70 µm constantan layer glued to a 60 µm Kapton film and was heated by the Joule effect. 196 chromel-alumel thermocouples of 80 µm bead measured the temperature on the Kapton side of the film. To reduce heat loss from the back of the model wall, it was insulated with phenolic foam (k = 0.02 W/mK) in which eight thermocouples were implanted. Measurements showed that heat losses by conduction from the back of the model were 6 to 8 W/m², that is, between 3 and 4% of the imposed heat flux on the wall. Thermocouples were used to measure the wall temperature along three streamwise paths and several spanwise lines. Fluid temperature in the boundary layer was measured by a thermocouple probe similar to that used by Blackwelder and Moffat (Blackwelder and Moffat [6]). Experiments were carried out at n U = 2, 3, 4.8, 7 m/s and the wall heat flux was ϕ = 200 W/m². The nominal flow velocity out of the boundary layer was measured upstream of the model's leading edge by a Pitot tube. In the present experiments pw U

  , which plots Stanton number versus Reynolds number for five nominal velocities. This figure shows the results for the reference case of n U = 3 m/s without a vortex-triggering grid. The results are also compared with Blasius heat transfer correlation for a laminar boundary layer (laminar flat-plate) and the correlation for a turbulent boundary layer on a flat plate (Kays and Crawford [11]) with constant wall heat flux. On Figure 3 we have shown the evolution of the Stanton number as a function of Reynolds number for U n =3m/s in the absence of the vortex-triggering grid (for simplicity). Note that in this Figure the experimental points follow the laminar flat-plate curve up to certain Reynolds numbers, confirming the initial very weak or zero effect of curvature on heat transfer. The measurements then diverge from the laminar flat-plate curve and after a constant-Stanton-number "plateau", heat transfer increases drastically before being saturated by turbulence. On Figure 3 laminar flat plate (Blasius) and turbulent flat plate boundary layer curves are used as reference for the evolution of heat transfer along the concave surface. Four different stages in the longitudinal evolution of the boundary layer can be distinguished in Figure4: • zone OA: this zone corresponds to the leading edge and its junction with the concave part. In this zone St follows approximately the curves corresponding to the laminar boundary layer on the flat plate (2D steady-laminar flow). Deviation of the first few points is due to flow acceleration at the leading edge.

  higher on the concave wall than on a flat plate. This heat-transfer enhancement in the turbulent regime can be explained by the mechanisms responsible for heattransfer enhancement in the laminar flow: the primary Görtler instability induces longitudinal vortices that are not completely broken down by transition and thus are still weakly present. The secondary instability induces large-scale coherent vortices in the turbulent flow that are more efficient for transport than the smaller fine grain turbulent scales. An attempt to correlate the variations of Stanton number with Görtler or Reynolds number has not given a conclusive result. In fact, Momayez et al.[5] plotted experimentally obtained values of the Stanton number versus Reynolds and Görtler number as shown on Figures 3 and 5. Instead of collapsing together the curves show a certain dispersion, implying that probably neither Gr nor Re is the proper control parameter of the problem. Question rises as to whether an other combination of the problem variables ( n

∞(,

  Görtler number based on momentum thickness) is the principal parameter and Λ is the non-dimensional wavelength in the Görtler instability. For this case we assume where a, b, c, d are the unknown coefficients determined by the GA program. The coefficients are calculated separately for each zone (OA, AB, BC and CD). Table2presents the values of mean error (%) for this correlation calculated for different nominal velocities U n and perturbation strength (wire diameter).In a second try for correlation we consider the influence of the curvature through the dimensionless number R / x , which represents the local angle of rotation of flow streamlines. This choice gives more weight to the radius of curvature of the concave wall. The influence of n U is still taken into account by the Reynolds number or the Görtler number, and we naturally nondimensionalize the wavelength with the boundary-layer thickness. Then the expected relation becomes:

From the construction point

  of view, the third correlation is obviously the worst fit, clearly because the hypothesized law lacks enough parameters to follow the complex experimental results. The two other correlations both have four control parameters and give roughly the same uncertainty. The final decision was made by putting in the three correlations the (flow and geometrical) parameters corresponding to 44 experimental configurations and comparing the Stanton-number versus Görtler-number curves obtained from the correlations with the corresponding experimental curves. The best results were obtained with the second correlation, where the effect of axial distance, relative curvature and relative wavelength are explicitly taken into account. Figures 6 and 7 compare the results of the second correlation with the corresponding experimental results. We have chosen only a few cases among the 44 optimization cases and are compared with 44 corresponding experimental cases.

4 .,

 4 CONCLUSIONSHeat transfer enhancement by streamwise Görtler vortices is considered. Passive flow manipulation for heat transfer intensification is due to the centrifugal instability caused by wall curvature. Previous results of the authors show however, that this heat transfer enhancement is sensitive to the wavelength and strength of the Görtler vortices. Therefore, perturbation grids were fixed upstream of the concave wall to manipulate and to force the wavelength of the vortices. Following a parametric study a large number of experimental data (in the form of the Stanton number as a function of the different flow and geometrical parameters of the problem) is now available. However, it has not been possible to establish a correlation between these parameters which fits the experimental data properly. Searching for such optimal correlation is a nonlinear optimization problem. There exit many methods for this optimization, we have opted for Genetic Algorithm. It can rapidly give sufficiently satisfactory solutions.GA estimator used in search for global optimal estimation of the correlation coefficients. Given the very distinct behavior of the Stanton number in a large range of the Görtler number (which covers 2D laminar, 3D laminar, transitional and turbulent regimes) optimization and evaluation of the mean values of errors were carried out for each zone separately. The results of the correlations were then compared with experimental results. This optimization was also performed for two values of perturbation strength. It was concluded that the most optimal correlation is in the form of with values of constant coefficients as given in table 1.

Figure 1 :

 1 Figure 1 : Görtler vortices in the boundary layer at x=490 mm from the leading edge ( n U = 2 m/s).

Figure 2 .

 2 Figure 2. Schematic diagram of the concave-convex model.

Figure 3 .

 3 Figure 3. Evolution of the Stanton number as a function of Reynolds number for various free-stream velocities in the absence of the vortex triggering grid

Figure 4 .

 4 Figure 4. Evolution of Stanton number versus Reynolds number for U n =3m/s in the absence of the vortex-triggering grid.

Figure 5 .

 5 Figure 5. Evolution of the Stanton number as a function of Görtler number for various free-stream velocities in the absence of the vortex triggering grid.

Figure 6 .

 6 Figure 6. Comparison between experimental (filled) and simulated (empty) evolution of Stanton number versus Görtler number for different and U n =3m/s, d w =0.18mm. Simulation results are obtained from the second correlation.

Figure 7

 7 Figure 7 Comparison between experimental (filled) and simulated (empty) evolution of Stanton number versus Görtler number for different U and λ =30mm, d w =0.18mm. Simulation results are obtained from the second correlation.

2 :Figure 2 : 7 Figure 2 .Figure 3 .

 22723 Figure 2 : Görtler vortices in the boundary layer at x=490 mm from the leading edge ( n U = 2 m/s).

Figure 4 .Figure 5 .²Figure 6 .Figure 7

 4567 Figure 4. Evolution of Stanton number versus Reynolds number for

Table 1 :

 1 Values of different coefficients for second correlation:

Table 3 :

 3 Values of different coefficients for second

	•,ο U=2m/s
	, U=3m/s
	, U=4.8m/s ♦ U=7m/s

Table 4 :

 4 Mean values of errors (%) for different experimental configurations in the present work, for first correlation

Table 3 :

 3 The of parameters covered in the experiments and in simulations reported in this work .

Rechenberg, I., / Evolutionsstrategie: Optimierung technischer Systeme und Prinzipien der biologischen Evolution /, Frommann-Holzboog,Stuttgart, 1973