Dr Mohammed Adil Sbai
email: a.sbai@brgm.fr

F De Smedt

A Larabi

A Generalized Approach for Modeling 3D Transient Free and Moving Boundaries in Coastal Aquifers

Keywords: Free surface, seawater intrusion, numerical modeling, finite element method

Coastal aquifers involve varying conditions in time and space, owing to the occurrence of free and moving boundaries, such as the water table, seepage face and saltwater intrusion interface. A fast Updating Procedure (FUP) is developed for solving these interfaces in steady and transient conditions with the finite element method. Several test examples, which involve confined and phreatic groundwater flow in coastal aquifers, are studied to demonstrate the FUP capability of predicting accurate results. Comparisons of the FUP numerical results are made against available analytical solutions, laboratory experiments measurements and other computer codes.

To show also the capability of this technique to solve real field situations, it was applied to the coastal aquifer of Martil in Morocco to study and understand the aquifer response to changes in recharge and total rate of pumping water, and their effects on seawater intrusion.

INTRODUCTION

The problem of saltwater intrusion can be treated with two methods, one considers that saltwater and freshwater are miscible and that a transition zone exists [START_REF] Voss | AQUIFEM-SALT: A Finite-element model for aquifers containing a seawater interface[END_REF], in which situation miscible density dependent groundwater flow models are traditionally used. The other method is based on an abrupt approximation [START_REF] Bear | Dynamics of fluids in porous media[END_REF][START_REF] Huyakorn | Multiphase approach to the numerical solution of a sharp interface saltwater intrusion problem[END_REF][START_REF] Larabi | Numerical solution of a 3D groundwater flow involving free boundaries by a fixed finite element method[END_REF][START_REF] Sbai | Modeling three dimensional groundwater flow and transport by hexahedral finite elements[END_REF] and for which (i) the transition zone is thin relative to the thickness of the freshwater lens, (ii) freshwater and saltwater are immiscible, and (iii) the freshwater and saltwater are separated by a sharp interface. Additionally, the approach considering flow in both fresh and salt water regions, leads to a coupled system of equations to be solved. This is cumbersome for a fully three dimensional model owing to numerical constraints and higher computational costs involved. Instead, a new numerical procedure is developed, assuming a sharp interface between freshwater and saltwater, for which the location, shape and extent must be determined. A first steady state version of this model was presented by [START_REF] Larabi | Numerical solution of a 3D groundwater flow involving free boundaries by a fixed finite element method[END_REF], who elaborate a part of the numerical technique which will be described, and also validate the developed model by running several benchmarks. In this work an enhanced version of their steady state model is developed, several features were implemented as a robust solver based on an M-matrix preconditioned incomplete factorization. Also, special boundary conditions such as the seepage face, and sea outflow face are efficiently handled. But, the most important achievement is the development of a transient model for simulating the moving sharp interface and the moving water table in case of free surface flow conditions. Numerical algorithms used as a basis for this model are similar to those for the steady state model for updating unsaturated hydraulic conductivity, with a small exception. The main difficulties due to the strong non-linearly of the equations systems, and which are related to estimation of the capacitance matrix coefficients or the so-called cord slope approximation, are tackled for the case of moving interfaces by adopting an idealized soil characteristic curve whose parameters are internally estimated, and which are depending on the local saturation in the interfaces neighborhood, and the finite element mesh configuration. This strategy proves to be very efficient in combination with a modified Picard iteration featuring a variable relaxation technique. A full description of the theory and the methodology used for developing GEO-SWIM code is found in [START_REF] Sbai | Modeling three dimensional groundwater flow and transport by hexahedral finite elements[END_REF]. In this article the key developments are highlighted, followed by an exposition of theoretical problems for which the technique has been validated, and a realistic field application of the model to the Martil aquifer in the northern part of Morocco.

THEORETICAL DEVELOPMENTS

The nonlinear partial differential equation governing the variable saturated flow in 3D heterogeneous porous medium [START_REF] Huyakorn | Computational methods in subsurface flow[END_REF]] is given by

R h t h h S + ∇ ∇ = ) ( ) ( K K K K ∂ ∂ (1)
Where S(h) = S s + dθ/dh is the general storage term [L -1 ] the first term being the specific storage coefficient and the second being the water capacity, h is the groundwater potential [L], K is the unsaturated hydraulic conductivity tensor [LT -1 ] depending upon position and pressure or degree of saturation, R is the source/sink term [T -1 ], and ∇ is the del operator [L -1 ], (∂/∂x, ∂/∂y, ∂/∂z) T , where x = (x,y,z) T is the spatial vector [L], and t represents time [T]. In order to have a unique mathematical problem, boundary conditions and the initial potential distribution h(x,t) are necessary. The location of the water table and the fresh-saltwater interface are unknown a-priori and hence the nonlinearity of the problem rises in consequence.

The finite element method in space using hexahedral elements, coupled within a fully implicit finite difference method in time, are used to solve (1). This results in a nonlinear system of algebraic finite element equations. in matrix form at the (k+1) th time

step t k+1 , k = 0, 1, … } { } { ] [ ] [ 1 B h S G = ÷ ø ö ç è ae ∆ + + k t t (2) Where 1 } { + k t h
is the unknown vector containing the nodal potentials at time 1

+ k t ,
] [G is the global conductance matrix depending on the geometrical and conductive properties of the flow domain, ]

[S is the diagonalized storage matrix, and the right hand side vector } {B incorporates boundary conditions together with sources, sinks, and storage terms at time step t k .

In the present approach a fixed finite element mesh is used to discretize the whole domain, including saturated, unsaturated and saltwater regions. The groundwater flow is assumed to occur only in the saturated freshwater zone, which means that the flow in the unsaturated and saltwater regions are negligible compared to the freshwater flow. The saltwater interface and the water table are effectively considered as impervious boundaries, but not as boundary conditions of the problem.

The fast updating procedure technique

This technique involves implicit updating of the saturated flow field independently of the unsaturated and saltwater zones, while they are still included in the meshed domain. This is performed in two steps, the first concerns the approximation of the relative hydraulic conductivity between nodal points in the neighborhood of the water table and the fresh-saltwater interface, this method is described in [Sbai et al., 1998;Sbai and De Smedt, 1998a] and in more detail in [START_REF] Sbai | Modeling three dimensional groundwater flow and transport by hexahedral finite elements[END_REF]. The second part of the algorithm, approximates the nonlinear storage or the time dependent term in the right hand side of (2). The nodal storage depends on the water table and the saltwater interface positions, and is evaluated as

dp d V i θ (3)
Where V i is the control volume patch centered on node i, and p is the water pressure [L]. The method used to evaluate the derivative term in (3) affects significantly the convergence behavior of the iterative schemes, due to steep gradients and discontinuties or points of inflection in the soil curves. The functionals showed in Fig. 1 are used to achieve the numerical differentiation method adopted in this model. At the beginning of each time step the differential in function (a) is directly used to approximate the chord slope, in the next iterations nodal water content values are relaxed following the function (b) and the differential expression is calculated thereafter. The limiting depth parameter d used for the water table zone is evaluated as

i i i S V d = (4)
Where i S is the projected surface area attributed to node i over the horizontal plane. The storage variations due to changes of the saltwater zone displacement are evaluated to be equal to the saturated water content θ s . The newly distinguished pressure distribution around the fresh-saltwater interface is smaller in size in comparison with the equivalent pressure distribution existing around the water table, because small variations in the water table position involves greater displacement of the fresh-saltwater interface. And, hence much larger variations in the storage term. The coefficient of proportionality is taken equal to the density difference ratio δ. The merit of the FUP procedure results from the fact that only relative hydraulic conductivity and nodal storage values have to be adapted, while the finite element mesh and saturated conductance coefficients remain fixed, such that the method is computationally fast.

P, p

MODEL APPLICATIONS AND RESULTS

Free surface flow and seawater intrusion in coastal aquifers

Several test examples were studied in steady conditions analytical [START_REF] Larabi | Numerical solution of a 3D groundwater flow involving free boundaries by a fixed finite element method[END_REF]] and time dependent conditions [Sbai et al., 1998] to demonstrate the numerical technique capability of predicting accurate results efficiently. A comprehensive validation of the model predictions is achieved by comparing to laboratory experimental measurements of 3D moving free surface flow in earth dam model and a moving saltwater interface in a 3D laboratory sand box model of irregular shape [START_REF] Sbai | Modeling three dimensional groundwater flow and transport by hexahedral finite elements[END_REF]. These problems are rather complex because they allow for real 3D flow and for heterogeneous dam material. The predicted 3D results compare well with the measurements. More other applications dealing with free surface flow and saltwater intrusion, have been also studied for comparison with other numerical schemes, but not included in this paper.

The Martil Coastal aquifer (Morocco)

