domain: spi phys

text:
Valentin Joncquieres
email: joncquieres@cerfacs.fr

Francois Pechereau

Alejandro Alvarez Laguna

Anne Bourdon

Olivier Vermorel

A 10-moment fluid numerical solver of plasma with sheaths in a Hall Effect Thruster

Electric propulsion can reach higher exhaust velocities compared to chemical systems and thus result in lower propellant mass requirements. Among the different electric propulsion systems, Hall effect thrusters are used for spatial propulsion since the 1970s. However inside a Hall thruster, complex physical phenomena such as erosion or electron anomalous transport which may lower thruster efficiency and lifetime, are not yet fully understood. Thanks to high performance computing, numerical simulations are now considered for understanding the plasma behavior. With the renewed interest for such electric propulsion to supply small satellites, numerical solvers able to predict accurately the real thruster efficiency have become crucial for industry. This paper presents the approach used and first validation tests of such a solver. The AVIP code solves plasma equations in complex industrial geometries using an unstructured parallelefficient 3D fluid methodology. AVIP also includes a Particle-In-Cell (PIC) solver used as a reference for validation. While full 3D PIC simulations of a Hall thruster still require unaffordable computational time, fluid models provide in a reasonable time 3D results on the plasma behavior inside the discharge channel. In this category, standard drift-diffusion models [START_REF] Hagelaar | Role of anomalous electron transport in a stationary plasma thruster simulation[END_REF][START_REF] Bareilles | Critical assessment of a two-dimensional hybrid Hall thruster model : Comparisons with experiments[END_REF][START_REF] Lam | A 2-D Hybrid Hall Thruster Simulation That Resolves the ExB Electron Drift Direction[END_REF] are fast and robust but at the cost of strong hypotheses and simplifications. In particular such models do not describe explicitly the sheath formation in the vicinity of walls and often use analytical models instead. They are limited to simple configurations and only provide a first insight into plasma complex phenomena. The present approach includes a more detailed two-fluid plasma model without drift-diffusion approximation. After the description of the formulation and main features of the solver, the paper focuses on wall boundary conditions which are crucial for the formation of sheaths. It is demonstrated in particular that a vacuum boundary condition is not adapted to fit PIC results. A boundary condition based on wall thermal fluxes is more realistic. The mesh resolution is also found to be critical. The simulation methodology is finally applied to a 2D simulation of a typical Hall effect thruster in order to observe the plasma properties inside the discharge chamber.

I. Nomenclature

II. Introduction

H effect thrusters have been massively studied since the 1970s for satellite propulsion. They are still widely used to cover a range of applications from satellite station keeping to orbital transfer and positioning. However, the magnetised plasma physics inside these thrusters is particularly complex. Simple analytical models are used for the design of operating parameters suited to the requested mission. Once constructed, the device undergoes several experimental qualifications to demonstrate the thruster efficiency and life time of the device during more than 10,000 hours in a vacuum test facility. These costly and long tests are required because of the lack of insight in Hall thrusters plasma behavior. With the rise of high performance computing, numerical simulation has become attainable and critical for the design of new thrusters in the industrial context.

Accuracy and reliability are of course essential and different numerical plasma models with various orders of accuracy have emerged. On one side, Particle-In-Cell (PIC) models follow the trajectories of a representative number of physical particles called macro-particles in the presence of electromagnetic fields [START_REF] Birdsall | Particle-in-cell charged-particle simulations, plus Monte Carlo collisions with neutral atoms, PIC-MCC[END_REF][START_REF] Taccogna | Kinetic simulations of a plasma thruster[END_REF]. They are often coupled to a Monte Carlo module to deal with statistical collisions and a Maxwell or Poisson solver to solve electromagnetic fields. On the other side, fluid models describe the statistical behavior of each population of particles. The fluid equations are derived for the statistical moments from the Boltzmann equation supposing a Maxwellian distribution of particles.

Assumptions can be made such as quasineutrality and equality of ion and electron density, or stationary state neglecting temporal and convective terms, in the so-called drift-diffusion models.

In this paper, a two-fluid detailed model for electrons and ions is presented, without quasineutrality nor drift-diffusion hypotheses. Similarly to [START_REF] Shumlak | Approximate Riemann solver for the two-fluid plasma model[END_REF][START_REF] Loverich | A discontinuous Galerkin method for the full two-fluid plasma model[END_REF][START_REF] Hakim | Extended MHD Modelling with the Ten-Moment Equations[END_REF][START_REF] Futtersack | Modélisation fluide du transport magnétisé dans les plasmas froids[END_REF][START_REF] Alvarez Laguna | Fully-implicit finite volume method for the ideal two-fluid plasma model[END_REF], each species is ruled by a system of Euler equations with source terms representing electromagnetic forces, ionization and other collision processes. The computation of collision source terms follows the idea of Benilov [START_REF] Benilov | Momentum and energy exchange between species of a multicomponent gas mixture due to inelastic and reactive collisions[END_REF][START_REF] Benilov | A kinetic derivation of multifluid equations for multispecies nonequilibrium mixtures of reacting gases[END_REF]. A Poisson equation is solved for the electric potential. The set of equations has been implemented in the AVIP-Fluid solver, an unstructured parallel-efficient 3D code for the numerical resolution of PDEs. AVIP-Fluid has been validated against several benchmark cases reproducing the main characteristics of a cold magnetised plasma such as plasma oscillation or collision processes. A particular attention is paid on wall boundary conditions as the creation of stable sheaths with correct properties is a difficult challenge for non-quasineutral fluid models. The validated model is then applied to a simplified 2D r-z configuration of Hall thruster. The objective is to analyse the discharge inside a Hall thruster and to demonstrate the capacity of the code for this specific application. Section III describes the system of fluid equations solved in AVIP. Then, the numerical methods of AVIP are briefly presented in Section IV with particular focus on the Poisson solver which is a key parameter for good computational performances. In Section V the sheath formation in a bounded plasma is studied and a numerical strategy to capture them with AVIP-Fluid is detailed. A comparison with AVIP-PIC results is performed to examine the influence of the boundary condition on the sheath formation. Finally, the Section VI shows results about the representative case for Hall effect thrusters.

III. Physical modeling

In this section, the fluid model used in AVIP-Fluid is presented. It includes conservative equations for electrons, ions, and neutrals.

Starting from the Boltzmann equation, the evolution of the distribution function f α (ì v, ì x, t) of a plasma species α writes :

df α dt + ì v • ì ∇ ì x f α + ì F α m α ì ∇ ì v f α = df α dt collisions (1)
where ì F α are the external forces exerted on the species α of mass m α . A distinction is done between gradient operators ì ∇ ì x in physical space and ì ∇ ì v in phase space, the phases being the velocity components. The left hand side of Eq. 1 represents the collisionless behavior of the species distribution function with convective and time derivative terms and effect of external forces. The right hand side describes the collisional processes.

In a plasma under electromagnetic fields, ì F α is the Lorentz force :

