Fabien Lotte

Raphaëlle N Roy

Brain-Computer Interfaces' Contributions to Neuroergonomics

Keywords: Brain-Computer Interface, Neuroergonomics, EEG

Brain-Computer Interfaces (BCIs) are systems that can translate brain activity patterns into messages or commands for an interactive application. As such the technology used to design them, and in particular to design passive BCIs which are a new means to perform mental state monitoring, can greatly benefit the neuroergonomics field. Therefore, this chapter describes the classical structure of the brain signal processing chain employed in BCIs, notably presenting the typically used preprocessing (spatial and spectral filtering, artefact removal), feature extraction and classification algorithms. It also gives examples of the use of BCI technology for neuroergonomics applications, either offline for evaluation purposes (e.g. cockpit design or stereoscopic displays' assessment), or online for adaptation purposes (e.g. video game difficulty level or air traffic controller display adaptation).

Introduction

Brain-Computer Interfaces (BCIs) are communication and control systems that enable their users to send commands and messages to a computer application by using only their brain activity, this activity being measured and processed by the system (Clerc et al., 2016a). A typical example of a BCI would be an application in which the user can move a cursor on a computer screen towards the left or towards the right, by imagining left or right hand movements respectively. While there are various ways to measure brain activity in BCIs [START_REF] Wolpaw | BCI meeting 2005workshop on signals and recording methods[END_REF], portable brain imaging techniques are typically used for practical applications. In particular, ElectroEncephaloGraphy (EEG) and functional Near InfraRed Spectroscopy (fNIRS) have been used for practical BCI applications, outside laboratories. Nonetheless, EEG remains by far the most used measure of brain activity for BCI design, both in laboratories and in real-life applications. Therefore, in the rest of this chapter, we are going to focus only on EEG-based BCIs.

BCIs can be divided into 3 categories: active, reactive and passive BCIs [START_REF] Zander | Towards passive brain-computer interfaces: applying brain-computer interface technology to human-machine systems in general[END_REF]. With an active BCI, the user voluntarily imagines some specific mental tasks (e.g., imagining left or right hand movements), whose resulting EEG patterns are translated into specific commands, e.g., moving the cursor left when the BCI recognizes an imagine left hand movement in EEG signals.

With reactive BCIs, various stimuli (often visual ones) are presented to the user, each one associated to a different command. Each stimulus is designed to evoke a different brain response (Event Related Potential -ERP -or Evoked Response EP) when the user pays attention to it. This brain response can be detected in EEG signals and thus translated into the command associated to this stimulus. The most iconic example of reactive BCIs is the P300-speller, in which the user is presented with a matrix containing all letters of the alphabet, these letters being randomly flashing [START_REF] Clerc | Brain-Computer Interfaces 2: Technology and Applications[END_REF]. The user is asked to pay attention to flashes on the letter he wants to spell (the target letter), which will give rise to a P300 ERP (a positive increase in EEG signal amplitude appearing about 300 ms after a rare and relevant stimulus) in the user's EEG signals when the target letter is flashed. No such P300 will appear when other letters are flashed. This thus enables to identify the target letter by finding out which letter evokes a P300 when flashed.

Another very widespread type of reactive BCIs are Steady-State Visual

Evoked Potential (SSVEP)-based BCIs. With such BCIs, the user is presented with various flickering visual objects (e.g., buttons on screen), each object flickering at a different frequency and being associated to a different command. If the user pays attention to one of these flickering objects, it will give rise to an SSVEP in his/her EEG signals, i.e., in an increase of occipital EEG signal power at the same frequency as the flickering frequency of the object, and at its harmonics. For instance, paying attention to a button flickering at 10Hz, would lead to an increase in 10 Hz EEG power, and possibly 20 Hz EEG power as well. Detecting the SSVEP enables to identify the object the user is paying attention to, and thus to send the corresponding command.

Finally, passive BCIs are used to monitor the user's mental states in order to adapt the application accordingly, without the user sending any volontary command through EEG signals. For instance, a passive BCI can be used to continuously estimate mental workload levels in EEG signals in order to present the user with a human-computer interface, e.g., a plane cockpit interface, that is not too cognitively difficult too use, nor too boring.

Passive BCIs are typically the kind of BCIs that can be used for neuroergonomics research and applications. Nonetheless, many of the tools developed for active and reactive BCI, and in particular EEG signal processing tools, are the same as the ones used for passive BCIs, and thus that can be used for neuroergonomics as well. Therefore, in this chapter, we present a short overview of the tools developed for BCI research that can contribute to neuroergonomics.

In particular, we will first briefly present tools to process and classify EEG signals online, in order to estimate the user mental state. Then we will show how passive EEG-based BCIs can be used for neuroergonomics and illustrate this with existing works. We then present some brief perspectives for the field.

Signal Processing

In BCI research, various signal processing tools were developed to estimate in real-time the users' mental states from their EEG signals, and this despite the noisy, non-stationary and data scarce nature of those signals. Typically, EEG signal processing in BCI follows a pattern recognition pipeline, which consists in:

1. Preprocessing EEG signals, which mostly consists in filtering them to increase their Signal-to-Noise Ratio (SNR); 2. Extracting features to describe EEG data in a compact way; 3. Classifying these features [START_REF] Lotte | A review of classification algorithms for EEG-based brain-computer interfaces[END_REF][START_REF] Bashashati | A survey of signal processing algorithms in brain-computer interfaces based on electrical brain signals[END_REF][START_REF] Makeig | Evolving signal processing for brain-computer interfaces[END_REF].

We describe below the main approaches available to perform these different steps, with a focus on approaches that can be used online. In addition, all these steps can be dynamically adapted online as well, so we briefly mention how this can be done.

Preprocessing

Preprocessing EEG signals typically consists in filtering the signal in various ways, in order to reduce the influence of artefacts such eye movements (ElectroOculoGraphy -EOG) or muscle tension (ElectroMyoGraphy -EMG) [START_REF] Fatourechi | EMG and EOG artifacts in brain computer interface systems: A survey[END_REF], and to highlight the EEG patterns representative of the mental state of interest. The most basic filtering is spectral filtering, i.e., restricting the EEG signals to some specific oscillatory components, e.g., only the alpha (8-12 Hz) and theta (4-7 Hz) rhythms to estimate mental workload.

Interestingly enough, some algorithms were developed to automatically identify the best frequency band for each subject, see, e.g., [START_REF] Pregenzer | Frequency component selection for an EEG-based brain to computer interface[END_REF].

Another essential preprocessing step is spatial filtering, i.e., combining the signals from multiple EEG channels in order to obtain a new signal with higher SNR. Different algorithms were developed to optimize spatial filters from examples of EEG data, in order to obtain EEG features that are maximally different between mental states, so as to recognize them as well as possible [START_REF] Blankertz | Optimizing spatial filters for robust EEG single-trial analysis[END_REF]. For online applications, we can notably cite the Common Spatial Patterns (CSP) algorithm to estimate mental states based on oscillatory activity [START_REF] Blankertz | Optimizing spatial filters for robust EEG single-trial analysis[END_REF], or xDAWN to estimate states based on Event Related Potentials [START_REF] Rivet | xDAWN algorithm to enhance evoked potentials: Application to brain computer interface[END_REF][START_REF] Roy | A comparison of ERP spatial filtering methods for optimal mental workload estimation[END_REF]. The Source Power COmodulation (SPOC) algorithm also enables to find spatial filters such that the power of the spatially filtered signals maximally co-vary with a continuous target variable [START_REF] Dähne | SPoC: a novel framework for relating the amplitude of neuronal oscillations to behaviorally relevant parameters[END_REF]. As such it can be used to estimate continuous mental states, such as attention or workload levels. Some extensions of such algorithms were also proposed so as to be more robust to noise or limited training data [START_REF] Samek | Divergence-based framework for common spatial patterns algorithms[END_REF][START_REF] Lotte | Signal processing approaches to minimize or suppress calibration time in oscillatory activity-based brain-computer interfaces[END_REF].

Finally, reducing the influence of EOG or EMG artefacts is also desirable.

While there are many effective algorithms to remove artefacts offline based on Independent Component Analysis (ICA) [START_REF] Urigüen | EEG artifact removal: state-ofthe-art and guidelines[END_REF], these algorithms typically cannot be used online as they are not computationally efficient enough. For online application, simpler and faster, but nonetheless useful algorithms are used. Notably, to remove EOG, regression-based algorithms based on explicit measures of EOG are often used [START_REF] Schlögl | A fully automated correction method of EOG artifacts in EEG recordings[END_REF].

For removing more general types of artefacts online, an interesting recent development is the FORCe algorithm, which combines wavelet decompositions and heuristics to remove artefactual wavelet components [START_REF] Daly | FORCe: Fully Online and automated artifact Removal for brain-Computer interfacing[END_REF].

Feature extraction

Once the EEG processed, they can be described by features [START_REF] Bashashati | A survey of signal processing algorithms in brain-computer interfaces based on electrical brain signals[END_REF]. For BCIs based on oscillatory activity, the typically extracted feature are the Band Power (BP) of the EEG signals in various frequency bands and channels. For ERP-based BCIs, the used features are typically the amplitude of the preprocessed EEG time points, for each channel, after down-sampling [START_REF] Blankertz | Singletrial analysis and classification of ERP components: a tutorial[END_REF]. These two types of features are by far the most used and give good results. It should be mentioned though that other types of features are being explored, such as complexity features, describing signal regularities, or connectivity features, quantifying how synchronized signals from different channels or frequency bands are [START_REF] Lotte | A tutorial on EEG signal-processing techniques for mentalstate recognition in brain-computer interfaces[END_REF]. While these are not as efficient as BP or time points, they can improve the overall accuracy when used in combination with them.

Classification

Classifiers learn from data which feature values correspond to which class, i.e., to which mental state here. There are multiple variants of classifiers that have been explored for BCIs, see [START_REF] Lotte | A review of classification algorithms for EEG-based brain-computer interfaces[END_REF] for a review. When it comes to online use though, only a few classifiers are typically used, the main ones being Linear Discriminant Analysis (LDA) (and its variants such as shrinkage LDA [START_REF] Blankertz | Singletrial analysis and classification of ERP components: a tutorial[END_REF] or Step-wise LDA [START_REF] Lotte | A review of classification algorithms for EEG-based brain-computer interfaces[END_REF]) and Support Vector Machine (SVM). Both classifiers are linear, and are fast to train and to use. They can also be trained from rather little training data, which makes them ideal for practical online BCI use. In recent developments, Riemannian geometry-based classifiers, which classify covariance matrices rather than vector of features, also prove very promising, including online [START_REF] Yger | Riemannian approaches in braincomputer interfaces: a review[END_REF][START_REF] Barachant | Multiclass brain-computer interface classification by riemannian geometry[END_REF].

Adaptation

As previously mentioned, EEG signals tend to be non-stationary, and the environment in which the BCI is used also leads to varying amount of external noise. As such, to reach optimal performances, it is worth considering adaptive signal processing algorithms, whose parameters are dynamically changed and optimized during online use [START_REF] Shenoy | Towards adaptive classification for BCI[END_REF]. There are a number of variants of the above mentioned algorithms that were thus designed so as to optimize online, in an incremental way, the spectral filters, spatial filters, features and classifiers, as new EEG data become available (see (Mladenovic et al., 2017, in press) for a review). While most of these algorithms remain to be tested in ecological conditions outside laboratories, they seem promising to deal with EEG fluctuations due to variation in context, noise and recording conditions that are typically encountered in real-life applications.

Contributions to Neuroergonomics

As explained earlier, since the beginning of the twentieth century the Brain-Computer Interface technology has been transposed to monitoring mental states of users. Systems that take into account information about mental states extracted from neurophysiological measures have been called biocybernetic systems or passive BCIs [START_REF] Fairclough | Fundamentals of physiological computing[END_REF][START_REF] Zander | Towards passive brain-computer interfaces: applying brain-computer interface technology to human-machine systems in general[END_REF][START_REF] Roy | Neurophysiological markers for passive braincomputer interfaces[END_REF]. Such systems can be used either offline or online, with different applicative goals and for various mental states. Definitions of mental states critical for the neuroergonomics field are given below, as well as examples of applications.

Mental states

There is a diversity of mental states that are relevant for characterizing a user/operator's state. The mental states of interest can be separated into two categories, the mental states linked to the main characteristics of the task performed by the user/operator, e.g., fatigue, and mental states linked to critical states of the system the operator interacts with, e.g., inattentional blindness.

Those two types of mental states globally generalize to the level to which the subject/operator has recruited and engaged cognitivo-attentional resources. Therefore, one can consider that the mental states linked to the main characteristics of the task relate to global resource engagement, while mental states linked to critical states of the system relate to local resource engagement. Examples of mental states that fall into these categories and that are classically estimated in neuroergonomics applications are listed in Table 1 below. For more details on each of these mental states please refer to [START_REF] Roy | Neurophysiological markers for passive braincomputer interfaces[END_REF].

It is worth noting that in addition to the mental states listed above which mainly rely on cognitive processes, emotional/affective states are of tremendous importance, as they are inseparable from cognitive states. For instance, workload can induce stress and frustration. It is therefore important to realize that mental states are never measured separately, which may be why systems trained on a particular set of data acquired in a specific setting are generally difficult to apply on another set acquired in a different setting. What's more, all these mental states interact in real life settings and thus decrease the system's performance if the latter is not conceived accordingly [START_REF] Roy | Mental fatigue and working memory load estimation: interaction and implications for EEGbased passive BCI[END_REF]. For a review on affective BCIs please refer to (Mühl et al., 2014).

Offline use: Evaluation

Passive BCIs technologies can be used offline, and in fact to this day they mostly are. More particularly, they are used for different purposes, the primary one being the evaluation of a product, a work setting or a work task, in order to determine their usability, performance, and generally their impact on the user. Hence, passive BCIs can be used for the evaluation of the comfort of stereoscopic displays [START_REF] Frey | Classifying EEG signals during stereoscopic visualization to estimate visual comfort[END_REF], for the evaluation of the difficulty of a game [START_REF] Allison | Workload assessment of computer gaming using a single-stimulus event-related potential paradigm[END_REF], a multitasking environment (Roy et al., 2016), a flying task [START_REF] Dehais | Auditory alarm misperception in the cockpit: An EEG study of inattentional deafness[END_REF] or a surgical training procedure [START_REF] Zander | Automated task load detection with electroencephalography: Towards passive brain-computer interfacing in robotic[END_REF], and also for the evaluation of prolonged and monotonous tasks such as driving [START_REF] Yeo | Can {SVM} be used for automatic {EEG} detection of drowsiness during car driving?[END_REF]. Another way of using these systems is to perform an evaluation of the user him/herself, for instance to determine his/her fitness to perform a coming task, or to determine his/her learning type. This could be promising, and to our knowledge has not been done yet.

Online use: Adaptation

Passive BCIs aim at being used online, and ideally should be so. Although to this day the scientific literature on the subject is mostly speculative, it seems to be the goal of most researchers in the field. The online use of passive BCIs allows to "close the loop" between a user and the system, and also to include the user in a more global system and regard him/her as a sub-system herself.

In order to do so the system has to adapt to the measured and inferred mental states of the user using counter-measures -if the detected state has a negative impact on performance-or more generally implicit modifications of the system.

Most studies developed and presented in the offline use section actually intend to progress towards an online evaluation of the user's mental states.

Since a few years the technological developments and the increase in variety of origin of the community members have allowed the implementation of systems that perform the measurements online, with for instance using fNIRS an online in-flight workload monitoring [START_REF] Gateau | Real-time state estimation in a flight simulator using fNIRS[END_REF], and using EEG during classical human-computer interaction tasks [START_REF] Heger | Online workload recognition from EEG data during cognitive tests and human-machine interaction[END_REF] and numerical learning tasks [START_REF] Spüler | EEG-based prediction of cognitive workload induced by arithmetic: a step towards online adaptation in numerical learning[END_REF]. Also, a few studies have recently been published with an actual adaptation of the system to the user's mental state as inferred from neurophysiological measures. For instance, workload level detection through EEG can be used to adapt the level of difficulty of a multitasking environment (i.e. the Multi-Attribute Task Battery) [START_REF] Prinzel | A closed-loop system for examining psychophysiological measures for adaptive task allocation[END_REF] and more recently that of a Tetris video game [START_REF] Ewing | Evaluation of an adaptive game that uses EEG measures validated during the design process as inputs to a biocybernetic loop[END_REF] and of an air-traffic controller display [START_REF] Aricò | Adaptive automation triggered by EEG-based mental workload index: a passive brain-computer interface application in realistic air traffic control environment[END_REF]. Additionally, passive BCIs can also be used to improve active BCIs, i.e., in which the user controls an effector. This was demonstrated by Chavarriaga and collaborators who used the EEG responses to errors committed by the active BCI (a.k.a., Error Related Potentials) to adapt the whole system and increase its performance [START_REF] Chavarriaga | Learning from EEG error-related potentials in noninvasive brain-computer interfaces[END_REF].

Perspectives

Passive brain-computer interfaces offer a very promising means to achieve online objective mental state monitoring of users and operators. Therefore it is of great interest to pursue their development for neuroergonomics applications.

Even though the literature on passive BCIs has increased drastically these last few years most research is still conducted in laboratory and not in ecological settings. Yet portable recording devices that are quite robust to environmental noise have been released (e.g. dry EEG systems [START_REF] Nijboer | Usability of three electroencephalogram headsets for brain-computer interfaces: A within subject comparison[END_REF]). Therefore, this may be due to the lack of neural features and learning algorithms robust to changes in tasks, settings and subjects. Research should therefore focus on these matters, as well as try and develop systems that actually work online but also adapt both at the signal processing level and at the interface level.

  to time-on-task: fatigue, vigilance, boredom and mind wandering; Linked to mental workload: load in working memory, divided attention, social or temporal stress. Critical states of the system Inattentional blindness or deafness phenomena, automation surprise/confusion.

Table 1 :

 1 Mental states classically estimated in neuroergonomics applications and generated in response to either main characteristics of the task or critical states of the system.

surgery. J Med Robot Res, (pp. 1-10).


domain: phys.meca.mefl phys spi

text:
X Nicolas

P Traore

A Mojtabi

J P Caltagirone

AUGMENTED LAGRANGIAN METHOD AND OPEN BOUNDARY CONDITIONS IN 2D SIMULATION OF POISEUILLE±BE Â NARD CHANNEL FLOW

Keywords: open boundary conditions, Poiseuille±Be Ânard ¯ow, augmented Lagrangian, thermal instabilities

The main objective of this study is to compare the in¯uence of different boundary conditions upon the incompressible Poiseuille±Be Ânard channel ¯ow (PBCF) in a 2D rectangular duct heated from below. In a ®rst technical part the algorithm used to carry out this work, based on the augmented Lagrangian method, is presented. The implementation details of the ®ve different open boundary conditions (OBCs) and the periodic boundary conditions (PBCs) tested in the present paper are also given. The study is then carried out for 1800 ` Ra 4 10,000, 0 ` Re 4 10 and 0Á67 4 Pr 4 6Á4. The ®ve selected OBCs, applied at the outlet of the computational domain, respectively express the following conditions: a square pro®le for the velocity (OBC1), mass conservation (OBC2), zero second derivative of the horizontal velocity component (OBC3), a mixed boundary condition combining Dirichlet and Neumann conditions (OBC4) and an Orlanski-type boundary condition (OBC5). A good estimation of the perturbation amplitude and of the length of the perturbed zone at the outlet boundary is proposed. It is shown that OBC5 causes very little perturbation in the recirculating ¯ow compared with the other OBCs.

INTRODUCTION

Often in computational ¯uid dynamics (CFD), to compute an out¯ow, one encounters the problem of open boundary conditions (OBCs). Computational time costs and memory limitations do not permit simulations on large physical domains. To solve the problem mathematically, the domain must be truncated and con®ned between arti®cial boundaries onto which numerical conditions are introduced, depending on the nature of the governing equations. However, it is dif®cult to establish these boundary conditions when the studied ¯uid system is isolated from the effects of its environment; information about the dynamics of the ¯uid through the boundary is all the more necessary as the ¯ow is of an elliptic nature (i.e. is made of eddies and recirculations). The OBCs are the result of assumptions which do not re¯ect exactly the real physical phenomenon; if they are not properly chosen, they can have spurious effects such as numerical instabilities, re¯ections at the boundary and, generally, errors that can propagate and alter the results throughout all the computational domain. Thus the choice of good OBCs is crucial in CFD. This choice is closely linked to the physical behaviour of the ¯ow and to the method used to solve the problem.

The main objective of this paper is to study the in¯uence of ®ve different OBCs, at the outlet of the computational domain, upon the thermoconvective structures that are liable to develop when simulating incompressible Poiseuille±Be Ânard channel ¯ow (PBCF). The PBCF is a mixed convection ¯ow in a horizontal rectangular channel heated from below. If there is no heating and the channel is open, it is a simple Poiseuille ¯ow. If vertical walls close each end of the channel and if the vertical temperature gradient is suf®ciently high (Ra b 1708, where Ra is the Rayleigh number), it is a thermoconvective Be Ânard ¯ow. For three-dimensional channels, when the two phenomena are added, the ¯ow structure is quite complex. For a given Prandtl number Pr and a given aspect ratio of the duct, several types of thermoconvective structures, depending on Ra and Re (where Re is the Reynolds number), can be observed: there can be either transversal rolls, longitudinal rolls (i.e. thermoconvective rolls with their axes either orthogonal or parallel to the axis of the channel) or intermittent patterns.

The stability of these different con®gurations is a very interesting problem; this has given rise to many theoretical 1±7 and experimental 1,4,8±11 papers. To our knowledge, only two similar papers [START_REF] Schro Èder | Three-dimensional convection in rectangular domains with horizontal through¯ow[END_REF][START_REF] Chen | Laminar, buoyancy induced ¯ow structures in a bottom heated, aspect ratio 2 duct with through¯ow[END_REF] have recently presented results obtained by means of direct three-dimensional numerical simulation. Other works, 14±16 carried out in the frame of CVD (chemical vapour deposition), have mainly focused on heat transfer enhancement related to different thermoconvective structures. When studying their stability numerically, several dif®culties appear: it is necessary to ensure, ®rst, that the amplitude of the perturbation caused by the OBC does not cause bifurcations between the different types of thermoconvective structures (transversal or longitudinal rolls) and, second, that the length of the perturbed zone above the outlet is limited compared with the length of the computational domain.

In the present paper the analysis of the in¯uence of OBCs is limited to the two-dimensional PBCF. In this case the ¯ow structure consists of two-dimensional contra-rotative rolls with their axes transversal to the direction of the average channel ¯ow and moving away with the mean ¯ow. The purpose of this study is to determine the OBC which causes the smallest perturbation among the ®ve selected ones. With this aim in view, the ®ve OBCs are all compared with a benchmark solution of the PBCF that was computed by Evans and Paolucci 17 for Re 10, Ra 10,000 and Pr 2 3 . They are also compared with simulations involving periodic boundary conditions (PBCs). The Nusselt number, the variation of the transversal roll wavelength l and the time period t are presented for the different OBCs. Several criteria are also proposed for estimating the length of the perturbed zone and the amplitude of the perturbation; their evolution as a function of the Rayleigh and Reynolds numbers is studied for the following range of dimensionless parameters: 0 `Re 4 10, 1800 `Ra 4 10,000 and 0Á67 4 Pr 4 6Á4.

The 2D PBCF was already proposed by Sani and Gresho 18 as a test case for an OBC minisymposium that took place in Swansea, U.K. on 10 July 1989. Bottaro [START_REF] Bottaro | Note on open boundary conditions for elliptic ¯ows[END_REF] and Kobayashi et al. [START_REF] Kobayashi | Comparison of several open boundary numerical treatments for laminar recirculating ¯ows[END_REF] compared different OBCs with Evans and Paolucci's benchmark solution. Other test cases were proposed during the minisymposium; the summary and remarks of this event are described in Reference 18.

As mentioned before, the choice of the OBCs is partly linked to the method used to compute the solution of the problem. For the present work a ®nite volume numerical code based on the augmented Lagrangian method has been developed. This method consists of an optimization technique to determine a velocity±pressure saddle point under the incompressibility constraint H ? V 0; the saddle point is computed by an iterative algorithm of Uzawa type. [START_REF] Uzawa | Iterative methods for concave programming[END_REF][START_REF] Temam | Navier±Stokes equationsÐTheory and Numerical Analysis[END_REF] This method at present is not frequently used. However, our experience in CFD has shown that it is very ef®cient in comparison with the classical numerical schemes: it permits us to simulate ¯ows with important constraints or with strong variations in the control parameters and in the geometry, for which other methods diverge. The augmented Lagrangian method was formulated in the publications by Fortin and Glowinski [START_REF] Fortin | Me Âthodes de Lagrangien Augmente Â[END_REF] and Glowinski, [START_REF] Glowinski | Numerical Methods for Nonlinear Variational Problems[END_REF] in which a variety of problems in the ®elds of ¯uid mechanics and elasticity are treated. The ®rst part of this paper presents the main outlines of the method. For more details see References 23±26, in which numerous remarks on the implementation of the algorithm are available.

NUMERICAL PROCEDURE

Governing equations

The problem is formulated considering the usual hypotheses for a Newtonian and incompressible ¯uid and the Boussinesq approximation is assumed to be valid. The three conservation equations (continuity, momentum and energy equations), provided with adequate initial and boundary conditions, enable us to express the solution of the PBCF. In primitive variables (velocity V, pressure P, temperature T) and under its dimensionless form the system that has to be solved takes the following expression:

I H ? V 0Y 1 dV dt V ? HV ÀHP 1 Re DV Ra Re 2 Pr TkY 2 dT dt V ? HT 1 RePr DTX 3 V b b b b b b b b b b X
The reference length, velocity, temperature and pressure for scaling are the height H of the duct, the average velocity V o of the PBCF, the temperature difference T hot 7 T cold between the bottom and the top of the channel and r o V o (where r o is the mass per unit volume of the ¯uid) respectively. Thus Re V o Han is the Reynolds number, Ra gb(T hot 7 T cold )H 3 ana is the Rayleigh number and Pr naa is the Prandtl number. Here n, b and a are the kinematic viscosity, thermal expansion coef®cient and thermal diffusion coef®cient of the ¯uid respectively. g is the gravity constant and k is the vertical unit vector. The time scheme used to solve (I) is Gear's second-order backward implicit scheme ((a, b, c) (1Á5, 2, À 0Á5) in system (II) below). The time discretization being uniform, the time step is denoted Dt; the current time and the following time step (time at which the ®elds are unknown) are nDt and (n 1)Dt respectively. The superscripts of the primitive variables correspond to the time index. In the momentum equation the buoyancy term RaTkaRe 2 Pr is treated explicitly; therefore the two coupled equations ( 2) and ( 3) are solved separately. The advective term (V ? HV being linearized, the time-discretized form of (I) consists of computing (V n1 Y P n1 Y T n1 ) as follows:

II H ? V n1 0Y 4 A n V n1 HP n1 FY 5 BT n1 GY 6 V b b b b X
where

A n a Dt V n ? H À 1 Re DY F Ra Re 2 Pr T n k bV n cV nÀ1 Dt Y B a Dt V n1 ? H À 1 RePr DY G bT n cT nÀ1 Dt X

Augmented Lagrangian algorithm

In system (II) the main dif®culties occur in the computation of equations ( 4) and ( 5); these dif®culties are due to the velocity±pressure coupling, the treatment of the continuity constraint H ? V 0 and the determination of boundary conditions on pressure. It could be argued that the coupled linear system in (V, P) corresponding to equations ( 4) and ( 5) can be solved directly; however, the associated matrix is often very large and very ill-conditioned, making it dif®cult to ensure the stability and convergence of the iterative solver for the linear system. It is often preferred to separate the computation of the velocity from that of the pressure. Here this is achieved by using the augmented Lagrangian method. [START_REF] Fortin | Me Âthodes de Lagrangien Augmente Â[END_REF][START_REF] Glowinski | Numerical Methods for Nonlinear Variational Problems[END_REF] It has been shown [START_REF] Fortin | Me Âthodes de Lagrangien Augmente Â[END_REF][START_REF] Glowinski | Numerical Methods for Nonlinear Variational Problems[END_REF] that the variational formulation of the problem `to ®nd (V, P) in order that (4) and ( 5) are satis®ed' consists of minimizing a functional (derived from the momentum equation) under the incompressibility constraint H ? V 0. This minimization problem is equivalent to solving a saddle point search problem for (V, P) after having introduced the augmented Lagrangian associated with the functional and the constraint de®ned above. The determination of the saddle point (V, P) is obtained by the iterative Uzawa algorithm. [START_REF] Uzawa | Iterative methods for concave programming[END_REF][START_REF] Temam | Navier±Stokes equationsÐTheory and Numerical Analysis[END_REF] (V n , P n ) being known, (V n1 , P n1 ) are obtained using the following external loop with iteration index k (time index is n):

III V k0 Y P k0 V n Y P n Y A k r V k1 F n À HP k Y 7 P k1 P k À sH ? V k1 Y 8 when kH ? V k1 k `eY then V n1 Y P n1 V k1 Y P k1 X V b b b b b b b b b b X
In equation ( 7), A k r A k À rHH ? . The additional operator ÀrHH ? comes from the augmented Lagrangian and takes into account the incompressibility constraint; its function is to increase the rate of convergence of algorithm (III). Here r, s and e are three positive constants; e has a small value and is used to determine the convergence criterion (in this paper all the results were obtained with e 10 ±5 ). kH ? V k1 k represents the arithmetic mean of the absolute values of H ? V at each mesh point.

Fortin and Glowinski [START_REF] Fortin | Me Âthodes de Lagrangien Augmente Â[END_REF] show that the convergence condition for algorithm (III) is 0 `s `2(r 1aRe). They also prove that a greater rate of convergence of (III) is achieved as r increases; on the other hand, the bigger r is, the more ill-conditioned the matrix of the linear system (derived from A k r in equation ( 7)) will be. Thus the choice of r must be a compromise in relation to two con¯icting requirements: minimizing the number of iterations on k and minimizing the number of iterations necessary to solve the linear system. Note that when very small time steps are used (e.g. for unsteady problems), the matrix of the system becomes more diagonally dominant and it is possible to increase the value of r. Fortin and Glowinski experimentally show, [START_REF] Fortin | Me Âthodes de Lagrangien Augmente Â[END_REF] by several numerical tests, that the optimum rate of convergence is approximately obtained when s is slightly bigger than r, when r is large. Thus, in this research, after having tested several values of r and s between 10 and 5000, we chose r 400 and s 500 when Dt 0Á0005. The complete algorithm, to move from time nDt to time (n 1)Dt, is presented in Figure 1.

The augmented Lagrangian term ÀrHH ? V k1 couples the velocity components and subsequently adds eight diagonals (in 2D problems) in the matrix associated with A k r . Therefore a highly ef®cient method has to be used to solve the linear system derived from (7), since U and W must be calculated at the same time. Since the discretization of the advective term V ? HV makes the linear system matrix non-symmetric and since it is a regular matrix, the Bi-CGSTAB (biconjugate stabilized gradient) method [START_REF] Van Der | Bi-CGSTAB: a fast and smoothly converging variant of Bi-CG for the solution of nonsymmetric linear systems[END_REF] with a preconditioning based on a modi®ed and incomplete Gauss factorization MILU [START_REF] Gustafsson | On ®rst and second order symmetric factorization methods for the solution of elliptic difference equations[END_REF] was the chosen solver. The same solver is used to compute T from the linear system derived from (6).

The space discretization of the momentum equation ( 7) and energy equation ( 6) is generated by using a ®nite control volume method on a staggered grid. [START_REF] Patankar | Numerical Heat Transfer and Fluid Flow[END_REF] This grid is Cartesian and uniform in each space direction, with the space step Dx in the horizontal direction and Dz in the vertical direction. P and T are computed at the nodes of the control volumes; U and W, the horizontal and vertical velocity components, are computed at the centre of the horizontal (resp. vertical) edges of the control volumes (see Figure 2). A central difference scheme is used for the space discretization of the convective ¯uxes, while the diffusive ¯uxes are discretized with central second-order derivative approximations. We note that the ®nite volume method leads to the following discretized equation for each control volume: 29

a P ff P a E ff E a N ff N a W ff W a S ff S bfY 9
where f stands for one of the three variables T, U or W and a i (f) (with i P, E, N, W, S) and b(f) are coef®cients depending on the discretization scheme. The subscript P refers to the central control volume point (i, k) and subscripts E, N, W and S refer to its neighbouring points to the east, north, west and south respectively (see Figure 2). To be able to compute the linear system corresponding to ( 9) for all the control volumes, equation ( 9) must be veri®ed everywhere, even at the outlet boundary. However, this is not always possible when, for example, f E is unknown on this boundary. To solve the problem, OBCs are implemented. In this paper, two formulations have been used. The `weak formulation' consists of expressing f E from the discretized form of the OBC and introducing it in equation ( 9). For instance, in the case of 9), the coef®cients are modi®ed in the following way:

a W f a W f a E fY a E f 0Y a P fY a N fY a S f

and bf remain unchangedX

The `strong formulation' consists of directly identifying the discretized form of the OBC with equation (9). For instance, in the case of Figure 2(b) the Neumann condition gives f P f E . Then the coef®cients are

a P f a E f 1Y a N f a W f a S f bf 0X
So far, no hypothesis has been made about the nature and geometry of the ¯ow; it has only been assumed that the system was provided with adequate initial and boundary conditions in order to be solved. The conditions that have been implemented to compute the PBCF are presented in the following subsection. Note that one of the advantages of Uzawa's algorithm (III) (compared e.g. with the classical SIMPLE, SIMPLER, PISO, etc. methods) is that no boundary condition on pressure is required. Indeed, owing to the explicit treatment of P in (8), only the initial conditions for pressure (and the entire velocity ®eld) are needed to compute P at each time step.

Computational con®gurations and boundary conditions

Two main computational con®gurations have been used to compute the PBCF. The ®rst one is a rectangular channel with one of the ®ve tested OBCs at the outlet and a conductive Poiseuille ¯ow at the inlet; the second one is a rectangular channel with PBCs. The description of these con®gurations and the method used to implement the different OBCs or the PBCs are given below.

First con®guration. This is illustrated in Figure 3 and its initial and boundary conditions for (V, P, T) are given in Table I (except for the ®ve OBCs, whose de®nition and implementation are given in the next subsection).

This con®guration allows us to observe the space ampli®cation of the perturbation until non-linear saturation occurs. [START_REF] Mu Èller | Convective patterns in horizontal ¯ow[END_REF] When thermoconvection develops in the PBF, three zones can be distinguished (see Figure 3): (i) for 0 4 x 4 x in , the inlet zone in which the perturbation is growing; then, after its saturation, (ii) for x in 4 x 4 x out , a fully established periodic ¯ow of transversal rolls; and (iii) near Figure 3. First computational con®guration with conductive Poiseuille ¯ow at inlet and one of ®ve tested OBCs at outlet the outlet, for x out 4 x 4 LaH, the zone on which we are going to focus, where the rolls are more or less distorted by the OBC (the longer this zone is, the more spurious the effect of the OBC).

Note that the length of the inlet and of the outlet can vary according to the values of Ra, Re and Pr. Therefore computational domains with different aspect ratios LaH (where L and H are the length and height of the duct respectively) are used to obtain a fully developed periodic ¯ow for each computed case. For the test case proposed by Evans and Paolucci 17 (Ra 10,000, Re 10, Pr 2

3 ), x in % 2; therefore we chose LaH 5 and for the space and time steps we took (Dx, Dz, Dt) (0Á0725, 0Á05, 0Á0005). For Pr 6Á4, when (Ra, Re) (4700, 0Á18) and (4700, 3), we chose LaH 10; but, when (Ra, Re) (1804, 0Á18), we took LaH 20, because x in % 9. In all cases with Pr 6Á4, (Dx, Dz, Dt) (0Á1, 0Á05, 0Á0005) unless stated otherwise.

OBCs at outlet of ®rst con®guration. Many OBCs are used by the CFD research community. Frequently their physical interpretation is dif®cult, but they permit one to close a problem mathematically. When simulating a ¯ow, no clear criterion permits one to select one of these OBCs over the others; the only way to proceed is to analyse the behaviour of the ¯ow a posteriori. Subsequently we decided to test ®ve of the most commonly encountered OBCs. Their de®nitions and the ways in which they are implemented are given below. Note that the strong formulation is used to implement the OBCs for U, while the weak formulation is used for W (except for OBC1) and T.

OBC1 U V o Y W 0Y dT dx 0Y
where

V o z1 z0
À6z 2 À zdz is the average velocity of the ¯ow. To be more precise, in order to take into account both the no-slip condition on the horizontal plates of the channel and the mass conservation, we took U(i, 0) U(i, Nz) 0 and U(i, k) [Nza(Nz 7 1)]V o for 0 `k `Nz (Nz is the total number of meshes in the z-direction).

OBC1 is a coarse way to simulate the porous side walls [START_REF] Luijkx | On the existence of thermoconvective rolls, transverse to a superimposed mean Poiseuille ¯ow[END_REF][START_REF] Schro Èder | Three-dimensional convection in rectangular domains with horizontal through¯ow[END_REF] that are often placed at the inlet and outlet of experimental channels in order to avoid direct contact of the measurement zone with the feed pump. The discretization of the condition for U gives (see Figure 2(b))

U k1 Nx 1Y k À U k1 NxY k ÀDxaDzW n NxY k 1 À W n NxY kX
In the strong formulation this expression implies

a P U a W U 1Y bU ÀDxaDzW n NxY k 1 À W n NxY kY a E U a N U a S U 0X
OBC2 is directly obtained by applying the conservation of mass at the outlet boundary. We have veri®ed (not shown here) that OBC2 has exactly the same in¯uence on the numerical solution (same pro®les of U, W and T at the outlet, same amplitude of perturbations, same perturbed length, etc.) as a Neumann boundary condition on U, W and T. The only difference is that the computational time is shorter with OBC2. Indeed, with the augmented Lagrangian algorithm, a larger number of loops in k (see system (III) or Figure 1) is necessary to reach the convergence criterion kH ? Vk `e when using the simple Neumann condition. This criterion is more easily veri®ed by OBC2 since it ensures the continuity equation on the boundary. OBC2 must be seen here as a means to close the problem mathematically.

OBC3

d 2 U dx 2 0Y dW dx 0Y dT dx 0X
U k1 Nx 1Y k is obtained explicitly by discretizing the condition for U with a second-order upwind scheme in space:

U k1 Nx 1Y k 2Á5U n NxY k À 2U n Nx À 1Y k 0Á5U n Nx À 2Y kX Therefore a P U 1Y bU 2Á5U n NxY k À 2U n Nx À 1Y k 0Á5U n Nx À 2Y kY a E U a W U a N U a S U 0X
Note that the semi-implicit treatment of this OBC implies an ill-conditioning of the linear system matrix: a P U 1 and a W U 2Á5. OBC3 is less restrictive than the Neumann boundary condition: whereas the latter implies a constant extrapolation of U at the outlet boundary, OBC3 can be seen as a linear extrapolation of U on this boundary. Note that a second-order derivative was also tested for T, but the results were never as good as those obtained with the ®rst-order derivative.

OBC4 dU dx Àb U U À U I Y dW dx Àb W W À W I Y dT dx Àb T T À T I X
For U the discretization gives (see Figure 2(b))

U k1 Nx 1Y k À U k1 NxY k ÀDxb U fU k1 Nx 1Y k U k1 NxY ka2 À U I kg
and the strong formulation implies

a P U 1 Dxb U a2Y a W U 1 À Dxb U a2Y bU Dxb U U I kY a E U a N U a S U 0X
For W the discretization gives (see Figure 2(c))

W k1 Nx 1Y k À W k1 Nx À 1Y k À2Dxb W W k1 NxY k À W I k
and the weak formulation implies

a W W a W W a E WY a P W a P W 2Dxb W a E WY bW bW 2Dxb W W I ka E WY a E W 0Y a N W
OBC5 dU dt V o dU dx 0Y dW dt V o dW dx 0Y dT dt V o dT dx 0X
For U the discretization gives

U k1 Nx 1Y k U k1 NxY k À U n Nx 1Y k À U n NxY k 2Dt V o Dx U k1 Nx 1Y k À U k1 NxY k 0Y
which implies

a P U V o Dxa2DtY a W U V o À Dxa2DtY bU U n Nx 1Y k U n NxY kDxa2DtY a E U a N U a S U 0X
For W the discretization gives

W k1 NxY k À W n NxY k Dt V o 2Dx W k1 Nx 1Y k À W k1 Nx À 1Y k 0Y which implies a W W a W W a E WY a P W a P W 2DxaV o Dta E WY bW bW 2DxaV o DtW n NxY ka E WY a E W 0Y a N W and a S W remain unchangedX
This boundary condition is a derived and simpli®ed form of Orlanski's boundary condition. [START_REF] Orlanski | A simple boundary condition for unbounded hyperbolic ¯ows[END_REF] It was initially introduced by Bottaro [START_REF] Bottaro | Note on open boundary conditions for elliptic ¯ows[END_REF] and later tested by Kobayashi et al. [START_REF] Kobayashi | Comparison of several open boundary numerical treatments for laminar recirculating ¯ows[END_REF] on the same con®guration as the one used by us. OBC5 can be seen as an approximation of the momentum equation on the boundary; U, W and T are convected by a transport equation whose celerity is the average velocity V o of the ¯ow (V o being de®ned above). In this way the re¯ection of these quantities when passing through the boundary is expected to be avoided.

Second con®guration with PBCs (periodic boundary conditions).

To try to free ourselves from the problem of OBCs, we have computed the PBCF with PBCs since this ¯ow is spatially and temporally periodic when it is fully developed. When implementing PBCs, it is no longer possible to induce the ¯ow in the channel by enforcing a Poiseuille ¯ow at the inlet; now it is induced by imposing a constant pressure gradient DP o between the inlet and outlet boundaries during all the computation. Taking DP o as the new scaling factor for pressure (instead of r o V o ), the boundary conditions for pressure are P 1 at the inlet and P 0 at the outlet. Taking (DP o /r o ) 1a2 as the new scaling factor for velocity (instead of V o ), the Reynolds number in system (I) is now de®ned by Re (DP o ar o ) 1a2 Han.

The periodicity for V and T is assumed in the following way: at each time step (n 1)Dt, (V n , T n ) obtained at time step nDt in the centre of the duct (x La2H) is simultaneously enforced as a boundary condition at the inlet and outlet of the channel (see Figure 4); in this way, at any moment, Figure 4. Second computational con®gurationÐsimulation of periodic boundary conditions the convective roll which goes out from the numerical domain goes into it on the opposite side. The top and bottom boundary conditions as well as the initial conditions are the same as those of the ®rst con®guration (see Table I). Note that because of the special way of proceeding, the number of transversal rolls in the computational domain is necessarily a multiple of four. Therefore the transversal roll wavelength l is imposed by the aspect ratio LaH. In Evans and Paolucci's benchmark case, [START_REF] Evans | The thermoconvective instability of a plane Poiseuille ¯ow heated from below: a proposed benchmark solution for open boundary ¯ows[END_REF] characterized by Ra 10,000, Re 10 and Pr 2 3 , it was shown that l % 1Á44. Thus for the simulation with PBCs a channel with LaH 2Á88 has been chosen so that four transversal rolls with the right wavelength appear in the duct; in this case the steps in space are (Dx, Dz) (0Á0725, 0Á05).

RESULTS AND DISCUSSION

The results obtained with the ®ve OBCs are compared with a reference solution for each set of Rayleigh, Reynolds and Prandtl numbers. These solutions are called `Reference' in all the ®gures and tables. They are computed over domains twice as long as those used to test the OBCs and only their ®rst half is taken into account. Indeed, the perturbation due to the OBC at the outlet of the channel is supposed to be negligible on this ®rst half. For the Reference the grid is ®ner and the OBC used is OBC5.

In the ®rst step all the tests are carried out for the case proposed by Evans and Paolucci. The Reference is computed on a domain of aspect ratio 10 with (Dx, Dz) (0Á05, 0Á0385). Figure 5 permits us to compare the streamlines, isotherms, iso-values of horizontal velocity component and Figure 6. Field pro®les at outlet boundary at time t* for different OBCs and PBCs iso-vorticity contours obtained with the different OBCs and with the Reference. Each solution is recorded at the same time t*, de®ned as the time when the temperature reaches a local minimum at the ®xed position xY z L À [START_REF] Mu Èller | Convective patterns in horizontal ¯ow[END_REF] 4 lY Ha2, where l is the roll wavelength (l 1Á44). From the OBC point of view this is the most unfavourable time, since this is the moment when a thermoconvective roll sits astride the outlet boundary: nearly half of this boundary is submitted to an entrance ¯ow. This is also the most interesting time to evaluate the behaviour of the OBCs faced with return ¯ows. Figure 5 shows that the solution is very perturbed by OBC1, OBC2, OBC3 and OBC4, whereas OBC5 seems to allow a good undisturbed evacuation of the transversal rolls.

To focus on the perturbation at the outlet boundary, Figure 6 presents the c, T, U and W vertical pro®les at x LaH for each OBC and for the PBCs compared with the Reference (where c is the streamfunction). The recording time is t* as before. The pro®les which are furthest from the reference solution are obtained with OBC2; the behaviours of OBC1, OBC3 and OBC4 are not very good either. On the other hand, OBC5 behaves almost like the Reference (except for W). OBC5 is also the i and k are the horizontal and the vertical space index respectively, Nt is the total number of time iterations and Nz is the total number of meshes in the z-direction. In this way the envelopes that are obtained illustrate clearly the magnitude of the perturbation and the length of the perturbed zone for each OBC. The amplitude of the perturbation varies signi®cantly from one OBC to another; for instance, compared with the Reference, the amplitude of the perturbation is less than 7% for OBC5 and reaches 100% for OBC1 and OBC2. The length of the perturbed zone is more dif®cult to determine precisely, because the perturbations propagate upstream oscillating slightly. Nevertheless, for the present ¯ow the length of the perturbed zone varies approximately between H for OBC5 and 1Á5H for OBC2. Figure 7 is also a good means to evaluate the length of the entrance zone. It is approximately 2H long. Thus for x 5 2 the convective ¯ow of transversal rolls can be considered as being established. Table II presents several physical parameters computed in this zone and characterizing the thermoconvective ¯ow. These parameters are the wavelength l and time period t of the transversal rolls and the space and time average Nusselt number hNui on the bottom and top plates of the channel:

1Á4465 1Á3319 À 2Á6495 4Á3958 À 5Á0587 5Á0319 Reference 1Á439 1Á298 2Á643 À 2Á720 4Á380 À 5Á002 4Á929 OBC1 1Á413 1Á298 2Á578 0 1Á05 0 0 (100%) (76%) (100%) (100%) OBC2 1Á419 1Á296 2Á537 À 2Á348 5Á653 À 0Á469 0Á045 (14%) (29%) (91%) (99%) OBC3 1Á419 1Á298 2Á638 À 1Á474 4Á828 À 3Á741 3Á641 (46%) (10%) (25%) (26%) OBC4 1Á406 1Á298 2Á582 À 0Á558 3Á174 À 2Á559 2Á456 (80%) (28%) (49%) (50%) OBC5 1Á425 1Á296 2Á649 À 2Á616 4Á072 À 4Á909 4Á814 (4%) (7%) (2%) (2%) PBC 1Á440 1Á310 2Á662 À 2Á741 4Á358 À 4Á951 4Á859 (0Á8%) (0Á5%) (1%) (1Á4%)
where

Nut 1 LaH À 2 LaH 2 1 2 dT dz z0 xY t dT dz z1 À xY t ! dxX
The values of Umin(i), Umax(i), Wmin(i) and Wmax(i) the furthest from the Reference, for x 2 to LaH, are also given; these values are respectively denoted Umin, Umax, Wmin and Wmax and the discrepancies with respect to the Reference are given as a percentage below each value. All these results (except hNui) obtained for the ®ve OBCs and the PBCs are compared with the benchmark of Evans and Paolucci. 17 From Table II it appears that the internal structure of the ¯ow is not modi®ed by the choice of the OBCs. Indeed, taking into account the margin of error on the measurements ( AE 0Á1 on l and AE 0Á001 on t), the wavelength and time period of the transversal rolls do not vary as a function of the OBCs; the ¯uctuations of the Nusselt number are negligible too. On the other hand, the variations in the velocity can be locally very important; the discrepancies of Umin, Umax, Wmin and Wmax with respect to the Reference are small only with OBC5 and with the PBCs.

The in¯uence of the OBCs as a function of Ra and Re is studied for ¯ows corresponding to the experiments of Ouazzani et al., [START_REF] Ouazzani | Etude expe Ârimentale de la convection mixte entre deux plans horizontaux a Á tempe Âratures diffe ÂrentesÐ2[END_REF] i.e. for water ¯ows (Pr 6Á4) at an average temperature of 23 C (n 0Á93 6 10 ±6 m 2 s ±1 ) in a duct of height H 4Á15 mm. In this aim, Figure 8 gives Wmax(i) as a function of x for three combinations of Ra and Re; Wmax(i) is plotted dimensionally. Ra 1804 is close to the critical Rayleigh number Ra* characterizing the onset of thermoconvection in the Poiseuille±Be Ânard ¯ow. Note that the appearance of fully established thermoconvective rolls is rejected far from the inlet (x % 7) when Re 3 for Ra 4700 (see Figure 8(c)). When the Rayleigh number is smaller (e.g. Ra 1804) and Re 3, the thermoconvective rolls cannot appear in the channel at Pr 6Á4. [START_REF] Ouazzani | Etude expe Ârimentale de la convection mixte entre deux plans horizontaux a Á tempe Âratures diffe ÂrentesÐ2[END_REF] Whatever the combination of Ra and Re, the smallest amplitude of perturbation is always obtained with OBC5. The behaviour of the ¯ow with OBC1, OBC2 and OBC4 is approximately the same (especially at Ra 1804, Figure 8(a)). The behaviour of OBC3 is different: the perturbation oscillates at the smallest Reynolds number Re 0Á18 (Figures 8(a The fact that the length of the perturbed zone increases when Ra diminishes is in good agreement with the linear theory. The latter predicts, [START_REF] Platten | Convection in Liquids[END_REF] in the case of free convection (i.e. with rigid boundaries at the outlet), that the length necessary to obtain fully developed thermoconvective rolls varies as [(Ra 7 Ra*)aRa*] À1a2 : therefore near the boundaries the length in question tends to in®nity when Ra tends to Ra* (with Ra b Ra*). The length of the perturbed zone also increases when Re diminishes, because the perturbations caused by the OBCs can go upstream all the more easily as the forced convection is weak.

CONCLUSIONS

In the ®rst part of this paper the implementation details of the algorithm based on the augmented Lagrangian method are given; the implementation of the different open boundary conditions (OBCs) and of the periodic boundary conditions (PBCs) is also described.

The second part is dedicated to the analysis of the response of the Poiseuille±Be Ânard channel ¯ow (PBCF) solved with the ®ve different OBCs. In particular, comparison criteria are introduced to ®nely evaluate the amplitude of the perturbation and the length of the perturbed zone at the outlet boundary of the computational domain. The results show the very good behaviour of the ¯ow when Orlanskitype boundary conditions (OBC5) are used at the outlet; they also show that the ¯ow is comparatively very perturbed with the other four OBCs. On the other hand, the PBCs give excellent results; they also allow us to reduce the length of the computational domain (and subsequently to save on storage and computing time), since it is possible to work with only two space periods of the thermoconvective ¯ow.

It is shown that the improvement obtained with OBC5 is due to the important reduction of the amplitude of the perturbation without any signi®cant reduction of the length of the perturbed zone. The latter remains approximately constant whatever the OBC used. It is also shown that the length of the perturbed zone can vary from H to 5H with Ra and Re; the length increases when Ra and Re diminish.

Figure 1 .

 1 Figure 1. Augmented Lagrangian algorithm

Figure 2 .

 2 Figure 2. Grid system and positioning of variables ((a) T and P; (b) U; (c) W) near outlet of computational domain

Figure 2 (

 2 Figure 2(a) or 2(c), applying the Neumann boundary condition dfadx 0 gives f E f W . By introducing this expression in (9), the coef®cients are modi®ed in the following way:

  and a S W remain unchangedX OBC4 is a Fourier-type boundary condition: this is a mixed condition between Neumann and Dirichlet conditions. In the case of open ¯ows (e.g. ¯ows around an obstacle) it permits the connection of the recirculating local ¯ow with the fully established parallel ¯ow at in®nity, (U I , W I , T I ), by adjusting the three positive constants (b U , b W , b T ). In the case of channel ¯ows, (U I , W I , T I ) generally represent the upstream non-perturbed ¯ow. In the present paper they represent the thermally strati®ed Poiseuille ¯ow: U I À 6(z 2 7 z), W I 0 and T I 1 7 z. Neumann boundary conditions are obtained when (b U , b W , b T ) (0, 0, 0) and Dirichlet boundary conditions when (b U , b W , b T ) (I, I, I). By choosing (b U , b W , b T ) between these two limits, the outlet behaviour of the ¯ow can be considerably improved. Since the three constants depend on the nature of the ¯ow and on the chosen non-dimensionalization, the best compromise for (b U , b W , b T ) cannot be obtained by a straightforward theoretical principle; it can only be obtained by numerical experiments. For the present work, after having tested and compared several sets of (b U , b W , b T ), the optimal value of the triplet was found to be (b U , b W , b T ) % (5, 5, 0).

Figure 5 .

 5 Figure 5. Field plots for different OBC treatments at time t* (Ra 10,000, Re 10 and Pr 2 3 )

Figure 7 .

 7 Figure 7. Envelopes of U and W extrema along axis of channel; visualization of perturbed zone for the different tested OBCs (Ra 10000, Re 10 and Pr 2 3 )

Figure 8 .

 8 Figure 8. OBC effect upon envelope of vertical velocity component as a function of Rayleigh and Reynolds numbers

  ) and 8(b)), particularly when Ra 1804. As it has already been shown in the case of the ¯ow of Evans and Paolucci, the length of the perturbed zone varies very little from one OBC to another. On the other hand, it varies with Ra and Re: it is around 5H in Figure8(a) for Ra 1804 and Re 0Á18, around 3H in Figure 8(b) for Ra 4700 and Re 0Á18 (a bit more for OBC4) and around H in Figure 8(b) for Ra 4700 and Re 3.

Table I .

 I Inlet, bottom and top boundary conditions as well as initial conditions for ®rst computational con®guration (with OBC)

	OBC2
	Inlet BC (x 0)	Bottom BC (z 0) dU dx À dW dz Y	dW dx	Top BC (z 1) 0Y dT dx	0X	Initial condition (t 0)

Table II .

 II Comparison of ¯ow parameters obtained with different OBCs and with Reference values for Ra 10,000, Re 10 and Pr 2 3 . The percentages give the discrepancy with respect to the Reference

	l	t	hNui	Umin	Umax	Wmin	Wmax
	Evans and Paolucci 17

ACKNOWLEDGEMENTS

The numerical simulations for the present work were performed on the IBM-SP2 computer of the CNUSC (South University National Centre of Computation). The authors would like to thank the staff of the CNUSC for their technical support and training.


F Michaud

J Y Collot

G Ratzov

Jean-Noël N Proust

A Dano

J F Lebrun

M J Hernández

G Loayza

A Khaoulani

Y Stoll

H Pouderoux

L De Min

A honeycomb

INTRODUCTION

Only a few cases of kilometer-scale honeycomb pattern (HP) on the seafloor have been reported along continental margins and generally into siliceous or mudstonedominated sediments. The main examples are [START_REF] Berndt | Kilometre-scale polygonal seabed depressions in the Hatton Basin, NE Atlantic Ocean: Constraints on the origin of polygonal faulting[END_REF] kilometer-scale seabed high-centered polygons [START_REF] Berndt | Kilometre-scale polygonal seabed depressions in the Hatton Basin, NE Atlantic Ocean: Constraints on the origin of polygonal faulting[END_REF] outlined by 20-m-deep, elongated pockmarks, formed by sediment contraction related to fluid circulation along polygonal faults; and (2) polygonal crests separating hectometer-to kilometer-scale, 60-m-deep depressions formed by longterm unstable bottom currents [START_REF] Sun | Three-dimensional seismic characterization of a complex sediment drift in the South China Sea: Evidence for unsteady flow regime[END_REF]. Here, we used new multibeam data collected during the 2012 ATACAMES cruise (RV L'Atalante) on the Carnegie Ridge offshore Ecuador (http://campagnes.flotteoceanographique.fr/campagnes/12010010/; Fig. 1) to report, for the first time, a HP morphology in carbonate sediment.

The HP consists of a network of kilometer-scale, hectometer-deep, polygonal depressions (Figs. 2 and3). Based on a low-resolution bathymetric data set, the depressions were initially reported as densely packed subcircular features [START_REF] Michaud | Fields of multi-kilometer scale sub-circular depressions in the Carnegie ridge sedimentary blanket: effect of underwater carbonate dissolution?[END_REF]. The ATACAMES high-resolution multibeam data and seismic profiles (72 tracks, 50-450 Hz) allow the nature of the polygonal depressions to be investigated. We present a coherent model for their formation in relation to sedimentation, fluid circulation, and large-scale seafloor erosion. Additionally, we discuss the timing, the origin, and the widespread distribution of the HP in relation to the closing of the Isthmus of Panama and subduction of the Carnegie Ridge.

CARNEGIE RIDGE FRAMEWORK

The interplay of the Galapagos hotspot with the Cocos-Nazca spreading center led to the formation of the Cocos, Malpelo, and Carnegie Ridges (Fig. 1A; [START_REF] Lonsdale | Erosional troughs in deep-sea carbonates and their relationship to basement structure[END_REF]. Deep Sea Drilling Project (DSDP) Site 157 [START_REF] Van Andel | Deep Sea Drilling Project Site 157[END_REF], and Ocean Drilling Program (ODP) Sites 1238 and 1239 (Fig. 1B; Mix et al., 2003) reveal a 400-500-m-thick Miocene to Pleistocene carbonate sediment cover over the oceanic basement. Two lithologic units were defined at Sites 157 and 1238, an upper unit, consisting of ~300 m of oozes and chalk, and a lower unit, composed of ~100 m of lithified interbedded chalk and chert horizons (Fig. DR1 in the GSA Data Repository1).

Subcircular depressions fields are present on the Carnegie, Cocos, and Malpelo Ridges (Fig. 1A), and they have been globally attributed to regional submarine carbonate dissolution processes enhanced by bottom currents [START_REF] Lonsdale | Submarine Geology of the Malpelo Ridge[END_REF].

DEPRESSION GEOMETRIES

The seafloor depressions of the HP area are typically 0.8-2.0 km wide, separated by 100-250 m high, dominantly linear ridges (supplemental information Table 1) (Fig. 2, 3, and Fig. S2). The ridges give the depressions their polygonal and, by places, hexagonal shape, and enclose sub-circular depression bottoms (Fig. 3). Most polygonal depressions lay in 2100-2800 m of water depth together with isolated, sub-circular to ovoid, larger depressions (3-4 km-wide and 400 m-deep) (Fig. 2); the largest being located deeper than 2500 m.

ACOUSTIC BASEMENT AND FEATURES OF SEDIMENTARY COVER UNITS

Acoustic Basement

The acoustic basement is marked by a set of strong and chaotic reflections (Fig. 4) overlain by a thinly laminated sedimentary cover 0.4-0.8 s two-way traveltime (TWTT) in thickness. The basement shows a remarkable step that is 0.8 s TWTT high and trends east-west in map view (Figs. 2 and4). South of the step, the acoustic basement is smooth and lies at a 4.3-4.5 s TWTT in depth; the sediment is 0.8 s TWTT thick, locally carved by isolated depressions and moats (Fig. 2). North of the step, the acoustic basement is shallower and locally crops out at the seafloor (stars on Fig. 4B). When probed by seismic data (Fig. 2), the depressions were observed to have preferentially (73%) developed where the acoustic basement is rough (Table DR1). However, where the acoustic basement is smooth, some lows filled with sediment represent channel-fill structures (Fig. 4A, profile ATAC129, zoom1). In the sedimentary cover, we differentiated a lower unit, u1, from an upper unit, u2 (Fig. 4), separated by a major erosional unconformity.

Lower Unit u1

The lower unit u1 (Fig. 4A) returns subparallel continuous reflectors. Unit u1 extends over the entire region, including the acoustic basement step, but it is absent in some depressions. The base of u1 locally shows a thin transparent or poorly reflective layer (Fig. 4A, profile ATAC129, zoom1). South of the basement step, u1 is 0.6 s TWTT thick. There, the poorly reflective lower half of u1 is bounded at the top by a domeshaped reflector crossing the original subparallel reflectors (Fig. 4A, Common Depth Point [CDP] 5000-7000; Fig. DR3, profile ATAC131 CDP 1000-2500). This domeshaped reflector suggests that a diagenetic front rose in the sediment cover crossing the original strata. To the north of the step, the thickness of u1 varies from 0 to 0.6 s TWTT further north along profile ATAC126 (Fig. 4B). An undulating geometry showing locally internal cutand-fill structures characterizes the unit above the rough basement area (Fig. 4B, zoom2 and zoom3;Fig. 4A, profile ATAC130, CDP 2000-4000), suggesting a sediment drift pattern (Faug.res et al., 1999).

The major unconformity truncates u1 in the buried depressions (Fig. 4A,ATAC129, CDP 1000-3000, zoom1; Fig. 4B,zoom2), and at the walls of most of the seafloor depressions (Fig. 4A, CDP 3600-4000, zoom2; Fig. 4B, CDP 15800-16500). In cross section, this erosional surface outlines the triangular shape of the ridges between the polygonal depressions carved into u1 (Fig. 4B, zoom2 and zoom3).

Normal faults with 5-to 20-m-scale vertical offsets (red arrows in Figs. 4A and4B, zoom1 and zoom4) deform u1. Most of them show characteristics of polygonal faults [START_REF] Gay | Polygonal faults furrows system related to early stages of compaction Upper Miocene to recent sediments of the Lower Congo Basin[END_REF] as they disappear upward and downward, do not end at a specific horizon, show a low frequency in lateral distribution, and may link up in polygons in plan view, as demonstrated in other parts of the world where three-dimensional (3-D) seismic control is available [START_REF] Cartwright | The genesis of polygonal fault systems: A review[END_REF].

Although time migration with a constant 1500 m/s velocity produced clear seismic imaging, imperfect migration velocities led to poor seismic imaging, preventing local geological interpretation. Indeed, only vertical acoustic anomalies associated with wellcollapsed diffractions in u1 (Fig. 4A, zoom1 and zoom4) are interpreted as fluid-escape pipes [START_REF] Cartwright | Seismic characteristics of fluid escape pipes in sedimentary basins: Implications for pipe genesis[END_REF].

Upper Unit u2

The upper unit u2 is thinly stratified and discontinuous, and its thickness varies from 0 in some depressions to 0.4 s TWTT south of the basement step. Although locally concordant with u1 (Fig. 4A, CDP 5000-7000; Fig. 4B, CDP 2000-4000), unit u2 unconformably overlies (Fig. 4, zoom1, zoom2, and zoom3) the erosional topographies carved in u1. Unit u2 shows internal discontinuities and undulating depositional patterns pointing to sediment drift (Fig. 4A,zoom1) and exhibits complex migration and aggradation geometries, especially in the rough acoustic basement areas (Fig. 4B,zoom2, 3; Fig. DR3, ATAC128 CDP 7000-8500, ATAC131 CDP 5000-7500). These geometries resemble contourite drifts produced by the interaction between oceanic bottom currents and the seafloor topography.

DISCUSSION

Origin of The Honeycomb Seafloor Morphology

smaller counterparts developed in underlying sedimentary units as a result of a sediment drift depositional pattern [START_REF] Sun | Three-dimensional seismic characterization of a complex sediment drift in the South China Sea: Evidence for unsteady flow regime[END_REF]. This pattern is characterized by wave-like geometries corresponding to cut-and-fill erosional features with honeycomb planforms. It was attributed to the action of unsteady bottom currents initiated in the late Miocene in the South China Sea. In our study area, without 3-D seismic data, the HP seafloor cannot be demonstrated to mimic buried honeycomb structures. Nevertheless, this hypothesis is not ruled out, as unit u1 exhibits wave-like geometry where the acoustic basement is rough (zoom2, and basement step along profile ATAC130 in Fig. 4), and depression 3 (Fig. 4B) shows that unit u2 sediment drift filled up a polygonal depression carved in u1, reproducing the underlying honeycomb pattern. Many polygonal depressions (Fig. 3A) carved into u1 (Fig. 4B, zoom2 and zoom3) are, however, deeper (~200-300 m) than those (~60 m) described by [START_REF] Sun | Three-dimensional seismic characterization of a complex sediment drift in the South China Sea: Evidence for unsteady flow regime[END_REF], suggesting that other important processes are involved.

Polygonal faults generally form in mud-rich sediment [START_REF] Cartwright | The genesis of polygonal fault systems: A review[END_REF] by dewatering and, among other processes [START_REF] Goulty | Geomechanics of polygonal fault systems: a review[END_REF], by differential compaction of the host sediment [START_REF] Berndt | Kilometre-scale polygonal seabed depressions in the Hatton Basin, NE Atlantic Ocean: Constraints on the origin of polygonal faulting[END_REF] related, to either sediment thickness variations (zoom4, Fig. 4; [START_REF] Lonsdale | Erosional troughs in deep-sea carbonates and their relationship to basement structure[END_REF] or silica diagenesis [START_REF] Davies | Initiation and propagation of polygonal fault arrays by thermally triggered volume reduction reactions in siliceous sediment[END_REF]. On the Carnegie Ridge, in the areas where we interpret a diagenetic front (Fig. 4A, CDP 4800-5500; Fig. DR3,ATAC131, CDP 2000-2500), u1 reflectors and the seafloor draw an elongated depression controlled by normal faults, likely relating to diagenesis conversion. Indeed, according to DSDP Leg 16 and ODP Leg 202, the sedimentary cover of the Carnegie and Cocos Ridges is dominantly composed of carbonate, but the sediment close to the basement reveals chert and micrite, indicating significant opal and carbonate diagenesis (Mix et al., 2003;[START_REF] Moore | Chert in the Pacific: Biogenic silica and hydrothermal circulation: Palaeogeography[END_REF]. Consequently, a polygonal fault network is expected in the studied HP area. According to [START_REF] Gay | Polygonal faults furrows system related to early stages of compaction Upper Miocene to recent sediments of the Lower Congo Basin[END_REF], furrows are initially produced at the seafloor along polygonal faults, with pockmarks developing above highdrainage chimneys at the triple junction of three adjacent hexagonal cells (Figs. 5A and5C). In a carbonate environment, [START_REF] Moore | Hydrothermal pits in the biogenic sediments of the equatorial Pacific Ocean[END_REF] showed that hydrothermal discharge contributes to sediment fracture and damage, creating pits over oceanic basement highs (Figs. 5B and5D). Based on a mass balance analysis, [START_REF] Bekins | Dissolution of biogenic ooze over basement edifices in the equatorial Pacific with implications for hydrothermal ventilation of the oceanic crust[END_REF] indicated that the formation of pits in the equatorial Pacific could have been enhanced by dissolution of carbonate sediment by fluids undersaturated with calcite exiting basement vents. Fluids are therefore considered as a key factor in the initiation of seafloor depressions on the Carnegie Ridge.

Subsequently, the pockmarks and pits were enlarged and deepened in response to the action of strong bottom currents [START_REF] Sun | The morphologies and genesis of mega-pockmarks near the Xisha Uplift, South China Sea[END_REF], shaping the regional unconformity between u1 and u2 [START_REF] Lonsdale | Erosional troughs in deep-sea carbonates and their relationship to basement structure[END_REF]. The currents preferentially removed sediments weakened by fracturing and fluid circulation, scouring depressions down to the indurated base of u1, or to the top of the oceanic basement (Fig. 5E). In the pit case, the currents action allowed the pit to grow radially, leaving at the end a polygonal depression pattern. In the pockmark case, the impact of current scouring is uneven, because it is controlled by the current direction with respect to the orientation of the hexagonal cell pattern. Indeed, triple junctions of adjacent cells offer both convergent and divergent (Fig. 5A) geometries for current action. The current is funneled and accelerated by the convergent geometry, so that erosion is locally enhanced, initiating subcircular depressions at every other triple junction. This process allows the pockmarks to grow in diameter, leaving at the end an HP that is offset laterally by a half-hexagonal cell with respect to the initial polygonal fault network. Many depressions show elongated or free shapes locally crossed by minor ridges (Fig. 3B) due to irregular basement topography, and partial or full erosion of some interdepression ridges. At the end, u2 sediment drift is plastered against the wall of some polygonal depressions carved in u1 (Fig. 5F).

Possible Regional Cause for the Scouring of the HP

Seismic records collected near DSDP Site 157 (Fig. 1B; [START_REF] Heath | Tectonics and sedimentation in the Panama basin: Geological results of Leg 16[END_REF][START_REF] Van Andel | Deep Sea Drilling Project Site 157[END_REF] show several unconformities in the upper 110-125 m of Pleistocene sediment. DSDP Site 157 drilled u2 immediately west of our study zone (Fig. 1B) and recovered 116 m of Pleistocene sediment above the Pliocene sequence (Fig. DR1). On profile ATAC126, u2 is ~80 m thick (Fig. 4B, CDP 2000-4500) and could then include the entire Pleistocene sequence drilled at DSDP Site 157. Consequently, the highly erosive episode that we associate with new current conditions is likely to have started at the transition between Pliocene and Pleistocene times. At a regional scale, DSDP Site 158 (Fig. 1A), drilled on the Cocos Ridge, revealed a shallow hiatus covering most of the late Pliocene [START_REF] Heath | Tectonics and sedimentation in the Panama basin: Geological results of Leg 16[END_REF]. Further northeast along the Cocos Ridge, ODP Site 1242 (Fig. 1A) shows a major unconformity (Mix et al., 2003) and a hiatus that encompasses the interval from 13 to 2.5 Ma. On the Malpelo Ridge, kilometer-scale erosional depressions were initiated during the mid-Pliocene [START_REF] Lonsdale | Submarine Geology of the Malpelo Ridge[END_REF]. If the origin and the age of the depression fields on the Carnegie, Cocos, and Malpelo Ridges are identical to those proposed in this study, the erosive episode appears to be of regional extent.

Major changes in the configuration of ocean gateway and related deep currents are triggered by geodynamic events [START_REF] Moore | Paleogene tropical Pacific: Clues to circulation, productivity, and plate motion[END_REF]. The formation of the Isthmus of Panama stands as one of the greatest geological events driving profound transformations in ocean circulation. The collision of the Carnegie Ridge with the South American margin clogged the subduction trench, contributing, together with the gradually emerging land of Panama, to the enclosure of the Panama Basin, and to the change in the regional bottom current configuration [START_REF] Lonsdale | Submarine Geology of the Malpelo Ridge[END_REF].

In a review, O'Dea et al. ( 2016) dated the emergence of the Isthmus of Panama at 2.8 Ma, close to the inferred age for the beginning of the current-driven, erosive regional episode observed on the Carnegie Ridge that carved and partially filled the honeycomb polygonal depression network.

In conclusion, polygonal faults, which are globally scarce in the studied carbonate ooze environment, were not clearly imaged in the HP area, possibly because they disappeared during the scouring process, as suggested in our model. In contrast, seafloor polygonal depressions statistically correlate with the rough basement topography. Hydrothermal fluid circulation and dissolution affecting the carbonate sediment overlying the rough oceanic basement provide support for the pit origin of the depressions. Strong bottom currents acting at the Pliocene-Pleistocene boundary, possibly as a result of the closure of the Isthmus of Panama, were instrumental in scouring the polygonal depressions.

FIGURE CAPTIONS

Figure 1 .

 1 Figure 1. A: Geodynamic setting (CNS-Cocos-Nazca spreading center; HS-hotspot).

Figure 2 .

 2 Figure2. Map of studied area (same color scale as that in Fig.1B). In pale and dark

Figure 3 .

 3 Figure 3. A-B: Close-ups (zooms) of depressions field (locations in Fig. 2). Red line

Figure 4 .

 4 Figure 4. Seismic profiles (see location in Fig. 2). A: Profiles ATAC130-129 showing

Figure 5 .

 5 Figure 5. Diagram showing two scenarios for origin of honeycomb pattern (HP). A:

  FIGURE 1

Figure DR1 :

 DR1 Figure DR1: Correlation between our seismic units (u1, u2, acoustic basement) with a seismic line, drilled lithology and age from the DSDP 157 (from van Andel et al., 1973).

Figure

  Figure DR2: A) Bathymetric map southeast of the studied area along the southern flank of the Carnegie ridge where high resolution multibeam were recorded during the Atacames cruise B) 3D zoom of the high resolution multibeam, showing honeycomb pattern around the seamounts

  TableDR1: Depressions characteristics along the seismic profiles. The numbers correspond to the depressions along the seismic profiles and located on Figure2,3, 4 and S1). Semi-quantitative estimation based on all seismic lines crossing the depressions that we have acquired in the area (figure 2 for line location). Some statistic is done about 1) the link between the basement roughness and the depression location (73% of depressions are related to a basement irregularity);and 2) about the sedimentary filling (38% of the depression are filled by u2 sedimentary unit.

Figure DR1 :

 DR1 Figure DR1: Correlation between our seismic units (u1, u2, acoustic basement) with a seismic line, drilled lithology and age from the DSDP 157 (from Van Andel et al., 1973).

Figure

  Figure DR2: A) Bathymetric map southeast of the studied area along the southern flank of the Carnegie ridge where high resolution multibeam were recorded during the Atacames cruise (same color scale as that in Fig. 1B and 2). In pale color the multibeam before the Atacames cruise (grid size 150m). In dark shade, multibeam from the Atacames cruise (grid size 50m). Iso contours=25m. B) 3D zoom high resolution multibeam showing honeycomb pattern around the seamounts (location in Fig S2 A). Vertical exaggeration x3. Northward vertical view.

Figure DR3 :

 DR3 Figure DR3: Interpretation of the available seismic profiles (ATAC 131, ATAC 128, ATAC 127 and ATAC 124S, location on Figure 2) in the study area in addition to the seismic profiles of Figure 4. With this additional figure, all available seismic lines are shown. Short dotted black line=acoustic basement; large dotted black line=unconformity between lower and upper seismic units; blue arrows=inferred upper limit of diagenesis and lithification; red rectangle=location of the zooms; u1=lower unit (u1a=poorly reflective; u1b = sub parallel continuous reflectors); u2=upper unit. Se = Side echo related to possible irregularities of the basement. Numbers correspond to the depressions considered in the table 1. Dotted blue line = seafloor depth EM122 vertical beam allows to discriminate side echoes figured above the seafloor.

ACKNOWLEDGMENTS

We thank the crew of R/V L'Atalante and GENAVIR. This work was supported by Institut National des Sciences de l'Univers (INSU) du Centre National de la Recherche Scientifique (CNRS). Thanks to Institut de la Recherche et du Développement (IRD), to Instituto Oceanográfico de la Armada (INOCAR) and Ministry of Higher Education, Science, Technology and Innovation of Ecuador (SENESCYT).