The Martil aquifer system is located in the northwestern part of Morocco on the Mediterranean coast (Fig. 2). It is the important local aquifer system in the region used for public supply, irrigation and industry in Tetouan city. In recent years, the aquifer has become vulnerable to potential pollution, especially due to seawater intrusion from the Mediterranean sea. A three-dimensional finite element model of this aquifer is developed using the FUP procedure. First, a steady-state groundwater flow is simulated. Calibration of this model will establish natural conditions and hydraulic conductivity ranges for the different aquifers. Afterwards a transient simulation is performed to predict future lateral extension of the saltwater encroachment due to pumping of groundwater.

The aquifer system consists of two aquifer units, separated by a leaky aquitard. The upper aquifer is formed of Quaternary alluvial deposits of the Martil river; the lower aquifer unit is composed of sandstone-limestone Pliocene formations, while the aquitard is mainly marl and clay. Variations in thickness are shown to be significant from North to South, and also along the West-East direction. To better characterize the thickness of the three hydrogeological units a reinterpretation of data obtained from a total of 59 wells and boreholes is performed with support of GIS tools.

The conceptual model used in this study is a multi-layer aquifer system of variable thickness. This aquifer system is assumed isotropic and homogeneous with each subregion. A structured surface mesh of 93 columns and 121 rows is used to approximate the aquifer domain, and some sub-zones are set to be inactive to fit the remaining part at the domain boundaries, especially at the eastern boundary near Tetouan city. In total, the aquifer system is divided into 9 nodal layers, with 88320 hexahedral elements and 101277 nodes. The finite element mesh used for the numerical simulation of the Martil aquifer system which is adjusted to fit the structure and the extension of the hydrogeological layers and the wells locations. Domain boundaries are set impervious, except at the eastern part, where the aquifer is in direct contact with the sea (outflow face boundary condition). The extent of the fresh water outflow to the sea is automatically determined by the model as part of the results. Nodes along the rivers paths are taken as fixed head with the assumption that water levels are equal to the elevation. Another condition taken into account is the effective recharge, which is assumed uniformly distributed over the whole surface basin. Saturated hydraulic conductivity values are deduced from calculated transmissivities of pumping tests analysis, conducted by the Regional Hydraulic Department of Tetouan.

A steady state simulation is performed under natural conditions (i.e. no pumping wells). The objective of this calibration process is to reproduce a natural groundwater flow pattern of the aquifer system and at the same time provide a range of confidence limits for the model parameters. The obtained results for each test run were compared to piezometric levels measured in 1966, for which it is assumed that the aquifer was not yet heavily pumped at that time. A trial and error calibration procedure is used to estimate the hydraulic conductivity of different layers to best reproduce the reference piezometric levels (1966). The tests show that the conceptual model is more sensitive to changes in the effective recharge value. Computed groundwater potential heads versus observed values are compared in Fig. 3, the fit is satisfactory except at the center of the plain and near the coastline where differences exist.

In general, the model is able to reproduce the same flow pattern, indicating that the main groundwater flow is directed W-E, with some convergence tendency to the rivers. the fresh-salt water interface is computed iteratively in parallel with the groundwater potential heads. W-E cross-sections displayed in Fig. 4 show the shape of the steady interface and water heads distribution in different parts of the aquifer. By comparing cross-sections C-C' and E-E' which are respectively parallel to the Alila and Martil rivers, it follows that saltwater intrusion is more sensitive along the Martil river as the hydraulic gradient is the smallest due to the flat topography. The parameters obtained after calibration are: saturated hydraulic conductivity K s = 3.5m/d, and porosity n = 0.4 in the alluvial formation; saturated conductivity K s = 0.05m/d, and porosity n = 0.25 in the clay/marl aquitard; saturated conductivity K s = 4.5m/d, and porosity n = 0.25 in the Pliocene limestone-sandstone formation; a freshsaltwater density difference ratio δ = 0.03; and an estimated effective recharge value of 50mm/y. Future prediction scenario for seawater intrusion. A second set of simulations is performed, to predict the present and future groundwater flow and seawater intrusion in the aquifer. A long-term transient simulation is made of 40 years starting in 1966, and using the steady state conditions from the first set simulations as initial head values. The pumping of groundwater is assumed to decrease the natural recharge by half during this period. Fig. 5 shows the positions of the computed moving interface plotted along cross-section F-F'.

CONCLUSION

A fast updating numerical procedure is developed for solving 3D groundwater flow involving saltwater intrusion and phreatic flow in steady and transient conditions.

The model uses a fixed FE mesh, and iteratively adjusts the water table and the saltwater interface position by neglecting the flow in the unsaturated and the saltwater zone. This is achieved without fully reconstructing the conductance and capacitance coefficients matrix of the FE element equations. During the solution procedure, the saturated conductance and capacitance coefficients are constant. For the conductance matrix, only the relative conductivities have to be adjusted during each iterates, depending on the status of water at each node. The capacitance matrix coefficients are also adjusted in the unsaturated and saltwater zones depending upon the change in the storage. The proposed numerical approach has been verified using a series of test problems under steady and transient flow situations. The comparison study has shown that the technique is accurate when compared with available analytical solutions.

Validation of the model is also made with respect to measurements from a 3D laboratory sand box model and other numerical schemes.

This model is applied to the aquifer system of Martil situated in Morocco using the FUP numerical procedure, in order to study seawater intrusion effects in terms of the shape and lateral extension. This study enables to understand the aquifer response to changes in recharge and total rate of pumped water, and their effects on seawater intrusion. Different scenarios are investigated for the period of 1966 to 2006, to predict future situations and the salinization risk from seawater intrusion. The obtained results show that the interface will move fast and travel over considerable distances in forthcoming years, and will produce an irreversible degradation of the groundwater quality. An alarming optimal management scheme to adopt in the near future is necessary for its safeguard. Sine this is the first time that a simulation model for groundwater flow and seawater intrusion in the Martil aquifer system is performed, this model should be improved as further data are obtained from the field.

Figure 1 .

 1 Figure 1. Saturation curves used in the model (a) at time periods start up and (b) during the iterative process.

Figure 2 .

 2 Figure 2. Study area and cross-sections locations.

Figure 3 .

 3 Figure 3. Comparison between (a) computed and (b) observed steady state groundwater potentials (presented as meters above sea level).

Figure 4 .

 4 Figure 4. Cross sectional views of the computed groundwater heads and fresh-salt water interface position.

Figure 5 .

 5 Figure 5. Cross-sectional view of the seawater intrusion; simulated moving interface positions from 1974 to 2006 with 8 years interval are plotted.


Abdoulaye Fall

Anael Lemaitre

François Bertrand

Daniel Bonn

Guillaume Ovarlez

Anaël Lemaître

Shear thickening and migration in granular suspensions

Keywords: numbers: 83, 80, Hj, 83, 60, Fg

The realm of complex fluids encompasses biological liquids such as blood, many liquid foodstuffs, building materials, glasses and plastics, crude oil, etc. Despite their importance, the most basic, quintessential question about the flow of these fluids has remained unanswered: it is generally impossible to predict their flow resistance, and it is even unclear why most fluids shear thin, whereas only some shear thicken. Understanding shear thickening, i.e. the increase of the apparent viscosity of materials with increasing flow rate, is thus an important issue in complex fluids with in addition a strong impact on energy consumption in industrial processes [1]. It is observed in dense colloidal suspensions [1,2], where it has been related to the formation of dense clusters of particles [2][3][4]. The viscosity rise with the shear rate is then usually reversible (it is a steady-state property), continuous, and is sharper at higher volume fractions [2][3][4]. For colloids, the competition between shear-induced cluster formation and Brownian motion that homogenizes the suspensions then naturally determine a critical shear rate for the onset of shear thickening.