ì F α = q α ( ì E + ì v × ì B) (2)
with q α the species electric charge and ì E and ì B the electric and magnetic fields respectively. In Hall effect thrusters the magnetic field ì B induced by the system is considered negligible compared to the magnetic field imposed by the coils so that ì B is taken constant : ì B = ì B(x). The time evolution of the electric field ì E is no longer governed by the full Maxwell equations but reduces to :

ì ∇. ì E = 1 0 α q α n α (3)
with 0 the vacuum permittivity and n α the number density of the species α. In general the resolution of this equation is done through the electric potential φ in the domain :

ì E = -ì ∇φ (4)
which leads to the Poisson equation :

∆φ = 1 0 α q α n α (5)
Macroscopic properties of each species α correspond to the statistical velocity moments φγ α :

ρ α = m α n α = m α ∫ f α dì v = m α φ0 α (6) m α n α ì u α = m α ∫ ì v f α dì v = m α φ1 α (7) α = 3 2 k B n α T α + 1 2 m α n α ì u α 2 = m α 2 ∫ ì v 2 f α dì v = m α φ2 α (8)
Here ρ α is the mass density proportional to the number density n α coming from the 0 th order moment φ0 α = n α . The velocity ì u α is related to the 1 st order moment φ1 α and the total energy α is the second order moment φ2 α which includes both kinetic energy and internal species energy, expressed with the mean temperature T α . k B is the Boltzmann constant.

Integrating the Boltzmann equation in the phase space for each species α and omitting collision terms, the collisionless transport equations are obtained :

∂ t ρ α + ∇.(ρ α ì u α ) = 0 (9) ∂ t (ρ α ì u α ) + ì ∇.(ρ α ì u α ì u α + p α I) = q α n α ( ì E + ì u α × ì B) (10)
∂ t ( α ) + ì ∇. 1 2 ρ α ì u 2 α + γ γ -1 p α ì u α + ì Q α = q α n α ì E. ì u α (11)
These equations are valid for each type of particles, being electrons, ions or neutrals. The fluid pressure p α , and the heat flux ì Q α come from higher order moments appearing in the Boltzmann equation and need to be expressed in order to close the system.

The heat flux ì Q α is the third order moment of velocity. An exact closure can be found by integrating even higher order moments of the equations as in [START_REF] Dujko | High-order fluid model for streamer discharges : I. Derivation of model and transport data[END_REF][START_REF] Markosyan | High-order fluid model for streamer discharges : II. Numerical solution and investigation of planar fronts[END_REF]. However this only moves the closure problem to higher order moments. Some authors solve a more sophisticated time-dependent equation for the heat flux depending on the magnetic field anisotropy [START_REF] Hagelaar | Modeling methods for low-temperature plasmas[END_REF]. This approach is not considered in the present fluid model because it introduces three additional equations in 3D to be solved and consequently additional computing time. Instead the system is closed with a simple Fourier law for the heat flux. Previous works [START_REF] Cagas | Continuum Kinetic and multi-fluid simulations of classical sheaths[END_REF] have revealed the importance of the electron heat flux in the vicinity of sheaths while the ion heat flux can be considered negligible. The Fourier law for electrons writes :

ì Q e = γk 2
B n e T e m e n n K en ì ∇T e [START_REF] Benilov | A kinetic derivation of multifluid equations for multispecies nonequilibrium mixtures of reacting gases[END_REF] with γ the heat capacity ratio of the fluid and K en the electron-neutral elastic collision rate. Ions and electrons are considered as monoatomic particles for which γ = 5/3 in three dimensions.

For the pressure p α , a perfect gas assumption for each species is made :

p α = n α k B T α (13)
Finally the 10-moment (or 2-fluid) collisionless system for electrons and ions reads :

∂ t ρ e + ∇.(ρ e ì u e ) = 0 ( 14
)
∂ t (ρ e ì u e ) + ∇.(ρ e ì u e ì u e + k B T e n e I) = -en e ( ì E + ì u e × ì B) (15)
∂ t ( e ) + ∇. ( 1 2 ρ e ì u 2 e + γ γ -1 k B T e n e ). ì u e = -en e ì E. ì u e (16)
∂ t ρ i + ∇.(ρ i ì u i ) = 0 (17) ∂ t (ρ i ì u i ) + ∇.(ρ i ì u i ì u i + k B T i n i I) = en i ( ì E + ì u i × ì B) (18)
∂ t ( i ) + ∇. ( 1
2 ρ i ì u 2 i + γ γ -1 k B T i n i ). ì u i = en i ì E. ì u i (19)
Equations 14 and 17 represent the conservation of ion and electron mass density. Equations 15-16 and 18-19 deal with the conservation of momentum and energy for electrons and ions respectively.

Equation systems [START_REF] Markosyan | High-order fluid model for streamer discharges : II. Numerical solution and investigation of planar fronts[END_REF][START_REF] Hagelaar | Modeling methods for low-temperature plasmas[END_REF][START_REF] Cagas | Continuum Kinetic and multi-fluid simulations of classical sheaths[END_REF] and [START_REF] Chen | Fluid modeling of electron heaheat in low-pressure, high frequency cacapacitive coupled plasma discharges[END_REF][START_REF] Barral | Numerical studies of Hall thrusters based on fluid equations for plasma[END_REF][START_REF] Escobar | On conditions at the sheath boundaries of a quasineutral code for Hall thrusters[END_REF] without external forces are known as Euler fluid equations. In a plasma, the two sets of Euler equations are coupled through the resolution of the electric field via the Poisson equation and collisions must be added. Collisions between two particles of the same population are considered negligible compared to collisions between different particles. Ionization plays a major role in a Hall effect thruster discharge and needs to be modelled accurately. In the present framework, only elastic collisions of charged particles (ions and electrons) with neutrals are taken into account and Coulomb collisions between charged particles are neglected. Excitation processes are also taken into account to prevent electron energy losses. However we choose not to represent the influence of these processes on neutrals.

With these hypothesis, the system of equations (14-19) becomes :

∂ t ρ e + ∇.(ρ e ì u e ) = S 0 ioniz,e (20)
∂ t (ρ e ì u e ) + ∇.(ρ e ì u e ì u e + k B T e n e I) = -en e ( ì E + ì u e × ì B) + S 1 en,e (21)
∂ t ( e ) + ∇. ( 1 2 ρ e ì u 2 e + γ γ -1 k B T e n e ). ì u e = -en e ì E. ì u e + S 2 en,e + S 2 ioniz,e -S 2 exc,e (22)
∂ t ρ i + ∇.(ρ i ì u i ) = S 0 ioniz,i (23)
∂ t (ρ i ì u i ) + ∇.(ρ i ì u i ì u i + k B T i n i I) = en i ( ì E + ì u i × ì B) + S 1 in,i (24)
∂ t ( i ) + ∇. ( 1
2 ρ i ì u 2 i + γ γ -1 k B T i n i ). ì u i = en i ì E. ì u i + S 2 in,i + S 2 ioniz,i (25)
where S γ δ,α is the source term representing, for each species α, the collision process δ ( ioniz for ionization, en and in for electron-neutral, ion-neutral elastic collision, exc for excitation). Collision frequencies are generally assumed constant or derived from an analytical solution [START_REF] Chen | Fluid modeling of electron heaheat in low-pressure, high frequency cacapacitive coupled plasma discharges[END_REF]. In the present paper, collision source terms are computed following the work of Benilov [START_REF] Benilov | Momentum and energy exchange between species of a multicomponent gas mixture due to inelastic and reactive collisions[END_REF][START_REF] Benilov | A kinetic derivation of multifluid equations for multispecies nonequilibrium mixtures of reacting gases[END_REF], which integrates cross sections of each considered collision depending on the species temperature and mean velocity.

Finally, like in several Hall effect thrusters models [START_REF] Barral | Numerical studies of Hall thrusters based on fluid equations for plasma[END_REF][START_REF] Escobar | On conditions at the sheath boundaries of a quasineutral code for Hall thrusters[END_REF], neutrals are supposed to be only governed by a mass conservation equation :

∂ t ρ n + ì u 0,n ∇.(ρ n ) = S 0 ioniz,n (26)
where u 0,n = [200 -300] m.s -1 is a constant axial velocity . The electron-impact ionization is represented by the source term S 0 ioniz,α in mass continuity equations 20, 23 and 26. The corresponding source terms read :

S 0 ioniz,e = S 0 ioniz,i = -S 0 ioniz,n = K ioniz n n n e (27)
with K ioniz the ionization rate coefficient depending on the cross section and the electron energy distribution function. S 2 ioniz,e and S 2 ioniz,i represent the energy transfer between particles associated to the ionization :

S 2 ioniz,i = 3 2 K ioniz n n n e k B T n , S 2 ioniz,e = -n n n e K ioniz E * ioniz ( 28
)
with E * ioniz the energy of the ionization reaction. It assumes that only electrons give their energy for the reaction and neutral energy is transferred to ions. The particle velocity deviation due to ionization is considered negligible in front of elastic collisions.

S 2 exc,e in Eq. 22 is the electron energy used to excite neutrals. Several excitation processes are considered :

S 2 exc,e = - N j n n K exc, j E * exc, j (29)
with E * exc, j the energy required to excite a neutral at the excited state j and K exc, j the corresponding reaction rate. Elastic collisions impact the momentum of particles without creation or exchange of particles, so that S 0 elastic,α = 0 . S 1 en,e and S 1 in,i represent the deviation of particle velocity due to elastic collisions with neutrals :

S 1 en,e = -n n n e K en (u e -u n ) , S 1 in,i = -n n n i K in (u i -u n ) (30)
S 2 en,e and S 2 in,i express the exchange of energy between particles during an elastic collision :

S 2 en,e = - 3 2 K en n n n e k B (T e -T n ) + 1 2 (m e u 2 e -m n u 2 n ) , S 2 in,i = - 3 2 K in n n n i k B (T i -T n ) (31)

IV. Implementation and numerical strategy

Inspired by efficient CFD codes [START_REF] Gourdain | High performance parallel computing of flows in complex geometries : I. Methods[END_REF][START_REF] Moureau | Design of a massively parallel CFD code for complex geometries[END_REF], AVIP is a massively parallel code based on a double domain decomposition. First the domain is split in several subdomains using the external library PARMETIS [START_REF] Karypis | MeTis : Unstructured Graph Partitioning and Sparse Matrix Ordering System[END_REF]. This decomposition allows a parallel use on a multitude of processors using the Message Passing Interface (MPI). Then each subdomain is divided into cell groups for an easy optimization of the use of processor memory. The AVIP solver handles unstructured 2D triangles and 3D tetrahedral meshes. This ability allows to deal with complex geometries and facilitates mesh adaptation and grid generation. The numerical resolution of the fluid equations is based on an operator splitting methodology for unstructured meshes. Convection and source terms are integrated separately in time and ensure a first order accuracy resolution. A Strang splitting [START_REF] Abgrall | Robust Finite Volume Schemes for Two-Fluid Plasma Equations[END_REF] is also available to get a second-order accurate splitting method. Convective terms are integrated spatially with a finite volume method using a cell-vertex method, with dual cells adapted to unstructured basis. Fluxes at volume interfaces are computed using various Riemann solvers. The HLLC Riemann solver [START_REF] Toro | Restoration of the contact surface in the HLL-Riemann solver[END_REF] appears to be robust and accurate for the present application and is used in the following. Time integration is performed with a third-order Runge Kutta explicit time scheme in order to avoid large numerical time advancement errors. All source terms (Lorentz forces and collision operators) are solved at once with an implicit time resolution using an external open-source solver from the library ODEPACK [START_REF] Hindmarsh | ODEPACK, A Systematized Collection of ODE Solvers[END_REF]. The use of an implicit scheme allows larger time steps especially for the ionization source time integration. Fluid boundary conditions will be detailed in Section V. The Poisson equation has elliptic properties and is the most time-consuming step in AVIP. A parallel-efficient solver, the Maphys software [START_REF] Agullo | Robust coarse spaces for Abstract Schwarz preconditioners via generalized eigenproblems[END_REF][START_REF] Agullo | On the parallel design of coarse space corrections for sparse hybrid solvers[END_REF], developed by the HIEPACS team at INRIA Bordeaux, is used. Maphys is based on a hybrid approach using the same numerical partitioning as the fluid computation. It combines the accuracy of direct methods on each processor and the efficiency of iterative methods to treat the interface system between cores.

V. Adapted boundary conditions for plasma sheaths

The length of sheaths is typically a few Debye lengths in a Hall effect discharge chamber. The very thin layer (generally few hundreds of µm) is crucial for the global behavior of the plasma inside the whole Hall effect thruster. In particular, plasma instabilities triggered at this location, may lead to wall erosion [START_REF] Sydorenko | Plasma-sheath instability in Hall thrusters due to periodic modulation of the energy of secondary electrons in cyclotron motion[END_REF][START_REF] Minelli | How to Build PIC-MCC Models for Hall Microthrusters[END_REF] and have a significant impact on the Hall thruster life time. Therefore the evaluation of the electric field and plasma densities inside the sheath is essential to understand and predict the efficiency of a Hall thruster. The physical mechanism of sheaths is fully understood and explained in the literature [START_REF] Anderson | The formation of a plasma sheath[END_REF][START_REF] Takamura | Sheath formation in a plasma with energetic electrons[END_REF][START_REF] Scheiner | Theory of the Electron Sheath and Presheath[END_REF]. The loss of electrons at the boundary produces a pressure gradient which accelerates them with larger velocities than ions. The subsequent negative charge difference created at the wall, generates an electric field that accelerates electrons in the direction opposite to their motion, i.e., away from the wall. The balance between electric force and pressure gradient leads to a stationary sheath with a large electric field [START_REF] Langmuir | The effect of space charge and initial velocities on the potential distribution and thermionic current between plane electrodes[END_REF]. The characteristic length of this phenomenon is the Debye length λ D defined as :

λ D = 0 k B T e n e e 2 ( 32
)
where T e is the electron temperature, e the elementary charge and n e the electron density. The Debye length corresponds to the characteristic length of the electric field induced by a charged particle. The influence on the potential of a charged particle is exponentially dependent on the Debye length :