As Brownian motion is absent in pastes made of large particles, the sharp shear thickening transition observed in for instance cornstarch suspensions [5] is highly surprising. In fact, the conditions of emergence of shear thickening in non-Brownian suspensions remain ill-characterized: in some systems, thickening was observed at low shear rates [1,[5][6][7], while in others no shear thickening (only viscous behavior) is observed whatsoever, even close to jamming [8][9][10]. Up to now it is thus impossible to predict whether a given system will shear thicken or not.

In systems for which it is observed, a more pronounced shear thickening [6,7] is observed near jamming, similarly to colloidal suspensions, and is attributed to aggregation of hydroclusters into a percolating network [4,11]. However, one should question whether these observations of sharp and discontinuous shear thickening reflect an intrinsic (local, steady-state) property of materials. For example, an important effect of confinement on shear thickening was recently evidenced in both colloidal [12] and noncolloidal [5] suspensions. In channel flows [12], the macroscopic shear thickened state was shown to form a plug flow and was not observed in large channels, which shows that local observations are crucial to get a better insight into shear thickening.

In this Letter, we address these puzzles by studying the emergence of shear thickening in the simplest of systems: model density-matched suspensions of non-Brownian particles in water. We use a wide gap Couette geometry to avoid confinement effects, and we access the intrinsic material behavior by measuring the local flow properties and particle concentration using Magnetic Resonance Imaging (MRI). In steady state, we show that the material is heterogeneous, and that the local rheology presents a continuous transition, from a viscous to a shear thickening, Bagnoldian, behavior (shear stresses proportional to the shear rate squared) at any fixed volume fraction, as predicted by a scaling analysis. The heterogeneity is shown to result from an unexpectedly fast migration of grains during the transient, which is attributed to the emergence of the Bagnoldian rheology. The migration process is accompanied by macroscopic transient discontinuous shear thickening, which is thus not an intrinsic property of granular suspensions.

Materials and methods -We study dense suspensions of noncolloidal monodisperse spherical particles immersed in a Newtonian fluid. We use polystyrene beads (diameter 40 μm, polydispersity <5%, density 1.05 g.cm -3 ) suspended in aqueous solutions of NaI to match the solvent and particle densities; the solution viscosity is 1 mPa.s. The density matching ensures that there are no gravity-induced contacts [13] and that the only source of normal stresses is shear [14,15]. We focus on results obtained at a 59% mean volume fraction as experiments at other high volume fractions show similar features. The material behavior is studied with a wide-gap Couette rheometer (inner radius: 4.1 cm; outer radius: 6 cm; height of sheared fluid: 11 cm) inserted in a MRI scanner, allowing us to access local velocity and particle volume fraction profiles in the flowing sample [9,16,17]. Sandpaper is glued to the walls and there is no significant slip on the velocity profiles. The inner cylinder velocity is controlled, and we record the torque exerted by the material on the inner cylinder with a Bohlin rheometer.

Macroscopic vs. local behavior -We first focus on the macroscopic behavior (Fig. 1a). The torque T values measured during a slow ramp in rotational velocity Ω (logarithmic ramp, 45 s/decade of shear) are shown as black circles. At the beginning of the ramp, T increases linearly with Ω, as expected for a homogeneous, Newtonian suspension [9,10]. Around Ω c ≈ 2.5 rpm (corresponding to a low average shear rate of 0.6 s -1 ), T presents a sudden increase -by a factor of 20 -which is a usual signature of "discontinuous" shear thickening, following which T continues to increase with Ω, yet at a slower pace. At the end of the ramp, Ω is kept constant and T reaches a stationary value. Subsequent slow changes (down or up) of the rotational velocity then drive the system reversibly along the curve T(Ω) in open symbols, i.e. the system has reached steady state. This curve presents very smooth, moderate, "continuous" shear thickening. We now turn to local measurements. During the initial increasing velocity ramp, at a low rotational velocity, Ω=2 rpm (<Ω c ), the density ( ) R φ remains uniform (Fig. 1b) while the flow (Fig. 1d) is homogeneous-there is no jammed region [5]. The velocity profile V(R) closely matches that of a Newtonian fluid (Fig. 1d) [9,18], consistent with the initial linear behavior of the torque. Shortly after shear thickening occurs (here at Ω=9 rpm) we observe that the material has become strongly heterogeneous (Fig. 1b) while the velocity profiles present a jammed region near the outer cylinder. This change is irreversible: the density profiles subsequently remain similar, even when Ω is increased further and, later, decreased below Ω c .

Clearly, the discontinuous shear thickening observed during the initial up-ramp is a transient phenomenon associated with a large-scale reorganization of the material, which involves shear-induced migration from low to high shear zones. While migration is expected in dense suspensions [9,[19][20][21], it is particularly striking here that the change in φ(R) occurs over a very short time interval, corresponding to a small total strain of order 100. Such a rapid migration is a puzzle as it is not predicted by classical theories [9].

Constitutive behavior -We now analyze the steadystate behavior. We first note that density and velocity profiles are in steady state whenever the torque is. Moreover, while φ(R) is Ω-independent in steady state, the dimensionless velocity profiles V(R,Ω)/V(R i, Ω) measured at various Ω do not superpose, implying that the local behavior is not simply viscous [9,18]. Finally, the flow is always strictly localized: for all Ω, there is a jammed region beyond a critical radius R m =5.7 cm; this corresponds to density threshold φ(R)>φ m ≈60.5% above which the material is jammed [9].

The material and flow being heterogeneous, macroscopic torque measurements T(Ω) are not sufficient to infer the intrinsic constitutive behavior, in particular the stress/strain-rate relationship in the shear thickening regime. This intrinsic behavior can however be obtained using our local measurements. The key point [9] is that the steady-state density profile ( ) R φ is independent of Ω; a change of variables can then be performed between radius R and ( ) R φ . In addition, the stress profile is prescribed by momentum balance 2b). Such a transition has been predicted to be a generic property of noncolloidal suspensions on the basis of theoretical dimensional arguments [22,23]. The 2 γ scaling signals a regime where particle inertia dominates over viscous forces [8,23], leading to a behavior analogous to that of dry granular materials (it needs not be associated with collisiondominated flows as Bagnold suggested [START_REF] Bagnold | Proc. R. Soc. A[END_REF]). It is particularly striking that inertial scaling arises in our dense, highly damped, suspension, with particles of size of only ~10 μm. Moreover, the critical shear rate c γ (i) is rather low (of order 1 s -1 ), (ii) vanishes almost linearly as the volume fraction tends to φ c ≈60.5%, which (iii) is identicalwithin the experimental accuracy -to the threshold φ m at which the material jams. Viscous/inertial transition -To understand what controls these scaling regimes, following [22,23], we write Newton's equations for a set of rigid particles. For the particles' centers of mass r i these read:

2 2 d d visc i i j i j m r t F F = + ∑
where ij F denote rigid contact forces and visc i F hydrodynamic forces which we suppose linear in terms of all the velocities entering into the problem. No other force is supposed to be involved. The key of the analysis is to remark that rigid forces ij F do not introduce, by definition, any force or length scale [22,23]. Two limiting cases can then be identified: "viscous" (V) when viscous forces are dominant over grain inertia 0

visc ij i j F F = + ∑
, and "inertial" (I) when grain inertia is dominant to the linear vanishing of ( ) c γ φ , and hence permits this transition to take place at low strain rates. Accelerated shear-induced migration -We now show that this transition explains the sudden migration associated with the macroscopically observed transient discontinuous shear thickening (Fig. 1a). In viscous suspensions, shear-induced migration is usually thought to be negligible when small particles are involved. Indeed, the typical strainscale for migration is very large: it is expected to be rate independent and to scale as ( ) [19][20][21], leading to an expected strain of order 50000 [9], more than 500 times higher than what we observe here at the onset of shear thickening. Our observations may be understood as a strong enhancement of migration kinetics in the inertial regime.

o i R R a ∝ - [9,
Within the framework of two-phase models, migration is driven by gradients of internal normal stresses within the particle network p ii Σ (not the total stress), and requires the fluid to filter through the granular phase to compensate for the local changes of packing fraction [20,21]. This filtration process exerts an average hydrodynamic drag U ∝ on the particle network, with U the average filtration velocity. The balance between these two effects controls the migration/filtration rate, leading to p ii U ∝ ∇Σ . When injected in a mass conservation equation

( ) t U ρ ρ ∂ ∂ = -∇
, this leads to a diffusion equation for the particle density ρ [20,21]. The local particulate stress p ii Σ entering this analysis is expected to display local viscous or inertial scaling over very short strainscales [22,[START_REF] Peyneau | [END_REF][START_REF] Da Cruz | [END_REF] compared to those of the migration process. If p ii γ Σ ∝ in the whole system, the timescale of migration scales as 1/γ and migration is controlled by a (large) rate-independent strainscale, which is the classical result [19][20][21]. Strikingly, the same analysis performed in the inertial regime now yields an unexpected 1/γ strainscale for migration: this explains why migration is much accelerated, and manifests itself abruptly when entering the inertial regime.

To confirm this analysis, we have studied the migration kinetics at constant Ω 's, starting each time from a homogeneous state. We define the critical strain migr ( ) Γ Ω as the macroscopic strain above which the instantaneous volume fraction profile matches the steady state one within experimental uncertainty (0.2 %). Fig. 1c shows migr ( ) Γ Ω

vs. Ω : it decreases strongly with Ω . Although the exact kinetics results from a complex history (as both γ and φ change locally in time), the asymptotic 1/γ decay predicted by the above scaling analysis in the Bagnoldian regime is roughly consistent with our observations (see line in Fig. 1c). Migration theories based on normal stresses [20,21] are thus shown here to be more generally applicable than diffusive theories [19].

Macroscopic shear thickening -To summarize, we propose the following scenario: (i) the intrinsic behavior of dense noncolloidal suspension presents a continuous transition at low strain rates from a viscous to a shear thickening, Bagnoldian, rheology characterized by shear stresses 2 γ ∝ ; (ii) in the Bagnoldian regime, a very fast particle migration then occurs towards low shear zones; (iii) the interplay between flow and migration shows up as a sharp shear thickening of the transient macroscopic stress.

Note that we checked the robustness of the above scenario by verifying that it is not specific to wide gap Couette geometry. We have studied the behavior of the same material in a parallel plate geometry and have observed similar phenomenology, i.e. transient discontinuous shear thickening associated with fast migration due to the 2 γ ∝ scaling of stresses.

We here contrast our findings with those of Brown and Jaeger [7] in a similar system (spheres in a Newtonian fluid). On the basis of steady-state macroscopic measurements, they find a transition between a roughly viscous and a shear thickening regime Here, we find instead that the local (intrinsic) rheology shows only viscous or Bagnoldian scalings (i.e. n=0.5) even when m φ φ → . We insist that in steady state such systems are heterogeneous, and that consequently the macroscopic stress/strain-rate relationship cannot be directly related to the local constitutive behavior, in particular in the shear thickening regime.

Let us finally note that our mechanism may also be at work in Brownian suspensions, in competition with or as an alternative to hydrodynamic clustering [2]. It is compatible with the reversibility of the shear thickening transition usually observed in Brownian suspensions: migration is indeed expected to be reversible due to the osmotic pressure. It is therefore particularly striking that our mechanism leads to a constant critical shear stress, exactly what is observed for colloids [2].

Fig. 1 : 2 /

 12 Fig.1: a) Torque vs. rotational velocity Ω when shearing a 59% suspension: increasing velocity ramp (filled circles) and stationary state (empty circles); the line is a viscous law. b) Volume fraction vs. radius R, for various Ω, during an increasing velocity ramp. c) Critical strain needed to complete migration in a 59% suspension vs. Ω; the dotted line indicates the 2 / γ γ transition on the 59% suspension. The line is a 1/γ scaling. d) Dimensionless velocity vs. R (empty symbols: before migration, filled symbols: various Ω after migration, line: Newtonian velocity profile).

  shear rate can be extracted from the velocity profile V(R) via: stress/strain-rate curve ( , ) τ γ φ -at fixed and well-defined density φ -is then obtained by collecting all measurements of local stress ( ) R τ and shear rate ( ) R γ for a fixed R and varying Ω. The results of this local analysis (Fig. 2a) show that, for a fixed volume fraction φ, a clear transition from a occurs at a critical shear rate ( ) c γ φ (Fig.

  ), with identical scaling with γ for all components of the stress tensor as rigorously shown in[23]. The full problem then reduces to (V) (resp. (I)) at low (resp. high) γ , which explains the existence of a crossover between the two simple scaling regimes τ γ ∝ (viscous) be so low and vanishes precisely at φ m . In the viscous (V) and inertial (I) regimes, the stresses are respectively of the form ( ) ρ and d are the particle density and diameter, and o η is the interstitial fluid viscosity. The cross-over between the viscous and inertial regimes is found by equating the two expressions for the stress, in the literature[START_REF] Da Cruz | [END_REF]27], this equation explains, as observed, that ( ) c γ φ vanishes (i) linearly with φ, (ii) at the jamming packing fraction φ m . Moreover, the crossover that c τ should indeed be independent of volume fraction. Although our stress measurements are not sufficiently accurate to assert that ( ) cst c τ φ ∼ , we then note in Fig.2athat, indeed, in the experiments c τ does not vary much. We finally conclude that it is the difference in singular behavior of the inertial and viscous