φ = φ 0 e -x λ D (33)
As a consequence, the Debye length should be resolved to correctly represent the influence of each charged particle on the overall field of electric potential.

Another standard result of the sheath theory is the Bohm criterion [START_REF] Riemann | The Bohm criterion and sheath formation[END_REF] which states that the ion velocity u i inside the sheath must be equal or larger than the ion sound speed to obtain a stable sheath :

u i ≥ k B T e m i (34)
The equality between the ion velocity and the ion sound speed defines the limit between the sheath and the presheath. This description is valid as long as λ D << L 0 where L 0 is the characteristic length of the domain. This condition depends on the electron temperature which must be correctly predicted.

For a sheath with fixed potential, the Child-Langmuir model [START_REF] Lieberman | Principles of plasmas Discharges and materials processing[END_REF][START_REF] Chabert | Physics of Radio-Frequency plasmas[END_REF] relates the voltage φ 0 across the sheath with the sheath thickness s, the electron temperature T e and the Debye length λ D . For a collisionless sheath, s reads :

s = 4 √ 2 9 1 2 eφ 0 k B T e 3 4 λ D ( 35
)
This relation is based on fluid macroscopic properties of ions and electrons, i.e., assuming a Maxwellian distribution function, and neglects the electron space charge.

Contrary to fluid solvers, the boundary condition for creating sheaths near walls in a PIC solver is quite simple. Each particle (ion or electron) exiting the domain is lost at the wall boundary which naturally leads to a stationary sheath at the wall or electrode.

Creating stable sheaths in fluid models is less straightforward. By construction, quasineutral solvers do not solve the sheaths and replace the sheath by equivalent boundary conditions [START_REF] Escobar | On conditions at the sheath boundaries of a quasineutral code for Hall thrusters[END_REF]. Analytical models [START_REF] Ahedo | Presheath/sheath model with secondary electron emission from two parallel walls[END_REF][START_REF] Hofer | Wall sheath and electron mobility modeling in hybrid-PIC Hall thruster simulations[END_REF][START_REF] Hobbs | Heat flow through a Langmuir sheath in the presence of electron emission[END_REF][START_REF] Stephens | Sheath and presheath potentials for anode, cathode and floating plasma-facing surfaces[END_REF][START_REF] Robertson | Sheaths in laboratory and space plasmas[END_REF] for the sheath and the presheath, that impose potential drop and plasma densities are commonly used.

On the other hand, drift-diffusion solvers coupled to a Poisson equation are able in principle to represent sheaths [START_REF] Hagelaar | Boundary conditions in fluid models of gas discharges[END_REF][START_REF] Hagelaar | Model of an inductively coupled negative ion source : I. General model description[END_REF]. The wall flux can be expressed considering that a charged particle hitting a wall is instantaneously neutralized and assuming that no particles are coming from the wall. Particles may exit at the wall due to their thermal motion and not only their mean velocity. This corresponds to the particles of the Maxwellian distribution function that have a positive (or negative) velocity fluctuation.

The present model is a 10-moment fluid formulation with Poisson equation. Imposing boundary conditions on such equation system is a crucial issue to reproduce correctly the sheath formation. Following the work of Cagas and Hakim [START_REF] Cagas | Continuum Kinetic and multi-fluid simulations of classical sheaths[END_REF], it is possible to create sheaths near walls assuming a vacuum just outside the domain for ions and electrons. This means that each charged particle at the wall is instantaneously neutralized and lost.

In [START_REF] Chen | Fluid modeling of electron heaheat in low-pressure, high frequency cacapacitive coupled plasma discharges[END_REF][START_REF] Becker | Advanced fluid modeling and PIC/MCC simulations of low-pressure CCRF discharges[END_REF], wall boundary conditions are expressed using explicit fluxes crossing the boundary, including the flux due to particle thermal motion, as drift-diffusion solvers. In this way, the half of the Maxwellian distribution representing particles that exit the domain, is considered as an outgoing flux.

In the following, a PIC simulation is first performed with the AVIP-PIC solver in order to understand the formation of the sheath and the behavior of macroscopic plasma properties. Note that the AVIP-PIC solver uses a Boris scheme for the transport of particles without collisions and in the presence of an electric field. Secondly, PIC results are considered as the reference and two different boundary conditions with AVIP-Fluid are compared with. Finally the influence of the mesh resolution on AVIP-Fluid results is addressed.

A. Test case on a collisionless sheath : preliminary results with AVIP-PIC

In this section the formation of a 2D collisionless anode sheath in a Xenon gas is studied (m Xe = 2.18 10 -25 kg) with the AVIP-PIC solver. The AVIP-PIC code is used as a reference to highlight the balance between pressure gradient and Lorentz force. The y direction is considered periodic. The anode potential φ 0 at x = 0 is set to 200 V and the cathode potential to φ L = 0V at x = L. No secondary electron emission or other wall interaction are considered. Electron particles are initialized with a Maxwellian velocity distribution function at an initial temperature of T e,0 = 10 eV. Cold ions are initially at T i,0 = 0.1 eV. Initially, 50 macro-particles per cell are set to ensure statistical convergence. Macro-particles weights are initialized in order to get same electron and ion densities of n e = n i = 5.0 10 16 m -3 . A fixed timestep of ∆t = 5 10 -12 s is used for the simulation. The Debye length computed with these initial conditions is λ D = 105.2µm and is the reference length for all simulation cases. After 100 ns the sheath starts to form as evidenced by ions and electrons densities in Fig. 1a. The electron density decreases faster than the ion density due to the high mobility of the much lighter electrons. The thickness of the sheath represents the area where the plasma is not quasineutral and is around 8 Debye lengths at t = 0.1µs. Also shown in Fig. 1a, the electron temperature slightly increases due to the presence of a strong electric field in Fig. 1b. The electric field magnitude increases in the sheath due to the source term proportional to (n in e ) in the Poisson equation, Eq. 5. The potential, fixed at the anode, is shifted in the sheath by 65 V at t = 0.1µs. This drop of potential is higher than the expected value at steady state, as the low mobility of ions induces long transient states.

F 2 Balance between momentum equation forces inside the sheath at t = 0.1µs, AVIP-PIC simulation

Figure 2 represents the different contributions to the electron momentum conservation equation derived from the PIC simulation. The convection term is negligible inside the whole domain and particularly in the sheath. As mentioned in the literature, the sheath results from the balance between pressure gradient and electric forces. Electrons exit the domain and create a negative pressure gradient. This pressure gradient tends to impose a negative electron velocity meaning that the force expulses the electrons out of the domain. The negative electric field created by the charge difference induces a positive Lorentz force and a negative electron velocity. At t = 0.1µs, the electric force is not strong enough to neutralize the pressure gradient and to stabilize the sheath. In the presheath, the Lorentz force is even negative. The steady state is not reached.

At equilibrium the plasma sheath stabilizes with characteristic potential drop and thickness. In the present case this is obtained at t = 0.3µs as shown in Fig. 3. The sum of momentum forces near the wall is slightly negative meaning that electrons still leave the domain. The net force at the wall is around 5 N. In the presheath, the electric force is now positive and compensates the pressure gradient. For an anode sheath, this part of electrons exiting the domain through The electron temperature decreases in the vicinity of the wall due to particles that leave the domain. In Fig. 4b, the maximum of electric field is located at the wall and is around -12, 500 V/m. The value is lower than in the transient state and no longer compensates for the loss of electron temperature at the wall. A potential drop of 3 eV is also observed at the sheath edge.