We thank Daniel Lhuillier for enlightening discussions on shear-induced migration.


E Deletombe
email: eric.deletombe@onera.fr

J Fabis

J Dupas

9-13 mai 2011, Presqu'île de Giens (Var) Vulnerability of fuel tanks with respect to hydrodynamic ram pressure -Experiments and FE modellings

Keywords: HRAM, Vulnerability, fuel tanks, 7, 62 mm

A qualitative and quantitative analysis is proposed to define relevant equivalent load functions to be used for non linear structural analysis of aircraft fuel tank subject to 7,62 mm hydrodynamic ram effects. For the qualitative point of view, the main conclusion is that any load case that would only be based on peak pressure analysis would underestimate the final damage level in the structure : both the "acoustic" wave and the cavity in the wake of the 7,62 mm bullet potentially generate (at different time scales) more significant energy/momentum loads onto the structure walls.

Introduction

Reducing vulnerability of aircraft fuel tanks against the hydrodynamic ram pressure effects generated by a high speed/high energy 7,62 mm projectile is a complex fluid/structure mechanical problem. The projectile is meant to perforate the fuel tank wall with enough remaining energy to induce high energy and momentum transfer to the fluid, hence generating a short transient shock wave (1-3 ms), followed by a longer cavity bubble growth and collapse phase (10-30 ms) that would both possibly rupture the structure at different time scales [START_REF] Ankeney | Physical Vulnerability of Aircraft due to Fluid Dynamic Effects[END_REF]. The accurate non linear simulation of the full HRAM phenomenon would require the use of complex multiphasic fluids and nonlinear rupture material/assembly models : it first leads to very high computing costs, and also raises many numerical questions (code coupling, stability, mesh sensitivity, accuracy, …). Such a multi-physical and multiscale (time and space) modelling is still in the field of theoretical and numerical Research.

The engineering objective of the presented work is to rapidly and roughly discriminate what has to be modelled to estimate the induced level of HRAM damage in a tank structure, assess its residual strength, and possibly improve its design. The final goal is to use an explicit lagrangian FE code to study the structural response of the fuel tank, using generic equivalent "load" functions (mapping process) instead of modelling the fluid medium and its interaction with the projectile, first, and the structure, second. The first hypothesis that is made is that one would only consider minor changes of the fuel tank design, that would not affect the generic load functions in return (soft coupling of the fluid vs structure responses). There are two possibilities to do so : either we model the full complex problem then extract the load functions that are computed at the interface between the fluid and the structure, which is another ongoing Research work at ONERA/DADS/CRD, or we use experimental results (and perform different tests for each specific HRAM configuration : projectile speed and energy, fuel tank geometry and overall stiffness) to get such generic loads.

The present work concerns the second approach. The case of a 7,62 mm ballistic impact is studied where the HRAM scenario differs from cases when the projectile has a constant shape with no cavity "bubble" but only a wake being created in the liquid [START_REF] Nishida | Experimental study of perforation and cracking of water-filled alumi-num tubes impacted by steel spheres[END_REF]. The tests which are documented in the paper have been performed at ONERA/DADS/CRD in the frame of the RTP3.32 EUCLID project, funded by the French MoD (DGA). Note that the general experimental observation are very similar to those described in [START_REF] Varas | Experimental analysis of fluid filled aluminium tubes subjected to high velocity impact[END_REF] (to see for more bibliography), where a 12,7 mm bullet was used.

General considerations to derive Vulnerability requirements

In the field of dynamic strength of structures, the authors consider that the main point to care about is the energy (rate) transferred to the targeted structure that has to cope with. For instance, the structure "adaptation" might be kinetic (fragments projection), by heat or by irreversible solid deformation, damage and rupture. The studied mechanical system in which energy conservation principle is applied, is made of different "sub-systems" (7,32 mm projectile, cavity bubble, liquid medium, and solid structure). For each sub-system, different physical phenomena are involved, and ranked, according to the energy amount they would dissipate, store and release during the HRAM event. Some of them can be evaluated from the available test data, but others cannot. So this part of the paper focuses on general theoretical considerations. It is not here question of performing any formal thermodynamical analysis of the problem (are we considering an isenthalpic, isentropic, isothermal, adiabatic, etc, process when studying our HRAM problem ?). We only propose to perform a coarse energy evaluation exercise on the open mechanical system with is constituted by : the cavity bubble, the liquid mass, and the solid structure. We will neglect any cavitation at the interface between the structure and the liquid. We assume that our fuel tank is at equilibrium at initial time (before impact), and that the different initial potential energies in each sub-system will not be released into damage increase during the HRAM event. The energy which is added in the system (part of the initial kinetic energy of the bullet) will be distributed between different energy parts :

S E ∆ System Total Energy Increase ) ( B E f Bullet Energy (Source) sC E ∆ Sub-system/Cavitation Bubble Energy sL E ∆ Sub-system/Liquid Energy sT E ∆ Sub-system/Tank Energy
During the process we suppose that :

) ( B S E f E = ∆
, and

sT sL sC B S E E E E f E ∆ + ∆ + ∆ = = ∆ ) ( (1)
We consider in the present work that the way the transferred energy (and momentum) are dealt with by the structure is properly solved by the explicit FE code that will be used in the end (wave propagation, momentum, material behaviours, etc). Then the structural problem may be written in terms of energy as :

sL sC B sT E E E f E ∆ - ∆ - = ∆ ) ( (2)
Each sub-system energy variation is divided into different components, for instance : the kinetic energy of the sub-system, the potential (internal stored) energy of the sub-system, and the dissipated energy in the sub-system. A number of assumptions is made to deal with unknown quantities, hopefully in a conservative way with respect to the ultimate goal of the study : the damage prediction in the structure.

Assumptions concerning the 7,62 mm bullet sub-system energy

Different assumptions are made concerning ) ( B E f : 1) as we cannot measure precisely the bullet velocity and mass just after it has perforated the front wall of the fuel tank, we will assume the initial mass and velocity is unchanged,