In Fig. 5, the Bohm criterion (Eq. 34) is verified : the mean ion velocity is well equal to the ion sound speed at the beginning of the sheath and assures a stabilized sheath. At the boundary, the ion velocity is twice larger than the ion sound speed. As expected ions are supersonic and should be treated with caution in the fluid formalism.

B. Wall boundary conditions for 10-moment fluid solvers

Reproducing sheaths in fluid solvers is a known difficulty. The configuration previously studied with AVIP-PIC is now computed with AVIP-Fluid. A refined mesh with 100 cells per Debye length in the sheath area is used which F 5 Bohm criterion : ion velocity compared to ion sound speed, AVIP-PIC simulation corresponds to a cell size of approximately 1µm near the wall. The HLLC Riemann solver is applied with a CFL number of 0.1. In these conditions, we suppose that numerical errors are negligible and that the results depend only on the boundary methodology and the physical model.

Two different types of boundary conditions are used and tested at the anode. First, the vacuum condition is implemented as in the paper of Cagas [START_REF] Cagas | Continuum Kinetic and multi-fluid simulations of classical sheaths[END_REF]. Numerically, the Riemann solver at the boundary condition computes a self-consistent flux between two fluxes : a zero density and pressure flux corresponding to a vacuum on one side (called the "ghost cell"), and the flux predicted by the scheme on the other side. The solution at the boundary condition adjusts itself to compute the surface flux adapted to the predicted variables. To avoid numerical instabilities, the p = 0 condition is not imposed in a hard way, but through a relaxation procedure :

p ghostcell = p in 100 ( 36
)
where p in is the predicted pressure at the boundary before the boundary condition procedure. The factor η p = 1/100 is arbitrary in the model and will be discussed later. Figure 6 shows results with the vacuum condition at t = 0.3 µs. The transient phase is not presented here as it is similar to the PIC simulation. On Fig. 6b, the potential drop is overestimated and reaches 40 V at the exit of the sheath. Consequently the electric field is higher inside the sheath and reaches -110, 000 V/m, which is 9 times the value of the PIC simulation. This huge electric field induces a larger sheath. On Fig. 6a, the steady sheath thickness is equal to around 10 Debye lengths, i.e., largely overestimated when compared to AVIP-PIC results. The sheath is also larger than predicted by the Child-Langmuir theory. In fact, the electron density induced by the vacuum condition is lower at the wall and leads to a difference of charge densities.

The electron temperature decreases inside the sheath and is however correct according to PIC results of Fig. 4. A sensitivity study showed no impact of the relaxation ratio η p on these results.

The second boundary condition specifies the thermal flux at the boundary. The thermal velocity represents the mean thermal motion of particles α :

v th,α = 2k B T α m α (37)
The thermal flux does not impact the mean particle velocity, i.e., the convective flux, but applies to the mass conservation [START_REF] Becker | Advanced fluid modeling and PIC/MCC simulations of low-pressure CCRF discharges[END_REF] :

F 0 th,α = n α,in v th,α 2 √ π (38)
with n α,in the density predicted by the scheme inside the domain and v th,α the thermal velocity of the species α. This flux represents the particles that leave the domain because of their velocity fluctuation and can be derived from the Maxwellian distribution function, see [START_REF] Hagelaar | Model of an inductively coupled negative ion source : I. General model description[END_REF]. Considering that the mass flux leaving the domain is at the temperature computed inside the domain T α,in , the corresponding flux for the energy conservation reads :

F 2 th,α = F 0 th,α • k B T α,in m α (39)
The components F γ th,α are applied to the corresponding equation of ( 20)- [START_REF] Hindmarsh | ODEPACK, A Systematized Collection of ODE Solvers[END_REF]. Moreover the flux predicted by the scheme inside the domain is kept if the wall normal component of the velocity is positive. If the wall normal velocity component is negative, meaning an incoming flux, the flux is set to zero. This can be synthesized as :

F γ bndy,α • ì n ⊥ = max F γ in,α • ì n ⊥ , 0 + F γ th,α (40)
where ì n ⊥ is the boundary normal vector . Results obtained for the thermal flux condition at t = 0.3 µs, are shown in Fig. 7. The sheath thickness is around 6 Debye lengths which is in better agreement with AVIP-PIC results. The electron density at the wall boundary is much higher compared to the vacuum condition and corresponds to PIC results. This leads to a lower difference of charge density. Consequently, the potential drop is equal to 10 eV and is slightly overestimated but closer to PIC results compared to the vacuum condition. Note that this drop has not yet reached a steady value as the transient regime is a bit longer for fluid simulations. Here, the electric field maximum is around -38, 000 V/m, much lower than for the vacuum condition. Finally the temperature drops in the sheath thickness and is in agreement with both previous simulations.

As explained in Section A, the sheath reaches a steady state thanks to the balance between pressure gradient and Lorentz force. The analysis of these forces for the two different fluid boundary conditions is made at t = 0.3µs on Fig. 8. For the vacuum condition, Fig. 8a, we observe a zero value of all momentum forces at the wall. Imposing a zero pressure outside the domain leads to a zero pressure on the two first Debye lengths inside the domain by continuity. The pressure gradient is close to zero such as the electric contribution to respect the force balance inside the sheath. This also explains the low sensitivity of the results to the pressure ratio η p , which only drives the convergence rate to the steady state but does not change the final pressure inside the domain. For the thermal flux condition, Fig. 8b, the sum of momentum forces at the wall is negative, around 2 N. Electrons leave the domain as already shown in PIC results. The pressure gradient is not null but results from the predicted flux.

It is clear from these results that the vacuum condition properties are not conceptually similar to the PIC boundary conditions. The low electron density imposed at the boundary because of the vacuum pressure creates a larger sheath than the PIC simulation.The thermal flux condition, intuited from the microscopic particle behavior at the wall, is more suited to match PIC results.

In the present problem, the Debye length (Eq. 32) appears to be a key quantity which controls main sheath properties. On Fig. 9, the Debye length is plotted in the vicinity of the wall and compared to the AVIP-PIC results. The Debye length slightly increases inside the sheath because of a lower electron density. The thermal flux condition shows a very good evaluation of the Debye length in the sheath/presheath area, whereas the vacuum condition leads to large errors in the sheath region and close to the wall. This is to be related to the too low imposed electron pressure and consequently too low density at the wall.