2) though we know that the 7,62 mm will tear apart and probably heat when it tumbles across the fluid, we will not consider the associated dissipated energy because we do not know their value,

3) we will neglect other kind of potential energy (such as gravity).

Then we assume that all the (initial) kinetic energy of the bullet sub-system is progressively transferred into the other sub-systems energies, and that this approximation should maximise the energy transferred to the composite fuel tank structure to damage it. The amount of transferred energy that is considered in the present study is shown on fig. 1 : the initial velocity of the 9 g bullet is about 800 m/s before impact, and the tank dimensions about 300 mm in height (bullet trajectory), 550 mm in width and 650 mm in depth. It is fully filled of water. The bullet velocity was recorded/evaluated using high speed digital cameras (two walls of the caisson are PMMA windows).

0,00E+00 5,00E+02 1,00E+03 1,50E+03 2,00E+03 2,50E+03 3,00E+03
3,50E+03 0 0,0001 0,0002 0,0003 0,0004 0,0005 0,0006 0,0007 0,0008 0,0009 0,001

Aluminum Caisson

Assumptions concerning the wake and cavity behind the bullet

The cavity bubble (see fig. 2) develops long after the initial pressure pulse (shock wave) has vanished. We will consider that :

1) the kinetic and potential (gravity) energy of the water vapour in the cavity are unimportant compared to its latent energy (energy stored in physico-chemical change),

2) any calorific energy created during the HRAM process will never reverse back into a mechanical energy that could notably load the fuel tank structure (and we do not know how much energy it could be), so we will consider that temperature is constant (ambient) in the cavity during the HRAM process,

3) no viscous energy is dissipated in the liquid vapour. The main energy to be considered for the cavity bubble then corresponds to the latent heat. An important question is : what is its amount (in fact, some air is entrapped in the wake of the bullet, and the gaz is an air/water vapour mixture) ? If we consider that the cavity is only full of water vapour (pure cavitation phenomenon), then the latent heat is the amount of energy stored (and released later) by water during its change of state without changing of temperature, meaning a phase transition :

v C C L V Q ρ = (3)
where Q is the amount of energy stored (and released) during the change of phase of the water, V C the volume of the cavitation bubble, C ρ the vapour density, and L v the specific latent heat of vapourization. To calculate the latent heat, we need the volume of the cavity bubble V C . This volume is experimentaly estimated from high speed camera pictures. The point then is to choose values for L v and C ρ . Those values would depend on the thermodynamical state of the vapour in the cavitation bubble : temperature and pressure, that we do not know and which possibly vary during the test. Different possible values are given hereafter for the thermodynamical coefficients : In the end, latent heat evolution can be calculated according to thermodynamical state assumptions. If we consider an isothermic cavitation process in the liquid under 1 Patm, 293 K conditions, then :

        - = T T R ML sat v e P P 1 1 0 0 = 0.0204 Patm ( 4
)
with T 0 =373 K, M = 0,018 kg/mol, L v = 2,455×10+6 J/kg, R = 8,31447 J/K/mol, P 0 = 10+5 Pa. Then the maximum latent energy reaches about 200 J (0,02 Patm), which is the estimation that we will take (no experimental information concerning the air/vapour mixture, the temperature and pressure level in the cavitation bubble during the test). This latent energy will be stored (no dissipation assumed) then released during a 8 ms long process (observed time of collapse of the cavity bubble).

Assumptions concerning the liquid medium (water)

Concerning the liquid medium : 1) we assume that its temperature is constant, and that it behaves as a non viscous fluid. No energy is then supposed to be dissipated in this sub-system, 2) the potential (gravity) energy rise (part of the water is moved upward) is small compared to other components, 3) its mass and volume are supposed to be constant (incompressible fluid, vapour mass very small -x1000 times less -compared to liquid mass in the fuel tank), which means that the main energy components will consist in kinetic and acoustic contributions. The next paragraph is specifically dedicated to the evaluation of the energy transfer function by the liquid medium onto the structure.

Evaluation of kinetic energy and momentum in the liquid medium

A number of simplifications or assumptions are made (sometimes because of unknown quantities) to evaluate, in the one hand, the purely acoustic energy (compressible liquid) and, on the second hand, the kinetic energy (cavity expansion, uncompressible hypothesis) in the liquid medium (water). Two kinds of experimental data are available, that concern first the (assumed spherical [START_REF] Lecysyn | Experimental study of hydraulic ram effects on a liquid storage tank : Analysis of overpressure and cavitation induced by a high-speed projectile[END_REF]) pressure shock wave (to derive acoustic energy), and the cavity bubble expansion velocity (to derive liquid velocity fields, then momentum and kinetic energy under uncompressible hypothesis).

Acoustic energy

The pressure gauge is located somewhat in the centre of the tank, at about 0.2 m away from the bullet impact point.on the entry wall of the structure : For a compressible fluid, a pressure wave induces a density change in the liquid medium (hence and optical indices change) which also means that the fluid particles concentrate thus move when the transient pressure pulse comes through. In an adiabatic process, the acoustic energy corresponds to the internal energy stored in a volume element V when the initial pressure P 0 increases to pa = P 0 + p. For a perfect fluid, the internal energy (small perturbations hypothesis) that is transported by pressure waves is known to be made of two components : a kinetic and a purely acoustic one : [START_REF] Townsend | Failure of fluid filled structures due to high velocity fragment impact[END_REF] with p a the acoustic pressure, u the fluid velocity, c the sound celerity, and ρ the fluid density (initial and actual values). So we consider a spherical pressure wave radiating from the impact point of the bullet on the tank wall, and inducing an approximated triangular p(t) pulse shape (see fig. 4), that can be correlated in space p(r), knowing the sound celerity and the 0,2 m position of the pressure gage from the impact point. The purely acoustic energy can then be coarsely evaluated from this 5 MPa peak / 0.150 ms long pressure pulse, and reaches a maximum of about 100 J. If we now consider the kinetic component : the increase of density for an average 2,5 MPa pressure peak is only 0,1% of the initial density, but about a 0,3 m 3 volume (then 300 kg) of water is concerned by the pressure pulse, which means that a 3 m/s average fluid velocity increase (in the half-spherical water ring that comes trough the pressure gauge) would be enough to "store" the remaining energy of the bullet initial kinetic energy (once the latent heat and purely acoustic energy is taken out). Unfortunately, it was not possible to measure the fluid velocity field with the experimental techniques that were available at ONERA/DADS/CRD (the particle image velocimetry technique could help, that is used by ONERA/DAAP Department). Nevertheless, this point is of first importance, since this kinetic energy would act on the structural wall thanks to momentum principle.