Finally the Bohm criterion theory is evaluated in Fig. 10. Ion velocities for AVIP-Fluid with the two proposed boundary conditions, also applied for ions, are normalized by ion sound velocities and compared to the AVIP-PIC results. The ratio of ion velocity to ion sound speed of the AVIP-Fluid with the thermal flux condition is closer to the AVIP-PIC profile than the AVIP-Fluid with vacuum condition. Inside the presheath, the ion velocity is however slightly overestimated. The Bohm criterion is respected for the thermal flux condition : the sheath thickness (estimated at around 6 Debye lengths) corresponds to the value of ion velocity equality. However, the Bohm criterion is not valid for the AVIP-Fluid with the vacuum condition. The sheath thickness was estimated around 10 Debye lengths in this case and the ion velocity is not over the ion sound speed between 5 and 10 Debye lengths from the wall. The ion velocity in the sheath is largely overestimated due to the strong pressure condition at the boundary.

Table 1 summarizes the results of each simulation with characteristic properties of the sheaths.

C. Influence of mesh resolution

Meshing near-wall regions with 100 cells per Debye length is not a viable solution for the simulation of complex 3D configurations since this would lead to too small computational time steps. The spatial resolution in the sheath is therefore a key factor and is studied in this subsection. Two meshes with respectively 2 cells per Debye length and 8 cells per Debye length inside the sheath are tested. Note that the characteristic Debye length chosen here is estimated away from the sheath, and equal to 105.2 µm, see Fig. 9. The HLLC Riemann solver of order one is used with the thermal flux boundary condition. Results with these two meshes are plotted on Fig. 11.

A too low resolution of the near-wall region generates a peak of electron temperature. With two cells per Debye length, the electron temperature reaches 120 eV at t = 0.3 µs and continues to grow. The sheath is not stable and much larger (more than 25 Debye lengths). The reason is that the Bohm criterion stability is not fulfilled because the high peak of electron temperature induces a higher ion sound speed. With 8 cells per Debye length in the near-wall region, the sheath is created with a slight heating of about 8 eV and the sheath thickness is overpredicted (around 11 Debye . From these results, the minimal resolution can be estimated at around 8 to 10 cells per Debye length. Note that this value is clearly dependent on the numerical scheme used (order 1 here). It means that an effort on the mesh resolution need to be done in order to correctly reproduce the sheath behavior near walls and electrodes.

VI. A 2D rz test case

The physical model has been validated on simple test cases in order to assess of the plasma modeling, from standard plasma oscillation to detailed collision processes (not shown here for the sake of brevity). This final section presents the ability of the Fluid AVIP solver to simulate a Hall effect thruster configuration. For this purpose, the SPT 100 model which has been widely studied in the literature [START_REF] Sullivan | Kinetic Study of Wall Effects in SPT Hall Thrusters[END_REF][START_REF] Escobar | Global stability analysis of azimuthal oscillations in Hall Thrusters[END_REF][START_REF] Boeuf | Tutorial : Physics and modeling of Hall thrusters[END_REF] is considered. A 2D r-z configuration of this thruster is chosen in order to get the global behavior of the plasma in the axial direction. Since the plume plays a role in the optimal operation of a Hall thruster, it is included in the computational domain. The main geometrical and operating parameters of the SPT 100 are referenced in the paper presented by Escobar [START_REF] Escobar | Global stability analysis of azimuthal oscillations in Hall Thrusters[END_REF] and are listed in Table 2 m n is the Xenon mass flow rate at the anode and u 0,n is the constant neutral velocity. The anode inside the Hall thruster has a fixed potential of φ a = 300V. The cathode, set to φ c = 0V is positioned at a distance L AC at the right boundary of the domain in order to get discharge between two plates inside the chamber. A source term at 1 mm from the cathode injects electrons to compensate for the electronic current lost at the anode. The thermal flux boundary condition described in the previous section is used for ions and electrons. The axial magnetic field profile fits the SPT-100 magnetic configuration [START_REF] Escobar | Global stability analysis of azimuthal oscillations in Hall Thrusters[END_REF] and reaches a maximum B max = 237 G. Dimensions of the channel L channel and h channel are described in Table 1.

Figure 12 shows the axial profile of the initial neutral density corresponding to the background solution of the cited paper [START_REF] Escobar | Global stability analysis of azimuthal oscillations in Hall Thrusters[END_REF]. Ion and electron densities are initially constant and set to n e,0 = n i,0 = 10 17 m -3 . With such initial conditions, this leads to an initial Debye length of λ D,0 = 52.6 µm.

The 2D mesh is made of 2, 2M of triangles trying to respect the cell size recommendation made in Section C. In the vicinity of walls, the cell size is fixed to 6 µm in order to have around 9 cells per Debye length in the sheath. Inside the Figure 14 shows preliminary results of the anode sheath for the r-z simulation after 0.1µs. The electron density decreases near wall boundaries to form a sheath with a thickness around 500 µm which corresponds to 10 Debye lengths. The solution is not yet stationary due to the low ion mobility. However it is possible to compute a theoretical thickness using the drift-diffusion law. For a potential of 300 V and an electron temperature of 5 eV, the theoretical thickness is around s th = 17λ D . Here, the magnetic field and collision processes are not negligible inside the sheath and the Child-Langmuir law is not correct. Moreover, previous results in Section V have shown a slight overestimation of the analytical theory compared to PIC results. In Fig. 14, the electron temperature is overpredicted inside the sheath due to numerical thermal heating. Using a slightly more refined mesh at the wall should reduce this thermal heating at wall boundaries.

In Fig. 15, preliminary results of plasma properties inside the discharge chamber and the plume are shown at t = 0.1µs. The formation of sheaths near chamber walls is visible on the electron density, see Fig. 15a. At the exit plane, in Fig. 15b the axial electric field is maximal leading to a strong acceleration of Xenon ions. This phenomena is standard in Hall effect thrusters and is responsible for the engine thrust.

VII. Conclusion

AVIP is a 3D unstructured parallel-efficient solver developed to simulate bounded plasmas in Hall thruster configurations. The AVIP-Fluid model takes into account separated fluid entities for ions, electrons and neutrals with inertia and convection effects, elastic and inelastic collision processes. The equation system is closed using a Fourier law for the heat flux and coupled to the Poisson equation for the electric field resolution. The numerical resolution of these equations is done with a finite volume method coupled to a Riemann solver for convection terms adapted to unstructured meshes. Source term integration is done with an open source implicit solver. The resolution of the Poisson equation is adapted to highly parallel computation using the MAPHYS solver.

Inside Hall effect thrusters, the behavior of sheaths plays a crucial role in the thruster efficiency and lifetime. The study of sheaths raises the issue of explicit fluid boundary conditions. The vacuum condition presented in the paper of Cagas [START_REF] Cagas | Continuum Kinetic and multi-fluid simulations of classical sheaths[END_REF] fails to recover sheath properties compared to PIC results. The thermal flux condition is more in adequacy with PIC results, as it takes into account a truncated distribution function at the boundary condition. The mesh resolution is also a key point to obtain accurate sheath behavior without numerical thermal heating. The minimal mesh resolution for a numerical scheme of order one, in the vicinity of walls, is estimated to 8 cells per Debye length.

Finally the sheath theory is verified on a more realistic configuration of a Hall thruster. A 2D r-z simulation of a SPT 100 is performed on a mesh satisfying previous recommendations. Sheaths are formed near walls with a thickness around 10 Debye lengths. Preliminary results show the global behavior of the plasma characteristics at the exit plane, i.e., a strong electric field and a high ion velocity.

  of the species α η p = pressure ratio for the vacuum condition ì F α = external forces on the species α F γ th,α = thermal fluxes applied on the boundary condition for the species α f α = distribution function of the species α k B = Boltzmann constant K δ = rate coefficient of the collision δ L 0 = characteristic length of the domain λ D = Debye length m α = mass of the species α m α = mass rate of the species α n α = density of the species α p α = pressure of the species α φ = potential φγ α = Boltzmann moment of order γ of the species α q α = electric charge of the species α ì Q α = heat flux of the species α ρ α = mass density of the species α s = sheath thickness S γ δ,α = fluid source term of the species α representing the collision δ t = time variable T α = temperature of the species α ì u α = mean velocity of the species α ì v = velocity vector v th,α = thermal velocity of the species α ì x = space vector

F 1

 1 Characterization of the plasma inside the anode sheath at t = 0.1µs, AVIP-PIC simulation

F 3

 3 Balance between momentum equation forces inside the sheath at t = 0.3µs, AVIP-PIC simulation the electrode supplies the electric circuit and creates a current. Away from the sheath the balance between these forces still play a role but the electric field is much reduced as the charge difference is smaller.

F 4

 4 Characteristics of the plasma inside the anode sheath at t = 0.3µs, AVIP-PIC simulation At equilibrium, Fig. 4a shows that the sheath thickness is strongly reduced down to 4 Debye lengths. In the present configuration, it leads to a sheath thickness of around 0.4mm. Note that he collisionless Child-Langmuir model (Eq. 35) based on fluid approximations, with an electron temperature of T e = 10 eV and a potential of φ 0 = 200V, predicts a sheath thickness of s th = 7.5λ D .

F 6

 6 Characteristics of the plasma inside the anode sheath at t = 0.3µs, AVIP-Fluid simulation with vacuum condition

F 7

 7 Characteristics of the plasma inside the anode sheath at t = 0.3 µs, AVIP-Fluid simulation with thermal flux condition

F 8

 8 Balance between electron momentum equation forces inside the sheath at t = 0.3µs, AVIP-Fluid simulation

F 9

 9 Debye length inside the domain for the various simulation cases F 10 Bohm criterion : Ratio of the ion velocity to the ion sound speed for all simulation cases T 1 Characteristics of the sheath for all simulation cases at t = 0.3µs Simulation case sheath thickness s potential drop φ D Bohm criterion respected ? AVIP-PIC 4 λ D 2 V yes AVIP-Fluid with vacuum condition 10 λ D 40 V no AVIP-Fluid with thermal flux condition 6 λ D 10 V yes Child-Langmuir law 7.5 λ D / /

  electron temperature of the plasma inside the anode sheath at t = 0.3µs for two different meshes, AVIP-Fluid simulations lengths)