∫∫∫ + = V a dV c p u I 2 2 2 2ρ ρ , with : ) ( 0 2 ρ ρ - = c p a

Liquid medium kinetic energy due to wake/cavity growth

Another water kinetic energy component is due to the 0.125 m radius / 10 ms cavitation bubble growth and collapse, that can coarsely be calculated. Looking at the high speed camera pictures, we assume as a first order approximation that the volume of water is a sphere of radius

) ( 2 t R with a gaz bubble of radius ) ( 1 t R
inside, both in radial expansion (u= r& ), and that ρ is constant (uncompressible hypothesis). Then, the associated variation of kinetic energy in the liquid can be expressed as :

∫ ∫ ∫ = ∆ ) ( ) ( 2 2 2 0 0 2 1 sin 2 1 ) ( t R t R k sL dr r r d d t E & π π θ φ φ ρ (6)
The conservation of mass of the liquid (i.e. dm/dt=0) gives :

dt t dR t R dt t dR t R ) ( ) ( ) ( ) (
∫ ∫ ∫ = ∆ ) ( ) ( 2 2 2 0 0 2 1 sin 2 1 ) ( t R t R k sL dr r K d d t E π π θ φ φ ρ (9)         - = ∆ ) ( 1 ) ( 1 ) ( 2 ) ( 2 1 2 t R t R t K t E k sL πρ , with : 2 1 1 2 ) ( R R r r t K & & = = (10)
The conservation of volume (incompressible fluid assumption) is now written to discard R 2 from the above equation :

) ( 3 4 ) ( 3 4 3 1 0 3 2 t R V t R π π + = (11)
The kinetic energy due to the cavity bubble growth (and collapse) can then be calculated from the experiment : a polynomial interpolation of R 1 (from test data) is found and derivated to get The kinetic energy largely overpasses the initial kinetic energy of the bullet at the beginning of the penetration phase, which can be easily attributed to experimental "errors" on the cavity "bubble" (which is more a wake than a bubble) radius measurements (see fig. 2). Once the cavity bubble dimension turns to be more measurable (around 0,001 s) the energy amount gets back around the initial kinetic energy of the bullet (which is at that time almost integrally transferred to the system, see fig. 1). Note on fig. 7 that a maximum 3 m/s value would then be obtained as an order of magnitude at the composite wall distance.

Conclusions

HRAM is a very difficult phenomenon to observe and quantify. Most of the time, the structural response of containers is measured (strain gauges, accelerometers, post-mortem analysis) which gives an indirect access to the loading they have to cope with. To study the liquid behaviour, traditional pressure transducers are often immersed in the tanks, but recent works propose to use high speed digital image cameras to observe and measure the dynamic evolution of the geometry of the wake behind, and the cavity around, these projectiles. The originality of the Research compared with the few papers found in the open literature, is that most of them consider constant projectile in shape and then drag coefficients, and only the very first stage of the penetration process (less than 1 ms). The present paper studies two firing tests performed in "confined" water containers, using a deformable NATO 7.62 mm bullet which tumbles and is torn apart during the water penetration process.

A qualitative and quantitative analysis is then proposed to define which kind of simplified equivalent load functions should be applied for FE non linear structural analysis of aircraft fuel tank subject to 7,62 mm hydrodynamic ram effects (subsonic projectiles, large fuel tanks). For the qualitative point of view, the main conclusion seems to be that any load case that would only be based on peak pressure analysis would underestimate the damage level in the structure (except at the entry point where the local pressure might be higher than the material strength [START_REF] Lecysyn | Experimental study of hydraulic ram effects on a liquid storage tank : Analysis of overpressure and cavitation induced by a high-speed projectile[END_REF], particularly for supersonic projectiles [START_REF] Townsend | Failure of fluid filled structures due to high velocity fragment impact[END_REF]). Then, both the acoustic pressure and the cavity growth that forms in the wake of the projectile could potentially transport a most significant kinetic energy trough the fluid, then momentum onto the structure. The share of kinetic energy due to acoustic and cavity mechanisms was not possible to get from the present work. Both of these kinetic energies will act (impulse, momentum) onto the structure walls, but at different time scales : the pressure wave crosses the tank within 0,2 ms and the cavity bubble expansion lasts around 10 ms (with a complex growth/collapse/rebounce scenario), and would both possibly constitute serious threats for the fuel tank integrity.

As it is not possible to share the energy between these two mechanisms, the final proposal would then be to model two "limit" cases with all the energy transported either by the acoustic wave (compressible fluid), on the one hand, or by the "cavity growth", on the other hand. Both cases can be used to derive impulse or momentum fields (pressure or velocity) that could then be applied onto the tank walls without modelling the bullet/fluid systems : this possibility will be studied at ONERA/DADS in the frame of the EDA BaTolUS project.

FIG 1 :

 1 FIG 1 : Evolution of Kinetic Energy (J) vs Time (s) of the 7,62 mm bullet in the tank

FIG 2 :

 2 FIG 2 : Pictures (0/0,2/2,5/4,5 and 6,5 ms) in caisson, with high speed camera (10+3 im/s)

FIG 3 :

 3 FIG 3 : Order of magnitude of cavitation bubble latent energy (J) vs Time (s) for different level of external pressure (blue 0,02 Patm, pink 1 Patm, and yellow 30 Patm), wrt the observed volume of the cavity bubble

FIG 4 :

 4 FIG 4 : Relative pressure pulse function P-Patm (bar) vs time (µs)

FIG 5 :

 5 FIG 5 : Cavitation Bubble Radius R 1 (m) vs Time (s) in composite caisson.

  FIG 6 : Evaluation of water kinetic energy (test and interpolation) and latent heat

  FIG 7. : Evaluation of the liquid velocity at the tank wall location, based on the measurement (and interpolation) of the cavity bubble dimension (uncompressible assumption).

The authors would like to thank the French MoD and DGA, for their financial support in the frame of the EUCLID RTP3.32 and EDA BaTolUS projects. The caisson specimens were designed and manufactured by EADS-IWF in the frame of the EUCLID RTP3.32 project.


Marc Fournier

Benoît Combès

Neil Roberts

Simon S Keller

Timothy J Crow

William D Hopkins

Sylvain Prima

Surface-based method to evaluate global brain shape asymmetries in human and chimpanzee brains

Keywords: Chimps MRI brain segmentation, Brain global shape computation, Surfaced-based bilateral brain asymmetries, Axis-wise asymmetry components, Chimp and human asymmetries comparison

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

INTRODUCTION

Anatomical asymmetries of the human brain [START_REF] Geschwind | Human brain: left-right asymmetries in temporal speech region[END_REF][START_REF] Lemay | Functional and anatomical asymmetries of the human brain[END_REF][START_REF] Toga | Mapping brain asymmetry[END_REF] have been known since the descriptions of asymmetries of the lateral Sylvian fissure by Eberstaller and Cunningham in the nineteenth century, but were largely forgotten until the rediscovery of the asymmetry of the planum temporale by Geschwind and Levitsky in 1968. Left larger than right planum temporale, as well as protruding right frontal and left occipital lobes (often termed the petalia or Yakovlevian torque), are the most striking asymmetrical features in most human brains. These asymmetries have been hypothesized to be related with brain functions, and are most often met in right-handed people whose language areas are mainly located on the left side of the brain. These patterns of anatomical/functional asymmetries have long been thought to be typical of human beings (Homo sapiens) [START_REF] Holloway | Brain endocast asymmetry in pongids and hominids: some preliminary findings on the paleontology of cerebral dominance[END_REF]. There has been much debate recently on this puzzling question [START_REF] Crow | Directional asymmetry is the key to the origin of modern Homo sapiens (the Broca-Annett axiom): A reply to Rogers' review of The Speciation of Modern Homo Sapiens[END_REF], and some have claimed that these asymmetry patterns are actually shared by other great apes [START_REF] Gannon | Asymmetry of chimpanzee planum temporale: humanlike pattern of Wernicke's brain language area homolog[END_REF][START_REF] Cantalupo | Asymmetric Broca's area in great apes[END_REF]. Answering this key question could have important implications in the understanding of language origins and in what makes the humans unique. Some techniques have been devised for the assessment of brain asymmetries in brain MR images of humans. Three approaches have been especially used. Most asymmetry studies are based on voxel-based morphometry (VBM) and variants [START_REF] Good | Cerebral asymmetry and the effects of sex and handedness on brain structure: a voxel-based morphometric analysis of 465 normal adult human brains[END_REF][START_REF] Watkins | Structural asymmetries in the human brain: a voxel-based statistical analysis of 142 MRI scans[END_REF][START_REF] Herve | Handedness and cerebral anatomical asymmetries in young adult males[END_REF]. Two other alternative techniques have also shown promising results, based on deformation-based morphometry (DBM) [START_REF] Thirion | Statistical analysis of normal and abnormal dissymmetry in volumetric medical images[END_REF][START_REF] Prima | Computation of the midsagittal plane in 3D brain images[END_REF] and low-dimensional analysis of asymmetries (LowD) [START_REF] Barrick | Automatic analysis of cerebral asymmetry: an exploratory study of the relationship between brain torque and planum temporale asymmetry[END_REF]. However, to our knowledge, only voxel-based morphometry technique [START_REF] Hopkins | Gray matter asymmetries in chimpanzees as revealed by voxel-based morphometry[END_REF] and related volume-based analysis [START_REF] Avants | Lagrangian frame diffeomorphic image registration: Morphometric comparison of human and chimpanzee cortex[END_REF][START_REF] Keller | A comparative magnetic resonance imaging study of the anatomy, variability, and asymmetry of broca's area in the human and chimpanzee brain[END_REF][START_REF] Hopkins | Volumetric and lateralized differences in selected brain regions of chimpanzees (pan troglodytes) and bonobos (pan paniscus)[END_REF][START_REF] Schenker | Broca's area homologue in chimpanzees (pan troglodytes): probabilistic mapping, asymmetry, and comparison to humans[END_REF] have been applied to chimpanzees. In this paper, we propose a new automated surface-based technique for the analysis of global shape asymmetries in humans and chimpanzees. First, the brain is segmented (and its hemispheres are separated) from high-resolution T1weighted 3T MR images, and its envelope is computed to assess its global shape. Second, asymmetries of this global brain shape are computed on an individual basis. A mean brain shape is then computed over a population of subjects, on which individual asymmetries are mapped for statistical analysis of the population asymmetries. Finally the output of these analyses on humans and chimpanzees are compared. Section 2 of this paper describes the material and methods: the segmentation pipeline, the brain global shape computation, the brain bilateral asymmetries evaluation, and the brain mean shape estimation. In Section 3 we present preliminary results on humans and chimpanzees obtained using the new technique, before concluding in Section 4.

MATERIAL AND METHODS

This section describes the subject databases and the methods designed to process the datasets. We use two databases of MR images, the first one of humans and the second one of chimpanzees, composed of 30 subjects each. We focus our study on a subgroup of the population represented in the databases in order to limit the inter-subject variability. We choose to analyze and compare humans and chimpanzees asymmetries on the right-handed young adult male subjects of both populations. For the two populations, this group shows significant asymmetry results. For the chimpanzees this group is composed of 7 subjects from 12 to 26 years old and for the humans the group is composed of 11 healthy control subjects from 21 to 35 years old.

Segmentation pipeline

There are many softwares available to process MR images in order to obtain human brain segmentation. Most of these tools use a human brain atlas to register the data in the segmentation process. Methods relying on human brain atlases do not work for chimpanzees because their brain anatomy is quite different from the human one. Figure 1 shows some of these differences on MR image slices. To our knowledge no chimpanzee brain atlas is available yet to use standard atlas-based methods for brain segmentation. While working with chimpanzee brain images, we need to define a pipeline without the use of an atlas. First we use BrainMask software (http://www.nitrc.org/projects/brainmask/) for the bias correction step followed by the initial brain mask evaluation for the whole brain. BrainMask provides an automatic and robust algorithm without atlas which works with chimpanzee brains. Then we use BrainVisa software (http://brainvisa.info/) to perform hemisphere split, surface extraction and brain alignment to the coordinate system axes based on anterior-posterior commissures landmarks. We use BrainVisa standard atlas-based process for human brains and its atlas-free feature for chimpanzee brains. In the atlasfree case for chimpanzees, we also use BrainVisa advanced correction and validation tools which provide an iterative threshold feature to obtain relevant results of hemisphere split.

Brain global shape computation

To alleviate the problems due to cortical variability, we base our asymmetry measures on a simplified yet accurate representation of the brain which closely follows the gyri without entering the sulci. We call this global shape the brain hull which we compute for each hemisphere obtained from BrainVisa output.

We first compute the 3D convex hull of each brain hemisphere vertices. This algorithm generates two meshes with few large triangles corresponding to coarse representations of each brain hemisphere which is not sufficiently accurate for our purpose. We subdivide these meshes to obtain a uniform triangulation with enough vertices to accurately describe the initial hemispheric brain meshes gyri positions. Then for each vertex of the subdivided convex hulls, we compute the minimal sphere centered on this vertex which intersects the boundary of the initial detailed hemispheric brain meshes. Finally we move the vertex in its opposite normal direction (towards the brain surface) on the minimal sphere surface.

This procedure illustrated in Figure 2(a) attracts vertices to the gyri positions and prevents them to penetrate inside the sulci, between neighboring gyri. It produces accurate hemispheric brain hulls shown in Figure 2

Brain bilateral asymmetries evaluation

The typical Yakovlevian torque human brain asymmetry of right-handed subjects is illustrated in Figure 3 as a combination of patterns which are exaggerated here to highlight this phenomenon. To quantify the magnitude of this torque we present a surface-based method to compute individual asymmetry maps (Section 2.3.1). These maps are defined over a vector deformation field that is decomposed in the three orthogonal axes. Then we explain how to put the individual asymmetry maps in a common coordinate system to perform statistical analyses (Section 2.3.2). The next two subsections provide specific implementation procedures on this strategy and additional methodological details can be found in [START_REF] Combès | Automatic symmetry plane estimation of bilateral objects in point clouds[END_REF][START_REF] Combès | Setting priors and enforcing constraints on matches for nonlinear registration of meshes[END_REF][START_REF] Combès | EM-ICP strategies for joint mean shape and correspondences estimation: applications to statistical analysis of shape and of asymmetry[END_REF].

Generation of individual asymmetry maps

We consider n meshes X 1 ,...,X n representing the n brains under study, independently for each population (human and chimpanzee separated). For each of these n meshes of a population, we propose to perform the following process.

• First step: Finding the symmetry plane.

We define the approximate symmetry plane P of a mesh X as the one best superposing X and its reflection S P (X) about P. For that, we consider X as a noised version of S P (X). This allows to consider each point xi of X as the realization of a random variable whose distribution is a mixture model composed of card(X) Gaussian laws N(S P (x k ),

σ 2 I), x k ∈ X.
Then we define the optimal reflection S P as a ML estimate, which is computed using an EM algorithm. The use of a multiscale strategy and of a truncated Gaussian kernel allows a fast, accurate and robust estimation of the plane [START_REF] Combès | Automatic symmetry plane estimation of bilateral objects in point clouds[END_REF].

• Second step: Registering X on S P (X).

Once P is known, we define the asymmetry field of X as (t(x) = T(x) -x) x∈X where T is the non linear transformation best superposing X and S P (X). For that, we consider T (S P (X)) as a noised version of X (here S P is known). Then we design a MAP approach by specifying i) a prior on T (a global affine model + a local order-one Tikhonov regularization on the linear component (t(x))_x), ii) a mixture density for the set T(S P (X)) and iii) a priori probabilities of matching between points of the two meshes (using geometrical descriptors almost invariant to the unknown transformation). Then, this MAP problem is solved using an EM algorithm [START_REF] Combès | Setting priors and enforcing constraints on matches for nonlinear registration of meshes[END_REF].