F 12 InitialF 13

 1213 Xenon neutral density corresponding to the background solution of Escobar paper [46] Mesh and geometry of the r-z configuration discharge chamber and away from the walls, cell sizes are around 50 µm to get 1 cell per Debye length.

F 14

 14 Densities and electron temperature inside the anode sheath at t = 0.1µs for the r-z configuration (a) Electron density (b) Axial electric field and ion velocity F 15 Plasma characteristics inside the Hall thruster chamber at t = 0.1µs

. T 2 Main parameters of the SPT 100 for the 2D r-z simulation

	m n	4.85 mg/s	u 0,n	300 m/s
	B max	237 G	φ a	300 V
	L AC	33 mm	L channel	25 mm
	h channel	15 mm	n 0,e	10 17 m -3
	T 0,e	5 eV	T 0,i	0.1 eV

Acknowledgments

This work is financially supported by Safran Aircraft Engines (former SNECMA) as a part of a CIFRE convention and supervised by


Arnaud Tom Weber

Geoffroy Lamarche

Yves Le Gonidec

Vanessa Lucieer

Yoann Ladroit

Gas Bubble Forensics Team Surveils the New Zealand Ocean

An international research group recorded the acoustic signatures of gas bubbles rising from a hydrothermal vent field to gather clues about greenhouse gases escaping into the atmosphere.

A synthetic seep generator, used to calibrate echo sounding data from naturally emitted gas bubbles, is deployed into New Zealand's Bay of Plenty. A recent study of bubble streams emitted from hydrothermal vents beneath the bay's surface has provided a foundation for developing methodologies to identify and characterize gases (including greenhouse gases) emitted from the ocean floor using solely their acoustic response.

In July 2018, our research group set sail for the Calypso hydrothermal vent field in New Zealand's Bay of Plenty aboard the New Zealand flagship R/V Tangaroa. The voyage was dubbed Quantitative Ocean-Column Imaging (QUOI). Our mission: Find bubbles-lots of them-and record everything they told us about themselves.

The vent field (Figure 1) did not disappoint: It emitted diffuse streams of carbon dioxide and discrete seeps of methane. Our ship's echo sounders sent out acoustic waves, and we recorded the echoes that these bubbles generated in response. The ways that these echoes varied, depending on the depth of the bubbles and the angles and frequencies of the incoming acoustic waves, provided a wealth of information.

We wondered whether we could use acoustics to fingerprint bubbles in the ocean. The world's oceans are full of bubbles (https://eos.org/articles/audio-reveals-sizes-of-methane-bubbles-risingfrom-the-seafloor). They come from natural sources like dissolving organic matter and gases seeping from great depths below the seafloor and from human sources like leaky oil and natural gas wells and pipelines. Many of these bubbles rise to the surface and release greenhouse gases like carbon dioxide (CO 2 ) and methane (CH 4 ) into the atmosphere. Some research expeditions have collected preliminary information on bubble streams that they discovered by chance while using echo sounders to map the morphology of the seafloor or search for fish aggregates and even using seismic equipment for subseafloor investigations. Our group wanted to dive deeper into the information that ocean bubbles provide. We wondered whether we could use acoustics to fingerprint bubbles in the ocean.

Gathering Clues from Bubbles