• Third step: Mapping the asymmetry field In order to provide specific results, we choose to analyze the asymmetry field by looking at its three components separately. This simply consists in projecting the asymmetry vector on the three coordinate axes and this yields the:

-postero-anterior component (difference in protrusion); -left-right component (difference in width); -head-foot component.

In the following, we call s i X the scalar asymmetry mapping at point x i of mesh X and S X the asymmetry map of X. Each subject has thus three asymmetry maps.

Generation of a mean shape & projection of the mapping

• First step: Computing the mean shape Individual asymmetry maps S X (X ∈ {X 1 ,...,X n }) have to be put in a common geometry to be compared. For this purpose, we compute the mean mesh M, defined as the mesh closest to all the meshes in the dataset (in a sense to be defined). In practice, we formulate the mean shape estimation as an iterative scheme that consists in the successive estimation of i) the mean shape and ii) the optimal similarity transformations between each mesh and the mean mesh [START_REF] Combès | EM-ICP strategies for joint mean shape and correspondences estimation: applications to statistical analysis of shape and of asymmetry[END_REF].

• Second step: Individual projection of the asymmetry mapping Together with the mean shape M, the previous algorithm provides us with the a posteriori probability that each point m i of M is matched with a point x j of X; we call this probability A ij (∑ j A ij = 1). This probability is used to project each individual asymmetry map on M. We call s i X→M the mapping information of X projected on point m i of M and define it as

s i X→M = ∑ j A ij s i X .
• Third step: Population statistical asymmetry maps At this point we now have a collection of n (×3) scalar asymmetry maps S X 1 →M ,....,S X n →M projected on a common mean mesh M. We compute a point-wise mean and t-test (corrected for multiple comparisons) on each of the three components of the asymmetry field over the population. The correction for multiple comparisons is performed as described in [START_REF] Nichols | Nonparametric permutation tests for functional neuroimaging: A primer with examples[END_REF] by i) fixing a supra-threshold p = 0.05, ii) building the permutation distribution of the maximal suprathreshold cluster size (by shuffling p-values over the mesh and computing the size of the largest cluster of significant vertices, nbTest=20.000 permutations) iii) considering the critical supra-threshold cluster size as the [α*nbTest]+1 largest value over the sampling distribution and iv) removing from the statistical map the clusters having a size smaller than the critical supra-threshold cluster size.

RESULTS AND DISCUSSION

Asymmetry maps interpretation

For a correct interpretation of the directional asymmetry maps results computed and presented in this section, one needs to understand how to read these maps. Figure 4 illustrates an axial view example of the asymmetry maps direction definition. Let X be an object (in green) and S P (X) its symmetrical image (in orange) according to the symmetry plane P. We compute the deformation superposing X on S P (X) illustrated with the arrows in Figure 4. The left-right component (x axis) of the deformation field (arrows a, b, c, d) are all positive deformations because they all point towards the positive x axis. The anterior-posterior component (y axis) of the deformation field has some positive values (arrows f, g) which point towards the positive y axis and some negative values (arrows e, h) which point towards the negative y axis.

The amplitudes of the deformation are illustrated using a color map (cold to hot type) as shown in Figure 4 color bar example. Asymmetry results use this color map over the brain mean shape of a population. This three components asymmetry maps definition is useful to highlight specific asymmetries in a chosen direction. However one needs to look carefully at the axis direction for asymmetry interpretation. Positive (resp. negative) values which lead to hot (resp. cold) colors do not necessary mean expansion (resp. shrinkage). Indeed, some deformations in Figure 4 (arrows e, f) fit with the intuitive interpretation of positive expansion and negative shrinkage, but other deformations (arrows g, h) are in the opposite direction. The convention for all results considered for the third component of the deformation field (not shown in Figure 5) is from Foot to Head. One can note that inverting the direction definition of a specific axis would result in inverting its respective red-blue asymmetry map pattern.

Humans and chimpanzees asymmetry results