Many scientific, industrial, and environmental enterprises make use of the ability to detect, characterize, and quantify bubbles generated from liquid or gaseous features in the ocean, but observations (https://eos.org/science-updates/building-an-accessible-archive-for-water-column-sonar-data) enabling us to measure the volume or identify the types of gases that are released from the seafloor are scarce. Such observations are important as inputs into global carbon cycle models, in helping monitor deep-sea wells and pipelines for potential gas leaks, and in mitigating the associated environmental and economic risks.

Modern echo sounders have been able to image bubble streams in the ocean and trace their movement and behavior as they rise through the water column. To date, analysis of such acoustic data has been largely restricted to localizing and mapping sources of bubbles on the seafloor [e.g., Skarke et al. (https://doi.org/10.1038/ngeo2232), 2014]. Although there has been some recent progress, it has been more difficult to use echo sounders to quantify and characterize bubble streams (numbers of bubbles, geometry, and composition) and, ultimately, the associated volumetric and mass flow of gas [e.g., Greinert and Nützel (https://doi.org/10.1007/s00367-003-0165-7), 2004; Weber et al. Differentiating between gas types such as CH 4 and CO 2 , the two most common greenhouse gases released at the seafloor, solely from their acoustic response remains a major challenge. Recent advances in echo sounders, including in our ability to calibrate and use them over a wide range of frequencies and to visualize and analyze their data, have provided a technological arsenal that we can bring to bear on this challenge.

How Do You Characterize a Bubble?

One of the most remarkable of the numerous gas seepages over the New Zealand continental shelf is the Calypso hydrothermal vent field [Stoffers et al. (https://doi.org/10.1130[Stoffers et al. (https://doi.org/10. /0091-7613(1999))027%3c0931:EMASHV%3e2.3.CO;2), 1999]. Vents emit large bubble streams that can show up on ships' sonars as gas flares (also called gas plumes). Yet many previous surveys only serendipitously imaged the seeps. These images enable qualitative observations, but gas bubble streams cannot be characterized with conventional oceanographic surveys, which provide only partial physical information on individual bubbles.

A synthetic image of the Calypso hydrothermal vent field showing sites where streams of gas bubbles rise from the seafloor. Credit: Erin Heffron Our project enables us to move beyond locating and defining the shape of the flares to assess the bubble size distribution, differentiate between CH 4 and CO 2 gases (both of which are known to be released at the Calypso vents), and eventually develop rapid, cost-effective methods of estimating the bubbles' greenhouse gas flux through the water column.

The intensity of a bubble's acoustic signal-and whether it generates an echo on a ship's sonar screen-is related to the bubble's acoustic scattering cross section, which quantifies the ability of a bubble to reflect the acoustic energy of an incoming acoustic wave back to the echo sounder. This physical (measurable) parameter depends on the geometry, resonance frequency, and damping properties of the bubbles, as well as on the nominal frequency and incident angle of the incoming wave [Clay and Medwin, 1977].

The need to understand how these physical parameters influence acoustic measurements highlights specific requirements for quantitative description of water column gas bubbles. These requirements include calibration of the acoustic devices used to study the bubbles [Ladroit et al. (https://doi.org/10.1007/s11001-017-9340-5) 2002] and how the frequency response of a bubble changes with the ensonification angle (the angle between the sonar beam and the vertical at the point where the beam reaches the bubble), and visual observations for ground truthing. We carried out all of these experiments in 180-meter water depths during the QUOI voyage.

The design of the QUOI voyage included very large overlaps of the sounder seafloor footprints, with up to 95% overlay of adjacent swaths (Figure 3). This enabled us to generate high-resolution images of both the seafloor and the water column, as well as full ensonification of bubbles -immersing them in sound waves with incidence angles from the vertical to a grazing angle of 60°. We were thus able to observe bubble flares from a range of incidence angles. We conducted a robust and complete calibration of all systems at the start of the survey to enable cross correlation between the multiple echo sounder systems. The swath echograms revealed the presence of multiple acoustic flares, and the map of the sum of all echoes within the water column projected on the seafloor highlighted the importance of both the geometry and intensity of the flares and the spatial distribution of the vents on the seafloor (Figure 4).

Our Technological Arsenal

Fully characterizing the bubbles required us to explore different echo sounder deployment methods and data acquisition procedures. We used multiple synchronous echo sounders to assess the backscatter strength of targets at different frequencies: EM302 and EM2040 multibeam echo sounders collected echo amplitudes athwart and along track simultaneously at 30 and 200 kilohertz.

An EK60 vertical echo sounder operated simultaneously at 18, 38, 70, 120, and 200 kilohertz. A wideband EK80 transducer mounted on a rotatable system emitted pulses in the frequency range of 84-270 kilohertz at different incidence angles.

The voyage deployed three technologically and scientifically innovative systems.

The voyage deployed three technologically and scientifically innovative systems. A mobile unit with a split-beam echo sounder that operates at 38 or 120 kilohertz was deployed directly on the seafloor to ensonify bubble streams horizontally. This generated contrasting images from a different point of view-sideways rather than top down. For 5 days, an autonomous hydrophone placed in the center of a large flare field recorded bubble sounds. We also deployed a synthetic seep generator (or bubble maker) on the seafloor to establish a control data set for many of the bubble parameters. The seep generator creates bubbles with known properties and at known rates and sizes. It allows us to develop an empirical model to which data from natural sites can be compared and to develop a quantitative approach to estimate bubble flux rates.

We collected sediment and water samples for ground truthing of the acoustic signals across the entire Bay of Plenty survey site (Figure 2, brown dots). Video transits allowed us to directly observe most of the active vents identified from high-amplitude echo integration of water column acoustic data (Figure 5, large red spot). This original, and multifaceted oceanographic survey has provided us with a data set from which we can develop robust water column surveying methodologies. We anticipate that this will lead to new protocols in identifying, quantifying, and classifying gas bubble streams in the ocean water masses. The QUOI experiments will advance understanding of active geological systems and the relationship the seabed has with the atmosphere.

Fig. 1 .

 1 Fig. 1. The Quantitative Ocean-Column Imaging (QUOI) mission visited the Calypso hydrothermal vent field off the northern coast of New Zealand in the Bay of Plenty. The black and white inset shows a wider view of northern New Zealand, with the area of the larger figure outlined by the black box.

(

  https://doi.org/10.1002/2014GC005271), 2014; Veloso et al. (https://doi.org/10.1002/lom3.10024), 2015; Weidner et al. (https://doi.org/10.1016/j.csr.2018.12.005), 2019].

Fig. 2 .

 2 Fig. 2. The QUOI voyage deployed a variety of acoustic devices over

  , 2018], multifrequency measurements of the backscatter strength of individual bubbles within a plume [Le Gonidec et al. (https://doi.org/10.1029/2001JB000558),

Fig. 3 .

 3 Fig. 3. Sketch of the QUOI acquisition protocol over a bubble seep to assess the frequency and angular dependency of the bubble backscatter strength using multibeam and single-beam echo sounders simultaneously in the range of 18-200 kilohertz. The image demonstrates the high swath overlap on the seafloor.

Fig. 4 .

 4 Fig. 4. Swath echograms recorded by a multibeam echo sounder (illustrated in the inset) over a gray scale map of the seafloor reflectivity (backscatter). The elongated red patches on the seafloor correspond to the vertical echo integration of water column backscatter data: They highlight positions and extensions of flares on the seafloor (processed using the SonarScope software by Institut Français de Recherche pour l'Exploration de la Mer (Ifremer (https://wwz.ifremer.fr/en/))).

Fig. 5 .

 5 Fig. 5. Example of a vertical echo integration in the Calypso hydrothermal vent field. The red patches identify active vent fields.

of 11 27/09/2019 à 06:38

27/09/2019 à 06:38

of 11 27/09/2019 à 06:38

of 11 27/09/2019 à 06:38

Acknowledgments

This project is funded by the Royal Society of New Zealand Catalyst: Seeding fund. The QUOI voyage on board R/V Tangaroa was funded by the National Institute of Water and Atmospheric Research. V.L. was supported by the Marine Biodiversity Hub through funding from the Australian government's National Environmental Science Program.

Author Information

Geoffro Lamarche ( geoffro.lamarche@niwa.co.nz