Figure 6 shows the most significant views (in the right-left and posterior-anterior components) of the asymmetry maps and p-values for human results. . Each row shows one of the components and displays the region where the asymmetry is the most significant. The p-values maps are computed according to the t-test with the null hypothesis that there is no asymmetry at the significance level p=0.05 (corrected). Red color indicates with more regions with asymmetries. At the population level, the left occipital lobe appears to be wider and to have a posterior protrusion compared to the right occipital lobe. Asymmetry maps observations show a similar behavior of the right frontal lobe respectively to the left frontal lobe. However in frontal lobes the deformations are smaller and less equally distributed (Figures 6(a) and 6(c)) and the t-test (not shown) is not significant. These findings using our surfaced-based asymmetry method agree with the literature about the Yakovlevian torque. Figure 7 shows equivalent results as in Figure 6 but with the chimpanzees. Very similar observations about the width and protrusion patterns of the frontal and occipital lobes in humans can be seen throughout the chimpanzee results. However, two differences are observed mainly in the frontal lobes (Figures 7(a) and7(c)) where the greater values of deformations are not located at the tip of the lobes but moved laterally on each side. Moreover this time the most significant t-test results were observed in both frontal and occipital lobes of the right-left component (Figures 7(e) and 7(f)) compared to the test on the posterior-anterior component (not shown). An interesting observation is that the maximum deformation amplitude for each component is almost the same in chimpanzees and humans even if the human brain size and volume is greater. In that sense, this suggests the chimpanzee brain would be more asymmetrical than the human one. Nevertheless, according to our surfaced-based asymmetry method findings, it seems that the Yakovlevian torque pattern of the human brain is also present in the chimpanzee brain, at least in our study on this adult male right-handed subgroup. Figure 8 shows an interesting finding using our surfacebased asymmetry method on humans in the occipital lobe area with a significant foot-head (vertical) asymmetry. The left lobe seems to be greater in the lower vertical direction compared to its counterpart. This observation can be seen as an addition to the Yakovlevian torque where the left occipital lobe is already greater in the sense that it is wider and that it protrudes more in the posterior direction compared to the right lobe. This observation is not present at all in the chimpanzees and compared to the previous results, it is the first significant difference finding between the human and chimpanzee brains using our surfaced-based asymmetry method. We also notice that the maximum asymmetry (1.7mm) is smaller than the other components of the deformation field in humans. The equivalent maximum foot-head deformation in chimpanzees (1.3mm, not shown) is also significantly smaller than the human one, compared to the other components where maximum values were almost identical in both populations. Figure 9 shows how significant are the deformation fields in human and chimpanzee on the left hemisphere for the right-left component. It is interesting to see that our surface-based method results show very similar and significant specific asymmetry in the right-left component of the deformation field around Broca's area in both populations. However the amplitude of the deformation in the right-left component is small in that area, especially in humans where it is almost not visible on the asymmetry map when displayed at the regular full asymmetry range. The specific asymmetry observation in Figure 9 is only present in that component of the deformation field. For both populations there is no significant asymmetry in the other two components around Broca's area.

CONCLUSION

In this paper we proposed a new surface-based method to compute bilateral brain asymmetries with respect to the midsagittal plane. We used this method on both chimpanzee and human populations in order to compare their asymmetries. We adapted a human brain segmentation pipeline to work with chimpanzee MR images. We proposed a brain global shape mesh processing algorithm to compute a brain hull in order to alleviate the problems due to cortical variability. We designed a method for mid-sagittal plane computation and individual brain asymmetry evaluation and a method for the estimation of a brain mean shape used to map statistical asymmetries at the population level. Results are presented as specific asymmetry maps for each component of the deformation field between the brain and its mirrored image. Our method allowed recovering typical global human brain asymmetry patterns previously observed and to find another specific asymmetry component. Future work will involve more extensive studies on a larger database including comparison with other subgroups of both investigated populations and also with other great apes populations. We will design an automated comparison method between populations including an integration of asymmetry on specific anatomical regions for easier interpretation of the results.

Fig. 1 .

 1 Fig. 1. Human and chimpanzee brain differences. (a) and (b) axial views show the surrounding tissues of the brain are thicker in chimpanzee. (c) and (d) sagittal views show different relative orientations of the anterior-posterior commissures and the spinal cord. (a) and (b) show on the same scale that the chimpanzee brain size is smaller than the human one.

Fig. 2 .

 2 Fig. 2. Brain hull mesh processing illustration. (a) Brain hull algorithm computation with: V i the initial vertex position of the subdivided convex hull; V f the final vertex position of the brain hull; -n the opposite normal vertex displacement direction. (b) Brain hull meshes superposed on a human brain hemispheric surfaces. (c) A chimpanzee brain hemispheric meshes output from BrainVisa. (d) The chimpanzee hemispheric brain hulls computed.

Fig. 3 .

 3 Fig. 3. Global brain asymmetries in right-handed human subjects when viewed from above. (a) Anterior (resp. posterior) protrusion of the brain in the right as compared to the left frontal lobe (resp. in the left as compared to the right occipital lobe). (b) Leftward frontal lobe and rightward occipital lobe deviation of the inter-hemispheric fissure.

Fig. 4 .

 4 Fig. 4. Asymmetry maps direction definition. Orthogonal components of the deformation field which are in the same (resp. opposite) direction as their positive axis definition have positive (resp. negative) values. A color map is used to represent the signed scalar values of the deformation field.

Figure 5

 5 Figure 5 illustrates the typical asymmetry color map of the Yakovlevian torque in the case of right-handed human subjects. The right-left and posterior-anterior components of the deformation field are shown according to the axes positive direction definitions shown in Figure 5. These two positive axes direction conventions are also the ones considered for all human and chimpanzee results presented.The convention for all results considered for the third component of the deformation field (not shown in Figure5) is from Foot to Head. One can note that inverting the direction definition of a specific axis would result in inverting its respective red-blue asymmetry map pattern.

Fig. 5 .

 5 Fig. 5. Asymmetry color map schematization. Typical Yakovlevian torque asymmetry map for right-handed human subjects. (a) Right-left component of the deformation field positively defined from Right to Left. (b)Posterior-anterior component of the deformation field positively defined from Posterior to Anterior.

  Figures 6(a) to 6(d) have indication of the positive direction of the deformation field. Their legend color bars indicate the signed deformation for each component and green color means no asymmetry (deformation)

Fig. 6 .

 6 Fig. 6. Humans asymmetry maps and p-values. (a) Front view of right-left component. (b) Back view of right-left component. (c) Front view of posterior-anterior component. (d) Back view of posterior-anterior component. (e) p-values back view of right-left component (resp. to b). (f) p-values back view of posterior-anterior component (resp. to d).

Fig. 7 .

 7 Fig. 7. Chimpanzees asymmetry maps and p-values. (a) Front view of right-left component. (b) Back view of right-left component. (c) Front view of posterior-anterior component. (d) Back view of posterior-anterior component. (e) p-values front view of right-left component (resp. to a). (f) p-values back view of right-left component (resp. to b).

Fig. 8 .

 8 Fig. 8. Humans foot-head deformation field component. (a) Bottom view of asymmetry map. (b) Bottom view of pvalues. The color bar legend is for the asymmetry map. Pvalues have the same legend as previous figures (0; 0.05).

Fig. 9 .

 9 Fig. 9. Chimpanzees and humans right-left p-values. (a) Humans left side view. (b) Chimpanzees left side view. P-values have the same legend as previous figures (0; 0.05) for both populations.


Catherine Berho
email: c.berho@brgm.fr

Anne Togola

Benjamin Girardeau

Alain Saada

Jean-Philippe Philippe Ghestem

Monitoring of volatile organic pollutants in groundwater by passive diffusive bags (PDBs)

teaching and research institutions in France or abroad, or from public or private research centers.


Parallel simulation of three-dimensional convective mixing in long-term geological CO 2 storage in saline aquifers M. A. Sbai (1) (1) BRGM, Water Division, Orléans, France (a.sbai@brgm.fr / Fax: +33 2-38643446) Accurate numerical modelling of the effective dissolution rate of aqueous CO 2 solutions is important for analysis of long-term storage efficiency and risk assessments. Theoretical and field studies showed the possibility of switching from a diffusioncontrolled mode of dissolution to a natural convection mode which is much more efficient for accelerating the dissolution mechanism. However, most of the previous investigations have been restricted to two-dimensional cross-section models similar to Elder (1947) problem. In this work, we developed a high resolution parallel simulator to investigate these mechanisms with grid blocks needed for accurate modelling of convective mixing in realistic three-dimensional configurations. The simulator sequentially couples groundwater flow and CO 2 transport equations by density and viscosity contrasts. Density and viscosity of the aqueous CO 2 solution is predicted by accurate equations of state module validated against previously published experimental data for a wide range of temperatures and pressures of interest. Numerical experiments show that two-dimensional modelling over-predicts the time to onset natural convection and under-predicts the time to achieve ultimate dissolution of the CO 2 -free phase in a homogeneous porous medium. The role of aquifer heterogeneity in further enhancing CO 2 dissolution is illustrated for a generic case study. Heterogeneity leads to a disconnection of the otherwise continuous fingers to form CO 2 blobs and ganglia, increasing the effective surface area between CO 2 -rich blobs and formation waters. Results suggest that high resolution three-dimensional modelling should be guided in CO 2 storage projects to (a) estimate the onset of natural convection, (b) determine the size of convective instabilities when they occur, and (c) calculating the maximum mixing achieved.